Previous use of the 3 omega method has been limited to materials with thermal conductivity tensors that are either isotropic or have their principal axes aligned with the natural cartesian coordinate system defined by the heater line and sample surface. Here, we consider the more general case of an anisotropic thermal conductivity tensor with finite off-diagonal terms in this coordinate system. An exact closed form solution for surface temperature has been found for the case of an ideal 3 omega heater line of finite width and infinite length, and verified numerically. We find that the common slope method of data processing yields the determinant of the thermal conductivity tensor, which is invariant upon rotation about the heater line’s axis. Following this analytic result, an experimental scheme is proposed to isolate the thermal conductivity tensor elements. Using two heater lines and a known volumetric heat capacity, the arbitrary 2-dimensional anisotropic thermal conductivity tensor can be measured with a low frequency sweep. Four heater lines would be required to extend this method to measure all 6 unknown tensor elements in 3 dimensions. Experiments with anisotropic layered mica are carried out to demonstrate the analytical results.

Anisotropic materials are relevant for a wide range of applications such as thermoelectrics,1–3 high-temperature superconductors,4 and lightweight heat spreaders.5 Often the crystallographic orientation of these anisotropic materials can be controlled depending on the material synthesis process,6–9 while in other cases, the orientation is unknown or may vary. However, the vast majority of thermal conductivity measurement techniques are restricted to isotropic materials, or anisotropic materials with principal conductivities aligned parallel with a sample’s orthogonal surfaces. Very recently, Feser et al.10 developed a time-domain thermoreflectance (TDTR) method which can measure an arbitrarily aligned thermal conductivity tensor using offset pump and probe beams, although the absence of a closed form analytic solution necessitates a time-consuming iterative fitting process. No analogous scheme has been developed for the electrothermal-based 3ω method,11,12 which uses simpler hardware than thermoreflectance-based techniques. Also, the two methods have fundamentally different heating geometries (spot vs. line) making it non-trivial to adapt analytic solutions from the TDTR method.

The 3ω method11,12 has been extensively used to measure the thermal conductivity in bulk substrates as well as thin films. Briefly, an alternating current source at an angular frequency ω is applied along a heater line deposited on top of the substrate of interest. The resulting joule heating causes a temperature field fluctuating at the 2nd harmonic. The temperature effect on heater line resistance combined with the 1ω current results in a voltage at 3ω. For semi-infinite substrates in the low frequency limit, the slope of this 3ω voltage with respect to the logarithm of frequency is inversely proportional to the thermal conductivity of the substrate (the “slope method”). The 3ω method has been extended to find the thermal conductivity of thin films using the differential 3ω method.13,14

For multi-layered anisotropic substrates with the materials’ principal axes aligned with the heater line and free substrate surface, Borca-Tasciuc et al.15 developed a general solution using integral Fourier transforms.14,16 They proposed an experimental scheme with multiple heater linewidths to extract the thermal conductivity tensor elements for aligned anisotropic layers. They showed that the slope of temperature versus logarithm of frequency is related to the geometric mean of the 2 principal thermal conductivities characteristic to the plane being probed by the 3ω line. Tong and Majumdar17 developed a data reduction scheme to determine the aligned anisotropy ratio of thermal conductivity and the interfacial thermal conductance in film on substrate systems. Ramu and Bowers18 gave a closed form solution for surface temperature in 3ω experiments on an aligned anisotropic substrate with small anisotropy ratios. Finite element method (FEM) simulations and a parallel offset thermometer line were used together to find the cross and in-plane thermal conductivity using an iterative approach in a later work.19 

Thus, all of these previous anisotropic 3ω techniques have been restricted to the aligned case: samples with principal thermal conductivities aligned parallel and perpendicular to the sample surface and heater line. This restriction eliminates the off-diagonal terms from the thermal conductivity tensor and thus simplifies the analysis, but limits the scope of applicability. An experimental procedure which can also identify any off-diagonal terms can also give insight into the orientation and symmetries of the crystal structure as a consequence of the Neumann’s Principle, which states that the physical properties of a material are at least as symmetric as its underlying atomic structure.20 

In this work, we solve the 3ω heating problem for the general case of a substrate with an anisotropic thermal conductivity tensor, which may be arbitrarily oriented with respect to the natural cartesian coordinate system defined by the heater line and free substrate face. We obtain the solution in closed form rather than in an integral form, thereby avoiding the need for time-intensive and potentially non-intuitive numerical fitting to extract the thermal conductivity tensor.

We first derive a closed form solution for the surface temperature profile in the case of an infinitesimally narrow line heater on an arbitrary anisotropic substrate and show it to be of similar form as that for an isotropic substrate. The solution is extended for a finite-width heater line and verified using FEM software. Experimental schemes are devised to use measurements of the heater line temperature in the low frequency regime to recover the complete non-diagonal thermal conductivity tensor. We demonstrate these predictions with measurements on mica, a layered anisotropic material.

In a cartesian coordinate system, the heat flux vector Q is related to the thermal gradient vector by Fourier’s Law,

Q i = j k ij T x j ,
(1)

and the most general anisotropic thermal conductivity tensor is

k = k xx k xy k x z k yx k yy k y z k zx k z y k z z .
(2)

k is real and symmetric, and it is known21 that kii > 0 and kiikjj > kij2. Thus, there are at most 6 independent elements for the most general anisotropic k. We express k in cartesian coordinates and require it to be independent of position, so that the analysis applies to layered materials like graphite, but not curved Swiss-roll topologies like a tree trunk.

The basic 3ω setup is depicted in Figure 1(a). We define the z axis to lie along the centerline of the heater line, and the x axis as orthogonal to both the heater line and the plane of the sample surface. The alternating current flowing through the heater line causes sinusoidal joule heating at the second harmonic, which is uniform over the heater line length and width. The magnitude of this heater power per unit length is defined as P′. With the usual assumption of an infinitely long heater line, there cannot be any temperature gradients in the z direction which simplifies the problem to the 2-d cross section shown in Figure 1(b). While it is possible to use Feldman’s algorithm16 to solve for this heating problem on a multi-layered substrate in a complicated integral form, following a process similar to that used by Feser et al.,10 the focus here instead is on achieving a closed form solution for an arbitrarily anisotropic substrate.

FIG. 1.

(a) The 3ω setup and (b) 3ω heater line on an anisotropic semi-infinite substrate. The cross-hatching in (b) is meant to suggest angled layers of an anisotropic material such as graphite or mica.

FIG. 1.

(a) The 3ω setup and (b) 3ω heater line on an anisotropic semi-infinite substrate. The cross-hatching in (b) is meant to suggest angled layers of an anisotropic material such as graphite or mica.

Close modal

We first find the solution for an infinitesimally narrow line heater deposited on such a substrate by using integral transform techniques. This infinitesimal line source is represented as a δ-function heat generation at x = 0, y = 0. The time-periodic heat generation is represented using the complex Euler notation, whereby the real and imaginary parts of T represent the in-phase and out-of-phase temperatures, respectively.22 

The governing equation is

k xx 2 T 2 x + 2 k xy 2 T x y + k yy 2 T 2 y + P e i 2 ω t δ ( x , y ) = C T ,
(3)

where C is the volumetric heat capacity. Since the sample is very large in the ±y-direction, suitable boundary conditions are

T y = 0 and T = 0  at  y ± .
(4)

In the x-direction, the top surface (x = 0) is insulated while the bottom surface is far away from the heater line, and thus,

k xx T x + k xy T y = 0 at x = 0 ,
(5)
T x = 0 and T = 0  at  x .
(6)

Exploiting the known 2ω forcing of this linear system, the time dependence of the governing equation can be removed by introducing T ¯ such that

T x , y , t = T ¯ x , y e i 2 ω t ,
(7)

where i = 1 . The factor of 2 arises because we follow the usual convention of defining ω as the frequency of the electrical current, so that 2ω represents the joule heating. A series of integral and Fourier transforms23 are used to remove x and y derivatives from the governing equation. Deferring the mathematical details to Appendix  A, Eqs. (A1)(A10), the solution for T ¯ is

T ¯ ( x , y ) = P π k xx k yy k xy 2 K 0 q xx k xx k yy x 2 + k x x 2 y 2 2 k xx k xy x y k xx k yy k xy 2 ,
(8)

where K0 is the modified Bessel function of the second kind and

q xx = i 2 ω C / k xx ,
(9)

where q xx represents the inverse of the penetration depth in the x-direction. As a check, when kxx = kyy and kxy = 0, the familiar isotropic solution from Cahill11 is correctly recovered.

Figure 2 shows an example of the temperature field resulting from a δ-function line heater, calculated with Eq. (8). The argument of Eq. (8) also shows that all isotherms are elliptical in shape and obey the relation

k xx k yy x 2 + k x x 2 y 2 2 k xx k xy x y = const ,
(10)

where const is a constant determining the temperature of the isotherm. It can be shown that the major and minor axes of these ellipses are aligned along the two principal thermal conductivity directions (proof in Appendix  B). This result is intuitive since the isotherms are stretched in the direction of least resistance or maximum thermal conductivity. It is also interesting to note that the δ-function line heater is a center of inversion symmetry in the xy coordinate system. This can also be confirmed mathematically by replacing (x, y) by (−x, − y) in the governing equation and the boundary conditions, which shows that the surface temperature along x = 0 is symmetric about the line heater even with the presence of cross terms in the thermal conductivity tensor.

FIG. 2.

Isothermal contours for the in-phase temperature response to a δ-function line heater on an anisotropic substrate. This example uses an anisotropic material with principal thermal conductivities 1 W/m-K and 5 W/m-K, rotated by 30° in the clockwise direction so that k = 2 3 3 4 W / m - K in the xy coordinate system. We define an effective thermal conductivity, k det = k xx k yy k xy 2 , with its characteristic thermal penetration depth, L c = k det / 2 ω C , which contains the only ω dependence. kdet is the square root of the determinant of k. A generic heat capacity of C = 2 × 106 J/m3-K is also used.

FIG. 2.

Isothermal contours for the in-phase temperature response to a δ-function line heater on an anisotropic substrate. This example uses an anisotropic material with principal thermal conductivities 1 W/m-K and 5 W/m-K, rotated by 30° in the clockwise direction so that k = 2 3 3 4 W / m - K in the xy coordinate system. We define an effective thermal conductivity, k det = k xx k yy k xy 2 , with its characteristic thermal penetration depth, L c = k det / 2 ω C , which contains the only ω dependence. kdet is the square root of the determinant of k. A generic heat capacity of C = 2 × 106 J/m3-K is also used.

Close modal

The next step in this analysis is to find the surface temperature profile, Tsurf(y), caused by a finite width heater, as well as the average temperature of this heater line which is what the 3ω voltage measures. The finite width heater is represented as a superposition of δ-function line heaters over the heater linewidth, leading to a convolution integral,

T surf ( y ) = 1 2 b b b T ̄ ( 0 , y y ) d y = 1 2 b b b P π k xx k yy k xy 2 K 0 q xx k xx y y k xx k yy k xy 2 d y ,
(11)

where the latter form follows from Eq. (8) after setting x = 0 to represent the surface. As mentioned after Eq. (8), for a δ-function line heater, the temperature solutions for anisotropic and isotropic substrates have very similar forms. It turns out that one can map the integral in Eq. (11) onto the corresponding integral for an isotropic substrate using a change of variables. As detailed in Appendix  A, Eqs. (A11) and (A12), in this way, we can obtain the surface temperature for an anisotropic substrate analogous to isotropic solutions found by Duquesne et al.24 The result here is

T surf ( y ) = P 4 k xx k yy k xy 2 1 Y Ξ i Ω 1 Y 2 + 1 + Y Ξ i Ω 1 + Y 2 ,
(12)

where Ω = b 2 2 ω C k xx k yy k xy 2 , Y = y b , and Ξ(y) = K0(y) L−1(y) + K1(y) L0(y), where Kn and Ln are the n-th order modified Bessel function of the second kind and Struve function,25 respectively. Similarly, the average of Tsurf over the heater linewidth, T , can also be found by using a change of variables to adapt the known isotropic solution24 for this anisotropic problem (Appendix  A, Eqs. (A13) and (A14)). The solution can be represented in terms of the Meijer-G function G 24 22 as24,26

T = i P 4 π Ω k xx k yy k xy 2 G 24 22 i Ω | 1 , 1 , 1 / 2 , 0 1 , 3 / 2 .
(13)

Equation (13) is plotted in Figure 3, along with low- and high-frequency asymptotic forms. The high-frequency limit of Eq. (13) is

T = P 2 k xx b q xx ,
(14)

which is the well known 1-D planar heating result where temperature gradients are present only in the x-direction.

FIG. 3.

Amplitude of in-phase and out-of-phase temperature oscillations of a finite-width heater line, as a function of heating current frequency (Eq. (13)). k det = det ( k ) = k xx k yy k xy 2 . In the low frequency limit, the slope ξ of the in-phase temperature rise is inversely proportional to kdet. In the high frequency limit, the magnitude of both in-phase and out-of-phase temperature rise is inversely proportional to k xx . The dashed lines depict the asymptotic forms given in the main text (Eqs. (14) and (15)).

FIG. 3.

Amplitude of in-phase and out-of-phase temperature oscillations of a finite-width heater line, as a function of heating current frequency (Eq. (13)). k det = det ( k ) = k xx k yy k xy 2 . In the low frequency limit, the slope ξ of the in-phase temperature rise is inversely proportional to kdet. In the high frequency limit, the magnitude of both in-phase and out-of-phase temperature rise is inversely proportional to k xx . The dashed lines depict the asymptotic forms given in the main text (Eqs. (14) and (15)).

Close modal

Of particular importance is the low frequency limit, because the simplicity and accuracy of the low-frequency slope method of 3ω data analysis are already well established for isotropic and aligned anisotropic substrates. Using the small argument approximation, K0(z) ≈ − ln(z) − γ + ln(2), we find the low-frequency limit of Eq. (13) as

T = P π k xx k yy k xy 2 1 2 ln b 2 2 ω C k xx k xx k yy k xy 2 + 3 2 γ i π 4 ,
(15)

where γ ≈ 0.5772 is the Euler constant.

Equation (15) is one of the major results of this paper and is used extensively in the experimental section below. Its slope, ξ = T / ln ω , is inversely proportional to

k d e t = k xx k yy k xy 2 = det k ,
(16)

with the proportionality constant −P′/2π. Importantly, we recognize kdet as the square root of the determinant of this 2-d thermal conductivity tensor, and as such it is invariant under rotation in the xy plane. Hence, changing the orientation of heater line by rotating it around its longitudinal (z) axis must always give the same slope value, and thus cannot yield any additional information to separate out the values of kxx, kyy, and kxy. On the other hand, Eq. (15) shows that the magnitude of temperature fluctuations does depend on the orientation of heater line. This fact will be exploited later to extract the individual thermal conductivity tensor elements.

Finally, the RMS 3ω voltage is related to the RMS 2ω temperature by11,12

V 3 ω = 1 2 α R I 1 ω T .
(17)

Here, I1ω is the RMS value of the AC source, R is the heater line resistance at zero current, and α is the temperature coefficient of resistance, all of which are easily determined experimentally making it straightforward to switch between V3ω and T .

To numerically validate our key analytical results, Eqs. (12), (13), and (15), simulations of a 3ω heating problem in 2-d were carried out using a commercial FEM package (COMSOL). To set up the heat conduction equations in the frequency domain, we used user-defined equations in the coefficient form PDE (partial differential equation) module. As a test case, we use the same k and C described in the caption of Fig. 2, with additional simulation parameters b = 5 μm and P′ = 3 W/m. For these parameters, the low frequency regime q xx b < < 1 and the high frequency regime q xx b > > 1 correspond approximately to ω < 103 rad/s and ω > 4 × 105 rad/s, respectively.

The validation is carried out in two steps. First, to check Eq. (12), in the most general case, we choose a mid-frequency regime (penetration depth comparable to heater linewidth: q xx b = 1 , corresponding to ω = 2 × 104 rad/s and Ω = 0.8), rather than either extreme frequency limits. To eliminate the effects of the far-field boundary conditions of the original semi-infinite problem in a finite-sized FEM domain, care is taken to ensure the sample extents in +x and ±y directions are much larger than the thermal wavelengths, confirmed by verifying that the same simulation results are obtained for both adiabatic and isothermal far-field boundary conditions. Figure 4(a) shows the comparison of surface temperatures obtained from our analytical solution (Eq. (12): lines) and FEM simulations (points). The agreement is excellent, including the detailed profile near and within the heater line itself.

FIG. 4.

FEM validation of closed form analytic solutions (Eqs. (12), (13), and (15)) for a general anisotropic sample. The figures compare analytical (lines) and FEM (points) results for (a) surface temperature at an intermediate frequency, ω = 2 × 104 rad/s, and (b) average heater line temperature over a wide frequency range. The dashed lines in (b) represent the low frequency limit of Eq. (13). Joule heating frequency is 2ω. k and C as in Fig. 2, with other parameters b = 5 μm and P′ = 3 W/m.

FIG. 4.

FEM validation of closed form analytic solutions (Eqs. (12), (13), and (15)) for a general anisotropic sample. The figures compare analytical (lines) and FEM (points) results for (a) surface temperature at an intermediate frequency, ω = 2 × 104 rad/s, and (b) average heater line temperature over a wide frequency range. The dashed lines in (b) represent the low frequency limit of Eq. (13). Joule heating frequency is 2ω. k and C as in Fig. 2, with other parameters b = 5 μm and P′ = 3 W/m.

Close modal

Second, to validate our solutions for T at arbitrary frequency, Eq. (13), and in the low frequency limit, Eq. (15), we performed FEM simulations spanning a wide frequency range. As shown in Figure 4(b), the results from FEM simulations (points) again agree very well with both analytical results (solid and dashed lines), thus further validating the analytical solutions.

In this section, we present two experimental schemes to use the low frequency solution of Eq. (15) to measure the thermal conductivity tensor of a substrate. The underlying focus is on the closed form solution to avoid nonlinear numerical fitting processes. Both approaches begin by using the slope method (Eq. (15), Figure 3) to find the thermal conductivity determinant in the natural xy coordinate system defined by the heater line and sample surface. Since rotation of the heater line around its longitudinal axis will give exactly the same measured slope, the schemes described below also use the magnitude of T to separate out the components of kij.

The first scheme requires knowledge of the sample’s volumetric heat capacity, C, and is depicted in the flowchart of Figure 5(a). In Eq. (15), we see that the magnitude of the in-phase temperature response at low frequency is independently related to both kxx and kdet. Thus, kdet from the slope method can be used with C to calculate kxx from the low frequency temperature response. This relationship is further clarified in Figure 5(b): if kxx is varied while keeping the determinant constant, then the magnitude of the in-phase T varies while its low-frequency slope remains constant.

FIG. 5.

(a) Process to determine kxx from the in-phase 3ω voltage measurements in the low frequency regime. The slope method gives k det = k xx k yy k xy 2 , which is used along with the in-phase magnitude and the sample’s heat capacity to determine kxx. (b) Effect of changing kxx while keeping kdet constant. For the in-phase curve at low frequency, the fixed determinant ensures constant slope while the magnitude depends on kxx. Simulation parameters as in Fig. 4.

FIG. 5.

(a) Process to determine kxx from the in-phase 3ω voltage measurements in the low frequency regime. The slope method gives k det = k xx k yy k xy 2 , which is used along with the in-phase magnitude and the sample’s heat capacity to determine kxx. (b) Effect of changing kxx while keeping kdet constant. For the in-phase curve at low frequency, the fixed determinant ensures constant slope while the magnitude depends on kxx. Simulation parameters as in Fig. 4.

Close modal

Figure 6 shows an example layout of heater lines to determine the elements of k. Heater lines 1 and 2 can be used to characterize the xy plane. In the low frequency limit, T for both heater lines will give the same slope value (determinant) but different magnitudes. Combining this information with C allows kxx and kyy to be calculated. It is then possible to extract the magnitude of the off-diagonal term, |kxy| from the determinant.

FIG. 6.

(a) To measure the thermal conductivity tensor for the xy plane of a given substrate, low frequency measurements from heater lines 1 and 2 can be used. (b) As shown in the flowchart, both lines will give the same slope T / ln ( ω ) , while their T magnitudes depend on the orientation of the heater line as well as C. Similarly, low frequency 3ω measurements from heater lines 1-4 on orthogonal faces in (a) can give the full 3-d thermal conductivity tensor. These approaches determine the magnitude, but not the sign of the off-diagonal terms such as kxy (see text).

FIG. 6.

(a) To measure the thermal conductivity tensor for the xy plane of a given substrate, low frequency measurements from heater lines 1 and 2 can be used. (b) As shown in the flowchart, both lines will give the same slope T / ln ( ω ) , while their T magnitudes depend on the orientation of the heater line as well as C. Similarly, low frequency 3ω measurements from heater lines 1-4 on orthogonal faces in (a) can give the full 3-d thermal conductivity tensor. These approaches determine the magnitude, but not the sign of the off-diagonal terms such as kxy (see text).

Close modal

To determine the sign of kxy, a third non-co-planar heater line can be used. Such a line is labeled γ in Figure 7, deposited on a face tilted by some angle θ with respect to face 1. Low frequency 3ω measurements from line γ will yield kxx, which can be related to the tensor elements in the xy coordinate system through standard tensor rotation rules,

k x x = k xx cos 2 ( θ ) + k xy sin ( 2 θ ) + k yy sin 2 ( θ ) .
(18)

Using Eq. (18), kxy now can be determined from the values of kxx, kyy, and kxx which are known from heater lines 1, 2, and γ, respectively. Note that lines 1 and 2 need not be on orthogonal faces for this method. Any set of parallel (z-oriented) heater lines on 3 non-co-planar surfaces will give a complete set of equations to solve for the 3 elements of the 2-d k tensor. A line at θ = 45° would maximize the sensitivity to kxy. Finally, multiple heater lines 1-4 on orthogonal faces will be required to find the full 3-d thermal conductivity tensor (Figure 6(a)) if only the magnitudes of the off-diagonal terms are needed, while additional heater lines on non-co-planar surfaces similar to line γ in Figure 7 can be used to determine their signs.

FIG. 7.

To fix the sign of kxy after determining k xy from Fig. 6, an additional heater line γ is required on a surface at some other angle. θ is defined as the angle between the x′ and x axes, measured in the counter-clockwise direction as shown. Using transformation rules (Eq. (18)), the effect of kxy is seen in kxx in the new x′ − y′ coordinate system. kxx can be measured using this new heater line, following the approach outlined in Figure 5.

FIG. 7.

To fix the sign of kxy after determining k xy from Fig. 6, an additional heater line γ is required on a surface at some other angle. θ is defined as the angle between the x′ and x axes, measured in the counter-clockwise direction as shown. Using transformation rules (Eq. (18)), the effect of kxy is seen in kxx in the new x′ − y′ coordinate system. kxx can be measured using this new heater line, following the approach outlined in Figure 5.

Close modal

The second scheme is differential and eliminates the heater linewidth and substrate heat capacity from the calculation, as long as their values are identical from sample to sample. Manipulating Eq. (15) and taking b and C as constant give

T line 2 ω 2 ξ 2 T line 1 ω 1 ξ 1 = ln ω 2 ω 1 × k yy k xx
(19)

and

T line 3 ω 3 ξ 3 T line 1 ω 1 ξ 1 = ln ω 3 ω 1 × k z z k xx × k xx k yy k x y 2 k z z k xx k z x 2 ,
(20)

where ξ n = T line n / ln ( ω ) , with n identifying the heater line orientation as per Figure 6. Here, ωn refers to any particular frequency used for the corresponding heater line’s voltage, and for convenience, one might use ω1 = ω2 = ω3. Since the thermal conductivity determinants in Eq. (20) are already known from slope measurements, it is now a direct calculation to find the ratio kxx/kyy using lines 1 and 2, and kxx/kzz using lines 1 and 3.

Note that this second scheme does not directly yield the off-diagonal elements of k, but they can be found using additional heater lines at different orientations. For example, using the differential scheme between heater lines 1 and γ in Figure 7 measures the ratio

k x x k xx = cos 2 ( θ ) + k xy k xx sin ( 2 θ ) + k yy k xx sin 2 ( θ ) .
(21)

Combining this additional information with the known values of kdet and kyy/kxx, it is a straightforward algebra to obtain all three of kxx, kyy, and kxy.

To demonstrate the major analytical results above, a self-consistency check was carried out with experiments on muscovite mica sheets, which was chosen for its moderate thermal conductivity anisotropy ratio at room temperature.27 To facilitate sample fabrication, we chose 0.25″ thick sheets (McMaster-Carr, product ID 8779K52). According to the vendor, the mica layers used to make these sheets are shaved directly from the natural mineral, layered horizontally, and then pressed into the required thickness. Mica is a silicate with a monoclinic crystal structure that is characterized by layered basal planes stacked in the c-axis direction of the unit cell. Thermal conductivity experiments27,28 on muscovite mica sheets usually presume an isotropic thermal conductivity in the ab plane, which is consistent with elastic moduli29–31 measurements that show a very weak or non-existent anisotropy in the ab plane of the silicate layers. Therefore, we assume the principal thermal conductivities Kab, Kab, and Kc, along the a-, b-, and c-axis directions, respectively. As shown in Figure 8, we aim to characterize a cross section of the layered material and prepare two samples accordingly. The first sample is the aligned case (Figure 8(a)), with Kc and Kab aligned perpendicular and parallel to the surface, and the 3ω experiments involve measuring these two principal values.

FIG. 8.

Schematics of the (a) aligned and (b) angled samples used for experiments. kθ is an effective thermal conductivity normal to the heater line, as defined in (Eq. (22)). In our experiments, the four different heater line orientations are 1 0 ° , 2 90 ° , 3 30 ° , 4 120 ° .

FIG. 8.

Schematics of the (a) aligned and (b) angled samples used for experiments. kθ is an effective thermal conductivity normal to the heater line, as defined in (Eq. (22)). In our experiments, the four different heater line orientations are 1 0 ° , 2 90 ° , 3 30 ° , 4 120 ° .

Close modal

The second sample is angled, to represent an arbitrary orientation rotated around the z axis. We introduce off-diagonal elements by machining orthogonal faces at angles θ and 90° + θ to the c direction (Figure 8(b)). For simplicity in presenting results from multiple orientations, for the rest of this section, it shall prove convenient to refer to the new coordinate system (x′, y′) as (θ, 90° + θ). Since θ is measured with respect to the c-axis of the material, for a surface of any θ, the corresponding kxx value is given straightforwardly by Eq. (18) with kxy = 0. In the notation of the present section, we denote kxx by kθ and have

k θ = K c cos 2 θ + K a b sin 2 θ .
(22)

Similarly, we shall denote the off-diagonal term kxy as kθ 90°+θ, which is easily shown to be

k θ 90 ° + θ = K a b K c sin θ cos θ .
(23)

Thus, for any rotated (x′, y′) coordinate system, the new k tensor is expressed as k θ k θ 90 ° + θ k θ 90 ° + θ k 90 ° + θ . The measured kij values in this coordinate system should match those calculated from Eq. (22) and (23). In our experiments, θ = 30°. To determine the off-diagonal term in the (θ, 90° + θ) coordinate system, we exploit the heater lines in the aligned configuration of Figure 8(a) by following the scheme described in Figure 7. The heater lines 1 and 2 in Figure 8(a) now correspond to parallel but non-co-planar heater lines at angles −30 and 60 in the (x′, y′) system of Figure 8(b). We can use either of these 2 lines now to measure the full thermal conductivity tensor in the (x′, y′) or (θ, 90° + θ) coordinate system.

All mica surfaces were polished using diamond lapping tools. A total of 16 heater lines were then deposited on several different samples, always at one of the four angles indicated in Figure 8, with the labeling (1 → 0°, 2 → 90°, 3 → 30°, 4 → 120°). Gold heater lines nominally 1.5 mm long, 60 μm wide, and 200 nm thick, with a 10 nm chromium adhesion layer, were patterned by evaporation through a shadow mask. The temperature coefficient of resistance was calibrated separately for every line. The 3ω measurements were performed in ambient air at a temperature of 27 ± 3 °C.

Invariance of the determinant. We first check the major theoretical prediction from (Eq. (15)) that the slope method for all heater lines should give the square root of the thermal conductivity determinant, kdet, which is expected to be a constant for all heater orientations since determinants are invariant upon rotation. This prediction is checked in Figure 9. For each orientation, we measured between 3 and 5 independent heater lines, with the resulting slope-values of thermal conductivity shown by the gray columns. The average for each orientation is indicated by the four blue solid lines. These averages indeed give nearly the same value, with a range from 0.79 to 0.86 W/m-K, in good agreement with the average from all 16 independent experiments (dashed line) of 0.81 W/m-K. This confirms the determinant is independent of heater line rotations about its axis. We also note that the overall standard deviation of ±0.10 W/m-K and standard error of the mean of ±0.02 W/m-K (where SE-Mean = SD/ ( N ) ) are much smaller than the observed differences between Kab and Kc as discussed later and shown in Figure 10.

FIG. 9.

Experimental confirmation of the invariance of thermal conductivity determinant, obtained using the slope method. Each column represents the measured kdet of a different microfabricated heater line (N = 16 total), grouped into four heater orientations. The aligned sample had θ = 0° and 90° while the angled sample which exercises an arbitrary anisotropy had θ = 30° and 120°. The average value for each orientation is depicted by the solid blue line, with the mean ± SE-mean values (where SE-mean = SD/ ( N ) , N being the number of samples) given on the plot. Finally, the average of all 16 measurements is given by the dashed black line. The agreement between the four orientations (range 0.79–0.86 W/m-K) is within the overall scatter of ±0.10 W/m-K (SD), confirming the constancy of the determinant to within noise.

FIG. 9.

Experimental confirmation of the invariance of thermal conductivity determinant, obtained using the slope method. Each column represents the measured kdet of a different microfabricated heater line (N = 16 total), grouped into four heater orientations. The aligned sample had θ = 0° and 90° while the angled sample which exercises an arbitrary anisotropy had θ = 30° and 120°. The average value for each orientation is depicted by the solid blue line, with the mean ± SE-mean values (where SE-mean = SD/ ( N ) , N being the number of samples) given on the plot. Finally, the average of all 16 measurements is given by the dashed black line. The agreement between the four orientations (range 0.79–0.86 W/m-K) is within the overall scatter of ±0.10 W/m-K (SD), confirming the constancy of the determinant to within noise.

Close modal
FIG. 10.

Effective thermal conductivity as a function of direction in mica samples. Blue circles show directly measured k-tensor elements. Kc and Kab are the principal values measured in the aligned samples, while k30°, k120°, and k30° 120° are measured in the angled samples. The calculated values (orange diamonds) for the angled sample are found by applying the tensor transformation rules of Eqs. (22) and (23) to the measured Kc and Kab. Error bars represent standard error of the mean, evaluated using measurements from multiple heater lines.

FIG. 10.

Effective thermal conductivity as a function of direction in mica samples. Blue circles show directly measured k-tensor elements. Kc and Kab are the principal values measured in the aligned samples, while k30°, k120°, and k30° 120° are measured in the angled samples. The calculated values (orange diamonds) for the angled sample are found by applying the tensor transformation rules of Eqs. (22) and (23) to the measured Kc and Kab. Error bars represent standard error of the mean, evaluated using measurements from multiple heater lines.

Close modal

Scheme 1:kelements. We now evaluate the individual tensor elements for both the aligned and angled samples. The first scheme (depicted in Figures 5 and 6) requires the volumetric heat capacity C = ρc, which was calculated using the measured gravimetric c (using differential scanning calorimetry: TA Instruments, Model Q20) and density ρ (using a geometric calculation of volume, confirmed by a water displacement method). The measured ρ of 2200 kg/m3 ± 1.5% is 18% lower than the literature,32 which may be due to imperfect pressing of the individual shaved layers in these industrial samples. The measured c is 804 J/kg - K ± 2.2% (average of 5 samples ± SE-mean), which is around 8% lower than reported for a single crystal specimen33 (876 J/kg - K); since we see no reason why c should have been affected by the sample manufacturing process, we use the average of these two values. The resulting volumetric C is 1.85 × 106 J/m3 - K at room temperature, which is around 20%-25% lower than reported previously for a single crystal specimen32–34 due mainly to the lower density. The low-frequency 3ω voltage data give the thermal conductivity perpendicular to the heater line, kθ. The resulting measurements are shown in Figure 10.

The first two values in Figure 10 are the measured principal thermal conductivities, Kc and Kab, extracted from heater line orientations 1 and 2, respectively. Using these principal values and the transformation rules of Eqs. (22) and (23), the tensor elements for the angled sample can be estimated, as shown by the orange diamonds. These are checked against direct measurements using heater lines 3 and 4 which give the diagonal elements, k30° and k120°, respectively (blue circles). Lines 3 and 4 can only measure k 30 ° 120 ° following the flowchart in Figure 6(b), and removing the sign ambiguity requires another additional line as noted in Figure 7. Here, we report the average of values obtained from the 2 sets of 3 non-co-planar heater line configurations, namely, 1 , 3 , 4 and 2 , 3 , 4 . Figure 10 shows that the calculated values of all three tensor elements for the angled sample are in good agreement with the directly measured elements, well within the standard errors of the means.

We see significant scatter in the thermal conductivity values reported in Figure 10. One obvious potential source of uncertainty is sample-to-sample variability since these are industrial grade samples synthesized from the natural mineral. We also note some more fundamental sensitivity issues which are intrinsic to this method, where the sensitivity of V3ω to a parameter β is defined following Gundrum et al.,35 

S β = ln V 3 ω ln β .
(24)

Consider this sensitivity to the various experimental parameters in the low frequency regime. Evaluation of Eq. (15) shows that V3ω is just as sensitive to C as it is to kθ, and twice as sensitive to b as to kθ. These samples exhibit a b variation of around ±10% which may be attributed to the simple shadow masked patterning process. We also note that the sensitivity of the slope to kdet is always 1, but over the frequency range of interest the sensitivity of V3ω to kθ is only around 0.2; therefore, it would not be surprising if the scatter (SE-mean) in Figure 10 were up to ∼5 times larger than in Figure 9.

Scheme 2: Ratios. Finally, to bypass the uncertainties in heater linewidth and heat capacity, we consider the second experimental scheme which uses a differential approach and deals in conductivity ratios. From the 7 lines deposited on the 2 aligned surfaces (heater line orientations 1 and 2), Eq. (20) gives Kab/Kc = 4.64 ± 0.48. Following this, we do two consistency checks. We first compare k120°/k30° measured directly by the differential method on the angled sample (Figure 8(b)) with the expected value of k120°/k30° from applying tensor transformation rules to Kab/Kc which gives

k 120 ° / k 30 ° = cos 2 120 ° + K a b / K c sin 2 120 ° cos 2 30 ° + K a b / K c sin 2 30 ° ,
(25)

yielding the estimate k120°/k30° = 1.95 ± 0.07. This compares favorably (14% difference) with the ratio measured directly by the differential method (Eq. (20)) on heater lines 3 and 4, which is k120°/k30° = 2.22 ± 0.25.

Similarly, for the off-diagonal element of the angled sample, tensor algebra gives

k 30 ° 120 ° / k 30 ° = cos 30 ° sin 30 ° K a b / K c 1 cos 2 30 ° + K a b / K c sin 2 30 ° ,
(26)

which using the measured Kab/Kc predicts k30°120°/k30° = 0.82 ± 0.06. This again compares favorably with the directly measured value k30°120°/k30° = 0.81 ± 0.21, which is the average value obtained by applying Eq. (21) to the set of lines {1,3,4} and {2,3,4}.

Finally, using the measured k30°120°/k30° and k120°/k30° from the differential scheme, and the average kdet from Fig. 9, we can evaluate the complete thermal conductivity tensor in the 30 ° , 120 ° coordinate system,

k = 0 . 65 0 . 53 0 . 53 1 . 45 W / m-K .
(27)

These tensor elements are systematically smaller than those obtained from the first scheme (last 3 blue circles of Fig. 10). The square root of the thermal conductivity tensor determinant, k det = k 30 ° k 120 ° k 30 ° 120 ° 2 , in Figure 10 is 1.05 ± 0.29 W/m - K, which is higher than kdet predicted from the slope method in Figure 9. This indicates possible errors in b and C, and the relatively low sensitivity to the thermal conductivity tensor elements in the first scheme as compared to the sensitivity to kdet in the slope method. On the other hand, the differential scheme which is independent of b and C works around uncertainties in measurements, thus confirming its advantage over the first scheme. This leads to a tighter scatter in the ratios reported by the differential method, which has inherently better sensitivity. Measuring kdet from the slope method, and the k-tensor element ratios from the differential scheme, ensures good sensitivity at each step and should yield the most accurate final results for the individual kij.

We report an analytic closed form solution for the 3ω heating problem on a substrate with arbitrary anisotropy in cartesian coordinates. The major analytical results are Eq. (13) and especially its low-frequency limit, Eq. (15). The extension of the familiar slope method gives the square root of the determinant of the thermal conductivity tensor, which is invariant upon rotation around the axis of the heater line. The magnitude of the in-phase temperature rise depends on this determinant as well as the effective thermal conductivity perpendicular to heater line, kxx in Figure 1 or kθ in Figure 8(b). The analytical solutions have been validated using FEM.

Two experimental schemes were devised to isolate the thermal conductivity tensor elements using multiple heater lines on orthogonal faces and demonstrated with experiments on industrial mica sheets. The first scheme can measure all the tensor elements, though it requires the volumetric heat capacity and gives only the absolute value of the off-diagonal elements. The second scheme is differential and does not need the heat capacity or heater linewidth, and gives anisotropy ratios. One additional heater line deposited on a non-orthogonal face (as in Figure 7) is required to fix the sign of the off-diagonal elements in the first scheme, and the specific tensor values rather than ratios in the second scheme. For our data set of 16 heater lines at 4 orientations, the first scheme over-predicts the thermal conductivity tensor elements, attributed to uncertainties in C and b values, while the differential scheme works around these issues.

To best combine the strengths and sensitivities of the various schemes, we recommend using the slope method to measure kdet and the differential scheme to measure kij ratios, information which can be combined to specify the complete and arbitrary tensor k.

We gratefully acknowledge support from a multidisciplinary research initiative (MRI) from the High Energy Lasers—Joint Technology Office (HEL-JTO) administered by the Army Research Office (ARO). The authors thank Sean Lubner, Geoff Wehmeyer, and Fan Yang for helpful discussions during the experiments.

Delta-function heater line

We present the details of the derivation beginning from Eq. (7). The equations are simplified by introducing the ratios εyy = kyy/kxx and εxy = kxy/kxx. Following a standard approach,23 we remove the y derivative by using the Fourier transform pairs T ̃ ( x , γ ) = e i γ y T ̄ ( x , y ) d y and T ̄ ( x , y ) = 1 2 π e i γ y T ̃ ( x , γ ) d γ . After applying this transform and Eq. (7) to the governing Eq. (3), we have

2 T ̃ 2 x 2 i γ ε xy T ̃ x γ 2 ε yy T ̃ + P k xx e i γ y δ ( x , y ) d y = i 2 ω α xx T ̃ ,
(A1)

where αxx = kxx/C is the thermal diffusivity in the x − direction. Next, the first-order x-derivative in Eq. (A1) can be removed by introducing

T ̃ ( x , γ ) = w ( x , γ ) e i γ ε xy x ,
(A2)

leading to

2 w 2 x γ 2 ( ε yy ε xy 2 ) w + e i γ ε xy x P k xx e i γ y δ ( x , y ) d y = i 2 ω α xx w .
(A3)

Similarly transforming the boundary conditions leads to

w x = 0 at x = 0 ,
(A4)
w x = 0 and w = 0  at  x .
(A5)

Note that for the boundary condition at x = 0 for this transformed temperature, the eigenfunction for the homogenous heat conduction problem is cos(βx).

The next step is to remove the x-derivatives using another pair of integral transforms: ζ ( β , γ ) = 0 cos ( β x ) w ( x , γ ) d x and w ( x , γ ) = 2 π 0 cos ( β x ) ζ ( β , γ ) d β . Thus, the solution in terms of transform variables is

ζ ( β , γ ) = P / k xx i 2 ω / α xx + β 2 + γ 2 ( ε yy ε xy 2 ) .
(A6)

After performing the required inverse transforms, the integral form in terms of x and y is

T ¯ ( x , y ) = 2 P π 2 k xx o 0 cos ( β x ) cos ( γ ( y ε xy x ) ) q xx 2 + β 2 + γ 2 ( ε yy ε xy 2 ) d β d γ ,
(A7)

with q xx = i 2 ω C / k xx defined in Eq. (9).

An important check of Eq. (A7) is to simplify it for an isotropic sample and verify it recovers the well-known closed-form solution given by Cahill in 1990.11 Setting εyy = 1 and εxy = 0, we see that this is indeed true,

T ¯ ( x , y ) = 2 P π 2 k o 0 cos ( β x ) cos ( γ y ) q 2 + β 2 + γ 2 d β d γ , = P π k K 0 q x 2 + y 2 ,
(A8)

where K0 is the modified Bessel function of the second kind, and the integrals were evaluated using symbolic math software (MAPLE).

Returning to the full solution of Eq. (A7), it is fruitful to make the substitutions γ 2 = γ 2 ( ε yy ε xy 2 ) and y = y ε xy x ε yy ε xy 2 . We have

T ¯ ( x , y ) = 2 P π 2 k xx ε yy ε xy 2 o 0 cos ( β x ) cos ( γ y ) q x x 2 + β 2 + γ 2 d β d γ .
(A9)

Thus, after this transformation, the anisotropic problem’s integral of Eq. (A9) is now formally equivalent to the isotropic problem’s integral of Eq. (A8). Therefore, the solution to Eq. (A9) follows immediately, and returning to the original (untransformed) variables, we have the solution for the δ-function line heater as

T ¯ ( x , y ) = P π k xx ε yy ε xy 2 K 0 q xx ε yy x 2 + y 2 2 ε xy x y ε yy ε xy 2 ,
(A10)

which is equivalent to Eq. (8) in the main text.

Finite-width heater line

We now obtain solutions for the finite heater-width case, including the surface temperature field and average temperature of the heater line. Through suitable substitutions, we can again transform the anisotropic problem into an equivalent isotropic one, thereby exploiting the exact analytical solutions obtained by Duquesne et al.24 for the isotropic case.

We begin here from Eq. (11) of the main text. Scale b, y, and y′ such that for any variable η, its scaled counterpart is

η ̃ = η k xx k det ,
(A11)

resulting in

T surf ( y ̃ ) = P 2 b ̃ π k xx k yy k xy 2 b ̃ b ̃ K 0 q xx y ̃ y ̃ d y ̃ .
(A12)

We have again succeeded in converting the relevant integral from an anisotropic to isotropic formulation. An equivalent isotropic integral was formulated in Fourier space by Cahill in 1990,11 though it was 20 years before the analytical solution was first obtained by Duquesne et al.24 We can adapt their isotropic solutions directly. Thus, the closed form solution for surface temperature caused by a finite-width heater line on an anisotropic substrate is written down directly as Eq. (12).

Proceeding similarly for the average temperature of the heater line on the anisotropic substrate, we write

T = P 4 b 2 π k xx k yy k xy 2 0 2 b 0 2 b K 0 q xx k xx u v k xx k yy k xy 2 d u d v .
(A13)

After scaling u, v, and b using Eq. (A11), we again obtain an integral formally identical to one solved by Duquesne et al.,24 

T = P 4 b ̃ 2 π k xx k yy k xy 2 0 2 b ̃ 0 2 b ̃ K 0 q xx u ̃ v ̃ d u ̃ d v ̃ .
(A14)

Finally, after suitably rewriting the isotropic result,24 the average heater line temperature for our anisotropic substrate can be expressed in terms of the Meijer-G function, resulting in Eq. (13) of the main text.

We consider the temperature field caused by a δ-line heater on an arbitrarily anisotropic substrate, given by Eq. (A10). The locus for isothermal lines was recognized in Eq. (10), which is reminiscent of the general form of an ellipse inclined at an angle ϕ from the x axis, given by

x cos ( ϕ ) + y sin ( ϕ ) 2 a 2 + x sin ( ϕ ) y cos ( ϕ ) 2 b 2 = const .
(B1)

Matching the coefficients of x2, y2, and xy in Eqs. (10) and (B1) gives 3 equations,

cos 2 ( ϕ ) a 2 + sin 2 ( ϕ ) b 2 = ε yy ,
(B2)
sin 2 ( ϕ ) a 2 + cos 2 ( ϕ ) b 2 = 1 ,
(B3)
2 cos ( ϕ ) sin ( ϕ ) 1 a 2 1 b 2 = 2 ε xy ,
(B4)

which can be solved for a, b, and ϕ,

a 2 = cos 4 ( ϕ ) sin 4 ( ϕ ) ε yy cos 2 ( ϕ ) sin 2 ( ϕ ) ,
(B5)
b 2 = cos 4 ( ϕ ) sin 4 ( ϕ ) cos 2 ( ϕ ) ε yy sin 2 ( ϕ ) ,
(B6)
tan ( 2 ϕ ) = 2 ε xy ε yy 1 = 2 k xy k yy k xx .
(B7)

This expression for ϕ is exactly the same as the known result36 for the angle between the principal conductivity of a material and the surface normal (the x-direction in Fig. 1). Thus, this result proves that the isotherms for a δ-line heater on a semi-infinite substrate are ellipses whose major and minor axes are exactly aligned with the principal thermal conductivities of the material, regardless of the rotation of the free surface around the z axis.

1.
M.
Carle
,
P.
Pierrat
,
C.
Lahalle-Gravier
,
S.
Scherrer
, and
H.
Scherrer
,
J. Phys. Chem. Solids
56
,
201
(
1995
).
2.
X.
Yan
,
B.
Poudel
,
Y.
Ma
,
W. S.
Liu
,
G.
Joshi
,
H.
Wang
,
Y.
Lan
,
D.
Wang
,
G.
Chen
, and
Z. F.
Ren
,
Nano Lett.
10
,
3373
(
2010
).
3.
L.-D.
Zhao
,
S.-H.
Lo
,
Y.
Zhang
,
H.
Sun
,
G.
Tan
,
C.
Uher
,
C.
Wolverton
,
V. P.
Dravid
, and
M. G.
Kanatzidis
,
Nature
508
,
373
(
2014
).
4.
D.
Morelli
,
G.
Doll
,
J.
Heremans
,
M.
Dresselhaus
,
A.
Cassanho
,
D.
Gabbe
, and
H.
Jenssen
,
Phys. Rev. B
41
,
2520
(
1990
).
5.
S.
Shen
,
A.
Henry
,
J.
Tong
,
R.
Zheng
, and
G.
Chen
,
Nat. Nanotechnol.
5
,
251
(
2010
).
6.
S. C.
Glotzer
and
M. J.
Solomon
,
Nat. Mater.
6
,
557
(
2007
).
7.
Y.
Gaillard
,
V. J.
Rico
,
E.
Jimenez-Pique
, and
A. R.
González-Elipe
,
J. Phys. D: Appl. Phys.
42
,
145305
(
2009
).
8.
K.
Robbie
,
J. Vac. Sci. Technol., A
15
,
1460
(
1997
).
9.
J. L.
Janning
,
Appl. Phys. Lett.
21
,
173
(
1972
).
10.
J. P.
Feser
,
J.
Liu
, and
D. G.
Cahill
,
Rev. Sci. Instrum.
85
,
104903
(
2014
).
11.
D. G.
Cahill
,
Rev. Sci. Instrum.
61
,
802
(
1990
).
12.
C.
Dames
and
G.
Chen
,
Rev. Sci. Instrum.
76
,
124902
(
2005
).
13.
D.
Cahill
,
M.
Katiyar
, and
J.
Abelson
,
Phys. Rev. B
50
,
6077
(
1994
).
14.
J. H.
Kim
,
A.
Feldman
, and
D.
Novotny
,
J. Appl. Phys.
86
,
3959
(
1999
).
15.
T.
Borca-Tasciuc
,
A. R.
Kumar
, and
G.
Chen
,
Rev. Sci. Instrum.
72
,
2139
(
2001
).
16.
A.
Feldman
,
High Temp. - High Pressures
31
,
293
(
1999
).
17.
T.
Tong
and
A.
Majumdar
,
Rev. Sci. Instrum.
77
,
104902
(
2006
).
18.
A. T.
Ramu
and
J. E.
Bowers
,
J. Appl. Phys.
112
,
043516
(
2012
).
19.
A. T.
Ramu
and
J. E.
Bowers
,
Rev. Sci. Instrum.
83
,
124903
(
2012
).
20.
R. E.
Newnham
,
Properties of Materials: Anisotropy, Symmetry, Structure
(
Oxford University Press Inc.
,
New York
,
2005
), p.
34
.
21.
L.
Onsager
,
Phys. Rev.
37
,
405
(
1931
).
22.
H. S.
Carslaw
and
J. C.
Jaeger
,
Conduction of Heat in Solids
, 2nd ed. (
Oxford University Press
,
London
,
1959
), p.
261
.
23.
D. W.
Hahn
and
M. N.
Özişik
,
Heat Conduction
, 3rd ed. (John Wiley & Sons, Inc., 2012), p. 614.
24.
J.-Y.
Duquesne
,
D.
Fournier
, and
C.
Frétigny
,
J. Appl. Phys.
108
,
086104
(
2010
).
25.
A.
Milton
and
I. A.
Stegun
,
Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables
, 9th ed. (Dover Publications Inc., New York, 1970), pp. 374 and 498.
26.

The Meijer-G function can be evaluated in Matlab using the MeijerG Function in the Mupad Notebook interface and is also available in Maple and Mathematica.

27.
A. S.
Gray
and
C.
Uher
,
J. Mater. Sci.
12
,
959
(
1977
).
28.
V.
Cermark
and
L.
Rybach
, in
Physics Properties of Rocks
, edited by
G.
Angenheister
(
Springer-Verlag GmbH
,
Heigelberg
,
1982
), p.
314
.
29.
R. E.
Newnham
,
Mineral Mag.
39
,
78
(
1973
).
30.
M. T.
Vaughan
and
S.
Guggenheim
,
J. Geophys. Res.
91
,
4657
(
1986
).
31.
P. Y.
Cholach
and
D. R.
Schmitt
,
J. Geophys. Res.: Solid Earth
111
,
1
(
2006
).
32.
J.
Faust
and
E.
Knittle
,
J. Geophys. Res.
99
,
19785
(
1994
).
33.
K. M.
Krupka
,
R. A.
Robie
, and
B. S.
Hemingway
,
Am. Mineral.
64
,
86
(
1979
).
34.
W.-P.
Hsieh
,
B.
Chen
,
J.
Li
,
P.
Keblinski
, and
D. G.
Cahill
,
Phys. Rev. B
80
,
180302
(
2009
).
35.
B.
Gundrum
,
D.
Cahill
, and
R.
Averback
,
Phys. Rev. B
72
,
245426
(
2005
).
36.
J. F.
Nye
,
Physical Properties of Crystals Their Representations by Tensors and Matrices
(
Oxford University Press
,
1985
), p.
45
.