# Quadratic scaling path integral molecular dynamics for fictitious identical
particles and its application to fermion systems

Yunuo Xiong<EMAIL_ADDRESS>Center for Fundamental Physics and School
of Mathematics and Physics, Hubei Polytechnic University, Huangshi 435003,
China  Shujuan Liu  Center for Fundamental Physics and School of Mathematics
and Physics, Hubei Polytechnic University, Huangshi 435003, China  Hongwei
Xiong<EMAIL_ADDRESS>Center for Fundamental Physics and School of
Mathematics and Physics, Hubei Polytechnic University, Huangshi 435003, China

###### Abstract

Recently, fictitious identical particles have provided a promising way to
overcome the fermion sign problem and have been used in path integral Monte
Carlo (PIMC) to accurately simulate warm dense matter with up to 1000
electrons (T. Dornheim et al., arXiv:2311.08098 (2023)). The inclusion of
fictitious identical particles in path integral molecular dynamics (PIMD) can
provide another way to simulate fermion systems. In a recent paper (J. Chem.
Phys. 159, 154107 (2023)), Feldman and Hirshberg improved the recursive
formula for PIMD of N identical bosons, significantly reducing the
computational complexity from $O(PN^{3})$ to $O(N^{2}+PN)$. In this paper, we
extend this latest recursive formula for bosons to PIMD of fictitious
identical particles to improve the efficiency of simulating fermion systems.
We also provide the virial estimator for calculating energy by using the
recursive technique. As an example, we use the quadratic scaling PIMD for
fictitious identical particles to study the simulation of hundreds of fermions
in a two-dimensional periodic potential, in the hope of providing a simulation
tool for two-dimensional Fermi-Hubbard model and other strongly correlated
fermion systems, such as the simulation of ultracold fermionic gases in
optical lattices.

##  I Introduction

For a quantum many-body system with $N$ particles, the numerical expression of
the wave function $\Psi(\textbf{r}_{1},\cdots,\textbf{r}_{N})$ increases
exponentially with the increase of the number of particles. For example, in
three-dimensional space, if each dimension is divided into $M$ parts for
approximation, then the numerical expression of the wave function of the
quantum many-body system requires $M^{3N}$ orders of magnitude of complex
numbers CeperleyBook . When we use numerical simulation to study the
thermodynamic properties of quantum many-body systems, the path integral
formalism can well solve the difficulty of being unable to numerically express
the wave function of quantum many-body systems. In the path integral formalism
Feynman ; Tuckerman , the treatment of the partition function involves high-
dimensional integration, which can be handled by Monte Carlo or molecular
dynamics, resulting in path integral Monte Carlo (PIMC) and path integral
molecular dynamics (PIMD)barker ; Morita ; CeperleyRMP  for the exact
numerical simulation of the thermodynamic properties of quantum many-body
systems.

PIMC and PIMD have successfully solved the problem of exact numerical
simulation of the thermodynamic properties of large-scale bosonic systems
CeperleyRMP ; Burov1 ; Burov1b . Unfortunately, for PIMC/PIMD, when simulating
fermions, the fermion sign problem ceperley ; Alex ; troyer ; Dornheim ; WDM
poses a serious obstacle, and the computational cost increases exponentially
with the number of fermions or the decrease in temperature. Recently, we used
the idea of fictitious identical particles XiongFSP ; Xiong-xi  to try to
overcome the fermion sign problem, and found that it is indeed promising to
accurately simulate the thermodynamic properties of many important fermionic
systems from zero temperature to high temperature, and the computational
complexity is no different from the case of bosons. Afterwards, the idea of
fictitious identical particles for overcoming the fermion sign problem was
systematically verified and confirmed by Dornheim et al. Dornheim1 , and after
extending the applicability of fictitious identical particles, they highly
accurately simulated the static structure factor and static linear density
response function of warm dense matter with one thousand fermions Dornheim2 .
Previously, the computational cost of $O(10^{7})$ CPU hours could only
simulate at most 34 fermions Bohme  under equal conditions Dornheim2 . Based
on the extraordinary success of PIMC with fictitious identical particles for
large-scale fermionic systems, Dornheim et al. Dornheim1 ; Dornheim2  strongly
suggest that further research be conducted on the methods and applicability of
fictitious identical particles.

This work is based on the following two motivations:

1\. The computational complexity of simulating fermionic systems in the
fictitious identical particle PIMC framework is $O(N^{2})$ Dornheim1 ;
Dornheim2 . In the initial paper XiongFSP ; Xiong-xi  on fictitious identical
particle PIMD, we used the recursion formula for bosonic PIMD developed by
Hirsheberg et al. HirshPIMD , to simulate fermionic systems with a
computational complexity of $O(PN^{3})$. From the efficiency perspective,
fictitious identical particle PIMD and the corresponding PIMC are still not
competitive. Fortunately, Feldman and Hirshberg HirshImprove  have made
significant improvements to the recursion formula in their latest study of
bosonic PIMD, resulting in a computational complexity of $O(N^{2}+PN)$. It is
clear that it is urgent and important to incorporate the latest recursion
formula for bosonic PIMD into fictitious identical particle thermodynamics.
This work will successfully achieve this and provide open source code to
facilitate future researches.

2\. Despite being a simplified model, the numerical simulation of the two-
dimensional Fermi-Hubbard model has been a long-standing unsolved problem
LeBlanc ; Qin . In the two-dimensional Fermi-Hubbard model, exact numerical
simulations of more than 100 fermions at any temperature involve significant
computational cost. In this paper, we will show how to use a small server to
complete the simulation of more than 100 fermions in a two-dimensional
periodic potential. Although we analyze a continuous-space model in fictitious
identical particle PIMD, this analysis can clearly also promote the
understanding of the simplified Fermi-Hubbard model. The current work also
lays the foundation for the numerical simulation of the physics related to the
three-dimensional Fermi-Hubbard model in the future, as well as the numerical
simulation of large-scale strongly correlated fermion systems. One direct
application of this simulation can be expected to be the experimental study of
ultracold fermionic gases in optical lattices Fermilattice .

The organization of this paper is as follows. In Sec.II, we introduce
fictitious identical particle PIMD and the latest recursion technique that
greatly improves computational efficiency. In particular, we give the virial
estimator for energy. In Sec.III, we briefly analyze the reasons for the
success of fictitious identical particles in overcoming the fermion sign
problem. In Sec.IV, we introduce the two-component fermion model in a two-
dimensional periodic potential. In Sec.V, we use the isothermal extrapolation
method to present the simulation results of the energy of fermions at medium
temperatures. In Sec.VI, we use the constant energy semi-extrapolation method
to present the simulation results of the energy of fermions at any
temperature, and compare the energy results with those of the isothermal
extrapolation. In Sec. VII, we discuss general strategies for reliable
simulation of large-scale fermion systems. Finally, in Sec.VIII, we give a
summary and discussion.

##  II The thermodynamics of fictitious identical particles with PIMD

For any quantum system, the partition function is

$Z(\beta)=Tr(e^{-\beta\hat{H}}).$ (1)

Here $\beta=1/k_{B}T$ with $k_{B}$ the Boltzmann constant and $T$ the
temperature. The parametrized partition function for single-component
fictitious identical particles with a parameter $\xi$ can be written as

$Z(\xi,\beta)\sim\sum_{p\in S_{N}}\xi^{N_{p}}\int
d\textbf{r}_{1}d\textbf{r}_{2}\cdots
d\textbf{r}_{N}\left<p\\{\textbf{r}\\}|e^{-\Delta\beta\hat{H}}\cdots
e^{-\Delta\beta\hat{H}}|\\{\textbf{r}\\}\right>.$ (2)

We may also call the above parametrized partition function as partition
function of fictitious identical particles. $S_{N}$ represents the set of $N!$
permutation operations. The factor $\xi^{N_{p}}$ is due to the exchange effect
of fictitious identical particles. $\xi=+1$ for boson partition function,
while $\xi=-1$ for fermion partition function. In addition, $\\{\textbf{r}\\}$
denotes $\\{\textbf{r}_{1},\cdots,\textbf{r}_{N}\\}$. $N_{p}$ is a number
defined to be the minimum number of times for which pairs of indices must be
interchanged to recover the original order $\\{\textbf{r}\\}$ from
$p\\{\textbf{r}\\}$. With the technique of path integral, the parametrized
partition function $Z(\xi,\beta)$ can be mapped as a classical system of
interacting ring polymers. It is this idea that eliminates the exponential
complexity difficulty involved in numerically expressing the quantum many-body
wave function.

A recursion formula HirshPIMD ; HirshbergFermi  is found by Hirshberg et al.,
to calculate the partition function for both bosons ($\xi=1$) and fermions
($\xi=-1$). The computational complexity is $O(PN^{3})$. Here $P$ denotes the
number of beads for a single particle when path integral is used to express
the partition function. Following this recursion formula, there are many
researches based on PIMD for identical particles HirshbergFermi ; Deuterium ;
Xiong4 ; Xiong2 ; Xiong5 ; Xiong6 ; Xiong7 . We emphasize that in a recent
work by Feldman and Hirshberg HirshImprove , a significant step is made for
the recursion formula of identical bosons so that the computational complexity
becomes $O(N^{2}+PN)$. After this improvement, the computational complexity of
PIMD is at the same level as that of ordinary PIMC. It is clear that
incorporating the latest techniques of Feldman and Hirshberg into fictitious
identical particles is very valuable.

In this work, we incorporate the latest recursive method HirshImprove  for
bosons into the recursive expression for the parameterized partition function.
The specific results are as follows.

Firstly, based on the path integral formalism and the recursive method in Ref.
HirshImprove , the discretized partition function of Eq. (2) is found to be:

$Z(\xi,\beta)=\left(\frac{mP}{2\pi\hbar^{2}\beta}\right)^{PdN/2}\int
e^{-\beta(V_{\xi}^{[1,N]}+\frac{1}{P}U)}d\mathbf{R}_{1}...d\mathbf{R}_{N},$
(3)

where $\mathbf{R}_{i}$ represents the collection of ring polymer coordinates
$(\mathbf{r}_{i}^{1},...,\mathbf{r}_{i}^{P})$ corresponding to the $i$th
particle. $P$ denotes the number of beads for a single particle. The system
under consideration has $d$ spatial dimensions. $V_{\xi}^{[1,N]}$ considers
the exchange effects of the fictitious identical particles with parameter
$\xi$ by describing all the possible ring polymer configurations. $U$ is the
interaction between different particles, which is given by

$U=\sum_{l=1}^{P}V(\mathbf{r}_{1}^{l},...,\mathbf{r}_{N}^{l}).$ (4)

Here $V$ denotes the interaction potential. The expression of
$V_{\xi}^{[1,N]}$ is the key to consider the exchange effects of fictitious
identical particles, which is given by the following recursive relation:

$e^{-\beta
V_{\xi}^{[1,N]}}=\frac{1}{N}\sum_{k=1}^{N}\xi^{k-1}e^{-\beta(E^{[N-k+1,N]}+V_{\xi}^{[1,N-k]})},$
(5)

where $E^{[N-k+1,N]}$ is defined as

$E^{[N-k+1,N]}=\frac{1}{2}m\omega_{P}^{2}\sum_{l=N-k+1}^{N}\sum_{j=1}^{P}(\mathbf{r}_{l}^{j+1}-\mathbf{r}_{l}^{j})^{2},$
(6)

where $\mathbf{r}_{l}^{P+1}=\mathbf{r}_{l+1}^{1}$ unless $l=N$ where
$\mathbf{r}_{N}^{P+1}=\mathbf{r}_{N-k+1}^{1}$. The harmonic strings connecting
the ring polymers of different particles have the frequency
$\omega_{P}=\frac{\sqrt{P}}{\beta\hbar}$.

We note that evaluating all the $E^{[N-k+1,N]}$ directly would take
$O(N^{3}P)$ computational time, since there are $O(N^{2})$ of them in total
and evaluating each one takes $O(NP)$ time. Fortunately, Feldman and Hirshberg
HirshImprove  purposed an iterative scheme to evaluate all the $E^{[N-k+1,N]}$
by eliminating all the redundant calculations and the resulting computational
complexity is only $O(N^{2}+NP)$.

In order to carry out molecular dynamics simulations, we must also calculate
the force resulting from the potential $V_{\xi}^{[1,N]}$, which is its
gradient with respect to each bead coordinate
$\nabla_{\mathbf{r}_{l}^{j}}V_{\xi}^{[1,N]}$. There are $O(NP)$ forces in
total and the evaluation of each force takes $O(N^{2})$ time from the
recursive definition of $V_{\xi}^{[1,N]}$, so the total computational
complexity is $O(N^{3}P)$ as well. Feldman and Hirshberg HirshImprove  also
derived a formula for $\nabla_{\mathbf{r}_{l}^{j}}V_{\xi}^{[1,N]}$ that is
mathematically equivalent to the direct definition but takes only
$O(N^{2}+NP)$ time to calculate, we present the result as the direct
application of this technique to fictitious identical particles as follows.

Firstly, we define another set of potentials through

$e^{-\beta
V_{\xi}^{[u,N]}}=\sum_{l=u}^{N}\xi^{l-u}\frac{1}{l}e^{-\beta(E^{[u,l]}+V_{\xi}^{[l+1,N]})},$
(7)

from the new $V_{\xi}^{[u,N]}$ we define the connection probabilities
$Pr[G[\sigma](l)=l^{\prime}]$ between particles $l$ and $l^{\prime}$ as
follows. For $l^{\prime}\leq l$:

$Pr[G[\sigma](l)=l^{\prime}]=\xi^{l-l^{\prime}}\frac{1}{l}\frac{1}{e^{-\beta
V_{\xi}^{[1,N]}}}e^{-\beta(V_{\xi}^{[1,l^{\prime}-1]}+E^{[l^{\prime},l]}+V_{\xi}^{[l+1,N]})}.$
(8)

For $l^{\prime}=l+1$:

$Pr[G[\sigma](l)=l+1]=1-\frac{1}{e^{-\beta
V_{\xi}^{[1,N]}}}e^{-\beta(V_{\xi}^{[1,l]}+V_{\xi}^{[l+1,N]})}.$ (9)

For $l^{\prime}>l+1$:

$Pr[G[\sigma](l)=l^{\prime}]=0.$ (10)

Now we can calculate the gradients
$\nabla_{\mathbf{r}_{l}^{j}}V_{\xi}^{[1,N]}$ efficiently. First, for ”interior
beads” (beads with indices $1<j<P$ for all $l$), the gradient is trivial:

$\nabla_{\mathbf{r}_{l}^{j}}V_{\xi}^{[1,N]}=m\omega_{P}^{2}(2\mathbf{r}_{l}^{j}-\mathbf{r}_{l}^{j+1}-\mathbf{r}_{l}^{j-1}),$
(11)

while for ”exterior beads”, it has been shown that we can use the connection
probabilities to calculate the gradients as

$\nabla_{\mathbf{r}_{l}^{1}}V_{\xi}^{[1,N]}=\sum_{l^{\prime}=1}^{N}Pr[G[\sigma](l^{\prime})=l]\cdot
m\omega_{P}^{2}(2\mathbf{r}_{l}^{1}-\mathbf{r}_{l}^{2}-\mathbf{r}_{l^{\prime}}^{P}),$
(12)
$\nabla_{\mathbf{r}_{l}^{P}}V_{\xi}^{[1,N]}=\sum_{l^{\prime}=1}^{N}Pr[G[\sigma](l)=l^{\prime}]\cdot
m\omega_{P}^{2}(2\mathbf{r}_{l}^{P}-\mathbf{r}_{l}^{P-1}-\mathbf{r}_{l^{\prime}}^{1}).$
(13)

The overall complexity for the gradient calculation part is also
$O(N^{2}+NP)$.

The energy estimator is not given in the paper by Feldman and Hirshberg
HirshImprove . Here, we will give the result of the energy estimator for the
fictitious identical particles, with the application of the virial theorem.
Even for identical bosons with parameter $\xi=1$, we emphasize that the
present virial estimator for energy is an important improvement, compared with
the previous works HirshPIMD ; HirshImprove ; XiongFSP ; Xiong-xi , where the
virial theorem is not used.

Based on the partition function, the average energy is

$E=-\frac{\partial\ln Z(\xi,\beta)}{\partial\beta}.$ (14)

It is well known that a direct energy estimator derived from the above
relation has variance in the order of $O(\sqrt{P})$, so a different estimator
known as virial estimator Herman  is devised which will make the variance as
small as possible to minimize the error in the simulation of the energy.

When the periodic boundary condition is also considered, following the
previous worksHerman ; Spada , the virial estimator of energy in the context
of fictitious particle PIMD is :

$<E>=\frac{dN}{2\beta}+<\epsilon_{\xi}^{(N)}>+<\frac{1}{P}\sum_{j=1}^{P}V(\mathbf{r}_{1}^{j},...,\mathbf{r}_{N}^{j})>+<\frac{1}{2P}\sum_{j=1}^{P}\sum_{l=1}^{N}(\mathbf{r}_{l}^{j}-\mathbf{r}_{l}^{1})\cdot\nabla_{\mathbf{r}_{l}^{j}}V(\mathbf{r}_{1}^{j},...,\mathbf{r}_{N}^{j})>,$

where

$\epsilon_{\xi}^{(N)}=-\frac{1}{2}m\omega_{P}^{2}[p(\mathbf{r}_{1}^{1},...,\mathbf{r}_{N}^{1})-(\mathbf{r}_{1}^{P},...,\mathbf{r}_{N}^{P})]\cdot[p(\mathbf{r}_{1}^{1},...,\mathbf{r}_{N}^{1})-(\mathbf{r}_{1}^{1},...,\mathbf{r}_{N}^{1})],$
(15)

averaged over all permutations $p$ for $N$ particles weighted by
$\xi^{N_{p}}e^{-\beta E_{p}}$. Here
$(\mathbf{r}_{1}^{1},...,\mathbf{r}_{N}^{1})$ and
$(\mathbf{r}_{1}^{P},...,\mathbf{r}_{N}^{P})$ represent two dN-dimensional
vectors, while ”$\cdot$” denotes the inner product between two dN-dimensional
vectors. We note that for distinguishable particles ($\xi=0$), there is only
one permutation with no particle interchanges and $\epsilon_{\xi=0}^{(N)}=0$.
Following the original PIMD formulation HirshImprove , a recursive relation
for $\epsilon_{\xi}^{(N)}$ can be derived:

$\epsilon_{\xi}^{(N)}=\frac{\sum_{k=1}^{N}\xi^{k-1}(\epsilon_{\xi}^{(N-k)}-\tilde{E}_{N}^{(k)})e^{-\beta(V_{\xi}^{[1,N-k]}+E^{[N-k+1,N]})}}{\sum_{k=1}^{N}\xi^{k-1}e^{-\beta(V_{\xi}^{[1,N-k]}+E^{[N-k+1,N]})}},$
$\epsilon_{\xi}^{(0)}=0,$
$\tilde{E}_{N}^{(k)}=\frac{1}{2}m\omega_{P}^{2}\sum_{l=N-k+1}^{N}(\mathbf{r}_{l}^{P+1}-\mathbf{r}_{l}^{P})\cdot(\mathbf{r}_{l}^{P+1}-\mathbf{r}_{l}^{1}),$
(16)

where $\mathbf{r}_{l}^{P+1}=\mathbf{r}_{l+1}^{1}$ unless $l=N$ where
$\mathbf{r}_{N}^{P+1}=\mathbf{r}_{N-k+1}^{1}$.

Similarly to the original bosonic PIMD algorithm, evaluating all
$\tilde{E}_{N}^{(k)}$ takes $O(N^{3})$ time (there are $O(N^{2})$ of them and
calculating each one takes $O(N)$ time). Fortunately, following the quadratic
scaling bosonic PIMD algorithm purposed by Feldman and Hirshberg HirshImprove
, the computational complexity can be reduced to $O(N^{2})$ using the
following iterative equations:

$\tilde{E}_{l}^{(1)}=0,$
$\tilde{E}_{l}^{(u)}=\tilde{E}_{l}^{(u-1)}-\tilde{E}_{o}(\mathbf{r}_{l-u+2}^{1},\mathbf{r}_{l}^{P},\mathbf{r}_{l}^{1})+\tilde{E}_{o}(\mathbf{r}_{l-u+2}^{1},\mathbf{r}_{l-u+1}^{P},\mathbf{r}_{l-u+1}^{1})+\tilde{E}_{o}(\mathbf{r}_{l-u+1}^{1},\mathbf{r}_{l}^{P},\mathbf{r}_{l}^{1}),$
$\tilde{E}_{o}(\mathbf{r}_{1},\mathbf{r}_{2},\mathbf{r}_{3})=\frac{1}{2}m\omega_{P}^{2}(\mathbf{r}_{1}-\mathbf{r}_{2})\cdot(\mathbf{r}_{1}-\mathbf{r}_{3}).$
(17)

It is worth noting that the virial estimator presented here is also suitable
for systems without periodic boundary conditions, such as identical particles
trapped in a harmonic potential.

We can also apply the new recursive techniques to the case of two-component
fictitious identical particles. Assume there are $N_{\uparrow}$ particles in
state $|\uparrow>$ and $N_{\downarrow}$ particles in state $|\downarrow>$, the
parametrized partition function for two-component fictitious identical
particles is

$Z(\xi,\beta)\sim\sum_{p_{1}\in S_{N_{\uparrow}}}\sum_{p_{2}\in
S_{N_{\downarrow}}}\xi^{P_{1}}\xi^{P_{2}}\int
d\textbf{x}_{1}d\textbf{x}_{2}\cdots
d\textbf{x}_{N}\left<p_{1}p_{2}\\{\textbf{x}\\}|e^{-\Delta\beta\hat{H}}\cdots
e^{-\Delta\beta\hat{H}}|\\{\textbf{x}\\}\right>.$ (18)

Here $N=N_{\uparrow}+N_{\downarrow}$. ${\textbf{x}}_{j}$ denotes both the
spatial state and internal state $s_{j}$ ($|\uparrow>$ or $|\downarrow>$) of
the $j$th particle. $p_{1}$ ($p_{2}$) is the permutation of particles with
internal state $|\uparrow>$ ($|\downarrow>$). $P_{1}$ and $P_{2}$ are the
permutation counts for permutations $p_{1}$ and $p_{2}$, respectively.

The recursive formula for single-component fictitious identical particles can
be directly extended to the case of two-component fictitious identical
particles. The technical details of the implementation can be found in the
open-source code. The reliability of the new algorithm and code has been
verified by comparing the calculation results with those of the previous code.
We also confirmed the consistency of the results, compared with Fig. 6 in Ref.
Dornheim1  about the energy of spin-polarized electrons in a 2D harmonic trap
simulated by PIMC for identical fictitious particles.

##  III General considerations on the use of fictitious identical particles
for the simulations of fermion system

###  III.1 Fermion sign problem and fictitious identical particles

For fermion systems, in PIMD/PIMC, for identical fermions of the same spin
state, the energy of a fermion system at temperature $T$ can be generally
simulated as follows. Due to the antisymmetry of the wave function with
respect to the exchange of different fermions, based on Monte Carlo or
molecular dynamics sampling:

$E(T)=\frac{A(T)}{S(T)}.$ (19)

Here, $A$ and $S$ are both positive numbers.

Due to the antisymmetry of the fermion wave function, here

$S\sim e^{-N/T}.$ (20)

The average sign $S$ decreases exponentially with the increase of particles or
the decrease of temperature $T$. Since the energy $E$ of a fermion system is
finite, it must be that

$A\sim e^{-N/T}.$ (21)

Unfortunately, for large-scale fermion systems or small-scale fermion systems
at extremely low temperatures, $S$ and $A$ are too small to accurately
calculate the energy of fermion systems using direct PIMD/PIMC. Numerical
simulation becomes difficult Dornheim  when the average sign $S\sim 10^{-3}$,
and it is unlikely to simulate meaningful results when $S<10^{-4}$.

For fictitious identical particles, we set

$E(\xi,T)=\frac{A(\xi,T)}{S(\xi,T)}.$ (22)

In this case, we still have the following exponential behavior Dornheim1

$A\sim e^{\xi N/T},~{}~{}~{}S\sim e^{\xi N/T}.$ (23)

The exponential behavior of $A$ and $S$ can be explained by the expression of
the partition function (2), because the number of terms in the partition
function increases exponentially with $N$ when $\xi\neq 0$. A rigorous
mathematical proof of the exponential behavior of $A$ and $S$ is a
mathematical problem worth studying in the future. When $\xi<0$, we still face
the difficulties caused by the sign problem. However, from a physical point of
view, $\xi$ represents a type of exchange interaction of identical particles.
When $\xi>0$, it represents attractive exchange interaction, and the stronger
the attractive exchange interaction becomes as $\xi$ increases. When $\xi<0$,
it represents repulsive exchange interaction, and the stronger the repulsive
exchange interaction becomes as $\xi$ decreases. Therefore, for $\xi$, there
is a monotonically decreasing relationship between $E(\xi,T)$ and $\xi$; while
for $T$, there is a monotonically increasing relationship between $E(\xi,T)$
and $T$.

The simple and monotonic relationship of energy derived from physics above
must be reflected in Eq. (22). However, this embodiment seems to be disturbed
by the exponential behavior of $A$ and $S$. For example, we may set $A\sim
e^{c_{1}\xi N/T}$ and $S\sim e^{c_{2}\xi N/T}$ more specifically. Here $c_{1}$
and $c_{2}$ are two possible different coefficients. However, once the
coefficients $c_{1}$ and $c_{2}$ even have a slight difference, it will lead
to completely unreasonable energy results when $N$ is large. If the energy is
to be reasonable, then we must require $c_{1}$ and $c_{2}$ to be exactly equal
when $N$ is large. The larger $N$ or the smaller $T$, the more precise $c_{1}$
and $c_{2}$ need to be equal. This consideration leads to the expression of
energy can only be set to:

$E(\xi,T)=\frac{f_{1}(\xi,T)g(\xi,T)}{f_{2}(\xi,T)g(\xi,T)}.$ (24)

Here, $g(\xi,T)$ represents the part that decreases exponentially with the
increase of $N$ or the decrease of $\xi$, and decreases exponentially with the
decrease of T. We do not need to pay attention to the specific expression of
$g(\xi,T)$. $f_{1}(\xi,T)$ and $f_{2}(\xi,T)$ represent the non-exponential
parts. As we have discussed from the physical point of view, the exponential
behavior of the numerator and denominator must be canceled out, so we have:

$E(\xi,T)=\frac{f_{1}(\xi,T)}{f_{2}(\xi,T)}.$ (25)

###  III.2 Isothermal extrapolation

Since $f_{1}(\xi,T)$ and $f_{2}(\xi,T)$ in the above expression represent non-
exponential behavior, they can be assumed to be relatively simple functions.
Considering the monotonic behavior of $E(\xi,T)$ with respect to $\xi$, for
the same temperature, we can consider setting $E(\xi,T)$ to:

$E(\xi,T)\approx a+b\xi+c\xi^{2}.$ (26)

Based on the energy simulation results for $\xi\geq 0$, we can easily obtain
the three coefficients $a,b,c$. Then, we have the energy of the Fermi system
as $E(\xi=-1,T)=a-b+c$. Fortunately, the validity of this simple quadratic
function can be judged in the exact numerical simulation for $0\leq\xi\leq 1$.

The above isothermal extrapolation method XiongFSP  has been shown to be
highly efficient and accurate for simulating medium-scale and large-scale
Fermi systems in the case of weak quantum degeneracy Dornheim1 ; Dornheim2 . A
further improvement is the constant energy semi-extrapolation Xiong-xi , which
can be used to efficiently and accurately simulate the energy of fermions from
zero temperature (high quantum degeneracy) to medium and high temperature
(moderate and weak quantum degeneracy).

It is worth noting that in the above analysis, the exponential behavior of the
sign factor caused by the fermion sign problem is turned from a difficulty
into an advantage of the isothermal extrapolation with respect to $\xi$. Of
course, this advantage also comes from the monotonic behavior of $E(\xi,T)$
with respect to $\xi$ and $T$ from the physical point of view. We also note
that the above analysis is universally applicable due to the ubiquitous
existence of the fermion sign problem, which means that the fictitious
identical particles have application value in many important Fermi systems.
One of the purposes of the current work is to discuss the application value in
the strong correlated Fermi systems related to the Fermi-Hubbard model.

The above isothermal extrapolation can be applied in principle to density
distribution, static structure factor, static linear density response
function, and imaginary-time density correlation function, etc Dornheim1 . In
the work of Dornheim et al. Dornheim2 , the isothermal extrapolation method
and PIMC have been used to give highly accurate simulations of the static
structure factor and static linear density response function of 1,000
electrons in the warm dense matter.

###  III.3 Constant energy semi-extrapolation

In principle, the constant energy semi-extrapolation method is more accurate
and has a wider range of applicability than the isothermal extrapolation
method. One of the reasons is that in the constant energy semi-extrapolation,
we can also use the following exact properties Xiong-xi :

$\left.\frac{\partial E(\xi,T)}{\partial T}\right|_{T=0}=0.$ (27)

In Fig. 1, we show how to use the constant energy semi-extrapolation to infer
the energy of fermions. In the region where $\xi\geq 0$, we can obtain the
constant energy curve $\xi_{E}(T)$ for a given energy $E$ and the simulated
energy results. Compared with the isothermal curve $E_{T}(\xi)$, the constant
energy curve $\xi_{E}(T)$ has one more condition (27), so we think it is a
semi-extrapolation method. In Fig. 1, we show different constant energy
curves, where $E_{g}<E_{1}<E_{2}<E_{3}<E_{4}<E_{5}$. For the ground state
energy of fermions, the constant energy curve corresponding to $E_{g}$ just
passes through $\xi=-1$ and is perpendicular to the $\xi$ axis. For $E_{3}$,
the curve should also be perpendicular to the $\xi$ axis. The intersection
point of this constant energy curve and the horizontal line $\xi=-1$
determines the temperature at which the Fermi system has energy $E_{3}$.

Figure 1: The figure shows the basic method of constant energy semi-
extrapolation.

Constant energy semi-extrapolation also provides the reasons and applicability
of choosing the simple Eq. (25) in isothermal extrapolation to obtain the
accurate energy of a Fermi system. Considering the dependence of the parameter
$\xi$ on the exchange interaction in fictitious identical particles, there is
a monotonically increasing relationship between $\xi$ and $T$ for all constant
energy curves. Moreover, these constant energy curves $\xi_{E}(T)$ will not
intersect. When the energy is relatively large, the constant energy curve is
longer and has a smaller curvature, which makes it more linear in the region
of $-1\leq\xi\leq 1$. This leads to the relationship between $E$ and $\xi$
being more linear when the temperature is given, which in turn makes the use
of Eq. (25) in isothermal extrapolation reasonable. However, at low
temperatures, especially $T=0$, isothermal extrapolation is likely to cause
large deviations, because the curvature of the constant energy curve near
$E_{g}$ may be large at this time. In this case, constant energy semi-
extrapolation is a necessary method. We note that when $E<E_{g}$, since the
constant energy curve must always intersect the $\xi$-axis vertically at
$T=0$, the curvature of these low-energy constant energy curves will increase
as the energy decreases, resulting in the ”black hole” region described in the
paper Xiong-xi . In the ”black hole” region, the relationship between energy
and $\xi$ may have a turning point in the region of $\xi\leq 0$ for a given
temperature, which may lead to the failure of the isothermal extrapolation
method. Once the constant energy curve corresponding to $E_{g}$ is outside the
”black hole” region, the constant semi-extrapolation method can in principle
obtain the accurate energy of a Fermi system at any temperature.

##  IV Two-component Fermi systems in a two-dimensional periodic potential

We consider $N_{\uparrow}$ and $N_{\downarrow}$ fermions in different spin
states, respectively. The total number of fermions is
$N_{\uparrow}+N_{\downarrow}$. We consider the fermions in the following two-
dimensional periodic potential:

$V(x,y)=h\left(\cos^{2}(\pi x)+\cos^{2}(\pi y)\right).$ (28)

We choose periodic boundary condition to carry out numerical simulations. The
size of the box that satisfies the periodic boundary condition is $L\times L$
($L$ is a natural number greater than zero), and the number of grids contained
in the box is $L^{2}$ . We show the shape of this periodic potential in Fig.
2.

Figure 2:  Shown is the two-dimensional periodic potential considered in this
paper. More specifically, we also employ periodic boundary condition to study
the lattice system with size $L\times L$. In the schematic diagram, the size
of the lattice system is $5\times 5$, and the number of grids it contains is
$5^{2}=25$.

For $N$ fermions, the Hamiltonian operator is:

$\hat{H}=-\frac{1}{2}\sum_{l=1}^{N}\nabla_{l}^{2}+h\sum_{l=1}^{N}\left(\cos^{2}(\pi
x_{l})+\cos^{2}(\pi
y_{l})\right)+\sum_{l=1}^{N_{\uparrow}}\sum_{m=1}^{N_{\downarrow}}\frac{g}{\pi
s^{2}}e^{-\frac{(\mathbf{r}^{\uparrow}_{l}-\mathbf{r}^{\downarrow}_{m})^{2}}{s^{2}}}.$
(29)

Here $\mathbf{r}^{\uparrow}_{l}$ and $\mathbf{r}^{\downarrow}_{m}$ represent
the coordinate of the fermion with spin $|\uparrow>$ and $|\downarrow>$,
respectively. We adopt the usual convention $\hbar=k_{B}=m=1$ in this paper.

In the numerical simulation of this paper, the parameters we choose are
$h=1,g=1,s=0.5$. The value of $s$ we choose will make the interaction between
particles mainly limited to the scope of a grid, in order to be as close to
the Fermi-Hubbard model as possible; for this reason, we did not consider the
interaction between fermions of different spin states in this paper. Of
course, there is no problem in carrying out simulations after adding this
interaction without extra difficulty.

##  V The energy results of isothermal extrapolation

To verify the rationality of the algorithm, we first calculate the energy of a
two-dimensional periodic box. The Hamiltonian operator is as follows:

$\hat{H}=-\frac{1}{2}\sum_{l=1}^{N_{\uparrow}}\nabla_{l}^{2}-\frac{1}{2}\sum_{l=1}^{N_{\downarrow}}\nabla_{l}^{2}.$
(30)

We consider a box of size $4\times 4$ with periodic boundary condition, and
$N_{\uparrow}=N_{\downarrow}=8$. At the temperature $T=5$, based on the
fictitious identical particle PIMD, we show the energy at different $\xi$ when
$\xi\geq 0$ in Fig. 3 with red dots. The yellow dots are the energy of
$-1\leq\xi\leq 1$ directly calculated by the expression of the canonical
ensemble without using PIMD. When directly calculating the energy in the
canonical ensemble, we use the following recursive formula for the partition
function of fictitious identical particles:

$Z_{\xi}^{(N)}(\beta)=\frac{1}{N}\sum_{k=1}^{N}\xi^{k-1}z(k\beta)Z_{\xi}^{(N-k)}(\beta),$
(31) $z(\beta)=\sum_{i}e^{-\beta\epsilon_{i}}.$ (32)

$\epsilon_{i}$ is the single-particle energy level. We note that the results
obtained by the two different methods are consistent when $\xi\geq 0$, which
proves the correctness of the PIMD algorithm and code.

At medium or high temperatures, we can accurately determine the energy of
fermions through the method of isothermal extrapolation. Here we demonstrate
this. In isothermal extrapolation, we use the quadratic function
$E(\xi)=a+b\xi+c\xi^{2}$. After we determine the coefficients $a,b,c$ by
fitting the region $\xi\geq 0$, the energy of the fermion is $a-b+c$. The
black line in Fig. 3 is the curve of isothermal extrapolation. In the absence
of interaction, we can calculate it without using PIMD under the canonical
ensemble, and the result is represented by yellow points. We note that the
yellow points agree with the entire black line. This result confirms our
algorithm and code.

Figure 3: The red points (PIMD) and yellow points (directly based on the
expression of the canonical ensemble) are the calculated results of $16$
fermions in a $4\times 4$ box with periodic boundary condition. The black line
represents the curve of isothermal extrapolation. The blue points are the
energy of $100$ fermions in a $10\times 10$ lattice system under a two-
dimensional periodic potential, based on the PIMD simulation. At the same
temperature of $T=5$, the energy of different $\xi$ in the region of $\xi\geq
0$ is simulated by PIMD, and then the energy of $\xi=-1$, which is the energy
of the fermions, is obtained by extrapolation. The blue line is the curve of
isothermal extrapolation.

We now turn to the numerical simulation of the Hamiltonian operator (29) with
inter-particle interaction in the two-dimensional periodic potential. In the
case of a $10\times 10$ lattice system, we simulated 100 fermions with
$N=2N_{\uparrow}=100$ at a temperature of $T=5$. In the range of $0\leq\xi\leq
1$, we numerically simulated the energy for different $\xi$. The number of
beads we used was $P=40$, and the MD steps we used in the simulation of
massive Nosé-Hoover chain Nose1 ; Nose2 ; Hoover ; Martyna ; Jang  were
$3\times 10^{6}$. For the selection of the number of beads and MD steps, we
need to check their convergence in actual research. In this case, the
calculation of the energy for each parameter $\xi$ on a single CPU takes about
4 hours. For a multi-core computer, it is easy to break through the numerical
simulation of more than 100 fermions by calculating for each $\xi$
simultaneously. In the past, even with the use of supercomputers, direct
PIMD/PIMC for the Fermi-Hubbard model with hundreds of fermions was often
hopeless. In Fig. 3, the blue points (PIMD simulation results) and the blue
line (isothermal extrapolation curve) give the idea and results of isothermal
extrapolation. As with the PIMC simulation of the energy of fictitious
identical particles by Dornheim et al. Dornheim1 , we do not show the error
bar in the figure because the fluctuation of the simulated energy can be
neglected. We note that there are significant differences in the total energy
between fermions, distinguishable particles, and bosons, indicating that
quantum statistical effects or quantum degeneracy effects are very important
in this example. After fitting the quadratic polynomial about $E(\xi)$, we
have $E(\xi)=642.46-73.25\xi-2.47\xi^{2}$. We note that the linear term plays
a dominant role, so it can be expected that isothermal extrapolation can get a
reasonable fermion energy.

For the half-filled case, the isothermal extrapolation method gives an average
energy of 7.13 per fermion in the case of a $10\times 10$ lattice system. We
also calculated the half-filled case of $8\times 8$ lattices (64 fermions) and
$12\times 12$ lattices (144 fermions). We summarize the results in Table 1. We
note that the average energy per fermion slowly decreases as the size of the
box increases, as expected.

lattices | particle number | average energy per particle
---|---|---
$8\times 8$ | 64 | 7.27(6)
$10\times 10$ | 100 | 7.13(6)
$12\times 12$ | 144 | 7.02(6)
Table 1: Shown is the average energy per fermion for different numbers of
lattice sites at $T=5$.

##  VI Constant energy semi-extrapolation and energy at any temperature

The previous section (Sec.V) analyzed the case of medium temperatures. For the
cases of low temperature and zero temperature, we need to be especially
careful. In the case of low temperature, when the energy is given to
isothermal extrapolation, the real curve in the region of $\xi\leq 0$ may have
a turning point, which may lead to systematic bias in the extrapolation. This
problem can be solved by means of constant energy semi-extrapolation. In
constant energy semi-extrapolation, we establish the constant energy curve
$\xi_{E}(T)$ in the region of $\xi\geq 0$ for the given energy $E$, and then
infer the temperature of fermions based on the setting of $\xi=-1$. The reason
why we think this is a semi-extrapolation method is that in determining the
$\xi_{E}(T)$ curve, we not only use the simulation results of $\xi\geq 0$, but
also use the exact relationship $\left.\frac{\partial E(\xi,T)}{\partial
T}\right|_{T=0}=0$ to determine the $\xi_{E}(T)$ curve. This constraint
condition can play a key role in establishing the curve $\xi_{E}(T)$.

To help readers understand the methods in this paper, we first analyze the
case of a medium-sized $6\times 6$ lattice with half filling as an example. At
this time, $N_{\uparrow}=N_{\downarrow}=18$. In the four cases of
$\xi=0,0.25,0.5,1$, we calculate the energy in the temperature range of $2\leq
T\leq 6$. In Fig.4, we show the energy results of PIMD based on fictitious
identical particles for different $\xi$ at different temperatures.

Figure 4: The figure shows the energy simulation results of different $\xi$ at
different temperatures in the case of $N_{\uparrow}=N_{\downarrow}=18$ in the
two-dimensional periodic potential under a $6\times 6$ lattice system. From
top to bottom, they are $\xi=0,0.25,0.5,1$. The solid lines are obtained by
fitting with the quadratic function $E=a+bT+cT^{2}$.

The constant energy semi-extrapolation method in this paper is improved over
previous work Xiong-xi , which has better numerical stability for the energy
simulation fluctuations involved in $\xi\geq 0$. For a given energy $E$, we
obtain the corresponding temperatures for different $\xi$ according to the
results shown in Fig.4, which means that we obtain four sets of numbers
$\\{\xi=0,T_{1}(E)\\}$, $\\{\xi=0.25,T_{2}(E)\\}$, $\\{\xi=0.5,T_{3}(E)\\}$,
$\\{\xi=1,T_{4}(E)\\}$. For these four sets of numbers, we use the following
function to fit:

$\xi+d(E)\xi^{2}=a(E)+b(E)T^{2}+c(E)T^{3}.$ (33)

After fitting to obtain the coefficients $a(E),b(E),c(E),d(E)$, we can
determine the constant energy curve $\xi_{E}(T)$, and then infer the
temperature $T(E)$ of the fermion system at the given energy $E$ after taking
$\xi=-1$ and solving the equation above. Of course, in order to improve the
fitting, we can also add different $\xi$ to carry out the calculation and
verify the convergence.

In Fig.5, we show the fermion energy from zero temperature to high temperature
(blue points). For the case of medium and high temperature ($T\geq 2$), we
also compare it with the isothermal extrapolation method (red points) and get
consistent results. However, at low temperature, for example $T=1$, the energy
given by isothermal extrapolation has larger energy fluctuations due to the
relatively large curvature of the isothermal curve. Compared with isothermal
extrapolation, another advantage of constant energy semi-extrapolation is that
it uses more simulation energy information of $\xi\geq 0$ when inferring the
temperature for a given fermion system. Therefore, we notice that the
fluctuations of the energy given by constant energy semi-extrapolation can be
basically ignored. Another advantage of constant energy semi-extrapolation is
that it is much more efficient than isothermal extrapolation if we want to
calculate the energy of fermion systems at all different temperatures. We also
notice that in constant energy semi-extrapolation, the efficiency of
simulating fermion systems at any temperature is even higher than that of
boson systems, because for boson systems, we need to calculate the
thermodynamic properties at any temperature separately.

Figure 5: The constant energy semi-extrapolation method for 36 fermions in a
$6\times 6$ lattice system from low temperature to finite temperature (blue
points). The red points in the figure are the results at temperatures
$T=1,2,3,4,5$ under isothermal extrapolation. We note that when $T\geq 2$,
even though the quantum degeneracy leads to a significant difference in energy
between fermions and bosons, the results of isothermal extrapolation and
constant energy semi-extrapolation are consistent. The constant energy semi-
extrapolation is more accurate with smaller fluctuations. In the low
temperature case of $T=1$, the fluctuations of isothermal extrapolation are
very obvious, and the results of constant energy semi-extrapolation are more
reasonable. The error bar of the red points shows the energy fluctuations at
isothermal extrapolation.

In conventional methods, the farther away from half-filled case, the more
troublesome the problem of fermion sign problem becomes LeBlanc . However, the
calculation of the energy of a fermion system based on the thermodynamics of
fictitious identical particles will not be disturbed by this factor. We will
now calculate the case of $N=2N_{\uparrow}=144$ in a $14\times 14$ lattice
system. In this case, it is equivalent to the case of doping in the Fermi-
Hubbard model. In Fig.6, we use red points to show the average energy of each
fermion. The blue points are the average energy of each fermion in a $6\times
6$ lattice system under half filling. We note that the average energy of each
fermion is higher at half filling, and the two tend to be consistent as the
temperature tends to high, which agrees with physical expectations.

Figure 6: The blue dots in the figure represent the average energy of each
fermion in a $6\times 6$ lattice system at half filling. The red dots
represent the average energy of each fermion in a $14\times 14$ lattice system
with $N=2N_{\uparrow}=144$. The black dots represent the average of each
fermion in a $8\times 8$ lattice system at full filling.

For the case of full filling, we can also simulate it. Here we simulate the
case of $8\times 8$ lattices, $N=2N_{\uparrow}=128$ fermions, shown by black
dots in Fig. 6. We see that the average energy of each fermion is
significantly larger than the situation of half filling, which is reasonable
from physical consideration.

As a final example, we use the constant-energy semi-extrapolation method to
study the energy of different fermion numbers in a $4\times 4$ lattice. We
start from the half-filled situation of $N_{\uparrow}=N_{\downarrow}=8$ and
gradually reduce the number of fermions. In Fig.7, we show the average energy
per fermion for different fermion numbers at different temperatures. We note
that at $T=0$ (red dots), there is a sawtooth structure due to the
antiferromagnetic distribution. When $N_{\uparrow}\neq N_{\downarrow}$, the
antiferromagnetic effect is somewhat destroyed, leading to a lower average
energy per fermion. This sawtooth structure gradually disappears with
increasing temperature, as expected physically.

Figure 7:  For a $4\times 4$ lattice system at different temperatures, we
studied the average energy per fermion for
$N_{\uparrow}+N_{\downarrow}=8+8,8+7,7+7,7+6,6+6,6+5,5+5,5+4$. We note that at
$T=0$, the ground state of the fermion system exhibits a clear sawtooth
structure due to the antiferromagnetic effect. This sawtooth structure
gradually disappears with increasing temperature.

It is worth mentioning that in this paper, the medium temperature we mentioned
refers to the case where the quantum degeneracy effect is significant (the
energy of fermion systems and boson systems has a significant difference) but
the isothermal extrapolation is applicable. Low temperature refers to the case
where the isothermal extrapolation has obvious deviation and fluctuation, and
thus requires constant energy semi-extrapolation.

##  VII General strategies for reliable simulation of larger fermion systems

From the physical and mathematical perspectives, if the thermodynamics of
fictitious identical particles can handle dozens of fermions, then under the
same physical conditions, we have the hope of handling thousands of fermions.
We will discuss this under the condition of isothermal extrapolation, and
these discussions are also suitable for constant energy semi-extrapolation.

Suppose that for a $L\times L$ lattice system with $N$ fermions, we can
successfully infer the energy of the fermions through a quadratic function at
temperature $T$, which is the following expression:

$E(\xi,N)=a+b\xi+c\xi^{2}.$ (34)

Suppose that the above formula can accurately describe the case of
$-1\leq\xi\leq 1$. Now we increase the fermion system by four times, that is,
the number of lattices is four times, and the number of fermions is four
times. Under the condition that other physical conditions have not changed, we
naturally expect that the total energy will increase by about four times, that
is, there is:

$\frac{E(\xi,4N)}{4}\approx a+b\xi+c\xi^{2}.$ (35)

Therefore, we can assume that if the isothermal extrapolation can be
established in the case of $N$ fermions in a $L\times L$ lattice system, then
the isothermal extrapolation can also be established in the case of $4N$
fermions in a $2L\times 2L$ lattice system. Of course, the simple Eq. (35)
does not accurately know the fermion energy of the larger system. We need to
numerically simulate the energy of $4N$ fictitious identical particles in the
case of $\xi\leq 0$. In this case, we will get

$\frac{E(\xi,4N)}{4}\approx\tilde{a}+\tilde{b}\xi+\tilde{c}\xi^{2}.$ (36)

Here, the coefficients $\tilde{a},\tilde{b},\tilde{c}$ may have small
differences from $a$, $b$, and $c$. Therefore, we still need to conduct
numerical simulations on larger quantum systems to determine the new
coefficients $\tilde{a},\tilde{b},\tilde{c}$.

For large-scale fermion systems, it is difficult to have other numerical
simulation methods for independent comparison. A reliable strategy is to start
from a small system. For a small system of $L\times L$, we can use direct
PIMD/PIMC or other methods to numerically simulate $E(\xi)$ in the entire
range of $-1\leq\xi\leq 1$ and compare the results with the fictitious
identical particle thermodynamics. If this comparison is successful, we can
confidently use the fictitious identical particle method in a larger system of
$2L\times 2L$ with the particle number increased to four times. By
continuously increasing the scale of the system and monitoring the changes in
the coefficients $a,b,c$, we can gradually simulate larger physical systems
and assess their reliability. This strategy is also suitable for constant
energy semi-extrapolation.

In their latest work, Feldman and Hirschberg HirshImprove  have successfully
simulated over a thousand bosons using a new PIMD technique. Therefore, it is
completely realistic to simulate thousands of fermions using fictitious
identical particles and PIMD. This achievable goal will be an important
collaborative approach for PIMC, which has already achieved the simulation of
a thousand fermions using fictitious identical particles Dornheim2 .

##  VIII Summary and discussion

In summary, we have applied the latest recursive formula techniques of PIMD
for bosons to PIMD for fictitious identical particles. We found that the
computational efficiency for overcoming the fermion sign problem has been
greatly improved compared to our previous work XiongFSP ; Xiong-xi . Our study
shows that PIMD has the same important value as PIMC in simulating complex
fermion systems. In addition, we applied the new algorithm of fictitious
identical particle PIMD to a fermion system in a two-dimensional periodic
potential, demonstrating that we have the potential to efficiently and
accurately simulate the thermodynamic properties of fermions for this
important physical system. In this paper, we considered continuous space
fermion systems, but these studies still have the potential to promote the
understanding of the Fermi-Hubbard model, which approximates actual systems.
Additionally, how to directly apply fictitious identical particles to
approximate lattice models such as the Fermi-Hubbard model is also a valuable
research topic for the future.

The purpose of this paper is not to study the specific physical phenomena of
the Fermi-Hubbard model, but to provide algorithms and codes, and to use the
two-dimensional periodic potential as an example to illustrate our methods. It
is a bit unfortunate that due to the great difficulties brought by the fermion
sign problem, there are no other methods that can efficiently simulate
hundreds of fermions in a periodic potential at any temperature. Therefore, we
cannot compare our results with other theoretical results at this stage.
However, based on various physical considerations, we have reason to believe
that the methods and calculation results here are reliable, which provides a
benchmark for the development of other methods to overcome the fermion sign
problem nodes ; Helium ; Militzer ; Mak ; Blunt ; Malone ; Schoof1 ; Schoof2 ;
Schoof3 ; Yilmaz ; PB1 ; PB2 ; Joonho ; Groth ; SWZhang ; QinMP ; DornheimMod
, such as diagrammatic Monte Carlo diMonte ; diMonte2 ; Burov2 . One future
scenario for testing the methods and applications is the periodic optical
lattice of ultracold fermionic atoms. Based on the fictitious identical
particle thermodynamics, PIMD/PIMC and the capabilities of supercomputers, it
is realistic to simulate thousands of fermions from zero temperature to high
temperature. Therefore, we have hope to compare directly with the experiments
of ultracold fermionic atoms in the future Fermilattice ; Bloch ; Randeria ;
lattice .

It is worth noting that in our previous work and this work, we mainly studied
the energy of fermion systems based on fictitious identical particles. As
pointed out by Dornheim’s recent work Dornheim1 ; Dornheim2 , the $\xi$
extrapolation method can be applied to density distribution, static structure
factor, static linear density response function, and imaginary-time density
correlation function in addition to energy. How to make the $\xi$
extrapolation method more widely applicable to physical scenarios is also an
important issue to be studied in the future.

As a mathematically constructed hypothetical particle, although the original
purpose of fictitious identical particles was to overcome the fermion sign
problem in order to simulate large-scale fermion systems, fictitious identical
particles also have potential in the following two areas:

1\. The physics involved in small negative $|\xi|$. Even without the $\xi$
extrapolation method, we can also use direct PIMD/PIMC to perform accurate
numerical simulations when $\xi<0$ and is close to 0. In this case, we can
still observe the repulsive exchange force caused by $\xi<0$. Although this
$\xi$ value cannot directly tell us the exact properties of the fermion system
($\xi=-1$), it may bring us important insights about fermion systems.

2\. The duality between bosons and fermions. The $\xi$ parameter serves to
continuously connect boson systems and fermion systems. For example, if we
observe a phase transition of bosons, such as a specific heat peak, at
$\xi=1$, then this specific heat peak may also be observed in a fermion system
at $\xi=-1$. This duality between bosons and fermions may have potential
application in the connection between superfluidity of bosons and
superconductivity of fermions, and other phenomena.

In Fig.8, we show five possible phase transition cases. In the figure, we show
the curve resulting from the relationship between $\xi$ and the transition
temperature $T_{c}$, that is, for each $\xi$, we assume that there is a
corresponding transition temperature. For the first case, if we find that the
phase transition curve $\xi(T_{c})$ terminates at $\xi\geq 0$ in the region of
$\xi\geq 0$, then we can determine that this phase transition curve will not
reach the fermions. In Fig.8, the light yellow area represents the region that
PIMD/PIMC can directly simulate. Curve II represents the curve that terminates
in the light yellow area. Curve III represents the termination at $\xi\geq-1$.
Curve IV represents the termination at exactly $\xi=-1$. Curve V represents
the phase transition occurring at a finite temperature in the fermion system,
which is particularly interesting because it may have implications for low-
temperature superconductor, high-temperature superconductivity and even room-
temperature superconductivity.

Figure 8: Five possible phase transition situations shown in the $\xi-T_{c}$
plane.

The thermodynamics of fictitious identical particles provides a new dimension
for analyzing important physical properties such as phase transitions.
Considering that $\xi$ is a continuous real number, this new dimension has the
potential to play a significant role. In this work, we mainly use PIMD (which
is similar to PIMC) to analyze the thermodynamic properties of fictitious
identical particles. Of course, we cannot ignore the significance of other
important methods such as Hartree-Fock theory and density functional theory
applied to $\xi$ physics. In any case, the thermodynamics of fictitious
identical particles is a new field that is just beginning. It is beneficial to
explore the hidden value in it by using various existing mature methods. The
data that support the findings of this study are available from the
corresponding author upon reasonable request. The code of this study is openly
available in GitHub code .

###### Acknowledgements.

This work is partly supported by the National Natural Science Foundation of
China under grant numbers 11175246, and 11334001. This work has received
funding from Hubei Polytechnic University. We acknowledge the discussion with
Prof. Xianqing Lin.

## References

  * (1) R. M. Martin, L. Reining, and D. M. Ceperley, Interacting Electrons: Theory and Computational Approaches (Cambridge University Press, Cambridge, UK, 2016).
  * (2) R. P. Feynman and R. H. Albert, Quantum mechanics and path integrals (McGraw-Hill, New York, 1965).
  * (3) M. E. Tuckerman, Statistical mechanics: theory and molecular simulation (Oxford University, New York, 2010).
  * (4) J. A. Barker, A quantum-statistical Monte Carlo method; path integrals with boundary conditions, J. Chem. Phys. 70, 2914 (1979).
  * (5) T. Morita, Solution of the Bloch Equation for Many-Particle Systems in Terms of the Path Integral, J. Phys. Soc. Japan. 35, 980 (1973).
  * (6) D. M. Ceperley, Path integrals in the theory of condensed helium, Rev. Mod. Phys. 67, 279 (1995).
  * (7) M. Boninsegni, N. V. Prokof’ev, and B. V. Svistunov, Worm Algorithm for Continuous-Space Path Integral Monte Carlo Simulations, Phys. Rev. Lett. 96, 070601 (2006).
  * (8) M. Boninsegni, N. V. Prokof’ev, and B. V. Svistunov, Worm algorithm and diagrammatic Monte Carlo: A new approach to continuous-space path integral Monte Carlo simulations, Phys. Rev. E 74, 036701 (2006).
  * (9) D. M. Ceperley, Path Integral Monte Carlo Methods for Fermions, in Monte Carlo and Molecular Dynamics of Condensed Matter Systems, edited by K. Binder and G. Ciccotti (Editrice Compositori, Bologna, Italy, 1996).
  * (10) M. Troyer and U. J. Wiese, Computational Complexity and Fundamental Limitations to Fermionic Quantum Monte Carlo Simulations, Phys. Rev. Lett. 94, 170201 (2005).
  * (11) T. Dornheim, The Fermion sign problem in path integral Monte Carlo simulations: quantum dots, ultracold atoms, and warm dense matter, Phys. Rev. E 100, 023307 (2019).
  * (12) A. Alexandru, G. Basar, P. F. Bedaque, and N. C. Warrington, Complex paths around the sign problem, Rev. Mod. Phys. 94, 015006 (2022).
  * (13) T. Dornheim, S. Groth, and M. Bonitz, The uniform electron gas at warm dense matter conditions, Phys. Rep 744, 1 (2018).
  * (14) Y. Xiong and H. Xiong, On the thermodynamic properties of fictitious identical particles and the application to fermion sign problem, J. Chem. Phys. 157, 094112 (2022).
  * (15) Y. Xiong and H. Xiong, On the thermodynamics of fermions at any temperature based on parametrized partition function, Phys. Rev. E 107, 055308 (2023).
  * (16) T. Dornheim, P. Tolias, S. Groth, Z. A. Moldabekov, J. Vorberger, and B. Hirshberg, Fermionic physics from ab initio path integral Monte Carlo simulations of fictitious identical particles, J. Chem. Phys. 159, 164113 (2023).
  * (17) T. Dornheim, S. Schwalbe, Z. A. Moldabekov, J. Vorberger, and P. Tolias, Ab initio path integral Monte Carlo simulations of the uniform electron gas on large length scales, arXiv:2311.08098 (2023).
  * (18) M. Böhme, Z. A. Moldabekov, J. Vorberger, and T. Dornheim, Static electronic density response of warm dense hydrogen: Ab initio path integral Monte Carlo simulations, Phys. Rev. Lett. 129, 066402 (2022).
  * (19) B. Hirshberg, V. Rizzi, and M. Parrinello, Path integral molecular dynamics for bosons, Proc. Natl. Acad. Sci. USA 116, 21445 (2019).
  * (20) Y. M. Y. Feldman and B. Hirshberg, Quadratic Scaling Bosonic Path Integral Molecular Dynamics, J. Chem. Phys. 159, 154107 (2023).
  * (21) J. P. F. LeBlanc, et al., Solutions of the Two-Dimensional Hubbard Model: Benchmarks and Results from a Wide Range of Numerical Algorithms, Phys. Rev. X 5, 041041 (2015).
  * (22) M. Qin, T. Schäfer, S. Andergassen, P. Corboz, and E. Gull, The Hubbard model: A computational perspective, Annu. Rev. Condens. Matter Phys. 13, 275 (2022).
  * (23) T. Esslinger, Fermi-Hubbard physics with atoms in an optical lattice, Annu. Rev. Condens. Matter Phys. 1, 129 (2010).
  * (24) B. Hirshberg, M. Invernizzi, and M. Parrinello, Path integral molecular dynamics for fermions: Alleviating the sign problem with the Bogoliubov inequality, J. Chem. Phys. 152, 171102 (2020).
  * (25) C. W. Myung, B. Hirshberg, and M. Parrinello, Prediction of a supersolid phase in high-pressure deuterium, Phys. Rev. Lett. 128, 045301 (2022).
  * (26) Y. Yu, S. Liu, H. Xiong, and Y. Xiong, Path integral molecular dynamics for thermodynamics and Green’s function of ultracold spinor bosons, J. Chem. Phys. 157, 064110 (2022).
  * (27) Y. Xiong and H. Xiong, Numerical calculation of Green’s function and momentum distribution for spin-polarized fermions by path integral molecular dynamics, J. Chem. Phys. 156, 204117 (2022).
  * (28) Y. Xiong and H. Xiong, Path integral and winding number in singular magnetic field, Eur. Phys. J. Plus 137, 550 (2022).
  * (29) Y. Xiong and H. Xiong, Path integral molecular dynamics for anyons, bosons, and fermions, Phys. Rev. E 106, 025309 (2022).
  * (30) Y. Xiong and H. Xiong, Path integral molecular dynamics simulations for Green’s function in a system of identical bosons, J. Chem. Phys. 156, 134112 (2022).
  * (31) M. F. Herman, E. J. Bruskin, and B. J. Berne, On path integral Monte Carlo simulations, J. Chem. Phys. 76, 5150 (1982).
  * (32) G. Spada, S. Giorgini, and S. Pilati, Path-integral Monte Carlo worm algorithm for Bose systems with periodic boundary conditions, Condens. Matter 7, 30 (2022).
  * (33) S. Nosé, A molecular dynamics method for simulations in the canonical ensemble, Mol. Phys. 52, 255 (1984).
  * (34) S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. 81, 511 (1984).
  * (35) W. G. Hoover, Canonical dynamics: Equilibrium phase-space distributions, Phys. Rev. A 31, 1695 (1985).
  * (36) G. J. Martyna, M. L. Klein, and M. Tuckerman, Nosé-Hoover chains: The canonical ensemble via continuous dynamics, J. Chem. Phys. 97, 2635 (1992).
  * (37) S. Jang and G. A. Voth, Simple reversible molecular dynamics algorithms for Nosé-Hoover chain dynamics, J. Chem. Phys. 107, 9514 (1997).
  * (38) D. M. Ceperley, Fermion nodes, J. Stat. Phys. 63, 1237 (1991).
  * (39) D. M. Ceperley, Path-integral calculations of normal liquid 3He, Phys. Rev. Lett. 69, 331 (1992).
  * (40) B. Militzer, E. L. Pollock, and D. M. Ceperley, Path integral Monte Carlo calculation of the momentum distribution of the homogeneous electron gas at finite temperature, High Energy Dens. Phys. 30, 13 (2019).
  * (41) C. H. Mak, R. Egger, and H. Weber-Gottschick, Multilevel blocking approach to the fermion sign problem in path-integral Monte Carlo simulations, Phys. Rev. Lett. 81, 4533 (1998).
  * (42) N. S. Blunt, T. W. Rogers, J. S. Spencer, and W. M. Foulkes, Density-matrix quantum Monte Carlo method, Phys. Rev. B 89, 245124 (2014).
  * (43) F. D. Malone, N. S. Blunt, James J. Shepherd, D. K. K. Lee, J. S. Spencer, and W. M. C. Foulkes, Interaction Picture Density Matrix Quantum Monte Carlo, J. Chem. Phys. 143, 044116 (2015).
  * (44) T. Schoof, M. Bonitz, A. V. Filinov, D. Hochstuhl and J. W. Dufty, Configuration Path Integral Monte Carlo, Contrib. Plasma Phys. 51, 687 (2011).
  * (45) T. Schoof, S. Groth, and M. Bonitz, Towards ab initio thermodynamics of the electron gas at strong degeneracy, Contrib. Plasma Phys. 55, 136 (2015).
  * (46) T. Schoof, S. Groth, J. Vorberger, and M. Bonitz, Ab Initio Thermodynamic Results for the Degenerate Electron Gas at Finite Temperature, Phys. Rev. Lett. 115, 130402 (2015).
  * (47) A. Yilmaz, K. Hunger, T. Dornheim, S. Groth, and M. Bonitz, Restricted configuration path integral Monte Carlo, J. Chem. Phys. 153, 124114 (2020).
  * (48) T. Dornheim, S. Groth, A. Filinov, and M. Bonitz, Permutation blocking path integral Monte Carlo: a highly efficient approach to the simulation of strongly degenerate non-ideal fermions, New J. Phys. 17, 073017 (2015).
  * (49) T. Dornheim, T. Schoof, S. Groth, A. Filinov, and M. Bonitz, Permutation blocking path integral Monte Carlo approach to the uniform electron gas at finite temperature, J. Chem. Phys. 143, 204101 (2015).
  * (50) S. Groth, T. Dornheim, T. Sjostrom, F. D. Malone, W. M. C. Foulkes, and M. Bonitz, Ab initio exchange-correlation free energy of the uniform electron gas at warm dense matter conditions, Phys. Rev. Lett. 119, 135001 (2017).
  * (51) J. Lee, M. A. Morales, and F. D. Malone, A phaseless auxiliary-field quantum Monte Carlo perspective on the uniform electron gas at finite temperatures: Issues, observations, and benchmark study, J. Chem. Phys. 154, 064109 (2021).
  * (52) J. Carlson, S. Gandolfi, K. E. Schmidt, S. Zhang, Auxiliary-field quantum Monte Carlo method for strongly paired fermions, Phys. Rev. A 84, 061602 (2011).
  * (53) M. Qin, H. Shi, and S. Zhang, Benchmark study of the two-dimensional Hubbard model with auxiliary-field quantum Monte Carlo method, Phys. Rev. B 94, 085103 (2016).
  * (54) T. Dornheim, M. Invernizzi, J. Vorberger, and B. Hirshber, Attenuating the fermion sign problem in path integral Monte Carlo simulations using the Bogoliubov inequality and thermodynamic integration, J. Chem. Phys. 153, 234104 (2020).
  * (55) N. V. Prokof’ev and B. V. Svistunov, Bold diagrammatic Monte Carlo technique: When the sign problem is welcome, Phys. Rev. Lett. 99, 250201 (2007).
  * (56) E. Burovski, N. V. Prokof’ev, B. V. Svistunov, and M. Troyer, Critical Temperature Curve in BEC-BCS Crossover, Phys. Rev. Lett. 101, 090402 (2008).
  * (57) P. C. Hou, B. Z. Wang, K. Haule, Y. Deng, and K. Chen, Exchange-correlation effect in the charge response of a warm dense electron gas, Phys. Rev. B 106, L081126 (2022).
  * (58) I. Bloch, J. Dalibard, and W. Zwerger, Many-body physics with ultracold gases, Rev. Mod. Phys. 80, 885 (2008).
  * (59) M. Randeria and E. Taylor, Crossover from Bardeen-Cooper-Schrieffer to Bose-Einstein Condensation and the Unitary Fermi Gas, Annu. Rev. Condens. Matter Phys. 5, 209 (2014).
  * (60) I. Bloch, Ultracold quantum gases in optical lattices, Nat. Phys. 1, 23 (2005).
  * (61) https://github.com/xiongyunuo/MD-Lab