Electromagnetic waves are described by Maxwell’s equations together with the constitutive equation of the considered medium. The latter equation in general may introduce complicated operators. As an example, for electron cyclotron (EC) waves in a hot plasma, an integral operator is present. Moreover, the wavelength and computational domain may differ by orders of magnitude making a direct numerical solution unfeasible, with the available numerical techniques. On the other hand, given the scale separation between the free-space wavelength $$\lambda _0$$
λ
0
and the scale L of the medium inhomogeneity, an asymptotic solution for a wave beam can be constructed in the limit $$\kappa = 2\pi L / \lambda _0 \rightarrow \infty$$
κ
=
2
π
L
/
λ
0
→
∞
, which is referred to as the semiclassical limit. One example is the paraxial Wentzel-Kramer-Brillouin (pWKB) approximation. However, the semiclassical limit of the wave field may be inaccurate when random short-scale fluctuations of the medium are present. A phase-space description based on the statistically averaged Wigner function may solve this problem. The Wigner function in the semiclassical limit is determined by the wave kinetic equation (WKE), derived from Maxwell’s equations. We present a paraxial expansion of the Wigner function around the central ray and derive a set of ordinary differential equations (phase-space beam-tracing equations) for the Gaussian beam width along the central ray trajectory.