Heat transport considerations in the mathematical analysis of the photoacoustic and photothermal effects

In this work, we solve the problem of modeling the generation of an acoustic pulse produced by the incidence of a pulsed laser light upon an elastic material. Our concern is about the heat transport during the absorption of electromagnetic radiation. We assume that the pulse duration is of the order of nanoseconds, and asses if under these conditions the contribution of the heat transport in the sample is an essential consideration in the description of the phenomena or if we can ignore it in the model. We begin with the energy balance analysis over the initial interaction of radiation with matter in the context of the formulation Meixner-Prigonine which is called the linear irreversible thermodynamics to describe the induced temperature field. Then we carry a momentum balance which yields the macroscopic elasticity equations with a heat source for the induced pressure field. Once established the equations for temperature and displacement fields, we solve them for the one-dimensional case, showing that the induced pressure has two components, one fast component and one slow component which is due to heat transport in the sample, which is one of the main contributions of the paper.


Introduction
Historically photoacoustic (PA) pulses have been generated by an optical source in two different ways, one called modulated mode and the other pulse mode configurations, the name indicates if the excitation is a periodic function in time, or a single pulsed excitation respectively. It was in the modulated mode that the phenomena of sound generation by an optical source was first discovered. This was mainly due to a periodic accumulation of thermal energy.
The advent of laser technology introduced a change in the way the photoacoustic pulses can be generated. Now it is possible to produce a light pulse of duration in the order of nanoseconds with high power, that is equally capable of producing observable photoacoustic pulses, this joined with new sensors technology of piezoelectric capacitors has enabled one to measure small photoacoustic signals in a broad range of frequencies.
But the underlying physics of the nature of both kind of photoacoustics pulses produced in pulsed mode or in modulated mode seemed to be different in fundamental ways. The resulting pressure wave is due mainly to thermal expansion. Therefore, as the thermal response of any material is slower than its elastic mechanical response, it should exist both a fast pressure pulse due to the mechanical response and a slow sound pulse due to the thermal response of the material. It is in this context that we reexamine the basic equations of thermoelasticity, and suggest the proper origin of the signals that are observed in the photoacoustic pulse generation, dealing with the mathematical analysis of pulse generation in pulsed mode, to show the differences between the slow and fast mechanical perturbations, hereafter referred to as PT and PA components.
Due to the small duration of the pulses when are generated in pulsed mode, the mathematical analysis of this situation has relied on one fundamental assumption, called the stress and heat confinement hypothesis. This hypothesis permit us to analyze mathematically the acoustic pulse generation overlooking the thermal response of the material, and focusing only on the mechanical part. In this work, we also address the implications of this assumption and what information is lost or not under such an assumption. We show that the mathematical analysis can be carried out without using this hypothesis, but the thermal response of the material needs to be solved completely instead.
Since the time scale of these events is related to the time duration of the pulsed laser used as optical excitation source, then it is necessary to identify which are the main sources of heat transport during the phenomena at the interfaces to properly chose the boundary conditions of the modeling. A vast bibliography dealing with the mathematical analysis of these two modes of generating photoacoustic pulses can be found in the references [1][2][3][4][5][6][7][8][9][10][11][12].
A rather concise historical development of the photoacoustic effect, its sensing techniques, and outstanding applications can be found in the tutorial in [13].
We organize this paper as follows: Initially, we discuss the energy conversion from radiation into heat, the photothermal phenomenon. We describe the conversion of the local absorbed optical energy as the heat source that induces thermoelastic deformation. Then we analyze the balance equations for the heat and linear momentum, whose equations model the temperature field and the displacement field, that compose the deformed state. Then we solve the one-dimensional case as a testing example considering the case in which a beam of light is impinging on a flat metal slab. Then analyze how legitimate is the omission of heat transport in the model. It means that we discuss the conditions for ignoring the temperature field, which implies the incorporation of the converted electromagnetic radiation into heat directly into the momentum balance equation, which implicitly is neglecting the heat transport inside the absorber. This condition is the heat and stress confinement hypothesis.
It is shown that under the assumptions of the Meixner-Prigonine formulation of the linear irreversible thermodynamics, the heat is transferred outside the absorber mainly by electromagnetic radiation, modeled by the Stephan-Boltzmann law, which will establish Robin's boundary conditions of problem analyzed.
Once the solution for the thermal response is found, then we get the solution of the wave equation for the scalar velocity potential, which is related to the pressure field which corresponds to the mechanical part. This procedure is carried on through the method of the Green's functions.
For purposes of proof of consistency, we compute numerical calculations, based on previously reported cases [10,14], where the heat transport is neglected, and then we compare with numeric calculations where the PT component is introduced.

Photo-thermal effect and thermo-elasticity
In this section, we derive the equation for the temperature field considering the conversion of electromagnetic radiation into heat which we call the photothermal (PT) effect. The absorption of energy from the external optical field follows the Beer-Lambert law [15], section 16.6.
This law describes how the electromagnetic radiation is absorbed into the matter. Here ( ) q z r is the density of the flux of energy in the z-direction, α is the absorption coefficient and q 0 is the incident energy flux upon the surface of the absorber. The current analysis is for linear elastic materials, and we are using the notation as in [16] for example. The electromagnetic energy is absorbed in the sample and accounts for a local increase of the internal energy of the sample.
Due to the small time scale of a duration of the phenomena, we need to justify some thermodynamic considerations regarding the generalizations of the second law of the thermodynamics for a non-equilibrium initial state as in [17,18]. To avoid the complexity of using a microscopic description of the system, we use an approach called the Meixner-Prigonine formulation. This is called the linear irreversible thermodynamic, and consist of four hypothesis: (a) the extension of the local equilibrium; (b) the extension of the validity of the second law of thermodynamics to irreversible process; (c) the use of linear constitutive equations (Fourier law and other transport properties); and (d) the validity of the symmetry of the Onsager reciprocity relations. These conditions are assumed implicitly through the rest of the work [19]. These assumptions permit us to deal with thermodynamic quantities in a dynamical sense, i.e. that they depend on time.
For the current discussion, let us recall the Landau's [16] description for the heat propagation in an elastic solid, including his tensor notation. In this formalism, one can consider that the optically induced heat (per unit of time and unit volume), can be expressed as ¶ ¶ T ; S t where S is the entropy per unit volume. For the condensed phase, the energy lost by the radiation phase is gained by the condensed phase and assuming the emissivity E = 0 and absorbance to be A a = a q e z 0 , therefore Here  q is the heat flux which units are watts per square meter J/(s m 2 ), S has units of  ( ) J Km 3 , q 0 is the electromagnetic energy flux with units J/(s m 2 ), and α the absorption coefficient with units of m −1 , f (t) is a rectangular function related to the time duration of the pulse. Further, we assume that the wave vector of the electromagnetic field is parallel with the z-axis. We also assume that the field begins to interact with the sample at z=0, and the light absorption occurs from that point and up to L following the Beer-Lambert decay law. Being then that the sample's thickness L where the absorber ends abruptly.
In (2), q 0 defines the light's energy density at z=0. It means that (2) describes the relationship between the transferred optical energy a a q e z 0 , the heat transport   · q and the increment in sample's local temperature T. In the Landau's formalism [16], he gives the Helmholtz free energy for an elastic material when an increase in the temperature is present [16], equation (6.1), Here u ll is the trace of the strain tensor, K is the bulk modulus, μ is the shear modulus, β is the thermal expansion coefficient, and F 0 (T) is the initial free energy within the interaction volume, is the strain tensor, u i is the displacement field and, δ ij is the Krocnecker's delta.
From the well know thermodynamic relation, = - ¶ ¶ S F T we get the related entropy given by negative of the partial derivative of the free energy with respect to the temperature. After applying this to (3), we obtain that It implies that the increase of entropy is proportional to the increase in volume. We now require to describe the heat flux,  q , due to the temperature gradient, ∇T. In it we use the Fourier's law k = -  q T , where κ is the thermal conductivity. Along with it, we substitute (4) in (2) to get, We recall that the heat capacities C p and C v fulfill the relationship [20], pg. 53, b where C p and C v are the specific heats of the elastic material at constant pressure and constant volume respectively. Substituting (6) in (5), Assuming that for the initial undeformed state there is no dilatation or contraction then , then the temporal change of the initial entropy is after combining equations (7) and (8), Thus, from the balance of energy, we have arrived to the description of the optically induced temperature field which is the photothermal effect.

The momentum equation and the scalar velocity potential
Now independently we followed an approach similar to that introduced by Arias and D Achenbach [21] to obtain the equations that describe the displacement field of the thermoelastic expansion. This equation results from the application of the second Newton's law which is the linear momentum balance, and corresponds to the thermoelastic wave equation which accounts for the propagation of the elastic deformation field [11,21,22] namely Here σ ik is the internal stress tensor, ρ is the medium density and ü i are the second time derivatives of the displacement components. The induced stress and pressure are calculated by the derivatives of the Free Helmholtz's energy with respect to each component of the deformation tensor [16], therefore it is given by (10) can be written as Here equations (9) and (12), are central for the analysis. Notice that these form a set of coupled partial differential equations and with appropriate boundary conditions, describe the evolution of the optically induced temperature and displacement fields. The set of equations (9) and (12), correspond to the thermo-elasticity or PT equations [2,21,23]. We proceed to simplify the momentum equation, to be rewritten in terms of the scalar velocity potential function, f s . This is defined as u. This velocity potential is related to the longitudinal modes and fulfills the wave equation,

The small stress approximation: the photothermal equation
Prior to decouple equations (9) and (12), we notice that when no dissipative forces are included, the system is at a state of free thermal expansion (i.e. no external forces are applied) and no internal stresses can be induced. Therefore σ ik =0 which implies that the trace of the strain tensor is simplified as u ll =β (T−T 0 ) ; more details can be found in [16], §6 equation (6.3). Therefore, after substituting u ll =β (T−T 0 ) in (9), we obtain This is the explicit form of the heat equation, with the heat source included, and decoupled from equation (12).
Interestingly the heat source, related to the optically absorbed energy, preserves its analytic dependence with the Beer-Lambert law.

The thermo-elastic equation under the heat and stress confinement hypothesis
In this section, we derive the thermo-elastic equation under the hypothesis of heat and stress confinement, and compare this model with the complete description given by (13), which is our main contribution to remark their differences, namely that it is the time derivative of the temperature field which is the source of the acoustic waves. From equation (14), we replace the time partial derivative of the temperature field in equation (13), as a result, we get . where σ is the Poisson ratio. With these considerations in mind, (15) turn out to be Notice that the right hand term in equation (16) corresponds to the source of the longitudinal acoustic field represented by the speed potential f s , under propagation. In it, the first term is related to thermal diffusivity. After introducing the dimensionless variables ¢ = -T T T T 0 0 and ¢ = z z L where T 0 is the reference temperature, and L is the width of a one-dimensional slab, the second factor of the right hand side of equation (16) can be written as , therefore the stress and thermal confinement condition is given by If this is the condition is met then one can ignore the term  k T C 2 p in equation (16), assuming that heat diffusion is negligible during the excitation pulse [8]. For example in a aluminum sample we see that this condition is fulfilled, so we expect that in our calculation the non-heat transport approximation is a good approximation. In that case then equation (16) is reduced to This is the well known PA wave equation for the scalar velocity potential [7,10,21]. A practical approximation for applications in elastic solids and human tissue is to assume s » 1 2 . Additionally, if the sample does exhibit only hydrostatic compression and no resistance to shear then, μ=0 and therefore σ=1/2. A noticeable example is biological tissue which is assumed that since it is ∼70 % water, this approximation can be applied since no transversal components propagate [8,16,[24][25][26]. Hence, in the equation (18) and it can be simplified to obtain the general analytic expression for the PA wave as reported elsewhere in the literature [2,3,7,10,21].
is the volumetric density of energy per unit time that is absorbed in the sample through nonradiative decays; otherwise known as the PA source or heat function [2,[8][9][10].
In terms of pressure, p, such that r = -f ¶ ¶ p t and if in addition we introduce the Grüneisen number , which we assume that is independent of the frequency and temperature, but in general it is a function of frequency, ω, and temperature, [27], pg. 27, [11,22,28], then equation (16) can be rewritten as This expression is a good enough approximation for describing the presence of the PA signal alone. It works well for a large number of applications where the hypothesis of stress and heat confinement holds. In all-optical PA sensing techniques, the PA and PT components can be sensed at once, depending on their relative magnitude and time lag due to the heat transport speed. The PT velocity depends on the thermal parameters κ, C p . We will use a Robin's boundary condition that models the heat lost by electromagnetic radiation. Other boundary conditions will modifiy the rate of heat transport and the problem should be modified accordingly to the boundary conditions for heat transport.
In this way, we have obtained the equation (1) in the [14] and the equation (5) in [10]. The latter reference covers a review of the most common deduction (the non-heat-transport approximation). In [10] equation (19) is solved using the Green's function method, including a discussion where the pulses generated are interpreted when they are detected by piezoelectric capacitors.

Solution to the Heat equation
In this section, we proceed to solve the heat equation with a source associated with the incident light pulse. The setup of the simulation is as follows: at = t 0, we make a laser pulse to fall upon an aluminum slab of width L, during a time interval of Δt=80 ns which corresponds to the pulse duration. This pulse will heat the aluminum slab producing a volume expansion results in a series of pressure pulses rebounds which are ideally measured at the end of the slab at z=L and at the begin of the slab z=0.
First we solve the heat equation (14), which solution is needed to calculate the source term of equation (13). We assume that the absorption occurs within the region 0<z<L, so we have that (13) can be written as    We also need boundary conditions for this problem. The absorber can interchange heat with the water in three ways, by radiation, convection, and conduction. Across the flat interface, the heat radiated as thermal radiation follows the Stefan-Boltzman law Here dσ is the cross section of the absorber and dt the time elapsed which in principle and accordingly to the assumption of the Meixner-Prigonine formulation can be as small as we want, here σ B =5.670 367×10 −8 W m −2 K −4 is the Stefan-Boltzmann constant, also, the sign − refers to the left side and + to the right side of the boundary respectively. The heat crossing the boundary toward the absorber due to thermal conductivity is where κ 1 the thermal conductivity of the absorber and κ 2 that of the surrounding water.
The heat passing to water by convection is Hereh is the convection coefficient which depends on the fluid dynamics around the flat surface, nature of the fluid if it is a liquid or a gas, the orientation relative to the gravity or if changes of phase are occurring at the interface.
The first case we analyze is that in which we assume that there is no difference of temperature between the absorber and the water. This case happens when the water can change its temperature so fast as to follow the temperature of the absorber. Therefore by equations (24) and (26) there is no heat transport by radiation or convection. In that case, we have that = k dQ 0 implying the boundary condition is the continuity of k ∂ z T and T across the interface. And also we need The other case is that when there exists a difference of temperature between the absorber and the water. This is the case when water cannot change its temperature as fast as that of the absorber.
We ask that no heat is not accumulated at the boundary, so we have that + = where n is the direction of the normal vector to the surface interface pointing outside the absorber, namely to -z at the z=0 plane and +z, at the z=L plane. Also here σ B =5.670 367×10 −8 W m −2 K −4 is the Stefan-Boltzmann constant, T 0 the water temperature, κ the thermal conductivity of the absorber. In what follows we assume that Robin's boundary condition is the one that better describes the heat transport outside the absorber, this assumption is justified by reason of the large water's heat capacity. Also as a consequence, this assumption is consistent with the implicit supposition that wave pressure is generated inside the absorber and not outside. Also if no convective currents are formed, then we can take the convection coefficienth as zero.
Therefore the general solution to equation (23) is, , the temperature is written as T=T 1 +T 2 +T 3 . Here we defined Numerically it is necessary to include the negative frequencies in the inverse Fourier transform, we believe this has to do with the irreversibility of the diffusion process. In the literature [29], r is known as the complex thermal wave number and n w 2 is the thermal diffusion length. After substituting (28) in (23), it turns out that A q is given by, W s 1/2 m −2 K −1 , also known as, thermal admittance' or, contact coefficient'. We introduce this parameter in equation (29), and obtain  ak a k w w = -+ - We see therefore that the source amplitude of the temperature increment is a function of the parameter α κ and the square effusivity times the frequency acting as an imaginary part. The coefficients a, b are found by the Robin's boundary condition (see [31] p.74).
Applying the boundary condition (27) to (28) at z=0 and z=L we found out that the coefficients a, b are a a a = - . Finally we can find the time domain solution by applying the inverse Fourier transform, we used the method of discrete inverse fast Fourier transform. We note that it is the time derivative of the temperature field in the sample's volume which constitutes the source term of the pressure wave equation (13).

The induced photoacoustic pressure
In this next section, we solve equation (13) by using the Green's function method in this way we assure to fulfill the boundary conditions. In the Fourier domain, for the time coordinate, the Green's function for the scalar velocity potential fulfills the equation Here, we solve the one-dimensional case for the Green's function in a mixed Fourier representation w ¢ ( ) g z z , , , for a three layer system when the light source impinges on the sample at z=0. Therefore we need to solve the wave equation for each medium consisting of water, metal, water: The solutions for each interval are, , where ¢ + z means approaching ¢ z from the right-hand side, and ¢z , means to approach ¢ z from the left-hand side.
By means of these boundary conditions we get that the coefficients A, F are: for further details, see [9,10]. Therefore, if we use the relation between the scalar velocity potential and the pressure w where l=1, 2, 3 for each media respectively. The displacement field for z<0 is calculated with the relationship By representing the absorber's density as ρ 2 and since w ¢ ( ) T z , and w ¢ ( ) g z z , , are given by the equations (28) and (40), respectively, then the integral (41)can be evaluated analytically in closed form. It means that, for z>L, Here a j , represents the different amplitudes of the temperature field, with j=1, 2, 3, that in (28), a 1 =a, a 2 =b and a 3 =A q . In turn α j takes the values α 1 =r, α 2 =−r and α 3 =−α, that correspond to three pressures p=P 1 +P 2 +P 3 respectively. r is the self consistent value for the exponential coefficient in e rz , as condition to solve the heat equation in one dimension, in the frequency domain.

Numerical results
In this section we present the results for equations (28) and (42) which represent the temperature and the pressure increment fields for an 80 ns incident laser pulse upon a 3 mm aluminum slab. The parameters for the simulation are shown in table (1). For the Fourier inversion of the pressure field we use the fast Fourier transform algorithm. For the fast pulse we use an upper maximum frequency W=(2π) 60.0·10 7 rad s −1 , which corresponds to a discretization time of 1.66 ns . For the slow components of the acoustic field we use instead W=(2π) 60.0·10 2 rad s −1 which corresponds to the discretization time of 166.66 μs.
In figure 1(a) we show the temperature increment upon the aluminum slab as a function of time, the aluminum interface is assumed to be at the origin of the reference system z=0, the temperature shows an increase, that reverses with the laser's pulse turn off. There is a small offset in the temperature field figure 1(a) (red color), this due to the method we chose to obtain the inverse Fourier transform. This offset is an inherent problem associated with the fast Fourier algorithm applied to the inverse Fourier transform [32]. For the corresponding laser's fluency, we observe a temperature increase of nearly 2 Kelvin degrees in 80 ns. From the analytic solution (28) we observe that we can decompose the temperature field in three terms, that we call T 1 , T 2 , T 3 . Associated with each temperature term we have a pressure field P 1 , P 2 , P 3 respectively. These pressure fields occur at different scales of time; being that the pressure fields P 1 and P 3 are of the order of 80 ns then these corresponds to the fast pressure pulse, we show P 1 +P 3 at z=0 in figure (1(a)) (blue color) the fast component of the pressure field, and on figure 1(b) the slow component P 2 of the pressure field. We note the that the amplitude of the fast pulse is approximately 3.3×10 7 greater that the amplitude of the slow pulse.
In figure 2(a) we compare the fast component of the pressure field with the pressure field obtained with the non-heat transport model which for our parameters range verifies the condition for heat confinement as a good one; both models show multiple replicas of the initial pulse. The details of these results and the interpretation of their respective spectra figure 2(b) can be found in [10].
In figure (3(a)) we display the pressure term P 2 and its Fourier spectra in figure 3(b). This plot is the most important results of this work because it shows that the non-heat transport model ignores a slow thermal pulse. This calculation shows that the magnitude of the fast pulse is seven orders of magnitude higher than the slow pulse.  Finally figure 4 shows the comparison in logarithmic scale of the fast pulse and its replicas both normalized with respect their maximum pressures, which shows that in a single pulse experiment, both fast and slow components are well distinguished showing each pulse in their proper time scales.

Conclusions
In this work, we have reexamined the basic equations for explaining the thermoelastic expansion that occurs in the process of generating photoacoustic pulses. The central question we place is about the validity of the hypothesis of heat confinement. We find that when we include heat transport in the modelling, the solution to the heat equation contains an additional source term, that is associated with a slow pressure pulse. This slow pulse is in the scale of milliseconds and has an amplitude several orders of magnitude shorter than the fast pressure component; whose pulse-width is within the nanoseconds scale, and corresponds to the time duration of the excitation laser pulse.
We conclude that if we are not interested in modeling the slow thermal response, then we can ignore the heat transport at least in the regime of the example we studied, due to the fact that its amplitude is much smaller that the mechanical pulse. However, one has to keep in mind the existence of a background low-frequency pressurecomponent.
Minding the material's physical and chemical properties, and the sensing method and configuration (transmittance or reflectance), we asses that this component cannot always can be neglected [9, 12, 33-36].  (2), now only the heat transport, for the field P 2 at z=L. (b) The pressure Fourier spectra at the plane z=L for the slow pressure component P 2 . The thermal confinement condition, as commonly stated, is expressed by inequality (17). The example we deal with, falls in this regime. Materials with different optical absorption, and also different thermal and mechanical impedances as stated here, may require to include the slow heat transport within the model considerations.
The current analysis shows that it is the time derivative of the temperature field which is the source of the acoustic waves. This is a vision that is lost if we ignore the heat transport in the model.
Another observation of the analysis is that the most reasonable boundary condition in the Meixner-Prigonine formulation for the heat equation is that of a Robbin type. In it the main terms transporting heat are electromagnetic radiation given by the Stefan-Boltzmann law, and the thermal conductivity.
Also we interpret that the fast PA component is mainly related to the elastomechanical properties of the sample. These manifested when the sample is excited externally with the laser pulses. From the sample's transport properties it is known [2-4, 16, 37, 38] that the fast PA signal corresponds to the longitudinal propagation of the elastomechanical perturbation. In turn, the slow component, is related to the thermal diffusivity and to viscous and elastic properties of the sample.
The numerical solutions show that the pressure due to the low frequency thermal field is much smaller respect the PA component (7 orders of magnitude), and propagates diffusively at a lower speed than the PA component lagged four orders of magnitude with respect to the PA component. However, this is not a generalization, since this condition may change depending on the acoustic resonant coupling of the transversal mechanical components [9,10,13,16].
Besides the complexity of the problem, we tried to keep the physical and mathematical sophistication to the minimum required, to asses how well the photoacoustic phenomena can be described in classical terms, avoiding the introduction of non-linear terms or even more realistic transport considerations than the ordinary Fourier law, showing that indeed our model can capture most of the main characteristics of the photoacoustic phenomena. Although the model we present can serve as a basic formalism on which more physics can be overlaid, it is far from capturing physical features reported in recent developments where is clear the impact of temperature distribution over the pulsed PA (or laser ultrasound) signal [12,33,34,[39][40][41][42][43].