  We have extended the canonical tree tensor network (TTN) method, which was
initially introduced to simulate the zero-temperature properties of quantum
lattice models on the Bethe lattice, to finite temperature simulations. By
representing the thermal density matrix with a canonicalized tree tensor
product operator, we optimize the TTN and accurately evaluate the thermodynamic
quantities, including the internal energy, specific heat, and the spontaneous
magnetization, etc, at various temperatures. By varying the anisotropic
coupling constant $\Delta$, we obtain the phase diagram of the spin-1/2
Heisenberg XXZ model on the Bethe lattice, where three kinds of magnetic
ordered phases, namely the ferromagnetic, XY and antiferromagnetic ordered
phases, are found in low temperatures and separated from the high-$T$
paramagnetic phase by a continuous thermal phase transition at $T_c$. The XY
phase is separated from the other two phases by two first-order phase
transition lines at the symmetric coupling points $ \Delta=\pm 1$. We have also
carried out a linear spin wave calculation on the Bethe lattice, showing that
the low-energy magnetic excitations are always gapped, and find the obtained
magnon gaps in very good agreement with those estimated from the TTN
simulations. Despite the gapped excitation spectrum, Goldstone-like transverse
fluctuation modes, as a manifestation of spontaneous continuous symmetry
breaking, are observed in the ordered magnetic phases with $|\Delta|\le 1$. One
remarkable feature there is that the prominent transverse correlation length
reaches $\xi_c=1/\ln{(z-1)}$ for $T\leq T_c$, the maximal value allowed on a
$z$-coordinated Bethe lattice.
