Solution of the spherically symmetric linear thermoviscoelastic problem in the inertia-free limit

The coupling between mechanical and thermal properties due to thermal expansion complicates the problem of measuring frequency-dependent thermoviscoelastic properties, in particular for highly viscous liquids. A simplification arises if there is spherical symmetry where - as detailed in the present paper - the thermoviscoelastic problem may be solved analytically in the inertia-free limit, i.e., the limit where the sample is much smaller than the wavelength of sound waves at the frequencies of interest. As for the one-dimensional thermoviscoelastic problem [Christensen et al., Phys. Rev. E 75, 041502 (2007)], the solution is conveniently formulated in terms of the so-called transfer matrix, which directly links to the boundary conditions that can be experimentally controlled. Once the transfer matrix has been calculated, it is fairly easy to deduce the equations describing various experimentally relevant special cases (boundary conditions that are adiabatic, isothermal, isochoric, etc.). In most situations the relevant frequency-dependent specific heat is the longitudinal specific heat, a quantity that is in between the isochoric and isobaric frequency-dependent specific heats.


Solution of the spherically symmetric linear thermoviscoelastic problem in the inertia-free limit
The coupling between mechanical and thermal properties due to thermal expansion complicates the problem of measuring frequency-dependent thermoviscoelastic properties, in particular for highly viscous liquids. A simplification arises if there is spherical symmetry where-as detailed in the present paper-the thermoviscoelastic problem may be solved analytically in the inertia-free limit, i.e., the limit where the sample is much smaller than the wavelength of sound waves at the frequencies of interest. As for the one-dimensional thermoviscoelastic problem ͓Christensen et al., Phys. Rev. E 75, 041502 ͑2007͔͒, the solution is conveniently formulated in terms of the so-called transfer matrix, which directly links to the boundary conditions that can be experimentally controlled. Once the transfer matrix has been calculated, it is fairly easy to deduce the equations describing various experimentally relevant special cases ͑boundary conditions that are adiabatic, isothermal, isochoric, etc.͒. In most situations the relevant frequency-dependent specific heat is the longitudinal specific heat, a quantity that is in between the isochoric and isobaric frequency-dependent specific heats. DOI: 10.1103/PhysRevE.78.021501 PACS number͑s͒: 64.70.PϪ

I. INTRODUCTION
Linear thermoviscoelasticity is the well-established discipline dealing with the irreversible thermodynamics of slightly perturbed systems where mechanical and thermodynamic degrees of freedom couple to each other ͓1-4͔. For glass-forming liquids cooled toward the calorimetric glass transition, relaxation times become very long compared to phonon times, approaching and eventually exceeding seconds. For such "ultraviscous" liquids, not only the mechanical moduli become complex and frequency dependent ͓5-7͔, but so do standard thermodynamic linear-response properties like the specific heat or the thermal expansion coefficient ͓8-26͔.
Of the 12 basic complex, frequency-dependent thermodynamic linear-response coefficients only three are independent ͑see, e.g., Ref. ͓25͔ and its references͒. If one assumes stochastic dynamics, which is believed to be realistic for viscous liquids on time scales much longer than phonon times, in fact only two thermoviscoelastic response functions are truly independent ͓1,8,27,28͔. No reliable measurements of a full set of ͑three͒ thermoviscoelastic response function appear to exist for any highly viscous liquid. Part of the reason for this may be the traditional focus in physics on phenomena on the molecular scale, but part of the problem most likely also comes from the considerable challenges associated with reliably measuring the thermoviscoelastic response functions. The problem is that the coupling between mechanics and thermodynamics caused by thermal expansion is nontrivial when the mechanical shear modulus is a significant fraction of the bulk modulus ͓29͔. This is true for solids in general, as well as for highly viscous liquids at frequencies of order the inverse relaxation time. For solids, however, because the isobaric and isochoric specific heats are almost identical, the problem is not serious. For ultraviscous liquids, on the other hand, the problem cannot be ignored. For such systems it was recently shown that the thermomechanical coupling implies that conventional methods fail to measure the isobaric, frequency-dependent specific heat ͓23͔. This is because truly isobaric conditions are difficult to establish experimentally due to the thermal expansion upon heating; in most experiments attempting to measure the frequencydependent specific heat the stress tensor is not proportional to the unit tensor ͑i.e., there is not hydrostatic conditions͒, and shear stresses relax on the very time scale that one wishes to monitor.
The purpose of this paper is to establish the theoretical framework for measuring a complete set of thermoviscoelastic response functions utilizing spherical symmetry. The thermomechanical equations describing the coupling of mechanics to thermodynamics have been well known for many years ͓1,2͔. We recently ͑with Olsen͒ presented the full analytic solution of the one-dimensional inertia-free case, i.e., where all motion is restricted to one direction and the wavelength of sound waves at the relevant frequencies is much larger than the sample size ͓23͔. The reader is referred to that paper as an introduction to the reasoning and techniques used in the present paper. Below, the spherically symmetric inertiafree case is treated, where the effective one-dimensional nature again allows the problem to be solved analytically. The solution is cumbersome, but once it has been arrived at a number of experimentally relevant special cases are fairly easy to work out.
The solution is formulated in terms of the so-called transfer matrix ͓23,30͔ that links infinitesimal variations at the boundaries for the following four quantities: entropy ͑or equivalently heat͒ input, volume displacement, pressure change, and temperature change. The transfer matrix-not to be confused with the transfer matrix of statistical mechanical models-is useful because it directly describes how the system interacts with its surroundings. These interactions are conveniently pictured and described via the energy bond technique ͓31-34͔. An energy bond is characterized by a displacement variable and an effort variable, generalizing the concepts of charge and voltage. The product of an effort and a differential displacement variable is a generalized work that gives the ͑free͒ energy transferred into the system from its surroundings.
In standard thermodynamics there are two energy bonds, one thermal and ͑if the shear modulus is negligible͒ one mechanical; see Fig. 1͑a͒. The thermal energy bond is characterized by entropy S as the displacement variable and temperature T as the effort variable. For the mechanical energy bond the displacement is the volume V, and the effort is the negative pressure −p. The time derivatives of displacement variables define flow variables, signaled by arrows in Fig. 1. The two energy bonds of Fig. 1͑a͒ represent the fundamental thermodynamic identity In linear irreversible thermodynamics, small perturbations in temperature and pressure are linear functionals of the volume and entropy flows, or vice versa. The temperature perturbation dT = T − T 0 is around a reference temperature T 0 , and the pressure perturbation dp = p − p 0 is around a reference pres-sure p 0 . The supplied heat dQ is related to the externally supplied entropy, dS ext = dQ / T. When the volume V 0 is so small that temperature and pressure can be assumed homogeneous throughout the volume, the relaxation is described in the time domain by memory kernels of the form where T is the isothermal compressibility, ␣ p the isobaric expansion coefficient, and c p the isobaric specific heat per unit volume. Instead of Eq. ͑1͒ we now have The function E − T 0 S ext + p 0 V is the imparted free energy, which is not a state function since dissipation in the system degrades this energy: During a cyclic process one has where ⌺ is the entropy production in the system. Note that the entropy production is quadratic in the perturbations, implying that entropy is conserved to first order ͓3͔, a fact that is utilized below. It follows from the above that the generic conjugated variables for a small system that relaxes mechanically and thermally with no shear forces are ͑dT , dS ext ͒ and ͑−dp , dV͒. From an experimental point of view, however, it is not convenient only to consider an infinitesimal volume V 0 . If one wishes to study relaxation on a time scale expt , the volume V 0 may only be considered small if the heat diffusion time D across the volume is much smaller than expt . This is a restriction that can be coped with only by analyzing the influence of heat diffusion on the response-an important purpose of this paper. A further complication arises when the shear modulus becomes comparable to the bulk modulus, which is the case in the relaxation region of viscous liquids. In this case heat diffusion and mechanical stresses couple nontrivially. To keep the discussion below as simple as possible, we look at the situation with highest symmetry, that of a sphere of inner radius r 1 and outer radius r 2 , posing the question: What is the relation between the thermal and mechanical variables at the boundaries of the system? In order to take the shear forces properly into account, it is shown below that pressure must be replaced by the "radial pressure" ␦p r,1 at r 1 and ␦p r,2 at r 2 . ␦V 1 is the volume swept by the surface at r 1 as a consequence of a radial small displacement u͑r 1 ͒. Correspondingly, ␦V 2 is the volume swept at r 2 ͑both FIG. 1. Energy bond graphs ͓31-34͔, a useful tool for modeling linear thermoviscoelasticity. ͑a͒ The standard thermodynamic energy bonds, one thermal and one mechanical. In energy bond terminology, the temperature and negative pressure, T and −p, are so-called effort variables analogous to voltage, and the entropy flux Ṡ = dS / dt and volume flux V = dV / dt are so-called flow variables analogous to electrical current. In quasiequilibrium the product of an effort variable and its flow variable is the energy flux into the system from its surroundings; more generally it gives the flux of free energy-in the sense of available work-into the system. ͑b͒ Symbolic figure of the energy bonds of linear irreversible thermodynamics for the spherically symmetric situation. For the thermal energy bonds the efforts are the temperature variations ␦T 1 and ␦T 2 at radius r 1 and r 2 , respectively, and the flows are the entropy fluxes ␦Ṡ 1 and −␦Ṡ 2 into the system ͑fluxes are in the positive radial di-rection͒. For the mechanical energy bonds the flows are the volume fluxes −␦V 1 and ␦V 2 into the system, respectively, whereas pressure is now the "radial pressure," ␦p r,1 and ␦p r,2 , respectively, defined in Eq. ͑44͒. The four energy bonds define eight variables. As shown in this paper, the fundamental physical equations provide four constraints among these eight variables. Thus there is a linear relation between the four "outer" and the four "inner" energy bond variables. This relationship is expressed in terms of the 4 ϫ 4 transfer matrix calculated below that provides all information needed to interpret any experimental situation characterized by specific boundary conditions. in the positive radial direction͒. The net volume change is dV =−␦V 1 + ␦V 2 . The temperature perturbations at the two surfaces are denoted by ␦T 1 and ␦T 2 , respectively, and the entropy fluxes in the positive radial direction by ␦S 1 and ␦S 2 , respectively. The net entropy influx is dS ext = ␦S 1 − ␦S 2 . Now Eq. ͑5͒ becomes d͑E − T 0 S ext + p 0 V͒ = ␦T 1 ␦S 1 − ␦T 2 ␦S 2 + ␦p r,1 ␦V 1 − ␦p r,2 ␦V 2 .

͑7͒
As illustrated in Fig. 1͑b͒, this gives rise to four energy bonds, two referring to the outer radius r 2 and two to the inner radius r 1 .
The transfer matrix T͑r 2 , r 1 ͒ is by definition the 4 ϫ 4 matrix that links the ͑generally complex and frequencydependent, see below͒ infinitesimal variations of the four energy bond variables at radius r 1 with the four energy bond variables at radius r 2 . Switching notation such that ␦p r,1 is denoted by ␦p r ͑r 1 ͒, etc., the transfer matrix is thus defined by ␦p r ͑r 2 ͒ ␦T͑r 2 ͒ ␦V͑r 2 ͒ ␦S͑r 2 ͒ = T͑r 2 ,r 1 ͒ ␦p r ͑r 1 ͒ ␦T͑r 1 ͒ ␦V͑r 1 ͒ ␦S͑r 1 ͒ . ͑8͒ The transfer matrix is related to the response matrix ⌫͑r 2 , r 1 ͒ that by definition links the four effort variables to the four displacement variables: ␦p r ͑r 1 ͒ ␦T͑r 1 ͒ ␦p r ͑r 2 ͒ ␦T͑r 2 ͒ = ⌫͑r 2 ,r 1 ͒ From the fluctuation-dissipation theorem the response matrix is known to be symmetric, a fact that is explicitly confirmed below.
In the next section fundamentals are summarized. In Sec. III the full dynamic equations are formulated and brought into dimensionless form by scaling with complex units. In Sec. IV the equations are solved, and the transfer and response matrices are calculated. In Sec. V several experimentally relevant special cases are considered. Finally, Sec. VI gives a brief discussion.

II. DEFINITIONS AND CONSTITUTIVE RELATIONS
The energy bond variables of standard thermodynamics give rise to a number of dc ͑i.e., static͒ linear-response coefficients as follows. If the variables of interest are those of the two thermodynamic energy bonds of Fig. 1͑a͒, ͑T , p , S , V͒, there are altogether 24 thermodynamic coefficients of the form ͑‫ץ‬a / ‫ץ‬b͒ c with a, b, and c chosen among T, p, S, and V ͓25,35,36͔. These coefficients form 12 pairs that are trivially related by inversion ͓͑‫ץ‬a / ‫ץ‬b͒ c =1/ ͑‫ץ‬b / ‫ץ‬a͒ c , etc.͔. As is well known, the 12 coefficients are not independent, but related by Maxwell relations. This leaves the following eight basic linear-response coefficients ͑where the specific heats here and throughout the paper are per unit volume͒: Isobaric specific heat c p ϵ T Isobaric expansion coefficient Adiabatic contraction coefficient

͑17͒
Some well-known relations between the thermodynamic coefficients are summarized in the Appendix for reference, where relations are simplified somewhat by letting the heat capacities be represented by the related variables V ϵ 1 V ͑ ‫ץ‬S ‫ץ‬T ͒ V and p ϵ 1 V ͑ ‫ץ‬S ‫ץ‬T ͒ p . In systems with relaxing degrees of freedom the thermodynamic coefficients generally become complex and frequency dependent. Suppose, for instance, that the system is subjected to an infinitesimal periodic pressure variation with angular frequency described as p͑t͒ = p 0 +Re͓␦p exp͑st͔͒, where s = Ϯ i, depending on convention ͓23͔, is the socalled Laplace frequency. The volume then varies periodically as V͑t͒ = V 0 +Re͓␦V exp͑st͔͒, where both ␦p and ␦V are generally complex and frequency dependent. If this takes place at constant temperature, the complex frequencydependent isothermal compressibility is defined via T ϵ −␦V / ͑␦pV 0 ͒, where V 0 is the average volume. According to a basic theorem of linear irreversible thermodynamics, the Maxwell relations among the dc linear-response quantities in Eqs. ͑14͒-͑17͒ translate into Onsager relations for the frequency-dependent coefficients, which reflect time reversibility. This is a special case of the so-called correspondence principle ͓1,2,6͔: Any ͑dc͒ thermodynamic relation or equation involving linear thermodynamic and/or mechanical quantities applies unchanged when constitutive properties are replaced by the corresponding complex, frequency-dependent quantities. These frequency domain functions are related to the corresponding memory kernels like, e.g., where we adhere to the ordinary ͑sloppy͒ notation of physics that it is to be read from the context whether T refers to a time or frequency domain function. In the following we stay in the frequency domain, though. Thermoviscoelasticity describes the coupling between thermal and mechanical deviations from equilibrium. This paper deals with the linear case that is well understood as regards fundamentals. Linearity means that the system is assumed to be infinitesimally close to equilibrium. Deviations from equilibrium are quantified in terms of the infinitesimal displacement field u = u͑r , t͒, temperature variation field ␦T͑r , t͒ = T͑r , t͒ − T 0 , etc. In this approximation the radial heat displacement is given by ␦Q͑r͒ = ␦S͑r͒ / T 0 ; we switch to Q simply because heat capacity traditionally is defined via this quantity.
The isothermal bulk modulus K T ͑inverse isothermal com-pressibility͒ is defined by If ‫ץ‬ i is the derivative with respect to the ith spatial coordinate The relative volume change is given by the trace of the strain tensor: Denoting the stress tensor by គ គ = ij , the shear modulus G is defined ͓3͔ as follows ͑the isothermal and adiabatic shear moduli are always identical͒: If there are both infinitesimal displacements and spatial temperature variations, this is generalized into the so-called Duhamel-Neumann relation ͓3,4͔ that is the following constitutive relation linking mechanical and thermodynamic properties: For relaxing systems, K T , G, and ␤ V are complex and frequency dependent, and for periodically varying boundary conditions both strain and stress tensors are generally complex. By the correspondence principle the Duhamel-Neumann relation applies also for the frequency-dependent case.
The isothermal longitudinal modulus M T is defined ͓3͔ by Similarly, the adiabatic longitudinal modulus M S is defined by ͓3͔͒, a useful quantity for the description of one-dimensional thermoviscoelasticity is the "longitudinal specific heat" defined by An important result of the present paper is that the longitudinal specific heat also plays a significant role for the spherically symmetric case. Again, c l is frequency dependent for relaxing systems. A useful identity ͓23͔ follows from Eqs. ͑A1͒ and ͑A7͒: From this one finds that c l may be interpreted as a generally complex convex combination of c p and c V : This shows that the longitudinal heat capacity is effectively in between the isobaric and isochoric heat capacities. In analogy to the standard abbreviation we define ␥ l , a quantity that is also generally complex and frequency dependent, as follows ͓23͔: If the heat-current density is denoted by j, the heat conductivity is defined via Fourier's law, The heat conductivity is generally assumed to be frequency independent, an assumption that was recently confirmed ͓37,38͔. The solution below, however, applies also if were to depend on frequency.

A. The frequency-independent case
As mentioned, entropy conservation always applies to first order, a fact that in Ref. ͓3͔ is referred to as "the equation of continuity for heat." In terms of the infinitesimal displacement field u͑r , t͒ and the infinitesimal temperature field ␦T͑r , t͒, for a volume element dV at average mass den-sity 0 , the basic thermoviscoelastic equations of motion reflecting entropy conservation and Newton's second law

͑32͒
We shall be concerned only with the inertia-free limit, i.e., the limit where the sample is much smaller than the wavelength of sound waves at the frequencies of interest. In this limit the acceleration term is negligible; thus

͑34͒
Before proceeding we note that the identity ␤ V = ␣ p K T ͓Eq. ͑A5͔͒ implies that, if there is no thermal expansion upon heating, then ␤ V = 0. In this case the two equations decouple and reduce to the ordinary elastic equation of motion in the inertia-free limit and the heat-conduction equation, respectively. Thus the coupling between mechanics and thermodynamics arises only when the thermal expansion coefficient is nonzero, as is indeed intuitively obvious.
The pressure variation ␦p is defined by ͑assuming here and henceforth that the average stress tensor is zero͒ This definition applies also for nonhydrostatic conditions, i.e., when គ គ is not proportional to the unit tensor. When the specific heat or the expansion coefficient is termed isobaric, it refers to a situation where the trace of the stress tensor is constant, not necessarily its individual diagonal components. The relative volume change of Eq. ͑21͒ is related to changes in pressure, ␦p, and temperature, ␦T, by the following equation ͑see, e.g., Ref. ͓23͔͒:

͑36͒
Using Eqs. ͑36͒ and ͑A1͒ one finds that under isobaric con- where the heat-diffusion constant D p involves the isobaric specific heat D p = / c p . In general, isobaric conditions do not apply and the full coupled system Eqs. ͑33͒ and ͑34͒ must be solved.

͑39͒
Here a 1 is an integration constant that is a function of time determined by the boundary conditions. Substituting Eq. ͑39͒ into Eq. ͑34͒ yields

B. The case of spherical symmetry and periodically varying fields
We now assume that all fields depend only on radius r and that their time dependence is harmonic, i.e., proportional to exp͑st͒ ͑s = Ϯ i͒. Henceforth u and ␦T refer to the complex frequency-dependent amplitudes of the infinitesimal radial displacement and temperature fields, respectively.
With these assumptions, if differentiation with respect to r is denoted by a prime, using the identity ١ · u = r −2 ͑r 2 u͒Ј Eqs. ͑39͒ and ͑40͒ become and c l s␦T + T 0 ␤ V sa 1 = r −2 ͑r 2 ␦TЈ͒Ј.

͑42͒
Throughout the paper it is important to remember that not only are the field amplitudes u and ␦T generally complex functions of s, but so are all constitutive properties like ␤ V , M T , etc. For simplicity of notation, however, the frequency dependence will usually not be explicitly indicated. We need two further equations for two auxiliary fields. One is the "radial pressure" i.e., the normal force per unit area in the negative radial direction exerted on a spherical surface of radius r from the outside. The radial pressure is generally different from the pressure ␦p of Eq. ͑35͒; the relation between the two follows from the definition of the shear modulus Eq. ͑22͒: If G = 0-in particular for any liquid at zero frequency-the radial pressure equals the pressure. Since ⑀ rr = uЈ, where ⑀ rr is the rrth component of the strain tensor in spherical coordinates ͓3͔, Eq. ͑44͒ and the Duhamel-Neumann relation Eq. ͑23͒ imply The other auxiliary field to be introduced is the timeintegrated heat-current density ␦q. Invoking Fourier's law of heat conduction Eq. ͑30͒, we get

C. Scaling to dimensionless variables
The solution of the thermoviscoelastic equations in spherical symmetry is extremely involved when expressed in terms of the original variables. We solve the equations below by proceeding in analogy to the solution of the onedimensional problem detailed in Ref. ͓23͔, where considerable simplification was obtained by scaling with suitable variables. One is free to scale by a characteristic length, time, temperature, and mass, or other combinations of these four fundamental dimensions. We choose pressure instead of mass as the fourth fundamental dimension. There is no mathematical requirement that scaling variables must be real; in the frequency domain they may well be complex and frequency dependent.
The scaling units are chosen as follows. The unit of time is the inverse complex Laplace frequency, s −1 . The unit of length is the ͑complex and frequency-dependent͒ heatdiffusion length It is occasionally convenient to use the ͑complex and frequency-dependent͒ wave number k ϵ 1 / l D , The unit of temperature is the average temperature T 0 , and the unit of pressure is the ͑complex and frequency-dependent͒ isothermal bulk modulus K T . With these units we define the following dimensionless variables: It should once again be emphasized that all thermodynamic and mechanical constitutive quantities are generally complex and frequency dependent, although for simplicity of notation this fact is rarely explicitly indicated below. The scaled problem involves only the three independent frequency-dependent response functions c, g, and ␣ . Nevertheless, it is convenient to introduce two further dimensionless frequency-dependent parameters ␥ and ␥ l of Eqs. ͑28͒ and ͑29͒, which are functions of c, g, and ␣ . Before proceeding we note that via the above definitions and Eqs. ͑A1͒, ͑A4͒, and ͑A5͒, the parameters ␥ and ␥ l obey These equations imply the following identity that turns out to be useful: When rewritten in terms of the above dimensionless variables, Eqs. ͑42͒, ͑41͒, ͑45͒, and ͑46͒ become ͑where the prime now implies differentiation with respect to r͒

͑64͒
These are the fundamental dimensionless thermoviscoelastic equations of spherical symmetry to be solved.

IV. SOLUTION IN TERMS OF THE TRANSFER MATRIX
Equation ͑62͒ leads to by which Eq. ͑63͒ is simplified as follows: ␦p r = − ͑1 + g͒a 1 + 3gr −1 ũ .

͑101͒
These relations follow from Eq. ͑82͒ and the connection between ⌫ and T . On the other hand, the symmetry relation

͑102͒
and that

A. A massive sphere
In this first application of the formalism ͑Fig. 2͒ we inquire into how a solid sphere responds to a compression −␦V and a heat supply −␦Q ͓39͔ applied at radius r 2 ͑the transferred heat is positive when ␦Q Ͻ 0 because ␦Q refers to the heat flow in the positive radial direction͒. This is calculated from the transfer matrix by putting r 1 =0. If r 2 Ӷ ͉l D ͉, we expect the response matrix G 0 to be given by the constitutive equations ͑13͒, ͑15͒, and ͑10͒ generalized to frequencydependent coefficients ͑where all variables refer to the outer radius r 2 ͒: FIG. 2. Compression −␦V and heat input −␦Q at the surface of a sphere give rise to changes ␦p r and ␦T in radial pressure and temperature, respectively, at the surface. For a finite sphere, the response ͓Eq. ͑107͔͒ is given by the constitutive properties K S , ␣ S , and c V , as well as by heat diffusion. Eight kinds of input-output relations are discussed in the text.
In general, the relation is This follows from Eqs. ͑76͒ and ͑87͒ via the boundary conditions ␦Ṽ ͑r 1 =0͒ = 0 and ␦Q ͑r 1 =0͒ = 0. Returning to dimensional variables, Eqs. ͑105͒ and ͑106͒ yield where f D is a function of frequency defined by

͑108͒
The frequency dependence of f D derives primarily from that of l D . Note that f D → 1 for → 0 and f D → ϱ for → ϱ.
Asymptotically, one has f D = r 2 / ͑3␥ l l D ͒ for → ϱ. The function f D describes how heat diffusion affects the measurement of the thermal and mechanical properties of a massive sphere that can only, of course, be accessed at the surface. Below, several different thermoviscoelastic experiments on a massive sphere are considered. In principle there are 24 such experiments ͓25͔, corresponding to the 24 coefficients discussed at the beginning of Sec. II: One may choose any of the four variables ␦p r , ␦T, −␦V, or −␦Q at the outer radius r 2 as the controlled input, any of the three remaining variables as the measured output, and either of the last two to be fixed. Since there are only three independent functions of frequency in the matrix of Eq. ͑107͒, it is superfluous to discuss all 24 experiments. It does make sense, however, to discuss more than just three cases, because the experimental challenges may vary from case to case.
We discuss eight cases, corresponding in the lowfrequency limit to the eight frequency-dependent thermodynamic response functions T , S , c V , c p , ␣ p ␣ S , ␤ V , and ␤ S . These cases are detailed below where, occasionally, relations from the Appendix are utilized.

(a) Compression with isothermal boundary
At low frequencies this converges to K T / V. This result shows how one is limited upward in frequency when attempting to do isothermal bulk modulus measurements. Because f D → ϱ at large frequencies, the equation also describes the transition to the adiabatic bulk modulus above the characteristic heatdiffusion frequency D defined by D ϵ D / r 2 2 , where D is the heat diffusion constant and r 2 the sample size.
The bulk modulus can be measured in the frequency range 1 Hz-50 kHz by the so-called piezoelectric bulk modulus gauge ͑PBG͒ ͓16͔. This is a piezoelectric ceramic hollow sphere that may be filled with liquid. The ceramic shell has electrodes on the inside and outside and thus constitutes an electrical capacitor. Due to the piezoelectric effect the frequency-dependent capacitance-which can readily be measured-depends on the bulk modulus of the liquid which can be found after a calibration of the PBG. This device has a radius of 10 mm, and since typical liquid heat diffusivities are of order 0.1 mm 2 / s, the characteristic heat diffusion frequency becomes 10 −3 s −1 . Thus experiments performed with the PBG above 0.1 Hz can safely be said to be adiabatic, despite the fact that no special measures are taken to make the boundary conditions adiabatic.
(b) Adiabatic compression Ironically, this boundary condition is difficult to achieve experimentally, whereas the isothermal experiment gives K S at most frequencies. Thus K S is easier to measure than K T ͓16͔.
(c) Isochoric heating The low-frequency limit is 1 / ͑Vc V ͒ giving the isochoric frequency-dependent specific heat. The high-frequency limit, however, involves the longitudinal specific heat. Note that in this limit-even though the overall volume is constant-it is c l that appears, not c V .

͑113͒
The low-frequency limit is 1 / ͑Vc p ͒, giving the isobaric frequency-dependent specific heat. The high-frequency limit, is identical to the isochoric high-frequency limit. Note that it is the longitudinal specific heat that enters into Eq. ͑114͒, not c p . This is similar to the fact that the frequency-dependent specific heat obtained from plane-wave effusivity measurements is not the isobaric specific heat, but the longitudinal ͓23͔. In that case nonisotropic stresses could be conceived as arising from the special kind of mechanical boundary conditions needed in order to keep the model of the plane-plate setup one dimensional. Here we see, however, that nonisotropic stresses may arise in the liquid itself, not necessarily coming from clamping the boundaries. This substantiates a conclusion of Ref. ͓23͔, namely, that it is not possible to probe the isobaric specific heat directly by effusivity measurements. Note also that the thermal admittance per unit area is the same as for the planar geometry ͓23͔:

͑115͒
Recently, radial heat effusion from the surface of a spherical cavity inside an infinite medium was shown also to involve the longitudinal specific heat ͓40͔.

͑116͒
At low frequencies this approaches ␣ p , whereas its highfrequency asymptotic form is given by This is a radial version of the Onsager relation that corresponds to the Maxwell relation Eq. ͑14͒. If an experiment were conceived where one measures the heat flux −␦Q needed to keep the temperature constant at the surface while applying a periodically varying radial pressure, one would find the same response function Eq. ͑116͒, including the diffusion dependence. Similar radial versions of Onsager relations corresponding to the Maxwell relations Eqs. ͑15͒-͑17͒ hold for the last three response functions.
(f) Radial pressure in response to heating at constant volume Thus ␣ S can be measured without interference from heat diffusion. This is not trivial, since the penetration depth ͉l D ͉ of the temperature field into the sphere is frequency dependent, and it is the temperature field that creates the pressure variation.
(g) Radial pressure in response to a controlled temperature oscillation at constant volume. This case leads to In contrast to case ͑f͒, this response function is diffusion influenced. It approaches ␤ V at low frequencies, whereas (h) Volume expansion in response to heating for a free surface As for case ͑f͒ we get a simple result that is independent of heat diffusion. In summary, the important role played by the longitudinal frequency-dependent specific heat is evident. Moreover, there is now an exact description of the transition from the adiabatic to the isothermal regime of bulk modulus measurements utilizing the PBG ͓16͔.
B. The "thermally thick limit" ͦl D ͦ ™ r 2 when r 1 ™ r 2 We now proceed to discuss the case where the inner radius is small and the sample is much larger than l D . This section prepares the theoretical basis of ongoing experiments where a small spherical thermistor placed in the center of the PBG should make it possible to simultaneously measure the frequency dependences of ␣ S , K S , and c l on the same sample. Supplemented by shear modulus measurements ͓24͔, this provides a complete set of thermoviscoelastic response functions of a liquid. Below we determine the reduced transfer matrix X for a situation where a mechanical boundary condition at r l and a thermal boundary condition at r 2 are stipulated. That is, X gives the linear relationship There are four possibilities for X , denoted below by Ã , B , C , and D , depending on the different boundary conditions: ͑125͒ X = C for ␦p r ͑r 1 ͒ = 0, ␦Q ͑r 2 ͒ = 0, ͑126͒ X = D for ␦p r ͑r 1 ͒ = 0, ␦T͑r 2 ͒ = 0.

͑141͒
͑3͒ The relation between the heat displacement at radius r 1 and the volume and negative normal-stress variations at radius r 2 does not involve any "delay" caused by thermal diffusion:

C. Mechanical boundary conditions
If we control the boundary conditions solely via isobaric or isochoric constraints, there is a thermal transfer matrix Ỹ connecting heat and temperature variations at the inner and outer radii: Depending on boundary conditions, we define ͓42͔ Ỹ ͑r 2 ,r 1 ͒ = H ͑r 2 ,r 1 ͒ for ␦Ṽ ͑r 1 ͒ = 0 and ␦Ṽ ͑r 2 ͒ = 0,
FIG. 4. Energy bond graph diagram ͓33͔ of the physical interactions between the thermistor, liquid, and thermal leak to the cryostat, modeling the setup of Fig. 3. The thermistor acts as both heat source and thermometer. The liquid is described by the transfer matrix K for the case with mechanical clamping at the thermistor and a free outer surface.

͑157͒
We shall not explicitly give the components for the general case, but limit ourselves to the thermally thin limit r 2 Ӷ ͉l D ͉ ͑i.e., ͉kr 2 ͉ Ӷ 1͒ where the results simplify considerably. In this limit one finds ͓20͔ after Taylor expanding the matrix elements of T that play the role of thermal resistance and capacitance, respectively. These results imply that in the thermally thin limit c p is measured if r 2 ӷ r 1 , whereas c l is measured if r 2 Х r 1 . We finally note that the above formulation is easily incorporated into a model taking into account the thermal heat loss to the surroundings ͓20͔.

D. No thermomechanical coupling
As mentioned in the beginning of Sec. III, if the isobaric thermal expansion coefficient ␣ p is zero, one has ␤ V = 0 and there is no thermomechanical coupling. In this case, the components T 21 , T 23 , T 41 , and T 43 vanish, implying that heat and temperature variations at r 2 depend only on heat and temperature variations at r 1 . Also, T 12 , T 14 , T 32 , and T 34 vanish, implying that pressure and volume variations at r 2 depend only on pressure and temperature variations at r 1 . Note that when there is no thermomechanical coupling, all specific heats are identical: This is often a good approximation for solids, but rarely for liquids. When there is no thermomechanical coupling, heat diffusion is described by a 2 ϫ 2 thermal transfer matrix T th defined ͓43͔ as follows: ͩ ␦T͑r 2 ͒ ␦Q͑r 2 ͒ ͪ=T th ͑r 2 ,r 1 ͒ͩ ␦T͑r 1 ͒ ␦Q͑r 1 ͒ ͪ. ͑162͒ The components of T th are found by substituting ␣ p = 0 into Eq. ͑77͒. The results are as follows ͑where k =1/ l D ͒ T 11 th = r 1 r 2 cosh͓k͑r 2 − r 1 ͔͒ + 1 kr 2 sinh͓k͑r 2 − r 1 ͔͒, ͑163͒ T 12 th = − s 4 sinh͓k͑r 2 − r 1 ͔͒ kr 1 r 2 , ͑164͒ T 21 th = 4c l k 3 ͕͑1 − k 2 r 1 r 2 ͒sinh͓k͑r 2 − r 1 ͔͒ − k͑r 2 − r 1 ͒cosh͓k͑r 2 − r 1 ͔͖͒, ͑165͒ T 22 th = r 2 r 1 cosh͓k͑r 2 − r 1 ͔͒ − 1 kr 1 sinh͓k͑r 2 − r 1 ͔͒. ͑166͒ Interestingly, the same purely thermal 2 ϫ 2 transfer matrix describes a low-viscosity liquid ͑g → 0͒ even when ␣ p 0, if either the inner or the outer surface is free, i.e., if ␦p r ͑r 1 ͒ =0 or ␦p r ͑r 2 ͒ = 0. In these two cases, however, not all three specific heats are identical, only c p = c l applies. This is consistent with the above remark regarding the validity of the standard heat diffusion equation ͑37͒.

VI. CONCLUDING REMARKS
Thermoviscoelastic response functions are notoriously difficult to measure. This paper establishes the theoretical framework necessary for developing experimental methods that utilize spherical symmetry for measuring such response functions. From the complete solution of the problem in the form of the transfer matrix the equations describing any realistic experimental situation may be derived, as exemplified in the last section.
The thermoviscoelastic response functions are important to determine for liquids approaching the glass transition ͑still in metastable equilibrium above the transition͒. For such ultraviscous liquids all thermodynamic coefficients become complex and frequency dependent for frequencies in the range of the inverse Maxwell relaxation time. To the best of our knowledge there are yet no reliable measurements of a complete set ͑i.e., three ͓1,2,4,25,28͔͒ of thermoviscoelastic response functions for any such liquid. The determination of such complete sets, from which all other thermoviscoelastic response functions are easily calculated, serves the obvious purpose of elucidating the macroscopic dynamics and thermodynamics of ultraviscous liquids. Very recent theoretical developments even further stress the importance of developing reliable methods for measuring thermoviscoelastic response functions. It now appears that the class of van der Waals liquids ͑possibly supplemented by some liquids forming bulk metallic glasses͒ have particularly simple properties: Liquids like these with nondirectional chemical bonds exhibit strong correlations between equilibrium pressure and energy fluctuations ͓28,44͔. For such "strongly correlating viscous liquids" it has been shown that there is basically only one independent thermoviscoelastic response function ͓25͔. It would be interesting to have this prediction subjected to experimental tests. Moreover, for strongly correlating viscous liquids there are indications from computer simulations that thermoviscoelastic measurements can determine the exponent of the so-called density scaling that collapses the relaxation time's pressure and temperature dependence onto a master curve ͓45͔.
As regards the above results, it is notable that the longitudinal specific heat c l plays a dominant role. It is perhaps not surprising that c l enters repeatedly into the equations describing the one-dimensional case ͓23͔-after all, this is what it was defined to do-but it is less obvious that c l also plays a dominant role for the case of spherical symmetry.