It is worth to note that Multiwfn is able to analytically calculate 3-order derivative tensor of orbital wavefunction, see "subroutine orbderv" in function.f90. Because electron density can be easily evaluated based on wavefunctions of occupied orbitals, realizing 3-order derivative of electron density in Multiwfn is straightforward.

In fact, in the "subroutine calchessmat_lapl" in function.f90, the 3-order derivative tensor of electron density (the array "rhotens3") is just the intermediate quantity for evaluating gradient of Laplacian of electron density, you can directly extract out this part of code.