The derivation of diabatic spin-orbit (SO) Hamiltonians is presented, which are expanded in terms of nuclear coordinates to arbitrary order including the treatment of multi-mode systems, having more than one mode of the same symmetry. The derivation is based on the microscopic Breit-Pauli SO operator and the consequent utilization of time reversal and spatial symmetry transformation properties of basis functions and coordinates. The method is demonstrated for a set of 2E and 2A1 states in C3v∗ (double group) symmetry, once for a 3D case of one a1 and one e mode and once for a 9D case of three a1 and three e coordinates. It is shown that the general structure of the diabatic SO Hamiltonian only depends on the basis states and is strictly imposed by time reversal symmetry. The resulting matrix can be expressed easily by a power series using six parametrized structure matrices as expansion coefficients multiplied by the associated monomials in terms of symmetrized coordinates. The explicit example presented here provides a full-dimensional diabatic SO model for methyl halide cations, which will be studied in the future.