Transient Multiscale Workflow for Thermal Analysis of 3DHI Chip Stack

Abstract: Modern package designs make use of technologies such as backside power delivery (BSPD) and 3D stacked chiplets that require accounting for the heterogeneity in back end of the line (BEOL) structures in hot-spot prediction. Multiscale homogenization strategies have been demonstrated to be effective for steady-state simulations, however accurate 3D transient simulations that include BEOL structures remain an open challenge. In this work, we demonstrate a transient thermal workflow that accounts for the 3D heterogeneous structures in the BEOL for problems with strong- and weak- temporal scale separation under the assumption of temperature independent constitutive properties. Our workflow, based on Bloomfield et. al. 2025, automatically extracts, meshes, and homogenizes thermal properties from GDSII and OASIS files to construct thermal property maps. Property maps (heat capacity and conductivity) have been generated for a 1 mm by 1 mm SoC-style model die that was constructed with LibreLane for 100 by 100 grids with 5 micron by 5 micron representative volume elements (RVEs), and 50 by 50 grids with 10 micron by 10 micron RVEs. The expressions for a transient effective conductivity are provided and a demonstration of the impact of the transient effects are provided for a single RVE. Finally, transient conductivity maps have been provided for a time integration timestep of dt=0.001.

• The paper introduces a rigorous transient multiscale homogenization workflow that integrates direct microstructure extraction from EDA tools for realistic thermal simulations.
• The paper employs a variational finite-element approach to capture transient microstructural inertia and anisotropic BEOL effects in chip stacks.
• The paper validates the method against fully resolved BEOL models, achieving less than 0.4% error while significantly reducing computational complexity.

Transient Multiscale Workflow for Thermal Analysis of 3DHI Chip Stack

Introduction and Context

Thermal management in 2.5D and 3D heterogeneous integration (3DHI) has become increasingly challenging due to extreme spatial and temporal scale disparities intrinsic to modern stacked chip architectures. Capturing the influence of intricate back-end-of-line (BEOL) metallization and spatially-variable power densities typically demands either gross simplification or prohibitive computational expense. Traditional approaches—such as volume-averaged or 1D analytical homogenization, RC-network reduction, or rapid but less interpretable ML surrogates—often fall short when resolving both spatially local microstructure effects and transient, non-steady-state thermal phenomena. The presented work addresses these deficiencies by introducing a rigorous, variationally consistent transient multiscale homogenization workflow integrated with direct extraction and meshing of microstructure from hierarchical EDA layout files (GDSII/OASIS).

Methodological Foundation

Variational Multiscale Homogenization

The framework employs first-principles finite-element (FE) discretization at both macro and micro scales. Governing equations are enforced in strong and weak forms with no distinction between macro/micro variables at the level of presentation, facilitating variational consistency and avoiding the limitations of coordinate-split or 1D flow assumptions. Microstructural temperature fields are decomposed into mean, gradient, and fluctuation components. Scale transition, crucially, is derived using the Hill–Mandel principle, ensuring energy consistency between scales and accommodating the inertia induced by finite RVE-size—a critical aspect when time-scale separation is not asymptotically strong.

A persistent innovation is the imposition that $\langle \dot{\epsilon} \rangle$ and $$\langle \mathbf{q} - \dot{\epsilon} (\mathbf{x} - \bar{\mathbf{x}}) \rangle$$ are returned from micro to macro at every time step, allowing for exact transfer of transient inertia effects into the macroscale finite-element solution.

Automated RVE Extraction

The methodology makes extensive use of EDA toolchains for realistic microstructure modeling. Here, RVEs are algorithmically extracted from open-source SoC designs (LibreLane with Sky130A PDK) using only a standard GDSII and process stackup description. This enables direct mapping from design intent to simulation geometry, representing CPU, crypto, peripheral, and memory logic at a high spatial resolution.

Figure 1: Tile-based dominance map of the LibreLane-generated Sky130A SoC test vehicle, showing the spatial distribution of CPU, cryptographic, peripheral, and memory logic by standard-cell density.

Multiscale Transient Property Computation

Homogenized Properties and RVE Grid Effects

Homogenized thermal conductivity tensors and local volumetric heat capacities are computed at each RVE from the extracted mesh. Two RVE grid sizes are compared (5×5 and 10×10 μm), each completely spanning the 5.4 μm BEOL stack, and distributed uniformly over the 1 mm² die in dense spatial lattices. These are used to generate high-resolution, directionally dependent property maps for transient simulation.

Spatial Distribution of Properties

The resultant spatial maps reflect the anisotropic structure of the BEOL. As shown in the axis-aligned conductivity maps, in-plane (x) conductivities consistently exceed out-of-plane (z) values—a feature faithfully captured at both RVE grid scales due to the macroscopic alignment of high aspect ratio horizontal power-delivery lines. Volumetric heat capacity exhibits less spatial variation but is affected by the distribution of high-density metal regions.

Transient Integration and Weak Scale Separation

Strong Temporal Scale Separation

In the limit where the microstructure equilibrates instantaneously relative to the macroscale (strong timescale separation), precomputed steady-state effective conductivities are directly interpolated in the macroscale simulation. This regime is computationally efficient and compatible with commercial FE solvers.

Weak Temporal Scale Separation and Microstructural Inertia

When loading timescales approach microstructural response times, transient effects in homogenized flux become non-negligible. The exact variationally consistent transient homogenization incorporates additional terms involving both the spatial coordinate deviation and the time derivative of internal energy, introducing cross-coupling of macro temperature and gradients. For fast load ramps, this introduces an inertia-driven lag in heat flux response not captured by simpler models.

Figure 2: Impact of transient contributions to the homogenized conductivity for a single 10×10 μm RVE, illustrating when RVE thermal inertia manifests for increasing loading rates.

Macroscale Demonstration and Model Validation

Die-Scale Demonstration

The approach is demonstrated by applying a spatially varying surface heat-flux map as a Neumann boundary on the bottom of the die. The temperature field resulting from this heterogeneous loading correlates spatially with the logical block distribution, showing realistic hot-spot formation and propagation.

Figure 3: Applied heat flux map driving the transient thermal response in the macroscale simulation.

Validation Against Explicit BEOL Model

For model validation, steady-state results from homogenized meshes are benchmarked—at multiple RVE sizes—against a fully resolved 50×50 μm explicit BEOL finite-element model, using matched boundary conditions. Homogenized models systematically converge toward the explicit solution as RVE size increases, with the 15×15 μm RVE yielding less than 0.4% error in average bottom surface temperature, despite a dramatic reduction in computational complexity.

Figure 4: Fully resolved 50×50 μm BEOL reference geometry utilized for validation, highlighting the nontrivial complexity of the explicit metallization structure.

Figure 5: Comparison of top-surface temperature field between the fully resolved explicit BEOL model and homogenized RVE-based models.

Implications and Future Developments

The demonstrated workflow marks a critical advance in chip-scale transient thermal analysis by removing restrictive 1D heat flow assumptions and permitting the inclusion of full microstructural BEOL heterogeneity extracted directly from standard design flows. Strong validation against explicit-resolved models confirms the efficacy of RVE-based transient property maps under the assumption of temperature independent properties and fixed timestep integration. The approach is immediately compatible with commercial solvers, enabling enhanced design-time thermal evaluation without recourse to full explicit simulation—a major practical benefit for 3DI/2.5D IC stack layout optimization and reliability assessment.

A clear theoretical implication is the controlled treatment of RVE inertia and the emergence of additional macro-micro coupling terms in the transient regime; these will become relevant for even higher fidelity as integration density and power densities further increase. Future developments should include incorporation of temperature-dependent constitutive relations, quantification of the additional coupling vector terms, and adaptive time-step or RVE size selection to optimize computational efficiency while maintaining accuracy. Extension to coupled electro-thermal simulation domains and integration with physical-design automation tools is also foreseeable.

Conclusion

This work establishes an authoritative reference workflow for transient multiscale modeling of thermal phenomena in realistic, heterogeneous, 3D-stacked IC substrates. By rigorously bridging geometric hierarchy from standard design files to macroscale FE simulation, and by enabling accurate transient analysis without reductionist assumptions, it substantially improves both the theoretical and practical landscape for the design and reliability assessment of next-generation chip stacks.