  We present the first successful application of the matrix product state (MPS)
representing a thermal quantum pure state (TPQ) in equilibrium in two spatial
dimensions over almost the entire temperature range. We use the Kitaev
honeycomb model as a prominent example hosting a quantum spin liquid (QSL)
ground state to target the two specific-heat peaks previously solved nearly
exactly using the free Majorana fermionic description. Starting from the
high-temperature random state, our TPQ-MPS framework on a cylinder precisely
reproduces these peaks, showing that the quantum many-body description based on
spins can still capture the emergent itinerant Majorana fermions in a ${\mathbb
Z}_2$ gauge field. The truncation process efficiently discards the high-energy
states, eventually reaching the long-range entangled topological state
approaching the exact ground state for a given finite size cluster. An
advantage of TPQ-MPS over exact diagonalization or purification-based methods
is its lowered numerical cost coming from a reduced effective Hilbert space
even at finite temperature.
