  The low symmetry of monoclinic $\beta$-Ga$_2$O$_3$ leads to elaborate
intrinsic defects, such as Ga vacancies split amongst multiple lattice sites.
These defects contribute to fast, anisotropic Ga diffusion, yet their
complexity makes it challenging to understand dominant diffusion mechanisms.
Here, we predict the 3D diffusivity tensors for Ga interstitials
(Ga${_i^{3+}}$) and vacancies (V${_{Ga}^{3-}}$) via first principles and direct
solution of the master diffusion equations. We first explore the maximum extent
of configurationally complex ''$N$-split'' Ga interstitials and vacancies. With
dominant low-energy defects identified, we enumerate all possible elementary
hops connecting defect configurations to each other, including interstitialcy
hops. Hopping barriers are obtained from nudged elastic band simulations.
Finally, the comprehensive sets of (i) defect configurations and their energies
and (ii) the hopping barriers that connect them are used to construct the
master diffusion equations for both Ga${_i^{3+}}$ and V${_{Ga}^{3-}}$. The
solution to these equations yields the Onsager transport coefficients, i.e. the
components of the 3D diffusivity tensors $D_{{Ga}_i}$ and $D_{V_{Ga}}$ for
Ga${_i^{3+}}$ and V${_{Ga}^{3-}}$, respectively. It further reveals the active
diffusion paths along all crystallographic directions. We find that both
Ga${_i^{3+}}$ and V${_{{Ga}}^{3-}}$ diffusion are fastest along the $c$-axis,
due to 3-split defects that bridge neighboring unit cells along the $c$-axis
and divert diffusing species around high-energy bottlenecks. Although isolated
Ga${_i^{3+}}$ diffuse faster than isolated V${_{Ga}^{3-}}$, self-diffusion of
Ga is predominantly mediated by V$_{Ga}^{3-}$ due to the higher V$_{Ga}^{3-}$
defect concentration under most thermodynamic environments.
