Proving Sonoluminescence in Zero Knowledge: How We Built a 1,664-Byte Proof for Real Physics

A Bubble, a Collapse, and 1,664 Bytes

A gas bubble in water, driven by a 26.5 kHz acoustic field, collapses so violently that it reaches 12,348 Kelvin -- hotter than the surface of the sun. For roughly 100 picoseconds, it emits light. This phenomenon, sonoluminescence, was discovered in 1934 and remains one of the most extreme energy-focusing mechanisms in nature.

We built a zero-knowledge proof that verifies the entire simulation of this collapse -- 2,359 time steps of a nonlinear differential equation, gas pressure calculations, temperature evolution, and Stefan-Boltzmann radiation -- in a constant-size proof of 1,664 bytes. The verifier learns three numbers: the bubble's equilibrium radius, its final temperature, and total radiated emission. Everything else -- surface tension, viscosity, pressure traces, intermediate states -- stays hidden.

The proof system is halo2 PLONK with KZG commitments on bn256. The codebase is roughly 3,500 lines of Rust. And the hardest part of building it had nothing to do with cryptography.

The hardest part was arithmetic.

The Gap Between Physics and Fields

The Rayleigh-Plesset equation governs how a bubble wall moves in a liquid:

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

Every term in this equation involves real numbers. The equilibrium radius R₀ is 5 micrometers -- 5 × 10⁻⁶ meters. Surface tension σ is 0.0728 N/m. The time step dt is 1 nanosecond. During collapse, the bubble radius drops from 25 μm to 0.77 μm, gas pressure swings by a factor of 200×, and temperature varies by 170×.

A ZK circuit operates in a completely different world. Every value is an element of a prime field F_p, where p is the bn256 scalar field modulus -- roughly 2²⁵⁵. There are no decimals. There are no fractions. There is no floating point. Multiplication doesn't overflow; it wraps modulo p. Division isn't truncation; it's multiplication by a modular inverse.

These two computational models -- real-number physics and modular field arithmetic -- are fundamentally incompatible. Bridging them is the entire engineering challenge of zk-physics, and the technique we used generalizes to any scientific computation you might want to prove in zero knowledge.

Scaling Into the Field: Why 10³⁰

The standard approach to representing real numbers in a finite field is fixed-point scaling: pick a large integer S, multiply every physical value by S, and work with the resulting integers. Your "real number" 0.0728 becomes the integer 72,800,000,000,000,000,000,000,000,000 -- which is 0.0728 × 10³⁰.

The choice of S = 10³⁰ isn't arbitrary. It's constrained from both sides.

From below: the smallest meaningful physical value in the simulation is dt² = 10⁻¹⁸ seconds². Scaled by 10³⁰, that becomes 10¹². Any smaller scaling factor and this value rounds to zero, killing the integration. The dynamic range of the simulation -- 200× in pressure, 170× in temperature, plus intermediate products -- demands at least 10²⁸ of scaling to preserve enough significant digits through chained multiplications.

From above: when you multiply two scaled values, the product is on the order of S² = 10⁶⁰. This product must fit in the field without wrapping. The bn256 field modulus p is approximately 2²⁵⁵, which is roughly 5.8 × 10⁷⁶. That leaves 10⁷⁶ / 10⁶⁰ = 10¹⁶ -- sixteen orders of magnitude of headroom. Enough to accumulate error across thousands of steps without any product accidentally wrapping around the field modulus and producing garbage.

You can see the scaling factor and the reasoning behind it in src/field_utils.rs:

// src/field_utils.rs:10
/// Scaling factor: 10^30
/// Chosen to accommodate sonoluminescence's 200x pressure / 170x temperature
/// variation while keeping products (10^60) well below the ~2^255 field modulus.
pub const SCALE: u128 = 1_000_000_000_000_000_000_000_000_000_000; // 10^30

This single constant governs the entire system. Every physical value, every intermediate computation, every constraint in the circuit is denominated in units of 10⁻³⁰.

Multiplication That Proves Itself

Here's where the circuit design diverges sharply from normal programming. In Rust, you multiply two numbers and get a result. In a ZK circuit, you must constrain the relationship between inputs and outputs so that a cheating prover can't substitute a false result.

For scaled multiplication, the naive constraint would be: given scaled values a and b, the product c satisfies a × b = c × S. But this is a field equation -- it holds modulo p. A cheating prover could find a c' such that a × b = c' × S + k × p for some integer k, and the constraint would pass. The result would be mathematically valid in the field but physically meaningless.

The fix is to constrain the remainder explicitly. The scaled multiplication gate in src/chips/field_arith.rs enforces:

a × b - c × S - r_mul = 0  (mod p),  where 0 ≤ r_mul < S

The prover must supply both the quotient c and the remainder rmul, and the circuit range-checks rmul to ensure it's less than S. This transforms a modular equation into an integer equation: the prover is proving that a × b = c × S + rmul over the integers, not just mod p. The range check on rmul is what makes this work -- without it, the constraint has p/S ≈ 10⁴⁶ valid solutions for any given a and b. With it, there's exactly one.

Division uses the same principle, inverted:

a × S - c × b - r_div = 0  (mod p),  where 0 ≤ r_div < b

This is not modular inverse division. It's integer division with a proven remainder. The distinction matters: modular inverse division gives you an exact field element, but that element has no relationship to the truncated integer quotient that physics requires. When you divide a pressure by a radius to get an acceleration, you need truncation semantics, not field semantics.

Every multiplication and division in the physics simulation -- and there are dozens per time step -- carries this quotient-remainder constraint. The circuit doesn't trust the witness values; it algebraically proves that each arithmetic operation was performed correctly under integer semantics.

Witness Generation: The Other Half of the Problem

The circuit constrains relationships between values. The witness generator computes the actual values. These two computations must agree exactly -- not approximately, not to 15 significant digits, but bit-for-bit in every intermediate result.

This is harder than it sounds. The witness generator runs in Rust, outside the circuit. The obvious implementation would use f64 arithmetic:

// DON'T DO THIS
fn to_scaled(v: f64) -> i128 {
    (v * 1e30) as i128
}

An f64 has a 53-bit mantissa -- about 16 decimal digits of precision. When you multiply 0.0728 (surface tension) by 10³⁰, the exact result is 72,800,000,000,000,000,000,000,000,000. The f64 representation gives you 72,800,000,000,000,004,194,304 -- correct to 16 digits, then garbage. That discrepancy, roughly 10⁻¹⁶ in relative terms, is negligible for any numerical simulation. For a ZK proof, it's fatal. The circuit computes the exact integer product using field arithmetic. The witness supplies the approximate f64 product. They disagree. The proof fails.

The witness generator in zk-physics uses num-bigint for arbitrary-precision integer arithmetic, avoiding floating-point entirely for all intermediate computations. The scaling function in src/field_utils.rs uses a 3-way split to handle large products:

// src/field_utils.rs:56
/// Split factor for wide multiplication. 10^10, so SPLIT^3 = 10^30 = SCALE.
const SPLIT: i128 = 10_000_000_000; // 10^10

pub fn scaled_mul(a: i128, b: i128) -> i128 {
    let sign = if (a < 0) != (b < 0) { -1i128 } else { 1 };
    let a_abs = a.unsigned_abs();
    let b_abs = b.unsigned_abs();
    let sp = SPLIT as u128;

    let a0 = a_abs % sp;
    let a1 = (a_abs / sp) % sp;
    let a2 = a_abs / (sp * sp);

    let b0 = b_abs % sp;
    let b1 = (b_abs / sp) % sp;
    let b2 = b_abs / (sp * sp);

    // Products with i+j >= 2 contribute to the integer part after dividing by SCALE
    let sum_2 = a0 * b2 + a1 * b1 + a2 * b0;
    let sum_3 = a1 * b2 + a2 * b1;
    let sum_4 = a2 * b2;

    let mut result: u128 = sum_2 / sp;
    result = result.wrapping_add(sum_3);
    if let Some(s4_scaled) = sum_4.checked_mul(sp) {
        result = result.wrapping_add(s4_scaled);
    } else {
        result = result.wrapping_add(sum_4.wrapping_mul(sp));
    }
    sign * result as i128
}

Each operand is decomposed into three 10-digit limbs (a = a₂ × 10²⁰ + a₁ × 10¹⁰ + a₀). The partial products are accumulated by their power of SPLIT, and the final sum is divided by SCALE = SPLIT³. The maximum intermediate value in any partial product is (10¹⁰ - 1)² ≈ 10²⁰, which fits comfortably in u128 (max ~3.4 × 10³⁸).

Division is even trickier. Since a × SCALE can exceed u128 range, the witness generator uses iterative long division with a chunk factor of 10⁶, performing five rounds (since (10⁶)⁵ = 10³⁰ = SCALE):

// src/field_utils.rs:100
pub fn scaled_div(a: i128, b: i128) -> i128 {
    // ...
    const CHUNK: u128 = 1_000_000;
    const ROUNDS: usize = 5;

    let mut acc: u128 = a_abs / b_abs;
    let mut rem: u128 = a_abs % b_abs;

    for _ in 0..ROUNDS {
        let wide = rem * CHUNK;
        let digit = wide / b_abs;
        rem = wide % b_abs;
        acc = acc * CHUNK + digit;
    }

    sign * acc as i128
}

Each round multiplies the running remainder by 10⁶, divides by b, and appends the quotient digit. The maximum intermediate product is rem × CHUNK < b × 10⁶. For physical values up to ~200 × SCALE ≈ 10³², this gives ~10³² × 10⁶ = 10³⁸ -- just within u128 range.

These two functions -- scaledmul and scaleddiv -- are the foundation of the entire system. Every physical computation in both the witness generator and the circuit reduces to sequences of these operations.

Chaining Physics Steps Through the Circuit

With reliable fixed-point arithmetic in place, encoding the Rayleigh-Plesset equation becomes mechanical -- tedious, but straightforward. Each physics step in src/chips/physics_step.rs chains together roughly 20 constrained operations:

1. Velocity from the Störmer-Verlet scheme: Ṙ = (Rcurr - Rprev) × S / dt
2. Gas pressure via the polytropic relation: Pgas = Pinitial × (R₀/R)⁵
3. Temperature from adiabatic compression: T = T₀ × (R₀/R)²
4. Net pressure combining gas pressure, ambient pressure, surface tension, and viscous damping
5. Acceleration from the full Rayleigh-Plesset equation
6. Next radius via Verlet update: Rnext = 2Rcurr - R_prev + R'' × dt²
7. Emission via Stefan-Boltzmann: E = σ_SB × 4π × R² × T⁴ × dt (accumulated only when T > 5,000 K)

The choice of Störmer-Verlet integration over forward Euler is deliberate. Verlet is a symplectic integrator -- it conserves energy over long simulations, even with large time steps. The bubble collapse involves a 30× change in radius over ~2,359 steps, with enormous accelerations near minimum radius. Forward Euler would accumulate energy drift that corrupts the physics. Verlet's update formula, R{n+1} = 2Rn - R{n-1} + R''n × dt², costs the same number of constraints as Euler (one multiply, two additions) but stays stable through the violent collapse phase.

These steps are chained N times in src/circuits/sonoluminescence.rs, with each step's output radius feeding into the next step's input. The circuit has k = 18, meaning 2¹⁸ = 262,144 rows of constraints.

What the Verifier Sees

The proof exposes exactly three public inputs:

1. R₀ -- the equilibrium bubble radius (5 μm, scaled)
2. Final temperature at the last simulation step
3. Cumulative radiated emission -- the total Stefan-Boltzmann energy output

Everything else is private: the 2,359 intermediate radius values, all pressure traces, velocity history, surface tension, viscosity, ambient pressure, and the acoustic driving signal. The verifier can confirm that a valid sonoluminescence simulation was executed, producing the claimed temperature and emission, without learning anything about the physical parameters that produced it.

Verification takes 2.56 milliseconds on Apple Silicon. The proof is 1,664 bytes -- constant regardless of how many time steps the simulation runs. This is a property of KZG-based polynomial commitments: proof size depends on the commitment scheme, not the circuit size. A 10-step simulation and a 10,000-step simulation produce identically sized proofs.

Why This Matters Beyond Bubbles

The zk-physics codebase proves something broader than sonoluminescence: that ZK circuits can encode real scientific computation. Not hash preimages, not Sudoku puzzles, not financial transactions -- actual differential equations with physical units, extreme dynamic range, and nonlinear dynamics.

The fixed-point arithmetic technique -- scaled integers with explicit quotient-remainder constraints -- generalizes to any computation where you need real-number semantics inside a prime field. Climate models, orbital mechanics, fluid dynamics, structural engineering: any simulation that runs as a numerical ODE integrator can, in principle, be constrained the same way.

The practical barriers are constraint count (262,144 rows for ~3,000 physics steps) and witness generation complexity (every intermediate value must be computed at arbitrary precision). But the proof size stays constant, and verification stays in the low milliseconds. For applications where you need to prove a simulation was executed honestly -- regulatory compliance, scientific reproducibility, competitive modeling -- the 1,664-byte proof is a compelling primitive.

The complete eight-part series on Jacobian walks through every layer of this system, from the field arithmetic gates to the privacy tests that verify information leakage. The deep-dives on custom halo2 chip design and KZG proof mechanics are available to Pro subscribers -- if you're building circuits that go beyond toy examples, they'll save you weeks of debugging the exact witness-circuit mismatches we already found and fixed.

The Constraint That Haunts You

There's a pattern that emerges when you build systems like this. The cryptography works. The proof system works. The polynomial commitments and KZG pairing checks are battle-tested by teams with far more resources than ours. What breaks -- every time, on every project -- is the boundary between the witness and the circuit. The moment your Rust code computes a slightly different integer than your constraint system expects, the proof is unprovable. Not wrong. Not insecure. Unprovable.

The entire 707-line field_arith.rs chip exists to make that boundary airtight. Every multiply, every divide, every range check is there because we hit a case where the witness and the circuit disagreed by one part in 10¹⁶ -- and one part in 10¹⁶ is enough to make the prover fail silently.

If you're encoding real-world computation into a ZK circuit, the question isn't whether you'll hit this problem. The question is whether you'll find it before your users do.