The mean direction of waves in a directional distribution
It is helpful to understand the meaning of Fourier coefficients in the directional wave spectrum analysis.

A pitch-and-roll buoy measures time series of water elevation $\eta$, and slopes in both Cartesian directions, $\eta_x$ and $\eta_y$. The Fourier components $a_0$, $a_1$, $b_1$, $a_2$, $b_2$, are related to cross-spectra of elevation and slope time series (introduced by Longuet-Higgins 1963):

$$
a_0 = \dfrac{C_{11}}{\pi}
$$

$$
a_1 = \dfrac{Q_{12}}{\pi}
$$

$$
b_1 = \dfrac{Q_{13}}{\pi}
$$

$$
a_2 = \dfrac{C_{22}-C_{33}}{k^2\pi}
$$

$$
b_2 = \dfrac{2C_{23}}{k^2\pi}
$$

where the cross-spectra are:

$$
C_{11}(f) = C[\eta(f)\eta(f)] = \int_0^{2\pi} E(f,\theta)d\theta
$$

$$
Q_{12}(f) = -Q[\eta(f)\eta_x(f)]k^{-1} = -\int_0^{2\pi} E(f,\theta)\cos{\theta}d\theta
$$

$$
Q_{13}(f) = -Q[\eta(f)\eta_y(f)]k^{-1} = -\int_0^{2\pi} E(f,\theta)\sin{\theta}d\theta
$$

$$
C_{22}(f) = C[\eta_x(f)\eta_x(f)]k^{-2} = \int_0^{2\pi} E(f,\theta)\cos^2{\theta} d\theta
$$

$$
C_{23}(f) = C[\eta_x(f)\eta_y(f)]k^{-2} = \int_0^{2\pi} E(f,\theta)\sin{\theta}\cos{\theta} d\theta
$$

$$
C_{33}(f) = C[\eta_y(f)\eta_y(f)]k^{-2} = \int_0^{2\pi} E(f,\theta)\sin^2{\theta} d\theta
$$

From here you can see that the mean direction, as defined in the paper that you cite, is:

$$
\theta_0 = \arctan\left({\dfrac{b_1}{a_1}}\right) = \arctan\left({\dfrac{Q_{13}}{Q_{12}}}\right) = \arctan\left({\dfrac{\int_0^{2\pi}E(f,\theta)\sin{\theta}d\theta}{\int_0^{2\pi}E(f,\theta)\cos{\theta}d\theta}}\right)
$$

Because $Q_{12}$ and $Q_{13}$ are proportional to the integrated energy in the zonal and meridional directions, respectively, the arctangent of their ratio yields mean wave direction.

The authors could use higher moments to calculate direction, but this would have different physical meaning. In this specific case, using $\theta = \arctan\left({b_2/a_2}\right)$ would yield peak (dominant) wave direction in the context of pitch-and-roll buoys. 

You see that the above analysis is inherently limited by the quantities that are measured, specifically $\eta$, $\eta_x$, and $\eta_y$. If higher order quantities, say $\eta_{xx}$, $\eta_{xy}$, $\eta_{yy}$ were available, the directional spectrum could be described at a higher accuracy even using just Fourier decomposition.

Further reading:


Oltman-Shay and Guza (1984)
Nondirectional and directional wave data analysis procedures at NDBC