Critical review of the M.Sc. thesis and its code

Purpose: record what the thesis established, what it left open, and what in its method or code must change before the results can go into a paper. Every finding cites the source so it can be checked. Paths are relative to the project root (related/ holds the symlinks).

Sources reviewed in full: thesis TeX (related/msc-thesis/src/chapters/*.tex), thesis Python (related/msc-thesis.py/src), Dullemond’s coag_py (related/msc-thesis.py/lib/coag_py), Dullemond’s notes Dust evolution with binning methods (19 pp.).

1. What the thesis did

flowchart LR
    A[Disk model<br/>irradiated power law,<br/>r0 = 10 au midplane] --> B[Dust kinematics<br/>BR, RD, AZ, TU]
    B --> C[Collision rates<br/>sigma dv]
    C --> D[Outcome probabilities<br/>MB-smoothed v_frag = 1 m/s,<br/>no bouncing]
    D --> E[Kernel K_kij<br/>Kovetz-Olund + cancellation,<br/>MRN fragmentation]
    E --> F[Integration<br/>scipy Radau per macro step]
    E --> G[Sampled kernel<br/>pairs without replacement]
    G --> F
    F --> H[Diagnostics<br/>mass error, accuracy error,<br/>RMS dN/dt, sparsity]

2. What is solid and carries over

1. The kernel construction (Kovetz and Olund plus cancellation handling) matches coag_py.create_coag_kernel and conserves mass per collision to about 10⁻¹². It becomes the reference for M4 parity tests.
2. The physical model is a reasonable, documented “standard model” of local dust evolution at one radius (close to Birnstiel 2010 and DustPy). Good enough as the testbed for a numerical-methods paper.
3. The diagnostics (mass error, RMS dN/dt, sparsity, sampled-pair maps) are the right ones to keep; some need better definitions (F9).
4. The main conclusion (sampling pays off when the particle state has more than one property) is plausible and is exactly the argument Dullemond’s notes build toward. It has not been tested yet. Testing it is the core of the paper.

3. Findings that change the method

Severity: A invalidates a quantitative result, B biases or weakens a result, C presentation or citation.

F1 (A): the sampling estimator is biased low; this explains the “slower evolution”

Thesis ch. 4, “Sampling of the Total Kernel”, Observation 3, attributes the slower approach to equilibrium at low ρ_sample to “effectively lower collision rates”. That is a direct consequence of the estimator, not a property of sampling in general:

• Pairs are drawn without replacement (np.random.choice(..., replace=False), related/msc-thesis.py/src/models/kernel/sampled.py:94-95, config allow_duplicate_sampling = false).
• The sampled kernel is the full kernel restricted to the drawn pairs, with no reweighting (sampled.py:69-70, Kernel(..., ijs=ijs)).
• Hence E[F̂_k] = Σ_ij π_ij K_kij N_i N_j with inclusion probability π_ij ≤ 1. Every term is scaled down by its inclusion probability, so the rate of every process is underestimated. A Horvitz-Thompson weight 1/π_ij, or importance sampling with replacement and weight 1/(S p_ij), removes the bias.
• The with-replacement branch (sampled.py:71-72) divides by P_ij but not by S, and it is applied to a kernel built from the set of unique pairs, so multiplicities are lost. It is not an unbiased estimator either, and the thesis runs did not use it.
• Dullemond’s notes (Eq. 1.70) and his 2021 prototype (coag_py/mc_sampling/coag.py, sample_kernel_sum) already use the unbiased form K/(F_i F_j) with draws with replacement. The thesis diverged from it.

Consequence: the accuracy error in thesis Fig. 106 mixes bias (systematic lag) and variance (noise). Only the variance part is intrinsic to sampling. The headline negative result of the thesis is partly an artifact and should be re-measured with an unbiased estimator. This is the strongest single argument for a follow-up paper.

F2 (A): the documented sampling probability is not what the runs used

The thesis defines P_ij ∝ ΔM^net_ij N_i N_j m_i m_j with ΔM^net_ij = Σ_k m_k |K_kij|. In the code path used for all sampling figures (examples 104, 106, 107 import SolverV2):

• SolverV2.__init__ sets self.W_ij = weights(self.mg, self.kernels[0]) (related/msc-thesis.py/src/models/solver/solver_v2.py:91).
• weights() (solver_v2.py:223-226) calls np.sum([...]) without axis=0, which reduces to a scalar.
• SampledKernel uses this scalar (sampled.py:39-44), so P_ij ∝ N_i N_j m_i m_j. The ΔM^net factor is a constant and drops out after normalization.

So the runs sampled by N_i N_j m_i m_j, i.e. by the product of bin masses ρ_i ρ_j, not by the documented relevance measure. The per-pair throughput (which vanishes for pairs that are cut at the grid edge, and is tiny where Δv is small) never entered. This may explain the “peaks at the upper edge” (Observation 4), since pairs near m_max carry most mass but partially have zero kernel.

F3 (B): the documented proposal is not the variance-optimal one

Even as documented, P_ij ∝ ΔM^net N_i N_j m_i m_j has an extra m_i m_j factor justified only heuristically. For the unbiased IS estimator of the vector F = Σ_ij K_·ij N_i N_j, the proposal that minimizes E‖F̂ − F‖ in a weighted norm ‖·‖w is p_ij ∝ ‖K·ij‖_w N_i N_j (standard result for vector-valued IS; derivation in 04-experiments.md, E7). With the mass-weighted L1 norm this is exactly p_ij ∝ ΔM^net_ij N_i N_j, i.e. the thesis’ relevance measure without m_i m_j. The m_i m_j factor double counts mass and oversamples big-big pairs.

F4 (B): one frozen sample per macro step, with long steps

Each of the 200 macro steps calls scipy Radau once over [t_n, t_{n+1}] with a kernel sampled once at t_n (solver_v2.py:126-155, N_subst = 1). On a log time grid the late macro steps span about 10⁸ s, so a single random draw drives the solution for a long interval and MC noise does not average out within that interval. Noise averaging needs either more, shorter steps, or several independent samples per step (averaged), or resampling per integrator step. This is a design variable to study (E9), not a fixed choice.

F5 (B): integrator tolerances and order

• scipy Radau is the 3-stage Radau IIA method of order 5 (also stated in the react_0d.py docstring). The thesis abstract, ch. 3, and ch. 5 say “4th order”. Fix in any write-up.
• solve_react_0d_equation passes no rtol/atol, so scipy defaults apply (rtol = 10⁻³, atol = 10⁻⁶). N_k is in m⁻³ and spans many decades (rough estimate: about 3·10² m⁻³ for monomers at 10 au, many orders less for the largest bins). An absolute tolerance of 10⁻⁶ m⁻³ means bins below that are not error-controlled. Needs a scaled state or per-component atol (ADR 0003). To verify in M6.
• The professor’s react_0d.py is a wrapper around scipy Radau with an analytic Jacobian plus a Newton-iterated backward Euler. It is not a custom Radau implementation. Our own Radau IIA (ADR 0008) is new work.

F6 (B): molecular viscosity in the radial gas velocity

The thesis (ch. 2, “Gas Viscosity” and “Radial Gas Velocity”) and the code (related/msc-thesis.py/src/models/disk/disk.py:54,72, physics.viscosity(u_th, lambda_mfp)) use the molecular viscosity ν = ½ u_th λ_mfp in the Lynden-Bell and Pringle formula for u_g. Birnstiel et al. (2010) use the turbulent viscosity ν_t = α c_s H_p. ν_mol is smaller than ν_t by the Reynolds number (a rough estimate for this model gives Re ≈ 10⁵ at 10 au), so u_g is underestimated by about five orders of magnitude. Impact on Δv is small (the drag term of radial drift scales with u_g and differs between grains only through St), but it is physically wrong and must be fixed in M2 with a toggle for both choices so the effect can be shown.

F7 (B): negative densities are clipped and the mass is redistributed silently

fix_negative_densities (related/msc-thesis.py/src/functions/solver/__init__.py, called at solver_v2.py:157) zeroes negative bins and rescales all positive bins to restore total mass. This conserves mass by construction, so the reported mass error after such a fix does not measure the integrator’s conservation. Negative values are a symptom (step too large, noisy estimator, tolerance) and must be reported, not hidden (global rule: no fallbacks that hide bugs). New code: count and report negativity events, and fail in tests.

F8 (B): kernel boundary and fragmentation choices

• Pairs with m_i + m_j beyond the grid are skipped entirely for coagulation (kernel.py, near_upper_bound, k_h >= N_m), including their loss term. Particles in the top bins become inert, which creates the “traffic jam” the thesis describes. Make the policy explicit and configurable (ADR 0005).
• MRN fragments extend up to m_i + m_j. Dullemond’s notes argue for independent shattering of each body up to max(i, j), and add cratering for extreme mass ratios. Both are in scope for M4.
• Bouncing is absent (P_bounce = 0). It is part of the target model (Dominik and Dullemond 2024) and changes the equilibrium substantially.

F9 (C): diagnostic definitions

• Δ_stab is defined as a global relative mass error in ch. 3 but as a sum of per-bin relative errors in ch. 4. Use the global one.
• Δ_acc divides by ρ^c_i(t), which is zero in many bins (the thesis footnote mentions gaps). Replace by a mass-weighted norm that is well defined (see 03-verification.md).
• assert N.all() >= 0 (sampled.py:49) compares a boolean with 0 and is always true.

F10 (C): text and citation errors to fix before reuse

4. What the thesis left open (and Dullemond’s notes flag as TODO)

5. A gap nobody has addressed yet: the linear solve

For p particle properties the flattened state has N = N_m N_φ ⋯ entries, the full kernel sum costs O(N³) with fragmentation, and so does a dense LU of the N × N Jacobian in any implicit method. Sampling reduces the RHS cost to O(S · L) (S draws, L recipe entries per draw) but not the linear algebra. An implicit sampled scheme only wins if the Jacobian is also treated cheaply:

• Sparse Jacobian assembled from the same sampled recipes (O(S · L) non-zeros), solved by sparse LU or by Krylov (GMRES/BiCGStab) with a cheap preconditioner.
• Approximate Jacobians, which W-methods (Rosenbrock-W) and simplified-Newton Radau tolerate by design.
• Jacobian-free Newton-Krylov using sampled RHS evaluations.

This turns “does sampling pay off” from an RHS question into a whole-scheme question, and it is a natural, novel section of the paper (E13).

6. Summary for the plan

1. Re-implement the thesis model faithfully (with F6 to F8 fixed behind toggles) and reproduce its figures as a regression baseline.
2. Replace the biased estimator by unbiased IS / HT estimators, use the variance-optimal proposal, and re-measure the cost-accuracy trade-off with work-normalized efficiency.
3. Add variance reduction, the implicit sampled form, and cheap Jacobian strategies.
4. Move to two particle properties, where the payoff is expected, and compare against full summation, tensor-train or other fast deterministic methods where feasible, and the group’s RPMC.