I'm showing how to emulate PM using additive oscillators and distortion/waveshaping.
The math behind this matches how PM works.
Phase modulation adds the output of the modulating signal to the phase (input) of the carrier signal, like this:
S_output(time) = A_carrier * S_carrier(time * ω_carrier + ɸ_carrier + ModIdx * S_modulator(time * ω_modulator + ɸ_modulator))
The input to the carrier is a sum of: constant carrier phase, time * frequency of the carrier, which can be thought of as a saw wave of said frequency, and modulator output, which can be any signal we want to modulate with.
On the other hand, waveshaping distortion works like this:
S_output(time) = A_post * S_curve(A_pre * S_input(time))
The input signal is multiplied by the pre-gain, passed through the curve, and amplified by post-gain. If we match the waveshaping curve with the shape of the carrier signal, we can emulate Phase Modulation. If we want multi-stage PM (O1 modulates O2 modulates O3, etc.), we can just attach multiple waveshapers in series, combining inputs to match modulations.