A vertical mesh replaces a spatial derivative with a difference between neighbouring samples. That replacement slightly changes the magnetic coupling of a wave, so a real grid calculation need not grow at exactly the continuum rate. Refining the mesh should reduce the difference in a predictable way. The saved channel calculations let us check that prediction without generating a new trajectory.
Begin with the four real channel equations and the dispersion relation developed in the preceding chapters. You will need sine and cosine, derivatives as local slopes, and a sum as an average over samples. This chapter derives the grid's effective wave number, reads an actual spatial refinement, and explains why energy and stress should still be calculated from the saved fields.
Keep the fields as the evolving objects
The local coordinates are radial outward x, orbital y, and vertical z. The positive rotation rate is Ω, and the background orbital velocity is −qΩx. The density ρ and vertical magnetic field are uniform. Write magnetic disturbances in velocity units, b=δB/μ0ρ, with μ0 the magnetic permeability. Both b and the horizontal fluid velocity have units of length per time.
Distinguish a spatial field from its amplitude. Let Ux(z,t) and Uy(z,t) be the actual sampled horizontal velocity perturbations, and let Bx(z,t) and By(z,t) denote the horizontal magnetic disturbances in those velocity units. The upper-case magnetic symbols here do not denote SI magnetic-field strengths. The total orbital velocity in this rotating frame is
Here vA is the constant vertical Alfvén speed. Substituting Ux=ux(t)sin(kz) and Bx=bx(t)cos(kz), and similarly for the orbital components, recovers the four-amplitude equations. The magnetic force has a minus sign in the amplitude equation because the derivative of cosine is minus sine.
The independent offline calculation advances four arrays of sampled fields. Its spatial derivative reads neighbouring array values. A sine projection is made only when saving an observation; the projected amplitudes never replace the arrays or feed the next step. A single-mode experiment can have a four-dimensional description and still test a real-space derivative. These are two ways of representing the same restricted linear problem, with different numerical errors to examine.
A centred difference has its own wave number
Divide the periodic vertical length L into N equally spaced samples. Their spacing is h=L/N, with zj=jh and periodic indices: the neighbour after the last sample is the first one. A centred difference uses the slope across both neighbours,
(Dhf)j=2hfj+1−fj−1.
For a sinusoid, write the local phase as θj=kzj and the one-step phase change as η=kh. Use the addition formulas before taking any small-spacing approximation:
The quantity kh, called the modified or effective wave number, has the same inverse-length unit as k. The grid differentiates this sampled wave exactly as though its derivative coupling used kh. For a well-resolved fundamental wave, 0<kh<π, it is smaller than k. Expanding sine gives
kkh=khsin(kh)=1−6(kh)2+O((kh)4).
The leading derivative error therefore falls by about four when the spacing halves. For this mode, the spatially differenced equations, with time still continuous, use Kh=vAkh in place of the continuum coupling K=vAk. This statement describes the derivative. It does not assert that every possible diagnostic has exactly the same leading error.
The upper panel is an analytic reference: the fractional wave-number deficit for the fundamental periodic wave. The lower panel uses the actual saved four-component projection errors at the same final time. The dashed lower curve is a second-order reference anchored to the coarsest saved result. The two panels measure different quantities. Both axes showing error are dimensionless; the horizontal coordinate is the number of vertical samples.
Worked example: the same seed on three meshes
Use box length L0, time unit T0, and the saved parameters q=3/2, Ω=1/T0, K=0.6/T0. The fundamental wave has k=2π/L0, so the background speed is vA=0.6L0/(2πT0). It is not 0.6L0/T0.
At N=64, the angle between neighbouring samples is 2π/64. Thus
kkhKhT0≃0.998394393,≃0.599036636.
The continuum growing rate is sT0=0.44≃0.663324958. The same continuum initial vector is used on all three meshes. With radial amplitude A=10−4L0/T0, its four amplitudes are
Aa(0)=10.4/0.440.6/0.44−1.5.
To read the saved radial amplitude, multiply every actual velocity sample by the known sine pattern and average:
ax(t)=N2j=0∑N−1Ux,j(t)sj,sj=sin(N2πj).
The factor two compensates for the mean squared sine, which is one-half on these meshes. The other velocity uses the same sine, and the two magnetic components use cosine. Compare this four-component observation aN with the continuum vector ac at the same time:
eN=∥ac∥2∥aN−ac∥2.
The norm squares the four differences, adds them and takes a square root. This is the saved modal-vector error; it is not a separately measured maximum error at every grid point. At t=3T0, with the same time step 0.0025T0, the actual results are
Vertical samples
Projected vector error
Measured radial rate times the time unit
64
0.00128771789
0.662936179
128
0.000321942474
0.663227886
256
0.0000804864182
0.663300698
The error ratios are 3.999838 and 3.999960. They support second-order spatial convergence for this experiment. The measured rate is sN=t−1ln(∣ax(t)∣/A). Inserting Kh into the quartic gives shT0≃0.662855038 on the 64-sample mesh, close to but distinct from the measured finite-time rate. The original seed is a continuum eigenvector, not an exact eigenvector of the grid equations. Its small additional mode components help explain that distinction.
For a prediction at a new resolution, suppose the same second-order range continues. Halving the finest spacing suggests an error near 2.01×10−5. This fourth resolution is a prediction for a future comparison.
Sum the actual fields when measuring transport
For arbitrary sampled fields, the vertically averaged perturbation energy and the two outward momentum fluxes are
The two energy terms are the kinetic and magnetic pieces. Energy is energy per volume and each stress is force per area. When the fields are exactly the stated sine/cosine mode, averaging the squared patterns recovers the amplitude formulas with factors one-quarter and one-half. For a field containing other modes, products of four projected amplitudes need not equal these full-array sums. The saved diagnostics use the actual arrays.
Periodic centred differences also obey a useful cancellation. Shifting the summation index gives ∑jfj(Dhg)j=−∑jgj(Dhf)j. Magnetic tension therefore exchanges kinetic and magnetic energy without creating their sum in the continuous-time grid equations. The remaining source is qΩ(TRe,h+TMax,h). The next guide examines how a finite time step approximates that work.
Exercises with explanations
1. A wave alternates sign at every sample. Does the centred difference measure its steepest possible slope? No. For the alternating pattern, the two neighbours have the same value, so the difference is zero. This is the mesh's shortest represented wavelength. Increasing the wave number toward that limit does not keep increasing the effective coupling. A resolved fundamental mode and a near-grid-scale disturbance are different tests.
2. A calculation near the fastest continuum mode shows an unusually small growth-rate error. Does that prove fourth-order spatial accuracy? No. The growth curve is stationary at its maximum: its first derivative with respect to coupling vanishes there. A coupling error of second order can then produce a fourth-order change in the ideal growing eigenvalue. Finite-time projections can also contain other components. Check a nonstationary coupling and a field-vector error; the saved spatial series deliberately uses K=0.6/T0.
3. Why use the same continuum seed on all three spatial meshes? It exposes the effect of changing the spatial derivative while holding the physical initial disturbance fixed. Seeding each mesh with its own grid eigenvector answers a different useful question, used to isolate temporal error in the next guide.
4. Could a small projected error hide a grid-scale disturbance? Yes. A component orthogonal to the fundamental sine/cosine pattern can disappear from these four projections. Inspect actual fields and their energy/stress sums as well. Convergence here qualifies the observed fundamental; other disturbances require their own tests.
Sources
Balbus and Hawley, A Powerful Local Shear Instability in Weakly Magnetized Disks. I. Linear Analysis, ApJ 376, 214–222 (1991), printed pages 216–218, equations 2.2 and 2.5–2.9, supply the continuum component equations and dispersion relation. Primary scan. Stone and Gardiner, Implementation of the Shearing Box Approximation in Athena, arXiv:1006.0139v1, Sections II–III, define the local shear and discuss numerical source terms. Pinned primary source.
The spatial tables and field observations come from the eleven-case linear-channel-replay-public.zip, specifically records/spatial-64, records/spatial-128 and records/spatial-256. Athena uses separate finite-volume and constrained-transport algorithms. The centred-difference derivation here describes the bundle's simple periodic channel calculation. Reading these existing arrays adds zero new evolution steps.