Calculations at finite temperatures are fundamental in different scientific fields, from nuclear physics to condensed matter. Evolution in imaginary time is a prominent classical technique for preparing thermal states of quantum systems. We propose a new quantum algorithm that prepares thermal states based on the quantum imaginary time propagation method, using a diluted operator with ancilla qubits to overcome the non-unitarity nature of the imaginary time operator. The presented method is the first that allows us to obtain the correct thermal density matrix on a general quantum processors for a generic Hamiltonian. We prove its reliability in the actual quantum hardware computing thermal properties for two and three neutron systems.