Reproduce a streaming growth rate and a bounded patch

Freeze the experiment and install the small tools

Download the scientific scripts, inputs, frozen outputs and provenance below. The main calculation uses authored uniform gas/dust backgrounds and a small seeded perturbation. It is a numerical benchmark, not a fit to an observed disk. The separately supplied DSHARP brightness profile belongs to the last paper guide and never initializes a streaming mode.

Create a fresh environment and install the versions used for the frozen run. The input gives the two named published test cases, units, wave numbers, grid sizes, step and seed amplitude. The output records a SHA-256 hash of the unchanged input. Keep a copy before modifying parameters.

python3 -m venv .venv
source .venv/bin/activate
python3 -m pip install -r requirements.txt
python3 -W error replay-streaming.py --input streaming-input.json --out streaming-replay-rerun.json

Reading: Youdin & Goodman 2005, Streaming Instabilities, arXiv v1 · Benítez-Llambay, Krapp & Pessah 2019, FARGO3D multifluid tests, v3

Assemble two independent six-root problems

The first implementation uses the individual gas and dust equations from the previous chapter. It projects gas velocity onto two transverse basis vectors, then combines them with three dust velocities and one relative dust-density perturbation. NumPy diagonalizes this six-by-six complex matrix.

The second implementation directly transcribes the original Youdin–Goodman equations 26–31 in center-of-mass velocity, relative velocity, total-density perturbation and pressure. It retains the algebraic incompressibility constraint instead of eliminating it first. SciPy solves the resulting generalized eigenvalue problem; its finite eigenvalues are the six physical roots. The two infinite algebraic roots are constraint variables, not growing physical modes.

Compare all six finite roots. In the fixed examples the independently assembled systems agree to floating-point precision. This is stronger than computing a characteristic polynomial from a matrix and checking that same matrix again. It tests the variable change and the density feedback terms.

Reading: Youdin & Goodman 2005, Streaming Instabilities, arXiv v1

Compare a published full eight-state reference

For a weakly compressible gas, include its density and all three gas velocities, plus dust density and three dust velocities. Use $c_s/(\eta v_{\rm K})=20$ . The gas drag coefficient varies as $\rho_{\rm d}/\rho_{\rm g}$ , so its first-order feedback contains $\epsilon(\delta_{\rm d}-\delta_{\rm g})\mathbf W/t_{\rm stop}$ .

The full system gives LinA growth $0.4190091323\Omega$ and LinB growth $0.0154862262\Omega$ , matching Benítez-Llambay and colleagues 2019 Table 4 at its printed precision. Their Appendix E supplies the equations, and Section 3.5.3 equation 86 explains the small shift from the earlier Youdin–Johansen Table 1. The frozen output retains both published numbers and the independently calculated roots.

Frequency conventions need care. This tutorial uses $e^{\sigma t+ikx}$ , so positive real part means growth. The 2019 paper uses $e^{-\omega t+ikx}$ , making its tabulated complex eigenvalue the negative of ours. Its density amplitudes are absolute while the six-state calculation uses a relative dust amplitude; rescale an eigenvector before comparing its entries.

Reading: Benítez-Llambay, Krapp & Pessah 2019, FARGO3D multifluid tests, v3 · Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1

Let time integration measure the growth

Normalize the growing six-state eigenvector so the relative dust-density coefficient is one, then multiply it by $10^{-6}$ . The independent DOP853 ODE run integrates the actual complex amplitudes for $4/\Omega$ . It does not multiply a drawn pattern by a prescribed exponential. Fit the logarithm of the integrated dust amplitude against time to measure the growth rate.

The default incompressible rates are 0.41903111 and 0.015456678 in inverse-angular-time units. The ODE fit recovers them. Keep the seed small enough that velocities and density remain perturbations about their uniform backgrounds. A large-amplitude image produced from linear equations would still be a linear extrapolation, even if it looked like clumps.

Reading: Youdin & Goodman 2005, Streaming Instabilities, arXiv v1 · Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1

Run an actual radial–vertical grid

The grid helper evolves gas and dust velocity fields plus relative dust density at every cell of a periodic square containing one wavelength in each direction. It uses centered spatial differences, a discrete Fourier pressure projection to preserve incompressible gas motion, and fourth-order Runge–Kutta time stepping. All three velocity components are retained. It solves linearized field equations, so it does not include nonlinear advection or saturation.

A centered derivative represents a resolved wave with an effective wave number different from the continuum value. For each grid, seed the corresponding discrete eigenvector and compare its measured growth with that discrete expectation; no growth rate is supplied to the integrator. The same physical wavelength remains in the box. Increasing cells then moves the discrete rate toward the continuum result.

k_{\rm eff}=\frac{\sin(k\Delta x)}{\Delta x}

This is the exact derivative symbol of the centered stencil for a Fourier wave. It explains the spatial convergence before any physical interpretation.

Reading: Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1 · Benítez-Llambay, Krapp & Pessah 2019, FARGO3D multifluid tests, v3

Read the actual convergence and controls

The plot and table below read the frozen numerical output. For LinA, 16, 32 and 64 cells per wavelength give growth 0.41247010, 0.41739996 and 0.41862389, approaching the incompressible 0.41903111. For LinB they give 0.015472384, 0.015465589 and 0.015459218, approaching 0.015456678. The 32-cell half-step run separately checks temporal error.

A continuum eigenvector seeded on a coarse grid can excite several discrete modes. Early dust-density growth can then differ from the dominant discrete eigenvalue. Matching the seed to the discrete operator isolates the integrator test; refining the grid addresses the approximation to the continuum. Those are different questions.

The output also removes feedback, vertical variation and pressure-driven drift in separate controls. Their full incompressible systems do not retain positive growth. The equal-diffusivity example is explicitly a toy closure on every perturbation. The terminal-velocity cubic is stored as another comparison, with its stopping-time-times-wave-number measure: the short-wavelength LinA test lies outside its tight-coupling domain.

Actual synthetic LinA amplification. Solid cyan: independent temporal ODE. Dashed amber: measured RMS of the32-cell-per-wavelength radial–vertical grid, normalized by its initial RMS. Horizontal time uses inverse angular frequency; the two curves use their respective continuum/discrete eigenvector seeds.

Actual independent root and integrated-grid results, growth per inverse angular time. Six-state incompressible and full eight-state weakly compressible results are separate systems.
Case	Six-state root	Full eight-state root	32-cell grid	64-cell grid
linA	$0.4190311$	$0.4190091$	$0.4174000$	$0.4186239$
linB	$0.01545668$	$0.01548623$	$0.01546559$	$0.01545922$

Reading: Youdin & Goodman 2005, Streaming Instabilities, arXiv v1 · Umurhan, Estrada & Cuzzi 2020, Streaming Instability in Turbulent Disks, v4

Inspect the separate nonlinear radial calculation

The radial-results file stores actual finite-volume gas/dust snapshots and spatial/time-step comparisons for the pressure-trap experiment. Its mass sums and gas response can be checked without access to the interactive engine. The 128-cell default differs from the 1024-cell reference in both overall profile and peak height, as the drift chapter quantifies.

That radial model is nonlinear in density and velocity but has no vertical modes. The grid benchmark here is two-dimensional but linear. Neither is a nonlinear streaming-saturation calculation. To study clumping thresholds, use the stratified primary simulations in the next guide and retain their resolution, particle treatment and mixing assumptions.

Reading: Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1 · Li & Youdin 2021, Thresholds for Particle Clumping, v2

Exercises and explanations

Why can a coarse-grid seed constructed from a continuum eigenvector produce a transient growth rate different from the expected eigenvalue?
Read the explanation
The spatial stencil changes the linear operator. A continuum eigenvector is generally a mixture of discrete eigenvectors. Its components evolve at different rates and phases, producing a transient. A matched discrete seed isolates the time integrator; grid refinement tests the continuum approximation.
The full weakly compressible rate differs from the incompressible one in the fifth digit. Should the equations be altered to make them identical?
Read the explanation
They are different physical approximations. Compare each with its own independent reference, state the sound-speed and feedback conventions, and retain the actual difference.

Try the prediction

Connect the numerical benchmark to the playable model

Select the linear wave and export it before its small-amplitude stop. The stored complex coefficients and history describe the same density field.
Compare the integrated growth with the six-state incompressible benchmark, not the separate full eight-state value.
Return to the radial trap and explain which equation and geometry changes prevent its peak from being used as a streaming growth measurement.

Save these scientific files together and install the pinned NumPy/SciPy versions in requirements.txt. They run without a private checkout or interactive Lab engine.

Sources and calculation scope

Youdin & Goodman 2005, Streaming Instabilities, arXiv v1
10 September 2004 version 1; ApJ 620,459–469(2005). Full text read, Section 2 equations 1–7/12–17/24–31, Figures 1–3 and Section 4 mode geometry; Section 5 equation 39 and its restricted terminal-velocity domain; Section 6 distinguishes arbitrary linear amplitude from subsequent collapse. The PDF itself prints a later Draft Modified July 5, 2018 header; the archive version is explicitly frozen.
Youdin & Johansen 2007, Linear Evolution and Numerical Methods, v1
23 February 2007 version 1; ApJ 662,613–626. Sections 2.1–2.3 equations 1–14, Table 1 and Section 4/Figures 4–5. Fixed stopping time, two-way drag, three velocity components with radial–vertical dependence, finite-resolution eigenmode tests and weakly compressible published roots. Section 2.3 explicitly explains why steady radial traffic-jam arguments alone do not produce linear SI growth.
Benítez-Llambay, Krapp & Pessah 2019, FARGO3D multifluid tests, v3
Version 3, 10 June 2019; ApJS 241, 25. Read Sections 3.4–3.5, Table 4, equation 86 and Appendix E. The temporal exponential uses the negative of the paper's eigenvalue; the complex growth eigenvalue used here has the opposite sign. Table 4 gives the full eight-state LinA growth rate as 0.4190091323 and LinB as 0.0154862262, matching the independent matrix. Section 3.5.3 attributes the small difference from the 2007 table to gas-density drag-feedback terms; this distinction is already established in that paper.
Li & Youdin 2021, Thresholds for Particle Clumping, v2
14 October 2021 version 2; DOI 10.3847/1538-4357/ac0e9f. Full-text Sections 2,3.1–3.5 and 4: stratified simulation setup, surface metallicity versus volume loading, threshold map Figure 1, linear-growth maps Figure 7 and clumping comparison Figure 8. Equation 21 diffusion diagnostic fails to predict all nonlinear stratified clumping, explicitly cautioning against universal linear thresholds.
Umurhan, Estrada & Cuzzi 2020, Streaming Instability in Turbulent Disks, v4
25 May 2020 version 4, PDF header 26 May. Sections 3.1–3.3, equations 4–18 and Appendix A: explicit viscosity, dust diffusion/pressure and equilibrium forcing; Section 5 interprets growth with a turbulent dilution model. The Lab’s equal damping of all linear fields is a much simpler declared teaching closure, not a reproduction of this alpha-disk model.