Emeric Bourasseau
email: emeric.bourasseau@cea.fr

Vincent Dubois

Nicolas Desbiens

Jean-Bernard Maillet

Molecular simulations of Hugoniots of detonation products mixtures at chemical equilibrium: Microscopic calculation of the Chapman-Jouguet State

In this work, we used simultaneously the Reaction Ensemble Monte Carlo (ReMC) method and the Adaptive Erpenbeck Equation Of State (AE-EOS) method to directly calculate the thermodynamical and chemical equilibrium of mixtures of detonation products on the Hugoniot curve. The ReMC method (W. R. Smith and B. Triska, J. Chem. Phys. 100, pp 3019-3027 (1994)) allows to reach the chemical equilibrium of a reacting mixture, and the AE-EOS method (J. J. Erpenbeck, Phys. Rev. A, 46, p 6406 (1992)) constrains the system to satisfy the Hugoniot relation. Once the Hugoniot curve of the detonation products mixture is established, the CJ state of the explosive can be determined. Performing a NPT simulation at P CJ , TCJ , we then calculate the direct thermodynamic properties and the following derivative properties of the system using a fluctuation method: calorific capacities, sound velocity and Grüneisen coefficient. As the composition fluctuates, and the number of particles is not necessarily constant in this ensemble, a fluctuation formula has been developed to take into account the fluctuations of mole number and composition. This type of calculation has been applied to several usual energetic materials: nitromethane, tetranitromethane, hexanitroethane, PETN and RDX.

I. INTRODUCTION

In the classical description of ideal onedimensional planar detonation waves using the ZND model, the propagation of the wave in the explosive is associated with a sudden increase of pressure and temperature which brings the system onto a particular point of its hugoniot curve: the ZND state. At this point (high pressure and temperature), chemical reactions start to occur, and bring the system from the ZND state (on the hugoniot of the inert explosive) to a point on the hugoniot of the mixture of detonation products, the Crussard curve. Only one point of this curve could support a stationary behavior of the reactive wave, the theory using molecular interactions described by exponential-6 potential models. It is also possible to treat a mixture of several different molecules by reducing the mixture to an equivalent pure fluid using the Ree's mixture law in the so-called one-fluid approximation [3]. Nevertheless, the use of such methods necessarily implies several approximations concerning theoretical background as well as potential models, which make the method less accurate and transferable than molecular simulations.

In this work, we then propose a method to calculate the Crussard curve and the CJ state through molecular Monte Carlo simulations. This presents two main advantages. In the first hand, the use of statistical thermodynamics ensures that the exact EOS of the model material is obtained from atomistic simulations without any approximation. In a second hand, the Monte Carlo method allows to perform simulations using a large variety of potential models from simple, mono center Lennard-Jones materials, to complex, fully atomistic molecular systems (including all atoms models, fluctuating charges, ...), leading to more accuracy and transferability. In this paper, the method is tested with very simple potential models, where molecules are represented by single exponential-6 force centers, in the same fashion as in our thermochemical code [2]. This allows a direct comparison of our results with those coming from thermochemical codes. Results concerning more accurate simulations performed with more complex potential models will be published later.

In this paper, we describe the two methods (ReMC and AE-EOS) used simultaneously to obtain the Crussard curve of any given energetic material. The first method is the Reaction Ensemble Monte Carlo (ReMC) method. This method allows to reach the thermodynamical and chemical equilibrium of a molecular system constrained by chemical equations. This method has been proposed at the same time by Smith and Triska [4] and Johnson et al. [5]. The second method is the Adaptive Erpenbeck Equation Of State (AE-EOS) method, proposed by Erpenbeck [6] and adapted by Brennan and Rice [7]. This method constrains the system to converge toward a point on its hugoniot curve where the Rankine-Hugoniot relations are satisfied. The simultaneous use of both methods has been applied to simple systems constrained by only one chemical equation [7]. In this work, we have coupled the ReMC method implying several chemical equations and the AE-EOS method to perform simulations of realistic mixtures of detonation products composed of up to 7 different molecular species implied in up to 3 simultaneous chemical equilibrium.

Using this method, we have calculated the Crussard curve of 5 usual energetic materials: tetranitromethane (TNM), hexanitroethane (HNE), PETN, RDX and nitromethane (NiMe). Then, the CJ states of those five systems have been determined and compared to experimental CJ properties and results obtained from thermochemical calculations.

In a third time, we have performed Reaction Ensemble NPT simulations at the CJ conditions, in order to obtain thermodynamic derivative properties of the detonation products mixtures at the CJ point. Indeed, those quantities are known to play an important role in the behavior of the detonating system. More precisely, the calorific capacity at constant pressure, the sound velocity and the Grüneisen coefficient are needed for both inert explosive and detonation products mixture in order to feed hydrodynamic codes. Unfortunately, it is also very difficult to obtain those quantities by experimental measurements. Indeed, experiments needed to obtain derivative properties in suitable thermodynamic conditions are often too expensive and too complicated to be performed. Monte Carlo methods appear particularly well suited to obtain derivative properties because fluctuations formula can be used to obtain the main derivative properties through a single simulation. Lagache et al., for instance, have derived the expressions to obtain several thermodynamics properties in the isobaric-isothermal ensemble [8]. This method can be used to calculate derivative properties of inert explosives, but this can not be applied to the detonation products mixture. Indeed, the system after the detonation is under chemical equilibrium, and the derivative properties should be calculated in the appropriate statistical ensemble: the Reaction Ensemble. In fact, although this ensemble has already been used in several conditions (see [4,5,7,[9][10][11][12][13][14][15][16] for examples), no derivative properties calculation has already been performed in this ensemble to our knowledge. In this work, we derive the equations to propose a fluctuation formula that allows the calculation of the main derivative properties under chemical equilibrium conditions, performing a single ReMC simulation. This paper is organized as follow. In a first part, we present the two particular methods used to obtain our results: the ReMC method and the AE-EOS method, and we present a way to establish the fluctuation formula used to calculate derivative properties following the example of Lagache et al. [8]. In a second part, simulation details are presented together with the results obtained for the five systems of interest. Then, discussion and conclusions are provided at the end.

II. METHODS

A. AE-EOS Method

In order to compute the hugoniot curve of a system, we employed the Adaptive Erpenbeck -Equation Of State (AE-EOS) method proposed by Brennan and Rice [7]. As discussed in this paper, it is possible to implement this method following several algorithms. We did not follow the idea of Brennan, who chose a succession of NPT simulations to converge to the Hugoniot pressure at a given temperature. Indeed the hugoniot curve is usually plotted in the (us,up) -shock velocity vs particle velocity -diagram because it corresponds to direct observable of the system. Then, it is possible to plot the Hugoniot in the (P,V) diagram, using the Hugoniot relations. The temperature is not a direct measure: its calculation is not possible from the Hugoniot relations, and it generally implies the use of particular approximations. In these conditions, temperature does not seem to be the most relevant quantity to use as a constrain. Instead, we proposed here a way to implement the AE-EOS method based on a succession of NVT simulations. In our simulations, the volume is constrained and the temperature is adjusted in order to fulfill the hugoniot conditions. Finally, this method appears more efficient considering that NVT simulations converge more quickly than NPT simulations.

Starting from an initial configuration, the system is simulated in the canonical ensemble (NVT 1 ) and the following expression, corresponding to the hugoniot difference H g (1) is evaluated by averaging over several hundred of thousand iterations:

H g = E -E 0 - 1 2 (P + P 0 )(V 0 -V ) (1)
where E 0 , P 0 and V 0 are the energy, pressure and volume of the inert explosive before the detonation (the pole of the hugoniot). In fact, this term measures the gap between the simulated thermodynamic state and the real hugoniot state.

After several hundred of thousand iterations (typically 500000 in our case), the temperature of the simulation is slightly modified from T 1 to T 2 (± 10 K). During the following 500000 iterations, H g (2) is evaluated. At the end of this second step, the following derivative is evaluated:

dH g dT (2) = H g (1) -H g (2) T 1 -T 2 (2)
From here, the new temperature of the simulation T 3 is calculated through:

T 3 = T 2 - H g (2) dHg dT (2)
(3) and the simulations NVT 3 is performed for 500000 iterations, during which H g (3) is evaluated.

This process is automatically iterated every 5.10 5 steps, evaluating each time:

dH g dT (n) = H g (n -1) -H g (n) T n-1 -T n (4)
and a new temperature is calculated using:

T n+1 = T n - H g (n) dHg dT (n) (5)
until the hugoniot difference H g (n) has converged to the required accuracy. Once the Hugoniot temperature T hug is reached, an additional N V T hug simulation could be performed in order to compute accurately the pressure P hug . The advantage of this method is its rapidity. The convergence on the hugoniot curve is usually obtained after 5•10 6 steps only performed in the NVT ensemble (this corresponds to 10 temperature changes). This remains accurate even when simulations are performed in the Reaction NVT ensemble, where chemical equilibrium is reached within the first 200000 iterations for each NVT simulations.

B. ReMC method

The goal of the ReMC method is to compute thermodynamic properties of a multi component system at chemical equilibrium. In order to fulfill this condition, a particular statistical ensemble is defined: the Reaction Ensemble. The complete definition of this ensemble, and the rigourous way to obtain its density probability have been done firstly by Smith and Triska [4]. Here, we will only give a brief description of this ensemble and we invite interested readers to see [4] for more details. We consider a system of s different chemical species a i , and begin by writing the usual canonical ensemble density probability for a non reacting system [START_REF] Mcquarrie | Statistical Mechanics[END_REF]:

ρ ens = s i=1 V Ni N i !Λ 3Ni i • exp -βU + s i=1 βN i µ i • Q -1 ens (6)
where N i is the number of molecules of species i, V is the volume of the system, Λ i is the de Broglie wavelength of species i, U is the configurational energy of the system, and µ i is the chemical potential of species i. Q ens is the partition function of the canonical ensemble.

In the Reaction Ensemble, the temperature and the volume (or pressure) are fixed. In comparison with the canonical ensemble, two more constraints are applied to satisfy the chemical equilibrium:

• the number of atoms is fixed for each atom type, so the number of different molecules in the system is controlled by the chemical equation which defines the chemical equilibrium:

s i=1 ν i a i = 0 (7)
for example:

2N H 3 ⇄ N 2 + 3H 2 (8)
• the sum of chemical potentials over the different molecule species implied in the chemical reaction, weighted by stoechiometric coefficients, is equal to zero:

s i=1 ν i µ i = 0 (9)
for exemple:

µ N2 + 3µ H2 = 2µ N H3 (10)
The key point is to establish a Metropolis algorithm that satisfies those constraints over the simulation. To obtain the mechanical equilibrium of the system, the usual Monte Carlo moves can be used: translation, rotation and internal relaxation of molecules. The probability to accept a move is given by [START_REF] Frenkel | Understanding Molecular Simulation -From Algorithms to Applications[END_REF]:

P acc = min 1, exp ρ ens (new) ρ ens (old) (11)
where ρ ens (new) (respectively ρ ens (old)) is the density probability of the new (respectively old) configuration in the statistical ensemble. So, following equation 6, it can be shown that in the case of usual moves, only the total energy of the system varies between the two configurations. We obtain:

P acc = min (1, exp (-β∆U )) (12)
where ∆U = U (new) -U (old).

To obtain the chemical equilibrium, we used an additional move: the so-called reaction move, proposed at the same time by Smith and al. [4] and Johnson and al. [5]. This move consists first in choosing a direction to perform the reaction, secondly in deleting a set of reactant molecules randomly chosen in the system, and finally inserting product molecules. Following the example given before, the reaction move can be:

2N H 3 → N 2 + 3H 2 (13)
where two randomly chosen N H 3 molecules are deleted and one N 2 and three H 2 molecules are inserted, or:

N 2 + 3H 2 → 2N H 3 ( 14
)
where one N 2 and three H 2 molecules randomly selected are deleted, and two N H 3 molecules are inserted.

During this move, not only the energy is modified, but also the respective number of molecules involved in the reaction. Introducing the parameter ξ, positive (respectively negative) if the reaction move is performed in the forward (backward) direction, it is possible to show that the ratio between ρ ens (new) and ρ ens (old) is:

ρ ens (new) ρ ens (old) = s i=1 V ξνi N i ! (N i + ξν i )!Λ 3ξνi i exp -β∆U + ξβ s i=1 ν i µ i (15)
Following equation ( 9), the second term of the exponential should vanish. Nevertheless, as recommended by Ungerer et al. [START_REF] Ungerer | Applications of Molecular Simulation in the Oil and Gas Industry[END_REF], it is important to clarify the definition of the chemical potential in our equations before such simplification.

In equation ( 9), the chemical potential is the one commonly used in chemistry. The zero of the chemical potential scale (the reference state) in chemistry is usually the chemical potential of elements of the periodic table taken in the standard state (pure element at P 0 and T). In the following, this chemical potential scale is noted μ. Equation (9) could then be rewritten:

s i=1 ν i μi = 0 (16)
In equations ( 6) and ( 15), the chemical potential is the one calculated during the simulation through the Widom test insertion:

µ i = - 1 β ln V (N i + 1)Λ 3 i exp -β∆U + (17)
where ∆U + is the change in energy due to the insertion of a molecule i. The zero of this chemical potential scale is obtained for µ i = 0, i.e. for a molecule i whose insertion energy is equal to zero in a system of density equals to Λ -3 i . These two scales have a different reference state and the simplification in equation ( 15) cannot be done straightforwardly. Indeed, it appears necessary to first link these two scales. This is achieved by expressing the chemical potential of a perfect gas in both cases. In the chemistry scale, the chemical potential of a perfect gas is equal to the standard formation free enthalpy of one mole of perfect gas at P 0 and T:

μi(GP ) = ∆ f G 0 i (T ) N A (18)
In the simulation scale, the chemical potential of a perfect gas is equal to:

µ i(GP ) = - 1 β ln 1 βP 0 Λ 3 i exp -β∆U + int (19)
And we can write:

µ i -µ i(GP ) = μi -μi(GP ) (20)
and

s i=1 ν i µ i = s i=1 ν i μi -μi(GP ) + µ i(GP ) (21)
This becomes, following equation (9):

s i=1 ν i µ i = s i=1 ν i μi(GP ) + µ i(GP ) (22)
Finally, introducing ν = s i=1 ν i , and using equations ( 15), ( 18) and [START_REF] Ungerer | Applications of Molecular Simulation in the Oil and Gas Industry[END_REF], simple algebra leads to the following result: [START_REF] Van Thiel | [END_REF] In our case, this expression simplifies because we used only molecules represented by a single exponential-6 force centre (i.e. without internal energy contribution):

ρ ens (new) ρ ens (old) = (P 0 βV ) ξν exp -ξ ∆ r G 0 (T ) RT • s i=1 N i ! (N i + ξν i )! exp (-β∆U ) exp -β∆U + int ξνi
ρ ens (new) ρ ens (old) = (P 0 βV ) ξν exp -ξ ∆ r G 0 (T ) RT • s i=1 N i ! (N i + ξν i )! exp (-β∆U ext ) (24)
Finally, the acceptance probability along a reaction move is given by:

P acc = min 1, (P 0 βV ) ξν exp -ξ ∆ r G 0 (T ) RT • s i=1 N i ! (N i + ξν i )! exp (-β∆U ext ) (25)
Practically, when a reaction move is chosen, the algorithm is the following:

• The direction of the move is randomly chosen (ξ = 0 or 1). This determines the type R of reactant molecules to be deleted and P product molecules to be inserted.

• ν i molecules of each reactant type R are randomly chosen and deleted from the configuration.

• ν j molecules of each product type P are randomly inserted in the configuration (i.e. the place of insertion is randomly chosen).

• The energy of the new configuration is calculated, and the acceptance probability of the move is obtained using equation [START_REF] Van Thiel | report LLNL[END_REF].

• According to the Metropolis algorithm, the new configuration is accepted or rejected in the markovian chain depending on the comparison between the value of P acc with a real randomly taken between 0 and 1.

Combining translation, rotation and reaction moves, it is possible to simulate a chemical equilibrium at a given temperature and density. Nevertheless, the algorithm presented above is limited to low density systems only. Indeed, as density increases, the insertion of product molecules becomes difficult. This limitation should be overcome if one wants to compute hugoniot curves at high density.

The question regarding the insertion of molecule in dense phases is well known, and appears similarly in the Gibbs Ensemble, or in the Grand Canonical ensemble. This problem is usually solved by the use of a statistical bias in the algorithm that enhances the insertion probability. As an example, John Brennan has proposed the use of the cavity bias sampling to improve the insertion of molecules during a reaction move [12]. We rather choose to use the pre-insertion bias [START_REF] Bourasseau | [END_REF], already implemented in our MC code, but the two algorithms remain quite similar. This bias applied to a reaction move consists first of all on inserting the first product molecules in the empty space left by the deleted reactant molecules. If more product molecules have to be inserted (i.e. if ν > 0), the insertion is performed at a preselected location using the following criterium: for each insertion, a set of k locations are tested, and one position x is selected following the acceptance criterion:

P (x) = exp (-βU LJ (x)) W ( 26
)
with W , the Rosenbluth factor:

W = k j=1 exp (-βU LJ (j)) (27)
In our simulation, this bias has been used for the insertion of the P molecules of product, with k=1 for the R first molecules of product, inserted at the place released by the R deleted reactant molecules, and k=20 for the (P-R) last molecules if P > R.

Two important points have to be considered when using a statistical bias. The first one is that, in order to preserve microreversibility of the Metropolis algorithm, it is necessary to always perform the move in the same way, in the forward or backward direction. More precisely, if the move A + B → C + D + E is performed by deleting A and B, then replacing A by C and B by D and finally inserting E using the bias, the move C + D + E → A + B must be performed by first deleting C, D and E, and then by replacing C by A and D by B. The second important point is to take into account in the acceptance probability the fact that the move is not random anymore. In fact, a new term should be introduced in the acceptance probability, and for a biased move, we obtain:

P acc = min 1, exp ρ ens (new)P gen (n → o) ρ ens (old)P gen (o → n) (28
) where P gen (n → o) (respectively P gen (o → n)) is the probability to generate the old (new) configuration from the new (old) one using the bias. In the case of the pre-insertion bias applied to our reaction move, we obtain (see [START_REF] Bourasseau | [END_REF] for detailed description of the bias):

P gen (o → n) = R p=1 1 • P p=R+1 exp (-βU LJ (x p )) W (x p ) (29
) where x P are the P locations where the product molecules are inserted. Note that the second term vanishes if R P.

In the same way:

P gen (n → o) = P r=1 1 • R r=P+1 exp (-βU LJ (x r )) W (x r ) (30)
where x r are the R locations where the reactant molecules were placed before the move. Note that the second term vanishes if P R.

Finally, if the move is performed in the sens for which P > R, the acceptance probability is:

P acc = min 1, (P 0 βV ) ξν exp -ξ ∆ r G 0 (T ) RT • s i=1 N i ! (N i + ξν i )! exp (-β∆U ext ) • P p=R+1 exp (-βU LJ (x p )) W (x p ) (31)
And, if the move is performed in the sens for which R > P, the acceptance probability is:

P acc = min 1, (P 0 βV ) ξν exp -ξ ∆ r G 0 (T ) RT • s i=1 N i ! (N i + ξν i )! exp (-β∆U ext ) • R r=P+1 exp (-βU LJ (x r )) W (x r ) (32)
It is interesting to note that the only input data needed is the ∆ r G 0 (T ) of the simulated reaction.

This standard free enthalpy of reaction is easily obtained from experimental database.

It is also possible to compute a chemical equilibrium combining several chemical reactions. To do this, we used the same algorithm as presented above, but we added a preliminary step consisting on randomly choosing the reaction to perform before each reaction move.

Finally, note that it is also possible to simulate a chemical equilibrium at constant pressure. To do this, an other move has to be used during the simulation: the usual volume change Monte Carlo movement. As in the NPT ensemble, it consists on modifying the volume of the simulation box, performing an homothetic transformation of the system. The acceptance probability of this move is given by [START_REF] Mcquarrie | Statistical Mechanics[END_REF]:

P acc = min 1, exp -β ∆U + P ∆V + N ln V new V old (33)
C. Fluctuation Formula in the Reaction Ensemble

In this section, starting with the Reaction Ensemble partition function defined by Smith [4], we will develop equations following the example of Lagache et al. [8] in order to establish the fluctuation formula for thermodynamic derivative properties. We will use capital letters X for extensive properties, and lower case letters x for the associated molar properties.

As discussed previously, the system is composed of s molecular species assuming the following chemical equilibrium:

s i=1 ν i µ i = 0 (34)
where µ i is the chemical equilibrium of molecular type i and ν i is the stoechiometric coefficient of i. This condition can be considered as an external constrain imposed on the thermodynamic system. As a consequence, Smith and Triska have established the expression of the Reaction Ensemble partition function under chemical equilibrium from a Legendre transformation [4]:

Q = ξ s i=1 (V (ξ) q i ) Ni(ξ) Λ 3Ni(ξ) i N i (ξ)! × exp -β Ĥ (z 1 , ...) dz 1 ... ( 35
)
where ξ is the extent of the chemical reaction used to define all the accessible states of the system: for a given molecular type i, N i (ξ) = N i

• + ξν i , with N i (ξ) the number of molecules of type i for a given state, and N i

• the initial number of molecules of type i.

In equation (35), V (ξ) is the volume of the simulation box, q i is the part of the partition function corresponding to an isolated molecule of type i, Λ i is the de Broglie thermal wavelength, and Ĥ is the configurational enthalpy, Ĥ = U ext + U int + P V , expressed in terms of the relative coordinates z i of the molecules.

As done in previous parts, we set:

q i Λ 3 i = P 0 β × exp -β ∆ f G i • (T ) N A (36)
As a consequence, the Reaction Ensemble partition function can be written:

Q = ξ s i=1 (P 0 βV (ξ)) Ni(ξ) N i (ξ)! × exp -βN i (ξ) ∆ f G i • (T ) N A × exp -β Ĥ (z 1 , ...) dz 1 ... (37)
Then, an average property X in this ensemble is given by:

X = Q -1 ξ s i=1 (P 0 βV (ξ)) Ni(ξ) N i (ξ)! × exp -βN i (ξ) ∆ f G i • (T ) N A × Xexp -β Ĥ (z 1 , ...) dz 1 ... (38)
To obtain the expression of calorific capacity, it is necessary to derive this expression with respect to temperature. Two quantities depend directly on T in equation ( 38): β and ∆ f G i • (T ). Thus, through simple algebra and using properties of perfect gas concerning ∆ f G i • (T ), it can be shown that:

∂ X ∂β P = X Ĥ -X Ĥ - 1 β X N tot -XN tot + X s i=1 N i N A ∆ f H i • (T ) -X s i=1 N i N A ∆ f H i • (T ) (39)
where N tot = s i=1 N i is the total number of molecules which can fluctuate during the simulation, and ∆ f H i • (T ) is the standard entalpy of formation of molecular type i at temperature T.

In fact, in comparison with the fluctuation formula obtained in the isobare-isotherm ensemble [8], the equation (39) takes into account not only the fluctuations of energy, but also the fluctuations of composition and number of molecules. From this equation, taking X = V , the coefficient of thermal expansion α P can be calculated:

α P = 1 V ∂ V ∂T P ( 40
)
The derivation of calorific capacity is a little more complicated. Indeed, it needs the derivation of total energy in regard with temperature, including kinetic energy. But in a Monte Carlo simulation, the kinetic energy K is not taken into account, and as a consequence, only the residual calorific capacity can be calculated during the simulation. Thus, Lagache et al.

have proposed to clearly separate the two contributions, ideal and residual:

C P (T, P ) = C id P (T ) + C res P (T, P ) (41)
with

C id P (T ) = ∂ H id ∂T P ( 42
)
and

C res P (T, P ) = ∂ H res ∂T P ( 43
)
Considering that H id can not be evaluated in regard with K, but that it is possible to give an expression of it in regard with standard ther-

modynamical data (H id = i N i ∆ f H • i
NA ), we obtain:

C id P (T ) = i ∂ N i ∂T ∆ f H • i N A + i N i C • Pi
(44) and with H res = U ext + P V -N tot kT , we can thus obtain C res P (T, P ), and then C P (T, P ), using equation (39) with X = H res and X = N i .

To obtain the other derivative properties of interest, it is necessary to calculate the isothermal compressibility β T whose is the same as in the isothermal-isobaric ensemble:

β T = - 1 V ∂ V ∂P T = 1 V kT V 2 -V 2 ( 45
)
Then it is possible to obtain the calorific capacity at constant volume C V , the sound velocity C S and the Grüneisen coefficient γ:

C V = C P - α 2 P V T β T C S = V C P C V β T γ = V α P C v β T (46)
It is important to note that the quantity X (for example C P (T, P )) is extensive. To obtain the exact molar quantity, it is necessary to derive the molar quantity of X (for example the molar enthalpies h id and h res ), introducing again the N i 's in the derivative, for which fluctuations have to be taken into account. Nevertheless, the only constant quantity during a detonation is the mass of the system, and the use is to employ massic quantities, so we did not evaluate molar quantity in this work.

III. RESULTS

A. Simulation Details

All the simulations performed in this work have been done with the program GIBBS, owned by the Institut Francais du Petrole, the Université Paris Sud, and the CNRS, and developped in collaboration between those three owners and the CEA [START_REF] Ungerer | Applications of Molecular Simulation in the Oil and Gas Industry[END_REF].

In order to reach the chemical equilibrium of a system using the ReMC method, it is necessary to define a priori the set of chemical equations driving the system. To do this, the first step is to choose which type of molecule is expected at equilibrium. In our case, it is known that the detonation products of NiMe, HNE, TNM, PETN and RDX are the following: CO 2 , H 2 O, CO, N 2 , H 2 , O 2 , N O, N H 3 and CH 4 in a fluid phase, and eventually carbon in a solid phase. More precisely, it is possible to neglect some of those molecules, because their molar fractions are known to be very small. So on, table I shows the detonation products taken into account for each studied explosives.

The second step consists on finding a set of linearly independent chemical equations involving those molecules that satisfy the stoechiometry and the chemical specificity of the system. To do this, we applied the formula-vector matrix stoechiometric algorithm devised by Smith

NiMe TNM HNE PETN RDX CO2 x x x x x H2O x x x CO x x x x N2 x x x x x H2 x x x O 2 x x N O x x x x N H3 x x x CH4 x
Csolid x TABLE I: Molecular species taken into account in the detonation products mixture of the studied explosives.

and Missen [21]. It is important to note that the nitromethane case is particular due to the presence of solid carbon. We will show later that the solid phase is not included in the chemical equilibrium, and this is the reason why it has not been included in the formula-vector matrix stoechiometric algorithm used here. This method leads to the following set of reactions:

• NiMe: 2N H 3 ⇄ N 2 + 3H 2 CO + 2N H 3 ⇄ N 2 + CH 4 + H 2 O CO 2 + H 2 ⇄ CO + H 2 O
• TNM (All the carbon atoms are in the CO 2 molecules):

2N O ⇄ O 2 + N 2 • HNE: 2N O ⇄ O 2 + N 2 2CO 2 ⇄ 2CO + O 2
• PETN and RDX:

2N

H 3 ⇄ N 2 + 3H 2 2CO 2 + N 2 ⇄ 2CO + 2N O CO 2 + H 2 ⇄ CO + H 2 O
The potential models used in this study are the exponential-6 potential proposed by Fried et al. [START_REF] Fried | proceedings of the 12th International Detonation Symposium[END_REF] that reproduce the thermodynamic behavior of pure compounds under high temperatures and pressures. A cut-off equal to the half of the box length has been used to reduce the computing time. Periodic boundary conditions together with long range corrections have been used.

All the simulations performed in this work were initiated with a number of molecules corresponding to a system of several hundred of explosive molecules. During ReMC NVT simulations, the respective probabilities of choosing a translation move and a reaction move were 0.8 and 0.2. During ReMC NPT simulations, the respective probabilities of choosing a translation move, a volume change and a reaction move were 0.8, 0.02 and 0.18. The pole conditions used to compute H g in the AE-EOS method are given in table II. The convergence of AE-EOS ReMC simulations were generally obtained after 10 7 iterations, changing the temperature every 0.5•10 6 iterations and calculating the H g value on 0.25•10 6 iterations before evaluating the new temperature. The NPT simulations performed to obtain derivative properties should be much longer: around 20•10 6 iterations to converge, and averages were calculated on 30•10 6 iterations after convergence.

V0 (cm 3 .g -1 ) E0 (J.kg

B. Crussard curves

To calculate the Crussard curves of the four systems without solid phase of carbon (TNM, HNE, PETN, RDX), we have directly applied the AE-EOS ReMC method described previously. Five simulations have been performed for each system, at five different specific volumes between 0.4 and 0.95 cm.g -1 . Monte Carlo results are displayed in figure 1 (TNM), 2 (HNE), 3 (PETN) and 4 (RDX), and compared with results obtained with our thermochemical code using the same potential models [2]. For those four systems, we can see that Monte Carlo results are in good agreement with thermochemical ones. The small discrepancies seen in figures 1, 2, 3 and 4, are probably due to the various approximations used in the thermochemical code, for example in the MCRSR theory [1] and the one-fluid hypothesis. Those results, obtained with the same potential models in both cases, validate our microscopic method to calculate Crussard curves.

Concerning Nitromethane, the simulation process is a little more complicated due to the fact that we are not actually able to explic- itly simulate a chemical equilibrium involving a solid phase. To obtain the Crussard curve of nitromethane, we then applied the following two steps process:

• In a first step, we performed a calculation of the Crussard curve with our thermochemical code, in order to determine the quantity, the molar volume and the molar energy of the solid phase of carbon along the Crussard curve, using the Reevan Thiel equation of state [START_REF] Van Thiel | [END_REF]. • Then, we performed the five AE-EOS ReMC simulations at five different volumes, but with the stoechiometry corresponding to nitromethane minus the carbon trapped into the solid phase. During those simulations, the AE-EOS method is applied with a modified Rankine-Hugoniot relation:

H g = E tot -E 0 - 1 2 (P + P 0 )(V 0 -V tot ) (47)
with

E tot = E f luid + E solid (48)
V tot = V f luid + V solid (49)
where E f luid and V f luid are the energy and volume of the fluid calculated during the Monte Carlo simulation, whereas E solid and V solid are the energy and volume of the solid phase of carbon taken from the thermochemical calculation.

Results obtained concerning the nitromethane appear in figure 5. As far as nitromethane is concerned, Monte Carlo and thermochemical results are also consistent. Nevertheless, in this case, Monte Carlo results depend on the thermochemical ones through the Ree-van Thiel equation of State for solid phase of carbon. This is not satisfactory, because this means that the accuracy of a complex microscopic Monte Carlo simulation is bound to the quality of an empirical equation of state. Moreover, in our simulations, the quantity of solid carbon is determined a priori.

As a consequence, although we propose in this work a simple and successful way to take into account the solid phase in the fluid chemical equilibrium, this may be improved to better simulate the equilibrium between the solid and fluid phases. Some perspectives of this work are presented in the conclusion.

C. Chapman-Jouguet State

Once the Crussard curve has been calculated for each system, the five CJ points could be determined. To do this, we have fitted the calculated points by the following function, which gives the pressure on the Crussard curve in re-gard with η = V /V 0 :

P (η) = a • e -bη + c • η -d (50)
The CJ point is then the tangential point between the Rayleigh line and the Crussard curve. Thus, knowing the pole conditions, we are able to give a value of the derivative of the pressure along the Crussard curve with respect to the volume at the CJ point:

dP (η) dV (CJ) = P (η CJ ) -P 0 V (η CJ ) -V 0 (51)
which is equal to:

dP (η) dη (CJ) = P (η CJ ) -P 0 η CJ -1 (52)
Introducing equation (50), and assuming that P 0 << P (η CJ ), we obtain the CJ point where η satisfies:

a • e -bη + c • η -d (1 -η) = a • b • e -bη + d • c • η -(d+1)
(53) In this method, the CJ points calculated for the five systems are presented in table III. From the calculated P CJ and V CJ , we have determined the detonation velocity D CJ , using the following equation derived from the Hugoniot relations:

D CJ = V 0 (P CJ -P 0 ) (V 0 -V CJ ) ( 54
)
This property is interesting because it can directly be measured experimentally. a Reference [24] for NiMe Reference [START_REF] Van Thiel | report LLNL[END_REF] for TNM Reference [START_REF] Meyer | Explosives[END_REF] for HNE Reference [START_REF] Price | report Naval Surface Weapons Center[END_REF] for PETN Reference [START_REF] Aleksandrov | Combustion[END_REF] for RDX TABLE III: CJ states of the studied explosives calculated using our Monte Carlo method, using our thermochemical code [2], and, considering the detonation velocity, compared with experimental measurements [24][START_REF] Van Thiel | report LLNL[END_REF][START_REF] Meyer | Explosives[END_REF][START_REF] Price | report Naval Surface Weapons Center[END_REF][START_REF] Aleksandrov | Combustion[END_REF].

temperatures and below 1.8 % concerning detonation velocities.

We can also note that the calculated values of D CJ are in good agreement with the experimental values in each of the five cases. The difference between Monte Carlo simulation results and experimental ones is about 1.8 % for nitromethane, 1.6 % for TNM, 1.3 % fo HNE, 0.4 % for PETN and 1.2 % for RDX, which is very satisfactory considering that simple potential models have been used, and considering that experimental uncertainties concerning the positions of the CJ state are in the same order of magnitude. It is already anticipated that MC results would be more significant as realistic potential models become used. It is also important to notice that, considering detonation velocity, thermochemical results are sometimes closer to experiment than Monte Carlo ones. This can be explained by the fact that po-tential models are effective potentials. Indeed, they have been optimized through the use of thermochemical codes, and thus integrate the different sources of errors coming from the different approximations. Table IV shows that the derivative properties calculated at the CJ points using our Monte Carlo method are in good agreement with thermochemical results.

The discrepancies between the two methods are below 3 % concerning heat capacity at constant pressure (C P ), below 1.5 % concerning heat capacity at constant volume (C V ), below 1.8 % concerning sound velocity (C S ) and around 5 % concerning Grüneisen coefficient (Γ). It is important to note that a part of those discrepancies is due to the fact that the (P CJ , T CJ ) are not exactly the same.

The calculation of the sound velocity at the predicted CJ state appears really interesting because it allows to verify that the considered thermodynamic state is really the CJ state. Indeed, the CJ state is the only point in the Crussard curve where the following relation is verified: D CJ = C S + u P (CJ). u P is the particle velocity, obtained from the Hugoniot relation: u P = P -P0 ρ0D [29]. Using the Monte Carlo results presented in tables III and IV, it appears that this relation is verified for the 4 systems, with an accuracy of 0.5 %. This is very satisfactory because it validates in the same time the way to determine the CJ state from the Crussard curve, and the fluctuation formula used to calculate C S .

Those results concerning derivative properties are important because to our known, it represents the first results of microscopic calculations of such properties for systems under chemical equilibrium.

IV. CONCLUSION

We have proposed in this work a new microscopic method to calculate the crussard curve and the CJ state of detonation products mixtures of usual explosives. This method, based on Monte Carlo simulations employing ReMC and AE-EOS methods, has been validated by comparison with thermochemical calculations. The simulations presented here have been performed with simple potential models, and it can be expected that simulations performed with complex potentials would give more accurate results. This is the first perspective of our work.

We also have proposed a fluctuation formula to calculate derivative properties in the Reaction statistical ensemble. This formula takes into account the energy fluctuations, and the fluctuations of composition and number of molecules. This formula has allowed us to calculate heat capacities, sound velocities and Grüneisen coefficients through ReMC NPT simulations, in the 4 cases where no carbon phase was included. Concerning the last case (NiMe), the calculation of derivative properties was not possible with this formula, because of the way the solid phase has been taken into account.

Finally, we have proposed a simple way to take into account the solid phase of carbon which could appear in the detonation products. In this work, the properties of the solid phase are determined by thermochemical calculations via the Ree-van Thiel equation of state, and included in the Monte Carlo simulation. This leads to a satisfactory restitution of the Crussard curve and the CJ properties, but still presents unconsistency in the treatment of the equilibrium between the fluid mixture and the solid phase. This constitutes the second perspective of this work. We could think about several ways to improve this equilibrium. A first idea should be to explicitly treat the car-bon phase in the simulation. Nevertheless, this would strongly increase the computing time. Moreover, the intrinsic properties of the solid phase are not the main goal of this type of study. So, it should be sufficient to take into account the solid phase in a more implicit way. For example, considering that experimentally, the solid phase of carbon seems to appear as agregates, it could be interesting to simulate those agregates via mesoparticules with more or less interactions with the fluid mixture. Moreover, those improvements should allow us to calculate the derivative properties of system under a chemical equilibrium including a solid phase.

FIG. 1 :

 1 FIG. 1: Calculated Crussard Curve of Tetranitromethane (TNM).

FIG. 2 :

 2 FIG. 2: Calculated Crussard Curve of Hexanitroethane (HNE).

FIG. 4 :

 4 FIG. 4: Calculated Crussard Curve of RDX.

FIG. 5 :

 5 FIG. 5:Calculated Crussard Curve of Nitromethane (NiMe).

TABLE II :

 II Pole conditions used to compute Hg in the AE-EOS method.

  Table III shows a comparison between the calculated and measured values of D CJ .

		NiMe TNM HNE PETN RDX
	Monte Carlo Results
	PCJ (kbar) 118.7 149.7 61.42 83.17 93.83
	VCJ (cm 3 .g -1 ) 0.658 0.470 0.780 0.721 0.742
	TCJ (K)	3463 2250 6017 4677 4325
	DCJ (m.s -1 ) 6414 6364 4828 5599 6028
	Thermochemical Results
	PCJ (kbar) 123.1 155.0 62.57 85.85 97.22
	VCJ (cm 3 .g -1 ) 0.653 0.470 0.784 0.721 0.742
	TCJ (K)	3513 2229 6042 4685 4350
	DCJ (m.s -1 ) 6493 6464 4897 5696 6137
	Experimental Results a
	DCJ (m.s -1 ) 6370 6450 4950 5620 6050
	Table III shows that the CJ states calculated
	using our Monte Carlo method are in good
	agreement with thermochemical results. The
	discrepancies between the two methods are
	below 4 % concerning pressures, below 1 %
	concerning volumes, below 1.5 % concerning

  D. Derivative properties at the CJ pointFour ReMC NPT simulations have been performed in order to calculate the derivative properties of TNM, HNE, PETN and RDX. Because of the contribution of the carbon solid state in the CJ point of nitromethane, this system has not been considered. Results are presented in tableIV.First of all, it is important to note that the converged volume obtained at the end of those NPT simulations is consistent with the predicted CJ volume presented in table III.

		TNM HNE PETN RDX
	Monte Carlo Results
	CP (J.kg -1 .K -1 ) 1660 1826 2199 2024
	C V (J.kg -1 .K -1 ) 1568 1593 1977 1800
	CS (m.s -1 )	4927 3434 4177 4494
	Γ	0.619 0.396 0.436 0.536
	Thermochemical Results
	CP (J.kg -1 .K -1 ) 1660 1824 2233 2085
	C V (J.kg -1 .K -1 ) 1575 1592 1998 1827
	CS (m.s -1 )	5008 3492 4233 4551
	Γ	0.59 0.40 0.45 0.57
	TABLE IV: Derivative properties (heat capacities,
	sound velocities and Grüneisen coefficients) at CJ
	point for the four studied explosives.

V. ACKNOWLEDGEMENTS Pr P. Ungerer from IFP is gratefully acknowledged for fruitful discussions about the theory of the Reaction Ensemble and the way to obtain the acceptance probability. E. Auroux from CEA is gratefully acknowledged for his help on thermochemical calculations and for his knowledge about experimental data. All Monte Carlo simulations have been performed with the Gibbs code from IFP, CNRS and the Université Paris-Sud [19].


Stéphane Berciaud

David Lasne

Gerhard A Blab

Laurent Cognet

Brahim Lounis

Photothermal Heterodyne Imaging of Individual Metallic Nanoparticles: Theory versus Experiments

We present the theoretical and detailed experimental characterizations of Photothermal Heterodyne Imaging. An analytical expression of the photothermal heterodyne signal is derived using the theory of light scattering from a fluctuating medium. The amplitudes of the signals detected in the backward and forward configurations are compared and their frequency dependences are studied. The application of the Photothermal Heterodyne detection technique to the absorption spectroscopy of individual gold nanoparticles is discussed and the detection of small individual silver nanoparticles is demonstrated.

I Introduction

Several optical schemes have been used to perform the detection of nano-objects at the single entity level. Together with constant improvements in synthesis and characterization of nanosized materials, those highly sensitive methods make possible the development of new nano-components such as plasmonic devices 1,2 or single photon sources 3 . Until recently most optical detection methods of single nanometer sized-objects were based on luminescence.

Single fluorescent molecules or semiconductor nanocrystals have been extensively studied and are now widely implemented in various research domains ranging from quantum optics [START_REF]Special-Issue[END_REF] to life science 5 . However, luminescence based techniques suffer for some shortcomings, mainly associated with the photostability of the luminescent nano-object itself. Concurrently, for relatively large nanoparticles, Rayleigh scattering based methods have recently demonstrated a great applicability, especially for single metal nanoparticles (NPs) spectroscopy 6,7 or biomolecules imaging 8 . However, as the scattering cross-sections of the particles decrease with the sixth power of the diameter, these methods are limited to the study of rather large NPs (diameter > 20 nm).

Since the absorption cross-section of these NPs scales with the volume, an interesting alternative to Rayleigh-scattering relies solely on absorption properties. Indeed, excited near their plasmon resonance, metal NPs have a relatively large absorption cross section (~6×10P - 14 P cmP 2 P for a 5 nm diameter gold NP) and exhibit a fast electron-phonon relaxation time in the picosecond range 9 , which makes them very efficient light absorbers. The luminescence yield of these particles being extremely weak 10,11 , almost all the absorbed energy is converted into heat. The increase of temperature induced by this absorption gives rise to a local variation of the refraction index. This photothermal effect was first used to detect gold NPs as small as 5 nm in diameter by a Photothermal Interference Contrast (PIC) method 12 . In that case, the separated beams of an interferometer, only one of which propagating through the heated region. The sensitivity of this technique, though high, is limited. In particular, when high NA objectives are used, depolarization effects degrade the quality of the overlap between the two arms of the interferometer.

We recently developed another photothermal method, called Photothermal Heterodyne Imaging (PHI) 13 . It uses a combination of a time-modulated heating beam and a non-resonant probe beam. The probe beam produces a frequency shifted scattered field as it interacts with the time modulated variations of the refractive index around the absorbing NP. The scattered field is then detected through its beatnote with the probe field which plays the role of a local oscillator as in any heterodyne technique. This signal is extracted by lock-in detection. The sensitivity of PHI lies two orders of magnitude above earlier methods 12 and it allowed for the unprecedented detection of individual 1.4 nm (67 atoms) gold NPs, as well as CdSe/ZnS semiconductor nanocrystals. In addition, since the PHI signal is directly proportional to the power absorbed by the nano-object, this method could be used to perform absorption spectroscopy studies of individual gold NPs down to diameters of 5 nm 14 The goal of this paper is to give the theoretical framework of the PHI method and to compare the expected signals with the experimental results. In the following section, after a qualitative description of the principle of the PHI method, we will present an analytical derivation of the PHI signal. For this purpose the theory of light scattering by a fluctuating dielectric medium is used as in our case, the photothermal effect occurs on dimensions much smaller than the optical wavelength of the probe beam and previous derivations of photothermal techniques do not apply. The experimental study of the PHI signal is presented in Section III. Its variations with the modulation frequency are detailed both theoretically and experimentally in section IV. In section V, we present the results of the implementation of PHI spectroscopy of individual gold NPs followed by preliminary results obtained with silver NPs. In the last section, we briefly discuss further implementations of PHI.

II Theoretical model for the PHI signal

Throughout this article, we consider an absorbing nanosphere with radius a much smaller than the optical wavelengths embedded in a homogeneous medium whose thermal diffusivity is

C D κ =
(with κ the thermal conductivity of the medium and C its heat capacity per unit volume). When illuminated with an intensity modulated laser beam with average intensity

IB heat B , the NP absorbs ( ) [ ] t P abs Ω + cos 1
, where Ω is the modulation frequency, heat abs I P σ = with σ the absorption cross section of the particle. At distance ρ from the center of the NP (Fig. 1), the temperature rise can be derived using the heat point source model for heat diffusion [START_REF] Carslaw | Conduction of heat in solids[END_REF] :

( ) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - Ω ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - + = Δ th th abs r t r P t T ρ ρ ρ πκ ρ cos exp 1 4 , (1)
with Ω = D r th 2 the characteristic length for heat diffusion. The corresponding index of refraction profile is

T T n n Δ ∂ ∂ = Δ
, where T n ∂ ∂ are the variations of the refractive index with temperature (typically, 10P -4

P KP -1 P).
In the following, we will rather use the induced susceptibility profile:

T T n n Δ ∂ ∂ = Δ 2 χ
, where we denote n the non-perturbed refractive index of the medium.

The geometry of the problem is represented in Figure 1. A focused circularly polarized probe beam interacts with the time modulated susceptibility profile. At the beam waist, a plane wave-front is considered. Thus, the incident probe field can be written ( ) In practice, rB th B is smaller than the probe beam's wavelength (see below). Hence, a full electromagnetic derivation of the scattering field is necessary to further evaluate the detected signal. We used the model introduced by Lastovka 16 as a starting point for our derivation [START_REF] Chu | Laser Light Scattering[END_REF] .

( ) + - = e ρ E ρ i t i E t
First, the interaction of the incident field with the local susceptibility fluctuations give rise to a polarization ( )

t , ~ρ P : ( ) ( ) ( ) t n t t i , , , ~2 0 ρ E ρ P ρ χ ε Δ = (2)
The expression of the scattered field at a point M (with R OM = ) is derived by introducing the Hertz potential [START_REF] Nisbet | Proc. Roy. Soc. Lond. A[END_REF]

( ) ( ) ( ) ∫ ∫ - + - = ρ - R ρ R ρ P ρ R Π m c t t t dt d t ' ' , ' 4 1 ,
( ) ( ) ( ) t t , , R Π R E × ∇ × ∇ = (4)
At this point, it should be stressed that only components

( ) t , R E Ω
of the electric field which are frequency shifted by Ω from the incident field frequency will contribute to the PHI signal.

Moreover, in the far-field, an analytical expression of ( )

t , R E Ω
can be derived. Indeed, considering that the spatial extension of ( )

t , ρ χ Δ
is microscopic (of the order of rB th B and thus small compared to ρ R -), points M in the far-field region will be such that

R R r e ρ ρ R . - ≈ - =
. The zeroth-order of this expression will be taken in the denominator of the integrand of Eq. ( 3). Furthermore, considering that

Ω >> th m r c
, retardation effects can be neglected in the temporal integral of Eq. (3). Then, using Eq. ( 4) and introducing the vector

( ) R z s i n e e k k k - = - = Δ λ π 2
, the electric field scattered at frequency shifted by ± Ω with respect to incident field frequency ω writes in the far-field region:

( ) ( ) ( ) ( ) t t t , , 2
1 , R E R E R E - Ω + Ω Ω + = (5)
with:

( ) ( ) ( ) ( ) ( ) [ ] t R k i abs s e E I R n P T n n t Ω ± - ± + ± Ω Ω × ∇ × ∇ ∂ ∂ = ω θ ε κ π 0 2 0 2 1 4 2 , , e R E (6)
and

( ) ρ ρ ρ θ ρ ρ k 3 . exp d r i r I th th ∫ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ± Δ + - = Ω ± , ( 7
)
After integration, one obtains:

( ) ( ) ( ) ( ) Ω ± Ω = Ω ± , , , θ θ π θ ig f r I th 2 2 (8)
The functions f and g write after the introduction of the parameter ( ) ( )

2 sin 2 4 , θ λ π θ Ω = Δ = Ω D n r u th k : ( ) ( ) ( ) ( ) ( ) ( ) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + - - + + + = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + - - + + + + = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 u u u u g u u u u u u f (9)
The far-field component of the scattered field (which varies as R 1 ) has the form:

( ) ( ) ( ) ( ) ( ) ( ) ( ) + - Ω × × ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ Ω Ω Ω + Ω Ω ∂ ∂ - = e e e R E R R t R k i abs s e R E t g t f C P T n n t ω θ θ λ π 0 2 sin , cos , 2 , (10)
The detected signal will now be evaluated from this expression. It originates from the interference between

( ) t , R E Ω
and a local oscillator field proportional to the incident probe field. Assuming that the wave front of this local oscillator is spherical in the far-field region, the power distribution of the resulting beatnote per unit of solid angle in the direction ( ) φ θ, writes:

( ) ( ) ( ) ( ) ( ) [ ] θ θ θ λ φ θ θ θ cos 1 sin , cos , 2 sin , 2 2 + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ Ω Ω Ω + Ω Ω ∂ ∂ = Ω t g t f P P w C P T n n d d P d LO i
abs PHI (11) with PB i B the incident power of the probe, w the waist of the probe beam in the sample plane and PB LO B the power of the local oscillator. Experimentally we deal with two configurations where we detect either the backward or forward contribution of the scattered field. In the backward configuration, PB LO B is the reflected probe power at the interface between the glass slide and the medium surrounding the NPs. In the forward configuration, it is the transmitted probe power through the sample. More precisely, Integration of Eq. ( 11) leads to the beatnote power arriving on the detector and oscillating at frequency Ω. Ιn the backward/forward configuration, it writes:

( ) ( ) ( ) [ ] ) sin( ) cos( . 2 2 , / / 2 / / / t G t F w C P T n n P t P F B F B abs i F B F B F B Ω Ω + Ω Ω ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ = Ω λ π α η (12)
where ηB B B and ηB F B are the transmission factors of the optical path (in pratice, ηB B B ~ ηB F B ) and:

( ) ( )[ ] ∫ + Ω Ω = Ω max min sin cos 1 θ θ θ θ d f F F B (13) ( ) ( )[ ] ∫ + Ω Ω = Ω max min sin cos 1 , 1 / θ θ θ θ θ θ d g G F B
If one assumes that the collection solid angle is 2π in both experimental configurations,

0 min = θ (resp. 2
π ) and 2 max π θ = (resp. π ) should be used for the forward (resp. backward) configuration (see Fig. 1). Note that

( ) Ω F B F /
is in phase with the modulation of the heating and

( ) Ω F B G /
is in quadrature. Finally, demodulation of the signal power by the lock-in amplifier leads to the PHI signal which magnitude is simply proportional to ( )

Ω F B P / : ( ) ( ) ( ) 2 / 2 / 2 / / / . 2 2 Ω + Ω ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ = Ω F B F B abs i F B F B F B G F w C P T n n P P λ π α η (14)

III Experimental results and characterization of the signals

In the following, the experimental details will be given. A schematic of the setup is presented in Figure 2. A non resonant probe beam (632.8 nm, HeNe, or single frequency Ti:Sa laser) and an absorbed heating beam (532 nm, frequency doubled Nd:YAG laser or tunable cw dye laser) are overlaid and focused on the sample by means of a high NA microscope objective (100x, NA=1.4). The intensity of the heating beam is modulated at Ω by an acousto-optic modulator.

As mentioned previously, the PHI signal can be detected using two different configurations.

In the case of the detection of the backward signal, a combination of a polarizing cube and a quarter wave plate is used to extract the interfering probe-reflected and backward-scattered fields. In order to detect the forward signal, a second microscope objective (80x, NA=0.8) is employed to efficiently collect the interfering probe-transmitted and forward-scattered fields.

The intensity of the heating beam sent on the NPs ranges from less than Application of silicon oil on the sample insures homogeneity of the heat diffusion. The index of refraction of that fluid and its thermal conductivity are close to those of common glasses. Thus, there is no sharp discontinuity neither for the thermal parameters nor for the refractive indices at the glass-silicon oil interface and we can consider that the NPs are embedded in a homogeneous medium.

Figure 3 shows a two dimensional PHI image of individual 10 nm gold NPs corresponding to the backward (Fig. 3a) and forward (Fig. 3b) signals. Both images show no background signal from the substrate, indicating that the signal arises from the only absorbing objects in the sample, namely the gold aggregates. In both cases the heating intensities were the same (~500

kW/cmP 2 P)
and the NPs are detected with high signal-to-noise ratio (SNR) greater than one hundred.

According to Eq. ( 14), the resolution of the PHI method depends on the probe and heating beam profiles and also on the dielectric susceptibility profile. Since the spatial extension of the latter is much smaller than the size of the probe beam, the transverse resolution is simply given by the product of the two beams profiles. In figure 4, we study the resolution by imaging a single gold NP with two sets of beam sizes. In the first case, we used low aperture beams with Gaussian profiles measured at the objective focal plane (Fig. 4(a)). In the second case, higher aperture beams are used and their profiles contain diffraction rings arising from the limited aperture of the microscope objective (Fig. 4(b)). In both cases, the transverse profile of the PHI signal is in very good agreement with the product of the two beams profiles. In the first case (Fig. 4 , the PHI signal shows no saturation behavior, rather it is accompanied by fluctuations in the signal amplitude and eventually irreversible damage on the particle [START_REF] Link | [END_REF]20 .We further confirmed that the peaks stem from single particles by generating the histogram of the signal height for 321 imaged peaks (see Fig. 5(b)). We find a monomodal distribution with a dispersion of 30 % around the mean signal. At a given heating wavelength heat λ , the PHI signal is proportional to the absorption cross-section σ. According to the dipolar approximation of Mie theory [START_REF] Kreibig | Optical properties of metal clusters[END_REF] and with a SNR > 10 13 .

To clearly demonstrate the linearity of the PHI signals with respect to the volume of the NPs, the size dependence of the absorption cross section of gold NPs (at 532 nm, close to the maximum of the SPR) was studied for NP diameters ranging from 1.4 nm to 75 nm. As expected, a good agreement with a third-order law of the absorption cross-section vs the radius of the particles was found 13 .

IV Frequency dependence

We will now compare the amplitude of PHI backward and forward signals and discuss their dependences with respect to the modulation frequency. respectively (Eq. 12). In the calculation, we considered a 10 nm gold NP absorbing PB abs B = 500 nW, and we used

1 4 10 ~-- ∂ ∂ K T n , 3 1 6 10 . 2 ~--m K J C , 1 2 8 10 . 8 . 1 - - = s m D
, nm w 520 = in Eq. 14.

Two main features can be seen from Figure 6: first, and as expected in any thermal process, the PHI signals exhibit a low-pass behaviour. Second, a clear difference between the forward and the backward normalized signal magnitudes is found at low frequencies, whereas both signals decrease identically in the high frequency domain. Qualitatively, one can understand these features by considering the susceptibility profile as a scattering object with a characteristic size th r .

We introduce the size parameter

λ π ξ th r n 2 =
of the scattering theory [START_REF] Bohren | [END_REF] and the cut-off

frequency 2 2 2 ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = Ω λ π n D C corresponding to 1 = ξ .
On the one hand, for high frequencies (

C Ω >> Ω
) when the scattering object is small compared to the incident probe beam wavelength ( 1 << ξ

), the angular distribution of the scattered light is symmetric with respect to the focal plane [START_REF] Bohren | [END_REF] On the other hand, for low frequencies (

C Ω ≤ Ω , 1 ≥ ξ
) , the scattering object is large enough so that forward scattering becomes significantly more efficient than backward scattering [START_REF] Bohren | [END_REF] , as shown in figure 6 where

( ) i F F P P α
Ω is greater than

( ) i B B P P α Ω .
Fig. 6 (a-b) also present experimental data points from measurements performed on a single 10 nm gold NP. The modulation frequency ranged from 100 kHz to 15 MHz. As mentioned above, experimental data were adjusted by the model resulting from Eq.( 14) with D as the only adjustable parameter. For both configurations, limitations of the collection efficiency by the objective lenses were neglected. A good quantitative agreement between the experimental data and the theoretical forms of ( ) Ω is obtained for the same value of ( )

1 2 8 10 1 . 0 8 . 1 - - ± = s m D
in both configurations. This value is in good agreement with the technical specifications of the embedding medium. We also reported the values of the typical sizes of the scattering objects rB th B mentioned above and reported them on the top axis of Fig. 6 (a-b).

To further test our theoretical description of the PHI signals, Eq. ( 14) can be used to estimate numerically the PHI power detected at the detector in both configurations. A single 10 nm gold having an absorption cross section of individual 10 nm gold NPs). A qualitative agreement is found between the modeling and the experimental results. The discrepancy is certainly due to the assumption of a spherical reference wave in the calculation of the beatnote power. This optimal case is not reached in our experiments even with high NA microscope objectives.

In any case, the forward PHI signal is significantly larger than the backward one, due to a more efficient heterodyne amplification in the forward direction ( B F α α >> ). However, in the case of shot noise limited detection, the two configurations are expected to give identical SNRs, except at frequencies lower than Ω c for which the transmission scheme is expected to give better SNRs as forward scattering becomes significantly more efficient than backward scattering. Experimentally, shot noise limited detection is difficult to obtain for low frequencies, and it is only in the forward detection scheme and for sufficiently high incident probe powers ( )

mW P i 1
> that it is achieved. Nevertheless the excess noise was low enough to obtain close-to-optimal SNR in the backward detection scheme (see Fig. 3).

V Absorption spectroscopy of individual metallic NPs

PHI not only allows for highly sensitive detection of gold NPs but also opens new pathways towards absorption spectroscopy of single metallic NP and more generally of non-luminescent nano-objects. Several theoretical models predicted the existence of so-called intrinsic size effects in the optical response of metallic NPs with sizes significantly smaller than the electron mean free path [START_REF] Kreibig | Optical properties of metal clusters[END_REF]23 . Limitation of the electron mean free path as well as chemical interface damping [START_REF] Kreibig | Optical properties of metal clusters[END_REF]24 increase the damping rate of the Surface Plasmon Resonance (SPR), thus leading in the time domain to shorter dephasing times (down to a few fs). Because of such intrinsic size effects, the dielectric constant of the NPs differs from that of the bulk value and must include an additional surface damping contribution. Experimental studies on ensembles of metallic NPs revealed the existence of such effects [START_REF] Kreibig | Optical properties of metal clusters[END_REF]25 . However a quantitative description of those effects was made difficult mostly because of inhomogeneous broadening.

In order to overcome this shortcoming, PHI was used to record absorption spectra of individual gold NPs with diameters down to 5 nm 14 . Figure 7 represents the absorption spectra of 2 single gold NPs with diameters of 33 nm and 5 nm respectively, performed in the wavelength range 515-580 nm (i.e. photon energy range 2.41-2.14 eV). The values of peak resonance energies are not particularly affected by intrinsic size effects. On the contrary, a significant increase in width of the resonance Γ is clearly visible that cannot be described by Mie theory using the bulk values of the gold dielectric constant 26 . We found a good agreement between the experimental widths and Mie simulations by introducing the size dependent correction to the bulk dielectric constant that accounts for the presence of intrinsic size effects in the absorption spectra of small gold NPs 14 . Although the existence of intrinsic size effects in the optical response of gold NPs was unambiguously revealed, part of the damping processes are due to interband transition, which makes it difficult to connect the widths of the Plasmon resonances to the damping rate. Indeed, for gold, the energy threshold for interband transitions lies at ~2.4 eVB B and is preceded by an absorption tail starting at about 1.8 eV.

Consequently, the SPR spectra are asymmetric and defining a full-width-at-half maximum of absorption spectra for gold NPs is delicate, and even impossible for very small particles. A way to circumvent this additional damping would consist in red-shiftingB B resonant energies towards photon energy smaller than the onset of interband transition. This can be done either by embedding spherical gold NPs in a matrix with high refractive index or by studying the SPR of the long axis mode in gold nanorods 6,27 or even by using core-shell NP 28 .

Another possibility consists in using spherical silver NPs instead of gold NPs 29 . Indeed for silver, the resonant energy of spherical NPs embedded in a medium with refractive index 5 . 1 = n is at ~3 eV, whereas interband transitions start at ~3.7 eV 26 . Thus silver NPs are well adapted for studying the size dependence of ultrafast dephasing of the SPR at the single particle level. A limitation of silver NP is however their weak photostability due to photooxydation 30 . This limitation can be surpassed by using encapsulated silver NPs 31 e.g. with PEG (poly(ethylene glycol)) 32 . Figure 8

VI. Conclusion and Perspectives

In conclusion, we have presented the theoretical and experimental characterizations of the Photothermal Heterodyne Imaging technique which allows for the detection of individual non fluorescent nano-objects such as metallic nanoparticles (gold or silver) or even semiconductor nanocrystals 13 . The photothermal heterodyne signal was derived through a model based on classical scattering theory. Experiments on individual gold nanoparticles provided results in good agreement with our calculations. Photothermal Heterodyne Imaging, which is not limited by most shortcomings of fluorescence microscopy such as scattering background or photobleaching allows for long observation times with high signal to noise ratios.

The high sensitivity provided by PHI, coupled to a better control of the NPs synthesis should allow original studies of light-NP interaction processes and we foresee diverse applications in different fields such as in plasmonics or bioimaging. For instance, gold NPs have been used as labels for biological applications 8 , 33 . In this context, PHI is potentially the most sensitive method to detect gold NPs in live cells by optical means. For such applications, the local increase of temperature in the vicinity of the NP must not perturb the biological function of the specimens under study. Since the thermal conductivity of metals is much higher than that of the surrounding medium, the temperature inside a spherical NP can be considered as uniform and equal to the temperature at its surface. This temperature writes: a P T

  frequency of the incident field).

  cB m B being the speed of light in the medium. As the polarization variations P ~ are localized in the vicinity of the particle, at the observation point M, the total electric field is simply related to this Hertz potential by:

  transmission coefficients at the glass/sample interface.

  the desired signal-to-noise ratio and the NP size to be imaged). Backward or forward interfering fields are collected on fast photodiodes and fed into a lock-in amplifier in order to extract the beat signal at Ω. Integration time of 10 ms are typically used. Images are formed by moving the sample over the fixed laser spots by means of a 2D piezo-scanner. The size distributions of the gold NPs were checked by transmission electron microscopy (data not shown) and are in agreement with the manufacturer's specification (typically, 10% dispersion in size). The samples were prepared by spin-coating a solution of gold NPs diluted into a polyvinyl-alcohol matrix, (~1% mass PVA) onto clean microscope cover slips. The dilution and spinning speed were chosen such that the final density of NPs in the sample was less than 1 µmP -2 P.

  (a)), the full-width-at-half-maximum (FWHM) of the PHI profile is equal to 365 ± 5nm and reduces to 235 ± 5nm in the second case (Fig.4(b)), in accordance with the products of the beams profiles which are respectively equal to 360 ± 5nm (FWHM) and 213 ± 5nm (FWHM).The linearity of the PHI signals with respect to IB heat B was checked on individual 10 nm gold NPs (see Fig.5(a)). Even at high intensities ( )

Figures 6 (

 6 Figures 6 (a-b) present the theoretical dependence of the normalized signal magnitude ( ) i B

1 .

 1 leading to identical backward and forward PHI signals. Furthermore, the in-phase components This is the typical response of a driven system in a dissipative medium. As a consequence, one finds that the high frequency domain and decrease as Ω 1 .

  (a) shows a three dimensional representation of PHI image of individual PEG-coated silver NPs with average diameter of 5.3 nm. For this image, a modulated diode laser (Coherent Compass) emitting at 405 nm was used as heating laser with heating intensity of only 2 50 ~cm kW. As expected, we found that the SNR obtained with silver NPs is about 10 times higher than that of gold NPs of the same size, when identical heating intensities are used to excite the NPs at the peak of their SPR resonance. The PHI signal measured on PEG-coated silver NPs was stable during more than one minute (see Fig.8(b)). For those silver nanoparticles, the size distribution checked by TEM (data not shown) was poor: 20%. This translates in an dispersion of 90% in the distribution of the measured PHI signals (653 imaged peaks (see Fig.8(b)). For such experiments, sensitivity at the single particle level is all the more necessary as the heterogeneity in size and shape is obviously larger for silver NPs than for gold NPs. In addition, expected progresses in the synthesis of silver NPs such that narrower size distributions NPs can be obtain should permit a quantitative study of the SPR of individual silver NPs much smaller than 5 nm in diameter.

  inverse of the distance from the NP surface. For 5 nm gold NPs in aqueous medium, a heating intensity of 500 kW/cmP 2 Pleads to a surface temperature of ~ 1.5 K, sufficiently low for live cell integrities.

Figure captions Figure 1 (

 1 Figure captions

Figure 2 (

 2 Figure 2

Figure 3 (

 3 Figure 3

Figure 4

 4 Figure 4

Figure 5 (

 5 Figure 5

Figure 6 Theoretical

 6 Figure 6

Figure 7

 7 Figure 7

Figure 8 (

 8 Figure 8

Figure 1 (Figure 2 (Figure 3 (Figure 4 Figure 5 Figure 6

 123456 Figure 1 (color online)

Figure 7 Figure 8 (

 78 Figure 7

  ( )

	Π	, R	t	which obeys to an inhomogeneous wave equation with the local
	polarization variations ( ) t , ~ρ P	as a source term:

Acknowledgments

We wish to thank D. Ferning and C. Brust for the synthesis of silver nanoparticles and Philippe Tamarat for helpful discussions. This research was funded by CNRS and the French Ministry for Education and Research (ACI Nanoscience and DRAB) and the Région Aquitaine.