Guest bloggers Enrico Di Lucente, a postdoctoral research scientist at Columbia University, and Michele Simoncelli, an assistant professor at Columbia University, discuss how they used the Physics Builder in COMSOL Multiphysics® to create a phonon hydrodynamics custom physics interface for exploring nondiffusive heat transport in realistic geometries.
In dielectric materials with ultrahigh thermal conductivity, heat can violate diffusion and behave fluid-like, forming vortices, displaying temperature waves, and even locally backflowing against the temperature gradient. In this blog post, we introduce a phonon hydrodynamics custom physics interface that we developed for the COMSOL Multiphysics® software and show how continuum modeling can capture viscous heat transport phenomena beyond Fourier’s law.
Why Model Heat Beyond Fourier’s Law?
Fourier’s law has long been the standard framework for describing heat conduction as it successfully captures thermal transport in common materials and devices at room temperature and at millimeter or larger scales. However, recent experimental and theoretical studies have shown that Fourier’s picture can break down in microscale devices made of ultrapure, high–thermal-conductivity materials such as diamond, graphite, and hexagonal boron nitride.
In these systems, phonons — the primary heat carriers — undergo frequent momentum-conserving collisions within certain temperature ranges, which can extend from cryogenic (~70 K) up to near room temperature. This regime is called phonon hydrodynamics. As a result, heat transport becomes collective and fluid-like, leading to phenomena such as temperature waves (“second sound”), Poiseuille-like heat flow (i.e., faster in the center of a channel and slower at the boundaries), thermal backflow, and heat vortices.
Capturing these effects requires a model that goes beyond standard Fourier diffusion while remaining practical for realistic geometries. This is precisely the motivation behind the viscous heat equations (VHE) (Refs. 1 and 2) and their implementation in COMSOL Multiphysics®.
Origin of the Viscous Heat Equations
At a deeper level, the VHE originate from a systematic coarse-graining of the linearized phonon Boltzmann transport equation (LBTE), where the complex microscopic phonon dynamics are projected onto a small set of local-equilibrium fields: the temperature T(r,t) and a drift velocity u(r,t). This drift velocity emerges because, in the hydrodynamic regime, momentum-conserving (normal) phonon collisions dominate over momentum-relaxing (Umklapp) processes, allowing phonons to collectively carry crystal momentum. As a result, heat transport is no longer governed solely by temperature gradients but also by the evolution of this momentum field. The VHE therefore consist of two coupled equations:
- An energy balance equation for T(r,t)
- A momentum balance equation for u(r,t), analogous to the linearized Navier–Stokes equations
Crucially, two transport coefficients arise from the underlying LBTE symmetries: the thermal conductivity (associated with the odd-parity part of the phonon distribution) and the thermal viscosity (associated with the even-parity part). While the former describes diffusive heat flow, the latter accounts for momentum diffusion and viscous stresses in the phonon fluid, enabling phenomena such as heat vortices, backflow, and nonlocal transport. This framework naturally interpolates between regimes — it reduces to Fourier’s law when momentum is strongly dissipated and recovers previously proposed hydrodynamic models such as the Guyer–Krumhansl and dual-phase-lag equations as limiting cases, while remaining fully applicable to realistic materials with complex phonon dispersions.
From Microscopic Theory to a Continuum Model
At the microscopic level, phonon transport is described by the Boltzmann transport equation (BTE). While predictive, the BTE is extremely expensive to solve in complex geometries and is therefore ill-suited for device-scale modeling.
The viscous heat equations provide a mesoscopic alternative. They extend Fourier’s law by introducing an additional field: the phonon drift velocity. In this sense, they represent the phonon analogue of the linear Navier–Stokes equations for classical laminar fluids but applied to quantum phonon fluids. In this framework, temperature evolution is coupled to momentum balance equations for phonons, allowing the model to naturally interpolate between diffusive and hydrodynamic regimes.
A key feature of this approach is that all material parameters — such as thermal conductivity, viscosity, and relaxation rates — are determined from a single solution of the linearized BTE using first-principles calculations. These parameters are then incorporated into the continuum model, enabling efficient simulations in arbitrary geometries without compromising physical fidelity while retaining full quantum mechanical and ab initio accuracy.
Implementing the Viscous Heat Equations of Phonon Hydrodynamics
To make this framework accessible, we developed a custom physics interface using the Physics Builder in the COMSOL Multiphysics® software. Once installed, the interface appears as Viscous Heat Equations (VHE).
Features Available After Implementing
- A coupled temperature–velocity formulation of heat transport
- Steady-state and time-dependent simulations
- Boundary conditions tailored to hydrodynamic heat flow
- Compatibility with complex geometries and mesoscopic devices
Installing the Interface
The installation follows the standard Physics Builder workflow:
- Enabling the Physics Builder in COMSOL preferences
- Importing the provided builder file
- Adding the Viscous Heat Equations (VHE) interface from My Physics Interfaces
Get step-by-step information on how to install this interface by expanding the section below.
If the button is present, proceed to Step 3; otherwise, continue to Step 2.
Step 2. From the File menu, select Preferences. In the Preferences window, choose Physics Builder from the list on the left and then select the Enable Physics Builder checkbox:
Click OK and restart COMSOL Multiphysics® for the changes to take effect.
Step 3. Once the Physics Builder button appears, indicating that the interface is ready to use, open an existing model or create a new blank model. From the Windows menu, select Physics Builder Manager:
In the Physics Builder Manager, under Archive Browser, right-click Development Files and select Add Builder File:
Then import the vhe.mphphb model file.
Step 4. Open the Add Physics window and navigate to the My physics interfaces section where the interface is listed as Viscous Heat Equations (VHE):
Then select Viscous Heat Equations (VHE) to add it to the model.
Governing Equations and Boundary Conditions
The interface directly implements the viscous heat equations as a coupled system:
- One energy conservation equation for temperature
- Three momentum balance equations for the phonon drift velocity components
From the user’s perspective, these equations behave like any other COMSOL physics interface and can be combined with standard meshing, solvers, and postprocessing tools.
Boundary Conditions Tailored to Heat Flow Physics
Several boundary conditions are provided to model different physical regimes:
- Fixed temperature boundaries for ideal thermal reservoirs
- Temperature-gradient flux boundaries that constrain only the diffusive component
- Velocity constraints to control phonon momentum flow
- Slip boundaries, representing specular phonon reflection
- No-slip boundaries, modeling fully diffusive phonon scattering
These options allow users to continuously tune the simulation from purely diffusive to strongly hydrodynamic behavior.
The temperature boundary condition setup and slip and no-slip boundary conditions appear in the Viscous Heat Equations (VHE) interface as follows:
Example 1: Fourier Diffusion vs. Hydrodynamic Heat Flow
To illustrate the qualitative differences between diffusive and viscous heat transport, consider a simple two-rectangle geometry. A vertical temperature gradient is applied across the main domain, while a smaller side region is laterally connected.
Under Fourier’s law, heat flows directly from hot to cold, producing smooth isotherms. When the viscous heat equations are used instead, a very different picture emerges: Heat recirculates in the smaller domain, forming vortex-like patterns and producing a small but finite thermal backflow. Below we show the comparison of steady-state temperature profiles obtained with Fourier’s law (left) and the viscous heat equations (right), showing hydrodynamic recirculation and thermal backflow.
Example 2: Thermal Vortices in the Incompressible Limit
The hydrodynamic nature of heat transport becomes even clearer in the incompressible flow limit of the VHE. In this regime, the divergence of the phonon drift velocity vanishes, and heat transport is dominated by momentum flow rather than diffusion.
Simulations show the spontaneous formation of steady-state thermal vortices (Ref. 3) — a behavior entirely absent in purely diffusion models.
Below we show the steady-state temperature profiles in the diffusive regime (top) and incompressible hydrodynamic limit (bottom), highlighting vortex formation and heat backflow (Ref. 3).
Time-Dependent Simulations: Hydrodynamic Heat in Motion
The interface also supports fully time-dependent simulations. By applying transient boundary conditions — such as oscillating drift velocities or alternating temperature gradients — it is possible to observe the formation, motion, and reversal of heat vortices over time.
These simulations reveal heat behaving much like a driven viscous fluid, responding dynamically to external forcing.
Why This Matters!
The interface provides a practical bridge between microscopic transport theory and continuum modeling. It enables the exploration of nondiffusive heat transport in realistic geometries that are inaccessible to fully microscopic methods.
By making hydrodynamic heat transport accessible within COMSOL Multiphysics®, the interface opens new possibilities for studying thermal phenomena in next-generation materials and devices.
To learn more about the interface and download it, visit its Application Exchange entry: Phonon Hydrodynamics Interface.
About the Guest Authors
Enrico Di Lucente is a postdoctoral research scientist in the Department of Applied Physics and Applied Mathematics at Columbia University in the City of New York, with a joint affiliation in materials science and engineering. He joined the research group of Prof. Michele Simoncelli in December 2025, immediately after completing his PhD in the Theory and Simulation of Materials (THEOS) group at EPFL, under the supervision of Prof. Nicola Marzari. He defended his doctoral thesis titled Theoretical and Computational Advances in Quantum and Hydrodynamic Thermal Transport.
His research focuses on thermal transport beyond Fourier’s law, phonon hydrodynamics, and first-principles modeling of heat transport and other fundamental and coupled excitations in condensed matter, ranging from magnons to light–matter interactions. His work bridges microscopic transport theory and continuum modeling, with the goal of enabling predictive simulations of nondiffusive heat transport in realistic materials and device geometries. It also aims to guide experimental efforts toward the design of innovative devices capable of detecting nonstandard quantum transport phenomena, with potential impact across emerging technologies including electronics, energy storage, fusion shielding, spintronics, and hypersonics. His broader research interests include condensed matter physics, computational physics, quantum transport, and magnetism.
Michele Simoncelli has been an assistant professor in the Department of Applied Physics and Applied Mathematics at Columbia University since January 2025. His group develops the theoretical and computational framework to understand, quantitatively describe, and control quantum transport phenomena in materials involving, e.g., charge, heat, light and spin, their possible synergies or conflicts, and related macroscopic signatures. Prior to joining Columbia, he held the Crone Research Fellowship in the Physics Department at the University of Cambridge (2021-2024). There, he worked on fundamental quantum theory and computational methods to describe the emergence of hybrid crystal–glass properties in materials with controlled degrees of atomistic disorder in, for example, chemical composition, bond network topology, or geometry. He received his PhD from EPFL (Switzerland) in 2021 under the supervision of Nicola Marzari, presenting in his thesis novel microscopic and mesoscopic theories of thermal transport in solids: the Wigner transport equation, generalizing the semiclassical Peierls–Boltzmann equation, and the viscous heat equations, generalizing Fourier’s law.
Acknowledgments
The development of the custom physics interface and the underlying theoretical work benefited from discussions and collaborations with Zhiyi Wang (Université Grenoble Alpes; THEOS, EPFL) and Prof. Nicola Marzari (THEOS, EPFL; Theory of Condensed Matter, Cavendish Laboratory, University of Cambridge).
References
- M. Simoncelli, N. Marzari, and A. Cepellotti, Generalization of Fourier’s law into viscous heat equations, Physical Review X 10, 011019 (2020); https://journals.aps.org/prx/abstract/10.1103/PhysRevX.10.011019
- Dragašević, B. Rajkov, and M. Simoncelli, Viscous heat backflow and temperature resonances in extreme thermal conductors. Physical Review Letters 136, 186302 (2026); https://doi.org/10.1103/nbbn-56hr
- Di Lucente, F. Libbi, and N. Marzari, Vortices and backflow in hydrodynamic heat transport, Physical Review Letters 136(5), 056307 (2026); https://journals.aps.org/prl/abstract/10.1103/g9dx-hjyn
Comments (0)