With the second quantum revolution on its way, there is a rising interest in the development of programmable quantum devices capable of outperforming classical computers in certain tasks. Since, in reality, a physical system can never be cooled down to the absolute zero temperature, the complete description of the physics behind the quantum computing devices has to account for the finite temperature effects. However, strategies for analyzing the properties of quantum systems at finite temperature, i.e. mixed quantum states, require handling of spaces where the number of parameters scales exponentially with the system sizes, and therefore, their application is restricted to small physical systems. The class of numerical algorithms based on tensor network methods has turned out to be successful in overcoming this difficulty, extending possible simulations to large number of particles. Within this thesis, we develop and implement an efficient Tree Tensor Network based algorithm for computing the finite temperature density matrix of the many-body quantum systems. We present the numerical techniques for computing the purity, Von Neumann entropy, Rényi entropies of any order, negativity, and entanglement of formation for the mixed state systems. We apply the method on systems composed of N = 8, 16, and 32 qubits. We apply our approach to the one-dimensional quantum Ising model and successfully distinguish the three quantum phases - ferromagnetic, paramagnetic, and the disordered finite-temperature phase, as well as verify the critical point as the point of the largest entanglement. Moreover, we use our method to study the systems of neutral Rydberg atoms trapped in the optical tweezer arrays, representing a physical quantum computing and simulation platform. On these two examples of physical systems, we show that the temperature range accessible to our numerical simulations depends on the gaps between the lower energy levels of the system’s spectrum.