# Thermal properties of nanoporous materials, large scale modelling with the use of Monte Carlo phonon transport autocorrelation

David Lacroix, M. I. Nkenfack, Gilles Pernot, Mykola Isaiev

► To cite this version:

David Lacroix, M. I. Nkenfack, Gilles Pernot, Mykola Isaiev. Thermal properties of nanoporous materials, large scale modelling with the use of Monte Carlo phonon transport autocorrelation. Journal of Applied Physics, 2023, 134 (2), pp.025101. 10.1063/5.0155582. hal-04282245

HAL Id: hal-04282245

https://hal.science/hal-04282245v1

Submitted on 13 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

# Thermal properties of nanoporous materials, large scale modelling with the use of Monte Carlo phonon transport autocorrelation

D. Lacroix*¹, M.I. Nkenfack¹, G. Pernot¹, and M. Isaiev¹

¹Université de Lorraine, CNRS, LEMTA, Nancy F-54000, France

June 6, 2023

## Abstract

In the present work, we demonstrate the ability of a technique based on Monte Carlo resolution of the Boltzmann Transport Equation associated to the Green-Kubo autocorrelation of phonon heat flux to predict, at thermal equilibrium, the thermal conductivity tensor of nanoporous structures. This methodology, which is derived from a former work [D. Lacroix, M. Isaiev, G. Pernot, Phys. Rev. B, 104, 165202 (2021)] developed in the case of bulk systems, is used to predict thermal transport properties of Si porous matrices and Si phononic membranes at room temperature. A broad range of porosities and different pore network organisations are considered. Our results are compared to available experimental data and former modelling techniques. In addition, analytical models based on phonon mean free path are detailed and compared to numerical simulations.

## 1 Introduction

The study of the thermal properties of nano-structured materials and more specifically the control of these properties through small-scale nano-systems is of prime interest in various communities, for both the physical questions that are raised by this research and for the wide engineering applications that they promise. Over the last decade, significant efforts have been made to develop new semiconductor devices at sub-micron scales, such as nanoporous membranes [1] or phononic crystals [2, 3, 4], or suprlattices with embedded nanoparticles [5] where those particular designs strongly impact the thermal conductivity and even to control the anisotropy of the thermal transport in periodic structures.

For systems where the porosity is geometrically “organized”, commonly called Phononic Crystals (PNC), many studies have focused on the possibility of achieving interference effects in the transport of phonons, similar to what can be observed in optics with photons. In practice, coherence effects are not observable at temperatures close to ambient. It would only be possible at ultra-low temperatures, where the frequency spectrum of phonons is reduced and their wavelength is large compared to the size of the structures. However, if coherence effects are difficult to reach, nano-structuring allows to control the phonon mean free path and thus tailor the thermal conductivity of the fabricated systems. These considerations have led to different experimental studies associating elaboration and characterization of PNC, in several research groups, for different geometries (pore size and shape, film thickness, etc.), and it remains a very active research field [6, 7, 8, 9, 10, 11, 12].

On the numerical side, different models have been developed to predict the thermal transport properties of these objects. These include analytical approaches based on the evaluation of the mean free path (mfp) of energy carriers as a function of the geometry considered. One can also find works based on “finite element” tools that solve the classical equations of elasticity and give access to the band structure of these objects; these calculations are however limited to low temperatures [13]. Other approaches are nevertheless possible to evaluate the transport properties at higher temperatures. This is the case of molecular dynamics (MD) which allows a more realistic approach of atomic vibrations and induced thermal transport. However, if this tool is well mastered and

<footnote>*Corresponding author, david.lacroix@univ-lorraine.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Studied porous networks: left) packed spherical pores; center) aligned PNC membrane; right) staggered PNC membrane.</figcaption></figure>

classically used to describe the thermal conductivity of nano-structured semiconductors (nano-films, nano-wires, etc.), it remains limited to small structures (a few nanometers to a few tens of nanometers). In these studies, the compulsory requirement of a “scale change” is often an underlying assumption to explain difference that can be observed with experimental behavior at larger scales. Hence why, this assumption can be challenged in complex systems. In addition, boundary conditions which are often used, such as the “periodicity” of the system studied, may be far from the experimental reality or simply inappropriate to consider phonons with wavelengths greater than typical simulation domain size (tens of nanometers). To overcome these drawbacks, other simulation methods must be considered.

In the present paper, we discuss a general methodology based on the use of the Monte Carlo (MC) solution of the Boltzmann Transport Equation (BTE) combined to the equilibrium Green-Kubo formalism in complex porous networks like those plotted in Fig. 1. This work is in the continuity of a previous study made on bulk materials [14]. We demonstrate that a similar methodology can be extended to evaluate the thermal conductivity tensor of several micron length devices with included nanoscale features.

In the following section, the simulation methodology will be briefly recalled. Then, the geometrical features of the simulated devices will be presented. In the latter section, parametric models based on the energy carriers “mean free path” and Matthiessen’s rule will be described in order to assess general trends related to thermal conductivity variations with respect to the device porosities. In result section, different configurations of organized porous systems are investigated, at room temperature and their thermal conductivity tensors compared to the literature available experimental data. Developments in progress and perspectives to this work end the study.

# 2 Monte Carlo - Green-Kubo formulation of Boltzmann Transport Equation for phonons (BTE)

## 2.1 General formalism

At mesoscales, phonon dynamics is governed by the BTE. In the classical assumption where no external forces are applied on the considered system, it reads:

$$ \frac{\partial f}{\partial t} + \nabla_{\mathbf{K},p}\omega \cdot \nabla_{\mathbf{r}}f = \frac{\partial f}{\partial t} \Big|_{\text{scat}} \qquad (1) $$

with $f(\mathbf{K}, p, \mathbf{r}, t)$ the energy carrier distribution function. left hand side (lhs) of Eq. 1 refers to phonon transport or drift, while right hand side (rhs) is related to phonon scattering. In this
framework, the BTE can be solved with stochastic techniques known as “MC formalism” where phonons are considered as particles, with a given frequency and polarization. Such MC modeling relies on the following stages.

* First, their displacement (drift) during a given time step $\delta t$ is calculated according to their group velocity, which is frequency and polarisation dependant, along a randomly selected propagation direction.

* Then, scattering process involving 3-phonons interactions and phonon-impurity / defect are considered. The latter are known as intrinsic mechanisms which depends on the material nature. These scattering mechanisms tend to restore a thermodynamic equilibrium in the material in response to the drift of phonons. “Normal”, “Umklapp” and “Impurity” scattering processes are here described in the frame of the relaxation time approximation which allows to linearize the rhs of Eq. 1. In the MC procedure, the intrinsic phonon scattering is addressed through the calculation of a collision probability $P_{scat}$ (Eq. 2), the latter being a function of calculated phonon lifetime and simulation time step $\delta t$ with $\tau_{tot}(\omega, p, T)$ the total relaxation time of phonon mode (K, p). This relaxation time can be expressed according to the Matthiessen rule about individual relaxation times inverse summation.

$$P_{scat}(\omega, p, T) = 1 - \exp \left[ \frac{-\delta t}{\tau_{tot}(\omega, p, T)} \right] \quad (2)$$

* In addition to those intrinsic scattering mechanisms, boundary scattering on the edges of simulation domain can be considered in the case of nanostrucures. The latter one is different from the “intrinsic” ones as it does not alter the phonon spectra but modify the propagation direction and mostly reduce the phonon mean free path. Such boundary scattering is a major mechanism in porous systems as it dictates heat transport at the considered length-scales, i.e. lowering the mfp lowers thermal conductivity. Boundary scattering is a part of the drift process and is treated as a particular feature of phonon displacement. When an energy carrier “hits” a physical boundary of the simulation domain (pore or wall), the phonon is reflected specularly or diffusely according to the nature of the interface.

More details regarding the MC simulation procedure algorithm are given in [14]. In addition, phonon dispersion properties and phonon intrinsic scattering lifetime are given in appendices for silicon which is the test material considered in this work.

At this stage, it should be recalled that these MC simulations are done for a prescribed temperature, i.e. no thermal gradient is imposed in the studied systems as it is usually the case in MC modelling of phonon transport [15, 16, 17, 18]. As a consequence, thermal conductivity of the studied systems is evaluated by another approach than the one based on the assumption of the Fourier’s law applicability. In this study, we proceed as in our previous work [14], with the calculation of heat flux autocorrelation with Green-Kubo formulation to determine the thermal conductivity tensor of the system. The latter reads:

$$k_{\alpha,\beta} = \frac{V}{k_B T^2} \int_0^\infty \langle q_\alpha(0)q_\beta(t) \rangle dt \quad (3)$$

where $\alpha, \beta$ are the Cartesian coordinates of $\mathbf{q}$ ($q_x, q_y$ and $q_z$). For each time step, the heat flux carried by each “particle” is calculated and stored. The resulting heat flux history is then split into $M = N_t/10$ samples to compute discrete values of the heat flux autocorrelation (HFAC). In this work, we used the same formalism as the one detailed by P.K. Schelling et al. [19]

$$k_{x,y} = \frac{V}{k_B T^2} \frac{\delta t}{N_p} \sum_{p=1}^{N_{ph}} \sum_{m=1}^{M} \frac{1}{N_t - m} \sum_{n=1}^{N_t - m} q_x(m + n)q_y(n) \quad (4)$$

The computed thermal conductivities are averaged on the number of launched particles ($N_{ph}$) in the MC process. Time steps are typically of one picosecond in these porous systems where boundary scatterings dominate even at low temperatures (no ballistic transport and thus short correlation time). More details about convergence and accuracy of the method can be found in [14].
## 2.2 Analytic models for the studied systems

Thermal transport properties like thermal conductivity $k$ in several bulk semiconductors can be modelled with the “kinetic theory model” [20, 21] with a reasonably good accuracy. This is in particular true in the systems where phonon transport is diffusive. With this formalism, the bulk thermal conductivity $k_{bulk}$ integrated on the whole phonon spectrum reads:

$$k_b = \frac{1}{3} C_{ph} v_g \Lambda_b \qquad (5)$$

with $C_{ph}$ the phonon heat capacity, $v_g$ the energy carriers average group velocity and $\Lambda_b$ the phonon mean free path in bulk. In the latter expression, for a given material, it can be assumed that the velocity and the heat capacity remain constant. In the specific case of porous devices, assuming that pores induce no contribution to heat propagation (ideal case of “empty pores”), the phonon mean free path becomes $\Lambda_{eff}$ depending on the system porosity $\phi$, the bulk $\Lambda_b$ and the porous network $\Lambda_{pn}$ mfps. If independent phonon scattering mechanisms are assumed, the Matthiessen’s law can be applied and $\Lambda_{eff}^{-1} = \Lambda_b^{-1} + \Lambda_{pn}^{-1}$. The effective porous thermal conductivity $k_{eff}$ in the general framework of kinetic theory is:

$$k_{eff} = \frac{1}{3} C_{ph} v_g \Lambda_{eff} \qquad (6)$$

According to the above assumption, the ratio of Eq. 6 and Eq. 5 leads to the general dimensionless thermal conductivity $\kappa$:

$$\kappa = \frac{k_{eff}}{k_b} = \frac{1}{1 + \Lambda_b / \Lambda_{pn}} \qquad (7)$$

Considering the above points, the dimensionless thermal conductivity $\kappa$ of a porous system is characterized by the porosity of the device and the geometric features of the pore network; i.e. pore radius and spacing. As suggested in [16] and [22], a connection between the pore mfp and the porosity can be directly addressed as a simplified configuration of the average mfp outside and inside an inclusion: $\Lambda_{pn} = (1 - \phi)\Lambda_p^{out} + \phi\Lambda_p^{in}$. In the case of empty inclusions (pores), the latter becomes: $\Lambda_{pn} = (1 - \phi)\Lambda_p$. Eventually, Eq. 7 is rewritten under the general form:

$$\kappa(\phi) = \frac{1}{1 + \frac{\Lambda_b}{(1-\phi)\Lambda_p^{out}}} \qquad (8)$$

with $\Lambda_p^{out}$ the phonon mean free path in the solid of the porous system. In the limit case of a null porosity $\Lambda_p^{out} = \Lambda_b$ and $\kappa$ tends to 1. At this stage, the derived formalism remains general and does not depend on the shape and organization of pores in the device. However, having some dedicated analytic models that can help to evaluate thermal transport properties in porous systems from pristine cases (i.e. bulk or membrane without any pore) using average carriers mean free path and geometric considerations (i.e. pore shape and their organisation) is of interest for many applications. In the following subsections $\Lambda_p^{out}$ will be defined for the three considered cases described in Figs. 1 and 2 having in mind such considerations. That being said, those models are not proposed to give a validation of the MC-GK calculations detailed in the result section, but rather to try to give some insights about pore boundary scattering mechanisms at play in tailored nanoporous devices.

### 2.2.1 Organized spherical pores

The first case considered here is a theoretical one and corresponds to a stack of spherical inclusions set along the three main Cartesian directions, see Fig. 1-(left). The elementary cell that can be periodically reproduced by translation along $x$, $y$ and $z$ axes is plotted in Fig. 2-(left).

These elementary cells are used to calculate the phonon mfp around a pore $\Lambda_p^{out}$. To do so one can consider the invariant property of diffusive random walks that states that “the average length $\langle\Lambda\rangle$ of the random walk trajectories from entry point to first exit point is independent of the characteristics of the diffusion process and therefore depends only on the geometry of the system” [23]. This property is known as “Cauchy formula” (chord length distribution), it reads:

$$\langle\Lambda\rangle = \frac{4V}{S} \qquad (9)$$

where $V$ and $S$ are the volume and the collision surface of the considered domain that does not include any inclusion. This relation can be straightforwardly applied in the case of convex bodies.
<figure><img src="image_3.png" /><figcaption>Figure 2: Elementary cell for nanoporous mean free path calculation; left) Organized spherical pores, center) Aligned PNC membrane, right) Staggered PNC membrane.</figcaption></figure>

Here, for the specific cases discussed above (non-convex bodies), the unit cell includes a pore (convex hole) with given volume (sphere or cylinder) that needs to be removed from the unit cell (convex body). For such geometries, generalization of the Eq. 9 can be derived [24] in frame of Euler functions $\Gamma$. This derivation is complex, thus a simple derivation of mfp can be assessed through simple geometric considerations. Here, a particle travelling within the cell can diffusely scattered on pore edge surface, this surface can be limited to half of the pore volume as forward scattering through the pore is not possible. With this assumption, the modified Cauchy formula is supposed to be:

$$ \langle \Lambda \rangle = \frac{4(V_{cell} - 1/2V_{pore})}{S_{pore}} \qquad (10) $$

In the latter expression, cell volume and pore volume are calculated considering periodic unit cells of the Fig. 2. On the other hand, scattering surface remain limited to the pore surface (convex hole) as no scattering is expected from the periodic boundary edges of the unit cell (convex body). Thus, one can compute, in the present case, the mfp around the pore in the elementary cell shown in Fig. 2-(left) as:

$$ \Lambda_{p-sph}^{out} = \frac{[a^3 - (4\pi R^3)/6]}{\pi R^2} \qquad (11) $$

where $a$ is the periodic side length of the cell. Using the same formalism, one recovers straightforwardly the inner mfp within an inclusion as $\Lambda_p^{in} = 4R/3$. Eventually, the reduced thermal conductivity for a porous device made of aligned spheres is:

$$ \kappa_{sph}(\phi) = \left\{ 1 + \frac{\Lambda_b \pi R^2}{(1 - \phi)[a^3 - (4\pi R^3)/6]} \right\}^{-1} \qquad (12) $$

In the following, the same formalism is used to study PNC membranes.

### 2.2.2 PNC membranes

In the case of aligned and staggered PNC membranes, the procedure is applied with respect to the elementary cells plotted in Fig. 2-(center & right). For both cases, the collision's surface is the inner cylinder ($S = 2\pi Rh$), while the elementary cell volumes $V$ are obtained using a square or an equilateral triangle base. Applying the above considerations, one can obtained the mfp around holes as (with subscripts $A$ and $S$ for aligned and staggered PNC configurations respectively):

$$ \Lambda_{p-A-PNC}^{out} = \frac{[2a^2 - \pi R^2]}{\pi R} \\ \Lambda_{p-S-PNC}^{out} = \frac{[\sqrt{3}a^2 - \pi R^2]}{\pi R} \qquad (13) $$

On this basis the dimensionless thermal conductivities in both geometries are given below. However, it shall be noted that bulk mfp is no longer the relevant parameter in the case of a PNC membrane that has a small thickness, hence $\Lambda_f$ is considered for film.

$$ \kappa_{p-A-PNC}^{out} = \left\{ 1 + \frac{\Lambda_f \pi R}{(1 - \phi)[2a^2 - \pi R^2]} \right\}^{-1} \\ \kappa_{p-S-PNC}^{out} = \left\{ 1 + \frac{\Lambda_f \pi R}{(1 - \phi)[\sqrt{3}a^2 - \pi R^2]} \right\}^{-1} \qquad (14) $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Thermal conductivity of aligned spherical pores versus material porosity. Comparison to analytic modelling based on phonon mfp evaluation.</figcaption></figure>

These models (Eq. 12 and Eq. 14) will be used, in the next section, to evaluate thermal transport properties obtained using the stochastic model based on the combined use of Monte Carlo modelling of phonon transport and Green-Kubo phonon heat flux autocorrelation. Besides, in appendice section we provide complementary information about the reliability of the pore boundary scattering mfp models depicted by Eq. 11 and Eq. 13 through direct evaluation of the latter by using Monte Carlo ray-tracing techniques as it was done in former studies [25, 16].

## 3 Results

Results section is divided in three parts. The first one is related to porous systems with spherical pores that are organized along a 3D-Cartesian grid (this is our theoretical model). The second and the third parts deal with PNC membranes where the pores are hollow cylinders that can be aligned (section 3.2) or staggered (section 3.3). In both cases, our numerical simulations aim at reproducing experimental characterization of the thermal conductivity in Si PNC membranes elaborated in the group of M. Nomura and which were investigated in a former work using “classical” MC method [26, 27].

### 3.1 Organized spherical pores

Thermal conductivity of devices with spherical pores is given on Fig. 3. In this work, both number of pores $N_p$ and their radius are changed. The simulation domain is a cube of 1 µm length with periodic reflections on its faces (“bulk”-like system). For such boundary condition, a phonon that leaves the domain through a given plane defined by external normal $n_x$, $n_y$ or $n_z$ reenters the domain through the plane defined by opposite normal vector. Such boundary condition allows us to preserve the heat flux propagation direction. On the other hand, reflection on the edges of the pore is purely diffuse and back-scattering propagation direction, after a collision event, is randomly sampled from the collision point. Three cases are considered in order to investigate porosity and mean free path impact in such systems with an isotropic pore’s distribution. In this study, $N_p$ is set to: 5 × 5 × 5, 8 × 8 × 8 and 10 × 10 × 10 pores distributed regularly along x, y and z directions. The porosity varies in the range 1% < $\phi$ < 50%, meaning that the pore radius is adjusted accordingly with respect to $N_p$ (from 6 nm to 98 nm according to $N_p$ and $\phi$ values).

First, the model given by Eq. 12 is used to recover bulk mean free path $\Lambda_b$ as it is the only unknown parameter in this equation. A regression analysis of the MC-GK thermal conductivity data, with $k_b$=127.8 W m⁻¹ K⁻¹ (obtained with MC-GK simulation on bulk case) leads to an average bulk mean free path of $\Lambda_b$=376 nm. The latter is coherent with usually observed value of 300 nm at room temperature. The “overestimation” compared to the literature of $\Lambda_b$ can be reduced while considering porosities larger than 10%. As a matter of fact, when the pore radius
<figure><img src="image_5.png" /><figcaption>Figure 4: Thermal conductivity versus phonon scattering surface per unit volume ($S_{scat}/V$) in porous silicon with aligned spherical pores. Power law fit for large pore diameters (red line).</figcaption></figure>

becomes small, $\Lambda_{p-sph}^{out}$ becomes large (see Eq. 11) and does not verify the assumption of “closed domain” where phonons are confined. Yet, in this formalism, the dimensionless $\kappa_{sph}(\phi = 0)$ is set to 1 and thus ensures a correct fitting of MC-GK data. When $\Lambda_b$ is evaluated with porosities $\phi > 10\%$, we obtain $\Lambda_b=350$ nm; which is close to commonly admitted value for silicon.

As expected, it can be seen, that increasing the number of spheres decreases the thermal conductivity for each tested porosities. This is logical and can be simply explained as reducing the pore size for a given porosity increases the scattering surface per unit volume ratio $S_{scat}/V$, leading to enhance the phonon scattering and thus to a thermal conductivity reduction. Here, the material volume $V$ is the volume of the box minus the volume of the pore. In the case of organized spherical pores, this ratio can be expressed considering the unit cell defined in Fig. 2-(left).

$$ \frac{S_{scat}}{V} = \frac{3\phi}{R(1 - \phi)} \qquad (15) $$

Porous samples thermal conductivity can be plotted in this frame. Such representation can help to understand when the pore scattering becomes dominant in the thermal transport process. In Fig. 4, for all the considered MC-GK simulations where the porosity is larger than 10%, thermal conductivity decrease, according to a power-law with the surface per unit volume ratio. Here, when $\phi > 10\%$ $k \propto (S_{scat}/V)^\alpha$ with $\alpha = -0.6$. Below this “threshold” porosity, thermal conductivity is weakly modified by pore occurrence. Such observation is consistent with the one done for the bulk mfp evaluation validity range. Such behaviour was previously reported for MD simulations in porous silicon where the length scales of the elementary cell where much smaller [28] (side length $a = 10 \times a_0 = 5.43$ nm and $3.3 < \phi < 52.5\%$). In the latter calculation, power-law dependence was also noted with similar $\alpha$ exponent close to -0.5. More recently, in MD based study on porous silicon, similar trends were observed and expressed in the framework of Santaló’s formula [29].

In addition to the classic porosity or scattering surface to volume dependence of thermal conductivity, one can wonder about the phonon “mean free path” (or interaction length) distribution in porous systems. Here, we take advantage of phonon tracking all along the simulation duration to evaluate such mfp distribution. In all our simulations, a set of selected particles (phonons) are followed from their initial location up to their final location and their “jumping length” properties are calculated each time a phonon undergoes a scattering event changing its propagation direction (resistive event). The latter can be scattering with a pore, umklapp scattering or lattice defect scattering. From the thermal transport viewpoint, it corresponds to interactions which decrease energy packet correlation.

In the inset of Fig. 5 is plotted a typical mfp distribution histogram in silicon at room temperature for $\phi = 5\%$. It shows a continuously decreasing distribution with mfp that ranges between tens and hundreds of nanometers. This kind of distribution can be adjusted according to a “Generalized Pareto model” with two parameters: $k$ which is the shape parameter and $\sigma$ the scale parameter.
<figure><img src="image_6.png" /><figcaption>Figure 5: Mean free path probability plot according to generalized Pareto model for $\phi = 5\%$, $N_p=1000$ and $R=22.85$ nm; the inset shows mfp distribution.</figcaption></figure>

Increasing $k$ enlarge the tail of the distribution (i.e. more events characterized by long mfp), while increasing $\sigma$ shift the maximum of distribution toward larger mfp. In the present case we use a probability plot to show the good agreement between the model and the MC-GK extracted data. In this case, less than 1% of the computed mfp deviate from the Pareto distribution due to the finite size of the domain; i.e. simulation box of $1 \times 1 \times 1 \mu m^3$. With the generalized Pareto formalism, one can extract mean $\overline{mfp}$ and median $\mathcal{M}(mfp)$ mean free paths for the considered distribution knowing $k$ and $\sigma$.

$$ \overline{mfp} = \theta + \frac{\sigma}{1-k} \qquad (16) $$

$$ \mathcal{M}(mfp) = \theta + \frac{\sigma(2^k - 1)}{k} \qquad (17) $$

In previous expression, $\theta$ parameter is the location parameter always equal to 0 in our adjustments. With this formalism the $\overline{mfp}$ is found to be close to 114 nm in this system made of 1000 spherical pores with a radius of 22.85 nm ($\phi = 5\%$). Depending on the porosity of the system the $\overline{mfp}$ range from 130 nm to less than 70 nm when porosity reach 50% in this silicon cubic domain. The above discussion describes the fitting procedure applied to all the simulation data set, i.e. for $N_p=125$, 512 and 1000 pores, considering 10 distinct porosities in the range $0 < \phi < 50\%$. This allows to plot the evolution of porous sample thermal conductivity's versus $\overline{mfp}$ (see Fig. 6). From the log-scale plot given in Fig. 6, one can see that thermal conductivity scales with the mfp derived from Pareto fitting of phonon "jumping length" between resistive scattering events giving a general trend $k \propto (\overline{mfp})^\beta$. Here, $\beta = 2.84$, close to 3 which can make sense in this ideal case where the symmetry is respected along the three Cartesian coordinates.

## 3.2 Aligned phononic membranes

In the case of PNC membranes, we consider silicon thin membranes with the following characteristic dimensions (length $L_z=25 \mu m$, width $L_y=5 \mu m$ and height $L_x=145$ nm. Here, a schematic view of such PNC membrane corresponds to the "center case" given in Fig. 1. For all simulations and in the corresponding experiments, the membranes characteristic lengths ($L_x$, $L_y$, and $L_z$) remain constant. The only parameters that change are: the pitch of the membrane $a$ (which is the same in $y$ and $z$ directions) and the radius of the holes. Five distinct pitches were considered in aligned PNC membranes ([160, 200, 300, 350, 500] nm). the hole radius vary in the range of 15 nm to more than 200 nm, depending on the considered pitch. It allows to tune the porosity of the membranes between 0% (plain membrane) and roughly 60%.
<figure><img src="image_7.png" /><figcaption>Figure 6: Thermal conductivity as a function of the average phonon mean free path in porous system made of spherical organized pores. The red line corresponds to a power law fit.</figcaption></figure>

In those devices, the MC-GK simulation method was applied to evaluate the thermal conductivity tensor as a function of the membrane's patterns as well as the characteristic phonon mean free path. A typical example of phonon energy carrier tracking with such calculation's procedure is given in Appendices in Figs. 14 and 15. In the latter, five distinct phonon tracks are plotted. In *x* and *y* directions phonon packet reflections are diffuse while in *z* direction, we assume “periodic” boundary conditions which allow a free flow of phonons. The *z* direction in our calculation is the heat flow propagation direction considered in micro-TDTR experiments achieved by R. Anufriev [26, 27].

First, the thermal conductivity tensor of a typical PNC membrane with a pitch *a*=300 nm is given in Fig. 7. Autocorrelation of the heat flux is computed along the three main directions. As expected, in *x* direction TC value $k_x$ is much smaller than bulk value appraised with Si Callaway-Debye lifetimes ($k_b$=127.8 W m⁻¹ K⁻¹). Here $k_x$=20 W m⁻¹ K⁻¹ for all considered porosities. This makes sense as hollow cylinder axis is along *x* direction and the thickness of the membrane ($L_x$=145 nm) is smaller than average phonon mean free path at room temperature. Here, the $k_x$ component is only limited by the film thickness. Along *z* direction, the energy carriers are free to flow (periodic BC) but their paths are limited by cylindrical pore occurrence. As the pore radius gets larger, the “neck” between pores is reduced as well as the phonon flow. The thermal conductivity value in this direction $k_z$, computed with the MC-GK method, is in good agreement with micro-TDTR measurements carried out on membranes with the same geometric features, giving a additional confirmation of the methodology reliability. Eventually, the TC in *y* direction remains lower than $k_z$ as the system width is smaller and finite ($L_y$=5 µm). For the largest pore diameter (here *R*=135 nm), the neck is equal to 30 nm and phonon confinement is even more pronounced than in the film thickness direction leading to $k_y$ values below 20 W m⁻¹ K⁻¹. Similar calculations were performed for the four other pore spacing's, showing similar trends.

On this general ground, thermal conductivity of PNC aligned membranes can be calculated as a function of porosity for the five distinct pitch values and, as it was done for spherical pores, compared to the above established models (see Eq. 14. Here *a* values and pore radii *R* are prescribed and set to correspond to the experimental values reported in [27], except for the smallest pore radius (*R* < 40 nm) that were only numerically considered to assess thermal continuous modelling of thermal properties from pristine membranes to highly porous ones. All those calculations (66 different cases) are plotted in Fig. 8.

Hereafter, only *z*-component of the thermal conductivity tensor is discussed; it corresponds to the one experimentally measured. A first observation is the continuous decrease of the membrane TC's as the porosity increases. In this case, a three-time reduction of $k_z$ is observed between pristine Si membrane ($k_f$=63.18 W m⁻¹ K⁻¹) and the PNC-membrane with highest porosities. A second remark concerns the impact of the PNC-membrane's pitch size *a* on the variation of TC. As expected, for a given porosity, reducing the pitch *a* implies to reduce pore diameter and
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_8.png" /><figcaption>Figure 7: Thermal conductivity tensor in a PNC membrane as a function of the pore radius; Membrane geometric characteristics $L_x=145$ nm, $L_y=5$ µm, $L_z=25$ µm, pitch of the PNC $a=300$ nm, $T=300$ nm. Experimental value are taken from M. Verdier et al.[27]</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Thermal conductivity of PNC membranes with aligned cylindrical pores versus porosity. Comparison to analytic modelling based on phonon mfp evaluation in thin films.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9: Thermal conductivity of PNC membranes with aligned cylindrical pores, comparison of MC-GK calculations to micro-TDTR measurements [27]</figcaption></figure>

thus to increase the number of pores in the membrane. As an example, when $a=160$ nm there are $30 \times 155$ cylindrical holes in the membrane with radius up to 63 nm, while when $a=500$ nm there are $9 \times 49$ holes in the membrane with radius up to 213 nm. Increasing the number of inner phonon scattering boundaries with large surface to volume ratio naturally (small $a$ case) improve phonon scattering efficiency and thus lead to a thermal conductivity lowering. Finally, we can notice that the proposed model for $k$ evolution with respect to geometrical parameters ($a$ and $R$) and phonon mfp in the membrane $\Lambda_f$ is consistent. It leads to an average $\Lambda_f=243$ nm, which is logically smaller than the $\Lambda_b$ computed before in the case of bulk system with spherical organized pores. We will show in the next section that similar value is found for staggered PNC membranes.

To conclude on this first set of calculations, we have compared our MC-GK simulations to micro-TDTR measurements. Results are reported in Fig. 9. In such a plot, perfect agreement supposes that all data-points are aligned on the first bisector line. Such agreement is good for pitches in the range $200 < a < 350$ nm, while light overestimation is found for the smallest pitch ($a = 160$ nm) and a light underestimation is observed for the largest one ($a = 500$ nm). Such trends were already observed in our former work based on the used of classical MC simulations (based on Fourier's formalism) with a temperature difference set along thr $z$ direction [27].

Eventually, it was shown that the mfp distribution can be extracted from the MC simulation of phonon packet transport and, in the case of porous systems, the distribution can be fitted by generalized Pareto model. For PNC membranes, such approach is still valid and is used to evaluate the mfp related to intrinsic and boundary scattering events. Computing the average phonon mean free path for the 66 different simulation cases, we provide variation of simulated TC's with respect to their assessed mfp, see Fig. 10. From log-scale plot given in Fig. 10, one can see that the thermal conductivity also scales with the mfp derived from Pareto fitting of phonon "jumping length" between each scattering events. The general trend is also $k \propto (mfp)^\beta$ with $\beta = 2.57$, close to 2.5, which is smaller than the previously considered test case (bulk with spherical pore). This variation can be related to the strong downsizing of the structure along $x$-direction ($L_x \ll L_y, L_z$).

## 3.3 Staggered phononic membranes

Before concluding, we demonstrate the ability of the method to deal with other types of geometry, staggered PNC membranes have been considered. Thin film characteristic remain similar to the previously investigated membrane ($L_z=25$ µm, $L_y=5$ µm and $L_x=145$ nm). Pitches for hollow cylinders spacing are also similar ([160, 200, 350, 500] nm), but they are organized in staggered mode with a unit cell characterized by an equilateral triangle (see Fig. 2-right). For this kind of membranes, 52 distinct cases were investigated with the MC-GK method at room temperature and compared to micro-TDTR experimental data reported in [26, 27].

In Fig. 11 are plotted the TC values calculated with MC-GK method for all the considered
<figure><img src="image_11.png" /><figcaption>Figure 10: Thermal conductivity as a function of average phonon mean free path in PNC membranes with aligned cylindrical pores. Power law fit for large pore diameters (red line).</figcaption></figure>

samples. Like in the case of the aligned PNC membrane $k_z$ values are decreasing as the porosity increases and this effect is more pronounced in membranes with smallest pitches. This is consistent with theory, experimental observations and our former calculations (see Fig. 8). In addition, the modelling procedure, detailed in section 2.2.2, has been also used to evaluate the phonon mfp in this membrane. Here, the average value of the thin film mfp is $\Lambda_f=240$ nm, a value close to the one obtained previously for aligned PNC.

In Fig. 12 are reported the $k_z$ values numerically predicted as a function of their experimental counterparts. Agreement between both approaches is reasonably good, except for the largest pitch ($a=500$ nm) where the simulation globally underestimates the experimental measurement. This trend was also noticed in the case of aligned PNC membranes and therefore probably need other simulations and experiments to assess possible limitation of the modelling and/or of the measurement techniques.

## 4 Conclusion

This work aims at demonstrating the ability of a new modelling procedure (Monte Carlo-Green Kubo) to assess thermal properties of complex porous systems. In comparison to former studies also based on the Monte Carlo solution of the Boltzmann transport equation developed in the group [16, 27], this method proposes improvements at different levels that are useful for nano and microscale thermal engineering. Among them, we can cite: the possibility to directly extract thermal conductivity tensor of complex systems, the ability to perform calculations at a prescribed temperature, the modelling efficiency in systems where phonon scattering is dominant. The latter point is a major advantage compared to classical MC simulations; i.e. calculations where thermal properties are extracted from steady state behaviour of a differentially heated structure. In classical MC approach, highly porous membranes are scarcely considered as convergence calculations, even on computing grids, are long to achieve. In addition to that, these techniques based on "quanta" tracking into the system of interest allow to extract other parameters like energy carriers mean free path or angular distribution of the latter. Here, generalized Pareto distribution of mfp fits the phonon travelling distance between resistive scattering events. It thus brings insights about phonon engineering in nanosystems. In this work, we also try to provide analytical models which can be used to evaluate thermal transport properties in "organized porous systems". Those models based on geometrical considerations and phonon boundary scattering allow to catch general tendencies and are consistent with phonon mean free path evaluation carried on bulk or thin films. Further developments of the simulation tool are currently in progress. Among them, the implementation of methods to deals with multi-component materials is considered in order to study sample with nanoinclusions that are also of big interest (e.g. superlattices, quantum dots, ...).
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_12.png" /><figcaption>Figure 11: Thermal conductivity of PNC membranes with staggered cylindrical pores versus porosity. Comparison to analytic modelling based on phonon mfp evaluation in thin films.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Thermal conductivity of PNC membrane with staggered cylindrical pores, comparison of MC-GK calculations to micro-TDTR measurements [27]</figcaption></figure>
<table><caption>Table 1: Data used to fit Si dispersion curves.</caption><thead><tr><th></th><th>units</th><th>Si</th></tr></thead><tbody><tr><td>c<sub>TA</sub></td><td>×10<sup>-7</sup>(m<sup>2</sup>s<sup>-1</sup>)</td><td>-2.28</td></tr><tr><td>v<sub>TA</sub></td><td>×10<sup>3</sup>(m s<sup>-1</sup>)</td><td>5.24</td></tr><tr><td>c<sub>LA</sub></td><td>×10<sup>-7</sup>(m<sup>2</sup>s<sup>-1</sup>)</td><td>-2.22</td></tr><tr><td>v<sub>LA</sub></td><td>×10<sup>3</sup>(m s<sup>-1</sup>)</td><td>9.26</td></tr></tbody></table>

# 5 Acknowledgements

This paper contains results obtained in the frames of the project “Spider-Man” (ANR-18-CE42-0006). This work was performed using HPC resources from GENCI-TGCC and GENCI-IDRIS (A0130913052), in addition HPC resources were partially provided by the EXPLOR center hosted by the Université de Lorraine. Thanks to “STOCK NRJ” that is co-funded by the European Union within the framework of the Program FEDER-FSE Lorraine and Massif des Vosges 2014–2020.

# A Appendices

## Dispersion properties

Parabolic fits $\omega = c_p K^2 + v_p K$ are used to describe the dispersion properties of Si, for LA and TA polarizations, as suggested by E. Pop et al [30]. Here, isotropic dispersion properties are considered. Parameters $c_p$ and $v_p$ are detailed in Tab. (1).

## Relaxation times

Relaxation time formulations are given below (Eq. A.1); each scattering process has its own expression which depends on the frequency and the temperature [31]. $B^i$ constants are adjusted to fit the experimental TC curve of bulk material conductivity as a function of the temperature $k = f(T)$, we cannot directly derive them from analytical expressions as phonon group velocities are not constant and calculated from phonon dispersion properties (quadratic fit). Tab. (2) lists them for “Normal”, “Umklapp” and “Impurity” scattering.

$$ 
\begin{align*}
(\tau_I)^{-1} &= B_I^i \omega^4 = \frac{V_0 \Gamma}{4 \pi v_{g,i}^3} \omega^4, \\
(\tau_N^L)^{-1} &= B_N^L \omega^2 T^3, \\
(\tau_N^T)^{-1} &= B_N^T \omega T^4, \\
(\tau_U^i)^{-1} &= B_U^i \omega^2 T e^{-\theta_i / (3T)}
\end{align*}
\tag{A.1} $$

<table><caption>Table 2: Relaxation time parameters for Si, Ge, GaN and C (diamond) for the Debye-Callaway model.</caption><thead><tr><th></th><th>units</th><th>Si</th></tr></thead><tbody><tr><td>B<sub>N</sub><sup>T</sup></td><td>×10<sup>-13</sup>(s K<sup>-3</sup>)</td><td>2.0</td></tr><tr><td>B<sub>U</sub><sup>T</sup></td><td>×10<sup>-20</sup>(s)</td><td>8.3</td></tr><tr><td>B<sub>N</sub><sup>L</sup></td><td>×10<sup>-24</sup>(K<sup>-4</sup>)</td><td>0.8</td></tr><tr><td>B<sub>U</sub><sup>L</sup></td><td>×10<sup>-20</sup>(K<sup>-4</sup>)</td><td>3.0</td></tr><tr><td>Γ</td><td>×10<sup>-5</sup>(-)</td><td>7.3</td></tr><tr><td>V<sub>0</sub></td><td>×10<sup>-30</sup>(m<sup>3</sup>)</td><td>40.9</td></tr><tr><td>θ<sub>T</sub></td><td>K</td><td>230.0</td></tr><tr><td>θ<sub>L</sub></td><td>K</td><td>583.8</td></tr></tbody></table>

## Analytic model of mfp, evaluation with ray-tracing methodology

Here, in addition to the comparison to available experimental data, we wanted to see if the proposed analytic modeling can be expected from direct considerations about phonon collision with boundaries. The proposed models given by equations (Eq. 11 & Eq. 13) only take into account the geometrical features of the porous media (pitch $a$ and pore radius $R$) through a “porous mfp” (labeled $\Lambda_p$). Such mean free path does not depend on the phonon spectral properties as it is purely
<figure><img src="image_14.png" /><figcaption>Figure 13: Pore scattering mean free path; comparison of MC evaluation (Eq. A.2) with analytic models (Eq. 11 and Eq. 13).</figcaption></figure>

geometric considerations that have been used to evaluate it. In order to demonstrate this, we have performed similar calculations to the ones presented in references [25, 16]. In these studies, the boundary scattering probability is calculated as the number of phonon packets that collide with a pore boundary $N_{pore}$ divided by the total number of phonon packets launched in the system $N_{tot}$. This collision probability can be related to a boundary lifetime and thus to a boundary mfp. It reads (see Eq. 2 in [25]):

$$P_{scat}(\omega, p) = 1 - \exp\left[\frac{-\delta t}{\tau_{pores}(\omega, p)}\right] = \frac{N_{pore}(\omega, p)}{N_{tot}(\omega, p)} \quad (A.2)$$

In the present work simulations parameters were: 10⁷ launched phonon packets, sampled over 500 frequency bands over the whole phonon spectra of silicon. Here, only one time step of $\delta t = 10$ ps is used. All calculations were performed for the organized spherical pore sample and for the two PNC membrane cases studied in the manuscript. Geometrical parameters of these cases were: i) 125 spheres in 1µm³ square box, ii) aligned PNC with pitch a=300nm, iii) staggered PNC with pitch a=350 nm. In both PNC cases pore radius was varied between 20 nm (smallest porosity, $\phi=1.2\%$) and 135 nm ($\phi=59\%$) or 150 nm ($\phi=65\%$) respectively for case ii) and iii). In case i) porosity range between 1% and 50%.

Hereafter, in Fig. 13, is plotted the pore boundary mfp for case i), ii) and iii) discussed above. On x axis is reported the theoretical $\Lambda_p$ computed with Eq. 11 and Eq. 13 (as a function of the pitch a, the radius R and the pore organization, organized, aligned or staggered), while on y axis is plotted the MC evaluated mfp $\Lambda_p$ obtained with Eq. A.2 which does not involve any adjustment parameter.

In Fig. 13 blue squares stand for spherical pore system, while orange diamond represent aligned PnC and green circle staggered one. In all cases MC calculation uncertainties are reported also. Perfect matching between model and MC simulations would result in all points aligned along the black dashed line ($x = y$). Naturally, this is not the case, especially for large porosities (small $\Lambda$ values), however the agreement between both approaches used to evaluate pore collision mfp is very good. This confirms the adequacy of the model to evaluate mfp due to pore scattering.

In addition, when porosity becomes large, the analytic model reaches its limit as pores are getting closer and thus the mathematical model based on a convex inclusion in a convex body ($\Lambda = 4V/S$) should take into account possible reflection at inner convex body boundary. In such case, surface scattering in the analytic model increases and thus theoretical collision mfp decreases. Such trend is consistent with our MC simulations that show smaller $\Lambda_{pMC}$ than $\Lambda_{pth}$ for large porosities. Those calculations confirm adequacy of the analytic modelling and give insights about its limitations.
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

# Phonon trajectories in PNC membranes

## References

[1] V Lysenko, S Perichon, B Remaki, and D Barbier. Thermal isolation in microsystems with porous silicon. *Sensors and Actuators A: Physical*, 99(1-2):13-24, 2002.

[2] David Song and Gang Chen. Thermal conductivity of periodic microporous silicon films. *Applied physics letters*, 84(5):687-689, 2004.

[3] Junki Nakagawa, Yuta Kage, Takuma Hori, Junichiro Shiomi, and Masahiro Nomura. Crystal structure dependent thermal conductivity in two-dimensional phononic crystal nanostructures. *Applied Physics Letters*, 107(2):023104, 2015.

[4] Maxime Verdier, David Lacroix, Stanislav Didenko, Jean-François Robillard, Evelyne Lampin, Thierno-Moussa Bah, and Konstantinos Termentzidis. Influence of amorphous layers on the thermal conductivity of phononic crystals. *Physical Review B*, 97(11):115435, 2018.

[5] Gilles Pernot, Mathieu Stoffel, Ivana Savic, Fabio Pezzoli, Peixuan Chen, Guillaume Savelli, A. Jacquot, Joachim Schumann, U. Denker, Ingolf Mönch, Ch. Deneke, Oliver G. Schmidt, Jean-Michel Rampnoux, Shidong Wang, M. Plissonnier, Armando Rastelli, Stefan Dilhaire, and Natalio Mingo. Precise control of thermal conductivity at the nanoscale through individual phonon-scattering barriers. *Nature Materials*, 9(6):491-495, 2010.

[6] Yaolan Tian, Tuomas A Puurtinen, Zhuoran Geng, and Ilari J Maasilta. Minimizing coherent thermal conductance by controlling the periodicity of two-dimensional phononic crystals. *Physical Review Applied*, 12(1):014008, 2019.

[7] Sergei Gluchko, Roman Anufriev, Ryoto Yanagisawa, Sebastian Volz, and Masahiro Nomura. On the reduction and rectification of thermal conduction using phononic crystals with pump-shaped holes. *Applied Physics Letters*, 114(2):023102, 2019.

[8] Marianna Sledzinska, Bartlomiej Graczykowski, Jeremie Maire, Emigdio Chavez-Angel, Clivia M Sotomayor-Torres, and Francesc Alzina. 2d phononic crystals: Progress and prospects in hypersound and thermal transport engineering. *Advanced Functional Materials*, 30(8):1904434, 2020.

[9] Thomas Vasileiadis, Jeena Varghese, Visnja Babacic, Jordi Gomis-Bresco, Daniel Navarro Uríos, and Bartlomiej Graczykowski. Progress and perspectives on phononic crystals. *Journal of Applied Physics*, 129(16):160901, 2021.

[10] Masahiro Nomura, Roman Anufriev, Zhongwei Zhang, Jeremie Maire, Yangyu Guo, Ryoto Yanagisawa, and Sebastian Volz. Review of thermal transport in phononic crystals. *Materials Today Physics*, page 100613, 2022.

[11] Roman Anufriev, Yunhui Wu, Jose Ordonez-Miranda, and Masahiro Nomura. Nanoscale limit of the thermal conductivity in crystalline silicon carbide membranes, nanowires, and phononic crystals. *NPG Asia Materials*, 14(1):35, 2022.

[12] Antonin M Massoud, Valeria Lacatena, Maciej Haras, Emmanuel Dubois, Stéphane Monfray, Jean-Marie Bluet, Pierre-Olivier Chapuis, and Jean-François Robillard. Heat dissipation in partially perforated phononic nano-membranes with periodicities below 100 nm. *APL Materials*, 10(5):051113, 2022.

[13] Nobuyuki Zen, Tuomas A Puurtinen, Tero J Isotalo, Saumyadip Chaudhuri, and Ilari J Maasilta. Engineering thermal conductance using a two-dimensional phononic crystal. *Nature communications*, 5(1):1-9, 2014.

[14] David Lacroix, Mykola Isaiev, and Gilles Pernot. Thermal transport in semiconductors studied by monte carlo simulations combined with the green-kubo formalism. *Phys. Rev. B*, 104:165202, Oct 2021.

[15] David Lacroix, Karl Joulain, and Denis Lemonnier. Monte carlo transient phonon transport in silicon and germanium at nanoscales. *Phys. Rev. B*, 72:064305, Aug 2005.
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_15.png" /><figcaption>Figure 14: Phonon paths in an aligned PNC membrane, trajectories of 5 distinct energy carriers; Membrane geometric characteristics $L_x=145$ nm, $L_y=5$ µm, $L_z=25$ µm, pitch of the PNC $a=500$ nm, hole radius $r=110$ nm.</figcaption></figure>
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_16.png" /><figcaption>Figure 15: Phonon paths in an staggered PNC membrane, trajectories of 5 distinct energy carriers; Membrane geometric characteristics $L_x=145$ nm, $L_y=5$ µm, $L_z=25$ µm, pitch of the PNC $a=500$ nm, hole radius $r=110$ nm.</figcaption></figure>
[16] V. Jean, S. Fumeron, K. Termentzidis, S. Tutashkonko, and D. Lacroix. Monte carlo simulations of phonon transport in nanoporous silicon and germanium. *Journal of Applied Physics*, 115(2):024304, 2014.

[17] Jean-Philippe M Péraud, Colin D Landon, and Nicolas G Hadjiconstantinou. Monte carlo methods for solving the boltzmann transport equation. *Annual Review of Heat Transfer*, 17, 2014.

[18] Marianna Sledzinska, Bartlomiej Graczykowski, Francesc Alzina, Umberto Melia, Konstantinos Termentzidis, David Lacroix, and Clivia M Sotomayor Torres. Thermal conductivity in disordered porous nanomembranes. *Nanotechnology*, 30(26):265401, apr 2019.

[19] Patrick K. Schelling, Simon R. Phillpot, and Pawel Keblinski. Comparison of atomic-level simulation methods for computing thermal conductivity. *Phys. Rev. B*, 65:144306, Apr 2002.

[20] Joseph Callaway. Model for lattice thermal conductivity at low temperatures. *Phys. Rev.*, 113:1046-1051, Feb 1959.

[21] M. G. Holland. Analysis of lattice thermal conductivity. *Phys. Rev.*, 132:2461-2471, Dec 1963.

[22] Liang-Chun Liu and Mei-Jiau Huang. Thermal conductivity modeling of micro- and nanoporous silicon. *International Journal of Thermal Sciences*, 49(9):1547-1554, 2010.

[23] S. Blanco and R. Fournier. An invariance property of diffusive random walks. *Europhys. Lett.*, 61(2):168, jan 2003.

[24] Alain Mazzolo, Benoit Roesslinger, and Wilfried Gille. Properties of chord length distributions of nonconvex bodies. *Journal of Mathematical Physics*, 44(12):6195-6208, 2003.

[25] David Lacroix, Karl Joulain, Damian Terris, and Denis Lemonnier. Monte Carlo simulation of phonon confinement in silicon nanostructures: Application to the determination of the thermal conductivity of silicon nanowires. *Applied Physics Letters*, 89(10), 09 2006. 103104.

[26] Roman Anufriev, Jeremie Maire, and Masahiro Nomura. Reduction of thermal conductivity by surface scattering of phonons in periodic silicon nanostructures. *Phys. Rev. B*, 93:045411, Jan 2016.

[27] Maxime Verdier, Roman Anufriev, Aymeric Ramiere, Konstantinos Termentzidis, and David Lacroix. Thermal conductivity of phononic membranes with aligned and staggered lattices of holes at room and low temperatures. *Phys. Rev. B*, 95:205438, May 2017.

[28] Maxime Verdier, Konstantinos Termentzidis, and David Lacroix. Crystalline-amorphous silicon nano-composites: Nano-pores and nano-inclusions impact on the thermal conductivity. *Journal of Applied Physics*, 119(17):175104, 2016.

[29] M. Isaiev, Y. Mankovska, V. Kuryliuk, and D. Lacroix. Thermal transport properties of nanoporous silicon with significant specific surface area. *Applied Physics Letters*, page to be published, 2023.

[30] Eric Pop, Robert W Dutton, and Kenneth E Goodson. Monte carlo simulation of joule heating in bulk and strained silicon. *Applied Physics Letters*, 86(8):082101, 2005.

[31] D. T. Morelli, J. P. Heremans, and G. A. Slack. Estimation of the isotope effect on the lattice thermal conductivity of group iv and group iii-v semiconductors. *Phys. Rev. B*, 66:195304, Nov 2002.
<figure><img src="image_17.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.55582

+

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>
<figure><img src="image_23.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_24.png" /></figure>
<figure><img src="image_25.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_26.png" /></figure>
<figure><img src="image_27.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

Probability

<figure><img src="image_28.png" /><figcaption>Probability Plot - mfp distribution @300K in : Si</figcaption></figure>
<figure><img src="image_29.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

k (W/m.K)

<figure><img src="image_30.png" /></figure>
<figure><img src="image_31.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_32.png" /></figure>
<figure><img src="image_33.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_34.png" /></figure>
<figure><img src="image_35.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_36.png" /></figure>
<figure><img src="image_37.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_38.png" /></figure>
<figure><img src="image_39.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_40.png" /></figure>
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_41.png" /></figure>
<figure><img src="image_42.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_43.png" /></figure>
<figure><img src="image_44.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_45.png" /></figure>
<figure><img src="image_46.png" /></figure>

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0155582

<figure><img src="image_47.png" /></figure>

×10⁻⁵