Appendix B: Computational Methods and Phase Portrait Blueprint

This appendix provides the computational tools for numerically integrating the model systems, computing equilibria, and producing phase portraits for the models developed in these notes. At the present stage this appendix is a computational blueprint: it records the numerical protocols and selected sweep values, but it does not include a complete set of figures. The code extends the attractor-figure script used to generate the Chapter 0 phase portrait (submission/rebuild_attractor_figure.py) to handle the general -state, -action template of Chapter 2.

No new theorems appear here. The purpose is to give the reader a self-contained computational pipeline for reproducing the qualitative attractor phenomena described in Chapters 0, 10, and 11, and to specify which phase portraits, basin slices, and parameter sweeps should be generated before making numerical basin claims.

The ODE has the form

A general-purpose integrator takes the model data and a generator family as input, and produces trajectories. Here the law equation carries the timescale parameter of Chapter 2, ; the shipped Chapter 0 figure of Section B.3 uses .

RK4 integration. The standard fourth-order Runge-Kutta method applies without modification to the -components. The distribution component requires care: after each integration step, must remain on the simplex . Because the generator equation preserves the simplex exactly in continuous time, the RK4 approximation stays close to the simplex for small time steps. For robustness, we project back onto after each step: set any negative component to zero and renormalize.

The existing rebuild_attractor_figure.py implements this for the Chapter 0 model with and . The extension to general requires:

1. Replacing the hard-coded rhs function with one that computes the softmax policy, occupancy measure, critic coefficients, actor drift, and generator dynamics from the data tuple and the current state.
2. Replacing the scalar damping with the diagonal matrix .
3. Adding the simplex projection after each RK4 step.

The computational cost per step is for the softmax and drift evaluation, plus for the generator dynamics. For the models in these notes (, , ), a single trajectory of length with step size requires roughly evaluations and runs in milliseconds.

Implementation note. The integrator should accept the model data as a structured input (a dictionary or data class with fields for , , , , , , , and the generator matrices ) and return time series for all state components. This makes it straightforward to switch between the Chapter 0, Chapter 10, and Chapter 11 models by changing only the data input.

Equilibrium finding. An equilibrium satisfies , which is a system of equations. For the models in these notes, this system has dimension at most 8 (the Chapter 10 model). Newton’s method applies: starting from an initial guess , iterate

where is the Jacobian of the vector field. The Jacobian can be computed analytically (from the derivatives of the softmax, the critic coefficients, and the generator) or numerically (by finite differences). For the models here, numerical differentiation is adequate.

Boundary equilibria. With a multi-dimensional actor the damping matrix vanishes only at a corner, where every coordinate is saturated. At a corner the actor equation is automatic, and the remaining system–the critic equation and the distribution equation –is linear and can be solved directly. On an open face, where only some coordinates satisfy , the damping kills only the saturated coordinates’ drift; each unsaturated coordinate still requires its raw drift to vanish, in the notation of Chapter 10, a nonlinear root-finding problem of the same kind as the interior case (Section 10.6 locates the Chapter 10 A-bubble this way, at with ). For (the Chapter 0 and Chapter 11 models) every boundary point of the actor box is a corner, so the linear shortcut applies.

Interior equilibria require solving the full nonlinear system. A practical approach: first solve the reduced system on the invariant-law graph (where and ), which is a system in alone (dimension ). For (the routing model), this is a scalar root-finding problem. For (the recommendation model), it is a two-dimensional root-finding problem. Then lift the solution to the full phase space by computing and .

Nullcline computation. A nullcline of the -th component is the surface . In the -plane at fixed , nullclines are curves that can be computed by sweeping: for each value of on a grid, solve for (or vice versa). The intersections of the -nullcline and the -nullcline are the equilibria.

For the Chapter 0 model, the nullclines in the -plane are already familiar from Section 0.14. For the Chapter 10 and 11 models, the nullclines live in higher-dimensional spaces and must be visualized as slices.

Chapter 0 asymmetric model (the existing figure). The attractor figure already shipped with the paper plots the asymmetric two-state model (, ): three trajectories, two of them starting outside the absorbing set , all enter and converge numerically to a single interior equilibrium at . A completed extended figure set should include:

• The -phase portrait with nullclines, multiple initial conditions, and the interior equilibrium marked, together with the two boundary equilibria at where the damping vanishes.
• Time series of , , and along one trajectory, showing the initial transient and the approach to the interior equilibrium.
• The -plane at , showing how the distribution coordinate co-evolves with the actor.

Chapter 0 symmetric model (planned figure item). The symmetric model (), whose interior equilibrium sits at , is generated by adapting the asymmetric figure script (Exercise 0.9). A useful plot is a set of time series launched from several initial conditions, used to explore numerically whether trajectories from different regions all approach the interior equilibrium or are instead captured by the boundary equilibria at –the boundary-stability question that Chapter 0, Section 0.14 leaves open. The plot records numerical observations, not a proved basin classification.

Chapter 0 asymmetric model (). Breaking the reward symmetry by varying from 0 to 2 produces a parameter sweep:

• For small (e.g., ), the asymmetry shifts the interior equilibrium and changes the relative sizes of the basins.
• For , the symmetric case is recovered.
• For , the opposite asymmetry appears.

The sweep illustrates how the attractor’s internal structure responds continuously to parameter changes.

Section 2.8 routing model. The three-state model with scalar actor and critic produces a -phase portrait qualitatively similar to Chapter 0 but with different equilibrium locations driven by the routing rewards. The computation plan includes:

• The -portrait at the frozen invariant law .
• A comparison of the three-state and five-state (Chapter 11) phase portraits, showing how the additional states change the equilibrium structure.

All basin and stability statements in this blueprint and in the parameter sweeps of Section B.5 are numerical observations in the sense of Section B.3, not proved classifications.

Chapter 10 recommendation model. The eight-dimensional phase space requires slicing. The most informative views are:

• -plane at the critic and distribution equilibrium: set and , and plot the reduced actor drift as a vector field in the -square . The four equilibria (diverse, A-bubble, B-bubble, dual bubble) should be marked as fixed points. Any basin boundaries shown in this slice are numerical evidence, not a theorem.

• Basin boundary map: Color each point by which equilibrium the trajectory from converges to. The resulting colored plot would show the basins of the four equilibria in that chosen slice.

• Parameter sweep: varying (the diversity reward). For ranging from 0.1 to 1.0, recompute the interior equilibrium, the saturated-face transverse signs, and the basin boundaries in the chosen slice.

• Parameter sweep: varying the deep-engagement exit rate. Increasing the exit rate from states 3 and 5 in both generators (making it easier for users to leave deep engagement) tests whether the boundary configurations become more or less transversely attractive.

Chapter 11 routing model. The six-dimensional phase space is easier to slice because :

• -phase portrait at : the reduced system is two-dimensional and can be plotted directly. The three equilibria (left-spoke, right-spoke, interior) should be marked.

• -plane: how the hub traffic fraction co-evolves with the routing policy. At the interior equilibrium, is substantial; at the spoke-heavy equilibria, is small.

• Time series comparison: overlay trajectories starting from hub-heavy, left-heavy, and right-heavy initial distributions to test which equilibria the trajectories approach in the chosen regime.

• Comparison with the three-state model: side-by-side phase portraits of the Section 2.8 and Chapter 11 models, highlighting how the additional states create competing spoke-heavy attractors.

The attractor’s sensitivity to model parameters is best seen through continuation: track how the equilibria and their basins change as one parameter varies continuously.

Method. For each parameter value on a grid, compute all equilibria (boundary and interior) using the methods of Section B.2. Plot the equilibrium coordinates as functions of the parameter. Mark stability changes (where an equilibrium gains or loses stability) and equilibrium creation or destruction (where two equilibria collide and annihilate).

This is not a formal bifurcation analysis — we do not compute normal forms or classify the bifurcation type. For those tools, a good entry point is Kuznetsov, Elements of Applied Bifurcation Theory, Chapters 2—3, or, for the one-parameter pitchfork picture, Strogatz, Nonlinear Dynamics and Chaos, Chapter 3. The purpose here is visual: to show the reader how the attractor landscape responds to design choices.

Chapter 10: deep-engagement reward sweep. Vary (the reward for curating in deep-niche states) from 0.5 to 2.0. As the deep-niche reward increases:

• The saturated-face transverse sign should be recomputed at each reward value. Only after that sign reverses should the boundary configuration be described as attracting from the interior.
• The interior equilibrium should be tracked as the reward changes; symmetry breaking may move it off the diagonal.
• Basin claims should be made only from the corresponding colored slice or a higher-dimensional sampling table.

Chapter 11: hub reward sweep. Vary (the hub reward for both actions) from 0.1 to 1.0. This sweep is instructive because it is a lever that does not work: the common hub reward cancels in the reward-advantage term, for every (Exercise 11.7(c)), so it enters only through the critic forcing . Numerically, the balanced equilibrium’s attraction weakens slightly as the hub reward rises: the slope of the bracket function of Section 11.5 moves from at to at (Section 11.6). The spoke-heavy boundary equilibria persist (they are created by damping, not by rewards). The plot shows the equilibria barely move under this sweep, in contrast with the secondary matching reward sweep, which moves the regime decisively.

Chapter 11: secondary matching reward sweep. Vary on a grid from 1.0 to 4.0 and recompute the bracket function at each value. Plot and the boundary values , against , marking the pitchfork threshold where changes sign (the balanced equilibrium loses stability, and a mirror pair of attracting interior equilibria splits off) and the lock-in threshold where the boundary values change sign (spoke lock-in becomes stable). This is the sweep cited in Sections 11.5—11.6, and it identifies the secondary reward asymmetry as the lever that actually moves the regime.

Chapter 11: topology sweep. Vary the direct hub-to-secondary rate by adding a shortcut edge from H to C with increasing rate. This changes the network topology continuously and shows how a shortcut that bypasses the primary node affects the attractor structure.

These computational experiments complement the analytical results of Chapters 10 and 11. The analysis identifies the equilibria and their character (boundary versus interior, damping-created versus gradient-zero). The computations described here are the route for quantifying basin structure and its response to design choices.