# Weak decays of heavy-quark baryons

Janina Nicolini
On behalf of the LHCb collaboration  TU Dortmund University, Dortmund, Germany
Université Paris-Saclay, Orsay, France

###### Abstract

Weak decays of heavy-quark baryons offer an attractive laboratory to search
for effects beyond the Standard Model (SM), complementary to searches in meson
decays. However, the properties of these baryons are not well understood. In
particular, the $\Omega^{-}_{b}$ baryon is the least studied weakly decaying
$b$ baryon. Therefore, the most precise determination of the $\Omega^{-}_{b}$
baryon mass is reported, using $\Omega^{-}_{b}\to J/\psi\Omega^{-}$ decays at
the LHCb experiment. In addition, the relative production rate of the
$\Omega^{-}_{b}$ baryon in $pp$ collisions at the LHC is determined for the
first time. These quantities are crucial to predict and measure branching
fractions of $\Omega^{-}_{b}$ baryon decays.

##  1 Precision measurements of weakly decaying $b$ baryons

Weak decays of heavy-quark baryons provide a complementary laboratory to
search for effects beyond the SM and offer rich angular structures compared to
meson decays. These are particularly interesting as they allow to test spin
dependencies of possible new physics scenarios. The quark model predicts four
weakly decaying $b$ baryons with one heavy quark: the lightest one, the
$\Lambda_{b}^{0}=|bud\rangle$, the isospin-doublet of
$\Xi_{b}^{0}=|bsu\rangle$ and $\Xi_{b}^{-}=|bsd\rangle$ and the heaviest one,
the $\Omega_{b}^{-}=|bss\rangle$. The experimental precision of the measured
$b$ baryon properties scales largely with the production fraction $f_{i}$ of
each $b$ hadron and therefore the available statistics. While the probability
$f_{i}$ of a $b\bar{b}$ pair hadronising into an $\Lambda_{b}^{0}$ baryon is
about $f_{\Lambda_{b}^{0}}\sim 18\%$ [1], for each $\Xi_{b}$ baryon the
production fraction is only approximately $f_{\Xi_{b}}\sim 1.5\%$ [1]. For the
$\Omega_{b}^{-}$ baryon, the production fraction is assumed to be
$f_{\Omega_{b}^{-}}\sim 0.5\%$ [1], but this assumption is based on a
measurement by the CDF experiment [2] with only 16 $\Omega_{b}^{-}$ events.
The low expected production probability makes the $\Omega_{b}^{-}$ difficult
to study with respect to other weakly decaying $b$ baryons.
In order to search for effects beyond the SM, the properties of the $b$
baryons have to be known with high precision. The production fractions of
these baryons need to be known to measure the absolute branching fractions of
their decays. As a first step, the relative production fraction of the
$\Omega_{b}^{-}$ baryon has been determined relatively to the $\Xi_{b}^{-}$
baryon for the first time at the Large Hadron Collider (LHC) by the LHCb
experiment:

$R=\frac{f_{\Omega_{b}^{-}}}{f_{\Xi_{b}^{-}}}\times\frac{\mathcal{B}(\Omega_{b}^{-}\to
J/\psi\Omega^{-})}{\mathcal{B}(\Xi_{b}^{-}\to J/\psi\Xi^{-})}.$ (1)

These decay modes were chosen as they both provide a decay chain of weakly
decaying particles that cannot be mimicked by any mesonic decays, leading to
very low background contamination. The two final states of the
$\Omega_{b}^{-}$ and $\Xi_{b}^{-}$ decay only differ by one particle since the
$\Omega^{-}$ and the $\Xi^{-}$ are reconstructed in their decays to, $\Lambda
K^{-})$ and ${\Lambda\pi^{-}}$, respectively. Given this particular topology,
many systematic uncertainties cancel in the relative production fraction ratio
$R$ (Eq. 1). In addition, both decays have rather low energy releases making
them excellent candidates to measure the $\Omega_{b}^{-}$ baryon mass through
the mass difference with the $\Xi_{b}^{-}$ baryon. Since the $\Omega^{-}$
baryon has shorter lifetime compared to that of the $\Xi^{-}$, more
$\Omega_{b}^{-}$ decays happen within the LHCb acceptance. Therefore, it has a
slightly better detection efficiency compared to the $\Xi_{b}^{-}$ baryon.
The measurement [3] to determine the relative production fraction ratio $R$
and the $\Omega_{b}^{-}$ mass has been performed with a fully cut-based
selection in a fiducial phase-space.

##  2 The $\Omega_{b}^{-}$ mass measurement

The full Run 1 and Run 2 LHCb dataset corresponding to an integrated
luminosity of $9\;\mathrm{fb}^{-1}$ is used to determine the $\Omega_{b}^{-}$
baryon mass. Instead of directly extracting the mass, the mass difference
${m(\Omega_{b}^{-})-m(\Xi_{b}^{-})}$ is measured. This choice allows to cancel
the dominant systematic uncertainty due to the knowledge of the absolute
momentum scale. The mass difference is kept as a free parameter and extracted
from a simultaneous likelihood fit shown in Fig. 1.

x

Figure 1: Invariant-mass distributions in the (top) $\Xi_{b}^{-}\to
J/\psi\Xi^{-}$ and (bottom) $\Omega_{b}^{-}\to J/\psi\Omega^{-}$ datasets, in
(left) Run 1 and (right) Run 2 data-taking periods used for the mass-
difference measurement. (Taken from Ref.[3])

The resulting value is

$m(\Omega_{b}^{-})-m(\Xi_{b}^{-})=248.54\pm 0.51(\mathrm{stat})\pm
0.38(\mathrm{syst})\;\mathrm{MeV}/c^{2}.$ (2)

The dominant source for the systematic uncertainty is related to the poor
knowledge of the hyperon masses. Together with the most precise determination
of the $\Xi_{b}^{-}$ baryon mass [4] using the decay
$\Xi_{b}^{-}\to\Xi_{c}^{0}\pi^{-}$ from the LHCb experiment, with a value of
$m(\Xi_{b}^{-})=5797.33\pm 0.24\;(\mathrm{stat.})\pm
0.29\;(\mathrm{syst.})\;\mathrm{MeV}/c^{2}$, the $\Omega_{b}^{-}$ mass is
determined to be

$m(\Omega_{b}^{-})=6045.9\pm 0.5\;(\mathrm{stat.})\pm
0.6\;(\mathrm{syst.})\;\mathrm{MeV}/c^{2}.$ (3)

This measurement provides an improvement of more than a factor of two in
precision with respect to the previous best measurement [5].

##  3 The relative production fraction $R$

Due to the energy dependence of the production fractions and the limited
statistical power of the Run 1 dataset, the relative production fraction $R$
(see Eq. 1) is only determined through a simultaneous likelihood fit to the
Run 2 dataset. This sample corresponds to an integrated luminosity of
$6\;\mathrm{fb}^{-1}$ recorded at a centre-of-mass energy of
$13\;\mathrm{TeV}$. The selection is optimised to improve the sensitivity for
the production fraction, and the corresponding mass fits are shown in Fig. 2.

Figure 2: Invariant-mass distributions in the (left) $\Xi_{b}^{-}\to
J/\psi\Xi^{-}$ and (right) $\Omega_{b}^{-}\to J/\psi\Omega^{-}$ Run 2 datasets
used for the relative production estimation. (Taken from Ref.[3])

The relative production fraction is kept as a free fit parameter and measured
to be

$R=0.120\pm 0.008(\mathrm{stat})\pm 0.008(\mathrm{syst}).$ (4)

The dominant systematic uncertainties are due to simulation calibration linked
to the small size of the $\Omega_{b}^{-}$ dataset and the poor knowledge of
the $b$ baryon lifetimes. The measurement is in agreement within $1.2\sigma$
with the CDF measurement [6] but the LHCb measurement suggests a larger
relative production. In addition to the statistical limitation of the CDF
measurement with $16$ $\Omega_{b}^{-}$ events, the larger nominal value can be
related to the different production environments at the Tevatron and the LHC,
leading to differing production cross-sections.

##  4 Conclusion

These proceedings present the most precise measurement of the $\Omega_{b}^{-}$
baryon mass using the LHCb dataset corresponding to an integrated luminosity
of $9\;\mathrm{fb}^{-1}$. The measurement agrees with the world average and
previous measurements, which are displayed in Fig. 3.

Figure 3: Overview of the $\Omega_{b}^{-}$ mass measurements to date (red
points) [2, 5, 6, 7, 8], with corresponding statistical and systematic
uncertainties marked in red and blue, respectively, as well as the current
world average by the Particle Data Group [9] (light-blue band). The superseded
measurements are shown with lighter colours. The new LHCb average is presented
in orange, showing the total uncertainty. (Taken from Ref. [3])

In addition, the first determination of the relative production fraction of
the $\Omega_{b}^{-}$ baryon at the LHC at a centre-of-mass energy of
$\sqrt{s}=13\;\mathrm{TeV}$ was presented

$\frac{f_{\Omega_{b}^{-}}}{f_{\Xi_{b}^{-}}}\times\frac{\mathcal{B}(\Omega_{b}^{-}\to\Omega_{b}^{-}J/\psi)}{\mathcal{B}(\Xi_{b}^{-}\to\Xi_{b}^{-}J/\psi)}=0.120\pm
0.008(\mathrm{stat.})\pm 0.008(\mathrm{syst.}),$ (5)

using a dataset of $6\;\mathrm{fb}^{-1}$. Further theoretical input on the
branching fractions is needed to disentangle the ratio of production fractions
from the ratio of branching fractions and enable to perform a first absolute
branching fraction measurement of $\Omega_{b}^{-}$ baryon decays.

## Acknowledgments

J.N. acknowledges support from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme under grant
agreement No. 714536: PRECISION, the German Federal Ministry of Education and
Research (BMBF, grant no. 05H21PECL1) within ErUM-FSP T04, and the German
Academic Scholarship Foundation.

## References

## References

  * [1] Aaij, R. et al. (LHCb collaboration), Measurement of the $b$-Quark Production Cross Section in 7 and 13 TeV $pp$ Collisions, Phys. Rev. Lett. 118, 052002 (2017).
  * [2] Aaltonen, T. et al. (CDF1 collaboration), Observation of the $\Omega_{b}^{-}$ baryon and measurement of the properties of the $\Xi_{b}^{-}$ and $\Omega_{b}^{-}$ baryons. Phys. Rev. D 80, 7 (072003)2009.
  * [3] Aaij, R. et al. (LHCb collaboration), Measurement of the mass difference and relative production rate of the $\Omega_{b}^{-}$ and $\Xi_{b}^{-}$ baryons. LHCb-PAPER-2022-053 (in preparation).
  * [4] Aaij, R. et al. (LHCb collaboration), Observation of a new ${\mathrm{\Xi}}_{b}^{0}$ state. Phys. Rev. D 103, 012004 (2021).
  * [5] Aaij, R. et al. (LHCb collaboration), Observation of excited ${\mathrm{\Omega}}_{c}^{0}$ baryons in ${\mathrm{\Omega}}_{b}^{-}\rightarrow{\mathrm{\Xi}}_{c}^{+}{K}^{-}{\pi}^{-}$ decays. Phys. Rev. D 104, L091102 (2021).
  * [6] Aaltonen, T. A. et al. (CDF2 collaboration), Mass and lifetime measurements of bottom and charm baryons in $p\bar{p}$ collisions at $\sqrt{s}=1.96$ TeV. Phys. Rev. D 89, 7 (2014).
  * [7] Aaij, R. et al. (LHCb collaboration), Measurements of the $\Lambda_{b}^{0}$, $\Xi_{b}^{-}$, and $\Omega_{b}^{-}$ baryon masses. Phys. Rev. Lett. 110, 182001 (2013).
  * [8] Aaij, R. et al. (LHCb collaboration), Measurements of the mass and lifetime of the $\Omega_{b}^{-}$ baryon. Phys. Rev. D 93, 092007 (2016).
  * [9] Workman, R. L. et al., PDG - Review of particle physics. Prog. Theor. Exp. Phys. 2022, 8 (2022).


# Differentially Encoded Observation Spaces
for Perceptive Reinforcement Learning

Lev Grossman1 and Brian Plancher2 1Lev Grossman is with Berkshire Grey,
Bedford, MA, USA<EMAIL_ADDRESS>Plancher is with
Barnard College, Columbia University, New York, NY, USA<EMAIL_ADDRESS>
###### Abstract

Perceptive deep reinforcement learning (DRL) has lead to many recent
breakthroughs for complex AI systems leveraging image-based input data.
Applications of these results range from super-human level video game agents
to dexterous, physically intelligent robots. However, training these
perceptive DRL-enabled systems remains incredibly compute and memory
intensive, often requiring huge training datasets and large experience replay
buffers. This poses a challenge for the next generation of field robots that
will need to be able to learn on the edge in order to adapt to their
environments. In this paper, we begin to address this issue through
differentially encoded observation spaces. By reinterpreting stored image-
based observations as a video, we leverage lossless differential video
encoding schemes to compress the replay buffer without impacting training
performance. We evaluate our approach with three state-of-the-art DRL
algorithms and find that differential image encoding reduces the memory
footprint by as much as 14.2$\times$ and 16.7$\times$ across tasks from the
Atari 2600 benchmark and the DeepMind Control Suite (DMC) respectively. These
savings also enable large-scale perceptive DRL that previously required paging
between flash and RAM to be run entirely in RAM, improving the latency of DMC
tasks by as much as 32%.

##  I INTRODUCTION

With its ability to synthesize complex behaviors in both simulated and real
environments, deep reinforcement learning (DRL) has been applied to solve a
host of robotic-specific problems ranging from dexterous manipulation [1], to
quadrupedal locomotion [2], to high-speed drone racing [3], as well as more
general artificial intelligence (AI) feats of mastering tabletop [4] and
competitive video games such as Dota 2 [5] and Minecraft [6].

Despite these achievements, DRL remains very sample inefficient, often
requiring massive amounts of training data to learn. Because much of this
persistent data is loaded into experience replay buffers during training, DRL
is extremely memory intensive, limiting the number of computational platforms
that can support such operations, and largely confining state-of-the-art model
training to the cloud.

Interestingly, it is the observations in the replay buffer which often consume
the majority of the memory, reaching over 90% for proprioceptive models [7].
The recent interest in perceptive DRL models [6, 8] only further exacerbates
this problem. For example, as compared to the 224 bytes needed to store the
proprioceptive observation space for the DeepMind Control Suite (DMC) [9]
quadruped, a single 84x84 grayscale image requires over 7kB of memory.
Consequently, deployment to the edge is unrealistic for such approaches.
Still, many robots, especially those involved in tasks as consequential as
search-and-rescue and space exploration [10], will have to adapt to ever-
changing environmental conditions and continue to optimize and update their
internal policies over the course of their lifetime [11], often in remote
areas without access to fast network connections. As such, approaches that
reduce the overall memory footprint of perceptive DRL are needed to enable
edge deployments.

In this paper, we begin to address this issue through differentially encoded
observation spaces. That is, by reinterpreting stored image-based observations
as a video, we leverage lossless differential video encoding schemes to
compress the replay buffer without impacting training performance. We evaluate
our approach across ten Atari 2600 benchmark tasks [12] as well as two robotic
control tasks from within the DMC, using three state-of-the-art DRL
algorithms. We find that differential image encoding reduces the memory
footprint by as much as 14.2$\times$ and 16.7$\times$ for the Atari and DMC
tasks respectively. These savings also enable large-scale DRL that previously
required paging between flash and RAM to be run entirely in RAM, improving
latency for the DMC tasks by as much as 32%. We release our software and
experiments open-source at: github.com/A2R-Lab/DiffCompressDRL.

##  II RELATED WORK

Perceptive DRL: Perceptive or image-based control and deep reinforcement
learning (DRL) share common roots. First proposed a decade ago, Deep
Q-Networks (DQNs) [13] bridged the gap between deep learning and RL, achieving
super-human performance on many Atari 2600 video games using only visual
input. Many model-free, DQN-based variants have followed [14, 15], further
improving game-playing performance directly from pixel observations. Still,
these early model-free methods remained quite data expensive. To address this,
some turned to model-based or world-model learning in both the Atari [16] as
well as robotic control [6, 17] domains. More recently, model-free methods
have improved sample efficiency further through contrastive learning [18] and
image augmentation [19]. Despite these improvements, perceptive DRL remains
memory intensive, owing especially to the large size of image observations.
For memory constrained systems and more complex tasks that require larger
replay buffers, it has become necessary to store observations directly on
disk, increasing the average memory access time and thus overall model
training times [8].

Compressed Perceptive DRL: While recent image-based, model-free DRL methods
have improved performance and sample efficiency through data augmentation
[20], contrastive learning [18], and more robust latent observation space
representations [21], little work has been done (aside from moving from RGB to
grayscale [22]) to reduce the size of the stored image observations. Model-
based methods offer an alternative by learning a world-model directly, thereby
reducing the dependence on large amounts of replay experience [23, 24]. While
these methods do show promise in highly complex simulated environments [6] as
well as in some real laboratory settings [25], world-model learning tends to
still struggle when the underlying world dynamics become highly complex [26].
For this reason, model-free methods remain popular for robotics applications
leveraging learning in the real world [27].

Video Encoding for Learning: Traditionally, video encoding schemes aim to
reduce the visual redundancy in digital video files without impacting human-
perceived quality [28]. Modern video encoding standards such as MPEG-4 [29]
and H.264 [30] use both intra- and inter-picture compression techniques to
achieve this result. With the recent proliferation of deep learning for
computer vision, some have begun optimizing traditional video encoding
standards for learning frameworks [31]. While these methods have been applied
to classification [32], object detection [33], and instance segmentation [34]
tasks among others, little has been done to incorporate video compression
within DRL.

##  III BACKGROUND

###  III-A Reinforcement Learning

Reinforcement learning poses problems as Markov decision processes (MDPs),
where an MDP is defined by a set of observed and hidden states, $\mathcal{S}$,
actions, $\mathcal{A}$, stochastic dynamics, $p(s_{t+1}|s_{t},a_{t})$, a
reward function $r(s,a)$, and a discount factor, $\gamma$. The RL objective is
to compute the policy, $\pi^{*}(s,a)$, that maximizes the expected discounted
sum of future rewards, $\mathbb{E}_{s,a}\left(\sum_{t}\gamma^{t}r_{t}\right)$.

Perceptive RL, where states $s$ are represented by images, adds an additional
caveat to the RL problem formulation, as representing a state using a single
image rendering of the system is often not sufficiently descriptive [8]. To
alleviate this issue, states are approximated by concatenating together the
$f$ immediately preceding image frames $s_{t}=\\{o_{t-f+1},\dots,o_{t}\\}$,
where $o_{t}$ is the image observation at timestep $t$ [13].

DRL algorithms parameterize the solution to these large MDPs with neural
networks. In this work, we leverage three different state-of-the-art DRL
algorithms: Proximal Policy Optimization [35], Quantile Regression DQN [15],
and Data-regularized Q-v2 [8]. In the remainder of this section, we provide
more detail on each algorithm.

####  III-A1 Proximal Policy Optimization (PPO)

Proximal Policy Optimization [35] is an on-policy, policy gradient [36]
algorithm that employs an actor-critic framework to learn both the optimal
policy, $\pi^{*}(s,a)$, as well as the optimal value function, $V^{*}(s)$.
Both the policy $\pi_{\theta}$ and value function $V_{\phi}$ are parameterized
by neural networks with weights $\theta$ and $\phi$ respectively. During
training, PPO needs to store the parameters of both its value and policy
networks111Many PPO implementations save space by using one central MLP with
two additional single-layer policy and value function model heads. We focus on
the standard PPO model approach in this section for clarity.–usually shallow
multilayer perceptrons (MLPs)–as well as its on-policy rollout buffer
$\mathcal{D}_{k}$. $\mathcal{D}_{k}$ stores $(s,a,r)$ tuples, which are
refreshed during each iteration of the algorithm.

####  III-A2 Quantile Regression DQN (QR-DQN)

Quantile Regression DQN [15] improves upon the highly influential Deep
Q-Networks (DQN) [13] algorithm, which uses experience replay to learn an
optimal Q-function approximated by a neural network $Q_{\theta}$. While in the
original DQN paper the authors train using a Huber loss [37] to maximize mean
return, QR-DQN uses a custom quantile Huber loss in order to learn a
distribution over expected returns. With this change, QR-DQN is able to
significantly outperform DQN on a benchmark of Atari 2600 games. Unlike PPO,
QR-DQN stores collected experience across all iterations in a much larger
replay buffer $\mathcal{D}$, functioning as a size limited queue.

####  III-A3 Data-regularized Q-v2 (DrQ-v2)

Data-regularized Q-v2 [8] improves upon the Data-regularized Q [19] model-free
RL algorithm. Both algorithms use data augmentation in conjunction with an
actor-critic learning method in order to perform image-based continuous
control. DrQ-v2 applies image augmentation through random shifts and then
passes the augmented images through an encoder network $f_{\xi}$ to produce
the final observations, which are used by an actor-critic DDPG [38] setup to
learn the solution policy. DrQ-v2 achieves state-of-the-art model-free results
on the DeepMind Control Suite, rivaling the performance of the popular model-
based DreamerV2 [24], while doing so 4$\times$ faster (in terms of wall-clock
training time). Like QR-DQN, DrQ-v2 stores its collected experience data in a
large replay buffer $\mathcal{D}$.

Figure 1: A graphical example of a differential image encoding taken from the
Breakout Atari 2600 environment [12]. The two images, A and B shown on the
left and middle, are very similar. The only difference in the resulting B-A
image, shown on the right, occurs for the two different locations of the
falling object. We can therefore encode this small difference in a highly
compressed sparse matrix format.

###  III-B Memory Requirements for Perceptive DRL

Previous work showed that for proprioceptive DRL, the memory requirements for
PPO were dominated by the observation space, consuming over 90% of the total
memory [7]. Perceptive DRL does not change this and actually only makes
matters worse. For example, the DMC quadruped requires only 224 bytes for a
single proprioceptive observation, while a single 84x84 grayscale image
requires over 7kB of memory. For the off-policy QR-DQN as well as DrQ-v2, an
even greater percentage of the overall memory requirement is dominated by the
observation space since, as mentioned above, their replay buffers are
generally orders of magnitude larger than PPO’s rollout buffer.

###  III-C Differential Video Encoding

Modern video coding formats and codecs, like the H.264 [30] standard, make use
of a number of compression and quantization techniques in order to reduce the
overall size of video streams while preserving quality. One of the compression
techniques employed is differential encoding (DE) [39], which can be used to
compress both single images as well as multiple subsequent image frames by
encoding the difference in inter-pixel or pixel-block values/intensities. In
this work, we make use of a simplified version of differential video
encoding–with inspiration from standard text-based DE used by most version
control systems [40]–in order to compress the pixel-based observations found
in rollout and replay buffers. We present our encoding scheme in detail in the
following section.

##  IV METHOD

In this section, we detail a custom, lossless compression technique for image-
based observation spaces. While we later show empirical results using PPO, QR-
DQN, and DrQ-v2, we note that this differential encoding-based scheme can be
applied to any image-based RL method in order to reduce the size of its stored
observations.

###  IV-A Differential Encoder

During DRL training, there tends to be a level of temporal sortedness in
experience replay buffers [41, 42], owing to the fact that many model-free RL
algorithms add experience $(s,a,r)$ sequentially as they train. Therefore,
while two states $s_{i},s_{j}$ drawn at random from a replay buffer may differ
dramatically, two adjacent states $s_{i-1},s_{i}$ have a high likelihood of
being similar, as they may originate from the same episode, just one control
step apart. Applied to image-based observations, this means that images
$o_{i-1}$ and $o_{i}$ will tend to be visually quite similar.

This similarity between adjacent image frames allows for a natural
differential encoding-based compression scheme (Figure 1), which takes the
difference between $o_{i}-o_{i-1}$ and stores the result in sparse matrix
format $S_{i}=\texttt{SparseArray}(o_{i}-o_{i-1})$. To decode $o_{i}$, it is
then only necessary to store $o_{i-1}$ and $S_{i}$, as:

$o_{i}=o_{i-1}+S_{i}.$ (1)

In order to efficiently store and access these image differences, we define a
custom sparse matrix implementation SparseArray that stores all nonzero pixels
using two arrays inds and vals, where inds stores pixel locations and vals
stores pixel values. Because image observations in both the Atari 2600
benchmark and DeepMind Control Suite are stored as 84x84, 8-bit unsigned
integer arrays (assuming grayscale images), we let inds store 8-bit unsigned
integers and vals store 16-bit signed integers. Assuming $o_{i}-o_{i-1}$
results in $n$ nonzero pixel values, we can store a 7kB input image in $4n$
bytes.

While it is possible that $o_{i}-o_{i-1}$ results in an image with a majority
nonzero pixel values, in practice we find the number of nonzero pixels $n$ to
be small. Nevertheless, we handle this corner case explicitly in SparseArray.
If $4n>84\times 84$, we store the full image array in order to cap the maximum
image size at 7kB.

###  IV-B Observation Indexing

As noted in the previous section, image-based learning methods generally store
states $s$ as a stack of the $f$ prior image observations. Letting $f=4$ (as
is the case for our Atari 2600 experiments in Section V), this results in each
state taking up $84\times 84\times 4=28$kB of memory.

Instead of storing states in this full image stack format, we define two
separate arrays obs, which stores the individual image observations, and
obs_inds, which stores pointers into obs for each state. To reconstruct the
full state $s_{i}$, we simply concatenate the individual observations
referenced by $\texttt{obs\\_inds}_{i}$. In most cases
$\texttt{obs\\_inds}_{i}=\\{(i-f+1),(i-f+2),\dots,i\\}$ resulting in:

$s_{i}=\\{\texttt{obs}_{(i-f+1)},\texttt{obs}_{(i-f+2)},\dots,\texttt{obs}_{i}\\}.$
(2)

However, in various corner cases, this is not the case, necessitating the use
of the pointer array. For example, for the first few observations, at the
start of an episode, individual observations need to be repeated to fill a
history of size $f$.

Through observation indexing alone, we can reduce the memory required to store
all states by almost $f\times$, as $f$ pointers are a fraction of the size of
the $(f-1)$ image observations we no longer need to save.

###  IV-C Observation Compressor

With both the differential image encoder and smart observation indexing, we
can define an observation compressor to handle the storage and retrieval of
image-based observations. Following best practices from the video encoding
literature, we apply our approach to each set of $f$ images.

Our compressor, ObsComp, internally stores three arrays: obs, sparse_obs, and
obs_inds. Using these arrays we can then define the following two operations:
ObsComp.get(i) and ObsComp.set(s, i), which get and set the uncompressed state
$s_{i}$ respectively.

To set a state $s_{i}$, we first compute the observation index array
$\texttt{obs\\_inds}_{i}=\\{(i-f+1),\dots,i\\}$. Next, we check if $i\mod
f=0$, and if so store the raw image frame $\texttt{obs}_{i/f}=o_{i}$. If not,
we compute and store the SparseArray, $\texttt{sparse\\_obs}_{i/f,(i\mod
f)-1}=o_{i}-\texttt{obs}_{i-i\mod f}$.

To get a state $s_{i}$, we retrieve the observation indices
$\texttt{obs\\_inds}_{i}$ and concatenate the individual frames referenced
either directly from obs or from decoding compressed frames stored in
sparse_obs (Equation 1).

###  IV-D Theoretical Compression Factor

With these definitions in mind, we can evaluate the theoretical compression
factor for our approach, that is the reduction ratio in overall memory
consumption.

Let obs store uncompressed image observations $o_{i}$ and be of shape
$(d,I,I)$, where $d=|\mathcal{D}|/f$ is the size of the uncompressed replay
buffer divided by the frame stack size $f$, and $(I,I)$ is the size of a
single grayscale input image.

Let sparse_obs store compressed image observations $S_{i}$ in SparseArray
format and be of shape $(d,f-1)$. Note that each element of this array
contains the overhead to store a pointer to a SparseArray object, which can
then be accessed to perform observation decompression (Equation 1).

Finally, let obs_inds store frame stack indices (as explained in the previous
sub-section) and be of shape $(|\mathcal{D}|,f)$.

Assuming each $S_{i}$ incurs an overhead of an 8-byte pointer and $N$ is the
number of total nonzero pixel values for all $S_{i}\in\texttt{sparse\\_obs}$,
we can define the memory size of ObsComp as the following (in bytes):

$I^{2}d+8d(f-1)+4N+4|\mathcal{D}|f.$ (3)

In comparison, to store all observations without any form of compression would
take $|\mathcal{D}|\times I\times I\times f$ bytes. Overall, our method yields
a theoretical compression factor of:

$\frac{I^{2}|\mathcal{D}|f}{I^{2}d+8d(f-1)+4N+4|\mathcal{D}|f}.$ (4)

Letting $N=d(f-1)n$, where $n$ is the average number of pixels that need to be
stored per compressed image, and $I=84$ (as used in our experiments in Section
V), we can further simplify the compression factor in terms of $f$ and $n$:

$\frac{1764f}{(1762-n)/f+2+n}.$ (5)

Noting that $n=I^{2}\phi=7056\phi$, where $\phi$ is the average _percentage_
of pixels that need to be stored per compressed image, we plot the theoretical
compression factor resulting from our approach for varying values of $f$ and
$\phi$ in Figure 2. We observe that even with a relatively small frame stack
length of $f=4$, an order of magnitude compression factor can still be
achieved as long as at most an average of 5% of the pixels remain in the
encoded images. Similarly, with a larger frame stack length of $f=10$, we can
also still achieve an order of magnitude compression factor, even with up to
25% of the pixels remaining.

Figure 2: Theoretical compression factor for different values of frame stack
length $f$ and average percent of total pixels remaining per compressed image
$\phi$. Higher values correspond to more compression achieved.

##  V EXPERIMENTS AND RESULTS

We evaluate the effectiveness of our differential encoder-based compression
framework using tasks from both the Atari 2600 benchmark [12] as well as the
DeepMind Control Suite (DMC) [9]. We analyze our results in terms of training
speed, memory consumption (and associated compression factor), and convergence
performance. Our full implementation, including values for all hyperparameters
used in our experiments, can be found in our open-source GitHub repository at:
github.com/A2R-Lab/DiffCompressDRL.

###  V-A Methodology

For all Atari benchmark tasks, we set $f$ to the default value of $4$ and
train using two popular on- and off-policy reinforcement learning algorithms:
Proximal Policy Optimization (PPO) [35] and Quantile Regression DQN (QR-DQN)
[15]. We use publicly available implementations of both PPO and QR-DQN
(credit: stable-baselines3 [43]), modifying only the buffer logic in order to
accommodate observation compression. For the DeepMind Control Suite tasks, we
set $f$ to the default value of $3$ and train on a state-of-the-art off-
policy, data augmentation-based reinforcement learning algorithm: Data-
regularized Q-v2 (DrQ-v2) [8]. We use Meta Research’s official implementation
of DrQ-v2, again only modifying logic pertaining to the storage and retrieval
of image observations. Unless otherwise noted, we run all experiments using
the default PPO, QR-DQN, and DrQ-v2 hyperparameters. Exact values for all
hyperparameters can be found on our GitHub repository.222Training was done
using a high-performance workstation with a $3.2$GHz 16-core Intel i9-12900K
and a $2.2$GHz NVIDIA GeForce RTX 4090 GPU running Ubuntu 22.04 and CUDA 12.1.

| Comp
---
Type
Env |  | FPS
---
(base)
| FPS
---
(ours)
| GB
---
(base)
| GB
---
(ours)
| Total
---
Comp
| Walker | 119 | 142 | 63.5 | 7.1 | 8.9$\times$
Half | Quadruped | 184 | 243 | 63.6 | 7.1 | 9.0$\times$
| AVG | 152 | 193 | 63.6 | 7.1 | 9.0$\times$
| Walker | 119 | 73 | 63.5 | 3.8 | 16.7$\times$
Full | Quadruped | 184 | 135 | 63.6 | 4.1 | 15.5$\times$
| AVG | 152 | 104 | 63.6 | 4.0 | 15.9$\times$
TABLE I: DrQ-v2 training speed (FPS), replay buffer size (GB), and total
memory reduction due to compression (Total Comp) of both speed-optimized
(Half) and memory-optimized (Full) compression for the DMC [9] tasks of Walker
Walk and Quadruped Walk.

###  V-B Robotic Control Tasks

In this section, we evaluate our approach on the Walker Walk and Quadruped
Walk tasks from the DMC robotic control suite [9].333We make two minor
modifications to the default environments: switching from RGB to grayscale
image observations and removing the default checkerboard floor pattern. We
define two separate levels of image-based compression: Half and Full. Half, or
speed-optimized compression uses a fully vectorized implementation of
observation indexing (Section IV-B) but does not leverage observation
compression (Section IV-C), while Full uses the entire compression stack
detailed in Section IV, but as a result, does not support vectorization.444We
note that this can be added through future work and as such speed results from
Full compression represent a conservative underestimate of future fully-
optimized performance. In addition, both implementations store all experience
values in our replay buffer directly in memory, unlike DrQ-v2, which by
default needs to store such values to disk.

Table I compares the training speed, replay buffer size, and total compression
factor (Equation 4) achieved with Half and Full compression using DrQ-v2
across the Walker Walk and Quadruped Walk DMC tasks.555We set $N$ to be the
maximum number of nonzero pixels stored at any time during training, across
all trials, to provide a conservative bound. We find that on average, Half and
Full compression reduce the replay buffer’s memory footprint by 9$\times$ and
15.9$\times$ respectively, with Full achieving a 16.7$\times$ compression
factor for Walker Walk.

In terms of training speed, our Half compression approach, which is fully
vectorized, improves latency over the baseline DrQ-v2 implementation by as
much as 32% while still achieving the aforementioned 9$\times$ compression
factor. As noted earlier, our Full compression implementation does not yet
leverage vectorization and so results in an average slow-down of 32% as
compared to the baseline DrQ-v2 implementation, which, despite storing its
replay buffer on disk, leverages multi-threaded data loaders and batch pre-
fetching in order to be as low latency as possible.

Figure 3 plots the learning curves of DrQ-v2 across the Walker Walk and
Quadruped Walk DMC environments using no compression (black), Half compression
(light blue), and Full compression (dark blue). We average our results across
five random seeds and display the standard deviation with the accompanying
shaded region. As expected, given that our compression is lossless, we find
convergence is not affected by our approach, as all final rewards both with
and without differential encoder-based compression are well within a standard
deviation of each other.

Figure 3: Learning curves of DrQ-v2 with Half (light blue), Full (dark blue),
and no (black) compression for the DMC [9] tasks of Walker Walk and Quadruped
Walk.

###  V-C Atari Environments

Table II compares the training speed and replay buffer size with and without
differential encoding-based image compression across 10 Atari 2600 benchmark
environments [12], as well as the compression factor, for both PPO and QR-DQN
using the Full compression approach.666We do not run the Atari experiments
with Half compression, as due to their low memory cost relative to the DMC
tasks, the Atari baselines are already fully vectorized in RAM. We find that
on average, Full compression is able to reduce the replay buffer’s memory
footprint by 8.9$\times$ and 9.9$\times$ for PPO and QR-DQN respectively.
While the end compression factor does vary by environment–as the number $N$ of
nonzero pixel values stored in the observation compressor (Equation 3) differs
depending on the level of inter-frame dissimilarity within the environment–we
find that 5/10 PPO and 6/10 QR-DQN environments achieve over a 10$\times$
compression factor, with our best-case environments achieving as much as a
14.2$\times$ reduction. As previously noted, Full compression is not fully
vectorized and can be optimized through future work, as such we find an
average 5% and 29% slow-down for PPO and QR-DQN respectively as compared to
the default, vectorized, stable-baselines3 implementation.

Figure 4 plots the learning curves of PPO and QR-DQN across the ten Atari
environments for the first 10M steps of training. These figures report the
average of five random seeds and also display the standard deviation with the
accompanying shaded regions. As in the DMC environments, we find that
convergence is not affected, with all final rewards being well within a
standard deviation of one another.

Alg | Env |  | FPS
---
(base)
| FPS
---
(ours)
| MB
---
(base)
| MB
---
(ours)
| Total
---
Comp
| Asteroids | 1605 | 1535 | 27.6 | 2.6 | 10.8$\times$
| BeamRider | 1339 | 1267 | 27.6 | 7.6 | 3.6$\times$
| Breakout | 1437 | 1383 | 27.6 | 2.0 | 14.2$\times$
| Enduro | 1117 | 1062 | 27.6 | 3.5 | 7.9$\times$
| MsPacman | 1432 | 1389 | 27.6 | 2.3 | 12.1$\times$
PPO | Pong | 1582 | 1512 | 27.6 | 2.0 | 13.9$\times$
| Qbert | 1507 | 1452 | 27.6 | 3.0 | 9.1$\times$
| RoadRunner | 1347 | 1267 | 27.6 | 3.0 | 9.0$\times$
| Seaquest | 1466 | 1371 | 27.6 | 2.9 | 9.6$\times$
| SpaceInvaders | 1468 | 1395 | 27.6 | 2.4 | 11.3$\times$
| AVG | 1430 | 1363 | 27.6 | 3.1 | 8.9$\times$
| Asteroids | 897 | 616 | 2694 | 265 | 10.2$\times$
| BeamRider | 792 | 547 | 2694 | 552 | 4.9$\times$
| Breakout | 833 | 587 | 2694 | 190 | 14.2$\times$
| Enduro | 722 | 536 | 2694 | 306 | 8.8$\times$
| MsPacman | 839 | 597 | 2694 | 219 | 12.3$\times$
DQN | Pong | 881 | 627 | 2694 | 192 | 14.0$\times$
| Qbert | 849 | 608 | 2694 | 213 | 12.7$\times$
| RoadRunner | 808 | 576 | 2694 | 290 | 9.3$\times$
| Seaquest | 732 | 484 | 2694 | 278 | 9.7$\times$
| SpaceInvaders | 861 | 607 | 2694 | 227 | 11.9$\times$
| AVG | 821 | 579 | 2694 | 273 | 9.9$\times$
TABLE II: Training speed in frames/steps per second (FPS), replay buffer size
(MB), and total memory reduction due to compression (Total Comp) of both PPO
and QR-DQN (DQN) across ten Atari 2600 [12] environments.

##  VI CONCLUSION AND FUTURE WORK

Figure 4: Learning curves of PPO (top) and QR-DQN (bottom) with (blue) and
without (black) pixel-based compression across ten Atari 2600 [12]
environments. Rewards are averaged across five trials and standard deviations
are shaded.

This paper presents a novel, differential encoding-based method for
observation compression, reducing the overall memory requirements of
perceptive DRL without impacting training performance.

We evaluate the compression factor, training speed, and learning performance
across ten Atari and two DMC tasks using three state-of-the-art on- and off-
policy perceptive DRL algorithms. We find that differential image encoding
reduces the memory footprint by as much as 14.2$\times$ and 16.7$\times$ for
the Atari and DMC tasks respectively. These savings also enable large-scale
DRL that previously required paging between flash and RAM to be run entirely
in RAM, improving latency for the DMC tasks by as much as 27%.

Admittedly, not all learning-based approaches will benefit equally from
observation compression. For instance, newer model-based techniques may trade
off large replay buffers for more expressive world models [6]. However, if we
are to realize lifelong, practical learning on the edge, curbing memory usage,
wherever it may be, is essential.

In future work, we hope to further optimize our method through enhanced
vectorization and parallelization in order to speed up training. Finally, we
hope to deploy our compression technique onto physical robot hardware and test
it in the context of both real-world edge RL and tiny robot learning [44]. We
hope that this effort will aid in reducing cost and compute barriers for
state-of-the-art RL across all robot platforms.

## References

  * [1] OpenAI, M. Andrychowicz, B. Baker, M. Chociej, R. Jozefowicz, B. McGrew, J. Pachocki, A. Petron, M. Plappert, G. Powell, A. Ray, J. Schneider, S. Sidor, J. Tobin, P. Welinder, L. Weng, and W. Zaremba, “Learning dexterous in-hand manipulation,” 2018. [Online]. Available: https://arxiv.org/abs/1808.00177
  * [2] T. Miki, J. Lee, J. Hwangbo, L. Wellhausen, V. Koltun, and M. Hutter, “Learning robust perceptive locomotion for quadrupedal robots in the wild,” _Science Robotics_ , vol. 7, no. 62, jan 2022. [Online]. Available: https://doi.org/10.1126%2Fscirobotics.abk2822
  * [3] Y. Song, M. Steinweg, E. Kaufmann, and D. Scaramuzza, “Autonomous drone racing with deep reinforcement learning,” 2021. [Online]. Available: https://arxiv.org/abs/2103.08624
  * [4] D. Silver, T. Hubert, J. Schrittwieser, I. Antonoglou, M. Lai, A. Guez, M. Lanctot, L. Sifre, D. Kumaran, T. Graepel, T. Lillicrap, K. Simonyan, and D. Hassabis, “Mastering chess and shogi by self-play with a general reinforcement learning algorithm,” 2017.
  * [5] OpenAI, :, C. Berner, G. Brockman, B. Chan, V. Cheung, P. Dębiak, C. Dennison, D. Farhi, Q. Fischer, S. Hashme, C. Hesse, R. Józefowicz, S. Gray, C. Olsson, J. Pachocki, M. Petrov, H. P. d. O. Pinto, J. Raiman, T. Salimans, J. Schlatter, J. Schneider, S. Sidor, I. Sutskever, J. Tang, F. Wolski, and S. Zhang, “Dota 2 with large scale deep reinforcement learning,” 2019.
  * [6] D. Hafner, J. Pasukonis, J. Ba, and T. Lillicrap, “Mastering diverse domains through world models,” 2023.
  * [7] L. Grossman and B. Plancher, “Just round: Quantized observation spaces enable memory efficient learning of dynamic locomotion,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , London, UK, May. 2023\.
  * [8] D. Yarats, R. Fergus, A. Lazaric, and L. Pinto, “Mastering visual continuous control: Improved data-augmented reinforcement learning,” 2021.
  * [9] Y. Tassa, Y. Doron, A. Muldal, T. Erez, Y. Li, D. de Las Casas, D. Budden, A. Abdolmaleki, J. Merel, A. Lefrancq, T. Lillicrap, and M. Riedmiller, “Deepmind control suite,” 2018.
  * [10] P. R. Gankidi and J. Thangavelautham, “Fpga architecture for deep learning and its application to planetary robotics,” in _2017 IEEE Aerospace Conference_ , 2017, pp. 1–9.
  * [11] S. Thrun and T. M. Mitchell, “Lifelong robot learning,” _Robotics and Autonomous Systems_ , vol. 15, no. 1, pp. 25–46, 1995, the Biology and Technology of Intelligent Autonomous Agents. [Online]. Available: https://www.sciencedirect.com/science/article/pii/092188909500004Y
  * [12] M. G. Bellemare, Y. Naddaf, J. Veness, and M. Bowling, “The arcade learning environment: An evaluation platform for general agents,” _Journal of Artificial Intelligence Research_ , vol. 47, pp. 253–279, jun 2013. [Online]. Available: https://doi.org/10.1613%2Fjair.3912
  * [13] V. Mnih, K. Kavukcuoglu, D. Silver, A. Graves, I. Antonoglou, D. Wierstra, and M. Riedmiller, “Playing atari with deep reinforcement learning,” 2013.
  * [14] M. Hessel, J. Modayil, H. van Hasselt, T. Schaul, G. Ostrovski, W. Dabney, D. Horgan, B. Piot, M. Azar, and D. Silver, “Rainbow: Combining improvements in deep reinforcement learning,” 2017.
  * [15] W. Dabney, M. Rowland, M. G. Bellemare, and R. Munos, “Distributional reinforcement learning with quantile regression,” 2017.
  * [16] L. Kaiser, M. Babaeizadeh, P. Milos, B. Osinski, R. H. Campbell, K. Czechowski, D. Erhan, C. Finn, P. Kozakowski, S. Levine, A. Mohiuddin, R. Sepassi, G. Tucker, and H. Michalewski, “Model-based reinforcement learning for atari,” 2020.
  * [17] D. Hafner, T. Lillicrap, I. Fischer, R. Villegas, D. Ha, H. Lee, and J. Davidson, “Learning latent dynamics for planning from pixels,” 2019.
  * [18] A. Srinivas, M. Laskin, and P. Abbeel, “Curl: Contrastive unsupervised representations for reinforcement learning,” 2020.
  * [19] I. Kostrikov, D. Yarats, and R. Fergus, “Image augmentation is all you need: Regularizing deep reinforcement learning from pixels,” 2021.
  * [20] M. Laskin, K. Lee, A. Stooke, L. Pinto, P. Abbeel, and A. Srinivas, “Reinforcement learning with augmented data,” 2020.
  * [21] D. Yarats, R. Fergus, A. Lazaric, and L. Pinto, “Reinforcement learning with prototypical representations,” 2021.
  * [22] M. Laskin, K. Lee, A. Stooke, L. Pinto, P. Abbeel, and A. Srinivas, “Reinforcement learning with augmented data,” 2020.
  * [23] D. Hafner, T. Lillicrap, I. Fischer, R. Villegas, D. Ha, H. Lee, and J. Davidson, “Learning latent dynamics for planning from pixels,” 2018. [Online]. Available: https://arxiv.org/abs/1811.04551
  * [24] D. Hafner, T. Lillicrap, M. Norouzi, and J. Ba, “Mastering atari with discrete world models,” 2022.
  * [25] P. Wu, A. Escontrela, D. Hafner, K. Goldberg, and P. Abbeel, “Daydreamer: World models for physical robot learning,” 2022.
  * [26] M. Tomar, U. A. Mishra, A. Zhang, and M. E. Taylor, “Learning representations for pixel-based control: What matters and why?” 2021.
  * [27] W. Yu, D. Jain, A. Escontrela, A. Iscen, P. Xu, E. Coumans, S. Ha, J. Tan, and T. Zhang, “Visual-locomotion: Learning to walk on complex terrains with vision,” in _Proceedings of the 5th Conference on Robot Learning_ , ser. Proceedings of Machine Learning Research, A. Faust, D. Hsu, and G. Neumann, Eds., vol. 164. PMLR, 08–11 Nov 2022, pp. 1291–1302. [Online]. Available: https://proceedings.mlr.press/v164/yu22a.html
  * [28] D. Wood, “Task oriented video coding: A survey,” 2022.
  * [29] T. Ebrahimi and C. Horne, “Mpeg-4 natural video coding – an overview,” _Signal Processing: Image Communication_ , vol. 15, no. 4, pp. 365–385, 2000\. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0923596599000545
  * [30] J.-W. Chen, C.-Y. Kao, and Y.-L. Lin, “Introduction to h. 264 advanced video coding,” in _Proceedings of the 2006 Asia and South Pacific Design Automation Conference_ , 2006, pp. 736–741.
  * [31] W. Yang, H. Huang, Y. Hu, L.-Y. Duan, and J. Liu, “Video coding for machine: Compact visual representation compression for intelligent collaborative analytics,” 2021.
  * [32] Z. Huang, C. Jia, S. Wang, and S. Ma, “Visual analysis motivated rate-distortion model for image coding,” 2021.
  * [33] Z. Yang, Y. Wang, C. Xu, P. Du, C. Xu, C. Xu, and Q. Tian, “Discernible image compression,” in _Proceedings of the 28th ACM International Conference on Multimedia_. ACM, oct 2020\. [Online]. Available: https://doi.org/10.1145%2F3394171.3413968
  * [34] N. Le, H. Zhang, F. Cricri, R. Ghaznavi-Youvalari, and E. Rahtu, “Image coding for machines: an end-to-end learned approach,” in _ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, jun 2021. [Online]. Available: https://doi.org/10.1109%2Ficassp39728.2021.9414465
  * [35] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal Policy Optimization Algorithms,” _arXiv:1707.06347 [cs]_ , July 2017\.
  * [36] R. S. Sutton, D. M. A. llester, S. Singh, and Y. Mansour, “Policy Gradient Methods for Reinforcement Learning with Function Approximation,” 2000, pp. 1057–1063.
  * [37] P. J. Huber, “Robust Estimation of a Location Parameter,” _The Annals of Mathematical Statistics_ , vol. 35, no. 1, pp. 73 – 101, 1964. [Online]. Available: https://doi.org/10.1214/aoms/1177703732
  * [38] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” 2015. [Online]. Available: https://arxiv.org/abs/1509.02971
  * [39] K. Sayood, “Chapter 11 - differential encoding,” in _Introduction to Data Compression (Fifth Edition)_ , fifth edition ed., ser. The Morgan Kaufmann Series in Multimedia Information and Systems, K. Sayood, Ed. Morgan Kaufmann, 2018, pp. 351–378. [Online]. Available: https://www.sciencedirect.com/science/article/pii/B9780128094747000112
  * [40] N. B. Ruparelia, “The history of version control,” _ACM SIGSOFT Software Engineering Notes_ , vol. 35, no. 1, pp. 5–9, 2010.
  * [41] L.-J. Lin, “Self-improving reactive agents based on reinforcement learning, planning and teaching,” _Machine Learning_ , vol. 8, no. 3, pp. 293–321, 1992. [Online]. Available: https://doi.org/10.1007/BF00992699
  * [42] T. Schaul, J. Quan, I. Antonoglou, and D. Silver, “Prioritized experience replay,” 2016.
  * [43] A. Raffin, A. Hill, A. Gleave, A. Kanervisto, M. Ernestus, and N. Dormann, “Stable-baselines3: Reliable reinforcement learning implementations,” _Journal of Machine Learning Research_ , vol. 22, no. 268, pp. 1–8, 2021\. [Online]. Available: http://jmlr.org/papers/v22/20-1364.html
  * [44] S. M. Neuman, B. Plancher, B. P. Duisterhof, S. Krishnan, C. Banbury, M. Mazumder, S. Prakash, J. Jabbour, A. Faust, G. C. de Croon, and V. Janapa Reddi, “Tiny robot learning: Challenges and directions for machine learning in resource-constrained robots,” in _2022 IEEE International Conference on Artificial Intelligence Circuits and Systems (AICAS)_ , Incheon, Korea, June 2022.


# OpenBTE: a Solver for ab-initio Phonon Transport in Multidimensional
Structures

Giuseppe Romano  Massachusetts Institute of Technology, 77 Massachusetts
Avenue, Cambridge, Massachusetts 02139, USA

###### Abstract

Controlling heat flow at the nanoscales is pivotal to several applications,
including thermal energy harvesting and heat management. However, engineering
nanostructures is challenging because phonon-boundary interaction, not
contemplated by Fourier’s law, must be taken into account. Nondiffusive
models, such as the Boltzmann transport equation (BTE), have been successfully
employed to capture size effects in complex structures; however, their
widespread has been hindered by the limited offer of open-source solvers in
this space. We fill this void by introducing `OpenBTE`, an efficient solver
for the steady-state phonon BTE in multidimensional structures. This tool is
interfaced to first-principles calculations, thus it unlocks the calculations
of thermal-related properties with no fitting-parameters. As an example, we
employ `OpenBTE` to compute the temperature and flux maps, as well as the
mode-resolved, effective thermal conductivity of Si membranes with infinite
and finite thickness. By unlocking fast nanoscale heat transport simulations,
`OpenBTE` may help accelerate the design of nanomaterials for thermal energy
applications.

###### keywords:

Thermal transport, nanostructures

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: `OpenBTE`
Developer’s repository link: https://github.com/romanodev/openbte.git
Licensing provisions: MIT
Programming language: `Python`
Nature of problem: Ab-initio calculation of thermal transport properties,
including temperature and heat flux maps, as well as the effective thermal
conductivity of multidimensional nanostructures.
Solution method: We implement the space-dependent, ab-initio Boltzmann
transport equation (BTE) for phonons. The first-principles calculation of
phonon group velocities, frequencies and scattering times, is delegated to
external packages to which `OpenBTE` is interfaced. The current implementation
focuses on a flavor of the BTE that interpolates phonon populations onto their
vectorial mean-free-paths, allowing for fast and accurate thermal transport
simulations.
Additional comments: Vectorization via `NUMPY` and parallelization (including
with shared memory) with `MPI4Py`.

##  1 Introduction

The increasing accuracy of first-principles calculations of heat transport,
along with the growing availability of computing resources, has enabled the
systematic prediction of lattice thermal conductivity in semiconductors and
two-dimensional (2D) materials [1, 2, 3]. If key theoretical developments have
solidified our understanding of phonon dynamics at the nanoscales, the fast
progress of this field has been unarguably aided by the release of open-source
packages. For example, `ShengBTE` [4] implements the iterative solution of the
Boltzmann transport equation (BTE) [5], with force constants computed by
first-principles [6, 7]. The tool can be primarily used for either bulk
materials or nanowires [8]. Later iterations include `AlmaBTE` [9], which,
among several improvements, features a Monte Carlo solver for superlattices,
and `FourPhonons` [10, 11], a tool that, as the name suggests, can handle
four-phonon scattering. Another notable package is `Phono3Py` [12], where the
BTE is solved both within the relaxation time approximation (RTA) and
including the full scattering operator [13]. Albeit not strictly required,
these tools are mostly used in tandem with super-cell approaches. When a
density functional perturbation theory (DFPT) is available for a given
material, it is often convenient to compute directly the Fourier-transformed
force constants, without resorting to supercells; a prominent package in this
space is `d3q` [14], implemented on top of `QUANTUM ESPRESSO` [15]. This
software also includes the possibility of simulating thin films [16] and
coherence effects [17]. In a similar effort, the recently released `Kaldo`[18]
combines the Green-Kubo model and the BTE to take into account material
disorder [19].

While the dissemination of these software manifests the increasing excitement
around phonon related applications, their focus is primarily on bulk materials
or simple geometries. In fact, most phonon dynamics simulations in complex
geometries are currently performed by MonteCarlo solvers [20, 21], with a
recent implementation provided by `MCBTE` [22]. Thus, an open-source tool that
solves space-dependent heat transport deterministically, i.e. in the same
spirit as in bulk calculations, is still pending. We fill this gap by
introducing `OpenBTE`, a solver for the first-principles, multidimensional
phonon BTE. Key capabilities include the calculation of the temperature and
heat flux maps, as well as the mode-resolve effective thermal conductivity.
Ab-initio harmonic and anharmonic properties are delegated to external
software. While flexible and straightforward to expand, `OpenBTE` currently
implements the anisotropic-mean-free-path BTE (aMFP-BTE) [23], a method based
on the interpolation of the phonon populations onto the vectorial MFP space.
This method has shown a 50x speed up with respect to the mode-resolve BTE for
thick Si membranes, while not compromising on the accuracy. Furthermore, it
scales with constant time with respect to the number of phonon branches,
therefore enabling fast simulations of size effects in complex-unit-cell
materials, such as bismuth telluride.

The paper is organized as follows. First, we review the model and the
assumptions therein, then we provide details on the implemented algorithm and
the software’s main modules. An example of a porous Si membrane will follow,
and final remarks conclude the paper. We also provides technical
implementation details in the appendices. Enabling fast simulations of heat
transport in nanostructures, `OpenBTE` sets out to accelerate the development
of thermal-related applications, such as thermal management and
thermoelectrics. `OpenBTE` is released under the `GPL-3.0` license and
currently hosted on GitHub [24].

##  2 Model

The steady-state, linearized phonon BTE in the temperature formulation reads
[23]

$-C_{\mu}\mathbf{v}_{\mu}\cdot\nabla\Delta
T_{\mu}(\mathbf{r})=\sum_{\nu}W_{\mu\nu}\Delta T_{\nu}(\mathbf{r}),$ (1)

where $\mu,\nu$ collectively indicate phonon branch and wave vector, running
up to $N_{p}$ and $N_{q}$, respectively. The unknowns $\Delta
T_{\mu}(\mathbf{r})$ are the phonon pseudo-temperatures (referred to hereafter
simply as temperatures), connected to the deviational non-equilibrium phonon
population via $\Delta T_{\mu}(\mathbf{r})=\Delta
n_{\mu}(\mathbf{r})\hbar\omega_{\mu}C_{\mu}^{-1}$. The mode-specific heat
capacity is $C_{\mu}=k_{B}\left(\eta_{\mu}/\sinh{\eta_{\mu}}\right)^{2}$, with
$\eta_{\mu}=\hbar\omega_{\mu}/k_{B}/T/2$. The terms $\hbar\omega_{\mu}$ and
$\mathbf{v}_{\mu}$ are the phonon frequencies and the group velocities,
respectively. In this formulation, the thermal flux is given by
$\mathbf{J}(\mathbf{r})=\left(N_{q}\mathcal{V}\right)^{-1}\sum_{\mu}\Delta
T_{\mu}(\mathbf{r})\mathbf{S}_{\mu}$, where $\mathcal{V}$ is the volume of the
unit cell, and $\mathbf{S}_{\mu}=C_{\mu}\mathbf{v}_{\mu}$. The term
$W_{\mu\nu}$ is the energy form of the scattering operator [25]. While we plan
to make available the implementation of Eq. 1, currently only the relaxation-
time-approximation (RTA) is supported; within the RTA, the RHS of Eq. 1
becomes $C_{\mu}/\tau_{\mu}\left[\Delta T_{\mu}(\mathbf{r})-\Delta
T(\mathbf{r})\right]$, leading to

$-\mathbf{v}_{\mu}\cdot\nabla\Delta T_{\mu}(\mathbf{r})=\frac{\Delta
T_{\mu}(\mathbf{r})-\Delta T(\mathbf{r})}{\tau_{\mu}},$ (2)

where $\tau_{\mu}$ is the relaxation time, and $\Delta T(\mathbf{r})$ is the
pseudo lattice temperature; this term is evaluated so that energy conservation
is satisfied, i.e. $\nabla\cdot\mathbf{J}(\mathbf{r})=0$, yielding $\Delta
T(\mathbf{r})=\left[\sum_{\nu}C_{\nu}/\tau_{\nu}\right]^{-1}\sum_{\mu}C_{\mu}/\tau_{\mu}\Delta
T_{\mu}(\mathbf{r})$ [26].

The simulation domain is a cuboid of size $L_{x}$, $L_{y}$ and $L_{z}$ with
periodic boundary conditions applied throughout by default. A difference of
temperature $\Delta T_{\mathrm{ext}}$ is also enforced along a chosen
Cartesian axis. At the walls of the pores, energy conservation leads to
$\sum_{\nu}\Delta T_{\nu}\mathbf{S}_{\mu}\cdot\mathbf{\hat{n}}=0$, which can
be split into incoming and outgoing flux, i.e. $\sum_{\mu}\Delta
T_{\mu}\mathrm{ReLu}(\mathbf{S}_{\mu}\cdot\mathbf{\hat{n}})=-\sum_{\mu}\Delta
T_{\mu}\mathrm{ReLu}(-\mathbf{S}_{\mu}\cdot\mathbf{\hat{n}})$; note that we
define $\mathrm{ReLu}(x)=x\Theta(x)$, with $\Theta(x)$ being the Heaviside
function. The boundary condition is imposed to phonons leaving the surface;
generally, we have $\Delta T_{\mu}^{-}=\sum_{\nu}R_{\mu\nu}\Delta
T_{\nu}^{+}$, where $-$ ($+$) are for outgoing (incoming) phonons, and
$R_{\mu\nu}$ is a mode-resolved bouncing matrix. However, currently only a
simplified model is available, with which we assume that phonons thermalize to
a single temperature [20]; this choice corresponds to
$R_{\mu\nu}=\left[\sum_{k}\mathrm{ReLu}(\mathbf{S}_{k}\cdot\mathbf{\hat{n}})\right]^{-1}\mathrm{ReLu}(\mathbf{S}_{\mu}\cdot\mathbf{\hat{n}})$.
Lastly, the effective thermal conductivity, assuming we have imposed the
temperature gradient along $x$, reads

$\kappa^{\mathrm{eff}}=-\frac{L_{x}}{\Delta
T_{\mathrm{ext}}}\sum_{\mu}\int_{A}\frac{d\mathbf{S}}{A}\mathbf{J}_{\mu}\cdot\mathbf{\hat{n}}.$
(3)

Equation 2, when solved iteratively, requires inverting as many linear systems
as the number of phonons modes, an endeavour that may easily become
prohibitive. To ameliorate the computational effort, we instead implement the
recently introduced anisotropic-MFP-BTE (aMFP-BTE), a model based on
interpolating the phonon temperatures onto the vectorial MFP space [23]. In
practice, instead of solving Eq. 2 $N_{p}N_{q}$ times at each iteration, we
solve it for a set of vectorial MFPs, $\mathbf{F}_{ml}$, located on a
spherical grid. The indices $m$ and $l$ label the MFP $\Lambda_{m}$ (up to
$N_{\Lambda}$) and phonon directions $\mathbf{\hat{s}}_{l}$ (up to
$N_{\Omega}$), respectively. As $N_{\Lambda}N_{\Omega}<<N_{p}N_{q}$, while not
compromising on accuracy, this scheme results in a much faster runtime with
respect to the case with mode-resolved resolution, i.e. no interpolation.

##  3 Implementation

As derived in A, the discretization of aMFP-BTE model in real and momentum
space yields the following iterative linear system

$\displaystyle\mathbf{A}_{ml}\mathbf{\Delta
T}_{ml}^{(n)}=\sum_{m^{\prime}l^{\prime}}\mathbf{S}_{ml}^{m^{\prime}l^{\prime}}\mathbf{\Delta
T}_{m^{\prime}l^{\prime}}^{(n-1)}+\mathbf{P}_{ml},$ (4)

where $\mathbf{A}_{ml}$ is the stiffness matrix for a given MFP and direction,
$\mathbf{S}_{ml}^{m^{\prime}l^{\prime}}$ is associated to the lattice
temperature and adiabatic boundary conditions, and $\mathbf{P}_{ml}$
represents the perturbation. Once Eq. 4 converges, the effective thermal
conductivity is computed as
$\kappa_{\mathrm{eff}}=\kappa_{\mathrm{bal}}/2+\sum_{ml}\mathbf{K}_{ml}^{T}\cdot\mathbf{\Delta
T}_{ml}$, where $\kappa_{\mathrm{bal}}$ is the ballistic thermal conductivity,
and $\mathbf{K}_{ml}$ is associated to the thermal flux.

The first guess to Eq. 4 is given by the standard heat conduction equation.
While it is straightforward to discretize the diffusive equation in orthogonal
grids, it becomes convoluted for unstructured meshes. In `OpenBTE`, we use the
approach described in [27], where the non-orthogonal contribution of the flux
is treated explicitly, i.e. it depends on previous solutions. The resulting
system, analogously to Eq. 4, is an iterative linear system,

$\mathbf{A}\mathbf{\Delta T}^{(n)}=\mathbf{S}\mathbf{\Delta
T}^{(n-1)}+\mathbf{P},$ (5)

where $\mathbf{S}$ arises from the non-orthogonality of the mesh. Note that
$\mathbf{S}_{ml}^{m^{\prime}l^{\prime}}$ in Eq. 4 and $\mathbf{S}$ in Eq. 5
are not related. The same applies with the terms sharing the same basename in
these two equations. The derivation of Eq. 5 is given in B.

Equation 4 entails, for each iteration $n$, the inversion of
$N_{\Lambda}N_{\Omega}$ matrices, whose size is the number of volumes in the
computational domain. To enhance the computational efficiency, three distinct
strategies have been put in place. First, by noting that the stiffness matrix
does not depend on $n$, we are able to compute LU factorizations for the first
iteration, and reuse them until convergence. Another major speed up is
achieved by vectorizing the assembly of the stiffness matrices, while also
exploiting their sparsity. Lastly, for each $n$, Eq. 4 consists of a set of
independent problems, thus it can be trivially parallelized [28]; to this end,
we assign a set of phonon directions to each core ($\Omega_{c}$), while the
indices $m$ run serially. The implementations of these three methods rely on
`SCIPY`, `NUMPY` [29] and `MPI4PY` [30], respectively.

When vectorization and parallelization are both used, we may save memory by
noting that, as shown in A, the stiffness matrices have the form
$\mathbf{A}_{ml}=\mathbf{I}+\Lambda_{m}\mathbf{A}_{l}$; thus vectorization can
be done only for $\mathbf{A}_{l}$, while $\mathbf{A}_{ml}$ is assembled on the
fly right before LU factorization. Similar arguments hold for
$\mathbf{P}_{ml}$ and $\mathbf{S}_{ml}^{m^{\prime}l^{\prime}}$. The overall
methodology is summerized in Algorithm 1, where, for readability, only
$\mathbf{A}_{ml}$ has been expanded upon.

Algorithm 1 Workflow for thermal conductivity calculations. The term
$\Omega_{c}$ is the set of solid angles assigned to core $c$.

1:$\mathbf{\Delta T}_{ml}^{(0)}\leftarrow$ Fourier’s law.

2:while error $>1e^{-3}$ do

3: $\kappa^{\mathrm{eff},(n+1)}=\kappa_{\mathrm{bal}}/2$

4: for l in $\\{\Omega_{c}\\}$ do

5: for m = 1:$N_{\Lambda}$ do

6:
$\mathbf{B}=\sum_{m^{\prime}l^{\prime}}\mathbf{S}_{ml}^{m^{\prime}l^{\prime}}\mathbf{\Delta
T}_{m^{\prime}l^{\prime}}^{(n-1)}+\mathbf{P}_{ml}$

7: $\mathbf{\Delta
T}_{ml}^{(n)}\leftarrow\left(\mathbf{I}+\Lambda_{m}\mathbf{A}_{l}\right)\mathbf{\Delta
T}_{ml}^{(n)}=\mathbf{B}$

8: $\kappa^{\mathrm{eff},(n+1)}+=\mathbf{K}_{ml}^{T}\cdot\mathbf{\Delta
T}_{ml}$

9: end for

10: end for

11: error
$\leftarrow|(\kappa^{\mathrm{eff},(n+1)}-\kappa^{\mathrm{eff},(n)})|/\kappa^{\mathrm{eff},(n+1)}$

12: $n\leftarrow n+1$

13:end while

##  4 Workflow

The easiest way to run `OpenBTE` is to use it as a module in `Python`. A
sketch for the workflow is illustrated in Fig. 1. The core module is `Solver`,
which runs both the heat conduction equation as a first guess to Eq. 4 and all
subsequent iterations, outlined in Algorithm 1. Before invoking it, though, it
is necessary to define the material’s internal boundaries using the `Geometry`
module, and bulk-related data via the `Material` module. The latter is the
interface to external packages. Once the simulation is finalized, results can
be plotted using the `Plot` module. Note that `OpenBTE` follows a pure
functional approach, where data is simply stored as dictionaries. Each module
is expanded upon in the next sections.

Figure 1: The software architecture of ``OpenBTE!. The key components are the
``Material!, ``Geometry!, ``Solver! and ``Plot! modules. Currently, the code
is interfaced to ``AlmaBTE! and ``Phono3Py!, which computes first-principles
data. Mesh generation is delegated to ``GMSH!, and output maps can be
visualized with ``Paraview!.

##  5 Material

The `Material` module converts mode-resolved bulk data, $C_{\mu}$,
$\tau_{\mu}$ and $\mathbf{v}_{\mu}$, into the spherically-resolved quantities
$\mathbf{G}_{ml}$ and $C_{ml}$. The input files, `rta.npz`, are created by
postprocessing data obtained by external packages. Currently, `OpenBTE`
provides interfaces to `AlmbaBTE` [9] and `Phono3Py` [12]. A small set of
precomputed `rta.npz` is also provided. The file created by `Material` is
stored into `material.npz` and ready to be used by `Solver`. Depending on the
geometry, this module has two models:

  * 1.

`rta3D`: this model must be used for geometries with $l_{z}>0$. In this case,
the angular interpolation of $\mathbf{F}_{\mu}$ is performed in both polar and
azimuthal angles.

  * 2.

`rta2DSym`: when $l_{z}$ is not specified, the actual simulation domain is 2D,
and this model must be used. Note that when the base material is 3D, using
this option is equivalent to considering an infinite thickness. In fact, in
this case, the mode-resolved $\mathbf{F}_{\mu}$ are mapped into
$\mathbf{F}_{tk}=\Lambda_{t}\cos(\phi_{k})$, where $\Lambda_{t}$ implicitly
includes the MFP and the azimuthal angle, and $\phi_{k}$ is the polar angle on
the x-y plane [23].

The model as well as the discretization grid can be chosen with

⬇

from openbte import Material

\parMaterial(filename=’rta_Si_300’,model=’rta2DSym’,n_phi=48)

where, by default, the material file is taken from the database, whose entries
are updated online. In this case, we choose to discretize the polar angle into
48 slices.

##  6 Geometry module

The `Geometry` module allows the user to build an arbitrarily patterned
nanomaterial. It interfaces with `GMSH` [31] to create an unstructured Delaney
mesh, and creates the file `geometry.npz`. While we plan to handle user-
provided meshes, `OpenBTE` currently supports 2D membranes with finite and
infinite thickness. The size of the unit cell is specified by `lx` and `ly`,
while the thickness is defined by `lz`. The direction of the applied gradient
is indicated with `direction` (default is `x`). Periodic boundary conditions
are applied along all directions; however, adiabatic surfaces can be imposed
along specific directions with `Periodic=[True,False,True]`, where each
component declares whether a surface is periodic. Currently, there are two
geometry models, as expanded upon in the next sections.

###  6.1 Lattice model

This model aids the creation of 2D periodic porous materials with pores having
either predefined shapes or custom ones. In analogy to crystals, a porous
material is defined as a set of unit vectors, which in our case are
$\mathbf{a}_{1}=l_{x}\mathbf{\hat{x}}$ and
$\mathbf{a}_{2}=l_{y}\mathbf{\hat{y}}$, and a base. The latter is a list of
x-y pairs, each of them being one pore. Here is an example with a single pore
in the origin

⬇

from openbte import Geometry

\parGeometry(model=’lattice’,base=[[0,0]],shape=’circle’)

In the example above, we chose a simple shape, a circle; however, it is
possible also to supply a user-defined shape using `shape=custom`. In such a
case, the keyword `shape_functions` must be provided, optionally along with
`shape_options`. For example, as shown in Fig. 2a-b and c, we can define a
structure having two pores, each of them described by the same custom shape.

###  6.2 Custom model

Using the option `model=custom`, it is possible to define an arbitrary
patterning by providing a list of polygons, with the periodicity being
automatically guaranteed. This option significantly widens up the number of
possible geometries at the expense of more input by the user. In spirit of
subtractive manufacturing, each polygon can be seen as a region where the
material is carved out. An example of custom model is depicted in Fig.2-d.

Figure 2: a) An example of structure using the keyword-value pair
``shape=custom! for the ``lattice! model of the ``Geometry! module. b) the
unit cell comprising two pores at designed positions. Note that when a pore
crosses a boundary it gets repeated. c) singe-pore view of material. The
coordinates of the shape are defined via the ``shape_function! option. d) An
example of ``custom! shape obtained by carving out the base material with
user-defined polygons.

##  7 Plot Module

Simulations results can be plotted via the `Plot` module. Currently, the
following features are implemented:

  * 1.

`vtu`: the temperature and flux maps are stored in a `vtu` file and ready to
be opened by `Paraview` [32].

  * 2.

`maps`: the maps are shown as web-ready, interactive structures.

  * 3.

`kappa_mode`: this model computes the mode-resolved thermal conductivity from
Eq 21, and stores in a file.

  * 4.

`line`: this option, which currently is implemented only for 2D systems,
allows for line plots, using linear interpolation.

##  8 Thermal transport in porous Si

(a)

(b)

(c)

(d)

Figure 3: a) Thermal transport reduction when using Fourier’s law, for regular
circular pores and for a user-define shape; the Eucken-Maxwell model is also
plotted. In the inset, the magnitude of the thermal map for the user-defined
structure. The porosity is 0.1 and $L$ = 50 nm. b) a cut of the temperature
maps, computed from both the BTE and Fourier’s law, along the $x$-axis. c) The
effective thermal conductivity, $\kappa^{\mathrm{eff}}$, of the user-defined
structure with thickness $t$ = 10 nm, and for different porosities. In the
inset, the case with porosities = 0.1 and 0.17 are shown. d) The mode-resolved
thermal conductivity, $\kappa^{\mathrm{eff}}_{\mu}$, versus the MFP,
$|F_{\mu,x}|$, for the cases $t$ = 10 nm and $t=\infty$.

In some cases, it is convenient to use Fourier’s law to estimate macroscopic
geometric effects. In this first part, we thus aim to assess the reliability
of the diffusive solver, implemented in `OpenBTE`. To this end, we consider a
membrane with circular pores and infinite thickness, for which
$\kappa_{\mathrm{fourier}}/\kappa_{\mathrm{bulk}}=r$ can be evaluated
analytically. In fact, in this case, the thermal transport suppression is
given by the Eucken-Maxwell model $r=(1-\phi)/(1+\phi)$ [33], $\phi$ being the
porosity. In Fig. LABEL:fig1a, we compare this model with the results from
`OpenBTE` for different porosities, obtaining excellent agreement ($<1\%$). We
note that by using the option `only_fourier=True` it is possible to run only
the Fourier’s model. Next, we compute $r$ for a more complicated structure,
comprising two concatenated sub-lattices of pores with different user-provided
shapes; this structure is generated using the `shape=custom` option from the
`lattice` model, a two-pore base, and different options for each pore. As
shown in Fig. LABEL:fig1a, we obtain a much larger heat transport suppression
due to the longest path heat must travel through the material.

Let us now compute $\kappa^{\mathrm{eff}}$ using the BTE, encoded in Eq. 4.
For the chosen porosity of $0.1$ and $L$ = 50nm, we obtain 17.4 Wm-1K-1, well
below the corresponding macroscopic value 99.4 Wm-1K-1. In Fig. LABEL:fig1b,
we plot a cut of the temperature maps, both for the BTE and Fourier cases,
along the $x$-axis. As expected, they have different trends, the BTE one being
flatter due to ballistic transport. When a finite thickness is added, phonon
undergo stronger suppression; in fact, for thickness $t$=10 nm and porosity
0.1, we obtain $\kappa^{\mathrm{eff}}=8.9$ W m-1K-1, roughly a half of the
value with $t=\infty$. In Fig. LABEL:fig1c, $\kappa^{\mathrm{eff}}$ for
different porosities is plotted. From the insets, we note that, as the
porosity becomes larger, the heat flux peaks around the phonon bottleneck,
e.g. the areas between the pores. Lastly, the mode-resolved thermal
conductivities for $t$ = 10 nm and $\infty$ are plotted in Fig. LABEL:fig1d;
analogously to Ref. [23], we can appreciate the stronger suppression for large
MFPs for the case with finite thickness.

##  9 Conclusion

We have presented `OpenBTE`, a software for computing deterministically space-
dependent heat transport in 2D and 3D structures. The efficient and modern
implementation of the underlying algorithms enable accurate simulations both
in the cloud and on common laptops, thus democratizing this type of
simulations. Future work includes going beyond the RTA of the BTE and adding
transient dynamics. GPU support and differentiability are also on the roadmap.
Furthermore, we plan to develop a hybrid Fourier/BTE solver, which has the
potential to dramatically reduce runtimes for the cases with 3D structures.
Filling an important gap in the current offer of publicly available phonon
transport solvers, `OpenBTE` sets out to guide experiments on nanoscale heat
transport and facilitate high-throughput prediction of thermal properties of
nanostructures.

##  10 Acknowledgments

Research was partially supported by the Solid-State SolarThermal Energy
Conversion Center (S3TEC), an Energy Frontier Research Center funded by the
U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences
(BES), under Award No. DESC0001. The author thanks all the early users, who
helped improving `OpenBTE`.

##  Appendix A The aMFP-BTE

`OpenBTE` employs the upwind finite-volume discretization of Eq. 2, where the
flux outgoing from volume $i$ contributes to the diagonal of the stiffness
matrix at row $i$ [27]. A sketch depicting incoming and outgoing flux is
illustrated in Fig. 4a. The discretized mode-resolved BTE reads [23]

$\sum_{j\mu^{\prime}}\left[\delta_{cc^{\prime}}\delta_{\mu\mu^{\prime}}+A_{ij\mu^{\prime}}+B_{\mu\mu^{\prime}j}\delta_{ij}+D_{\mu^{\prime}}\delta_{ij}\right]\Delta
T_{j\mu^{\prime}}=P_{i\mu},$ (6)

where

$\displaystyle A_{ij\mu}$ $\displaystyle=$
$\displaystyle\sum_{k}\mathrm{ReLu}(\mathbf{F}_{\mu}\cdot\mathbf{R}_{ik})\delta_{ij}-\mathrm{ReLu}(-\mathbf{F}_{\mu}\cdot\mathbf{R}_{ik})\delta_{jk},$
(7) $\displaystyle B_{\mu\mu^{\prime}j}$ $\displaystyle=$
$\displaystyle-\sum_{s}\mathrm{ReLu}(-\mathbf{F}_{\mu}\cdot\mathbf{\hat{n}}_{s})g_{si}\frac{\mathrm{ReLu}(\mathbf{S}_{\nu}\cdot\mathbf{\hat{n}}_{s})}{\sum_{k}\mathrm{ReLu}(\mathbf{S}_{k}\cdot\mathbf{\hat{n}}_{s})},$
(8)

arise flux discretization and adiabatic boundary conditions, respectively. The
term $D_{\mu}=C_{\mu}/\tau_{\mu}\left[\sum_{k}C_{k}/\tau_{k}\right]^{-1}$
contributes to the lattice pseudo-temperature; lastly, the perturbation is

$P_{i\mu}=-\sum_{j}\mathrm{ReLu}(-\mathbf{F}_{\mu}\cdot\mathbf{R}_{ij}^{P})\Delta
T_{\mathrm{ext}}.$ (9)

The terms $\mathbf{R}_{ij}$, $\mathbf{R}_{ij}^{P}$ and $g_{sc}$ define the
mesh, and are given by

$\displaystyle\mathbf{R}_{ij}$ $\displaystyle=$
$\displaystyle\begin{cases}\frac{1}{V_{c}}A_{ij}\mathbf{\hat{n}}_{ij},&\text{if
}i\text{ and }j\text{ are neighbors}\\\ 0,&\text{otherwise},\end{cases}$ (10)
$\displaystyle\mathbf{R}_{ij}^{P}$ $\displaystyle=$
$\displaystyle\begin{cases}\frac{1}{V_{i}}A_{ij}\mathbf{\hat{n}}_{ij},&\text{if
}i\text{ and }j\text{ are periodic to each other}\\\
0,&\text{otherwise}.\end{cases}$ (11) $\displaystyle g_{si}$ $\displaystyle=$
$\displaystyle\begin{cases}\frac{1}{V_{i}}A_{s},&\text{if }s\text{ is a side
of }i\\\ 0,&\text{otherwise},\end{cases}$ (12)

where $A_{ij}$ is the area of the surface between the volumes $i$ and $j$, and
$\mathbf{\hat{n}}_{ij}$ is its normal, pointing toward the volume $j$. The
term $V_{i}$ is the volume of the element $i$. We solve 6 iteratively,

$\sum_{j}\left(\delta_{ij}+A_{ij\mu}\right)\Delta
T_{j\mu}^{(n)}=P_{i\mu}+\sum_{\mu^{\prime}}\left(B_{\mu\mu^{\prime}i}+D_{\mu^{\prime}}\right)\Delta
T_{i\mu^{\prime}}^{(n-1)}.$ (13)

Equation 13 requires solving as many linear systems as the number of phonons
modes. To overcome this issue, `OpenBTE` implements the aMFP-BTE [23], where
the phonon temperatures are interpolated over the vectorial MFPs,
$\mathbf{F}_{\mu}=\sum_{ml}c_{ml}^{\mu}\mathbf{F}_{ml}$; the chosen
$\mathbf{F}_{ml}=\Lambda_{m}\mathbf{\hat{s}}_{l}$ are located uniformnly on a
spherical grid. The labels $m$ and $l$ refer to MFP and solid angle,
respectively. Analogously, the temperatures are $\mathbf{\Delta
T}_{\mu}=\sum_{ml}c_{ml}^{\mu}\mathbf{\Delta T}_{ml}$. In the aMFP-BTE
formulation, Eq. 13 becomes

$\sum_{j}\left(\delta_{ij}+\Lambda_{m}A_{ijl}\right)\Delta
T_{jml}^{(n)}=\Lambda_{m}P_{il}+\sum_{m^{\prime}l^{\prime}}\left(\Lambda_{m}B_{lm^{\prime}l^{\prime}i}+D_{m^{\prime}l^{\prime}}\right)\Delta
T_{im^{\prime}l^{\prime}}^{(n-1)},$ (14)

where

$\displaystyle A_{ij^{\prime}l}$ $\displaystyle=$
$\displaystyle\sum_{k}\mathrm{ReLu}(\mathbf{\hat{s}}_{l}\cdot\mathbf{R}_{ik})\delta_{ij}-\mathrm{ReLu}(-\mathbf{\hat{s}}_{l}\cdot\mathbf{R}_{ik})\delta_{jk},$
(15) $\displaystyle B_{lm^{\prime}l^{\prime}j}$ $\displaystyle=$
$\displaystyle-\sum_{s}\mathrm{ReLu}(-\mathbf{\hat{s}}_{l}\cdot\mathbf{\hat{n}}_{s})g_{si}\frac{\mathrm{ReLu}(\mathbf{S}_{m^{\prime}l^{\prime}}\cdot\mathbf{\hat{n}}_{s})}{\sum_{m^{\prime\prime}l^{\prime\prime}}\mathrm{ReLu}(\mathbf{S}_{m^{\prime\prime}l^{\prime\prime}}\cdot\mathbf{\hat{n}}_{s})}.$
(16)

The term $P_{il}$ is the perturbation, given by

$P_{il}=-\sum_{j}\mathrm{ReLu}(-\mathbf{\hat{s}}_{l}\cdot\mathbf{R}_{ij}^{P})\Delta
T_{\mathrm{ext}}.$ (17)

The quantities $D_{ml}=\sum_{\mu}c_{ml}^{\mu}D_{\mu}$ and
$\mathbf{S}_{ml}=\sum_{\mu}c_{ml}^{\mu}\mathbf{S}_{\mu}$ contributes to the
lattice temperature and thermal flux, respectively. Finally, after defining
$S_{iml}^{m^{\prime}l^{\prime}}=\Lambda_{m}B_{lm^{\prime}l^{\prime}i}+D_{m^{\prime}l^{\prime}}$
in Eq. 14, we obtain Eq. 4. To compute the effective thermal conductivity from
Eq. 3, we need to integrate over the hot contact and the momentum space.
According to the upwind scheme, the volume from which we pick the temperature
depends on whether the flux is outgoing from or incoming to such volume;
accordingly, Eq. 3 translates into

$\displaystyle\kappa_{\mathrm{eff}}=$ $\displaystyle-$
$\displaystyle\frac{L}{A\Delta T_{\mathrm{ext}}}\sum_{ijml}\big{[}\Delta
T_{cml}\mathrm{ReLu}(\mathbf{S}_{ml}\cdot\mathbf{R}^{K}_{ij})-$ (18)
$\displaystyle-$ $\displaystyle\left(\Delta T_{iml}+\Delta
T_{\mathrm{ext}}\right)\mathrm{ReLu}(-\mathbf{S}_{ml}\cdot\mathbf{R}^{K}_{ij})\big{]}=$
$\displaystyle=$
$\displaystyle\frac{\kappa_{\mathrm{bal}}}{2}+\sum_{ml}\mathbf{K}_{ml}^{T}\mathbf{\Delta
T}_{ml}$

where
$\kappa_{\mathrm{bal}}=L\sum_{ml}\mathbf{G}_{ml}\cdot\mathbf{\hat{n}_{\mathrm{ext}}}$
is the ballistic thermal conductivity, $\mathbf{\hat{n}}_{\mathrm{ext}}$ being
the direction of the applied temperature gradient; the matrices
$\mathbf{K}_{ml}$ are

$K_{iml}=\frac{L}{A\Delta
T_{\mathrm{ext}}}\sum_{j}\mathbf{S}_{ml}\cdot\mathbf{R}^{K}_{ij},$ (19)

where we use $\mathrm{ReLu}(x)-\mathrm{ReLu}(-x)=x$, and

$\displaystyle\mathbf{R}_{ij}^{K}$ $\displaystyle=$
$\displaystyle\begin{cases}A_{ij}\mathbf{\hat{n}}_{ij},&\text{if }i\text{ and
}j\text{ share a hot contact's side}\\\ 0,&\text{otherwise}.\end{cases}$ (20)

Figure 4: a) A scheme for two adjacent elements. In this case, the vectorial
MFP $\mathbf{F}_{\mu}$ points inward the element $i$, thus it contributes to
the off-diagonal part of the stiffness matrix of the BTE model. b) Local
reference system, $\mathbf{\hat{e}}_{\xi}-\mathbf{\hat{e}}_{\eta}$, used to
derive the non-orthogonal contribution to the heat flux for the heat-
conduction equation.

Lastly, it is also possible to compute the “mode-resolved” thermal
conductivity, given by

$\kappa^{\mathrm{eff}}_{\mu}=\frac{\kappa_{\mathrm{bal}}}{2}+\sum_{ml}\mathbf{K}_{\mu}^{T}c_{\mu}^{ml}\mathbf{\Delta
T}_{ml},$ (21)

where $c_{\mu}^{ml}$ maps $\mathbf{F}_{ml}$ to $\mathbf{F}_{\mu}$, and

$K_{i\mu}=\frac{L}{A\Delta
T_{\mathrm{ext}}}\sum_{j}\mathbf{S}_{\mu}\cdot\mathbf{R}^{K}_{ij}.$ (22)

##  Appendix B Heat conduction equation

The first guess to Eq. 4 is given by the standard diffusive equation
$\nabla\kappa\nabla\Delta T(\mathbf{r})=0$, with $\kappa$ being the bulk
thermal conductivity tensor. In a finite-volume fashion, we integrate over the
control volume $\Omega_{c}$

$\int_{\partial\Omega_{c}}\sum_{j}\kappa_{ij}\frac{\partial\Delta
T(\mathbf{r})}{\partial x_{j}}\hat{n}_{i}=0,$ (23)

where we used Gauss’ law. Eq. 23 is then transformed into a sum of integrals
over individual faces of the volume, i.e.

$\sum_{j}\mathbf{g}_{ij}\cdot\mathbf{S}_{ij}=0$ (24)

where $S_{ij}=\sum_{\beta}\kappa_{\alpha\beta}R_{ij,\alpha}$ and
$\mathbf{g}_{ij}$ is the gradient evaluated at the face between the volume $i$
and $j$ (which hereafter will be referred to as $s$). For orthogonal grids,
flux discretization simply leads to $\nabla T=\left(\Delta T_{j}-\Delta
T_{i}\right)\mathbf{d_{ij}}|\mathbf{d_{ij}}|^{-2}$, where $\mathbf{d}_{ij}$ is
the distance between the centroids of the two adjacent volumes. However, for
unstructured mesh, a non-orthogonal contribution must be added. In `OpenBTE`,
we follow the approach described by Murthy [27], described as follows. With no
loss of generality, we consider a 2D case. As shown in Fig. 4b, we introduce a
reference system based on the vectors $\mathbf{e}_{\eta}$ and
$\mathbf{e}_{\xi}$. The former is aligned with $\mathbf{d}_{ij}$ and latter is
parallel to the face. In this new system, the gradient transforms as

$\frac{\partial\Delta
T}{\partial\mathbf{x}}=\mathcal{J}(\mathbf{x},\mathbf{e})\frac{\partial\Delta
T}{\partial\mathbf{e}},$ (25)

with $\mathcal{J}$ being the Jacobian. Here we will not elaborate on its
structure but will only provide the final formula for the temperature gradient
[27]

$\mathbf{g}_{ij}=\frac{\mathbf{\hat{n}}_{ij}}{\mathbf{\hat{n}}_{ij}\cdot\mathbf{\hat{e}}_{\xi}}\left[\frac{\partial\Delta
T}{\partial
e_{\xi}}\bigg{|}_{s}-\mathbf{\hat{e}}_{\xi}\cdot\mathbf{\hat{e}}_{\eta}\frac{\partial\Delta
T}{\partial e_{\eta}}\bigg{|}_{s}\right].$ (26)

The gradient along $\mathbf{e}_{\xi}$ at the face $s$ is simply

$\frac{\partial\Delta T}{\partial e_{\xi}}\bigg{|}_{s}=\frac{\Delta
T_{j}-\Delta T_{i}-R_{ij}^{P}}{|\mathbf{d_{ij}}|},$ (27)

where

$\displaystyle R_{ij}^{P}$ $\displaystyle=$ $\displaystyle\begin{cases}-\Delta
T_{\mathrm{ext}},&\text{if }i\in\mathrm{hot}\text{ and
}j\in\mathrm{cold}\text{ and they share a side}\\\ \Delta
T_{\mathrm{ext}},&\text{if }i\in\mathrm{cold}\text{ and
}j\in\mathrm{hot}\text{ and they share a side}\\\
0&\text{otherwise}\end{cases}$ (28)

Using Eq. 27, we rewrite Eq. 26 as

$\mathbf{g}_{ij}=\frac{\mathbf{\hat{n}}_{ij}}{\mathbf{\hat{n}}_{ij}\cdot\mathbf{d}_{ij}}\left[\Delta
T_{j}-\Delta T_{i}-R_{ij}^{P}\right]+\mathbf{F}_{ij},$ (29)

where $\mathbf{F}_{ij}$ is called “secondary flux.” In contrast, Eq. 27 is
termed “primary flux.” The secondary flux is treated explicitly, i.e. using
precomputed temperatures. To avoid unambiguities in the 3D case, we compute
$\mathbf{F}$ as the difference between the total flux and the primary flux,
i.e.

$\mathbf{F}_{ij}=\mathbf{g}_{ij}-\frac{\mathbf{\hat{n}}_{ij}}{\mathbf{\hat{n}}_{ij}\cdot\mathbf{d}_{ij}}\left(\mathbf{g}_{ij}\cdot\mathbf{d}_{ij}\right).$
(30)

Combining Eqs. 30-29, and explicitly indicating the iteration index, we have

$\mathbf{g}_{ij}^{(n)}=\mathbf{g}_{ij}^{(n-1)}+\frac{\mathbf{\hat{n}}_{ij}}{\mathbf{\hat{n}}_{ij}\cdot\mathbf{d}_{ij}}\left[\left(\Delta
T_{j}^{(n)}-\Delta
T_{i}^{(n)}\right)-\mathbf{g}_{ij}^{(n-1)}\cdot\mathbf{d}_{ij}\right].$ (31)

In light of these results, Eq. 24 translates into the following linear system

$\mathbf{A}\mathbf{\Delta T}=\mathbf{B}+\mathbf{P},$ (32)

with stiffness matrix and right hand side

$\displaystyle A_{ij}$ $\displaystyle=$
$\displaystyle\mathbf{S}_{ij}\cdot\frac{\mathbf{n}_{ij}}{\mathbf{\hat{n}}_{ij}\cdot\mathbf{d}_{ij}}-\delta_{ij}\sum_{k}\mathbf{S}_{ik}\cdot\frac{\mathbf{n}_{ik}}{\mathbf{\hat{n}}_{ij}\cdot\mathbf{d}_{ik}},$
$\displaystyle B_{i}$ $\displaystyle=$
$\displaystyle\sum_{j}\mathbf{g}_{ij}\cdot\left(\frac{\mathbf{d}_{ij}(\mathbf{\hat{n}}_{ij}\cdot\mathbf{S}_{ij})}{\mathbf{\hat{n}}_{ij}\cdot\mathbf{d}_{ij}}-\mathbf{S}_{ij}\right).$
(33)

The perturbation term is given by

$P_{i}=\sum_{j}\mathbf{S}_{ij}^{P}\cdot\frac{\mathbf{n}_{ij}}{\mathbf{\hat{n}}_{ij}\cdot\mathbf{d}_{ij}},$
(34)

where $\mathbf{S}_{ij}^{P}=\mathbf{S}_{ij}R_{ij}^{P}$. The gradients
$\mathbf{g}_{ij}$, appearing in Eq. B, depend on precomputed temperatures.
Their calculation if performed by interpolating the gradients from the
elements $i$ and $j$, i.e.

$\mathbf{g}_{ij}=\omega_{ij}\mathbf{g}_{i}+\left(1-\omega_{ij}\right)\mathbf{g}_{j},$
(35)

where

$\mathbf{g}_{i}=\frac{\partial\Delta T}{\partial\mathbf{x}}\bigg{|}_{i},$ (36)

and $\omega_{ij}$ is the distance between the centroid of the element $i$ and
the intersection between $\mathbf{d}_{ij}$ and the face $s$. The gradient at
the element’s centroid is evaluated assuming local linear variation of $\Delta
T$ around the element $i$, resulting in as many linear equations as the
neighbours of $i$,

$\Delta T_{i}+\mathbf{g}_{i}\cdot\mathbf{d}_{ij}=\Delta T_{j}+R_{ij}^{P},$
(37)

which, in matrix form, become

$\mathbf{M}_{i}\mathbf{g}_{i}=\mathbf{Q}_{i}.$ (38)

Equation 38 is a linear system to be solved for each element. Since
$\mathbf{M}_{i}$ has more rows (number of neighbors) that columns
(dimensions), we solve it using the least square method, i.e.

$\mathbf{g}_{i}=\left(\mathbf{M}_{i}^{T}\mathbf{M}_{i}\right)^{-1}\mathbf{M}_{i}^{T}\mathbf{Q}_{i}.$
(39)

Combining Eqs. 39-B and-32, we may formalize the heat conduction equation in
unstructured grids as the following iterative linear system

$\mathbf{A}\mathbf{\Delta T}^{(n)}=\mathbf{S}\mathbf{\Delta
T}^{(n-1)}+\mathbf{P},$ (40)

as reported in the main text. Lastly, after Eq. 40 is solved,
$\kappa_{\mathrm{fourier}}$ is computed by

$\kappa_{\mathrm{fourier}}=-\frac{L}{\Delta
T_{\mathrm{ext}}A}\sum_{s\in\mathrm{hot}}\mathbf{g}_{ij}\cdot\mathbf{S}_{ij}.$
(41)

## References

  * [1] L. Lindsay, A. Katre, A. Cepellotti, N. Mingo, Perspective on ab initio phonon thermal transport, Journal of Applied Physics 126 (5) (2019) 050902.
URL https://www.osti.gov/pages/biblio/1550742

  * [2] L. Lindsay, C. Hua, X. Ruan, S. Lee, Survey of ab initio phonon thermal transport, Materials Today Physics 7 (2018) 106–120.
URL https://www.sciencedirect.com/science/article/abs/pii/S254252931830141X

  * [3] N. Mingo, D. Stewart, D. Broido, L. Lindsay, W. Li, Ab initio thermal transport, Length-scale dependent phonon interactions (2014) 137–173.
URL https://link.springer.com/chapter/10.1007/978-1-4614-8651-0_5

  * [4] W. Li, J. Carrete, N. A. Katcho, N. Mingo, ShengBTE: A solver of the Boltzmann transport equation for phonons, Comput. Phys. Commun. 185 (6) (2014) 1747–1758.
URL http://www.sciencedirect.com/science/article/pii/S0010465514000484

  * [5] M. Omini, A. Sparavigna, An iterative approach to the phonon boltzmann equation in the theory of thermal conductivity, Physica B: Condensed Matter 212 (2) (1995) 101–112.
URL https://www.sciencedirect.com/science/article/abs/pii/0921452695000163

  * [6] D. A. Broido, M. Malorny, G. Birner, N. Mingo, D. Stewart, Intrinsic lattice thermal conductivity of semiconductors from first principles, Applied Physics Letters 91 (23) (2007) 231922.
  * [7] K. Esfarjani, G. Chen, H. T. Stokes, Heat transport in silicon from first-principles calculations, Phys. Rev. B 84 (8) (2011) 085204.
URL https://journals.aps.org/prb/abstract/10.1103/PhysRevB.84.085204

  * [8] W. Li, N. Mingo, L. Lindsay, D. A. Broido, D. A. Stewart, N. A. Katcho, Thermal conductivity of diamond nanowires from first principles, Physical Review B 85 (19) (2012) 195436.
URL https://journals.aps.org/prb/abstract/10.1103/PhysRevB.85.195436

  * [9] J. Carrete, B. Vermeersch, A. Katre, A. van Roekeghem, T. Wang, G. K. Madsen, N. Mingo, almaBTE : A solver of the space–time dependent Boltzmann transport equation for phonons in structured materials, Comput. Phys. Commun. 220 (2017) 351–362.  doi:10.1016/j.cpc.2017.06.023.
  * [10] Z. Han, X. Yang, W. Li, T. Feng, X. Ruan, Fourphonon: An extension module to shengbte for computing four-phonon scattering rates and thermal conductivity, arXiv preprint arXiv:2104.04895 (2021).
URL https://journals.aps.org/prb/abstract/10.1103/PhysRevB.93.045202

  * [11] T. Feng, X. Ruan, Quantum mechanical prediction of four-phonon scattering rates and reduced thermal conductivity of solids, Physical Review B 93 (4) (2016) 045202\.
  * [12] A. Togo, L. Chaput, I. Tanaka, Distributions of phonon lifetimes in brillouin zones, Phys. Rev. B 91 (2015) 094306.  doi:10.1103/PhysRevB.91.094306.
  * [13] L. Chaput, Direct solution to the linearized phonon boltzmann equation, Physical Review Letters 110 (26) (6 2013).  doi:10.1103/PhysRevLett.110.265506.
  * [14] L. Paulatto, F. Mauri, M. Lazzeri, Anharmonic properties from a generalized third-order ab initio approach: Theory and applications to graphite and graphene, Physical Review B 87 (21) (2013) 214303\.
URL https://journals.aps.org/prb/abstract/10.1103/PhysRevB.87.214303

  * [15] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, et al., Quantum espresso: a modular and open-source software project for quantum simulations of materials, Journal of physics: Condensed matter 21 (39) (2009) 395502.
URL https://iopscience.iop.org/article/10.1088/0953-8984/21/39/395502/meta

  * [16] L. Paulatto, D. Fournier, M. Marangolo, M. Eddrief, P. Atkinson, M. Calandra, Thermal conductivity of bi 2 se 3 from bulk to thin films: Theory and experiment, Physical Review B 101 (20) (2020) 205419.
URL https://journals.aps.org/prb/abstract/10.1103/PhysRevB.101.205419

  * [17] M. Simoncelli, N. Marzari, F. Mauri, Unified theory of thermal transport in crystals and glasses, Nature Physics 15 (8) (2019) 809–813.
URL https://www.nature.com/articles/s41567-019-0520-x

  * [18] G. Barbalinardo, Z. Chen, N. W. Lundgren, D. Donadio, Efficient anharmonic lattice dynamics calculations of thermal transport in crystalline and disordered solids, Journal of Applied Physics 128 (13) (2020) 135104.
URL https://aip.scitation.org/doi/full/10.1063/5.0020443

  * [19] L. Isaeva, G. Barbalinardo, D. Donadio, S. Baroni, Modeling heat transport in crystals and glasses from a unified lattice-dynamical approach, Nature communications 10 (1) (2019) 1–6.
URL https://www.nature.com/articles/s41467-019-11572-4

  * [20] C. D. Landon, N. G. Hadjiconstantinou, Deviational simulation of phonon transport in graphene ribbons with ab initio scattering, J. Appl. Phys. 116 (16) (2014) 163502.
URL https://aip.scitation.org/doi/full/10.1063/1.4898090

  * [21] J.-P. M. Péraud, N. G. Hadjiconstantinou, Efficient simulation of multidimensional phonon transport using energy-based variance-reduced monte carlo formulations, Physical Review B 84 (20) (2011) 205331\.
URL https://journals.aps.org/prb/abstract/10.1103/PhysRevB.84.205331

  * [22] A. Pathak, A. Pawnday, A. P. Roy, A. J. Aref, G. F. Dargush, D. Bansal, Mcbte: A variance-reduced monte carlo solution of the linearized boltzmann transport equation for phonons, Computer Physics Communications 265 (2021) 108003.
  * [23] G. Romano, Efficient calculations of the mode-resolvedab-initiothermal conductivity innanostructures, arXiv preprint arXiv:2002.08940 (2020).
URL https://arxiv.org/abs/2002.08940

  * [24] G. Romano, Openbte, https://github.com/romanodev/OpenBTE (2013).
  * [25] G. Romano, Phonon transport in patterned two-dimensional materials from first principles, arXiv preprint arXiv:2002.08940 (2020).
URL https://arxiv.org/abs/2002.08940

  * [26] C. Hua, A. J. Minnich, Transport regimes in quasiballistic heat conduction, Physical Review B 89 (9) (2014) 94302\.
URL https://journals.aps.org/prb/abstract/10.1103/PhysRevB.89.094302

  * [27] J. Y. Murthy, S. Mathur, Numerical methods in heat, mass, and momentum transfer, School of Mechanical Engineering Purdue University (2002).
URL
https://engineering.purdue.edu/~scalo/menu/teaching/me608/ME608_Notes_Murthy.pdf

  * [28] S. A. Ali, G. Kollu, S. Mazumder, P. Sadayappan, A. Mittal, Large-scale parallel computation of the phonon boltzmann transport equation, International journal of thermal sciences 86 (2014) 341–351.
URL https://www.sciencedirect.com/science/article/pii/S1290072914002233

  * [29] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, et al., Array programming with numpy, Nature 585 (7825) (2020) 357–362.
  * [30] L. D. Dalcin, R. R. Paz, P. A. Kler, A. Cosimo, Parallel distributed computing using python, Adv. Water Resour. 34 (9) (2011) 1124–1139.
URL https://www.sciencedirect.com/science/article/pii/S0309170811000777

  * [31] C. Geuzaine, J.-F. Remacle, Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities, Int. J. Numer. Meth. Eng. 79 (11) (2009) 1309–1331.
URL http://onlinelibrary.wiley.com/doi/10.1002/nme.2579/full

  * [32] J. Ahrens, B. Geveci, C. Law, Paraview: An end-user tool for large data visualization, The visualization handbook 717 (8) (2005).
URL https://www.paraview.org/

  * [33] D. Hasselman, L. F. Johnson, Effective thermal conductivity of composites with interfacial thermal barrier resistance, J. Compos. Mater. 21 (6) (1987) 508–515.
URL https://aip.scitation.org/doi/full/10.1063/1.4945776