Effective medium theory for the low-temperature heat capacity of a metasolid plate

Abstract

Nanopatterning can be used to strongly control the thermal properties of solids, but theoretical understanding relies often on complex numerical simulations. Here, an analytical theory is derived for the low temperature heat capacity of a nanopatterned phononic crystal plate, focusing on the geometry of a square lattice of cylindrical holes in an isotropic matrix material. Its quasistatic elastic properties were studied using an anisotropic effective medium theory, that is, considering it as a homogenized metasolid. The effective elastic parameters can then be used as an input for an anisotropic plate theory, yielding analytical expressions for the dispersion relations of the three lowest phonon modes that are dominant in the low temperature limit below 1K. Those results were then used to derive a simple analytical formula for the heat capacity, which was compared numerically with the exact results for an example material. The effects of material and geometric design parameters in the formula are also discussed, giving simple guidelines how to tune the heat capacity up to an order of magnitude or more.

Introduction

Recent decade has seen an expansion of the range of metamaterials — with properties not found in nature—to the mechanical realm, from acoustic wave control in fluids¹ and mechanics of solid materials^2,3 to elastic waves in solids⁴. Typically these metamaterials are engineered from a microstructure, which strongly modifies the overall elastic response of the material body, which can often be analyzed in terms of effective materials parameters such as the effective elastic moduli and mass density. In the area of wave motion, such homogenization is useful when the wavelengths in question are larger than the length scales of the microstructure. Development of effective medium theories (EMTs) has also been an important topic in material physics and engineering in general, for the understanding the elastic behavior of composites and layered materials⁵. Research has progressed from specific cases of materials with random or periodic solid or void inclusions, to general periodic three-dimensional (3D) structures of full anisotropic elasticity^6,7.

On the other hand, vibrations in materials are also excited thermally, with phonons as the quanta of such vibrations. In the full quantum mechanical picture, phonons consist of atomic level vibrational degrees of freedom, consisting of acoustic and in most cases also optical phonon modes. Continuum elasticity does not describe the high-frequency optical modes or the high-wave-vector dispersion of the acoustic modes, so that it cannot be applied to study thermal phonons in general. However, at low temperatures, only the low-frequency long-wavelength modes are excited, meaning that continuum elasticity is valid description even for thermal phonons. Therefore, at low temperatures it is possible to have a strong, orders of magnitude control of thermal conductance^8,9,10 or heat capacity¹¹ by modifying the acoustic phonon dispersion relations through micro and nanostructuring of the material. Such an approach has experimentally been shown to be possible for thermal conductance using periodic two-dimensional (2D) holey phononic crystals (PnCs) in the nanoscale^12,13,14, relying on strong interference effects due to Bragg scattering, which start to operate for wavelengths of the order of the periodicity of the structure. At higher temperatures, controlling thermal conductance with PnCs rely more on diffusive scattering from surfaces¹⁵. Another possibility is to use pillared structures, which add effects of local resonances of the pillars into the mix¹⁶.

In contrast, here we want to consider whether the thermal properties of such nanostructures can be understood from the perspective of an elastic metamaterial (Fig. 1). In other words, we consider to what extent homogenized, effective material parameters (elastic stiffness tensor, density) can be derived and used to predict thermal properties in the limit where the dominant thermal wavelength is larger than the size scale of the nanostructure. For the case study reported here, we concentrate on the heat capacity of a plate with isotropic material parameters, perforated with cylindrical circular holes (Fig. 1). In addition, we will try to answer the question to what extent such a thermal metamaterial exhibits properties not observed in a naturally occurring solid.

Fig. 1: Schematic of the process of homogenization of a holey phononic crystal plate.

The real holey material with an isotropic elastic matrix material (parameters λ, μ, ρ) will be considered as a continuous homogenized but anisotropic metasolid with effective elastic parameters ${c}_{ij}^{* }$,ρ*.

The driving motivation for developing such an effective medium theory is simplicity and physical insight. Without such a theory, the analysis of thermal properties of a PnC plate structure requires extensive numerical computations of millions of data points for the phonon dispersion data using for example the finite element method^8,11,12. Complexities arise because (i) for PnC plates the displacement fields are more complicated than bulk plane waves, and thus have to be solved by full 3D wave models; and (ii) computation of thermal properties require a fine sampling of the wave vector-space in all 2D directions in the plane of the plate. Thus, simultaneous fine tuning and optimization of geometric and physical design parameters by direct computation is challenging.

Here, we develop an analytical effective medium theory for the low temperature heat capacity of holey phononic crystal metasolid and show that it is accurate in the sub-wavelength limit where the dominant thermal wavelength is larger than the PnC period. We use a quasistatic homogenization procedure based on the average eigenstrain approximation¹⁷, which is shown to be valid in the low-temperature limit by comparing to the direct numerical solution. This approach allows us to study the emergent anisotropy of the effective elastic parameters of a square lattice PnC plate as a function of the hole filling factor, in contrast to previous effective medium theories for PnC plates, which demonstrated results only for hexagonal lattices exhibiting in-plane isotropy^18,19,20. Analytical expressions valid in the long wavelength limit for all three low frequency anisotropic plate (Lamb wave) modes (flexural, dilatational and shear) are also given as a function of both the magnitude and the azimuthal angle of the wave vector.

Results and discussion

Effective medium theory for a phononic crystal in the low-frequency limit

For the purpose of understanding how a phononic crystal plate geometry affects its heat capacity in the low temperature limit, it suffices to consider only the three lowest-frequency acoustic phonon modes. In the long wavelength limit, these can be categorized based on the standard plate (Lamb wave) theory²¹ into an in-plane horizontal shear mode (SH mode), a flexural or antisymmetric mode (A mode), and a symmetric mode (S mode), also known as longitudinal or dilatational mode. (The symmetry and antisymmetry refers to the reflection symmetry of the modes in the mid-plane of the plate.) We will use an analytic formalism to express the low-frequency dispersive properties of these modes which is simple enough so that a closed form expression for the heat capacity of holey phononic crystal plates can be obtained.

In the case of PnC plates with circular hole patterns in a 2D square lattice (Fig. 1), the obvious complication to the analytical treatment is the periodic pattern itself. The exact mathematical description for the low-frequency phonon modes in such a geometry depends on the elasticity tensor c and the mass density ρ of the plate material, and is based on the solution of the elastic wave equation (details in the Methods section), which only allows numerical solutions using e.g., the finite element method.

However, with the use of an elastic EMT, the effect of the hole pattern can be incorporated in the effective elasticity tensor of the material. One of the simplest models available for periodic media in the long wavelength limit is the average eigenstrain approximation²², which describes the material in the macroscopic scale, whereas the voids are assumed microscopic. In the average eigenstrain approximation the problem domain including the holes is replaced by a continuous domain without the holes (Fig. 1), with a new effective elasticity tensor c* and effective mass density ρ*. The new effective tensor c* (Eq. (11), Methods) depends only on c, the lattice geometry, the shape of the holes and their volume filling factor f, and the Poisson ratio of the matrix material ν, whereas ρ* = (1 − f)ρ. Often the effective material represented by c* has a higher anisotropy than the original material c, which is the case in our example as well.

Our homogenization approach is then the following: we first consider a periodic isotropic 3D material containing infinitely long cylindrical voids forming a 2D square lattice, for which the calculation of the effective tensor c* is particularly simple and can be done analytically (Methods). (For an anisotropic matrix material, only numerical solutions exist⁶.) After obtaining the effective tensor c* for such a 3D effective medium, we imagine cutting a slice out of it in the direction perpendicular to the axes of the voids, and use Lamb-wave theory for anisotropic materials²³ to write asymptotic approximations of the dispersion relations for the lowest plate modes.

In Fig. 2, we plot the effective elastic moduli as a function of the hole filling factor f, calculated based on the above effective medium theory for the case of infinitely long cylindrical holey PnC, with the cylinder axes aligned with the x₃-coordinate direction. The material of the matrix in all the examples in this article is isotropic silicon nitride (Young’s modulus E = 250 GPa, density ρ = 3100 kg/m³, Poisson ratio ν = 0.23). The infinite sums of Eq. (14) were calculated up to order ∣n_i∣ ≤ 2000. We see how the material naturally becomes more compliant (c* elements decrease) by a factor of 2-3 as the holes take up more volume. Anisotropy is seen to emerge as some of the components that are equal in the f → 0 limit split. Later, we will show that for the high-end range f > 0.5 the EMT theory used here becomes less accurate.

Fig. 2: Effective elastic stiffness tensor components of a SiN PnC metasolid as a function of the hole filling factor f.

We have used the average eigenstrain approximation for a square lattice of infinitely long cylindrical holes aligned in the x₃-direction.

The full boundary value problem (Methods) does not have closed form expressions for the solutions. However, asymptotic approximations can be given for all three modes in the low frequency limit. These approximations have been available for a while even for the most general anisotropic plates²³. In ref. ²³, using Stroh’s formalism, a plane wave ansatz for phonon displacement fields and normal tractions in the plate was converted to a system of six linear differential equations, from which approximations to the propagator of the coefficient matrix (Stroh matrix) led to asymptotic presentations of the dispersion relations. We apply the formulas from ref. ²³ to our effective medium and obtain for the quadratic antisymmetric Lamb (flexural) mode

$${\omega }_{{{{{{{{\rm{A}}}}}}}}}({{{{{{{\bf{k}}}}}}}})=\frac{{k}^{2}d}{2\sqrt{3{\rho }^{* }}}\sqrt{{{{{{{{\mathcal{B}}}}}}}}-\frac{1}{2}{{{{{{{\mathcal{A}}}}}}}}\,{\sin }^{2}\,2\phi },$$

(1)

while for the two remaining non-dispersive symmetric Lamb and the horizontal shear modes we get

$$\begin{array}{ccc}{\omega }_{{{{{{{{{\rm{S}}}}}}}}}_{+},{{{{{{{{\rm{SH}}}}}}}}}_{-}}({{{{{{{\bf{k}}}}}}}})=\frac{k}{\sqrt{2{\rho }^{* }}}\sqrt{{{{{{{{\mathcal{B}}}}}}}}+{c}_{66}^{* }\pm \sqrt{{({{{{{{{\mathcal{B}}}}}}}}-{c}_{66}^{* })}^{2}-[2{{{{{{{\mathcal{A}}}}}}}}({{{{{{{\mathcal{B}}}}}}}}-{c}_{66}^{* })-{{{{{{{{\mathcal{A}}}}}}}}}^{2}]\,{\sin }^{2}\,2\phi }},\end{array}$$

(2)

where the plus sign corresponds to to the symmetric mode (S) and the minus sign to the shear mode (SH).

In the above formulas we used the polar coordinates k = (k, ϕ) for the 2D in-plane wave vector, and notations

$${{{{{{{\mathcal{A}}}}}}}} ={c}_{11}^{* }-{c}_{12}^{* }-2{c}_{66}^{* }\\ {{{{{{{\mathcal{B}}}}}}}} ={c}_{11}^{* }-{({c}_{13}^{* })}^{2}/{c}_{33}^{* }.$$

(3)

Constant ${{{{{{{\mathcal{A}}}}}}}}$ is called the coefficient of elastic anisotropy which is zero for isotropic materials. From Eqs.(1) and (2) we clearly see that if ${{{{{{{\mathcal{A}}}}}}}}=0$, the ϕ-dependent terms drop out. The effective density is geometrically calculated from the amount of missing material, ρ* = (1 − f)ρ.

To demonstrate the effect of the filling factor on the lowest wave modes, we plot in Fig. 3 the slowness surfaces k/ω for the A, S and SH modes from Eqs. (1), (2), for three different f = 0, 0.25, 0.5 at ω = 10¹⁰ rad/s, which is a frequency dominantly populated at 0.1 K¹². (The slowness of the linear S and SH modes naturally do not depend on frequency, but for the flexural A-mode it does.) As the hole fraction increases from zero, the wave speeds for all modes become slower and anisotropic. The strongest effect happens for the SH mode, whose anisotropy (ratio of max to min slowness) reaches ~ 2. It is also notable that the A and S modes are predominantly slowed in the diagonal directions unlike the SH mode, which slows down mostly along the unit cell axes.

Fig. 3: Radial plot of the slowness surfaces (k/ω vs ϕ) of the three lowest modes of the holey SiN PnC metaplate.

The lowest modes (A, S, and SH) are shown for three filling factors f = 0, 0.25, 0.5. For the flexural A-mode, slowness is a function of frequency, here plotted for a thermally relevant frequency ω = 10¹⁰ rad/s. We have used the average eigenstrain approximation for a square lattice of infinitely long cylindrical holes aligned in the x₃-direction.

Moreover, we can compare the dispersion relations obtained from the EMT, Eqs. (1), (2), to the exact numerically computed dispersions (using the finite element method) for some subset of chosen PnC geometries. Figure 4 shows an example of this comparison for a PnC plate of thickness d = 10 nm and period a = 10 nm, for three different f = 0.1, 0.3, 0.5 in the Γ − X direction. We see that the EMT describes the lowest linear S and SH modes accurately over most of the Brillouin zone (BZ), whereas the flexural A-mode can be described by the EMT for a somewhat smaller range of k-vectors, with deviations starting around ω ~ 2 − 3 × 10¹¹ rad/s. This frequency starts to limit the accuracy of the EMT theory at temperatures around 1 K in this case, as will become clearer later. More dispersion relation comparisons can be found in Supplementary Figs., for the cases d = 10 nm, a = 10 nm (Supplementary Fig. 1), d = 50 nm, a = 50 nm (Supplementary Fig. 2), d = 100 nm, a = 50 nm (Supplementary Fig. 3), and d = 300 nm, a = 50 nm (Supplementary Fig. 4).

Fig. 4: Dispersion relations of the lowest modes.

The EMT results for A (yellow dot-dash), S (green dash), and SH (purple dot) mode, calculated based on Eqs. (1), (2), are compared with the exact numerically calculated dispersions (black solid), for three different filling factors f = 0.1, 0.3, 0.5, in Γ-X direction. SiN parameters E = 250 GPa, ρ = 3100 kg/m³, ν = 0.23 and a plate thickness d = 10 nm were used, and in the exact calculation a period a = 10 nm.

Low-temperature heat capacity of a metaplate

The phononic heat capacity of a solid is given by the temperature derivative of the internal energy of the solid

$$C=\frac{\partial U(T)}{\partial T}=\frac{\partial }{\partial T}\mathop{\large \sum}\limits_{{{{{{{{\bf{k}}}}}}}},j}\frac{\hslash {\omega }_{j}({{{{{{{\bf{k}}}}}}}})}{\exp (\hslash {\omega }_{j}/{k}_{{{{{{{{\rm{B}}}}}}}}}T)-1},$$

(4)

where the summation is over all phonon modes j over the first atomic Brillouin zone. The phonon dispersion relations ω_j(k), j = 1, …n are the only unknown quantities, and in the low-frequency limit they depend only on the elastic properties of the material and the geometry of the solid, as was discussed in the previous section. For bulk samples, the k-summation can be converted to an integral, which can be extended over the whole k-space in the low temperature limit, as the Bose-Einstein occupation factor goes to zero for high k-values anyways. Finally, in the low temperature limit only the lowest three linear bulk modes contribute, and by elementary integration we arrive at the well known cubic Debye-law²⁴ for the volume specific 3D heat capacity

$${C}_{{{{{{{{\rm{3D}}}}}}}},V}=\frac{2{\pi }^{2}}{15}\frac{{k}_{B}^{4}}{{\hslash }^{3}}\left(\frac{2}{{c}_{T}^{3}}+\frac{1}{{c}_{L}^{3}}\right){T}^{3},$$

(5)

where all of elasticity is hidden in the definitions of the bulk transverse c_T and longitudinal c_L speeds of sound. Note that in ref. 11 the numerical prefactor is in error for the 3D Debye formula.

For the case of thin plates, the quasi-continuous wave vector space exists only in the 2D plane of the plate, so that the k-space integral is two-dimensional, and the conversion from a sum to an integral gives a factor proportional to the area of the plate A instead of a volume. Another feature is of course that the dispersion relations ω_j(k) are much more complex even for isotropic plates^21,25, with an infinite number of branches (in the continuous limit). Performing the temperature differentiation in Eq. (4), we can write for the area-specific heat capacity C_A (units [J K⁻¹ m⁻²])

$${C}_{{{{{{{{\rm{2D}}}}}}}},A}=\frac{{\hslash }^{2}}{16{\pi }^{2}{k}_{{{{{{{{\rm{B}}}}}}}}}{T}^{2}}\mathop{\large\sum}\limits_{j}{\int}_{2D}{{{{{{{\rm{d}}}}}}}}{{{{{{{\bf{k}}}}}}}}\,\,{\omega }_{j}^{2}({{{{{{{\bf{k}}}}}}}})\,{\sinh }^{-2}\left(\frac{\hslash {\omega }_{j}({{{{{{{\bf{k}}}}}}}})}{2{k}_{{{{{{{{\rm{B}}}}}}}}}T}\right).$$

(6)

In the low-temperature limit only the three lowest branches will contribute, and by using the lowest order analytical approximations for their dispersion relations, Eq. (6) can be integrated. For the case of an isotropic plate it was done in ref. ²⁵, yielding a temperature dependence C ~ aT + bT², where the linear term is generated by the flexural A mode and the quadratic term by the S and SH modes.

We now calculate similarly the low-temperature area-specific heat capacity for our metasolid plate, using the asymptotic analytical dispersions of Eqs. (1), (2). First, consider the shear and symmetric modes of Eq. (2). Inserting the form of ω(k) into Eq. (6) and carrying out the integration in polar coordinates, we get for those two modes

$${C}_{{{{{{{{\rm{S}}}}}}}},{{{{{{{\rm{SH}}}}}}}}}= \,\frac{3\zeta (3){k}_{{{{{{{{\rm{B}}}}}}}}}^{3}{T}^{2}{\rho }^{* }}{{\hslash }^{2}{\pi }^{2}}\cdot \left\{2\pi \sqrt{\frac{{\alpha }^{2}+\beta }{\beta (\beta +\gamma )}}\pm \frac{4}{\alpha }K\!\left(\frac{\gamma }{{\alpha }^{2}}\right)\right.\\ \left.\mp \frac{4({\alpha }^{2}+\beta )}{\alpha \beta }\Pi\!\left(-\frac{\gamma }{\beta },\frac{\gamma }{{\alpha }^{2}}\right)\right\},$$

(7)

where

$$\alpha ={{{{{{{\mathcal{B}}}}}}}}-{c}_{66}^{* },\quad \beta =4{{{{{{{\mathcal{B}}}}}}}}{c}_{66}^{* },\quad \gamma =2{{{{{{{\mathcal{A}}}}}}}}\alpha -{{{{{{{{\mathcal{A}}}}}}}}}^{2},$$

(8)

ζ(x) is the Riemann zeta function, K(m) is the complete elliptic integral of the first kind (with m = k², where k is the elliptic modulus), Π(n, m) is the complete elliptic integral of the third kind, and the upper (lower) signs correspond to the S (SH) mode. Here we note that the dispersion relations in Eq. (2) are always real-valued, as ${{{{{{{\mathcal{B}}}}}}}}$ and ${c}_{66}^{* }$ are non-negative and real due to general physical constraints for the elasticity tensor that cause the matrix $({c}_{IJ}^{* })$ and all its principal submatrices to be positive definite (and thus have non-negative determinants). Interestingly, the S and the SH mode contributions produce opposite signs for the elliptic integral terms, causing them to cancel out in the summation for the total heat capacity.

When the antisymmetric mode given by Eq. (1) is plugged into Eq. (6) and integrated, only a single complete elliptic integral of the first kind remains. Combining the three mode contributions together, we finally end up with an analytical formula for the area-specific heat capacity of the phononic crystal metaplate in the low-temperature limit:

$${C}_{{{{{{{{\rm{EMT}}}}}}}},A} =\frac{{k}_{{{{{{{{\rm{B}}}}}}}}}^{2}}{3\hslash d}\sqrt{\frac{3{\rho }^{* }}{{{{{{{{\mathcal{B}}}}}}}}}}K\!\left(\frac{{{{{{{{\mathcal{A}}}}}}}}}{2{{{{{{{\mathcal{B}}}}}}}}}\right)T \\ \quad+\frac{3\zeta (3){k}_{{{{{{{{\rm{B}}}}}}}}}^{3}}{{\hslash }^{2}\pi }\left(\frac{{\rho }^{* }}{{{{{{{{\mathcal{B}}}}}}}}}+\frac{{\rho }^{* }}{{c}_{66}^{* }}\right)\sqrt{1-{{{{{{{\mathcal{A}}}}}}}}\frac{{{{{{{{\mathcal{A}}}}}}}}-2{{{{{{{\mathcal{B}}}}}}}}\,+\,2{c}_{66}^{* }}{({{{{{{{\mathcal{A}}}}}}}}-2{{{{{{{\mathcal{B}}}}}}}})({{{{{{{\mathcal{A}}}}}}}}+2{c}_{66}^{* })}}{T}^{2}.$$

(9)

The first T-linear part of the formula corresponds to the antisymmetric mode, while the quadratic part contains the contributions of the shear and symmetric modes. In the quadratic part we see two terms proportional to ${\rho }^{* }/{{{{{{{\mathcal{B}}}}}}}}$ and ${\rho }^{* }/{c}_{66}^{* }$ that resemble and converge to the full plate isotropic (inverse) wave speeds ${c}_{{{{{{{{\rm{S}}}}}}}}}^{-2}$ and ${c}_{{{{{{{{\rm{SH}}}}}}}}}^{-2}$, respectively, in the zero filling factor limit. It is then clear that in the f → 0 limit the above formula reduces to the previously reported analytical results for the full isotropic plates²⁵.

As a first look of the analysis of Eq. (9), in Fig. 5 we plot an example of the mode specific contributions for the heat capacity as a function of T (below 1 K) for a 50 nm thick SiN, which is a typical material from which suspended devices are made¹².

Fig. 5: Modal contributions (A, S, and SH) to the area-specific heat capacity of the holey SiN PnC metaplate.

Heat capacity is plotted vs. T, using Eq. (9). Three different filling factors f = 0, 0.25, 0.5 are shown, and isotropic SiN parameters E = 250 GPa, ρ = 3100 kg/m³, ν = 0.23 and a plate thickness d = 50 nm were used.

For this case, the flexural A mode dominates up to around 0.5 K. For all modes, interestingly, increasing f leads to an increase in the heat capacity, more strongly for the S and SH modes than for the A mode. The S mode is always seen to give about a factor of three smaller contribution than the SH mode. However, as can be directly seen from Eq. (9), the EMT theory predicts that the SH and S modes would contribute more with a denser material. A thicker plate would reduce the contribution of the A mode, but would also diminish the temperature range of validity of the lowest mode approximation. We turn to address the question of the range of accuracy of the EMT theory next.

Numerical analysis of the approximation error

The result of Eq. (9) includes several approximation steps for which it is difficult to give formal error estimates. It is clear, though, that in the f → 0 limit the average eigenstrain approximation reproduces the original isotropic elasticity tensor. Numerical comparisons, as well as experimental studies have been conducted in literature which indicate fairly good approximation for the average eigenstrain model over a wide range of the parameter f ∈ [0, 0.5]²². Another possible source of error is the number of elements chosen in the calculation of the infinite series S_I and S_IJ in Eq. (14). Nemat-Nasser et al.²² reported a typical difference of less than 1% between the truncation points n_i = ± 40 vs. n_i = ± 50. To effectively remove this source of error we have calculated all the results in this work using ∣n_i∣ ≤ 2000.

Apart from the homogenization error, the exact heat capacity of a plate at an arbitrary non-zero temperature includes contributions from higher frequency phonon branches (with a non-zero ω at k = 0), but in Eq. (9) we only considered the three acoustic modes (with ω → 0 for k → 0) that survive in the T → 0 limit. Eq. (9) is therefore exact only in the T → 0 K, f → 0 limit.

For the purpose of using Eq. (9) in practice, it is important to find some guidelines for which parameter ranges the error is not too high. The approximation error can be analyzed by comparing the heat capacity obtained from the effective medium theory with the exact numerically calculated heat capacity of the PnC. We use the finite element method to calculate phonon dispersion relations for several PnC plate designs. A second-order Gaussian quadrature is used for the numerical integration of Eq. (6). Depending on the lattice constant, either 1600 (a = 50 nm) or 7800 (a = 10 nm) quadrature points in the irreducible Brillouin zone were chosen, with an increasing sampling density at the k → 0 limit. A typical size of the finite element mesh was ~ 7000 elements. We remark that these seemingly exaggerated choices were actually important for the reduction of the error arising from numerical integration over a temperature range of many orders of magnitude.

Keeping isotropic silicon nitride as our reference material, in Fig. 6 we compare the relative error of the effective medium approximation to the exact numerically calculated heat capacity of a selection of PnC designs. Significant modifications to the phonon spectrum often appear near the first Brillouin zone (BZ) boundary (Bragg conditions), so we have chosen to compare designs with smallish lattice constants (a = 10 nm and a = 50 nm) to move the corresponding frequency ranges of the BZ edges higher. To ascertain that the excited thermal phonons do not yet satisfy the Bragg conditions in the temperature range considered, we also plot the dominant thermal phonon wavelength for the plate (Lamb) modes corresponding to the temperature in the upper x-axis scale. The plate thickness d was also varied over the range 10–300 nm.

Fig. 6: Relative error of the effective medium theory.

The plots show the ratio of the EMT heat capacity (C_EMT) with the numerically calculated exact specific heat capacity C_PnC, for SiN phononic crystal plates (PnC) of varying lattice constants a and thicknesses d [a a = 10 nm, d = 10 nm; b a = 50 nm, d = 50 nm; c a = 50 nm, d = 100 nm; d a = 50 nm, d = 300 nm]. The upper x-axis scales give the dominant thermal wavelengths of an uncut plate.

We see in Fig. 6 that for d = 10 nm, the effective medium theory (EMT) is very good over a wide range of temperatures: the error is roughly within a few percentage points for f ≤ 0.5 below T = 100 mK, and still only ~ 20% even at 1 K. (A typical low-temperature experiment has a base temperature of 10–100 mK). Similarly, for the other designs considered, our EMT is very good in the low temperature limit for the cases f ≤ 0.5, staying below ~ 20% for T < 0.1 K. The case f = 0.6 seems to deviate noticeably from the exact results even at the very low temperature regime, indicating that the average eigenstrain model fails to properly address the highest filling factor nanostructures. Based on the comparisons in Fig. 6, we can conclude that Eq. (9) gives a good estimate for the SiN PnC heat capacity for temperatures below 100 mK up to medium filling factors f ~ 0.5. Increased accuracy at higher temperatures would most likely require the inclusion of higher order plate modes, for which analytical formulas do not exist. Most importantly, though, being analytic, Eq. (9) can easily be used to predict how various design parameters affect the PnC plate heat capacity.

The effect of the geometric design and material parameters

Coming back to Eq.(9), we now discuss its predictions of the effect of several design parameters on the low-temperature heat capacity of PnC plates. First, it is worth mentioning that the plate thickness dependence only appears in the first, flexural wave term, which is linear in temperature. This results from the fact that the first order approximation is sufficient for the horizontal shear and symmetric modes. (Second order approximations for these modes would have small negative quadratic thickness dependencies which would have an impact only at higher temperatures²³.) Thus, the plate thickness has the biggest effect at the lowest temperatures where the linear term dominates (see Fig. 5). This observation also agrees with our previous exact calculations for SiN and Si¹¹, where the thickness parameter was seen to dominate at the lowest temperatures, while the lattice constant had a larger effect at ~ 10 mK and above. EMT theory thus confirms that the area-specific heat capacity can also be inversely proportional to the plate thickness for PnC plates, similar to what happens for unperforated plates^11,25.

The area-specific heat capacity also has simple linear (S and SH modes) and square root (A mode) dependencies on the effective density ρ* = (1 − f)ρ. Because of this linearity, the density of the material is a more sensitive design parameter for tuning the heat capacity at higher temperatures where the S and SH modes dominate.

The complete elliptic integral K(x) present in the T-linear term is real valued in the range $x\in \left[0,1\left[{{{{{{{\rm{with}}}}}}}}\,K(0)=\pi /2\right.\right.$, and it has a pole at x = 1. It seems therefore possible to strongly increase the heat capacity if a pattern/material combination with $x:= {{{{{{{\mathcal{A}}}}}}}}/(2{{{{{{{\mathcal{B}}}}}}}})\approx 1$ could be found. However, this is not straightforward, as we discuss in the following. The value of x can be analyzed as a function of f and Poisson ratio ν; Young’s modulus of the plate material does not affect this value. By increasing the hole filling fraction f, the value of x also increases slowly. It also increases when ν decreases. As a reference, for f = 0.5, ν = − 0.99 we have x = 0.734 and K(x) ≈ 2.1286. Increase in the value of K is therefore only ~ 36% higher compared to the isotropic unperforated plate case, so unfortunately the heat capacity increase by this route seems quite limited at best.

How the Young’s modulus, Poisson ratio and filling factor affect the heat capacity is a more complicated question and requires numerical comparison of the results. Taking again the d = 50 nm thick silicon nitride membrane as a reference, we see in Fig. 7 how the elastic parameters of the PnC matrix material affect the heat capacity, keeping ρ constant. Softer materials (E small) have a higher heat capacity, and over the range E = 50−400 GPa, the heat capacity can triple, and further improvement/reduction seems possible among natural materials. Increasing the filling factor does not change the heat capacity much between materials of different Young’s moduli: it mainly affects the shear modulus of the effective material, which only starts to play a role for the heat capacity when the SH mode becomes dominant in the softest materials at the higher end of the temperature range. Moreover, changing the Poisson ratio ν has only a small effect on the heat capacity. Over the range of ν for natural and artificial isotropic materials reported so far (including negative values of ν), only a ±5% overall difference is seen compared to the reference plate at the low temperature limit.

Fig. 7: Effect of the materials parameters on the EMT hear capacity.

The heat capacity of a 50 nm thick PnC EMT metasolid for various Young’s moduli (a) and Poisson ratios (b), compared to the EMT result of the same PnC structure with SiN parameters (E = 250 GPa, ρ = 3100 kg/m³, ν = 0.23).

In Fig. 8 we compare how changing the filling factor affects the EMT heat capacity in two different scenarios. On the left panel, the heat capacities of SiN PnC metaplates with different filling factors and various thicknesses are compared to full plates of the same thickness (represented by the color) by plotting the enhancement factor, i.e., the ratio C_EMT,A/C_full.

Fig. 8: Heat capacity enhancement factor.

The enhancement factor C_EMT,A/C_full is plotted vs. T for a changing filling factor f of SiN PnC metaplate with constant thickness (a), and constant area mass density [kg/m²] (b).

Material inside the holes would simply be removed in this case, as is done in practice. With an increasing filling factor, the density of the metaplate therefore decreases as more material is removed. As reported before¹¹, we see that removing material increases the heat capacity in this case. The relative increase does not seem to depend on the plate thickness at the lowest temperatures, as the plotted curves for various thicknesses overlap there. With the PnC patterning considered in this study, heat capacity can be roughly doubled at the lowest temperatures by choosing f = 0.5. Higher enhancements seem possible with increasing temperature, but the EMT theory starts to lose accuracy there.

If instead of removing the material we “mold” the plate by changing the plate thickness accordingly to maintain a constant area density, then an increasing filling factor decreases the heat capacity relative to a full plate with thickness d_ref, as seen in the right panel in Fig. 8. Similarly to the case where material was removed, the reference plate thickness d_ref does not seem to affect the heat capacity enhancement factor at the low temperature limit.

Thermal metamaterial

We also want to make comments on to what extent the thermal metamaterial considered here exhibits properties not found in naturally occurring materials. For low temperatures, the natural material to compare to is the 3D Debye solid. Considering its specific heat capacity/unit mass, it behaves “normally” in the sense that it only depends on material parameters and constants, and by adding more mass, heat capacity increases proportionally. It also has the property that the denser the material, the larger the heat capacity increase/unit mass. Now, the PnC metamaterial we have discussed has some “strange” and in some sense negative properties. If we consider the limit where the flexural mode dominates, then the total heat capacity (J/K) is proportional to $\sqrt{\rho (1-f)}A/d$, from which we can derive that the specific heat capacity/unit mass is proportional to $1/({d}^{2}\sqrt{\rho (1-f)})$. It still depends on the extensive variables d and f. By removing mass either by thinning the plate or by increasing the hole sizes (filling factor), one can increase the specific heat, and where one removes the atoms (plate surfaces or hole surfaces) affects the specific heat differently. Also, we see that the lighter the matrix material, the more the specific heat increases, entirely opposite to the natural case.

Conclusions

Focusing on phononic crystal geometries with periodic square lattice hole patterns in isotropic plates, we derived an analytical effective medium theory for the heat capacity of a phononic metaplate, applicable at low temperatures. The purpose of the theory is to greatly simplify the calculation of heat capacity of such structures, and to aid in understanding which parameters should be tuned to control it. We used a quasistatic elastic effective medium theory, allowing us to include the emergent anisotropic elastic response of the metasolid. We then proceeded to describe the characteristics of the three dominant wave modes at low temperatures, and using those results, derived a simple to use formula, which directly exposes how the heat capacity depends on various material and geometrical parameters. We analyzed the error by comparing the formula to the exact numerically calculated heat capacity data, and found rough guidelines for the hole pattern filling fraction and temperature where the theory is most applicable for silicon nitride. Finally, we analyzed how the design parameters affect the metaplate heat capacity, concluding that either increase or decrease can be achieved by tuning the parameters. The derived formula shows that the effects of different parameters (such as the density of the matrix material, its elastic modulus, the filling factor, the plate thickness) would add, providing simple blueprints to tune the heat capacity by over an order of magnitude. In the future, we can also consider the possibility of an equivalent effective medium theory for low temperature thermal conductance, but that theory is more complex because the anisotropic group velocity has to be taken into account, as well.

Finally, we comment on the possible applications. Such tuning of heat capacity can be relevant for low-temperature radiation detectors, as the heat capacity of the absorber element in a bolometric sensor dictates the temperature rise when electromagnetic quanta are absorbed in the element, thus affecting the sensitivity of the device. The most sensitive detectors operate at the temperature range of our study, at or below 0.1 K.

Methods

Exact mathematical description for the holey PnC plates

The exact mathematical description for the low-frequency phonon modes in the geometry of a plate of thickness d with a circular hole pattern in a 2D square lattice can be written as

$$\nabla \cdot {{{{{{{\bf{c}}}}}}}}:\varepsilon = -\rho {\omega }^{2}{{{{{{{\bf{u}}}}}}}}\quad {{{{{{{\rm{in}}}}}}}}D\backslash {{\Omega }}\\ {\varepsilon }_{ij} = \frac{1}{2}\left({\partial }_{i}{u}_{j}+{\partial }_{j}{u}_{i}\right)\\ {{{{{{{\bf{u}}}}}}}}{| }_{{B}_{i}^{\pm }} = \exp (i({{{{{{{\bf{k}}}}}}}}\cdot {{{{{{{\bf{x}}}}}}}}))\,{{{{{{{\bf{u}}}}}}}}{| }_{{B}_{i}^{\mp }},\qquad i=1,2\\ \hat{{{{{{{{\bf{n}}}}}}}}}\cdot {{{{{{{\boldsymbol{\sigma }}}}}}}}({{{{{{{\bf{u}}}}}}}}) = 0\quad {{{{{{{\rm{on}}}}}}}}(\partial D\backslash {{\Omega }})\backslash {B}_{i}^{\pm },$$

(10)

where ε is the infinitesimal strain tensor related to the displacement field u, the elastic properties of the matrix material are given by the mass density ρ and the elasticity tensor c = c_ijkl, which is assumed to be isotropic. Volume $D=[-\frac{a}{2},\,\frac{a}{2}]\times [-\frac{a}{2},\,\frac{a}{2}]\times [0,\,d]\subset {{\mathbb{R}}}^{3}$ is the simple tetragonal unit cell of the periodic PnC plate that has side lengths a (in the x₁–x₂ plane) and height d, and where ${B}_{i}^{\pm }$ are the periodic boundaries with the Bloch-Floquet boundary conditions. The hole is centered within the cell D and is a cylinder Ω standing in the x₃-direction, with a radius r < a/2 and a height d. The stress-free boundary condition (the last equation in Eq. (10)) is applied on the remaining boundaries. In the above, the double dot symbol is the tensorial double contraction operator. As usual, we express the elasticity tensor as a 6 × 6 matrix following the Voigt convention, so that the contraction operation reduces to weighted matrix multiplication (see e.g.,^17,21).

Detailed methodology for the calculation of the effective elasticity tensor

The effective elasticity tensor in the average eigenstrain approximation is given by¹⁷

$${{{{{{{{\bf{c}}}}}}}}}^{* }={{{{{{{\bf{c}}}}}}}}:\left\{{{{{{{{{\bf{1}}}}}}}}}^{(4{{{{{{{\rm{s}}}}}}}})}-f{({{{{{{{{\bf{1}}}}}}}}}^{(4{{{{{{{\rm{s}}}}}}}})}-{{{{{{{{\bf{S}}}}}}}}}^{{{{{{{{\rm{P}}}}}}}}})}^{-1}\right\},$$

(11)

where ${{{{{{{{\bf{1}}}}}}}}}_{ijkl}^{(4{{{{{{{\rm{s}}}}}}}})}=\frac{1}{2}({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})$ is the fourth order symmetric identity tensor and f is the volume filling factor of the void Ω. The tensor S^P = S^P(f, ν) only depends on the hole shape, the lattice geometry, and the Poisson ratio of the original matrix material.

For infinitely long cylindrical holes in a square lattice surrounded by an isotopic matrix material, the tensor S^P is given by¹⁷

$${S}_{ijkl}^{{{{{{{{\rm{P}}}}}}}}}= \,\frac{1}{2}\left({\delta }_{il}{S}_{I(j,k)}+{\delta }_{ik}{S}_{I(j,l)}+{\delta }_{jl}{S}_{I(i,k)}+{\delta }_{jk}{S}_{I(i,l)}\right)\\ -\frac{1}{1-\nu }{S}_{I(i,j)I(k,l)}+\frac{\nu }{1-\nu }{\delta }_{kl}{S}_{I(i,j)}$$

(12)

where S_I and S_IJ, I, J = 1, …, 6 are infinite sums that have to be numerically evaluated, and where indices I(i, j) are given by the elements of the matrix

$$I=\left(\begin{array}{lll}1&6&5\\ 6&2&4\\ 5&4&3\end{array}\right).$$

(13)

The only non-zero terms S_I and S_IJ in this geometry are

$$\begin{array}{ll}{S}_{1}={S}_{2}=f\cdot {\mathop{\sum }\limits_{{n}_{i} = 0}^{\pm \infty }}^{{\prime} }{\left(\frac{2{J}_{1}(x)}{x}\right)}^{2}{({\bar{\xi }}_{1})}^{2}\hfill\\ {S}_{11}={S}_{22}=f\cdot {\mathop{\sum }\limits_{{n}_{i} = 0}^{\pm \infty }}^{{\prime} }{\left(\frac{2{J}_{1}(x)}{x}\right)}^{2}{({\bar{\xi }}_{1})}^{4}\hfill\\ {S}_{66}={S}_{12}={S}_{21}=f\cdot {\mathop{\sum }\limits_{{n}_{i} = 0}^{\pm \infty }}^{{\prime} }{\left(\frac{2{J}_{1}(x)}{x}\right)}^{2}{({\bar{\xi }}_{1}{\bar{\xi }}_{2})}^{2},\hfill\end{array}$$

(14)

where $x=\sqrt{4\pi f}\sqrt{{n}_{1}^{2}+{n}_{2}^{2}}$, ${\bar{\xi }}_{i}={n}_{i}4\pi f/x$, ${n}_{i}\in {\mathbb{Z}}$ for i = 1, 2, and J₁(x) is the first order Bessel function of the first kind. The primed sum symbol here means that the case n₁ = n₂ = 0 is excluded from the summation. Elements ${c}_{ijkl}^{* }$ can also be given closed form expressions in terms of the filling factor, however, they remain complicated and still contain numerically approximated parts related to the above infinite sums²⁶, and we therefore prefer to compute the sums directly up to a high truncation order.