## 2 Theoretical Part
### Path Integrals  
For the calculation of the diagonal and single-momentum density matrices by Equations (11) and (15), we use the method of path integrals [38]. Let us consider the density matrix \(\rho(q^{A},q^{B})=\langle q^{B}|\mathrm{e}^{-\beta\hat{H}}|q^{A}\rangle\). Decomposing the statistical operator into the product of \(M\) high-temperature operators and using \(M-1\) complete sets of q-states, one can represent the density matrix in the form of a multiple integral:  
\[\rho(q^{A},q^{B})=\int\cdots\int\mathrm{d}^{3N}q^{1}\ldots\mathrm{d}^{3N}q^{M- 1}\prod_{k=0}^{M-1}\langle q^{k+1}|\mathrm{e}^{-\epsilon\hat{H}}|q^{k}\rangle _{q^{\hat{H}}=q^{A}}\text{,} \tag{17}\]  
where \(\epsilon=\beta/M\). Using the symbol \(|\{\mathbf{q}_{a}\}\rangle\) for the non-antisymmetrized \(q\)-state \(|\mathbf{q}_{1}\rangle\)\(|\mathbf{q}_{2}\rangle\)\(\ldots\)\(|\mathbf{q}_{N}\rangle\) and the symbol \(Pq\) for the permutation, one can take into account the Fermi statistics and rewrite the expression (17) as follows [39]:  
\[\rho(q^{A},q^{B})=\int\cdots\int\mathrm{d}^{3N}q^{1}\ldots\mathrm{ d}^{3N}q^{M-1}\sum_{P_{u},P_{d}}(-1)^{P_{u}+P_{d}}\\ \times\prod_{k=0}^{M-1}\langle\{\mathbf{q}_{a}^{k+1}\}|\mathrm{e} ^{-\epsilon\hat{H}}|\{\mathbf{q}_{a}^{k}\}\rangle\Big{|}_{q^{\hat{H}}=q^{A} \atop q^{\hat{H}}=Pq^{B}}\text{.} \tag{18}\]  
The non-antisymmetrized high-temperature matrix elements are well-known [34]:  
\[\langle\{\mathbf{q}_{a}^{k+1}\}|\mathrm{e}^{-\epsilon\hat{H}}|\{ \mathbf{q}_{a}^{k}\}\rangle\rangle=\lambda_{\epsilon}^{-3N}\\ \times\exp\Biggl{\{}-\frac{\pi}{\lambda_{\epsilon}^{2}}\sum_{a=1 }^{N}(\mathbf{q}_{a}^{k+1}-\mathbf{q}_{a}^{k})^{2}-\frac{\epsilon}{2}\left(U(q ^{k+1})+U(q^{k})\right)\Biggr{\}}+O(M^{-2})\text{,} \tag{19}\]  
where \(\lambda_{\epsilon}=\sqrt{2\pi\hbar^{2}\epsilon/m}\). Thus, the approximation for the density matrix with an accuracy of \(O(M^{-1})\) is:\[\rho(q^{A},q^{B})\approx\lambda^{-3N\Lambda}\int\cdots\int\mathrm{d}^{ 3N}q^{1}\cdots\mathrm{d}^{3N}q^{M-1}\sum_{p_{u},p_{d}}(-1)^{P_{u}+P_{d}}\\ \times\exp\left\{-\sum_{k=0}^{M-1}\sum_{a=1}^{N}\frac{m}{2}\left( \frac{\mathbf{q}_{a}^{k+1}-\mathbf{q}_{a}^{k}}{e\hbar}\right)^{2}-\sum_{k=0}^{ M-1}\frac{e}{2}\Big{(}U(q^{k+1})+U(q^{k})\Big{)}\right\}\Bigg{|}_{ \begin{subarray}{c}q_{ u}^{0}=q^{A}\\ q_{u}^{M}=P_{q}{}^{B}\end{subarray} }, \tag{20}\]  
where \(\lambda=\sqrt{2\pi\hbar\beta/m}\) is the thermal wavelength of electrons and \(C(M)\) is a constant depending on the number of high-temperature terms. At \(M\to\infty\) the multiple integral turns into a path integral over all \(3N\)-dimensional trajectories, and the expression becomes exact:  
\[\rho(q^{A},q^{B})=\lambda^{-3N}\sum_{P_{u},P_{d}}(-1)^{P_{u}+P_{d }}\int_{ \begin{subarray}{c}q(0)=q^{A}\\ q(\beta\hbar)=P_{q}{}^{B}\end{subarray} }\mathrm{D}^{3N}q(t)\\ \times\exp\left\{-\frac{1}{\hbar}\int_{0}^{\beta\hbar}\mathrm{d}t \left[\frac{m}{2}\sum_{a=1}^{N}\dot{\mathbf{q}}_{a}^{2}(t)+U(q(t))\right] \right\}. \tag{21}\]  
(Note that for some singular attractive potentials the continuous limit must be considered more carefully [40]). The Formula (21) contains the fermionic sign problem (FSP) due the sign-alternating permutations and, thus, cannot be used in Monte Carlo simulations directly. To avoid this problem we substitute the variables:  
\[\mathbf{q}_{a}(t)=\mathbf{z}_{a}(t)+\left(1-\frac{t}{\beta\hbar}\right)\mathbf{ q}_{a}^{A}+\frac{t}{\beta\hbar}\mathbf{q}_{Pa}^{B}\quad(a=1,\ldots,N), \tag{22}\]  
obtaining the path integral with zero boundary conditions:  
\[\rho(q^{A},q^{B})=\lambda^{-3N}\sum_{P_{u},P_{d}}(-1)^{P_{u}+P_{d }}\int_{ \begin{subarray}{c}z(0)=0\\ z(\beta\hbar)=0\end{subarray} }\mathrm{D}^{3N}z(t)\\ \times\exp\left\{-\frac{1}{\hbar}\int_{0}^{\beta\hbar}\mathrm{d}t \left[\frac{m}{2}\sum_{a=1}^{N}\dot{\mathbf{z}}_{a}^{2}(t)+U(q(t))\right] \right\}\exp\left\{-\frac{\pi}{\lambda^{2}}\sum_{a=1}^{N}\left[(\mathbf{q}_{Pa }^{B})^{2}-(\mathbf{q}_{a}^{A})^{2}\right]\right\}.