Understand the
timing of heat.

Many nodes: the matrix model

Positive capacities and reciprocal passive conductances produce real, nonnegative decay rates. With each connected component linked to a fixed-temperature reservoir, all rates are positive. An isolated network has a zero mode: its mean temperature does not relax to a prescribed ambient.

There can be repeated rates. A slow mode may be weak or invisible for a particular forcing and sensor. The slowest network mode alone does not give every node’s rise time. Passive networks can contain geometric loops; reciprocal conduction does not prohibit loops.

Deriving Elmore delay and the higher moments

1. Choose the input vector b and output node j. Define H(s) = eⱼᵀ(K+sC)⁻¹b.
2. Let x = K⁻¹b, v = K⁻¹Cx, and w = K⁻¹Cv. Expanding around s = 0 gives H(s)/H(0) = 1 − (vⱼ/xⱼ)s + (wⱼ/xⱼ)s² − … .
3. The first moment is μ = vⱼ/xⱼ = ∫₀∞[1 − y(t)]dt for the normalized step response y. The second raw impulse-response moment is 2wⱼ/xⱼ; its variance is σ² = 2wⱼ/xⱼ − μ².
4. For one stepped reservoir and an insulated chain end, xᵢ = 1 for a unit step. Then μⱼ = Σᵢ Cᵢ Rshared(i,j), where Rshared is the resistance common to the two paths back to that reservoir.
5. For a uniform chain, μⱼ = rc[j(j+1)/2 + j(N−j)]. At the last node, μₙ = rcN(N+1)/2.

The moment is exact for the stated model; using it as a threshold time is an approximation. Neither t₅₀ = μ − 0.35σ nor t₁₀–₉₀ = 2.2σ is a universal identity. Compute the full response when thresholds matter. Gupta et al., Elmore delay bounds ↗

Delay locates the transition; rise time measures its duration

Here, threshold delay t₅₀ measures time from the applied step to 50% of the final temperature change. Elmore / first-moment delay μ is the area above the normalized step response. Rise time tᵣ = t₉₀ − t₁₀ measures the transition between 10% and 90%. A pure transport delay shifts both crossing times equally and leaves this rise time unchanged.

Exact formulas for a single lumped capacitance

Fractions are relative to the initial-to-final temperature change. These formulas also describe cooling when the response is normalized this way. A constant heat input changes the final temperature but not τ in this linear model.

A useful Elmore-style shortcut

This fits one exponential using the first moment. It is exact for a single RC model and an estimate for a multi-node network. It can overpredict rise time because part of μ represents a delayed onset. For internal heat input or multiple sinks, use the normalized matrix expression above for μ.

Two moments: include the response width

A practical refinement fits a gamma-shaped impulse response using the exact mean μ and variance σ². Define shape a = μ²/σ² and scale θ = σ²/μ. If P(a,z) is the regularized lower incomplete gamma function, the estimated normalized step is y(t) ≈ P(a,t/θ).

The browser evaluates the inverse numerically. With one lumped capacitance, a = 1 and θ = τ, recovering ln(9)τ exactly. For a narrow, nearly Gaussian impulse shape, a rough width estimate is 2.563σ; for a single exponential it is 2.197σ. Neither coefficient is universal. The gamma fit is a moment-matching approximation, not a guaranteed bound or an exact multi-node solution; inspect its error against time marching below.

FOUNDATIONS / HEAT BALANCE → TEMPERATURE RESPONSE → TIMING

From the heat-transfer equation to delay and rise time

A thermal response has no single universal “response time.” We first derive the temperature history, then define the timing quantity we want to extract. In this section time is measured from an ideal input step at t = 0. Temperatures refer to a stationary solid, with no bulk fluid advection or phase change.

A. Heat equation · B. Lumped model · C. Delay versus rise time · D. A true delay · E. Distributed networks

A. Start from local energy conservation and Fourier’s law

1. Balance a small fixed control volume. Let q be the conductive heat-flux vector (W/m²), q̇ᵥ the volumetric heat generation (W/m³), ρ the density (kg/m³), and cₚ the specific heat capacity (J/(kg·K)). Stored thermal energy increases through heat generation and net inward conduction:

For constant isotropic conductivity k and constant ρcₚ, this becomes ∂T/∂t = α∇²T + q̇ᵥ/(ρcₚ), where α = k/(ρcₚ). All terms in the energy equation have units W/m³. If k varies spatially, retain ∇·(k∇T); do not replace it by k∇²T without justification.

2. Specify initial and boundary conditions. The differential equation alone cannot determine a time response. Typical choices, with outward surface normal n, are:

A prescribed-temperature boundary and convection through a finite h are different models. For a slab with coordinate ξ = x/L and dimensionless time Fo = αt/L², the source-free equation becomes ∂T/∂Fo = ∂²T/∂ξ². This identifies L²/α as the diffusion time scale. A threshold time also depends on the dimensionless solution, boundary conditions, and sensor location.

B. Integrate the heat equation to derive lumped capacitance

1. Integrate over the body. Applying the divergence theorem gives an exact global balance for constant properties:

Here Q(t) = ∫ᵥq̇ᵥdV is total internal heating in watts; insulated surfaces contribute zero outward heat flow. The surface integral is over the heat-exchanging surface. A volume-averaged temperature alone does not close this exact equation because the surface temperature can differ from the average.

2. Make the lumped approximation. When internal gradients are small, approximate T(r,t) ≈ T(t), including at the surface. Define Cth = ∫ᵥρcₚdV = ρcₚV for a homogeneous body. For a common ambient temperature and constant h, G = hAₛ. More generally G = ∫ₛh dA if the same ambient acts over the surface:

Every term now has units W. The capacity Cth is in J/K, not kilograms: mass m must be multiplied by cₚ. Several independent conductances to fixed reservoirs contribute Gtotal = ΣGₗ and forcing ΣGₗTₗ + Q, giving τ = Cth/Gtotal. Resistances in a series path combine before taking its reciprocal.

3. Check when the approximation is sensible. For convection, choose Lc = V/Aₛ and define Bi = hLc/k. Comparing scales based on that same Lc:

Small Bi means internal equilibration is fast compared with exchange with the surroundings. Bi < 0.1 is a common engineering screening criterion, not an error guarantee for every geometry, local heat source, or material arrangement. The approximation is obtained from the integrated energy balance; setting ∇²T = 0 in the local equation and then discarding the boundary flux would incorrectly eliminate the mechanism that heats or cools the body. MIT: transient heat transfer and the lumped approximation ↗

4. Solve the post-step ordinary differential equation. For constant T∞, Q, and G > 0 after t = 0, set dT/dt = 0 to find Tss = T∞ + Q/G. Let z(t) = T(t)−Tss. Substitution removes the constant forcing:

For z(0) = 0, the solution is already at equilibrium and there is no transition. With a nonzero final change, normalize progress as y(t) = [T(t)−T₀]/[Tss−T₀]. This definition gives y(0) = 0 and y(∞) = 1 for both heating and cooling:

For a body with no heat-loss path (G = 0), constant nonzero Q instead gives T(t) = T₀ + Qt/Cth. There is no finite steady temperature, so the steady-state percentage rise-time formulas do not apply.

C. Derive separate definitions for delay and rise time

Threshold delay tₚ is the elapsed time from the input step to the first crossing y = p. In a monotone response the crossing is unique. A common convention is tdelay = t₅₀, but every report should name its threshold. Solving the lumped response gives:

Rise time is the duration between two selected crossings, not the elapsed time from the input step. For fractions p₁ < p₂:

Elmore / mean delay μ uses the complete response instead of one crossing. For the normalized lumped model, the impulse response is f(t) = dy/dt = e^−t/τ/τ. Its first moment, using integration by parts, is:

Equivalently μ = ∫₀∞[1−y(t)]dt. Thus μ equals τ for this one-pole response, while t₅₀ and t₁₀–₉₀ do not. The standard deviation of this impulse response is also τ, since ∫t²f(t)dt = 2τ² and σ² = 2τ²−τ².

The phrase “rise time to 63.2%” sometimes means an elapsed response time, but it is not the 10–90% rise time used here. A pure exponential never reaches exactly 100% at a finite time; always specify a settling tolerance. For nonmonotone systems, first-crossing times and settling time need additional care because a response can cross a threshold more than once.

D. Show mathematically why a true delay does not lengthen rise time

Suppose the applied input is postponed by a known time d before reaching the same lumped body. This is an imposed transport or actuation delay, not an automatic property of Fourier conduction. Solving the ODE before and after the delayed step gives:

The extra d appears in both thresholds and cancels in their difference. In contrast, the mean delay includes the rectangle of height one and width d before the response starts:

The variance stays the same because shifting a distribution changes its mean but not its width. If the response is known to be a delayed single exponential, its moments imply τ = σ, d = μ−σ (requiring μ ≥ σ), and tr = ln(9)σ. This is a different shape assumption from the gamma fit; it is not a general network identity, and μ−σ must not be interpreted automatically as physical transport delay.

A dimensional example: the same rise time, different delay

Take Cth = 100 J/K and G = 5 W/K, so τ = 20 s. The boundary steps from 20°C to 80°C, with Q = 0. The 10%, 50%, and 90% temperatures are 26°C, 50°C, and 74°C. For cooling from 80°C to 20°C, they would be 74°C, 50°C, and 26°C, with identical normalized timing.

E. Recover the multi-node model from the same heat equation

Divide the body into control volumes. Integrate the energy equation over each volume, approximate its temperature by Tᵢ, and replace each conductive face flux by a conductance times the adjacent temperature difference. For a uniform connecting segment, Gᵢⱼ = kA/ℓ; for unequal cell-center layers, add the two half-layer resistances and any contact resistance before inverting.

Internal fluxes cancel when all node balances are added, which preserves total energy. Kᵢᵢ is the sum of all conductances out of node i, including boundary links, and Kᵢⱼ = −Gᵢⱼ. With constant post-step forcing, solve KTss = b. The error z = T−Tss then satisfies C dz/dt = −Kz.

For a nonzero transition at node j, Aⱼₖ = −aₖφⱼₖ/(Tss,ⱼ−T₀,ⱼ), and ΣₖAⱼₖ = 1. These coefficients depend on the initial state, forcing, and observation point; individual coefficients need not all be positive even for a monotone response. Thresholds must solve ΣₖAⱼₖe^{−tₚ/τₖ} = 1−p, generally numerically. But the remaining-area definition integrates term by term:

This explains why neither the slowest τₖ nor the exact μⱼ automatically equals a threshold delay or rise time. The derivations below obtain these same moments by linear solves, then explicitly choose a one-pole or gamma approximation to estimate rise time.

Diffusion is not a literal waiting period. The classical Fourier heat equation is parabolic. Its ideal mathematical solution generally has an arbitrarily small response at a remote point for any t > 0; it does not create a sharp finite arrival time. A practical “diffusion delay” therefore needs a specified threshold or measurement resolution. Treating L²/α as an exact dead time and then adding a bulk RC rise time can count the same storage twice. Finite-speed or non-Fourier models are outside this page’s scope.

Continue to the exact network-moment derivation ↓

DERIVATION / FROM THE NETWORK TO RISE TIME

Detailed derivations of the two estimates

For a gentler introduction with every step worked out, start with Why the estimates work and its two- and ten-node examples.

Both estimates begin with exact moments of the selected input-to-temperature response. The approximation enters only when we choose a response shape from those moments. Below, C is the diagonal matrix of node heat capacities, K is the grounded conductance matrix, b is the forcing vector for a unit input, and eⱼ selects observation node j.

Start with the heat equation and timing definitions ↑

1. Exact network moments · 2. Elmore-based rise time · 3. Two-moment rise time · 4. Worked 10-node example

1. Derive the exact first and second moments

Step 1 — Transform the heat balance. With input u(t), the thermal network obeys C dT/dt + KT = bu(t). Taking a Laplace transform at zero initial conditions gives:

For a unit step, let x = K⁻¹b be the steady temperature-change vector. Thus Hⱼ(0) = xⱼ. Normalize the transfer function as Fⱼ(s) = Hⱼ(s)/xⱼ, and the unit-step response as yⱼ(t) = Tⱼ(t)/xⱼ. Then Fⱼ(0) = 1 and yⱼ(∞) = 1.

Step 2 — Expand the inverse without changing matrix order. Factor K + sC = K(I + sK⁻¹C), and use the geometric series near s = 0:

Define three successive linear solves. They use the same K, so a matrix factorization can be reused:

Step 3 — Identify the moments. Let fⱼ(t) = dyⱼ/dt be the normalized impulse response. It integrates to one. Expand its Laplace integral directly:

The factor of two is essential: the coefficient of s² is half the second raw moment, not the moment itself. μ has units seconds, m₂ and σ² have units seconds squared. The coefficient wⱼ/xⱼ belongs to the normalized transfer-function expansion; it is not generally a denominator coefficient.

Step 4 — Connect the first moment to step-response delay. Since f = y′ = −(1−y)′, integration by parts gives μ = [−t(1−y)]₀∞ + ∫₀∞(1−y)dt. The boundary term vanishes for a stable finite RC network. Similarly, m₂ = 2∫₀∞t(1−y)dt.

Thus μ measures the total area remaining above the normalized response. It is not defined as a particular crossing time. These integrals also provide a way to estimate moments from a sufficiently long measured or time-marched response, with care for the unobserved tail.

2. Derive the Elmore-based rise-time estimate

Step 1 — Reduce the matrix expression to shared paths. Write Z = K⁻¹. Then:

For a tree with one stepped reservoir and no other sink, the final temperature change is identical at every node, so xᵢ/xⱼ = 1. To interpret Zⱼᵢ, inject one watt at node i while the root reservoir is held at zero temperature change. That watt flows along the unique path from i to the root. The temperature at node j equals the sum of resistances on the part shared by the root-to-i and root-to-j paths. Therefore Zⱼᵢ = Rshared(i,j):

The second form follows by exchanging the sums: each edge carries the contribution of all capacities downstream of it, including branches. With several sinks or local heat input, retain Zⱼᵢ and the steady-state weights xᵢ/xⱼ; an unweighted shared-path sum is not generally valid.

Step 2 — Fit a single pole to the first moment. Choose F₁(s) = 1/(1+sτfit). Its expansion is 1−sτfit+s²τfit²−… . Matching the exact constant and first-order terms forces τfit = μ:

Step 3 — Solve for the thresholds. Set y₁(tₚ) = p, so e^−tₚ/μ = 1−p. Subtract the two threshold times:

This also predicts t₅₀,1 = μ ln(2). Calling μ itself a “delay estimate” does not make μ the exact t₅₀. The rise-time formula is exact for a single lumped RC; for a network it is the consequence of an explicitly chosen one-pole fit.

Why one moment is insufficient. A pure delay d followed by a single RC has F(s) = e^−sd/(1+sτ), μ = d+τ, and σ² = τ². Its actual 10–90% rise time is ln(9)τ because the delay cancels in the difference. The one-moment estimate ln(9)(d+τ) counts that delay as extra rise duration. This example explains the limitation; a diffusion network need not have a literal dead-time segment.

3. Derive the two-moment gamma estimate

Step 1 — Choose a positive two-parameter shape. Mean and variance do not uniquely determine a waveform, so a shape assumption is still needed. This page chooses the gamma impulse-response family:

Here a is dimensionless and θ is a time scale. Substituting z = t/θ verifies ∫f₂dt = 1. Applying the same substitution to ∫tᵏf₂dt gives θᵏΓ(a+k)/Γ(a). Using Γ(a+1) = aΓ(a):

Step 2 — Match the exact network moments. Equating aθ = μ and aθ² = σ², first divide the second equation by the first, then solve for a:

Equivalently, the gamma transfer function is F₂(s) = (1+sθ)^−a. Its expansion is 1−aθs+[a(a+1)θ²/2]s²−… . Substitution yields 1−μs+[(μ²+σ²)/2]s²−… , matching the exact network through s². This is a two-moment fit; it is not necessarily a two-pole model. A noninteger a generally does not correspond to a finite rational RC transfer function.

Step 3 — Integrate to obtain the step response. Define the regularized lower incomplete gamma function P(a,z). Then:

Step 4 — Invert the thresholds. P⁻¹(a,p) denotes the nonnegative z satisfying P(a,z) = p; it is an inverse in the second argument, not 1/P. Consequently:

The inverse is evaluated numerically. For an independent calculation, use gamma.ppf(p, a, scale=theta) in SciPy, with zero location. The SciPy gamma reference defines the distribution and its inverse CDF. Choosing this family is an approximation made here, not a universal Elmore rise-time theorem.

Step 5 — Check the single-node limit. A single RC gives F(s) = 1/(1+sτ), μ = τ, m₂ = 2τ², and σ² = τ². Therefore a = 1 and θ = τ. Since Γ(1) = 1 and P(1,z) = 1−e^−z, the two-moment formula reduces exactly to τ ln(9).

Step 6 — Understand the limits. If a is large, the gamma shape approaches a normal shape with mean μ and standard deviation σ, giving tᵣ ≈ [Φ⁻¹(0.9)−Φ⁻¹(0.1)]σ ≈ 2.563103σ. For a = 1, the factor is ln(9) ≈ 2.197225. For other shapes neither constant is exact. If a < 1, the fitted impulse response is singular at t = 0 even though a finite RC network has a finite initial derivative; this can limit early-time accuracy. Matching two moments does not fix skewness, a separate transport delay, or the long-time decay rate. The estimate can be high or low and has no general error bound.

4. Work through the 10-node example numerically

Use the default left-step / insulated-end ladder: N = 10, each r = 0.1 K/W, each c = 0.1 J/K, and observe node 10. A unit boundary step gives xᵢ = 1 K. Here Zⱼᵢ = r min(i,j), so the complete first-moment vector is:

At the last node, the first solve gives μ₁₀ = 0.01(1+2+…+10) = 0.55 s. For the second solve, use vᵢ = xᵢμᵢ and sum the same inverse-matrix row again:

One moment gives tᵣ,1 = ln(9)(0.55) = 1.208474 s. Two moments give a = 0.55²/0.2035 = 1.486486 and θ = 0.2035/0.55 = 0.370000 s. The gamma quantiles are:

The full modal solution gives t₁₀ = 0.141062 s and t₉₀ = 1.136826 s, hence tᵣ = 0.995764 s. Refined RK4 time marching agrees to the displayed six decimals. The one-moment estimate is about 21.36% high; the two-moment estimate is about 4.71% high. Notice that a reasonably close rise-time difference does not imply that both individual threshold times are equally accurate.

Compare the other examples and time-step refinement below ↓ · Inspect the selected node’s individual moment contributions ↑

PHYSICAL LAYERS → THERMAL CIRCUIT

Convective boundaries & non-uniform layers

Consider air heating a copper plate, followed by an insulating polymer and a sensor. The far face of the sensor is insulated. Each node stores the heat capacity of its entire layer, while its temperature represents the layer center. Assume one-dimensional conduction, common cross-sectional area A, constant properties, and perfect contact between layers.

1. Place the nodes at the layer centers

Fourier’s law gives Q̇ = kA(Tleft−Tright)/L, hence R = ΔT/Q̇ = L/(kA). A uniform half-layer has half that resistance. Newton’s convection law gives Q̇ = hA(Tamb−Tsurface), hence Rconv = 1/(hA).

2. Convert the geometry into actual circuit links

Add an interface contact resistance to its corresponding inter-node link if needed. A right convective boundary would add r₃b = R₃/2 + 1/(hright Aright); it would change the topology and require the general matrix moment formula instead of the single-sink sum below.

3. Accumulate the resistance back to ambient

Each complete intervening layer contributes its full resistance; the observed layer contributes only the half between its center and the previous interface. These paths describe a lumped approximation to the geometry, not an exact spatial temperature field inside each layer.

4. Derive the Elmore delay at the sensor

Apply a step in Tamb with no internal heat generation and no other sink. At steady state, every node has the same temperature change, so xᵢ/x₃ = 1. The common path from node i and node 3 back to ambient has resistance Sᵢ. The general moment formula therefore reduces to:

This expansion shows how convection loads all capacities, whereas each downstream resistance affects a smaller subset. Equivalently μ₃ = r₀₁(C₁+C₂+C₃) + r₁₂(C₂+C₃) + r₂₃C₃. At an upstream observation node j, the correct shared resistance is Smin(i,j), so μⱼ = Σᵢ Smin(i,j)Cᵢ.

5. Follow a numerical example

For continuity with the existing unequal-layer preset, use illustrative full-layer resistances R = [0.005, 2.5, 0.1] K/W, capacities C = [3.5, 0.8, 0.2] J/K, and Rconv = 0.5 K/W. These are circuit inputs for an example, not certified copper, polymer, or sensor material data.

For the same input and sensor, the two-moment estimate is about 6.71640 s; the full modal solution and RK4 time marching give about 6.47804 s. The one-moment estimate is about 28.00% high and the two-moment estimate about 3.68% high. The explorer shows the full transient, all node times, and individual moment contributions.