Modeling terahertz heating effects on water

: We apply Kirchhoff’s heat equation to model the inﬂuence of a CW terahertz beam on a sample of water, which is assumed to be static. We develop a generalized model, which easily can be applied to other liquids and solids by changing the material constants. If the terahertz light source is focused down to a spot with a diameter of 0.5 mm, we ﬁnd that the steady-state temperature increase per milliwatt of transmitted power is 1 . 8 ◦ C / mW. A quantum cascade laser can produce a CW beam in the order of several milliwatts and this motivates the need to estimate the effect of beam power on the sample temperature. For THz time domain systems, we indicate how to use our model as a worst-case approximation based on the beam average power. It turns out that THz pulses created from photoconductive antennas give a negligible increase in temperature. As biotissue contains a high water content, this leads to a discussion of worst-case predictions for THz heating of the human body in order to motivate future detailed study. An open source Matlab implementation of our model is freely available for use at www.eleceng.adelaide.edu.au/thz .


Introduction
Applications of terahertz technology are very promising for use in the everyday life, e.g. in security systems [1], for detection of cancer [2], for package inspection [3] and in spectroscopy [4,5].Furthermore there have been recent advances in techniques for performing THz-TDS on various liquids and aqueous samples [4][5][6] and driven by these techniques, there is much interest in exploiting THz-TDS for studies of solvation dynamics of fluids [7], hydration levels in substances [2], proteins in aqueous solution [8], relaxation processes in liquid [9], etc.As with all other electromagnetic waves, terahertz beams will create heat under interaction with a given medium.If a sample is not kept at a constant temperature by an external energy (C) 2010 OSA 1 March 2010 / Vol. 18, No. 5 / OPTICS EXPRESS 4728 source or sink during an experiment, the increase in temperature could induce uncertainties in the extracted parameters.We elucidate this problem by expounding an appropriate model to calculate a given temperature increase in a medium at steady-state.The model can easily be generalized to work for all types of fluid samples and even solids by choosing appropriate input parameters, but since water has been widely investigated within the terahertz bandwidth and its thermal conductivity and absorption coefficient is well defined we choose to use water as a case study.
As we will find in Section 5, this paper shows that in order to create a local increase in temperature of 1 • C, where a THz beam is focused down to a spot diameter of 0.5 mm, the transmitted THz CW power must be more than 0.56 mW.For comparison, quantum cascade lasers can produce CW terahertz beams with output powers up to 17 mW [10].If a sufficient time exposure of such a CW beam is maintained to allow thermal equilibrium, our results indicate a temperature increase of 30 • C for a spot diameter of 0.5 mm -thus possibly affecting the measured material parameters -and one should in this case consider using an external heat sink to keep the sample temperature steady.However, note that accelerated electron sources of terahertz radiation, such as synchrotrons and free-electron lasers, have average powers as high as 20 W [11][12][13].At these levels, THz radiation may possibly boil a water sample and thus this highlights the need to take appropriate precautions.
When specifying a THz source it is often convenient to quote figures in terms of power density, but since the result in this paper has a nonlinear dependence on the beam radius it is more useful to keep the dimensions in terms of watts and then state the beam radius usedone can then easily calculate a given power density afterwards.

The Method
All calculations in this paper are performed on a disc of water with thickness d and radius b.The THz beam with radius a is focused onto the center of the disc of water.An illustration of the boundary conditions are shown in Fig. 1.
The heat equation for a system where the temperature has reached a steady value is given by [14], where k (T ) is the thermal conductivity at temperature T and g (z) is the dissipated power density at the thickness z.
In order to gain physical insight we keep the model simple by considering only transmitted power through the sample.Later in Section 5.3, we will consider the effect of transmittance to find what value of source power corresponds to a given transmitted power.

Assumptions and boundary conditions
In order to solve the heat Eq. ( 1), appropriate boundary conditions for the system need to be found.It is a good approximation to assume that there is negligible loss of heat to the air on the top of the disc, (∂ T /∂ z = 0).We choose the disc to be relatively thick, which means the bottom will stay at a constant temperature (T = T 0 ).At the edge of the disc, the boundary conditions can be chosen as either those for a constant temperature or those for a good heat sink.Both cases will be explored, but it turns out that if the initial values for the different parameters are chosen properly, the final result is the same.We furthermore assume that the function is continuous within the medium.
To achieve an analytical solution to the heat Eq. ( 1) we use Kirchhoff's transformation, [a]: where k 0 = k (T 0 ), which transforms the nonlinear problem into a linear one [15].Due to the symmetry of the incident beam, we can reduce the heat equation to the following form in cylindrical coordinates: where g = αP exp (−αz) / πa 2 k 0 is the transformed power density with α as the absorption coefficient and P as the transmitted power.The corresponding boundary conditions for the transformed temperature, (U (r, z) = T (r, z) − T 0 ), are shown in Fig. 1.

The solution of U(r, z)
By using the two Boolean variables χ and θ , taking on the values 0 or 1, U (r, z) can be written as a single expression, which holds when either 0 ≤ r ≤ a or a < r ≤ b and for the two different boundary conditions at r = b, b.1 and b.2.The definitions of the two variables are given in Table 1.
Table 1.Definition of the Boolean variables χ and θ It turns out that the temperature change is directly proportional to the power, which makes it convenient to carry out the simple transformation: where P 0 = 1 mW and P is a dimensionless constant, which is chosen to be 1.This gives Ũ(r, z) the dimensions of Kelvin per milliwatt [K/mW] and all calculations will be based on this parameter.To find out how much a specific amount of transmitted terahertz average power P (in milliwatts) will heat up the water sample, only a simple multiplication is needed: This transformation will furthermore decrease the number of numerical errors, since the result will be around 0.5 − 2K/mW and not a small number, e.g.
The complete solution of the heat equation is where the variables p n , g n , A n , B n (r) and C n (r) are given by B n (r) = (−1) χ K χ( χ p n a + χ p n r) I χ (χ p n a + χ p n r) (10) where I m (x) and K m (x) denote the modified Bessel functions of order m of the first and second kind, respectively.Equation ( 6) is a generalized model that is applicable to any dielectric material, provided that its thermal conductivity and absorption coefficient are known.The full derivation for Eq. ( 6) can be found in [14].Note that it is shown that the simplification of a uniform absorption profile is unrealistic [14], therefore in this paper we adopt the correct exponential absorption profile.
The beam used in these calculations has a uniform beam profile, where normally a laser beam has a Gaussian profile.However, we argue that it is a sufficiently good approximation to look at a uniform profile especially when the spot size of the beam is much smaller than the diameter of the disc of water, and furthermore the added complication of modeling a Gaussian beam profile is not justified for a calculation of this type where we are finding the final equilibrium temperature.For very small sample sizes, where the spot size becomes significant, a Gaussian beam profile will need to be considered -however, this case is not within the scope of this paper.

Thermal conductivity and the absorption coefficient
The thermal conductivity, k 0 [W/mK], and the absorption coefficient, α [1/m], for liquid water are two quantities that previously have been extensively investigated.Ramires et al. [16] approximate, on the basis of experimental data, the temperature dependence (from 0 • C to 100 of the thermal conductivity of water as a second order polynomial: In 2007, Ellison performed experiments on water that led to an empirical expression for the complex permitivity as a function of both the temperature (from 0 • C to 100 • C) and the frequency (from 0 THz to 25 THz) [17].In order to convert the real part (ε ) and the imaginary part (ε ) of the complex permittivity to an expression for the absorption coefficient the following relation is used: where c 0 is the speed of light in vacuum.The expression can be derived from equations given in [18].The data from [16] is very well documented and the result from [17] is in excellent agreement with [4,5,19], so no further tests of the accuracy of the two formulas are needed.

Initial conditions
It turns out that in order to obey the boundary condition at the bottom of the disc, the disc itself needs to be relatively large (b = 50 mm, d = 15 mm), otherwise Ũ(r, z) evolves nearly linearly in the z direction and does not become zero at the very bottom of the disc.
All the initial values can be seen in Table 2, and are shown as vertical red dashed lines in most of the forthcoming plots.Due to the exponential absorption profile, it is expected that the solution also should behave exponentially.This is verified in Fig. 2 where an exponential tail is seen in the solutions both along the z-and radial coordinates.The summation in the solution, Ũ(r, z), turns out to converge to a final value, i.e. the more elements of the sum that are included the less influence the last term has.In order to obtain meaningful results a proper value for the number of terms to be included is found.In Fig. 3(a) this convergence is shown; after 1000 terms the difference in Ũ(r, z) is around 10 −6 K/W and after 4000 terms the value is around 10 −9 K/W, which means that the numerical error will occur in the 9th significant figure of the result.But since the error in the thermal conductivity occurs in the fourth significant figure [16] it makes no sense to use that many terms in the sum.So in all the following calculations the number of terms will be set to 1000.
In [14], the C n (r, z) term was neglected in the final solution since its influence became negligible.This is not necessarily the case in the calculation performed here, because the considered beam width is much larger than the one used in [14].However, it actually turns out that if the initial conditions given in Table 2 are used, the dependence of C n (r) is at maximum around 10 −6 in all calculations that means the error is about as small as the chosen numerical error.Therefore the term is neglected in all calculations.This also means that it does not matter whether boundary condition b.1 or b.2 is used.Even though the dependence is very small, C n (r) still needs to be taken into account if the radius of the disc is small.This is illustrated in Fig. 3(b), where the two different boundary conditions and the case where C n (r) = 0 are compared.But when the radius of the disc becomes small other uncertainties occur, which means the boundary condition at the bottom of the disc no longer will hold.This is also why the chosen initial value for the disc radius is b = 50 mm.If the chosen thickness of the water disc is too small, the profile along the z-axis loses its exponential behavior, which is unphysical.Thus we have chosen a thickness of 15 mm for the sake of example.(b) The temperature change vs. the beam radius, a, at (r, z) = (0, 0).It is hard to focus a THz beam with a center frequency of 1 THz down to a spot radius smaller than 0.25 mm [20] due to the diffraction limit, and most often the radius will be around 0.5 mm.Note that this quantity is investigated even further in Fig. 5.
Two other significant parameters are the thickness of the sample and the radius of the beam.As it can be seen in Fig. 4(a) the value of Ũ(r, z) increases when the thickness of the sample becomes larger and seems to saturate around 30 − 70 mm.When small values for d are chosen, the profile of the solution, along the z-axis (at r = 0), is almost linear and the behavior becomes unrealistic, i.e solutions with an exponential tail that converge to zero for large z are closer to physical reality.Hence the chosen value is d = 15 mm, and on Fig. 4 seen that the solution is almost saturated at this point, which means Ũ(r, z) will be maximal.In Fig. 2(b), the exponential profile of the solution along the z-axis is clearly seen.
According to [20], a terahertz beam can be focused down to a spot size around 0.5mm in diameter, which leads to the chosen value for a.Looking at Fig. 4(b) it is seen that the beam radius is very significant for the final result, e.g. if the spot size is hypothetically focused down to a radius of 0.1 mm instead of 0.25 mm the amount of heating will double.This motivates us to investigate the influence of the beam width even further in the next section.

Beam width dependence
In Fig. 5 the effect of different beam widths is seen, and even though the heating is large for small beam widths it appears the penetration depth of heating is greatly reduced.This means the transmitted beam mainly will heat up the top layer of water right where the beam hits.According to Fig. 5(g), temperature varies by 20-40% over the area where the beam is incident.Thus the center temperature can be used as the actual worst-case temperature of the sample, and is especially applicable when using a reflection setup where it only is the surface layer of the sample from which optical constants are extracted.

Frequency and temperature dependence
Building upon the work of and [17], the dependence of the initial temperature and the frequency on the temperature increase can now be investigated.
When the initial temperature of the water disc is increased it is expected that the temperature change per milliwatt will decrease, because the water does not need that much heating in order to reach the same total temperature as for lower initial temperatures.This is exactly what is seen in Fig. 6(a).
The frequency seems to have a relatively large effect on the final temperature change -in Fig. 6(b) it is seen that Ũ(r, z) increases almost 40 % from 1.4 K/mW to 2.0 K/mW when the frequency is varied from 0.1 THz to 10 THz.So reduction of the center frequency of the beam lowers the amount of heating.

Application of the model in various situations
In the contour plot of Fig. 7(a) it is clearly seen that the maximum heating takes place on the very top of the disc right in the middle of where the beam is incident.Here the temperature change per milliwatt is 1.8K/mW.Figure 7(b) shows the transmitted power required to give a temperature increase of 1 • C.
Future terahertz devices with applications in security analysis claim to produce powers up to 100 mW at a frequency of 0.65 THz [1], with a spot diameter of 5 mm this gives a temperature increase of 19 • C and if the beam is focused down to a smaller spot the increase in temperature can be up to 170 • C -with such systems it will not be possible to keep the sample at a constant temperature by using a heat sink.But since the exposure time in such cases normally will be fairly short, i.e the temperature will most probably not reach steady-state, the model will provide an overestimate or worst-case of the final heating, so a more advanced model where the time is taken into account is of interest for future work.
Furthermore recent techniques have made it possible using an electron accelerator to make THz beams with average powers over 20 W, which under the right conditions may be sufficient to boil water [11][12][13].The initial conditions given in Table 2 is used in both figures.

Application to terahertz time domain spectroscopy
Until now we have only considered CW sources, but since terahertz time domain spectroscopy, where a pulsed beam is emitted using a photoconductive antenna, is more commonly used today, a small discussion of how hazardous such systems are might be appropriate, and for comparison the average beam power of the system will be used.
As shown in Fig. 7, a power of 0.56 mW and a beam diameter of 0.5 mm is required to create a temperature increase of 1 • C, whereas the average power from a typical laser-based THz-TDS system is < 1 µW [13,20].It is furthermore interesting to note that from Planck's black body curve the normal ambient worst-case THz background power level at room temperature on a 0.5 mm diameter spot is 9.9 µW, which by itself indicates that these systems will barely affect the measurements by heating.The value is calculated by integrating the Planck's function [21] in the interval from 0.1 THz to 10 THz, and multiplying by the spot area.Where, c 0 is the speed of light in vacuum, h is Planck's constant, and k B is Boltzmann's constant.Thus the THz beam from photoconductive antennas will at maximum give a negligible temperature rise of 1.8 mK, and this is only if the water is held completely steady right at the focal point for a long time until the temperature reaches the steady-state value.This heating is barely noticeable and it can be argued that no external heat sink is needed when working with pulsed THz radiation from photoconductive antennas.
Though it is possible to amplify the terahertz signal up to an average power of 3 mW [20], here the temperature increase will be 5.3 K, but again one must keep the system in a steady position right in the focal point for long time to achieve this effect, if the sample is moved from the focal point to a place where the spot size is 10 times as large (a = 2.5 mm) the heating decreases with almost a factor 10 to 0.57 K.

Application to the human body
The basal layer of human skin contains 70 % of water [22], therefore our model may be exploited to give a very rough estimation of how a given THz source can affect the human body.However, note that all results are worst-case approximations, and that the temperature increase be much lower when the exposure time and the circulatory system under the skin surface (lymph, blood streams etc.) are properly taken into account.
Since the use of terahertz techniques within medical applications and other related areas is becoming commonplace [2,[23][24][25][26] it is still interesting to carry out the calculations to give an idea of how THz radiation might heat the human body in the worst-case.
Terahertz photons have significantly lower energies than do X-rays, and are unable to damage biological cells via ionization.Furthermore, in terms of photon energy terahertz radiation is noninvasive when considering biological samples, in spite of the nonlinear interaction between tissue and coherent terahertz radiation as verified by Fröhlich [27] and Grundler and Kaiser [28].This is evident from the low energy levels involved in the absorption process, when a THz beam is incident on a sample.For example, at a worse-case frequency, ν, of 10 THz, the molar energy calculated from E = N A hν, where N A is Avogadro's number and h is Planck's constant, equals 4 kJ/mol.In comparison, the ATP, ADP, and AMP processes that are fundamental to biological systems take up 31.8,31.8, and 14.2 kJ/mol, respectively [29].Thus terahertz photon energy is below the level for activating any biologically significant process.
However recent research by Alexaandrov et al. [30] has indicated that under certain conditions (high power, long exposure time and the right frequency) the THz field might be able unzip double-stranded DNA, creating bubbles in the double strand that could significantly interfere with processes such as gene expression and DNA replication.Note that these results are all numerical and have not yet been experimentally verified.In past years, several experiments have been performed investigating the influence of THz exposure.The conclusion is that high power exposure for a long time will in some way be hazardous for the human body.Bondar et al. [31] showed that, for a high power density of 15 kW/m 2 , which corresponds to an output power of 3.0 mW with a beam diameter of 0.5 mm, a 30 min exposure of a 3.6 THz beam to mice resulted in behavioral changes, but at an exposure time of 5 min no changes could be detected.The exact reason for these differences is not known, and it could be, as Alexaandrov et al. suggest, a result of nonlinear interference with the DNA, but it could also be caused by heating of tissue due to absorption of the THz radiation as discussed in this paper.However, in the worst-case, our results show the local temperature increase would at most be only 5 • C.
If we use water as a worst-case approximation to biotissue the result is that, for a THz beam with a diameter of 0.5 mm, 0.56 mW is required to create a temperature increase of 1 • C. Quantum cascade lasers can produce CW terahertz beams with output powers up to 17 mW [10], a long time exposure of such a beam will create a temperature increase of 30 • C, this means that a long time exposure could affect the human body -this motivates more accurate heating studies for predicting THz induced temperature increase for either therapeutic applications or for determining safety limits under specific conditions.However, note that accelerated electron sources of terahertz radiation, such as synchrotrons and free-electron lasers, have average powers as high as 20 W [11][12][13].At these levels THz radiation is clearly unsafe.

Accounting for the effect of transmittance
As per the usual procedure in physics, we start with a simple model that gives physical insight upon which further complexity can then be added as dictated by practical requirements.So far, throughout this paper, we have motivated the physical understanding of THz heating effects by considering the transmitted THz power as an input to the heating model.In order to then calculate the actual source THz power this corresponds to, in practice, requires us to now consider the transmittance, T given by [21]: Where nsample is the complex refractive index of water, which in terms of the complex permittivity reads [18]: We again use the complex permitivity formula by Ellison, which is valid from 0 THz to 25 THz and 0 • C to 100 • C [17].
In Section 5, our heating model thus far showed that for a 0.5 mm CW beam spot size at 1 THz, the steady-state temperature increase per milliwatt is 1.8K/mW in terms of transmitted power.Now by applying Eq. ( 15) and ( 16), at 1 THz, we obtain a transmittance T = 0.87, so we find that in terms of source power this becomes 1.6K/mW.This example, for a single interface, shows that whilst there is a small difference between transmitted and source power, the transmitted power may be used as a rough worst-case upper limit for the source power.Such an approximation errs on the side of underestimating maximum allowable source power, which provides a favorable safety margin.

Conclusion
In this paper we have developed the generalized model for estimating the heating effect, in terms of the steady-state temperature increase per milliwatt, in a sample exposed to terahertz radiation.We have exemplified this model using water as a topical case of interest in the community.From the model we learn that the heating effect is more intense at the sample's surface for a smaller beam width, and is more penetrative for a larger beam width.Furthermore, the temperature increase is roughly in direct proportion to the terahertz frequency and roughly inverse proportion to the initial temperature.An application of the model to terahertz spectroscopy shows that conventional pulsed THz-TDS causes a negligible rise in temperature in a water sample.However, it has been shown that the heating effect on a sample in a high-power CW terahertz setup is significant and could question the reliability of the results from the performed measurements if no care is taken.We have also highlighted that exposure to humans or animals can lead to heating of the body, and this motivates further study that builds upon this work.

Fig. 1 .
Fig. 1.Illustration of the cross-section of the disc of water, which is used in the calculations.Here, b and d are the radius and the thickness of the disc, respectively, a is the radius of the THz beam and [a], [b], [c] and [d] denote the chosen boundary conditions for the transformed temperature, U (r, z).The calculations are performed in cylindrical coordinates, which is why only half of the cross-section is needed by argument from symmetry.

Fig. 2 .Fig. 3 .
Fig.2.The temperature change per milliwatt, Ũ(r, z), (a) at z = 0 as a function of the radial distance, r, and (b) at r = 0 as a function of the thickness, z.It is seen that the maximal heating occurs just in the center of the beam and that the result has an exponential tail, but this may vary depending on the choice of dissipated power.

Fig. 4 .
Fig. 4. The dashed red line on both graphs indicates the chosen initial values.(a) The temperature change per milliwatt at (r, z) = (0, 0) as a function of the sample thickness, d.If the chosen thickness of the water disc is too small, the profile along the z-axis loses its exponential behavior, which is unphysical.Thus we have chosen a thickness of 15 mm for the sake of example.(b) The temperature change vs. the beam radius, a, at (r, z) = (0, 0).It is hard to focus a THz beam with a center frequency of 1 THz down to a spot radius smaller than 0.25 mm[20] due to the diffraction limit, and most often the radius will be around 0.5 mm.Note that this quantity is investigated even further in Fig.5.

aFig. 5 .Fig. 6 .
Fig. 5. (a)-(f) Contour plots, with the normalized radius, r/a, on the x-axis and the thickness, z on the y-axis showing the heating caused by the THz beam for six different beam radii a = 0.05 mm, a = 0.25 mm, a = 0.75 mm, a = 2.50 mm, a = 5.00 mm and a = 10.00 mm.The colors indicate the level of heating and the color-bar above each plot shows the conversion (in Kelvin per milliwatt), the vertical dashed lines show where the beam hits the water.(g)-(h) The normalized temperature change, Ũ(r, z)/ Ũmax , as a function of the normalized radius, r/a, and the thickness z, respectively.The six curves represent the six different beam radii as shown on the contour plots.

Fig. 7 .
Fig. 7. (a) Contour plot of the sample near the center, the colors indicate the temperature change per milliwatt, Ũ(r, z), and the vertical lines show where the terahertz beam hits the disc.(b) The power needed to create a temperature change of 1 K as a function of the beam radius.The red vertical dashed line indicates the chosen initial value for the beam radius and the horizontal line shows the power corresponding to this value (0.56 mW).Note that P and U are directly proportional P = U/ Ũ(r, z) , i.e with a beam radius of 0.25 mm the amount of transmitted power needed to heat the water for instance 2 • C is 2 × 0.56 mW = 1.2mW.The initial conditions given in Table 2 is used in both figures.

Table 2 .
Initial values for the calculations.