Proving Physics: Encoding Differential Equations in ZK

A differential equation describes how something changes. A zero-knowledge proof says "this computation happened correctly." Put them together, and you can prove that a physical simulation ran faithfully — without revealing what was simulated.

Here's how we encoded the Rayleigh-Plesset ODE into a halo2 circuit, step by step.

The Physics in 30 Seconds

A gas bubble in liquid oscillates under pressure. The governing equation (Rayleigh-Plesset) is:

R'' = [P_gas - P₀ + 2σ/R - 4μṘ/R] / (ρR) - 1.5(Ṙ²/R)

Where R is the bubble radius, P_gas is internal gas pressure, P₀ is ambient pressure, σ is surface tension, μ is viscosity, ρ is liquid density, and Ṙ and R'' are the velocity and acceleration of the bubble wall.

Under the right conditions, the bubble collapses from 5 μm to 0.77 μm and reaches 12,348 K. That's sonoluminescence — light from sound.

Our job: prove this simulation ran correctly in zero knowledge.

Circuit Architecture

The circuit has three layers:

SonoluminescenceCircuit
├── FieldArithChip (scaled multiply, divide, add, sub)
│ ├── Gate: s_mul → a × b = c × S
│ ├── Gate: s_div → c × b = a × S
│ ├── Gate: s_add → a + b = c
│ └── Gate: s_sub → a - b = c
└── PhysicsStep (chained N times)
├── Velocity: Ṙ = (R_curr - R_prev) / dt
├── Gas pressure: P_gas = P_initial × (R₀/R)^5
├── Temperature: T = T₀ × (R₀/R)²
├── RP acceleration: R'' = δP/(ρR) - 1.5Ṙ²/R
└── Verlet update: R_next = 2R - R_prev + R''×dt²

Every arithmetic operation is a constrained region using our FieldArithChip. The physics step chains ~30 of these regions together.

One Timestep: 30 Constrained Regions

Here's what happens in a single physics_step() call. Every line is a constraint that the verifier checks:

1. Velocity (2 regions)

r_diff = R_curr - R_prev // sub gate
rdot = r_diff × S / dt // div gate

Störmer-Verlet derives velocity from positions. No separate velocity variable is stored.

2. Gas pressure (5 regions)

ratio = R₀ × S / R_curr // div: R₀/R
ratio2 = ratio × ratio / S // mul: (R₀/R)²
ratio3 = ratio2 × ratio / S // mul: (R₀/R)³
ratio4 = ratio3 × ratio / S // mul: (R₀/R)⁴
ratio5 = ratio4 × ratio / S // mul: (R₀/R)⁵
p_gas = P_initial × ratio5 / S // mul: P_gas

Six constrained operations for gas pressure. Each intermediate result is its own verified cell.

3. Temperature (2 regions)

ratio2 = (already computed above)
temperature = T₀ × ratio2 / S // mul

We reuse ratio2 from the pressure computation — the circuit shares intermediate values through cell references.

4. Surface tension and viscosity (5 regions)

two_sigma_R = 2σ × S / R_curr // div: 2σ/R
four_mu_rdot = 4μ × rdot / S // mul: 4μṘ
viscous_term = four_mu_rdot × S / R_curr // div: 4μṘ/R

5. RP acceleration (4 regions)

delta_p = p_gas - P₀ + two_sigma_R - viscous_term  // add/sub chain
rho_R = ρ × R_curr / S             // mul: ρR
term1 = delta_p × S / rho_R        // div: δP/(ρR)
rdot_sq = rdot × rdot / S          // mul: Ṙ²
term2 = 1.5 × rdot_sq × S / R_curr // mul + div: 1.5Ṙ²/R
rddot = term1 - term2              // sub: R''

6. Verlet update (3 regions)

rddot_dt2 = rddot × dt² / S       // mul: R''×dt²
two_r = 2 × R_curr / S             // mul: 2R
r_next = two_r - R_prev + rddot_dt2 // add + sub

Total: ~30 constrained regions per timestep. For a 3,000-step simulation, that's 90,000 independently verified operations.

The PhysicsStep Function

In code, the step function signature reveals the structure:

pub fn physics_step<F: PrimeField>(
chip: &FieldArithChip<F>,
mut layouter: impl Layouter<F>,
params: &PhysicsParams<F>,
state: &StepState<F>,
sin_i: &AssignedCell<F, F>,
) -> Result<StepResult<F>, Error> {
// 1. Velocity
let r_diff = chip.sub(
layouter.namespace(|| "R_curr - R_prev"),
&state.r_curr,
&state.r_prev,
)?;
let rdot = chip.scaled_div(
layouter.namespace(|| "rdot = r_diff * S / dt"),
&r_diff,
&params.dt,
)?;
// ... 28 more constrained operations
Ok(StepResult { r_next, r_curr, rdot, p_gas, temperature, rddot, emission })
}

PhysicsParams holds all physical constants as assigned cells (loaded once via load_constant with Fixed column constraints). StepState holds R_curr and R_prev. The function returns StepResult with everything the next step needs.

Chaining Steps: The Full Circuit

The outer circuit loops N times, feeding each step's output into the next step's input:

Step 0: (R₀, R_start) → (R₁, R₀, T₀, ...)
Step 1: (R₁, R₀) → (R₂, R₁, T₁, ...)
Step 2: (R₂, R₁) → (R₃, R₂, T₂, ...)
...
Step N-1: (R_{N-1}, R_{N-2}) → (R_N, R_{N-1}, T_{N-1}, ...)

The chaining is enforced by halo2's copy constraints — the r_next cell from step n is constrained to equal the r_curr cell of step n+1. The prover can't skip steps or reorder them.

Three Public Inputs

Out of the entire simulation — thousands of intermediate values, all the physical parameters — only three values are exposed:

• R₀ — the equilibrium bubble radius (what physical system is being claimed)

• Final temperature — the result (did the bubble get hot enough?)

• Total emission count — how many steps exceeded the sonoluminescence threshold

These are assigned to the instance column. Everything else — pressure, viscosity, surface tension, every intermediate radius and velocity — stays in advice columns, hidden by the zero-knowledge property.

Circuit Sizing

The number of rows needed scales linearly with steps:

k = min k such that 2^k ≥ n_steps × 30 + 30

Steps
k
Rows (2^k)

10
9
512

100
13
8,192

500
15
32,768

1,000
16
65,536

3,000
18
262,144

At k=18, the circuit has 262,144 rows — and generates a proof of exactly 1,088 bytes. The same size as the 10-step proof.

Why This Structure Works

Three properties make this circuit design viable:

Linearity: Each step uses exactly the same gate configuration. The circuit "template" is defined once; it's instantiated N times.

Locality: Each step only references the current and previous radius. There's no global state, no random access into the simulation history. This makes the constraint system sparse and efficient.

Composability: The FieldArithChip's primitive operations (mul, div, add, sub) compose into arbitrarily complex physics. Adding new terms to the ODE — acoustic driving, different gas models — just means more chip method calls per step.

The constraint system doesn't care that we're simulating physics. It just sees polynomial identities over field elements. The physics is in the structure of which identities we check and how they chain together.

This is Part 3 of our 8-part series. Part 2: Building Custom halo2 Chips covers the constraint system in detail. Up next: Part 4: 1,088 Bytes to Prove a Star — why proof size is constant and what that means.

New here? Subscribe free for the full series.

---

Originally published at https://jacobian.ghost.io/proving-physics-encoding-differential-equations-in-zk/.

Subscribe to **Jacobian* for weekly ZK/privacy insights: jacobian.ghost.io*

New here? Grab the free **ZK Privacy Patterns* guide: Get it free*