Interactive companion · numerical experiments

Where to put the traps

Optimizing the mean first-passage time with many small traps in heterogeneous media — D. S. Grebenkov & T. Kolokolnikov

A diffusing particle starts from a distribution ω(x) in a medium of diffusivity D(x). You have N absorbing traps of radius ε. Where do they go? Two answers compete, and which one is right depends on a single number — the total trapping strength κ = νN, where ν = 2π ⁄ log(1/ε) is what one small trap is worth in two dimensions. Note how weakly ν depends on ε: that logarithm is the whole story.

Energy minimisers

Direct minimisation over the N trap positions, live in a rectangle (0,L)×(0,H) — Fig. 3's own domain — or in the unit disk, of the paper's pairwise energy (2.27) — whose minimiser is ½(α+ω) for every κ — or of the full linear system (2.22) solved to all orders in ν, where κ does enter, through ν = κ/N. The heterogeneity is a tanh step — across the width in the rectangle, radial in the disk — placed either in the medium D(x) or, with D constant, in the release density ω(x). In full-system mode each step solves the bordered symmetric system for (C₁…C_N, T) and differentiates ū_ω = T − ΣC_jV_ω(x_j) through it by an adjoint. Here ν = κ/N, so the ε implied by your κ and N is shown — watch it, because the whole construction assumes ε ≪ 1.

Background: whichever profile carries the step — D(x), or ω(x) in the constant-A mode. The high-contrast region is where particles spend their time. Cells: Voronoi partition used to estimate the local density as ρ ≈ 1/(N·S_j), exactly as in Fig. 1 of the paper. Areas come from a 560²-equivalent raster (reshaped to the domain's aspect) with an exact per-cell candidate set, renormalised so ΣS_j = |Ω|.

Emergent density vs. the prediction

Vertical axis: trap density ρ. Voronoi estimate at each trap, against the optimum at this κ and both asymptotic limits — Eq. (1.5) and Eq. (1.6).

Voronoi estimate optimal ρ at this κ weak limit ½(α+ω) strong limit ∝√(ωα)

Points are single traps; the scatter is finite-N noise, not disagreement — cell areas are accurate to 0.17% rms against an 1800² reference, so the spread you see is real. The estimator is biased low in the outermost ring, where cells are clipped by ∂Ω — a boundary layer the continuum limit does not resolve.

Flow

traps N100

Model being minimised

running

Convergence

steps0

max |∂E/∂x_k|—

energy E_κ/N²—

L¹ to weak ½(α+ω)—

L¹ to κ-optimal—

ν = κ/N—

implied ε = e^−2π/ν—

The paper stops at max|x′_k| < 10⁻⁵. Reseeding lands on a different configuration every time — the measured density barely moves. Every step is a backtracking line search, so E_κ decreases monotonically. At κ ≪ 1, Φ → ½(U+V_ω) and E_κ is exactly the paper's energy (2.25).

Energy descent

max |∂E_κ/∂x_k| against descent step, log scale

Plateaus are the descent crossing between rearrangements of the outer rings — the discrete problem has many local minima, and the plateau structure is where they trade places.

The κ dial

The optimality system for arbitrary κ — (−Δ+κρ)u = A, (−Δ+κρ)p = ω, up = const — solved on a radial grid and minimised over ρ. One slider takes you from the weak-trapping answer to the strong-trapping one, through everything in between that has no closed form.

weak ½(α+ω) strong ∝√(ωα) optimal at this κ

Vertical axis: trap density ρ. The dotted vertical line is the step radius r₀.

Mean capture time u(r)

Vertical axis: u, the expected time to capture from radius r.

under the optimal ρ under uniform ρ |A|/κ

With ω = a the optimum makes u exactly constant at |A|/κ — capture time stops depending on where the particle started. That is the paper's main result, and it is the only case where flattening is the answer: elsewhere the optimum minimises ∫ωu, which is a different thing. With A constant it is uniform placement that gives the flat profile, and the optimum deliberately tilts away from it.

Acknowledgements

This page was written by Claude (Opus 5) — the interface, the numerical methods, and the prose alike — in conversation with the authors, as a companion to the paper. It is not a substitute for the paper; where the two disagree, the paper is authoritative.

Everything on it was checked against something independent. The optimality-system solver reproduces both of the paper's published tables — Ω = (0,1), α = 1, ω = 2·1_x<1/2 to three decimals, and the (0,3)×(0,1) rectangle of Fig. 3 to four decimals at every κ from 0.1 to 5000. The rectangle's Neumann Green's function, which the paper cites rather than prints, is derived here and checked against an independent finite-difference solve of the same boundary-value problem. The quadrature for V_α and V_ω matches the closed form for a quadratic profile to 1e−6. The adjoint gradient of the full system (2.22) agrees with central differences to 3.8 × 10⁻⁵. Voronoi cell areas sit within 0.17% rms of an 1800² reference. Where a question could not be settled at the sizes a browser can reach — whether (2.22) converges to the homogenised optimum as ν grows — the page says so rather than implying an answer.