  Tensor network states are expected to be good representations of a large
class of interesting quantum many-body wave functions. In higher dimensions,
their utility is however severely limited by the difficulty of contracting the
tensor network, an operation needed to calculate quantum expectation values.
Here we introduce a method for the time-evolution of three-dimensional
isometric tensor networks which respects the isometric structure and therefore
renders contraction simple through a special canonical form. Our method
involves a tetrahedral site-splitting which allows to move the orthogonality
center of an embedded tree tensor network in a simple cubic lattice to any
position. Using imaginary time-evolution to find an isometric tensor network
representation of the ground state of the 3D transverse field Ising model
across the entire phase diagram, we perform a systematic benchmark study of
this method in comparison with exact Lanczos and quantum Monte Carlo results.
We show that the obtained energy matches the exact groundstate result
accurately deep in the ferromagnetic and polarized phases, while the regime
close to the critical point requires larger bond dimensions. This behavior is
in close analogy with the two-dimensional case, which we also discuss for
comparison.
