  Bilayer CrI$_3$ attracted much attention owing to peculiar switching between
the layered ferromagnetic and antiferromagnetic order upon stacking
alternation. This finding pointed out the importance of the apparently small
interlayer exchange, yet, existing literature addresses only its isotropic
part. To fill this gap, we combine the density functional theory with
Hamiltonian modeling to examine the anisotropic interlayer exchange in bilayer
CrI$_3$ - Dzyaloshinskii-Moriya (DMI) and the Kitaev interaction (KI). We
develop and apply a novel computational procedure that yields the off-diagonal
exchange matrix elements with $\mu$eV accuracy. Inspecting two types of bilayer
stacking, we found a weak interlayer KI and much stronger DMI between the
sublattices of monoclinic bilayer and their complete absence in rhombohedral
bilayer. We show how these anisotropic interactions depend on the interlayer
distance, stacking sequence, and the spin-orbit coupling strength and suggest
the dominant superexchange processes at play. In addition, we demonstrate that
the single-ion anisotropy largely depends on stacking, increasing by 50% from
monoclinic to rhombohedral structure. Remarkably, our findings prove that
iodines, owing to their spatially extended 5p orbitals featuring strong
spin-orbit coupling, are extremely efficient in mediating DMI across the van
der Waals gap in two-dimensional magnetic heterostructures. Given that similar
findings were previously demonstrated only in metallic multilayers where the
DMI shows a much longer range, our study gives promise that the chiral control
of spin textures can be achieved in two-dimensional semiconducting magnetic
bilayers whose ligands feature strong spin-orbit coupling.
