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.

## INTRODUCTION

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ω method^{11,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 Majumdar^{17} 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 Bowers^{18} 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.

## THE MODEL SETUP

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

and the most general anisotropic thermal conductivity tensor is

**k** is real and symmetric, and it is known^{21} that *k _{ii}* > 0 and

*k*>

_{ii}k_{jj}*k*

_{ij}^{2}. 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 algorithm^{16} 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.

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

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

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

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

where $i= \u2212 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 transforms^{23} 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 \xaf $ is

where *K*_{0} is the modified Bessel function of the second kind and

where $ q xx $ represents the inverse of the penetration depth in the *x*-direction. As a check, when *k _{xx}* =

*k*and

_{yy}*k*= 0, the familiar isotropic solution from Cahill

_{xy}^{11}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

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 *x* − *y* 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.

## FINITE HEATER WIDTH AND AVERAGE HEATER TEMPERATURE

The next step in this analysis is to find the surface temperature profile, *T _{surf}*(

*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,

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

where $\Omega = b 2 2 \omega C k xx k yy \u2212 k xy 2 $, $Y= y b $, and Ξ(*y*) = *K*_{0}(*y*) *L*_{−1}(*y*) + *K*_{1}(*y*) *L*_{0}(*y*), where *K _{n}* and

*L*are the n-th order modified Bessel function of the second kind and Struve function,

_{n}^{25}respectively. Similarly, the average of

*T*over the heater linewidth, $ T $, can also be found by using a change of variables to adapt the known isotropic solution

_{surf}^{24}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 $ as

^{24,26}

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

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

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, *K*_{0}(*z*) ≈ − ln(*z*) − *γ* + ln(2), we find the low-frequency limit of Eq. (13) as

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, $\xi =\u2202 T /\u2202ln \omega $, is inversely proportional to

with the proportionality constant −*P*′/2*π*. Importantly, we recognize *k _{det}* 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

*k*,

_{xx}*k*, and

_{yy}*k*. 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.

_{xy}Finally, the RMS 3ω voltage is related to the RMS 2ω temperature by^{11,12}

Here, *I*_{1ω} 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 *V*_{3ω} and $ T $.

## NUMERICAL VALIDATION

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 *ω* < 10^{3} rad/s and *ω* > 4 × 10^{5} 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 × 10^{4} 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.

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.

## EXPERIMENTAL SCHEMES AND PROOF-OF-CONCEPT

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 *x* − *y* 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 *k _{ij}*.

### Scheme 1

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 *k _{xx}* and

*k*

_{det}. Thus,

*k*

_{det}from the slope method can be used with

*C*to calculate

*k*from the low frequency temperature response. This relationship is further clarified in Figure 5(b): if

_{xx}*k*is varied while keeping the determinant constant, then the magnitude of the in-phase $ T $ varies while its low-frequency slope remains constant.

_{xx}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 *x* − *y* 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 *k _{xx}* and

*k*to be calculated. It is then possible to extract the magnitude of the off-diagonal term, |

_{yy}*k*| from the determinant.

_{xy}To determine the sign of *k _{xy}*, 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

*k*

_{x′x′}, which can be related to the tensor elements in the

*x*−

*y*coordinate system through standard tensor rotation rules,

Using Eq. (18), *k _{xy}* now can be determined from the values of

*k*,

_{xx}*k*, and

_{yy}*k*

_{x′x′}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

*k*. 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

_{xy}*γ*in Figure 7 can be used to determine their signs.

### Scheme 2

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

and

where $ \xi n =\u2202 T line \u2212 n /\u2202ln ( \omega ) $, 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

*k*/

_{xx}*k*using lines 1 and 2, and

_{yy}*k*/

_{xx}*k*using lines 1 and 3.

_{zz}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

Combining this additional information with the known values of *k*_{det} and *k _{yy}*/

*k*, it is a straightforward algebra to obtain all three of

_{xx}*k*,

_{xx}*k*, and

_{yy}*k*.

_{xy}## EXPERIMENTS WITH MICA

### Concept

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 experiments^{27,28} on muscovite mica sheets usually presume an isotropic thermal conductivity in the ab plane, which is consistent with elastic moduli^{29–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 *K _{ab}*,

*K*, and

_{ab}*K*, along the

_{c}*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

*K*and

_{c}*K*aligned perpendicular and parallel to the surface, and the 3ω experiments involve measuring these two principal values.

_{ab}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 *k*_{x′x′} value is given straightforwardly by Eq. (18) with *k _{xy}* = 0. In the notation of the present section, we denote

*k*

_{x′x′}by

*k*

_{θ}and have

Similarly, we shall denote the off-diagonal term *k*_{x′y′} as *k*_{θ 90°+θ}, which is easily shown to be

Thus, for any rotated (*x*′, *y*′) coordinate system, the new **k** tensor is expressed as $ k \theta k \theta 90 \xb0 + \theta k \theta 90 \xb0 + \theta k 90 \xb0 + \theta $. The measured *k _{ij}* 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.

## RESULTS

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, *k*_{det}, 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 *K _{ab}* and

*K*as discussed later and shown in Figure 10.

_{c}*Scheme 1:* *k* *elements*. 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/m^{3} ± 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 specimen^{33} (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 × 10^{6} J/m^{3} - K at room temperature, which is around 20%-25% lower than reported previously for a single crystal specimen^{32–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, *K _{c}* and

*K*, 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,

_{ab}*k*

_{30°}and

*k*

_{120°}, respectively (blue circles). Lines 3 and 4 can only measure $ k 30 \xb0 120 \xb0 $ 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 *V*_{3ω} to a parameter *β* is defined following Gundrum *et al.*,^{35}

Consider this sensitivity to the various experimental parameters in the low frequency regime. Evaluation of Eq. (15) shows that *V*_{3ω} 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 *k*_{det} is always 1, but over the frequency range of interest the sensitivity of *V*_{3ω} 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 *K _{ab}*/

*K*= 4.64 ± 0.48. Following this, we do two consistency checks. We first compare

_{c}*k*

_{120°}/

*k*

_{30°}measured directly by the differential method on the angled sample (Figure 8(b)) with the expected value of

*k*

_{120°}/

*k*

_{30°}from applying tensor transformation rules to

*K*/

_{ab}*K*which gives

_{c} yielding the estimate *k*_{120°}/*k*_{30°} = 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 *k*_{120°}/*k*_{30°} = 2.22 ± 0.25.

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

which using the measured *K _{ab}*/

*K*predicts

_{c}*k*

_{30°120°}/

*k*

_{30°}= 0.82 ± 0.06. This again compares favorably with the directly measured value

*k*

_{30°120°}/

*k*

_{30°}= 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 *k*_{30°120°}/*k*_{30°} and *k*_{120°}/*k*_{30°} from the differential scheme, and the average *k*_{det} from Fig. 9, we can evaluate the complete thermal conductivity tensor in the $ 30 \xb0 , 120 \xb0 $ coordinate system,

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 \xb0 k 120 \xb0 \u2212 k 30 \xb0 120 \xb0 2 $, in Figure 10 is 1.05 ± 0.29 W/m - K, which is higher than *k _{det}* 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

*k*

_{det}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

*k*

_{det}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

*k*.

_{ij}## SUMMARY

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, *k _{xx}* 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 *k*_{det} and the differential scheme to measure *k _{ij}* ratios, information which can be combined to specify the complete and arbitrary tensor

**k**.

## Acknowledgments

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.

### APPENDIX A: CLOSED FORM SOLUTION FOR THE 3ω PROBLEM ON AN ANISOTROPIC SUBSTRATE

#### Delta-function heater line

We present the details of the derivation beginning from Eq. (7). The equations are simplified by introducing the ratios ε_{yy} = *k _{yy}*/

*k*and ε

_{xx}_{xy}=

*k*/

_{xy}*k*. Following a standard approach,

_{xx}^{23}we remove the

*y*derivative by using the Fourier transform pairs $ T \u0303 ( x , \gamma ) = \u222b \u2212 \u221e \u221e e i \gamma y \u2032 T \u0304 ( x , y \u2032 ) d y \u2032 $ and $ T \u0304 ( x , y ) = 1 2 \pi \u222b \u2212 \u221e \u221e e \u2212 i \gamma y T \u0303 ( x , \gamma ) d\gamma $. After applying this transform and Eq. (7) to the governing Eq. (3), we have

where *α _{xx}* =

*k*/

_{xx}*C*is the thermal diffusivity in the

*x*− direction. Next, the first-order x-derivative in Eq. (A1) can be removed by introducing

leading to

Similarly transforming the boundary conditions leads to

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: $\zeta ( \beta , \gamma ) = \u222b 0 \u221e cos ( \beta x \u2032 ) w ( x \u2032 , \gamma ) d x \u2032 $ and $w ( x , \gamma ) = 2 \pi \u222b 0 \u221e cos ( \beta x ) \zeta ( \beta , \gamma ) d\beta $. Thus, the solution in terms of transform variables is

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

with $ q xx = i 2 \omega 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,

where *K*_{0} 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 $ \gamma \u2032 2 = \gamma 2 ( \epsilon yy \u2212 \epsilon xy 2 ) $ and $ y \u2032 = y \u2212 \epsilon xy x \epsilon yy \u2212 \epsilon xy 2 $. We have

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

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

resulting in

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

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

### APPENDIX B: MAJOR AXES ORIENTATION OF ELLIPTICAL ISOTHERMS

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

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

This expression for *ϕ* is exactly the same as the known result^{36} 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.

## REFERENCES

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.