Contribution to the improvement of thermal comfort in a bioclimatic building by integration of a phase change material

Ya couba Ca m a ra 1 ,2 * , Xa vie r Che sne a u 1 , Ce llou Ka nt e 2

. In addition, we did a comparative study between air temperatures for habitats with and without phase change material in the false ceiling to see the existence of thermal comfort in the habitat which is the goal of this study. The values of these air temperatures observed at 16 h in a modern habitat and a conventional habitat are respectively 28,14 °C and 29,02 °C. The values of storage and restitution heat efficiencies in PCM are respectively 25.92% and 19.34% for the first day.

Keyw ords: Therm al com fort , Heat st orage, Bioclim at ic habit at

I N TROD U CTI ON

Faced with the rising cost of energy and the negative consequences of greenhouse gas emissions, it has become essential to develop innovative systems that permit to improve energy efficiency, better control of energy demand and a reduction in CO 2 emissions (Bedecarrats, 2010).

The housing sector consumes about 40% of primary energy consumption and contributes up to 33% of greenhouse gas emissions. In addition, almost half of the building's energy consumption is devoted to the air conditioners in developed countries [START_REF] Ramakrishnan | Thermal Energy Storage Enhancement of Lightweight Cement Mortars With The Application of Phase Change Materials[END_REF][START_REF]A Review on Buildings Energy Consumption Information[END_REF]. Nevertheless, much of the space conditioning energy is ultimately lost through the building envelope due to the lack of thermal inertia in modern buildings. Thus, energy saving and green building design have become an important research area. A common theme among effective building design strategies is the improvement of the thermal mass of the building envelope for storing thermal energy during the day and restitute to the building during the cooler night. This can reduce maximum temperatures and temperature fluctuations inside the habitat [START_REF] Soares | Review of Passive Pcm Latent Heat Thermal Energy Storage Systems Towards Buildings' Energy Efficiency[END_REF]Baetens and al., 2010).

Bioclimatic architecture, for its part, is a technique whose cost is not too high, which allows to make spectacular savings. A large number of homes built by our ancestors already used this technology: spacious, using quality materials and a good finish. Then the concern for yield and profit has led to the loss of quality, the use of unhealthy and less expensive materials and the reduction of spaces, to obtain what we see today: rabbit cages and villas where it's not good to live (Boukli [START_REF] Hacene | Aspects Energétiques, Economiques et Environnementaux d'une Habitation Ecologique[END_REF]. The application of bioclimatic principles is a factor in reducing energy consumption and CO2 emissions in the construction sector [START_REF] Tzikopoulos | Modeling Energy Efficiency of Bioclimatic Buildings[END_REF][START_REF]Review of Bioclimatic Architecture Strategies for Achieving Thermal Comfort[END_REF].

Heating, ventilation and air conditioning (HVAC) are the largest energy consumers in the building sector (Hoy-Yen and al., 2010).

The design of buildings is based mainly on energy saving criteria [START_REF] Berger | Ambiances Radiatives et Confort Thermique[END_REF]; It associates the insulation, solar contributions, inertia with thermal comfort, because it is defined on the one hand as the satisfaction expressed by an individual with respect to the thermal environment of the medium in which it evolves (Boukli [START_REF] Hacene | Thermal Requirements and Temperatures Evolution in an Ecological House[END_REF].

The feeling of comfort is a difficult concept to quantify. It depends on a large number of criteria aimed at achieving an adequate thermal equilibrium between the human body and its environment [START_REF] Brun | Amélioration du confort d'été dans les bâtiments à ossature par ventilation de l'enveloppe et stockage thermique[END_REF].

Thermal comfort is an obvious need and all built structures (habitable sector) or constructed (transport sector) are equipped with systems to heat or cool the indoor climate (Ahmad, 2004;[START_REF] Candas | Confort Thermique[END_REF]. Numerous physical or psychophysical studies have shown relationships between the intensity of the sensation and that of the stimulus. The sensation increases according to the intensity of the stimulus, or the gap between the stimulus and a threshold.

In the last twenty years, several forms of bulk encapsulated phase change materials have been marketed for active and passive solar systems, including direct gain. However, the surface of 5 most encapsulated commercial products was insufficient to provide heat to the building after the phase change material was melted by direct sunlight. On the other hand, the walls and ceilings of a building offer large areas of passive heat transfer in each zone of the building [START_REF] Neeper | Thermal Dynamics of Wallboard with Latent Heat Storage[END_REF].

The integration of phase-change materials in building walls allows a large amount of energy to be stored by latent heat in a small volume of material [START_REF] Younsi | Etude expérimentale d'un mur trombe-composite à échelle réduite intégrant des matériaux à changement de phase[END_REF]. The special feature for phase change materials is that they can store heat energy by latent heat, as well as by sensible heat [START_REF] Kuznik | Optimization of Phase Change Material Wallboard for Building Use[END_REF]. Thermal energy storage is an attractive technology because it is the most appropriate method to correct the gap between energy demand and supply [START_REF] Oro | Review on Phase Change Materials (PCMs) for Cold Thermal Energy Storage Applications[END_REF][START_REF] Hussein | Review of Development Survey Of Phase Change Material Models In Building Applications[END_REF].

Belén Zalba et al. [START_REF] Zalba | Review on Thermal Energy Storage With Phase Change: Materials, Heat Transfer Analysis and Applications[END_REF] used phase-change materials in free-cooling systems. Free cooling is understood as a means of storing freshness outdoors during the night, to cool the interior of the building during the day. The use of phase change materials is appropriate because of the small difference in temperature between day and night indoors and outdoors.

Mohamed M. Farid et al. [START_REF] Khudhair | A Review on Energy Conservation in Building Applications With Thermal Storage By Latent Heat Using Phase Change Materials[END_REF] have defined that, the latent heat accumulator is one of the most efficient means of storing thermal energy. Unlike the sensible heat storage method, the latent heat storage method provides a much higher storage density, with a lower temperature difference between storage and heat release. Phase change materials can provide a large latent heat storage capacity over the range of temperature typically encountered in buildings, so they can improve the degree of thermal comfort. Lavinia Socaciu et al. [START_REF] Socaciu | Review on Phase Change Materials for Building Applications[END_REF] studied phase change materials in building construction, the study shows that the incorporation of appropriate phase change material into buildings can be an effective solution to reduce energy consumption and cool down the buildings. A small amount of phase change material can lift the thermal inertia substantially without increasing the mass of the structure.

The incorporation of appropriate phasechange materials into walls, ceilings and floors of buildings can directly capture solar energy and store large amounts of thermal energy in the building envelope without the associated large structural mass at sensible heat storage, which helps to reduce the frequency of changes in internal air temperature and keep the temperature close to the desired temperature for a longer period of time [START_REF] Fang | Study on Preparation of Montmorillonite-based Composite Phase Change Materials and Their Applications in Thermal Storage Building Materials[END_REF].

Chi-ming Lai et al. studied the thermal performance of a honeycomb aluminum wall panel incorporating microencapsulated phase change material. The honeycomb phase change material module has sufficient heat conduction enhancement. As a result, internal temperatures within the module were similar and did not show latent heat of stratification. This study applied in addition to honeycomb phase change material, has the ability to regulate the temperature of photovoltaic panels, and the results indicated that photovoltaic panels with nest phase change material module of bees have increased the efficiency of the photovoltaic panel generator compared to normal photovoltaic panels attached directly to the exterior wall of a building (Chi-ming and Shuichi, 2014). 6 Y. Zohir et al. [START_REF] Younsi | Etude expérimentale d'un mur trombe-composite à échelle réduite intégrant des matériaux à changement de phase[END_REF] have experimentally studied a reduced scale composite-wall containing phase-change materials. In this study, they found that the storage of latent heat thermal energy appears indeed very interesting in comparison with sensible heat storage. The main advantages are the storage of a large amount of heat in a reduced volume of phase change materials and the return to a temperature level close to thermal comfort temperatures.

Alvaro Gracia et al. made a study on phase change materials and thermal energy storage for buildings. In this study, they used the method of integrating phase change materials into the racks of clay bricks to increase energy efficiency in buildings. For this study, they discussed the notions of storage by sensible heat, latent heat and thermochemical reaction. Through this study, they have shown that it is possible to do heating and cooling with thermal energy storage in buildings (Gracia and Cabeza, 2015). B. Zivkovic et al. made a study on the isothermal phase analysis of phase change materials contained in rectangular and cylindrical containers. The calculation method for the isothermal phase change of encapsulated phase change material in a single container is based on the enthalpy method. The theoretical model was verified with a test problem and an experiment carried out to evaluate the validity of the hypotheses of the mathematical model.

The results obtained with the proposed method show a very satisfactory agreement with the results obtained from other calculation models. In addition, the comparison between the melting time of the phase change materials contained in the rectangular and cylindrical containers was performed and the results show that the rectangular container requires nearly half of the melting time obtained for the cylindrical container with the same volume and the same zone of heat transfer between the coolant and the wall of the container. It is therefore preferable to use rectangular containers to encapsulate the phase change material rather than to the cylindrical containers [START_REF] Zivkovic | An Analysis of Isothermal Phase Change Of Phase Change Material Within Rectangular and Cylindrical Containers[END_REF].

N U M ERI CAL M ETH OD 1. Climatic Data

The data we used for this research are the climatic data of the typical day of March with an overall radiation of 1000 W/m 2 and a maximum and minimum temperature respectively of 35°C and 25°C. These data allow us to find the hourly variations of the ambient temperature and the solar flux by considering a sinusoidal variation (figure 1).

Figure 1 shows the evolution of horizontal global solar flux and ambient air temperature for the typical day in March. Indeed, March is the hottest time of the year for Guinea. We therefore select the climatic data for the typical March day as input to our program as they allow us to analyze the thermal behavior of the habitat for extreme climatic conditions.

Description of the Physical Model of the Bioclimatic Habitat

The bioclimatic roof model we use is a flat solar collector whose roof is made of transparent tiles. Consider, as shown in Figures 2 and3, a bioclimatic type habitat that can be decomposed into a parallelepiped-shaped enclosure and a roof separated by a rectangular panel containing phase-change materials encapsulated as false ceilings (Figure 4). The roof is likened to a flat wall of rectangular section, inclined at an angle to the horizontal false ceiling, which is extended to a length sufficient to protect the facade of the It is thus composed of an inclined transparent tile plate forming an angle of 30 ° with respect to the horizontal made of a panel of phase-change material and the vertical walls in stabilized earth bricks. Under the roof air circulate in the enclosure by natural convection. The outer face of the roof has transparent tiles. The thermophysical properties of the materials constituting the roof, the false ceiling and the parallelepiped enclosure, assumed to be constant are reported in Table 1.

Mathematical Formulation of the Model

Simplifying Hypotheses

The methodology adopted for describing the thermal behavior of our building model is based on nodal analysis [START_REF] Kemajou | Real Impact of The Thermal Inertia on the Internal Ambient Temperature of the Building in the Hot Humid Climate : Simulation and Experimental Study in The City of Douala in Cameroun[END_REF].

The detailed study of the heat transfer phenomena involved in the functioning of the building (Figures 2, 3 and4) leads us to make a number of hypotheses, the main ones being:

Basic Equations

The establishment of transfer equations is based on the analogy between thermal and electrical transfers.

In a general way, the instantaneous variation of the energy within a component of the habitat is equal to the algebraic sum of the flux densities exchanged within this component.

It is written:

      n i x xij i i i i DFSA T s T Cp M 1  ...(1)
 xij : Heat flux density exchanged by the transfer mode x (conduction, convection and radiation) between the media (i) and (j), (W.m -2 )

S:

Wall section (m 2 )

DFSA i : Solar flux density absorbed by the material (i) (W.m -2 )

DFSA i =  i  i ...( 2
)
 i : Thermal absorption coefficient of the material (i)

 i : Solar flux density captured by the surface of the medium (i) (W.m-2).

By introducing an exchange coefficient h xij and by linearizing transfers, we can write:

 xij = h xij (T j -T i ) ...(3)
Thus, the equation 1 is written:

  i j n i x xij i i i i T T h DFSA T s T Cp M       1 ...(4)
We will then apply Equation 3 to the various medium of our system.

Exterior wall of building ... [START_REF] Hacene | Aspects Energétiques, Economiques et Environnementaux d'une Habitation Ecologique[END_REF] Internal wall of building

  pex pi pex pex pex pex pex pex pex T T Ep K t T s Cp M          ) ( . pex vc pex vc pex amb ex T T hr T T hc      
  pi pe pi pi pi pi pi T T Ep K t T s Cp M       ) ( 1 pi i i pi i n i pi airh i T T F hr T T hc        ...(7)

Air Renewal

Considering a volume of air V exchanged between the outside and the inside of the habitat, this corresponds to a quantity of heat Q in joules. By supposing that volume is exchanged every hour. We will thus have an air exchange flux  ra (J.h -1 ):

  airins airout ra T T Q c     ... (8)
With:

Q. the volumic flow rate in m 3 .h -1 c the heat of the air (c=1225J.m -3 .K -1 ).

T airout ambient air temperature outside the home T airins air temperature inside the enclosure of the building

Heat Equation in the Phase Change Material

Continuity of flux at the solid-liquid interface, by posing that the difference between the flux on either side of the interface is equal to the amount of instantaneous heat released or absorbed in

Boundary Conditions

Boundary conditions are defined along both axis X and Y.

Following the Y axis, we have:

... [START_REF] Neeper | Thermal Dynamics of Wallboard with Latent Heat Storage[END_REF] The development of this equation gives us:

... [START_REF] Oro | Review on Phase Change Materials (PCMs) for Cold Thermal Energy Storage Applications[END_REF] ...( 21)

Following the X axis, we have:

...( 22)

By analogy with equations 20 and 21, we find:

... [START_REF] Ramakrishnan | Thermal Energy Storage Enhancement of Lightweight Cement Mortars With The Application of Phase Change Materials[END_REF] ... [START_REF] Soares | Review of Passive Pcm Latent Heat Thermal Energy Storage Systems Towards Buildings' Energy Efficiency[END_REF] After the discretization of these equations, we find the following dimensionless equations:

Following Y, we have:

... [START_REF] Socaciu | Review on Phase Change Materials for Building Applications[END_REF] ... [START_REF] Tzikopoulos | Modeling Energy Efficiency of Bioclimatic Buildings[END_REF] Following X, we have:

. . . [START_REF] Vasilios | Mathematical Modeling of Melting and Freezing Processes[END_REF] ...( 28)

M ETH OD S

For the resolution of the equations in the building and in the phase change material, we used the Gauss and Thomas algorithms respectively. The enthalpy method is applied for the resolution of the heat equation, it is known by its robustness and results from the fact that it makes it possible to determine the temperature field without having recourse to the knowledge of the progression of the front of solidification in the time.

It thus makes it possible to solve phase change problems for complex geometries but also when there is a pasty zone. The principle of this method is that for both phases, a single variable, enthalpy H, is used as an unknown and thus reduces the equation system to a single heat transport equation (9 to 16). For this method, generally, the reference temperature is equal to the phase change temperature (melting temperature).

RESU LTS AN D D I SCUSSI ON Model Validation

Many validation methods of the model exist, for our case, we validated our model with the experimental study carried out by [START_REF] Zivkovic | An Analysis of Isothermal Phase Change Of Phase Change Material Within Rectangular and Cylindrical Containers[END_REF] who have made an experimental study of the heat transfers in a parallelepipedic tank filled with PCM (Hydrated salt) soaked in a tank water at a temperature of 60 °C. The walls of this tank are made of stainless steel (100x100x20 mm) and thermally insulated (adiabatic process). The temperature in the middle of the PCM is measured using a K-type thermocouple placed in the center of the tank. The curve in Figure 6 below allows us to say that there is indeed a close qualitative and quantitative relationship between our numerical simulation and the chosen experimental model. The fusion time of our numerical model is long compared to the experimental model, this is due to the fact that the phenomenon of convective heat transfer at the level of the PCM is neglected in the numerical model.

PCM STUD I ED

The following results concerning the phasechange material alone have been obtained with the same boundary conditions as the used for validation.

Figure 7 shows the evolution of the phase change for the hydrated salt and paraffin. According to this figure, we find that the hydrated salt has a melting point higher than that of paraffin RT27 with a longer melting time than for paraffin. This indicates that the heat stored in the hydrated salt is important compared to that stored in paraffin RT27. The melting temperatures of these two materials (hydrated salt and Paraffin RT27) are respectively 29.9 °C and 27 °C. Figure 9 shows the evolution of the heat stored at middle point of the phase change material as a function of the thickness after 8000 second. This figure shows that the stored heat increases according to the variation of the thickness from 5mm up to 20 mm, at 25 mm where the heat is constant and starts to decrease from 30 mm up to 50 mm where the stored heat is less important. So there is less heat stored in the center with an increase in thickness because the heat takes more time longer to reach the center.

H ABI TAT STU D Y Temperature Distribution

Figure 10 represents the evolution of the air temperature in the enclosure of the habitat with and without phase-change material. According to this figure, we find that the air temperature in modern habitat (with PCM) decreases compared to conventional habitat (without PCM) because the PCM stores the heat during the day and restitutes it during the night, while in the classic habitat, this heat penetrates inside the habitat and the air is overheated. The temperatures observed at 16 h are respectively for a classical habitat and a modern habitat 29.64 °C and 28.71 °C. The ambient air temperature is low compared to the indoor temperatures in the 6 h -13 h and 23 h -5 h intervals due to the effect of the sky. The ambient air overheats and exceeds the indoor temperature of the room from 14 h to 22 h for the modern habitat. than those of conventional habitat with a difference of about 2 °C. This demonstrates that integrating PCM into habitat provides significant thermal inertia compared to habitat without PCM. The temporal shift between the outdoor and indoor temperatures of the south wall for habitats with and without PCM are respectively 5 h and 4 h.

Figure 12 shows the wall thickness influence on air temperature inside the habitat enclosure with and without PCM. We find that the more the thickness of the wall varies, the greater the thermal inertia and air temperature inside the habitat. By comparing the temperature distributions of the two habitats, we observe that the air temperature decreases in the habitat integrating PCM. This shows that some of the energy provided by the solar flux is stored in the PCM and does not reach the habitat enclosure.

Figure 1 2 : Thick ne ss W a ll I nflue nce on H a bit a t Air Te m pe ra t ure w it h a nd W it hout PCM

We present in Figure 13 the influence of solar flux on the indoor air temperature with and without PCM for a thickness of the wall fixed at 15 cm. We find that for an important flux, the habitat integrating PCM stores a lot of heat during the day and restitutes it during the night into the enclosure. While the habitat without PCM overheats and the inside temperature rises. We can say that the PCM acts as a temperature damper. For a maximum global heat flux of 1300 W.m -2 , the maximum air temperature in the habitat without PCM is about 30°C and that of the habitat with PCM is 28°C. We present in Figure 15 the influence of the solar flux on the distribution of the south interior wall temperature of the habitat with and without PCM for a wall thickness fixed at 15 cm. This wall temperature have the same evolution that the solar heat flux. It increases from 6 h, passes through its maximum towards 13 h and then decreases until 18 h. We notice that the wall temperature is less important for the wall of the habitat with phase-change material it means that the latter stores the heat and so the thermal inertia is more important. For a maximum flux of 1600 W.m-2, the maximum temperature of the inner face of the south wall of habitat without PCM is approximately 30 °C and that of habitat with PCM is 28 °C. Figure 16 shows the influence of air change rate on air temperature in the enclosure of a habitat with and without PCM. On this figure, we can see that an increase of the air change rate induces an increase of the air temperature. This increase is also present during the night period but is less pronounced for the habitat with phasechange material. This is due to the exchange of warm air from the outside to the interior of the habitat during the day. It is there for necessary to limit the rate of air exchange in the habitat if we want to increase the efficiency of the energy storage in the phase-change material.

CON CLU SI ON

We presented a numerical modeling of thermal transfers within a building and heat storage by the enthalpy method in a phase change material. This enabled us to analyze the influence of certain parameters on the thermal behavior of the building. We validated our numerical model with the experimental results proposed by [START_REF] Zivkovic | An Analysis of Isothermal Phase Change Of Phase Change Material Within Rectangular and Cylindrical Containers[END_REF]. For show the utility to use phasechange material, we compared the results obtained for an habitat with and without phasechange material. We have shown the influence of different parameters on the evolution of the temperature inside the enclosure, such as the nature of the phase-change material, the wall thickness and the solar heat flux. The efficiency of storage and discharge of heat and the influence of the variation of certain physical parameters on the distribution of the temperature.

REFEREN CES

N OM EN CLATU RE

h r sol Wall-floor radiation coefficient [W.m -2 .K -1 ] h r vc Wall-celestial vault radiation coefficient [W.m -2 .K -1 ] T vc Celestial vault temperature [K] h c Convection transfer coefficient [W.m -2 .K -1 ] h r Radiation transfer coefficient [W.m -2 .K -1 ] h d Conduction transfer coefficient [W.m -2 .K -1 ] K Thermal conductivity [W.m -1 .K -1 ] h cex Convection transfer coefficient at a outer wall [W.m -2 .K -1 ]

Figure 1 :

 1 Figure 1 : H orizont a l Globa l Sola r Flux a nd Am bie nt Te m pe ra t ure of t he St a nda rd D a y of M a rch

Figure 2 :

 2 Figure 2 : D ia gra m of t he Roof of t he H a bit a t M ode l

Figure 3 :Figure 4 :

 34 Figure 3 : D ia gra m of t he Pa ra lle le pipe d Enclosure of t he Bioclim a t ic H a bit a t

  1. The heat transfers are unidirectional, perpendicular to the walls 2. The materials are assimilated to gray bodies 3. The thermophysical properties of materials are assumed to be constant 4. The rate of air change is imposed in the enclosure of the building 5. The convective heat transfer between the upper wall of the PCM and the interior air of the roof is neglected 6. The upper wall of the PCM receives a quantity of imposed solar flux 7. The heat transfer mode in the PCM is only conduction

Figure 5 :

 5 Figure 5 : Evolut ion of t he Te m pe ra t ure of t he PCM w it h Pha se Cha nge

B

  The systems of algebraic equations obtained by establishing energy balances on the different components of the habitat model are of the form (Kemajou and Mba, 2012): ...(29) with: T(t) Vector of temperatures at different timedependent nodes ()Tt t C Vector column consisting of thermal capacities at different nodes C K Square matrix composed of thermal conductance K Matrix coefficient for the different nodes (t) Vector column representing the inputs of the system (thermal excitations) The equations (5 to 7) are discretized by an implicit finite difference method, whose resolution requires iterative computation. The system is solved by the Gauss algorithm coupled to an iterative method because the transfer coefficients depend on the temperatures of the different building components.

Figure 6 :

 6 Figure 6 : M ode l Va lida t ion

Figure 7 :

 7 Figure 7 : Com pa rison Be t w e e n Tw o-pha se Cha nge M a t e ria ls

Figure 8

 8 Figure8represents the evolution of the heat stored in the phase-change material at the various points chosen for our study over time. On this figure, we find that the extreme points retain much more heat compared to the intermediate points. This is because the height of 16 mm receives the heat first, before returning it to the height of 2 mm which is the last point near the interior of the habitat.

Figure 8 :

 8 Figure 8 : Evolut ion of t he H e a t St ore d in t he PCM a t D iffe re nt Point s Ove r Tim e

Figure 11

 11 Figure 11 illustrates the temperature profiles of the exterior and interior faces of the south wall for habitat with and without PCM. Temperatures observed at the exterior face of the southern wall of conventional habitat dominate those of modern habitat with a difference of 6 °C between maximum temperatures. Temperatures of interior face of southern wall for modern habitat are lower

Figure 9 :

 9 Figure 9 : H e a t St ore d in t he PCM According t o t he Va ria t ion of t he Thick ne ss

Figure 1 0 :

 0 Figure 1 0 : Air Te m pe ra t ure Profile in t he H a bit a t w it h a nd W it hout PCM

Figure 1 3 :

 3 Figure 1 3 : Flux I nflue nce on I ndoor Air Te m pe ra t ure w it h a nd W it hout PCM

Figure 1 4 :Figure 1 5 :

 45 Figure 1 4 : Flux I nflue nce on I ndoor Roof Air Te m pe ra t ure of H a bit a t w it h a nd W it hout PCM

Figure 1 6 :

 6 Figure 1 6 : Air Cha nge Ra t e I nflue nce on I ndoor Air Te m pe ra t ure w it h a nd W it hout PCM

Figure 1 7 :

 7 Figure 1 7 : Sola r Flux I nflue nce on t he St ora ge a nd D ischa rge Efficie ncy

  h ripi Radiation transfer coefficient between two walls [W.m -2 .K -1 ]  te x Absorption coefficient of the outer cover  te x Heat flux density at the outer wall [W.m -2 ] Heat flux exchange by air renewal [J.h -1 ] H Enthalpy [J] M tex Mass of the outdoor roof [Kg] C p tex Thermal capacity of the outdoor roof [J.Kg -1 .K -1 ] T ti Temperature of the internal face of the roof [°C] T tex Temperature of the external face of the roof [°C] h r vc,pex Radiation transfer coefficient between the sky and the outer wall [W.m -2 .K -1 ] h r flor,pex Radiation transfer coefficient between the floor and outer wall [W.m -2 .K -1 ] h ci,pi Convection transfer coefficient between the air and the inner wall [W.m -2 .K -1 ] h ripi Radiation transfer coefficient between two internal walls [W.m -2 .K -1 ] T airh Air temperature inside the habitat [°C]

  temperature inside the habitat without PCM [°C] T AIPCM Air temperature inside the habitat with PCM [°C] T AMB Ambient temperature [K] T SEH Outside south wall temperature of habitat without PCM [°C] T SEPCM Outside south wall temperature of habitat with PCM [°C] T SIH Interior south wall temperature of habitat without PCM [°C] T SIPCM Interior south wall temperature of habitat with PCM [°C] T AIH10 Air temperature inside the home without PCM for wall thickness set at 10 cm [°C] T AIPCM15 Air temperature inside the home with PCM for wall thickness set at 15 cm [°C] T AFX7 Air temperature inside the home without PCM for a maximum flux of 700 W / m2 [°C] T AMFX10 Air temperature inside the home with PCM for a maximum flux of 1000 W/m 2 [°C] T ATFX7 Air temperature inside the roof of the home without PCM for a maximum flux of 700 W / m 2 [°C] T ATMFX16 Air temperature inside the roof of the home with PCM for a maximum flux of 1600 W / m 2 [°C] T SIFX7 Interior south wall temperature of the home without PCM for a maximum flux of 700 W / m 2 [°C] T SIMFX13 Interior south wall temperature of the home with PCM for a maximum flux of 1300 W / m 2 [°C] T AHTAU3 Air temperature inside the home without PCM for an air change rate set at 3 / 24h [°C] T AMTAU9 Air temperature inside the home with PCM for an air change rate set at 9 / 24h [°C] H 16 Enthalpy at the height of 16 mm [J] H 2 Enthalpy at the height of 2 mm [J] L PCM length[m]  ref Reference thermal conductivity of the PCM [W.m -1 .K -1 ]  s Thermal conductivity of the solid PCM [W.m -1 .K -1 ] T ref Reference temperature [K] T pcm(IN,JJ) Temperature at the point IN and JJ [K]  Convection coefficient [W.m -2 .K -1 ] X and Y andSpace steps [m]


K Pantzas

G Patriarche

David Troadec

S Gautier

T Moudakir

S Suresh

L Largeau

O Mauguin

P L Voss

A Ougazzaden

Nanometer

Keywords:

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Indium gallium nitride (InGaN) is used in commercial light-emitting diodes (LEDs), and lasers [START_REF] Nakamura | The Blue Laser Diode: The Complete Story[END_REF], and is being actively investigated for use in solar cells [2][3][4][5][6][7][8]. It is one of many scientifically and economically important materials [9][START_REF] Karpov | Modeling of III-nitride light-emitting diodes: progress, problems[END_REF][11] whose quality can potentially be improved if the material were better understood. Detailed chemical mappings of the alloy composition from transmission electron microscopy (TEM) and scanning transmission microscopy (STEM) images can be particularly interesting in such studies, providing useful insights into material growth.

While both TEM and high-angle annular dark field scanning transmission microscopy (HAADF-STEM) give qualitative chemical information, precise, quantitative maps cannot be readily obtained. Indeed, only relative intensity variations in TEM and HAADF-STEM reflect variations in composition. The difficulty lies in the fact that these relative intensity variations cannot be quantitatively interpreted directly from the TEM or STEM image [12]. In order to do 1 Georgia Institute of Technology, GT-Lorraine, 2-3 Rue Marconi, F-57070, Metz, France 2 UMI 2958, Georgia Tech-CNRS, 2-3 Rue Marconi, F-57070, Metz, France 3 CNRS-Laboratoire de Photonique et des Nanostructures (LPN), Route de Nozay, F-91460, Marcoussis, France 4 CNRS, Institut d'Electronique, de Microélectronique et de Nanotechnologie, UMR 8520, F-59652 Villeneuve Dascq, France 5 Laboratoire Matériaux Optiques, Photonique et Systèmes (LMOPS), EA4423, Université Paul Verlaine Metz-Supelec, 2 Rue Edouard Belin, F-57070, Metz, France that, one needs a means of establishing an absolute scale of composition.

Recently, Rosenauer et al [13] implemented a method, originally proposed by LeBeau and Stemmer [14], to obtain quantitative chemical maps from HAADF-STEM by fitting frozen lattice simulations to experimental images. In this paper we report an alternative for obtaining quantitative chemical maps from HAADF-STEM. In our case the STEM intensity variations are quantified experimentally, using a small set of energy dispersive x-ray spectroscopy (EDX) measurements as a chemical calibration of the observed intensity variations. This approach combines the chemical precision of EDX with the spatial resolution of HAADF-STEM, to obtain quantitative chemical mappings that cannot be obtained independently from either technique. It presents two advantages over direct EDX mappings: for the EDX signal to be reliably quantifiable, exposure times of more 60 s are required, making 1024 × 1024 mappings not only impractical but also subject to error due to small displacements of the sample during the acquisition time; moreover, the resolution of EDX is limited by beam-broadening to, at best, 1-2 nm.

In what follows, we first use elastic scattering theory to show how, within certain limits, the scattering law of electrons in HAADF-STEM can be evaluated experimentally through a small set of EDX measurements. We proceed to describe how to compute this scattering law in practice. Then we demonstrate how to use the resulting function to transform qualitative Z-contrast images of InGaN layers into precise, quantitative maps that reflect the variations of indium composition in these layers with nanometer-scale precision. Using the results from such maps we provide useful insights into the impact of fluctuations in the indium composition on the growth of epitaxial InGaN layers.

Experiment

The InGaN epilayers used for this study were all grown on commercial gallium nitride (GaN)/sapphire templates by metal-organic vapor-phase epitaxy (MOVPE), using nitrogen as the carrier gas and tri-methyl indium (TMIn), tri-methyl gallium (TMG) and ammonia (NH 3 ) as precursors for elementary indium, gallium and nitrogen, respectively. The growth temperature was 800 • C. Two samples among them are discussed below. Sample A was 70 nm thick, while sample B was 140 nm thick. The TMIn/III, the ratio of TMIn to the sum of TMIn and TMG, for sample A was 15%, while that of sample B was 25%.

X-ray diffraction (XRD) measurements were performed in a Panalytical X'Pert Pro diffractometer in a triple axis configuration. In order to perform accurate composition and strain measurements, reciprocal space mappings (RSMs) were taken along both the symmetric 0 0 0 6 and asymmetric 1 1 2 4 reflections. The effect of tilt was taken into account by averaging the results from measurements along all six asymmetric ϕ reflections of the wurtzite lattice, as suggested by Moram and Vickers [15]. The composition was then obtained by numerically solving the third-order polynomial linking the composition to the a and c parameters of wurtzite layers for any strain layers, introduced by Schuster et al [16].

Specimens from the samples were then prepared for STEM using focused ion beam (FIB) thinning and ion milling. The specimens were 80 nm thick. In order to preserve the sample surface during FIB preparation, a surface coating consisting of a 50 nm-thick layer of carbon, followed by 100 nm of silicon nitride (Si 3 N 4 ) was applied. Prior to STEM imaging, the specimens were cleaned using an argon plasma cleaner.

STEM and EDX were then performed in an aberrationcorrected JEOL 2200FS microscope, operating at 200 kV with a probe current of 150 pA, and a probe size of 0.12 nm at the full width at half maximum (FWHM). The convergence half-angle of the probe was 30 mrad and the detection inner and outer half-angles for the HAADF-STEM images were 100 and 170 mrad, respectively. The samples were imaged along the 1 1 2 0 zone axis.

Compared to TEM, which is known to cause image artifacts due to electron beam damage [17], currents of at most a few hundred pA are used in STEM and are expected to leave no residual damage on the specimen. Although the stability of InGaN in HAADF-STEM in particular has been demonstrated previously [13], we ensured that this was the case here also by taking several sets of images and comparing them. No change was observed in the sample between the different sets of measurements.

Quantitative measurements of the indium composition from EDX were obtained from the intensity ratio of the Lα line of indium (3.290 keV) to the Kα line of gallium (9.770 keV). The K line of elementary nitrogen (0.392 keV) was also taken into account and revealed that the alloy is stoichiometric. The acquisition time for each EDX spectrum was 60 s, during which no drift in the position of the electron beam was observed. The EDX spectra were acquired using a JEOL 2300D detector and the accompanying JEOL software. The k-factors used by the software had been previously re-calibrated using GaAs, InP, GaN, AlN, GaP and GaSb, as well as the ternary alloys In 0.48 Al 0.52 As and In 0.53 Ga 0.47 As. These alloys are lattice-matched to InP, allowing the composition to be precisely determined through XRD.

Quantification of the Z-contrast

In this section we explain how the quantification works. We first give a proof of principle from first-principles calculations. We then proceed to present an algorithm to implement it. Finally, we use the algorithm to quantify a HAADF-STEM image of sample A. This example is used to evaluate the spatial resolution and chemical precision of the algorithm.

Proof of principle

It has been shown that in HAADF-STEM, at detection inner half-angles greater than 60 mrad and for STEM specimens thinner than 100 nm the major contribution to the collected intensity is that of elastically scattered electrons [18]. In particular it has been shown that a good approximation of the fraction of the incident electron beam intensity scattered towards those angles by an atomic column of the STEM specimen is given by [19]

I n ∝ d Z α . (1)
Here, d is the length of the imaged atomic column (i.e. the STEM specimen thickness), Z is its average atomic number and α is an exponent that will be discussed in more detail later.

For an atomic column of In x Ga 1-x N the average atomic number is given by

Z InGaN = xZ In + (1 -x)Z Ga + Z N , (2)
or, equivalently,

Z InGaN = xZ InN + (1 -x)Z GaN , (3)
where

Z III-N = Z III + Z N .
Thus, one can relate the root α of I n to the concentration x of an atomic column of In x Ga 1-x N by the following equation:

α I n = α √ Kd(xZ InN + (1 -x)Z GaN ). (4)
K is a proportionality coefficient that depends on the imaging geometry and which is constant during the acquisition of the HAADF-STEM image. If one divides both sides of equation ( 4) by the reference intensity

I ref = KdZ α GaN , one finds an intensity ratio R that is linked to the composition by α √ R = (ζ -1)x + 1, (5)
where ζ = Z InN Z GaN is the ratio between the atomic numbers of pure indium nitride (InN) and pure GaN.

From equation ( 5) it is obvious that if one were to know the value of α, one could readily invert the equation and find

x = ξ( α √ R -1), (6)
where ξ = 1 ζ -1 is introduced to simplify the notation. Then one could proceed to compute a chemical composition map by applying equation ( 6) to the HAADF-STEM image.

In elastic scattering theory, the exponent α is equal to 2. However, it has been argued that a more appropriate modeling of Rutherford scattering would need to take into account electron screening. In such cases, α is given by

α = 2 -σ, (7)
where σ ≥ 0 is a factor modeling electron screening. This factor depends both on the scattering atom and the penetration depth of the incident electrons into the atom's electron cloud, which in turn depends on the beam acceleration voltage. As a result, the value of α depends on the material being imaged as well as the imaging conditions. Typical values of α range between 1.7 and 2 [19].

Although one cannot have an a priori knowledge of the value of α, one could estimate it using a set of EDX measurements of the composition, {x i } M i=1 , measured at the same locations as a set of intensity ratios {R i } M i=1 , and then using the following estimator α:

α = log R log ((ζ -1)x + 1) M i=1 . (8)
In section 3.2 we will discuss an algorithm that computes the estimator α and then applies equation ( 6) to obtain a concentration map.

Implementation

Here we will discuss how to practically implement the method outlined in section 3.1. The implementation proposed here was carried out in MATLAB, where the HAADF-STEM image is represented by a N × N matrix, of typical size N = 1024. For clarity, we will first introduce a few conventions on notation before giving the details of the implementation.

Notation.

As discussed above, each pixel of the STEM image corresponds to an intensity I collected at position (k, l) of the STEM probe. In what follows, the intensity I k,l , collected at position (k, l), will be referred to as I, whereas I will refer to the matrix of intensities [I] (k,l) . Furthermore, an N × N matrix of ones will referred to by 1 N .

Point-wise, or element-to-element, matrix division will be denoted by './', i.e. C = A./B will signify:

C k,l = A k,l /B k,l
for all couples (k, l) such that k = 1, . . . , N, l = 1, . . . , N.

(

Finally, a subset of M intensities I i , or a set of M EDX measurements x i will be denoted by {I i } M i=1 and {x i } M i=1 , respectively.

Algorithm.

A schematic diagram of the proposed algorithm is given in figure 1. In HAADF-STEM, the intensity is collected by an annular scintillation detector combined with a photomultiplier type (PMT), resulting in the useful signal being added on top of a background intensity from the polarization current of the detector. Thus, the first step consists of measuring the average background signal I b and subtracting it from the collected image intensity matrix I raw .

We then obtain a normalized intensity image I n that can be used for the quantification

I n = I raw -I b • 1 N . (10)
In the second step, we evaluate the reference intensity I ref defined previously. The reference is taken in the GaN buffer. Since specimen preparation by FIB results in a wedge shaped sample, shown in figure 2, d is a function of the position (k, l) in the image. In practice, one observes an intensity gradient, all the more evident at low magnifications. To account for this, I ref is obtained through a linear least-squares fit of the GaN part of I n , to a first-order, two-dimensional polynomial function. This function is then extrapolated to the whole image, and in particular to InGaN. Since InGaN is not too dissimilar from GaN and the gradient created from FIB is not too pronounced, we expect any residual thickness variations to be significantly smaller than the effect of the wedge.

A point-wise division between the background-free I n and the reference intensity I ref n results in the matrix of intensity ratios R: Then, a set {R i , x i } N i=1 of intensity ratios and the associated composition, independently evaluated by EDX measurements, is used to compute the estimator α given by equation (8). Since the EDX probe interacts with a volume of the specimen larger than the STEM probe, it would be inaccurate to use the intensity of the center of the probe as the value R i used in the estimation. Instead, an average ratio was computed for each position of the probe, obtained by using a moving average filter on the image. The mask for the moving average filter was obtained by modeling the EDX probe interaction volume by a truncated cone. The diameter was computed as the average between the entry and exit disks. This typically results in a 5 × 5 matrix for the mask. Finally, the function given in equation ( 6), combined with the resulting estimator α, is applied to the matrix of intensity ratios R to compute the quantitative composition map.

R = I n • /I ref n . (11)

Example

In this section we show an example of how the algorithm works on sample A. The spatial resolution and the chemical precision are estimated and limiting factors for both are discussed. Figure 3 Figure 4 shows an image of the intensity ratio R, obtained in stage three of the algorithm after subtraction of the background intensity and normalization by the reference intensity. The relative intensity variation between GaN and InGaN is preserved, proving that stages 1 and 2 of the algorithm induce no loss of information. A residual modulation of the intensity parallel to the GaN/InGaN interface is still apparent (highlights of the modulation are pointed out by black arrows in figure 4). Profiles taken 10 nm on either side of the interface are shown in the inset of figure 4. They reveal that this residual modulation is present in both GaN and InGaN and is related to the FIB preparation of the specimen. The impact of the modulation on the precision of the algorithm will be discussed at the end of this section.

In order to evaluate α, a set of 18 EDX measurements was performed. The distribution of the EDX measurements across the HAADF-STEM image is shown in figure 4 (red and blue dotted lines). As indicated in section 3.2.2, the corresponding set of intensity ratios {R i } M i=1 was computed by using a moving average filter with a 5 pixel × 5 pixel mask. This is the required mask size for a 2 nm EDX disk at the magnification used in for the HAADF-STEM image presented in figure 3. It should be noted that this averaging step does not affect resolution of the chemical mappings, since, at this stage, we are estimating the scattering law at places were the composition is known.

Equation ( 8) was used on the set {R i , x i } M i=1 , shown in table 1. A value of α = 2±0.1 was found to be the appropriate power model in this case. Similar values for the exponent α in the case of InGaN alloys were reported by Amari et al [20], using the exact 'sum-of-squares' scattering model:

R = xZ α In + (1 -x)Z α Ga + Z α N Z α Ga + Z α N . ( 12
)
To check the validity of our approach, we also evaluated the fit of our values to this model. A good fit could not be obtained unless we let α = 1.7. This calculation also revealed that, for indium compositions between 0 and 80%, the difference between our model and the sum-of-squares is negligible. Further Table 1. Summary of the EDX measurements and their corresponding intensity ratios, used for the map presented in figure 5. Numbers in italic were excluded from the calculation, since they were located in regions where FIB artifacts were suspected to be present. Figure 5 shows the concentration map obtained after application of equation ( 6). The colorbar indicates the computed indium composition. It has been extended below zero, to better show that the estimation error is randomly distributed around 0 in the GaN layer. The root mean square error (RMSE) between the measured and computed values of the composition was around ±0.5%, of the order of the chemical precision of EDX. The resolution of the map is that of the original HAADF-STEM image, since every intensity ratio in the original image was transformed using equation (6).

However, when interpreting the calculated composition, one needs to take into account artifacts introduced by small deviations from the planar thickness model used previously to compute I ref . Indeed, one can distinguish two main sources of artifacts, both related to the FIB preparation of the sample: residual large-scale thickness modulations, such as the ones pointed out by the black arrows in figure 4; and uncertainty due to the amorphous layer on the sidewalls of the specimen, created during the FIB etching. This second source of uncertainty induces an error in both the EDX measurement of the composition and a short-length scale error in the thickness model [21]. In what follows, we will attempt to estimate these errors. All errors are given in atomic %.

The short-scale error in thickness can be evaluated in the GaN part of the concentration map. Indeed, it is responsible for the specularity observed both in the HAADF-STEM image in figure 4 and in the concentration map in figure 5. The standard deviation of the concentration evaluated in the GaN part of the concentration map yields an average error of ±0.5%. Finally, the error due to large-scale thickness modulations can be evaluated in the composition map figure 5 by estimating the standard deviation of the composition along a line at the same position as the red dashed line of figure 4. The resulting error is computed at ±1%.

Thus, any variations in composition that are less than ±1.5% cannot be detected in this case. A similar precision was reported by Rosenauer et al in [13]. Further improving the sample preparation may contribute to removing these artifacts and increasing the precision of the map. For the purpose of the application presented in section 4, however, this value was sufficient.

Finally, the results obtained by the algorithm were compared to the concentration deduced from XRD RSMs. Figure 6 shows an average, tilt-free RSM for the asymmetric 1 1 2 4 (left) and the symmetric 0 0 0 6 (right) reflection, respectively. Table 2 compares the indium concentration deduced from XRD, EDX and the concentration map. It shows that the algorithm achieves similar results to what can be found by other techniques. Thus, the chemical mappings obtained can be used reliably in the analysis of InGaN layers.

Application to InGaN layers presenting double XRD diffraction peaks

In this section, the algorithm is applied to the study of thick InGaN layers presenting double diffraction peaks in XRD RSMs. Sample B is an example of such a case, where two InGaN diffraction peaks besides the intense GaN peak are apparent in both the asymmetric and the symmetric RSM (cf figure 7). The first diffraction spot corresponds to pseudomorphically strained 13% InGaN. This diffraction spot elongates towards its fully relaxed position, indicating an elastic relaxation of strain. The second diffraction spot corresponds to a fully relaxed 20% InGaN layer. Similar results have been observed for thick InGaN layers in [22][23][24][25][26][27][28][29].

A HAADF-STEM image along the 1 1 2 0 zone axis of sample B is shown in figure 8. It is similar to what has been reported in [26][27][28]. The GaN, InGaN and C/Si 3 N 4 cap (3D) growth mode, in agreement with results reported in [24,28,29].

If one inspects the InGaN layer closely, one can see the presence of two distinct regions: a region of homogeneous contrast near the GaN/InGaN interface, labeled InGaN 1, and a second region, where the contrast fluctuates laterally, concentrated at the top of the InGaN layer, labeled InGaN 2. There are two types of fluctuations in InGaN 2: short-scale variations, of the order of a few nanometers, which give InGaN 2 a 'veined' aspect, and long-scale fluctuations, of the order of a few tens of nanometers, concentrated essentially near the interface of the two InGaN sublayers. Figure 9 shows the corresponding chemical composition map. The boundary between InGaN 1 and InGaN 2 is clearly identified in the composition map. The average concentration of InGaN 1 computed by the map is 14%, while that of InGaN 2 is 23%. Thus, InGaN 1 and InGaN 2 can be unambiguously identified as the origin of the fully strained and fully relaxed diffraction spots in the RSM, respectively. It should be noted here that, when computing the average composition of InGaN 2, the region within 10 nm of the boundary between InGaN 2 and the carbon layer was excluded. Indeed, due to the pronounced roughness of the InGaN 2 sublayer, there are strong deviations from the first-order thickness model used in the algorithm, which result in erroneous composition estimates.

The map also reveals that both short-scale and long-scale lateral fluctuations in the HAADF intensity in InGaN 2 correspond to variations in the indium concentration of the order ±5%. This is significantly higher than the uncertainty due to FIB artifacts and also corresponds to the spread in composition determined using XRD. Due to the rapid variations in the HAADF intensity, it has been argued that the short-scale fluctuations in intensity may be linked to the presence of threading dislocations [27]. In order to check whether the algorithm does not misinterpret these fluctuations, atomic scale resolution images of the pyramids, where the threading dislocations are expected to emerge, were taken to verify if this is the case. An example of the results is given in figure 10(a). It shows an atomic scale resolution HAADF-STEM image of the tip outlined by a white dashed rectangle in figure 8. This particular pyramid was chosen because it is entirely within the FIB section and can be observed all the way to the tip. Moreover, it presents very strong short-scale fluctuations that can be easily misinterpreted. While a line of lower contrast is shown to run through the middle of the pyramid all the way to the top, no disruption of the atomic planes is observed. This is confirmed by the corresponding bright field (BF)-STEM image, shown in figure 10(b). Additional EDX measurements along the dotted red line shown in figure 10(a) (see inset) revealed a local decrease from 22% at the edges of the pyramid to 14% near the middle. We can thus conclude that the map accurately reflects fluctuations of the indium composition.

As discussed in [29], we believe the observed process, i.e. the elastic relaxation of InGaN 1, the transition to 3D growth, the increase in the average indium concentration in InGaN 2 and the lateral fluctuations of indium concentration in InGaN 2, to be the result of an accumulation of excess indium at the surface during early stages of growth. This in turn leads to a lateral modulation of the incorporation rate of indium into the lattice, which results in a lateral modulation of the growth rate. This eventually leads to the elastic relaxation of the initial InGaN, the appearance of additional facets favorable to growth and the transition to 3D growth. The strong faceting and eventual 3D growth are evidenced by the roughness of InGaN 2. No particular family of planes is favored in the faceting, as revealed by the presence of pyramids with different basal angles and heights. All pyramids, however, present lateral fluctuations of the composition to some extent. Finally, the stark difference between the average compositions of InGaN 1 and InGaN 2, as well as the significant lateral variations of indium composition in InGaN 2, can be explained by either, or both, of the two following mechanisms:

(i) The absence of compositional pulling [30,31] in the elastically relaxed InGaN 2 sublayer. (ii) A more favorable binding of indium atoms along the prismatic 1 0 1 l planes [32,[START_REF] Wunderer | Optimization of Semipolar GaInN/GaN blue/green Light Emitting Diode Structures on 1 101 GaN side facets[END_REF] that make up the facets of the pyramids.

Conclusion

In this paper we have shown that is possible to experimentally evaluate the scattering law of electrons in HAADF-STEM using EDX. In particular, we have evaluated it for InGaN epitaxial layers and used the resulting function to compute quantitative chemical maps of the indium composition. The resulting maps are shown to have nanometric resolution, are precise to 1.5% of the composition and are consistent with results from XRD. Based on these mappings, we have shown that fluctuations in the indium composition are responsible for the spontaneous formation of sublayers with distinct compositions and strain states during the growth of InGaN. Finally, several possible explanations for the origin of these fluctuations are discussed.

Figure 1 .

 1 Figure 1. Schematic diagram of the algorithm.

Figure 2 .

 2 Figure 2. Schematic cross-section of the STEM specimen prepared by FIB. The thickness variation has been exaggerated for visual clarity.

Figure 3 .

 3 Figure 3. HAADF-STEM image of sample A taken along the 1 1 2 0 zone axis. The inset shows the intensity profile from the top left to the bottom right corner of the image. The dashed blue line represents the linearly extrapolated I ref .

  shows an HAADF-STEM image of sample A. One clearly distinguishes the GaN template (bottom right), the GaN layer and the carbon-Si 3 N 4 cap layer. A close inspection of the GaN layer reveals a slight intensity gradient, roughly perpendicular to the GaN/InGaN interface and towards the bottom right of the image. The gradient becomes evident if one extracts the intensity profile from the top left corner of the image to the bottom right corner, shown in the inset of figure 3. As discussed in section 3.2.2, this gradient can be attributed to the gradient in the sample thickness. The dashed blue line represents the section of the extrapolated I ref along the profile line. It shows that a first-order polynomial is a reasonable approximation of I ref in the area of interest, i.e. the GaN and InGaN layers.

Figure 4 .

 4 Figure 4. HAADF-STEM image of sample A after subtraction of the background intensity I B and division by the reference intensity I ref . The inset shows profiles of the intensity ratio R along the dashed red and green lines, respectively. Fits in black are given as a guide to the eye. One can see similar oscillations in the intensity along the two lines, proof of a residual lateral modulation of thickness, associated with the FIB preparation.

Figure 5 .

 5 Figure 5. Concentration map of sample A obtained after the application of the algorithm; the colorbar indicates the calculated indium concentration in per cent.

Figure 6 .

 6 Figure 6. XRD reciprocal space mappings of sample A along the asymmetric 1 1 2 4 (left) and the symmetric 0 0 0 6 reflection (right), respectively.

Figure 7 .Figure 8 .

 78 Figure 7. XRD reciprocal space mappings of sample B along the asymmetric 1 1 2 4 (left) and the symmetric 0 0 0 6 reflection (right), respectively.

Figure 9 .

 9 Figure 9. Concentration map obtained of the HAADF presented in figure 8 after the application of the algorithm; the colorbar indicates the calculated indium concentration in per cent.

Figure 10 .

 10 Figure 10. Atomic scale resolution [HAADF (left) and BF (right)] of one of the pyramidal tips. They show no discontinuity in the atomic planes all the way to the top. Thus, the decrease in the HAADF intensity in the middle of the pyramid can only be attributed to a local decrease of the indium concentration in the middle of the pyramid. This is confirmed by EDX measurements (see the inset), revealing a drop in the indium composition along the middle of the pyramid. The red dots indicate the positions of the EDX measurements.

Table 2 .

 2 Average indium composition and the associated standard deviation for both samples, for three different techniques.

			XRD		EDX		EQHAADF
			Average	Standard	Average	Standard	Average	Standard
			composition	deviation	composition	deviation	composition	deviation
			%	%	%	%	%	%
	Sample A		12	1	12.7	0.9	14	1.4
	Sample B	InGaN 1 InGaN 2	13.5 22	1 5	13 23	1 6	14 23	1.4 5

Acknowledgment

This study has been funded by the ANR Habisol 2009, project NewPVonGlass.


Mario Barbatti
email: mario.barbatti@univ-amu.fr

Defining the Temperature of an Isolated Molecule

The microcanonical temperature of an isolated molecule is derived in terms of Boltzmann and Gibbs volume entropies within the quantum harmonic vibrational and equivalent degenerated model approximations. The effects of the entropy functional choice and various approximations are examined. The difference between Boltzmann and Gibbs volume temperatures is negligible for molecules bigger than ten atoms. However, it is significant for smaller systems, opening a way to probe them experimentally. A simple, analytical expression of the temperature as a function of the vibrational energy is provided, allowing predictions with a 3% margin of error compared to the exact harmonic estimate. The microcanonical temperature is discussed and exemplified with polycyclic aromatic hydrocarbon molecules and other molecules of astrophysical interest.

I. INTRODUCTION

Any excess of vibrational energy in a molecule above its zero-point level may induce physical-chemical reactions such as isomerization, dissociation, or rearrangement. The role of such vibrational energy excess is particularly relevant when a light-excited molecule transfers the photon energy to the ground state without emitting radiation. The classic example is the vision mechanism, where a complex sequence of events is triggered in the rhodopsin protein after the retinal chromophore converts the photon energy into vibrational energy. 1 We often say that a molecule with large vibrational energy is hot. However, if we ask what the molecule's temperature is, the question sounds odd. The fundamental thermochemistry of a canonical harmonic system tells that the environmental temperature determines the vibrational excess through

( ) / , 1 k B e k e ZP h k T k h E T E e   =+ -  ( 1
)
where EZP is the vibrational zero-point energy, Te is the environmental temperature, kB is the Boltzmann constant, h is the Planck constant, and the sum runs over the vibrational degrees k with harmonic frequency k.

Nevertheless, an isolated molecule cannot count on the environment to thermalize. The vibrational energy excess will be trapped there until it is irradiated or induces a chemical reaction.

A hot, isolated molecule is a microcanonical system with constant total energy, and, different from the canonical problem, we want to determine ( ) m TE , the microcanonical temperature as a function of the energy.

Defining the temperature of such an isolated molecule could be worth doing, 2 especially if we aim at describing large systems. Indeed, the temperature of isolated molecules and small finite systems has been under scrutiny for decades in diverse fields. It is needed for computing unimolecular reaction rates, [2][3][4] sampling initial conditions for dynamics, 5 studying fragmentation and radiative cooling of electrostatically trapped ions, [6][7][8] determining thermionic electron emission of clusters beams, 9 characterizing phase transitions in clusters, 10 or assigning molecular species in the interstellar medium. [START_REF] Hendecourt | Proceedings of the 135th Symposium of the International Astronomical Union[END_REF][START_REF] Leger | [END_REF][13] Defining such a temperature may be helpful even for understanding the photophysics of chromophores in vacuum-like hydrophobic cavities of proteins. 14 The microcanonical temperature can be used similarly to the canonical temperature to predict the vibrational population distribution when the molecule has an energy excess above its zero-point level. We will naturally use the microcanonical temperature for isolated molecules in a vacuum. We may also use it for assessing situations where the reaction time (for luminescence, isomerization, dissociation, ionization) is much shorter than the collisional time but longer than the vibrational equilibration time. For example, a photoexcited molecule equilibrates vibrationally in the first excited state within a few picoseconds. If it is in a vacuum or within a collisionless regime, its vibrational states are populated according to the microcanonical temperature corresponding to the energy excess. Thus, if we want to know whether the molecule will fluoresce before it does an intersystem crossing, we should consider that these processes' rates are functions of the microcanonical temperature (as opposed to rates computed at 0 K or at the canonical temperature).

The temperature of isolated systems has been addressed pragmatically using phenomenological approaches for applied purposes. Berry and Smirnov, 10 for instance, started from a reasonable estimate for the heat capacity of isolated clusters to derive their temperature. Another common approach has been to derive the canonical temperature for a canonical ensemble with kinetic averages equivalent to those of the microcanonical system. 3 Andersen et al. 15 followed this way to propose a high-energy approximation for the molecular temperature.

D'Hendecourt et al. [START_REF] Hendecourt | Proceedings of the 135th Symposium of the International Astronomical Union[END_REF] estimated the properties of molecules in the interstellar medium from the observed temperature, also using an equivalent-canonical approach.

To the best of my knowledge, none of these phenomenological-based works rest on a solid statistical basis and are valid over a broad energy domain.

On the other hand, more statistical-mechanics-oriented works have provided much knowledge on small, isolated systems. We know, for instance, that the microcanonical temperature can be defined independently of the probe, 16 the type of fluctuation we may expect, including the dependence on different entropy definitions, [17][18][19][20] and how to count microstates of finite coupled systems. [21][22][23][24][25] Nevertheless, these works usually do not allow for an everyday estimate of the temperature in applied cases of interest.

This paper aims to bridge this gap between theory and applications by deriving a microcanonical temperature for an isolated molecule that can be routinely used in diverse applied fields, has an extended validity domain, and is based on solid statistical mechanics grounds.

II. THEORETICAL MODEL

To define the microcanonical temperature of an isolated molecule, we must go back to the principles of statistical thermodynamics. For a system with total energy E, the temperature is defined as

1 , S T E -   =    (2)
where ( ) SE is the entropy, which is a function of the number of microstates at energy E. Later, we will discuss what functional of E is the most adequate when dealing with a small system far from the thermodynamical equilibrium.

If we approximate our isolated molecule by an ensemble of N quantum harmonic oscillators, the total molecular energy is

, ZP k k k E E n h   =+  ( 3
)
where nk is the number of quanta deposited in vibrational mode k, and  denotes a particular ensemble  

where k k Mn =  , is a very well-known statistical- mechanics textbook example. [START_REF] Callen | Thermodynamics and an Introduction to Thermostatistics[END_REF]27

A. Modeling an isolated molecule

For modeling the molecule, we assume (1) that it is rigid enough, so its vibrational modes are harmonic; (2) there are no external fields; and (3) it populates only a single electronic state.

The first assumption implies that we will not describe anharmonic modes, like intramolecular hydrogen bonds or internal rotations (like those methyl groups are prone to).

We can consider an isolated molecule as an ensemble of weakly coupled quantum harmonic oscillators. By "weakly coupled," I mean that energy can flow between vibrational modes, but the vibrational energy is still approximately given by Eq. (3). This energy flow is crucial for the molecule to ensure that all microstates sharing the same energy are equally probable. For instance, suppose the molecule has 1 = 300 cm -1 and 2 = 600 cm -1 . A microstate with n1 = 2 and n2 = 0 should be as likely to occur as another one with n1 = 0 and n2 = 1 (assuming that nk for k > 2 is the same in both microstates). Although a weak coupling is assumed, it is never explicitly considered here.

The second assumption implies that we will also not discuss the temperature associated with the energy allocated in the translational and rotational modes. Note, however, that these components are relevant when probing the isolated molecule. [START_REF] Fuchs | [END_REF] The third assumption means that we will not deal with the sudden vibrational energy change observed during internal conversion.

B. Reviewing the statistical mechanics' concepts

Before proceeding with the discussion, it is worth revisiting some fundamental concepts from statistical mechanics. 18 The basic quantities we need are the density of states ( ) ( )

Tr , E E H  =-  (5)
and the integrated density of states ( ) ( )

Tr . E E H  =  -   (6)
In these equations, E is the energy of a system with a Hamiltonian H . In Eq. ( 5),  is the Dirac function and, in Eq. ( 6),  is the unit-step function. The trace Tr denotes the integral over the Hilbert space.

The entropy can be defined in different ways. For regular systems in the thermodynamic limit ( )

N → , they all tend to converge to the same result. Nevertheless, when dealing with an isolated molecule, the number N of degrees of freedom is small, and we cannot assume such equivalence. Among several possibilities, we will discuss two entropy formulations, the Boltzmann entropy 29

( ) ( ) ln , BB S E k E  =   (7)
and the Gibbs volume entropy 18,29,[START_REF] Gibbs | Elementary Principles in Statistical Mechanics[END_REF] ( ) ( )

ln . GB S E k E =   (8)
In Eq. ( 7),  is a small energy constant.

With either entropy formulation, the microcanonical temperature can be calculated with Eq. ( 2) and the heat capacity with

1 . T C E -   =    (9)

C. Counting microstates

We can directly count the microstates contributing to either ( ) T temperatures. The superscript (n) denotes the numerical approach to distinguish from other strategies discussed below. To apply the numerical procedure, we note that the molecular energy given by Eq. (3) changes by no less than the energy of the slowest normal mode 1  . Thus, the number of microstates with energy between

E   or ( ) E   , compute
1 k ZP E E kh =+ and 1 k Eh  + ( ) 0,1, 2 k = is ( ) ( ) 1 1 1 , N N k n n k nn h E E E    =  (10)
where ( )

1 1 1 if , 0 otherwise. N k n n k k E E E h E     +   =    (11)
The integrated density of states is

( ) ( ) 1 0 . k kj j E h E  = =  (12)
For more advanced strategies to count microstates considering coupled anharmonic modes, see Refs. 21,23 and references therein.

D. Solving the degenerated problem

As mentioned, the degenerated problem is used below (Section II.E) as an approximation for the molecular problem. Besides that, it is also helpful to check some basic hypotheses and approximations underlying the temperature definition, namely,

• the equivalence between Gibbs volume and Boltzmann entropy,

• the extension of discrete into continuous variable needed to compute derivatives,

• the effect of Stirling's approximation, which is commonly used to simplify the results of the microstates' counting, and

•

The Boltzmann formulation

The number of microstates satisfying Eq. ( 4) is the number

of ways of allocating / / 2 M M E h N  =- quanta in N
vibrational modes. Thus, the number of microstates between E and Eh  + is the number of multisets of length M on N elements,

( ) ( ) ( ) 1! 1 , 1 ! ! M MN N M N E MMNM  +-   + -     = = =      -      (13)
where ( ) ( )

m n and ( ) m n
denote the number of multisets and the binomial coefficient, respectively. The dependence on EM is implicit in M. Taking h  = in Eq. ( 7) allows directly computing the Boltzmann entropy as

( ) ( ) (
)

1! ln . 1 ! ! B M B MN S E k NM  +- =  -   (14)
The extension of this equation from a discrete to a continuous variable ( M EE → ) allows using Eq. ( 2) to get

( ) 1 1, 22 e B B E N E N h T h h k   -      =  + - -+           ( 15
)
where  is the digamma function defined as ( ) ( )

ln d xx dx  =    ( 16
)
and  is the gamma function. The superscript (e) indicates this is an extended-to-continuous approximation.

Alternatively, we can employ Stirling's approximation

( ) ( ) ln ! ln n
n n n -in Eq. ( 14). Then, with Eq. ( 2), we get

( ) ( ) ( ) ( ) 1 22 ln , 2 s B B E N h h T E Nh k    -   +-  =   -    (17)
where the superscript (s) indicates Stirling's approximation. Finally, for large systems, we have 2

N

, and we get ( )

1 2 ln . 2 l B B E Nh h T E Nh k   -  +   =   -   (18)
Schwabl [START_REF] Schwabl | Statistical mechanics[END_REF] offers an alternative approach to computing B T . Instead of using combinatorial to estimate ( )

E



, he rewrites this quantity in terms of the Fourier transform of the Dirac function in Eq. ( 5). The resulting integral is evaluated via the Laplace method, which requires N → .

This approximation once more yields ( )

l B T .
The heat capacity is obtained by replacing B T into Eq. ( 9). With ( ) s B T , we have

( ) ( ) ( ) ( ) 2 2 22 2 41 22 ln , 2 s B B k EE C N N N hv hv E N hv E Nhv   = - - -    -   +-    -  (19)
while ( ) 22 1 ln . 42

l B T gives ( ) 22
l B B Nk E E Nhv C Nhv E Nhv  +     =-      -      (20)

The Gibbs volume formulation

To solve the degenerated problem in the Gibbs volume formulation, we proceed along the same lines as in Section II.D.1 but employing the integrated density of states in Eq. ( 6). For the quantum harmonic oscillator, this quantity is the number of microstates with energy between 0 and E, ( ) (  )

00 1 ! , !! MM M KK N K N E KK MN MN M NM ==   + -      = =           + +  ==    (21)
where the sum runs over the density of states defined in Eq. ( 13) and is evaluated using the hockey-stick identity. [START_REF] Jones | [END_REF] Thus, the Gibbs volume entropy is

( ) ( )! ln . !! G M B MN S E k NM +  =   (22)
Replacing it in Eq. ( 2) gives the extended-to-continuous ( ) M EE → approximation for the microcanonical temperature ( )

1 1 1 . 22 e G B E N E N h T h h k   -      =  + + - -+           (23)
Alternatively, with Stirling's approximation in Eq. ( 22), we have

( ) 1 2 ln , 2 s G B E Nh h T E Nh k   -  +   =   -   (24)
which is identical to ( ) l B T in Eq. ( 18) (we will explore this fact in Section II.D.3). Because of this equality, the heat capacity ( )

s G
C is also given by Eq. (20).

Comparing Boltzmann and Gibbs volume formulations

Let us now compare the previous results for the degenerated problem for a particular example of a system composed of six quantum harmonic oscillators, all with frequency 1015 cm -1 . Its number of states These values are much bigger than the vibrational frequency (0.13 eV in our example). Thus, we can compute the temperature by numerical differentiation of the entropy [Eq. ( 2)], using either ( ) N → , Boltzmann and Gibbs volume temperatures tend to the same values, 17 as we can see comparing Eqs. ( 18) and (24).

E  to get ( ) n B T or ( ) E  to get ( ) n G T .
We also can see in Figure 1-bottom that the extended-tocontinuous approximations ( )

( ) l B T
in Eq. (18). Although this approximation is fully justified for large systems, it has a major conceptual implication: it shifts the system from a Boltzmann to a Gibbs volume description, as made evident by Eqs. ( 18) and ( 24), where we can see that ( ) ( ) ls BG TT = .

Several authors, [START_REF] Callen | Thermodynamics and an Introduction to Thermostatistics[END_REF]33,[START_REF] Schroeder | An introduction to thermal physics[END_REF] starting from the Boltzmann formulation, already apply the large N approximation at the level of the number of microstates [Eq. ( 13)]. Because this approximation reduces Eq. ( 13) to Eq. ( 21), they inadvertently work in the Gibbs volume formulation.

E. The temperature of an isolated molecule ,

ZP Ak k E v NN  ==  (25)
and the harmonic mean 1 11 . The rationale for adopting A  is the direct analogy to the solution of the degenerated problem. Note, for instance, that Eq. ( 24) can be rewritten as

H k k N   -  =    (26)
( ) 1 2 ln . s ZP ZP G ZP B E E E T E E Nk -   + =    -   (27)
Thus, it is natural to search for an approximated solution to the non-degenerated problem with the same functional form. On the other hand, for a given E, there are many more microstates involving changes in the number of quanta in the low-frequency than in the high-frequency modes. Thus, it may be helpful to weigh the average frequency toward the low-frequency modes. If we take as weight the maximum number of quanta we can add to a specific vibrational mode, max /

kk n E h = corresponding to the microstate ( ) max 0, , , ,0 k n , we get max max , kk k H k k n n   =   (28)
providing a rationale for adopting the harmonic mean of the frequencies.

On comparing these two possibilities, we will see in Sections II.E.2 and II.E.3 that the arithmetic mean does better than the harmonic mean.

Test case: The temperature of H2O2

Let us now work out a specific example of a small, isolated molecule, H2O2, for which we can count the microstates and compute the temperature numerically. H2O2 has four atoms and N = 6 vibrational modes with harmonic frequencies 346, 956, 1330, 1469, 3704, and 3709 cm -1

(computed with the geometry, frequency, noncovalent, extended tight-binding (GFN1-xTB) method [START_REF] Grimme | [END_REF] ). The arithmetic mean frequency is

A  = 1919 cm -1
, and its harmonic mean frequency is behavior at low energies, delivering too high-temperature predictions.

H  = 1015 cm

Model validation

The equivalent degenerated model with arithmetic mean frequency gave an excellent approximation for the temperature of the isolated H2O2. Nevertheless, is this a general result or just a coincidence? To check this point, I have computed the Gibbs volume temperature for 30 arbitrary systems (10 with N = 6, 10 with N = 7, and 10

with N = 8) with frequencies randomly distributed between 300 and 3000 cm -1 . The systems with N equal 7 and 8 do not correspond to molecules because N must be 36 at N -. Nevertheless, they are useful to assess the errors. The systems with N = 8 have about 10 9 microstates up to 4 eV [see Eq. ( 21)]. For each energy value, the relative error

( ) ( ) ( ) ( ) / / / / ns B G B G H A n BG TT T T   - = (29)
between the numerical ( ) Most of the error is in the low-energy region, and they are independent of the system size. The mean relative errors in Table 1 T also with arithmetic mean frequency is a bit bigger, 8%. The mean relative errors with harmonic means are always larger than those with arithmetic mean. They are also dependent on the system size, reducing for bigger N. It is not shown in the table, but most of the error when using the harmonic mean frequency is in the low-energy region. Given the better performance of A  over H  , only the former will be used in Section III.

The better performance of the arithmetic over the harmonic mean is fortunate. It implies that low frequencies, commonly associated with anharmonic modes, are not more relevant than high frequencies in temperature determination. The microcanonical temperature is not an injective function of the energy. 18 This means we may have the same temperature for different energies. Therefore, for an energy E, the comparison between a temperature based on a continuous energy variable and that based on a discrete energy variable should be made considering the local mean value of the latter around E. Using the mean relative error as the margin of error is an effective way of considering the mean value of the discrete function.

III. DISCUSSION

F. Leading terms in the molecular temperature

The asymptotic expansion of ( ) s G T [Eq. (24)] at E → gives ( )

22 3 1 . 12 a G BB E Nh TO Nk k E E   = - +   ( 30
)
The leading term is linear in E and does not depend on the particularities of the molecule, only on its size through N.

In cluster studies, it is common to estimate the microcanonical temperature using the phenomenological formula 10,36 2 , where kin E is the mean kinetic energy of the cluster. For a harmonic system, 2 kin EE = and, therefore, Eq. ( 31) corresponds to the leading term in the Gibbs volume temperature expansion, Eq. ( 30).

kin cluster B E T Nk  (31)
If, instead of the Gibbs volume, we start from the Boltzmann temperature ( ) s B T (Eq. ( 17)), the asymptotic expansion at E → gives

( ) ( ) ( ) 22
1 2 2 1 1 . 1 12 a B BB Nh E hv TO N k k E E  - -  = - +  -  (32)
As expected, ( ) ( )

aa BG
TT → when N →. However, for small systems, the leading term depends on the molecular properties through  . Moreover, the ( )

1 N -term in the denominator causes ( ) a B T grows faster than ( ) a G
T with the energy.

Andersen et al. 15 derived the microcanonical temperature of a small isolated system starting from the Boltzmann entropy and searching for the average energy of a canonical system equilibrated at the microcanonical temperature. Their result for an ensemble of N (nondegenerated) quantum harmonic oscillators in the hightemperature limit is (

)

. 1 ABH B E T Nk = - (33)
Considering that in this limit E hv , their formula matches the leading order of the ( )

s B
T expansion in Eq. (32) .

A common assumption in phenomenological models is that the heat capacity (Eq. ( 9)) of the isolated molecule does not depend on the temperature, so a simple linear

relation 0 E CT E =+
holds between energy and temperature. 2,8 According to Eq. ( 30), such a hypothesis is satisfied for a molecule with N vibrational degrees if . 2 3 3

ZP E Nh E  = (34)
Note that the radiative cooling should bring the molecule to the nonlinear region even if this relation is initially satisfied.

G. Effect of the harmonic approximation

The isolated molecule model developed here is purely based on a harmonic approximation. Because the density of microstates as a function of the energy of a fully coupled anharmonic system grows faster than that of a harmonic system, 23 the harmonic microcanonical temperature should be an upper limit for the exact temperature. (Remember, the temperature is the inverse of the derivative of the density of states.) However, how far off the harmonic temperature is? I will not fully address this question here, but I will show two examples to provide insights into the errors we may expect when using the models derived in this paper.

Aieta et al. reported in Ref. 21 the integrated density of microstates ( )

E 
at energy E for fully coupled anharmonic systems, computed with the Wang-Landau Monte Carlo method using the Paradensum algorithm. 25 Among the systems they tested, hydrocarboxyl radical (HOCO) and anthracene are particularly interesting because they correspond to two molecular limits: the first is small and flexible, while the second is large and rigid.

The anharmonic temperature of these two molecules, computed by numerical differentiation of ( )

E  [see Eqs.
(2) and ( 8)], is shown in Figure 4. These temperatures can be directly compared to the Gibbs volume harmonic temperatures ( ) s G T , which are also shown in the figure. They were computed with the same harmonic frequencies used in Ref. 21 . , solid line) and with a fully coupled anharmonic model (dots). The fully coupled anharmonic model is calculated by numerical differentiation of the integrated density of states reported in Ref. 21 (see text).

As expected, the harmonic temperatures are higher than the anharmonic ones. Both approximations agree qualitatively, but the numerical differences are not negligible. For anthracene, the anharmonic temperature is about 150 K smaller than the harmonic one over the entire For HOCO, the anharmonic temperature is about 300 K lower than the harmonic one for intermediary excitation energies between 0.05 and 0.1 au (1.36 -2.72 eV). For lower energies, the comparison is unclear due to the strong oscillations in the anharmonic result. The anharmonic temperature is much smaller than the harmonic one at higher energies. The relative error reaches 20% at 0.2 au. This deviation may indicate the limitations of the harmonic approximation for small, flexible systems.

Assuming that the harmonic microcanonical temperature is an upper bound for the exact microcanonical temperature is reasonable for molecules, where the anharmonic modes fully couple. This hypothesis, however, may fail in systems with internal rotations and decoupled modes. For instance, Belega et al. 36 proposed that the microcanonical temperature of molecular clusters could be estimated by replacing N in Eq. ( 31) with an effective Neff factor computed from the reduced averaged kinetic energy in each mode. Their approach applied to the water pentamer yielded a temperature increase of up to 100 K compared to the harmonic approximation.

H. Experimental comparison between Boltzmann and Gibbs volume

Currently, there is an ongoing debate on which function is the correct way of computing the entropy for an isolated system. While Hilbert, Hänggi, and Dunkel 18,29 claimed that, among several entropy models, only Gibbs volume fulfills all three laws of thermodynamics simultaneously, Swendsen and Wang 19 argued that Gibbs volume fails to satisfy the postulates of thermodynamics and make incorrect predictions for systems with nonmonotonic energy densities. (See Ref. 37 for an account of this debate.)

Although this discussion goes beyond the scope of this paper, the present results create a way for experimentally testing which entropy, Boltzmann or Gibbs volume, delivers the best theoretical prediction. We could think of a setup where a small, isolated molecular ion is photoexcited in an electrostatic ion trap. 4,6,38 The laser pulse duration should be such as it allows thermal equilibration before radiative decay starts but without losing too much energy resolution. Then, monitoring its radiative cooling would inform on its peak temperature (before any cooling). As discussed in Section II.E. According to Eq. ( 20) and using  = 922 cm -1 , this temperature corresponds to E = 14.7 eV, yielding a theoretical estimate of ( ) s G C = 12.0 meV/K [Eq. (20)], in excellent agreement with the experiment. Note that this result is smaller than the high-energy limit, B Nk = 15.0 meV/K (N = 174), due to the nonlinear dependence of the temperature on the energy.

If no chemical process occurs, the hot, isolated molecule is expected to cool down due to blackbody radiation. [6][7][8] Since the seminal papers by Sellgren 40 and Leger and Puget, 12 some specific near-infrared (IR) spectral features of interstellar dust have been assigned to polycyclic aromatic hydrocarbons (PAHs). 41 These earlier works raised the hypothesis that those IR features arise from small, 5-10 Å molecules instantaneously heated to about 1000 K by UV radiation. Working on observational data, Leger et al. [START_REF] Hendecourt | Proceedings of the 135th Symposium of the International Astronomical Union[END_REF][START_REF] Leger | [END_REF] estimated a color temperature of 600 K, corresponding to a peak temperature of 950  150 K. Then, using the Einstein heat capacity model and assuming 6-eV photoexcitation, they predicted that the IR features should be due to carbon species with about 60 atoms.

We can revisit this assignment by using the microcanonical temperature derived here. Let us assume a peak temperature of T = 950  150 K and that the arithmetic mean frequency is  = 1100  100 cm -1 . Let us also allow a broader variation in the UV excitation, EUV = 5  1 eV. Using Eq. ( 24), the number of atoms is exp 1 2. 3 T delivers a good approximation for the microcanonical temperature within an error of 3% compared to the exact harmonic estimate.

UV at B E h N h k T      = -+      (35) Figure 5.

Boltzmann

Gibbs volume

Numerical counting

( ) We should also bear in mind that the microcanonical temperature models discussed here are strictly valid for molecules that are not subjected to external fields, populate a single electronic state, and whose vibrational modes can be well-represented by a harmonic approximation. There are two advantages of working within these boundaries. First, they provide simple analytical expressions that allow a straightforward temperature estimation beyond the classical regime.

Second, they enable the development of a conceptual understanding that is sometimes hidden when applying advanced numerical approaches. Nonetheless, the present modeling should be extended beyond these approximations in the future. The primary goal should be to tackle anharmonic effects to increase the prediction accuracy.

  number of microstates, we must determine the number of solutions of Eq. (3). This means the number of ways we can distribute different values of nk among the vibrational modes, keeping E constant. Despite the simplicity of this formulation, it is a formidably tricky problem to solve analytically. Fortunately, as I will also show later (Section II.E.3), we can get an approximated solution for the number of microstates by supposing the ensemble of vibrational frequencies   k  can be replaced by a single average frequency  . The degenerated problem,

the impact of the 2 N

 2 approximation, also commonly applied to simplify the results. In Subsections II.D.1 and II.D.2, we will discuss the solution of the degenerated problem in the Boltzmann and Gibbs volume formulations. Although the solution in the Boltzmann formulation is well known, 26, 27 I am not aware of any demonstration using the Gibbs volume formulation so far. Subsection II.D.3 compares both results and shows how the Boltzmann temperature tends to Gibbs volume temperature for large N.

  energy resolution is determined by , which is h for the degenerated problem. We may expect vibrational energy excess of about 2-4 eV for typical photoexcited systems.

Figure 1

 1 Figure 1.

1 .

 1 Defining the average frequency My hypothesis to compute the temperature of an isolated molecule is that there is an average frequency  for the N-degenerated problem that is a good approximation for the isolated molecule, composed of N different frequencies. I tested two possibilities for the equivalent degenerated model (EDM),  given as the arithmetic mean 2 1

- 1 .

 1 Up to 4 eV, it has 169,501 microstates, distributed as a function of E as shown in Figure 2-top. The number of states Gibbs volume temperatures of the isolated H2O2 computed by numerical differentiation of the entropy [Eq. (2)] are shown as dots in Figure 2, middle and bottom. The approximated results using the arithmetic mean frequency in the six-degenerated problem are shown in the middle and the harmonic mean in the bottom figure.

For a vibrational energyFigure 2 )

 2 Figure 2.

T

  temperatures and the temperatures of the equivalent degenerated model using Stirling's approximation was evaluated. These errors are shown in Figure3for temperatures computed with arithmetic mean frequencies.

Figure 3

 3 Figure 3. T  Relative error between the numerical temperature and Stirling's approximated temperature of the equivalent degenerated model with arithmetic mean frequencies computed for 30 random systems with N = 6, 7, and 8. Top: Boltzmann; bottom: Gibbs volume.

Figure 4 .

 4 Figure 4. Microcanonical temperature of HOCO (blue) and anthracene (red) computed with the harmonic model using Eq.

  ) temperatures for a few molecules (arithmetic mean frequency). N is the number of vibrational modes. Bottom: Density map of possible molecular systems with temperature 950  150 K, UV excitation 5  1 eV, and arithmetic mean frequency 1100  100 cm -1. The density grows from yellow to blue. The dots illustrate the absorption peak of a few PAH molecules with excitation energies taken from the literature: anthracene, 42 naphthalene (neutral and cation), phenanthrene, pyrene (neutral and cation),43 coronene, dicoronylene, ovalene, 44 and fullerene.45 Letting all three parameters vary according to a random normal distribution with standard deviation given by the error bars, we arrive at Nat = 61  38 for 10 6 random systems. This mean number of atoms agrees with the previous calculations. Nevertheless, as we can see in the density plot at the bottom of Figure 5, it also indicates that a large variety of possible molecular systems falls within the parameter constraints. The figure also shows a few PAH molecules ranging from small to large sizes. Except for naphthalene neutral and cation, they all fall in the highdensity region, as we would expect from PAH mixtures that should occur in the interstellar media. 44, 46 This analysis considers only molecules returning ground state via radiationless pathways. If they fluoresce, the peak temperature after re-equilibration in the ground state should correspond to a vibrational excess discounting the energy of the emitted photon. IV. CONCLUSIONS This paper explores the microcanonical temperature of an isolated molecule. The molecule is treated as an ensemble of quantum harmonic oscillators, and the temperature is derived in the frame of Boltzmann and Gibbs volume entropies. Although the direct solution to the problem seems not to be accessible analytically, we discuss how approximations based on an equivalent degenerated model using average frequencies can adequately describe the molecular temperature. The discussion tracks the effect of several hypotheses and approximations on the results, including the entropy functional, the type of average frequency, the extension of discrete into continuous variables, the use of Stirling's approximation, the assumption of a large number of degrees of freedom N, and the expansion in the leading orders. Table 2 summarizes the several expressions discussed in the text. The comparison between Boltzmann and Gibbs volume approaches shows that the temperatures converge to the same values for large N. The difference is already negligible for a ten-atom molecule. However, for smaller molecules, the difference is significant, and it should be possible to probe it experimentally. In particular, we discuss how a common numerical approximation corresponds, in fact, to a shift from the Boltzmann into the Gibbs volume formulation.We also see how the microcanonical temperature relates to previous phenomenological approaches. For a large vibrational energy excess, the common assumption of a linear relationship between temperature and energy (constant heat capacity) is valid. Nevertheless, dissociative and radiative cooling should bring the system to nonlinear regions.For practical purposes, for a molecule with more than ten atoms, I recommend using ( ) temperature for large N) and arithmetic mean of frequencies. For smaller molecules, the temperature can be computed with either ( ) point, it is unclear which one delivers the correct answer.

  confirm this picture.( )

	s T with arithmetic G
	mean frequency has a mean relative error of 3% no matter
	N. The mean relative error of ( )

s B

Table 1 .

 1 Mean relative error T

		between numerical and

Table 2 .

 2 Summary of the microcanonical temperatures derived under different approximations from Boltzmann and Gibbs volume entropies for a degenerated set of N quantum harmonic oscillators with frequency  . For the non-degenerated problem,  is the arithmetic mean frequency of the spectrum.  is the digamma function. For an isolated molecule well-represented by harmonic vibrational modes,( )

s G

ACKNOWLEDGMENTS M.B. thanks the funding provided by European Research Council (ERC) Advanced grant SubNano (Grant agreement 832237). The author is grateful to Professor A. Verga, who brought Schwabl's alternative derivation of B T to his attention, and Professor Michele Ceotto for discussions.