A deterministic solver for the phonon Boltzmann transport equation (BTE) that retains the full phonon–phonon scattering matrix W — rather than the diagonal relaxation-time approximation (RTA) — for thermal transport in geometries ranging from one-dimensional slabs to a three-dimensional FinFET-like device.

This repository contains the solver code accompanying the manuscript "Deterministic full-matrix phonon BTE solver for 3D device geometries" (Y. S. Ju). It reproduces the full-scattering-matrix construction, the hybrid mode-space solver, and the device geometries reported there.

The interatomic model used throughout is the Stillinger–Weber (SW) potential for silicon. Material-specific transport magnitudes (e.g., the bulk conductivity and the device temperature-rise correction) are properties of this SW model and should not be read as quantitative predictions for crystalline silicon; the structural results of the paper are independent of this choice.

---

What the solver does

The phonon BTE under study is

v_λ · ∇ f_λ(r) = Σ_μ W_λμ f_μ(r) + Q_λ(r),

where f_λ is the deviational phonon distribution for mode λ, v_λ the group velocity, W the full scattering operator, and Q a heat source. The defining feature of this solver is that it keeps the complete off-diagonal Wwhile remaining tractable in 3D, by a hybrid mode-space decomposition:

• Streaming (v · ∇) is handled by dense, upwind spatial sweeps, which are trivially mode-local and need no compression.
• Scattering (W) is applied as a compressed operator — either a truncated SVD of the in-scattering block (with the out-scattering diagonal carried exactly) or a tensor-train (TT) representation — with grid-independent rank.

The scattering matrix is constructed to satisfy three physical constraints that are load-bearing for a stable device solve:

1. Energy conservation, W c = 0, where c are the modal heat capacities.
2. Detailed balance, c_i W_ij = c_j W_ji.
3. Negative semi-definiteness on the non-equilibrium subspace.

Conservation is not optional: a non-conserving operator makes the device solver diverge. The construction enforces all three analytically and then projects out the residual, reaching machine precision on each (see "Verifying the install "
below).

---

Repository layout

__init__.py            Sets the version and its description
core/                  Tensor-train algebra (physics-agnostic)
  tensor_train.py        TensorTrain class, index bookkeeping
  tt_operators.py        TTOperator (MPO), add/scale/compress
  tt_solvers.py          AMEn linear solver

physics/               Phonon model and scattering operator
  sw_potential.py        Stillinger–Weber energy and forces
  si_physical_analysis_v2.py
                         Force constants, phonons, and the full-W
                         scattering-matrix construction (the core
                         physical contribution)
  sw_phonon_bridge.py    SWPhononComputation: end-to-end SW → (v, τ, C, W)
  phonon_data.py         PhononBandData container and reference models
  mode_space.py          ModeSpaceData: the mode-resolved model consumed
                         by every solver; factory mode_space_from_sw()
  svd_scattering.py      SVDScatteringOperator: low-rank W with exact γ
  tt_scattering.py       TTScatteringOperator: tensor-train W

solvers/               BTE solvers (hybrid mode-space)
  bte_slab_1d_hybrid.py  1D slab (validation against analytic reference)
  bte_slab_2d_hybrid.py  2D slab
  bte_monolithic_2d.py   2D monolithic domain
  dense_sweep_3d_hybrid.py
                         3D box dense-sweep reference
  bte_3d_finfet.py       3D FinFET-like device (fin + base, Schwarz
                         coupling) — the headline device geometry
  dsa_operator.py        Diffusion-synthetic acceleration (base)
  dsa_operator_1d.py     DSA, 1D
  dsa_operator_2d.py     DSA, 2D
  dsa_operator_3d.py     DSA, 3D (Fourier-cosine)

Diffusion-synthetic acceleration (DSA) is an optional convergence accelerator and is off by default; it is not part of the production methodology and is provided only for cross-checking.

All the code files are in Dryad_tt_bte_solver.zip.

---

Requirements

• Python ≥ 3.10
• NumPy
• SciPy

That is the complete runtime requirement. The solver runs entirely on
CPU/NumPy.

JAX is an optional dependency, used only by a single GPU-resident DSA correction method (DSAOperator3D.apply_correction_jax). It is not needed to run any solver or to reproduce any result in the paper; the relevant import is lazy, so the package loads and runs without JAX installed.

Install the requirements with:

pip install numpy scipy

---

Installing and importing

This is a flat package: put the directory that contains core/,  physics/,
and solvers/ on your Python path and import the subpackages directly. The simplest approach is to run Python from inside this directory, or to add it to
PYTHONPATH:

export PYTHONPATH=/path/to/this/directory:$PYTHONPATH

Then, in Python:

from physics import SWPhononComputation, SVDScatteringOperator, ModeSpaceData
from physics.mode_space import mode_space_from_sw
from solvers import BTESlab1D, DenseSweep3DHybrid
from solvers.bte_3d_finfet import DenseSweep3DFinFET

---

Verifying the install

The following builds the SW phonon model on a small grid, constructs the full scattering matrix, and confirms the conservation, detailed-balance, and semi-definiteness gates. It takes a couple of minutes (most of the time is the force-constant computation) and requires only NumPy and SciPy.

from physics import SWPhononComputation

comp = SWPhononComputation(q_grid_size=3)      # small 3x3x3 q-grid
comp.compute(build_W=True, verbose=True)       # builds v, tau, C, and W

# After construction the solver prints the matrix gates, e.g.:
#   Energy conservation ||W c|| / ||c||      ~ 1e-16   (Wc = 0)
#   Detailed balance     max|c_i W_ij - ...| ~ 1e-16
#   Eigenvalues          all <= 0            (semi-definite)
print("modes:", comp.W.shape[0])

A clean run prints corrected conservation and detailed-balance residuals at the level of machine precision (~1e-15 or smaller) and no positive eigenvalues.

---

Typical workflow

1. Build the phonon model

from physics import SWPhononComputation
from physics.mode_space import mode_space_from_sw

comp = SWPhononComputation(q_grid_size=7)      # production grids use larger N
comp.compute(build_W=True)

mode_data = mode_space_from_sw(comp)           # ModeSpaceData for the solvers
print(mode_data.summary())

ModeSpaceData exposes the mode-resolved velocities, heat capacities, the full
W, the RTA operator W_rta, the thermally active mask, and helpers such as
bulk_thermal_conductivity() and project_to_1d().

2. Compress the scattering operator (optional but recommended)

from physics import SVDScatteringOperator

svd_op = SVDScatteringOperator.from_mode_data(mode_data, rank=50)
svd_op.print_rank_summary()

The out-scattering diagonal γ is retained exactly; only the in-scattering
block is approximated, so the transport-critical relaxation rates are preserved.

3. Solve a geometry

1D slab (used for validation against the analytic Fuchs–Sondheimer reference):

from solvers import BTESlab1D

slab = BTESlab1D(N_x=100, L=200e-9, mode_data=mode_data.project_to_1d())
result = slab.solve(tol=1e-8, verbose=True)

3D FinFET-like device (the headline geometry):

from solvers.bte_3d_finfet import DenseSweep3DFinFET
# Construct with the fin/base spatial resolution and geometry, supply the
# mode_data and (optionally) the SVD-compressed operator, then solve.
# See the DenseSweep3DFinFET docstring for the full constructor signature
# and boundary-condition options (isothermal sink, diffuse side walls,
# y-uniform line heat source).

Each solver returns a result object carrying the converged temperature field, 
the heat flux, and solver diagnostics.

---

Reproducing the paper's headline quantities

The device result is the full-W correction to the peak temperature rise, defined as a same-grid difference between the full-W and RTA solutions on the FinFET-like geometry with a y-uniform line heat source on the top of the fin. For Brillouin-zone grids,s N ≥ 11 the correction is

δT_fW / (T_RTA − T_0) = 11.7 ± 0.3 %,

a bounded, grid-robust correction. The RTA baseline used for this comparison carries the same isotope channel as the full operator, so the two solutions differ only in the off-diagonal scattering.

To reproduce a single point: build the model at the desired grid. N
(q_grid_size=N), construct both the full and RTA operators from the same
ModeSpaceData solve the FinFET geometry with each, and take the difference of the peak temperatures. The bulk conductivity of the SW model on this grid is
k_bulk = 148.0 W m⁻¹ K⁻¹ (isotropic to better than 0.003%), which is a useful
cross-check that the phonon model is built correctly.

---

Numerical notes and known properties

• Γ-centred Monkhorst–Pack sampling. On Γ-centred grids, the absolute device temperatures converge non-monotonically at small odd N (the grid does not close under q → −q); monotone behaviour sets in by N ≈ 11–15. The
  same-grid correction is far less sensitive and is reported as a bounded band.
• SW group-velocity anisotropy. On a Γ-centred gri,d the SW group velocities carry a fixed [111] ellipticity, giving a non-negligible off-diagonal
  Cartesian conductivity (κ_xy/κ_xx ≈ 0.56). The scalar bulk conductivity is unaffected (the tensor is traceless-clean). The y-uniform line source usedfor thee device results makes∂T/∂y ≈ 0, so this term is geometrically inactive.
• Scattering-matrix scaling. The number of three-phonon channels, and hence the nonzeros ofW, scales as the square of the active mode count, so the operator is genuinely dense. Compression of the solution (not the operator)
  is what makes the solver efficient.

If you use this solver, we encourage you to cite the accompanying paper "Deterministic full-matrix phonon BTE solver for 3D device geometries," by Y. S. Ju in Scientific Reports (2026).