"Cooling by heating"- demonstrating the significance of the longitudinal specific heat

Heating a solid sphere at the surface induces mechanical stresses inside the sphere. If a finite amount of heat is supplied, the stresses gradually disappear as temperature becomes homogeneous throughout the sphere. We show that before this happens, there is a temporary lowering of pressure and density in the interior of the sphere, inducing a transient lowering of the temperature here. For ordinary solids this effect is small because c_p is almost equal to c_V. For fluent liquids the effect is negligible because their dynamic shear modulus vanishes. For a liquid at its glass transition, however, the effect is generally considerably larger than in solids. This paper presents analytical solutions of the relevant coupled thermoviscoelastic equations. In general, there is a difference between the isobaric specific heat, c_p, measured at constant isotropic pressure and the longitudinal specific heat, c_l, pertaining to mechanical boundary conditions that confine the associated expansion to be longitudinal. In the exact treatment of heat propagation the heat diffusion constant contains c_l rather than c_p. We show that the key parameter controlling the magnitude of the"cooling-by-heating"effect is the relative difference between these two specific heats. For a typical glass-forming liquid, when temperature at the surface is increased by 1 K, a lowering of the temperature in the sphere center of order 5 mK is expected if the experiment is performed at the glass transition. The cooling-by-heating effect is confirmed by measurements on a 19 mm diameter glucose sphere at the glass transition.


I. INTRODUCTION
Most solids and liquids expand when heated. Heat diffusion is a notoriously slow process, and heating a solid sample at its surface induces stresses in the sample that only disappear when temperature gradually becomes homogeneous again throughout. Heating a lightly fluent fluid that has a free surface (i.e., is free to expand), on the other hand, makes the entire sample expand on the time scale set by the sound velocity and sample dimensions. In this case, there are no transient stresses beyond the acoustic time scale. A liquid close to its glass transition provides an interesting case in between the solid and fluid behavior. Such a liquid behaves like a solid on time scales shorter than the Maxwell-relaxation time M ¼ =G 1 , where is the shear viscosity and G 1 the instantaneous shear modulus. The Maxwell-relaxation time becomes longer than 1 s when a liquid approaches its calorimetric glass transition, implying that induced stresses survive for seconds or more. This paper discusses the ''cooling-by-heating'' effect that arises when a sample is heated at a free surface. We show that this effect, which is present in all hard solids with a nonzero thermal-expansion coefficient, is generally magnified considerably for glass-forming liquids close to their glass-transition temperature T g . This is because close to T g the liquid is solidlike by having a large, nonzero dynamic shear modulus on short time scales and, at the same time, is liquidlike by having a fairly large thermal-expansion coefficient.
Returning to the case of a solid, what happens when heat is supplied at the (free) surface of a spherical sample? The outermost layers attempt to expand, obviously, but a priori one may imagine two different possibilities: (1) The expansion presses inwards, resulting in an increase of the pressure at the center of the sphere; or (2) the expansion turns outwards, thus transmitting a negative pressure into the sphere. Which of the two possibilities applies is answered by the application of standard thermoelasticity theory to the problem of calculating the stresses induced by the heating. This is done in the present paper. It turns out that case 2 applies-the sphere expands and pressure decreases in the interior of the sphere. This induces an adiabatic cooling inside the sphere. The phenomenon of cooling caused by heating at the surface is referred to below as the cooling-by-heating effect. The situation is not self-evident since if the shear modulus is small there is no change in pressure although there is an expansion. This is the case for a material like rubber.
The solution of the coupled thermomechanical equations detailed in Sec. III shows that the cooling-by-heating effect is proportional to the difference between the reciprocals of the isobaric specific heat c p and the longitudinal specific heat c l (all specific heats are per unit volume); the latter quantity is the heat needed to increase temperature by one unit if the associated expansion is confined to be longitudinal instead of isotropic [1,2]. The longitudinal specific heat is related to the isochoric specific heat c V by where M S and M T are the adiabatic and isothermal longitudinal moduli, respectively. This is analogous to the standard thermodynamic relation c p ¼ ðK S =K T Þc V relating the isobaric specific heat to the isochoric specific heats in terms of the adiabatic K S and isothermal K T bulk moduli.
where G is the shear modulus, one readily finds that c l is in between c V and c p . As we shall see, the relative difference a l ¼ ðc p À c l Þ=c p controls the strength of the cooling-by-heating effect, and we thus term this quantity the ''longitudinal thermomechanical coupling constant.'' Combining the equations above a l is found to be the product of two factors [1], The first factor shows that the cooling-by-heating effect is suppressed if the shear modulus is small compared to the longitudinal modulus. The other factor, ''the thermomechanical coupling,'' a ¼ ðc p À c V Þ=c p is a dimensionless measure of the coupling between thermal and mechanical perturbations. It can be expressed in terms of the expansivity, p ð1=VÞð@V=@TÞ p , as follows: where T 0 is the temperature. It follows that the coolingby-heating effect is quadratic in the thermal-expansion coefficient p . Since solids typically expand significantly less upon heating than do liquids, the cooling-by-heating effect is generally small in solids. As an example, for solid glucose the thermal-expansion coefficient [3] is 1:1 Â 10 À4 K À1 close to the glass transition whereas for liquid glucose it is 3:7 Â 10 À4 K À1 in the same temperature region. This potentially enhances the cooling-by-heating effect by a factor of 11 in the very viscous liquid compared to the solid glass. However, the change in c p [3] from 1:91 Â 10 6 JK À1 m À3 to 3:05 Â 10 6 JK À1 m À3 and in K T [4] from 10:75 Â 10 9 Pa to 6:49 Â 10 9 Pa reduces this to a factor of 8.
Here we have used a density of 1:52 Â 10 3 kg m À3 to convert specific-heat data from mass to volume. It is, however, not unusual for liquid expansivities to be near 10 À3 K À1 for which we would expect an enhancement of the thermomechanical coupling a ¼ ðc p À c V Þ=c p by a factor of 30. The shear modulus of glucose in the glassy state is G 1 ¼ 3:1 Â 10 9 Pa as deduced from the shearcompliance data of Meyer and Ferry [5].
The above relations all generalize to deal with cases of complex, frequency-dependent (dynamic) specific heats and moduli and expansivity, which are the relevant quantities when studying glass-forming liquids. Near the glass transition, the cooling-by-heating effect may be studied on time scales of seconds. Here, upon increasing the frequency, the factor G=M T in Eq. (2) increases while at the same time the factor ðc p À c V Þ=c p decreases. The enhancement of the cooling-by-heating effect thus critically depends on the relative time scales of the different relaxation processes at the glass transition. In the case that the shear modulus relaxes on a much faster time scale than the volume processes, the cooling-by-heating effect will be less pronounced. This situation is illustrated in Fig. 1. The model describing the relaxation behavior between highand low-frequency values is described in Sec. IV.
The present work discusses the basis of cooling by heating by referring to the equations of standard linear thermoviscoelasticity. Section II introduces the general framework of thermoelasticity and thermoviscoelasticity. It is shown that the heat-diffusion constant involves the longitudinal specific heat. Section III discusses the case when a finite amount of heat is fed into a sample at its surface at t ¼ 0, as well as the experimentally more easily realized case when temperature is suddenly increased at the surface. That section also presents analytical calculations of the ordinary solid case for which the constitutive properties do not undergo relaxation. Section IV gives calculations of a model glass-forming liquid, i.e., the case when the constitutive properties are frequency dependent. We estimate the effect to be of the order of 5 mK at The longitudinal coupling constant a l ¼ ðc l À c p Þ=c p (black line) is the product of those two frequency-dependent functions. Near the crossing of the two curves, the difference between c l and c p is generally largest. This determines the time scale of experiments on glass-forming liquids, where the cooling-by-heating effect is particularly large. the center of a sphere for a temperature increase at the sphere surface of 1 K. Section V confirms this prediction with measurements on a glucose sphere. Sections VI and VII briefly discuss and summarize the paper.

II. THERMOELASTICITY AND HEAT DIFFUSION
Thermoelasticity [6][7][8][9][10] deals with problems where displacement field uðr; tÞ and temperature field Tðr; tÞ in a material couple. It is a linear theory of small deformations given in terms of the strain tensor ij ¼ 1 2 ð @u i @x j þ @u j @x i Þ and small temperature increments T ¼ T À T 0 relative to a reference temperature T 0 . The material properties of a thermoelastic medium are given by a set of linear constitutive relations that express stress ij and increments in entropy density s in terms of ij and T. The hydrostatic pressure is related to the trace of the stress tensor p ¼ À1=3 P i ii and the relative compression is the trace of the strain tensor ¼ The following constitutive relations [11] define the shear modulus G, the isothermal bulk modulus K T , the isochoric specific heat c V , and the isochoric pressure coefficient V : We follow Biot [12] in assigning the symbol to the thermodynamic pressure coefficient, The material is furthermore characterized by its heat conductivity , which enters Fourier's law for the entropy current density j s : The interest in thermoelastic problems has, since Duhamel [13], mostly been focused on the calculation of thermal stresses deriving from an evolving temperature field. In the classical thermoelasticity theory the displacement and temperature fields are partially decoupled [6,7]. This comes from assuming that the development of the temperature can be found independently of the stresses by the conventional heat-diffusion equation: Here D is a heat-diffusion constant. After solving this equation the displacement field can be found from the quasistatic stress equilibrium equation: M T rðr Á uÞ À Gr Â ðr Â uÞ À V rT ¼ 0: (10) This approximate theory is referred to as the theory of thermal stresses [6]. According to many authors [6,[8][9][10], the correct treatment appeared remarkably late in the development of thermoelastic theory with Biot's paper [12] in 1956. Lessen [14] considered similar problems the same year. The heat-diffusion equation, Eq. (9), is replaced by which follows from entropy conservation, when this is combined with Eqs. (6) and (8). Entropy conservation may seem strange at first sight, but the entropy production per volume associated with heat conduction is Àj s Á 1 2 , i.e., a second-order effect disappearing in a linearized theory. In most cases the ordinary, decoupled heat-diffusion equation is a good approximation in the manner it is used in the theory of thermal stresses. However, this approximate theory is not able to describe the phenomenon of cooling by heating, which is the theme of this paper. It should be noted that the heat-diffusion equation, with the diffusion constant containing the isobaric specific heat c p , is exact for the nonviscous liquid state or soft matter with G ¼ 0, if part of the boundary (with normal vector n) is free to expand, i.e., P j ij n j ¼ 0. The proof runs as follows: The assumption G=M T ¼ 0 simplifies Eq. (10) to rðK T r Á u À V TÞ ¼ 0: However, the term under the gradient is according to Eq. (5), nothing but minus the pressure increment. Thus, we conclude that this pressure increment is uniform in space and only depends on time. Moreover, Eq. (4) ensures-when G ¼ 0-that all diagonal elements of the stress tensor are identical and equal to minus this pressure increment. If the normal component P j ij n j is zero on part of the boundary, it follows that the pressure is also zero there, but then it is zero throughout the body. Equation (5) then reduces to r Á u ¼ p T. Inserting this in the entropy equation, Eq. (11), one arrives at the ordinary decoupled heat-diffusion equation with D p ¼ =c p , when noticing that As we have seen, the temperature field, in general, does not exactly obey a diffusion equation. It does so when c p ¼ c V (, V ¼ 0) or for certain boundary conditions when G ¼ 0. However, as emphasized by Biot [12], the entropy density in fact does fulfill a diffusion equation and, moreover, with a diffusion constant containing the ubiquitous longitudinal specific heat: Applying the divergence operator to the inertia-free stress equilibrium equation, Eq. (10), gives Applying the Laplacian to the constitutive relation, (6), yields Fourier's law and the entropy conservation equation, (12), gives and thus @s @t ¼ c l r 2 s; with c l ¼ c V þ T 0 2 V M T being the longitudinal specific heat [1,2].
Note that this result is limited to the inertia-free cases. If one wishes to study coupled mechanical and thermal waves, the inertia-term @ 2 @t 2 u must be added on the right side of Eq. (10). Solutions of the equations in this case have been studied extensively [6] also in the spherically symmetric situation. Note, however, that acoustic wavelengths are much longer than thermal wavelengths. Thus, for a sample of a certain size there is an interesting time regime where acoustic waves have settled, but thermal diffusion has barely begun. Take as an example a sphere of radius 1 cm. For a typical sound velocity of 10 3 m=s and heatdiffusion constant of 10 À7 m 2 =s, the sound traveling time is 10 s while the diffusion time is 1 ks. It is within this time regime that we will find the cooling-by-heating phenomenon. Although the solution restricted to the inertiafree case that we present below is, in principle, contained in the coupled acoustic-thermal wave solutions including inertia, the phenomenon is obscured by the complicated structure of these solutions and seems not to have been recognized.
The thermoelastic theory that was originally developed for elastic solids without relaxation is easily extended to a thermoviscoelastic theory taking relaxation of all the constitutive parameters into account, as it is necessary for relaxing liquids near the glass transition. The most straightforward way of generalizing is to interpret the equations in the frequency domain, allowing all the constitutive parameters to be complex functions of the angular frequency !. The cases we study in the frequency domain cover thus both solids and thermoviscoelastic liquids, but can only be transformed analytically into the time domain for solids.
For relaxing liquids, one must do the transformation numerically.

III. ANALYTICAL SOLUTIONS OF THE SPHERE-HEATING PROBLEM
A. The case when the heat flow is controlled at a mechanically free boundary This subsection presents the analytical solution in the frequency domain to the situation when a sphere of a general viscoelastic material is subjected to a periodic heat input at the surface. The solution shows the temperature at the center varying, at high frequencies, 180 out of phase with the heat oscillation at the surface, indicating the cooling-by-heating effect. In order to give a more lucid and transparent understanding of the phenomenon, we translate the solution to the time domain. This can be done analytically by an inverse Laplace transformation if there is no frequency dependence of the constitutive properties. That is, we calculate the temperature and stress profile throughout the sphere following a stepped heat input at the surface at time zero.
Consider the case when a periodically varying heat QðtÞ ¼ RefQe i!t g is supplied at the surface of a sphere of radius R. The surface is assumed to be mechanically non-clamped, i.e., the sphere is free to expand. We wish to calculate how the periodically varying temperature and displacement fields vary throughout the sphere, i.e., to calculate the complex frequency-dependent amplitudes of temperature, Tðr; !Þ, and radial displacement field, uðr; !Þ. From these quantities the stress components rr ðr; !Þ, etc., may be calculated. Denoting the position vector by r and the complex frequency-dependent radial displacement field by uðr; !Þ ¼ uðr; !Þr=r, the coupled thermoelastic equations, (10) and (11), become @ @r M T r À2 @ @r The four boundary conditions are (1) no displacement at the center: uð0; !Þ ¼ 0; (2) no temperature gradient at the center: @T @r ð0; !Þ ¼ 0; (3) free surface, i.e., no radial stresses at the surface: rr ðR; !Þ ¼ 0; (4) heat supply boundary condition at the surface: @T @r ðR; !Þ ¼ i! Q 4R 2 . Denote the volume of the sphere by V 0 ¼ 4 3 R 3 and define the complex frequency-dependent thermal wave vector by k ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi i!c l ð!Þ= p . Furthermore, define the functions, sinhðkrÞ kr ðkRÞ 3 kR coshðkRÞ À sinhðkRÞ ; (21) ðkrÞ coshðkrÞ À sinhðkrÞ ðkRÞ coshðkRÞ À sinhðkRÞ : Introduce the characteristic heat-diffusion time The solutions to Eqs. (19) and (20) are These solutions were found by the transfer-matrix approach (see Ref. [2] and Appendix A 1), and can be verified by insertion, noticing that f 1 ð; sÞ ¼ Consider the low-and high-frequency limits of these expressions. The functions f 1 and f 2 both have the limits, 1 for ! ! 0 and 0 for ! ! 1. Thus, as expected, T ! Q=ðV 0 c p Þ in the low-frequency limit when heat has had time to distribute throughout the sphere. In the high-frequency limit, ! ! 1, the temperature amplitude becomes If we, for a moment, consider the no-relaxing case where the specific heats are real, we see that the temperature amplitude is in counterphase with the heat amplitude since c l < c p . For a propagating thermal wave, it would not be surprising that temperatures at points of some distance apart-e.g., at a half wavelength-had opposite phases. However, Eq. (26) holds throughout the sphere and is not associated with the diffusive thermal wave. This will be even more clear when we consider the response to a stepped heat input later on. We see that the longitudinal coupling constant a l controls the magnitude of the cooling-by-heating effect. The ratio of the amplitudes of the temperature at the center and the heat input at the surface is shown in Fig. 2. The phenomenon, cooling by heating, is seen to occur at high frequencies, albeit this is more conspicuous in the time domain.
For the displacement field, we find in the low-frequency limit the natural result, determined by the final temperature rise Q V 0 c p and the linear thermal expansion 1 3 p . However, at high frequencies we find Notice that this displacement, which is responsible for the cooling-by-heating effect, is only present when G Þ 0.
In the spherically symmetric case there are only two different components of the stress tensor, rr and ¼ . It follows from Eqs. (4) and (5) and the fact that rr ¼ @u @r and ¼ ¼ u r that rr ¼ ðK T þ 4 3 GÞ Â @u @r þ 2ðK T À 2 3 GÞ u r À V T, which by Eqs. (24) and (25) becomes Thus, at high frequencies there is an isotropic, uniform tensile stress in the interior of the sphere of the magnitude: On the other hand, all stresses vanish for ! ! 0 (as expected). The real part of the ratio between the complex amplitudes of temperature at the center and the heat supplied at the surface scaled with the isobaric heat capacity. At high frequencies the limit becomes the negative value Àa l =ð1 À a l Þ ¼ 1 À c p =c l .
In order to gain insight into the cooling-by-heating effect and show the significance of the longitudinal coupling constant a l , we transform Eqs. (24), (29), and (30) into the time domain, but only for a solid, i.e., in the case when all constitutive properties are frequency independent. If a delta-function heat flux is applied at t ¼ 0, the heat supplied at the surface is a Heaviside step function, QðR; tÞ ¼ Q 0 HðtÞ; in this case, calculating the inverse Laplace-Stieltjes transform leads to the following expressions for the temperature and stresses as functions of time after t ¼ 0 (see Appendix A 2): Here x 1 < x 2 < . . . are the positive roots of the transcendental equation x ¼ tanðxÞ. Note that F 1 and F 2 ! 0 for t ! 0 and F 1 and F 2 ! 1 for t ! 1. When a l ¼ 0, there is no cooling-by-heating effect according to Eq. (32). Furthermore, a l ¼ 0 implies either G ¼ 0 or V ¼ 0, and there is no immediate expansion and no induced stresses.
When a l Þ 0, the situation is quite different. In Fig. 3 we plot the scaled temperature change ðc p V 0 =QÞTðt=; r=RÞ for several radii r=R as given in Eq. (32). The longitudinal coupling constant is here fixed to a l ¼ 0:091 and time is given in units of the characteristic heat-diffusion time . The figure clearly shows the cooling-by-heating effect. Since a finite amount of heat was added at the surface at t ¼ 0, the surface temperature initially diverges. The interior of the sphere, even close to the surface, instantaneously cools to a uniform temperature. The expansion of the surface is immediately felt in the interior, and since no heat has yet arrived by diffusion, it cools adiabatically. This initial response is followed by an evolution in time where the temperatures of the different parts of the sphere converge and eventually equilibrate.
In order to understand better the physics of cooling by heating, we consider the components of the stresses given by Eqs. (33) and (34), respectively. In Fig. 4 the rr component of the stress tensor is plotted scaled with the initial uniform interior stress 0 ¼ 4 First, we note that the boundary condition rr ðR; tÞ ¼ 0 is fulfilled. As the surface receives heat and expands, an immediate traction is felt in the interior of the sphere. rr is positive, seeking to stretch a volume element in the radial direction under the entire evolution to thermal equilibrium.
The scaled stress component ðr; tÞ is shown in Fig. 5. One notices an immediate, uniform increase of this stress component throughout the sphere of the same size as rr . The initial stress is thus isotropic. Note that shifts sign during the thermal equilibration process, in contrast to rr . This can be understood in a physical picture: Consider the outer region that has been reached by the inflowing heat at a certain point in time. If that region were free it would expand, but it is kept in place by the inner unheated region that has not expanded thermally yet. This creates a negative stress on surfaces with a normal vector at right angles to the radius vector.
We conclude this section by a simple result. If one compares the instantaneous temperature drop, Eq. (26), and the instantaneous stress increase, Eq. (31), one finds that the ratio is given by the adiabatic temperature-pressure coefficient: B. The case when temperature is controlled at a mechanically free boundary The above studied case with heat-input control showed a rather simple cooling-by-heating behavior at short times or high frequencies. We now consider the case of controlling the temperature on the outer surface instead. There is still an effect, but it is not instantaneous. We only calculate the temperature at the center of the sphere. The surface is again mechanically free. Again, using the transfer-matrix technique in the frequency domain, one finds that the temperature amplitude Tð0; sÞ at the center is related to the temperature amplitude TðR; sÞ at the surface by Tð0; sÞ ¼ ÈðsÞTðR; sÞ, where ÈðsÞ ¼ 1 À x 3 À x 2 sinhðxÞ 3a l ½x coshðxÞ À sinhðxÞ À x 2 sinhðxÞ : Here x ¼ ffiffiffiffiffiffiffiffiffiffiffiffi ffi sð!Þ p , s ¼ i!. The characteristic diffusion time ð!Þ [Eq. (23)] and thermomechanical coupling constant a l ð!Þ [Eq. (2)] are, in the general thermoviscoelastic case, complex and frequency dependent and the temperature response can only be converted to the time domain numerically.
In order to calculate the temperature signal as a function of time, we again limit ourselves to the purely thermoelastic case, i.e., the case of a solid where and a l are real and frequency independent. For a Heaviside temperature step at the surface of the sphere, TðR; tÞ ¼ ÁTHðtÞ, the temperature at the sphere center is calculated via an inverse Laplace-Stieltjes transform of ÈðsÞ, where the residues are given by Here the y k 's denote the positive roots of the transcendental equation: In Fig. 6 we plot the solution, Eq. (39), for various values for the coupling constant a l . Time is given in units of the characteristic diffusion time and ÁT ¼ 1 K. We see that the cooling-by-heating effect is present also when a step in temperature (instead of heat) is applied to the surface. However, now the effect is not instantaneous but evolves gradually, reflecting the gradual heat diffusion at the surface mediated to the center by the stress field. Figure 6   furthermore shows that it is not enough to have a thermomechanical coupling ( p Þ 0) for the phenomenon to be present-only when c l Þ c p is there a cooling-by-heating effect. The next section studies the general, thermoviscoelastic case of the frequency dependence of the response functions, which describes supercooled liquids.

IV. THE THERMOVISCOELASTIC CASE
The above time-domain results apply for a thermoelastic solid only, whereas the frequency-domain results are general. The thermoelastic examples handled so far only involved frequency-independent constitutive parameters, corresponding to the high-frequency (low-temperature) limiting values of the curves sketched in Fig. 1. However, Fig. 1 indicates that the value of the coupling constant a l is larger at lower frequencies (higher temperature), thus suggesting that the effect of cooling by heating may be even larger in the very viscous liquid or simply at the glass transition. We investigate this issue in the time domain in this section. In the thermoelastic case the inversion of the problem to the time domain could be made analytically. This is not possible in the thermoviscoelastic case where the constitutive parameters are complex frequency-dependent functions.
In order to investigate the effect of going from solid to liquid we resort to numerical methods. Specifically, we transform ÈðsÞ of Eq. (38) into the time domain, accounting for the frequency dependence of the constitutive parameters via and a l . To do this we have to introduce a model of the constitutive parameters that enters via and a l . It is common in rheology to illustrate models like the Maxwell model by rheological networks or even their electrical analogue. We use this approach to model the thermoviscoelastic behavior. The purpose of the model is to interpolate between the thermodynamic coefficients at high frequencies, T;1 ,c p;1 , p;1 , and at low frequencies T;0 , c p;0 , p;0 . Network modeling assures internal consistency and agreement with the rules of linear irreversible thermodynamics. A one-parameter relaxation model implies that the Prigogine-Defay ratio is unity, which is not the case for glucose. Rather, with T g ¼ 300 K and Ác p ¼ 1:14 Â 10 6 JK À1 m À3 , Á T ¼ 6:1 Â 10 À11 Pa À1 and Á p ¼ 2:6 Â 10 À4 K À1 one finds We are thus forced to consider a model with two relaxation elements that cannot be lumped into one. The model of Fig. 7 is suited for this purpose. In order to still make it simple, the two relaxation elements are taken to be Debye-like. The model has a simple mathematical formulation in the frequency domain. We change the independent variables compared to Eqs. (5) and (6) and consider the complex amplitudes T and p to be the controlled stimuli creating a linear response in the amplitudes S and of the entropy density and dilation: In the model, the three measurable quantities, the isobaric specific heat, the isobaric expansivity, and the isothermal compressibility are related to the elements D, C, J A ð!Þ, and J B ð!Þ by The relaxation element J A is determined by three parameters J A;1 , ÁJ A ¼ J A;0 À J A;1 , and R A : and likewise for J B .
The parameters of the model can be established from the high-and low-frequency limits of p , T , and p . One finds D ¼ ÀÁ p =Á p , J A;0 ¼ À p;0 =D, J A;1 ¼ À p;1 =D, J B;0 ¼ T;0 À J A;0 , and J B;1 ¼ T;1 À J A;1 . The Prigogine-Defay ratio of the model is given by and the dynamic (frequency-dependent) Prigogine-Defay ratio [15] is given by 7. Electrical equivalent diagram of the interaction of a volume element with its surroundings through two gates. The thermal gate where entropy displacement S or temperature T can be controlled, and the mechanical gate where volume displacement V or pressure p can be controlled. S and V are generalized charge displacements and T and Àp are generalized voltages. The relaxational elements J A and J B are simple single-relaxation-time elements.
As expected, the Prigogine-Defay ratio is larger than unity, but becomes one if the element J B is nonrelaxing, in which case the model reduces to a single-parameter model. In the case of J A being non-relaxing there is no relaxation of p and p and a Prigine-Defay ratio cannot be defined. For glucose at 300 K, we calculate the parameters of the model to be D ¼ À1:46 Â 10 7 Pa K À1 , C ¼ 4:76 Â 10 3 JK À1 m À3 , J A;1 ¼ 7:53 Â 10 À12 Pa À1 , ÁJ A ¼ 1:78 Â 10 À11 Pa À1 , J B;1 ¼8:55Â10 À11 Pa À1 , and ÁJ B ¼ 4:33 Â 10 À11 Pa À1 .
The heat conductivity of glucose-needed to calculate the heat-diffusion time-is ¼ 0:35 WK À1 m À1 at 303 K according to Greene and Parks [16]. The shear relaxation of glucose is modeled by a Maxwell model Here the high-frequency shear modulus can be taken to be G 1 ¼ 3:1 Â 10 9 Pa [5]. The values of the thermodynamic parameters used to parametrize the model are given in Table I. The temperature dependence of the shear viscosity causes the shift in the loss peaks of the relaxations. Parks et al. [17] measured the viscosity of glucose in a wide temperature range from 295 to 418 K. We fitted their tabulated data by the expression ðTÞ ¼ 0:0125 exp½ð512:9 K=TÞ 6:42 Pa s. This holds within 20% over the entire temperature range except for the highest viscosity point of 9:1 Â 10 12 Pa s, which, however, is well beyond the glass transition and may be hard to measure reliably. A Vogel-Fulcher law is not as good a fit, deviating more than 40% in the measured temperature range. The rate parameters R A ¼ R A ðTÞ and R B ¼ R B ðTÞ are assumed to follow the temperature dependence of the viscosity [18]. It is found numerically that one should choose R A ðTÞ ¼ R B ðTÞ ¼ 30ðTÞ in order to get the loss-peak frequency of the shear modulus Gð!Þ and isothermal bulk modulus K T ð!Þ to coincide. The relaxation of the different response functions described by the model implies a frequency dependence of the longitudinal thermomechanical coupling via The modulus of this complex function was shown in Fig. 1. Also, the heat-diffusion time, now becomes complex. This makes the inversion to the time domain nontrivial and thus Eq. (38) was inverted numerically. The algorithm for the inverse Laplace transform is an improved version of de Hoog's quotient difference method [19] developed and implemented in Matlab by Hollenbeck [20]. The calculated temperature response at the center of the sphere to a step of 1 K at the surface is shown in Fig. 8. Time is now scaled by the fixed real-valued diffusion time 0 in the liquid regime, The figure shows that the effect of the thermomechanical coupling is absent at high temperatures. But as temperature is decreased and the liquid gets more and more viscous, a dip in temperature emerges. Going further down in temperature, the phenomenon of cooling by heating becomes most pronounced slightly above T g . Even further down in temperature, in the glassy state, the effect is still present, but small. One may ask what happens if the expansivity p;1 vanishes, so that the phenomenon is TABLE I. Literature data for glucose (at 300 K) used to parametrize the model depicted in the electrical equivalent diagram in Fig. 7. We have used a density of 1:52 Â 10 3 kg m À3 to convert specific-heat data from mass to volume.  FIG. 8. The change in temperature T at the center of a sphere, after a temperature step of ÁT ¼ 1K has been applied at the surface. A model of the thermoviscoelastic relaxation of glucose with realistic limiting thermodynamic parameters has been invoked. Time is scaled by the characteristic diffusion time 0 , which is 800 s for a glucose ball of radius 9.5 mm. The minimum of approximately À5 mK occurs for the 309 K curve just above T g at 0:017 0 , corresponding to 14 s. Notice that although the cooling-by-heating phenomenon is present in the glassy solid state, it becomes more pronounced at the glass transition. The squares marking the minima at each temperature are plotted against temperature in Fig. 10. absent in the glassy state: Will it still be present at the glass transition? The simulations shown in Fig. 9 confirm this expectation. Putting p;1 ¼ 0, but otherwise keeping the values of the rest of the parameters, we get a succession of temperature evolutions. Going down in temperature we see the cooling-by-heating phenomenon appearing at T g and afterward disappearing in the glassy state. Figure 10 shows the minimum temperature T MIN reached when p;1 Þ 0 as a function of temperature T 0 , emphasizing the phenomenon as characteristic of the glass transition.

V. EXPERIMENTAL VERIFICATION OF THE COOLING-BY-HEATING EFFECT
To prove the existence of cooling by heating, we mold glucose [-DðþÞ glucose, 98%, Sigma-Aldrich] into spherical samples with a thermistor placed at the center. Via the large negative-temperature coefficient (NTC) thermistor, we measure the temperature in the middle of the sphere. During measurements, the glucose sphere is placed in a cryostat, which makes it possible to change the temperature at the surface of the sphere quickly compared to the characteristic heat-diffusion time. A sketch of the setup together with a photo of one of the samples is shown in Fig. 11. The photo shows the wires that lead into the sphere, connecting the thermistor to the terminals on the peek plate shown in the photo. When mounted on a holder, the terminals get connected to the multimeter that does the resistance measurement.
The procedure in the experiments is the following. First, we bring down the temperature to the desired starting level. Then we wait for the temperature to equilibrate. This is monitored by measuring the resistance every fifth minute. Typical waiting time is 18 h. After the initial waiting time, we increase the sampling rate of the multimeter to about five data points per minute over a period of 1 h to get a baseline like the one shown in Fig. 12. Then we impose a 5-K temperature step and continue sampling data for another ten minutes with a sampling rate of 15 data points per minute. Figure 12 shows the temperature measured by the NTC thermistor during a measurement with a step from 298 to 303 K of the cryostat temperature. The baseline extends for about 20 min, and then a characteristic temperature dip appears. The magnitude of the dip is 7:3 AE 0:2 mK, which was reached 40 s after the temperature step was imposed.
The experiment has been repeated on three different samples; the inset in Fig. 12 shows the results of the first At high temperatures, the phenomenon of cooling by heating is absent. Going down in temperature, the liquid becomes more viscous and a dip in temperature emerges. Close to the glass transition, the cooling-by-heating phenomenon gets most pronounced. As temperature is decreased further and the liquid enters the glassy state, the effect is still present, but is reduced in strength.
measurement done on each sample, at the same temperature. Each marker represents the lowest temperature reached in one measurement. They are plotted against time after the temperature step is initiated. It has not been possible to reproduce the phenomenon on the same sample by recycling the temperature. We ascribe this to crystallization since the opacity of the spheres appeared to change as the experiment was lengthened. Although the cooling-by-heating effect should be present also in the solid state, it is here considerably smaller and not observable with our temperature resolution. Nevertheless, the phenomenon has been seen every time we repeat the experiment with a fresh supercooled sample. The sphere does not flow or deform to any appreciable degree even somewhat above the glass transition. This is supported by the following estimate. The characteristic flow time flow is proportional to the viscosity and inversely proportional to the gravitational force mg. By a dimensional argument, it follows that which means that flow $ 10 7 M . Thus, even at 310 K, where viscosity becomes 10 10 Pa and thereby the Maxwell time, 3 s, the flow time is 1 yr.

VI. DISCUSSION
The concept of a longitudinal specific heat has been identified in Ref. [1] as the relevant quantity within AC-calorimetric methods that utilizes heat effusion. The principle of the simplest of these techniques [21] is to measure the complex temperature response T ! at a planar surface to a heat-current density j Q;! generated at the same surface. The effusivity e ¼ ffiffiffiffiffiffi c p is found from the measured specific thermal impedance Z T ! =j Q;! ¼ 1= ffiffiffiffiffiffiffiffiffiffiffi ffi i!c p , and from the effusivity the specific heat can be calculated. A thorough analysis of the thermomechanical equations of this problem shows that the specific heat that comes into play in this situation is c l rather than c p . Effusivity measurements in spherical geometry have been shown also to involve the longitudinal specific heat [2,22]. The longitudinal specific heat is not a very well-known property, but it does appear in the textbook on elasticity by Landau and Lifshitz [23]. They show that the coupled thermoelastic equations decouple for certain boundary conditions of an infinite solid, namely, when temperature at infinity is constant and deformation there is zero. They showed further that the heat-diffusion equation is valid with a diffusion constant containing the effective specific heat ½ð1 þ Þc p þ 2ð1 À 2Þc V =½3ð1 À Þ, where is the isothermal Poisson ratio. Inserting ¼ ð3K T À 2GÞ=ð6K T þ 2GÞ, one readily finds that the effective specific heat is c l . The longitudinal specific heat also appeared in Biot's 1956 paper [12] in his diffusion equation for the entropy density. Although not very different from c p , there is a fundamental difference, and c l appears in many thermoelastic problems when they are treated exactly. In particular, as we have seen in this paper, there is only a cooling-by-heating effect if c l Þ c p . We originally proposed the name longitudinal specific heat because this is the heat needed to increase the temperature by 1 K if the associated expansion is confined to be longitudinal instead of isotropic.
Thermal relaxation and mechanical relaxation of supercooled liquids are mostly treated separately although they obviously are connected appearing as they are at more or less the same time scale. Recently, theory and simulation [24][25][26] have revealed a class of liquids, ''the strongly correlating liquids,'' that have strong correlations between equilibrium fluctuations of energy and pressure. The relaxation of thermal and mechanical properties of these liquids can be described in terms of a single relaxing parameter. This means that certain triples [15] of response functions as, e.g., T ð!Þ, c p ð!Þ, and p ð!Þ, are linearly related. This can also be formulated as whether the socalled Prigogine-Defay ratio is 1 or larger than 1. It has been conjectured that van-der-Waals-bonded liquids belong to this class but hydrogen bonded liquids do not. The temperature and density dependence of the relaxation time for many liquids can be described in terms of the single scaling variable =T. A striking prediction for strongly correlating liquids is that the density-scaling exponent can be calculated from the relaxation strengths of the response functions [27]. These developments call for accurate measurements of thermoviscoelastic response functions and a thorough treatment of the influence of mechanical boundary conditions in thermal experiments and thermal boundary conditions in mechanical experiments. The cooling-by-heating phenomenon evidently exemplifies this point. For convenience, we have chosen glucose to show the effect experimentally. This liquid has a Prigogine-Defay ratio significantly greater than one and thus we have had to describe its relaxation by a twoparameter model. However the cooling-by-heating phenomenon will also be present in a single parameter liquid so this is not a distinguishing feature between single parameter liquids and more complex liquids. More important is the issue of the relative time scales [18] of the shearmodulus relaxation compared to those of the scalar thermoelastic response functions that make up the thermoelastic coupling a of Eq. (3). If G relaxes on a shorter time scale than a-that is, if the green curve of Fig. 1 is moved to the right-then the maximum of the longitudinal thermomechanical coupling a l (black curve) becomes smaller, reducing the magnitude of the cooling-by-heating effect. Transient thermal stresses induced by surface heating of a sphere have been considered theoretically by Cheung et al. [28]. Their interest was fragmentation of brittle solids by surface heating. A heat current was applied uniformly within the polar angle regime 0 0 and the temperature and stress distributions calculated in time and space. This is a problem very similar to the one considered in this paper, although not generalized to situations with relaxation. However, the cooling-by-heating phenomenon was not seen since the standard decoupled heat-diffusion equation was used. The phenomenon thus seems apparently not to have been recognized in the literature, even though it belongs to classical continuum physics.

VII. SUMMARY
We have shown that cooling by heating occurs at the center of a solid spherical sample if it is heated at a mechanically free surface, reflecting a nontrivial thermomechanical coupling where the temperature initially decreases in the interior of the sphere. What happens is that, as heat diffuses into the outermost parts of the sphere, these parts expand and build up a negative pressure at the center of the sphere. This negative pressure couples to the temperature via the adiabatic pressure coefficient ð @T @p Þ S . The opposite effect also applies, of course: If the temperature of the surface is lowered, heating by cooling will be observed. In ordinary solids, the cooling-by-heating effect is almost negligible because their thermal expansion is generally small. The effect is particularly large for liquids close to their glass transition. The cooling-by-heating phenomenon establishes the difference between the longitudinal and isobaric specific heat, since the effect is only present when these two quantities differ. This is the case when the shear modulus is not vanishingly small compared to the bulk modulus and, simultaneously, the isobaric and isochoric specific heats differ significantly. Analytical results show that the phenomenon occurs in the elastic case (the solid), and numerical results based on a model of the glass transition with parameters determined by glucose data show that the effect is present dynamically also in the very viscous liquid. Even in the hypothetical case when the glassy state is assumed to have zero expansivity, the phenomenon will still appear at the glass transition.
The numerical results based on glucose data indicate that the drop in temperature at the center of the sphere is of the order of 5 mK with a duration of approximately 15 s when the temperature is increased by 1 K on the surface of a sphere of radius r ¼ 10 mm. This prediction has been confirmed experimentally.

ACKNOWLEDGMENTS
This work grew out of a suggestion by Niels Boye Olsen to whom we are obviously indebted. The center for viscous liquid dynamics ''Glass and Time'' is sponsored by the Danish National Research Foundation.

Solution from Sec. III: The transfer-matrix method
The solution to the problem given in Sec. III is based on the transfer-matrix formulation [1,2,29] of the general solution of the thermoviscoelastic problem in a spherically symmetric case. Including the radial stress field and the time-integrated heat-current density that relate to the temperature and displacement fields, the authors of Refs. [1,2] end up with four coupled equations to solve. Laplace transforming the equations relating the four fields and solving the resulting inhomogeneous system of four ordinary differential equations, the result is in the general form of a transfer matrixTðj; iÞ that links the dimensionless complex amplitudes of the fields at the boundary r i with those at r j : (A1) Here S, Ṽ, T, and p r are the complex amplitudes of entropy, volume, temperature, and the radial component of pressure (p r ¼ À rr ), respectively. The elements of the transfer matrix are given in Ref. [2]. From this general solution, one can work out different cases, like the ones in Sec. III. The boundary condition atr 1 ¼ 0 is no net flux of heat through the center of the sphere, i.e., S 1 ¼ 0 and no volume displacement Ṽ 1 ¼ 0. At the mechanically free outer boundaryr 3 the entropy supplied S 3 is given and p r;3 ¼ 0. Lettingr ¼r 2 be an intermediate variable radius betweenr 1 andr 3 one has