Universal Energy Functionals for Trapped Fermi Gases in Low Dimensions

The zero-range interacting systems are good models for many physical systems, such as the ultracold Bose gases [2, 3, 4, 5, 6], ultracold Fermi gases [7, 8, 9, 10], and few-nucleon systems [11, 12, 13]. If the mean inter-particle distance $d$ and the thermal de Brogile wavelength $\lambda$ are both much larger than the range $r_{e}$ of the interaction between the particles, the system may be approximated as a zero-range interacting system, and it has universal properties that do not depend on the details of the interaction. These universal properties depend on the interaction potential through the $s$-wave scattering length $a$, which characterizes the low-energy scattering properties. This universality exists in the Bose systems [14, 15, 16, 17, 18, 19], the Fermi systems [20, 21, 22, 23], and the mixtures [24, 25, 26, 27].

For the 3D two-component Fermi systems with $s$-wave contact interaction, it was found that there exists a universal parameter $\mathcal{I}_{3\mathrm{D}}$, called contact, characterizing the tail of the momentum distribution at large $\mathbf{k}$, where $\hbar\mathbf{k}$ is the single-particle momentum, and $\hbar$ is Planck’s constant over $2\pi$, and that this tail is related to many other physical properties of the system through some exact relations [28, 29, 30, 31]. The name contact comes from the fact that it is a measure of the number of pairs of fermions in two different internal states with small separations. These exact relations [28, 29, 30, 31] have been generalized to the 1D two-component Fermi system [32], the 2D two-component Fermi system [33, 34, 35, 36], the spin-orbit-coupled Fermi system [37, 38], the Bose system [39, 40], and the mixtures [34, 40].

As a zero-range interacting system, the 3D two-component Fermi gas trapped in a smooth potential has an elegant property: its energy can be expressed as a linear functional of the occupation probabilities of single-particle energy eigenstates, i.e. [1]

$$E=\frac{\hbar^{2}\mathcal{I}_{3\mathrm{D}}}{4\pi ma}+\lim_{\epsilon_{\text{M}}\rightarrow\infty}\left(\sum_{\epsilon_{\nu}<\epsilon_{\text{M}}}\epsilon_{\nu}n_{\nu}-\frac{\hbar\mathcal{I}_{3\mathrm{D}}}{\pi^{2}}\sqrt{\frac{\epsilon_{\text{M}}}{2m}}\right),$$

(1)

where $m$ is the mass of each fermion, $\epsilon_{\nu}$ is the single-particle energy of the $\nu$th single-particle level in the specified smooth external potential, $n_{\nu}=n_{\nu\uparrow}+n_{\nu\downarrow}$, and $n_{\nu\uparrow}$ ($n_{\nu\downarrow}$) is the occupation probability of the spin up (down) state in the $\nu$th level. This general functional can be regarded as a generalization of the energy of trapped non-interacting Fermi gases,

$$E=\sum_{\nu\sigma}\epsilon_{\nu}n_{\nu\sigma}.$$

(2)

Since the zero-range interaction model is valid for lower spatial dimensions, a straightforward idea is to generalize the energy functional Eq. (1) to lower dimensions. The 1D and 2D two-component Fermi gases have been studied for many years. Experimentally, one can realize them by confining the particles in some transverse directions and allowing the particles to move freely in the remaining dimensions [41, 42, 43].

In this paper, we follow the method used in Ref. [1]. We first study the one-particle density matrices of the 2D and 1D trapped two-component Fermi gases with contact interactions. We then generalize the linear energy functional Eq. (1) to 2D and 1D.

This paper is organized as follows. In Sec. II, we introduce the normalized $N$-body energy eigenstate and the 2D Bethe-Peierls boundary condition. Using the boundary condition, we expand the one-particle density matrix at small displacements. In Sec. III, we combine the one-particle density matrix with the single-particle imaginary time propagator to find the universal energy functional in 2D:

$\displaystyle E=\lim_{\epsilon_{\mathrm{M}}\rightarrow\infty}\left(\sum_{\epsilon_{\nu}<\epsilon_{\mathrm{M}}}\epsilon_{\nu}n_{\nu}-\dfrac{\hbar^{2}\mathcal{I}_{2\mathrm{D}}}{4\pi m}\ln{\dfrac{e^{2\gamma}ma_{2\mathrm{D}}^{2}\epsilon_{\mathrm{M}}}{2\hbar^{2}}}\right),$

(3)

where $\gamma=0.5772\cdots$ is Euler’s constant, $e=2.718\cdots$ is the base of natural logarithm, $\mathcal{I}_{2\mathrm{D}}$ is the 2D contact, $a_{2\mathrm{D}}$ is the 2D scattering length between two fermions in different spin states, $n_{\nu}=\sum_{\sigma}n_{\nu\sigma}$, $\sigma=\uparrow,\downarrow$, and $n_{\nu\sigma}$ is the occupation probability of the spin-$\sigma$ state of the $\nu$th single-particle level. If the external potential is zero, the single-particle levels reduce to plane-wave states and Eq. (3) reduces to the energy theorem in Refs. [33, 34, 35, 36]. One can extract the contact $\mathcal{I}_{2\mathrm{D}}$ from the asymptotic behavior of $\rho(\epsilon)$ at large $\epsilon$, where

$$\rho(\epsilon)\equiv\sum_{\nu\sigma}n_{\nu\sigma}\delta(\epsilon-\epsilon_{\nu}),$$

(4)

and the coarse-grained version of $\rho(\epsilon)$ has the following asymptotic expansion at large $\epsilon$:

$$\rho(\epsilon)|_{\text{cg}}=\frac{\hbar^{2}\mathcal{I}_{2\mathrm{D}}}{4\pi m}\epsilon^{-2}+O(\epsilon^{-3}).$$

(5)

We also derived the occupation probabilities of high energy states; see Eq. (50). In Sec. IV and Sec. V, we do analogous calculations for the 1D two-component Fermi system and find that

	$\displaystyle E$	$\displaystyle=$	$\displaystyle-\dfrac{\hbar^{2}a_{1\mathrm{D}}\mathcal{I}_{1\mathrm{D}}}{2m}+\sum_{\sigma\nu}\epsilon_{\nu}n_{\nu\sigma},$		(6)
	$\displaystyle\rho(\epsilon)|_{\text{cg}}$	$\displaystyle=$	$\displaystyle\frac{\hbar^{3}\mathcal{I}_{1\mathrm{D}}}{2\sqrt{2}\pi m^{3/2}}\epsilon^{-5/2}+O(\epsilon^{-7/2}),$		(7)

where $a_{1\mathrm{D}}$ is the 1D scattering length between two fermions in different spin states, and $\mathcal{I}_{1\mathrm{D}}$ is the 1D contact. If the external potential is zero, the single-particle levels reduce to plane-wave states and Eq. (6) reduces to the energy theorem in Ref. [32]. We also derived the occupation probabilities of high energy states in 1D; see Eq. (72). In Sec. VI, we summarize our results and discuss the utilities and generalizations of our energy functionals.

We consider a trapped two-component Fermi system in 2D, with $N_{\uparrow}$ spin-up fermions and $N_{\downarrow}$ spin-down fermions. The total number is $N=N_{\downarrow}+N_{\uparrow}$. Here the trapping potential $V(\mathbf{r})$ is assumed to be smooth. First we calculate the one-particle density matrix. Consider a normalized $N$-body energy eigenstate

	$\displaystyle\ket{\Phi}$	$\displaystyle=$	$\displaystyle(N_{\uparrow}!N_{\downarrow}!)^{-1/2}\int D_{1}^{\uparrow}D_{1}^{\downarrow}\Phi(\mathbf{r}_{1}\dots\mathbf{r}_{N_{\uparrow}}\mathbf{s}_{1}\dots\mathbf{s}_{N_{\downarrow}})$		(8)
			$\displaystyle\times\psi_{\uparrow}^{\dagger}(\mathbf{r}_{1})\dots\psi_{\uparrow}^{\dagger}(\mathbf{r}_{N_{\uparrow}})\psi_{\downarrow}^{\dagger}(\mathbf{s}_{1})\dots\psi_{\downarrow}^{\dagger}(\mathbf{s}_{N_{\downarrow}})\ket{0},$

where $\mathbf{r}_{1},\dots,\mathbf{r}_{N_{\uparrow}}$ are the position vectors of the spin-up fermions, $\mathbf{s}_{1},\dots,\mathbf{s}_{N_{\downarrow}}$ are the position vectors of the spin-down fermions, $\psi_{\uparrow}^{\dagger}(\mathbf{r})$ is the creation operator for a spin-up fermion at position $\mathbf{r}$, $\psi_{\downarrow}^{\dagger}(\mathbf{s})$ is the creation operator for a spin-down fermion at position $\mathbf{s}$, $D_{i}^{\uparrow}\equiv\prod_{\mu=i}^{N_{\uparrow}}\mathrm{d}^{2}r_{\mu}$, $D_{i}^{\downarrow}\equiv\prod_{\mu=i}^{N_{\downarrow}}\mathrm{d}^{2}s_{\mu}$, and $\Phi(\mathbf{r}_{1}\dots\mathbf{r}_{N_{\uparrow}}\mathbf{s}_{1}\dots\mathbf{s}_{N_{\downarrow}})$ is the $N$-body wave function which is antisymmetric under the interchange of the positions of any two spin-up (spin-down) fermions. When $\mathbf{r}_{1}$ and $\mathbf{s}_{1}$ are close, the wave function satisfies the 2D Bethe-Peierls boundary condition

	$\displaystyle\Phi$	$\displaystyle=$	$\displaystyle A\left(\frac{\mathbf{r}_{1}+\mathbf{s}_{1}}{2};\mathbf{r}_{2}\dots\mathbf{r}_{N_{\uparrow}}\mathbf{s}_{2}\dots\mathbf{s}_{N_{\downarrow}}\right)$		(9)
			$\displaystyle\times\frac{1}{2\pi}\ln{\frac{a_{2\mathrm{D}}}{|\mathbf{r}_{1}-\mathbf{s}_{1}|}+O(|\mathbf{r}_{1}-\mathbf{s}_{1}|)},$

where $a_{2\mathrm{D}}$ is the two-dimensional $s$-wave scattering length, and $A$ is a function of $(N-1)$ position vectors. The one-particle density matrix for the spin-$\sigma$ fermions is defined as

$$p_{\sigma}(\mathbf{r},\mathbf{r}+\mathbf{b})=\bra{\Phi}\psi_{\sigma}^{\dagger}(\mathbf{r})\psi_{\sigma}(\mathbf{r}+\mathbf{b})\ket{\Phi}.$$

(10)

In particular, by substituting Eq. (8) into the above definition, we find that

	$\displaystyle p_{\uparrow}(\mathbf{r},\mathbf{r}+\mathbf{b})$	$\displaystyle=$	$\displaystyle N_{\uparrow}\int D_{2}^{\uparrow}D_{1}^{\downarrow}\Phi^{*}(\mathbf{r},\mathbf{r}_{2}\dots\mathbf{r}_{N_{\uparrow}}\mathbf{s}_{1}\dots\mathbf{s}_{N_{\downarrow}})$		(11)
			$\displaystyle\times\Phi(\mathbf{r+b},\mathbf{r}_{2}\dots\mathbf{r}_{N_{\uparrow}}\mathbf{s}_{1}\dots\mathbf{s}_{N_{\downarrow}}).$

We will expand $p_{\uparrow}(\mathbf{r},\mathbf{r}+\mathbf{b})$ through order ${O(b^{3})}$ at small distance $b$. Since $\Phi$ is singular when two fermions in different spin states are close, we divide the $2(N-1)$-dimensional integration domain into two regions: $\mathcal{C}_{\eta}$ and $\mathcal{D}_{\eta}$. $\mathcal{D}_{\eta}$ is the region in which every spin-down fermion lies outside of the circle of radius $\eta$ centered at $\mathbf{r}$, that is, $|\mathbf{s}_{\mu}-\mathbf{r}|>\eta$ for $\mu=1,\dots,N_{\downarrow}$, which is shown in Fig. 1(a). $\mathcal{C}_{\eta}$ is the complement of $\mathcal{D}_{\eta}$. We set $\eta$ small but $\eta>b$. In $\mathcal{C}_{\eta}$ the cases that two or more spin-down fermions come inside the small circle of radius $\eta$ centered at $\mathbf{r}$ are possible, but the contributions from these cases are suppressed by Fermi statistics and are of higher order than $b^{4}$, see Fig. 1(c). Next, we calculate the integrals in these two regions and add them up, then the dependencies on $\eta$ will be canceled.

In $\mathcal{D}_{\eta}$, we expand $\Phi$ in powers of $\mathbf{b}$ as

	$\displaystyle\Phi(\mathbf{r+b},\mathbf{R})$	$\displaystyle=$	$\displaystyle\Phi(\mathbf{r},\mathbf{R})+\nabla_{\mathbf{r}}\Phi(\mathbf{r},\mathbf{R})\cdot\mathbf{b}$		(12)
			$\displaystyle+\dfrac{1}{2}\sum_{i,j=1}^{2}\dfrac{\partial^{2}}{\partial r_{i}\partial r_{j}}\Phi(\mathbf{r},\mathbf{R})b_{i}b_{j}$
			$\displaystyle+T_{\mathbf{b}}+O(b^{4}),$

where $\mathbf{R}\equiv(\mathbf{r}_{2}\dots\mathbf{r}_{N_{\uparrow}}\mathbf{s}_{1}\dots\mathbf{s}_{N_{\downarrow}})$ and $T_{\mathbf{b}}\equiv\dfrac{1}{3!}\sum_{i,j,k=1}^{2}\frac{\partial^{3}}{\partial r_{i}\partial r_{j}\partial r_{k}}\Phi(\mathbf{r},\mathbf{R})b_{i}b_{j}b_{k}$. Let $I_{\mathcal{D}}$ be the integral evaluated in $\mathcal{D}_{\eta}$, and $I_{\mathcal{C}}$ be the integral evaluated in $\mathcal{C}_{\eta}$. We find

	$\displaystyle I_{\mathcal{D}}$	$\displaystyle=$	$\displaystyle n_{\uparrow}^{\mathcal{D}}(\mathbf{r})+\mathbf{u}_{\uparrow}(\mathbf{r})\cdot\mathbf{b}+\dfrac{1}{2}\sum_{i,j=1}^{2}v_{\uparrow,ij}(\mathbf{r})b_{i}b_{j}$		(13)
			$\displaystyle+T^{\prime}_{\mathbf{b}}+O(b^{4}),$

where

	$\displaystyle n_{\uparrow}^{\mathcal{D}}(\mathbf{r})$	$\displaystyle=N_{\uparrow}\lim_{\eta\to 0}\int_{\mathcal{D}_{\eta}}D_{2}^{\uparrow}D_{1}^{\downarrow}|\Phi(\mathbf{r},\mathbf{R})|^{2},$		(14)
	$\displaystyle\mathbf{u}_{\uparrow}(\mathbf{r})$	$\displaystyle=N_{\uparrow}\lim_{\eta\to 0}\int_{\mathcal{D}_{\eta}}D_{2}^{\uparrow}D_{1}^{\downarrow}\Phi^{*}(\mathbf{r},\mathbf{R})\nabla_{\mathbf{r}}\Phi(\mathbf{r},\mathbf{R}),$		(15)
	$\displaystyle v_{\uparrow,ij}(\mathbf{r})$	$\displaystyle=N_{\uparrow}\lim_{\eta\to 0}\int_{\mathcal{D}_{\eta}}D_{2}^{\uparrow}D_{1}^{\downarrow}\Phi^{*}(\mathbf{r},\mathbf{R})\dfrac{\partial^{2}}{\partial r_{i}\partial r_{j}}\Phi(\mathbf{r},\mathbf{R}),$		(16)
	$\displaystyle T^{\prime}_{\mathbf{b}}$	$\displaystyle=N_{\uparrow}\lim_{\eta\to 0}\int_{\mathcal{D}_{\eta}}D_{2}^{\uparrow}D_{1}^{\downarrow}\Phi^{*}(\mathbf{r},\mathbf{R})T_{\mathbf{b}}.$		(17)

To calculate the contributions from the region $\mathcal{C}_{\eta}$, we use the Bethe-Peierls boundary condition (9). The region $\mathcal{C}_{\eta}$ can be approximately partitioned into $N_{\downarrow}$ subregions, and in the $\mu$th subregion ($\mu=1,\dots,N_{\downarrow}$) $\mathbf{s}_{\mu}$ is within the circle of radius $\eta$ centered at $\mathbf{r}$. The contributions to $I_{\mathcal{C}}$ from these subregions are equal due to Fermi statistics. In the first subregion (shown in Fig. 1(b)) we have

	$\displaystyle\Phi(\mathbf{r},\mathbf{R}^{\prime})$	$\displaystyle=$	$\displaystyle A\left(\frac{\mathbf{r}+\mathbf{s}}{2};\mathbf{R}^{\prime}\right)\times\dfrac{1}{2\pi}\ln{\dfrac{a_{2\mathrm{D}}}{|\mathbf{r}-\mathbf{s}|}}$		(18)
			$\displaystyle+O(|\mathbf{r}-\mathbf{s}|),$
	$\displaystyle\Phi(\mathbf{r}+\mathbf{b},\mathbf{R}^{\prime})$	$\displaystyle=$	$\displaystyle A\left(\frac{\mathbf{r}+\mathbf{s}+\mathbf{b}}{2};\mathbf{R}^{\prime}\right)\times\dfrac{1}{2\pi}\ln{\dfrac{a_{2\mathrm{D}}}{|\mathbf{r}+\mathbf{b}-\mathbf{s}|}}$		(19)
			$\displaystyle+O(|\mathbf{r}+\mathbf{b}-\mathbf{s}|),$

where $\mathbf{s}\equiv\mathbf{s}_{1}$ and $\mathbf{R}^{\prime}\equiv(\mathbf{r}_{2}\dots\mathbf{r}_{N_{\uparrow}}\mathbf{s}_{2}\dots\mathbf{s}_{N_{\downarrow}})$. We then do the following expansions:

$$A\left(\frac{\mathbf{r}+\mathbf{s}}{2};\mathbf{R}^{\prime}\right)=A(\mathbf{r};\mathbf{R}^{\prime})-\dfrac{\mathbf{q}}{2}\cdot\nabla_{\mathbf{r}}A+O(q^{2}),$$

(20)

	$\displaystyle A\left(\frac{\mathbf{r}+\mathbf{s}+\mathbf{b}}{2};\mathbf{R}^{\prime}\right)=\,$	$\displaystyle A(\mathbf{r};\mathbf{R}^{\prime})+\left(\dfrac{\mathbf{b}}{2}-\dfrac{\mathbf{q}}{2}\right)\cdot\nabla_{\mathbf{r}}A$
		$\displaystyle+O(|\mathbf{b}-\mathbf{q}|^{2}),$		(21)

where $\mathbf{q}=\mathbf{r}-\mathbf{s}$. So we have

$$I_{\mathcal{C}}=N_{\uparrow}N_{\downarrow}\int D_{2}^{\uparrow}D_{2}^{\downarrow}\int_{q<\eta}\mathrm{d}^{2}qF_{\mathbf{b}}+O(b^{4}),$$

(22)

where

	$\displaystyle F_{\mathbf{b}}$	$\displaystyle=$	$\displaystyle\dfrac{1}{4\pi^{2}}\ln{\dfrac{a_{2\mathrm{D}}}{q}}\ln{\dfrac{a_{2\mathrm{D}}}{|\mathbf{q}+\mathbf{b}|}}\left(A^{*}-\nabla_{\mathbf{r}}A^{*}\cdot\dfrac{\mathbf{q}}{2}\right)$		(23)
			$\displaystyle\times\left[A+\nabla_{\mathbf{r}}A\cdot\left(\dfrac{\mathbf{b}}{2}-\dfrac{\mathbf{q}}{2}\right)\right].$

Carrying out the integral $I_{\mathcal{C}}$ and adding it to $I_{\mathcal{D}}$, we get

	$\displaystyle p_{\uparrow}(\mathbf{r},\mathbf{r}+\mathbf{b})=$	$\displaystyle\,I_{\mathcal{C}}+I_{\mathcal{D}}$
	$\displaystyle=$	$\displaystyle\,n_{\uparrow}(\mathbf{r})+\mathbf{u}_{\uparrow}(\mathbf{r})\cdot\mathbf{b}+\frac{1}{8\pi}C_{2\mathrm{D}}(\mathbf{r})b^{2}\ln{\frac{b}{a_{2\mathrm{D}}e}}$
		$\displaystyle+\dfrac{1}{2}\sum_{i,j=1}^{2}v_{\uparrow,ij}(\mathbf{r})b_{i}b_{j}$
		$\displaystyle+\dfrac{1}{16\pi}b^{2}\left(\dfrac{1}{2}\ln{\dfrac{b}{a_{2\mathrm{D}}}}-\dfrac{3}{8}\right)\mathbf{w}^{*}\cdot\mathbf{b}$
		$\displaystyle+\dfrac{1}{16\pi}b^{2}\left(\dfrac{3}{2}\ln{\dfrac{b}{a_{2\mathrm{D}}}}-\dfrac{11}{8}\right)\mathbf{w}\cdot\mathbf{b}$
		$\displaystyle+T^{\prime}_{\mathbf{b}}+O(b^{4}),$		(24)

where

	$\displaystyle n_{\uparrow}(\mathbf{r})$	$\displaystyle=$	$\displaystyle N_{\uparrow}\int D_{2}^{\uparrow}D_{1}^{\downarrow}|\Phi(\mathbf{r},\mathbf{R})|^{2},$		(25)
	$\displaystyle C_{2\mathrm{D}}(\mathbf{r})$	$\displaystyle\equiv$	$\displaystyle N_{\uparrow}N_{\downarrow}\int D_{2}^{\uparrow}D_{2}^{\downarrow}\ |A(\mathbf{r};\mathbf{R}^{\prime})|^{2},$		(26)
	$\displaystyle\mathbf{w}(\mathbf{r})$	$\displaystyle\equiv$	$\displaystyle N_{\uparrow}N_{\downarrow}\int D_{2}^{\uparrow}D_{2}^{\downarrow}\ A^{*}(\mathbf{r};\mathbf{R}^{\prime})\nabla_{\mathbf{r}}A(\mathbf{r};\mathbf{R}^{\prime}).$		(27)

$n_{\uparrow}(\mathbf{r})$ is the spatial density of spin-up fermions at $\mathbf{r}$, $C_{2\mathrm{D}}(\mathbf{r})$ is the 2D contact density, and $\mathbf{w}(\mathbf{r})$ is related to the center-of-mass motion of small-distance pairs of fermions in different spin states. We can also find a similar expansion for $p_{\downarrow}(\mathbf{r},\mathbf{r}+\mathbf{b})$.

We define an absolutely convergent series

$$J_{\sigma}(\beta)\equiv\sum_{\nu}n_{\nu\sigma}e^{-\beta\epsilon_{\nu}}=\sum_{\nu}\expectationvalue{c^{\dagger}_{\nu\sigma}c_{\nu\sigma}}{\Phi}e^{-\beta\epsilon_{\nu}},$$

(28)

where $\beta$ satisfies Re$\beta\geq 0$, $\ket{\Phi}$ is an $N$-body energy eigenstate,

$$n_{\nu\sigma}=\bra{\Phi}c_{\nu\sigma}^{\dagger}c_{\nu\sigma}\ket{\Phi}$$

(29)

is the occupation probability of the spin-$\sigma$ state of the $\nu$th single-particle level,

$$c_{\nu\sigma}=\int\mathrm{d}^{2}r\phi^{*}_{\nu}(\mathbf{r})\psi_{\sigma}(\mathbf{r})$$

(30)

is the fermion annihilation operator of such a single-particle state, and $\phi_{\nu}(\mathbf{r})$ is the wave function of the $\nu$th single-particle level in the trapping potential $V(\mathbf{r})$ and satisfies the single-particle Schrödinger equation

$$\left[-\frac{\hbar^{2}}{2m}\nabla^{2}+V(\mathbf{r})\right]\phi_{\nu}(\mathbf{r})=\epsilon_{\nu}\phi_{\nu}(\mathbf{r})$$

(31)

and the normalization condition

$$\int|\phi_{\nu}(\mathbf{r})|^{2}\mathrm{d}^{2}r=1.$$

(32)

We rewrite $J_{\sigma}(\beta)$ as

$$J_{\sigma}(\beta)=\int\mathrm{d}^{2}r\mathrm{d}^{2}r^{\prime}U_{\beta}(\mathbf{r},\mathbf{r}^{\prime})p_{\sigma}(\mathbf{r},\mathbf{r}^{\prime}),$$

(33)

where $U_{\beta}(\mathbf{r},\mathbf{r}^{\prime})=\sum_{\nu}e^{-\beta\epsilon_{\nu}}\phi_{\nu}(\mathbf{r})\phi^{*}_{\nu}(\mathbf{r}^{\prime})$ is the propagator of a single particle moving in the potential $V(\mathbf{r})$ whthin a time $-i\hbar\beta$. For a small positive $\beta$, at $|\mathbf{r}-\mathbf{r}^{\prime}|\gg\hbar\sqrt{\beta/m}$ the propagator is exponentially suppressed, while at $|\mathbf{r}-\mathbf{r}^{\prime}|\sim\hbar\sqrt{\beta/m}$ we have a short imaginary-time expansion

	$\displaystyle U_{\beta}(\mathbf{r},\mathbf{r}^{\prime})$	$\displaystyle=$	$\displaystyle\dfrac{m}{2\pi\hbar^{2}\beta}\left[1-\frac{V(\mathbf{r})+V(\mathbf{r}^{\prime})}{2}\beta\right]$		(34)
			$\displaystyle\times\exp\left[-\dfrac{m(\mathbf{r}-\mathbf{r}^{\prime})^{2}}{2\hbar^{2}\beta}\right]+O(\beta).$

Recall that when $|\mathbf{r}-\mathbf{r}^{\prime}|$ is small, we also have an expansion of $p_{\sigma}(\mathbf{r},\mathbf{r}^{\prime})$. Defining $G_{\sigma}=U_{\beta}(\mathbf{r},\mathbf{r}^{\prime})p_{\sigma}(\mathbf{r},\mathbf{r}^{\prime})$, we have

	$\displaystyle G_{\sigma}$	$\displaystyle=$	$\displaystyle\dfrac{m}{2\pi\hbar^{2}\beta}\left[1-\frac{V(\mathbf{r})+V(\mathbf{r}^{\prime})}{2}\beta\right]\exp\left[-\dfrac{mb^{2}}{2\hbar^{2}\beta}\right]$		(35)
			$\displaystyle\times\Bigg{[}n_{\uparrow}(\mathbf{r})+\mathbf{u}_{\uparrow}(\mathbf{r})\cdot\mathbf{b}+\sum_{i,j=1}^{2}v_{\uparrow,ij}(\mathbf{r})\dfrac{b_{i}b_{j}}{2}$
			$\displaystyle+\dfrac{b^{2}C_{2\mathrm{D}}(\mathbf{r})}{8\pi}\ln{\dfrac{b}{a_{2\mathrm{D}}e}}+\dfrac{b^{2}\mathbf{w}\cdot\mathbf{b}}{16\pi}\left(\dfrac{3}{2}\ln{\dfrac{b}{a_{2\mathrm{D}}}}-\dfrac{11}{8}\right)$
			$\displaystyle+\dfrac{b^{2}\mathbf{w}^{*}\cdot\mathbf{b}}{16\pi}\left(\dfrac{1}{2}\ln{\dfrac{b}{a_{2\mathrm{D}}}}-\dfrac{3}{8}\right)\Bigg{]}+O(\beta),$

where $\mathbf{b}=\mathbf{r}^{\prime}-\mathbf{r}$. Substituting the above result into Eq. (33), we find

	$\displaystyle J_{\sigma}(\beta)$	$\displaystyle=$	$\displaystyle N_{\sigma}-\beta\int\mathrm{d}^{2}rV(\mathbf{r})n_{\sigma}(\mathrm{r})$		(36)
			$\displaystyle-\dfrac{\hbar^{2}\mathcal{I}_{2\mathrm{D}}\beta}{8\pi m}\left(1+\gamma+\ln{\dfrac{ma_{2\mathrm{D}}^{2}}{2\hbar^{2}\beta}}\right)$
			$\displaystyle+\dfrac{\hbar^{2}\beta}{2m}\int\mathrm{d}^{2}r\sum_{i=1}^{2}v_{\sigma,ii}(\mathbf{r})+O(\beta^{2}),$

where

	$\displaystyle N_{\sigma}$	$\displaystyle=$	$\displaystyle\int\mathrm{d}^{2}rn_{\sigma}(\mathbf{r}),$		(37)
	$\displaystyle\mathcal{I}_{2\mathrm{D}}$	$\displaystyle=$	$\displaystyle\int\mathrm{d}^{2}rC_{2\mathrm{D}}(\mathbf{r}).$		(38)

Outside of the tiny range of two-body interactions, the $N$-body Schrödinger equation is simplified as

	$\displaystyle E\Phi$	$\displaystyle=$	$\displaystyle\sum_{\mu=1}^{N_{\uparrow}}\left[-\dfrac{\hbar^{2}}{2m}\nabla^{2}_{\mathbf{r}_{\mu}}+V(\mathbf{r}_{\mu})\right]\Phi$		(39)
			$\displaystyle+\sum_{\mu^{\prime}=1}^{N_{\downarrow}}\left[-\dfrac{\hbar^{2}}{2m}\nabla^{2}_{\mathbf{s}_{\mu^{\prime}}}+V(\mathbf{s}_{\mu^{\prime}})\right]\Phi,$

where $\mathbf{r}_{\mu}\neq\mathbf{s}_{\mu^{\prime}}$ for all $\mu,\mu^{\prime}$. Multiplying both sides of Eq. (39) by $\Phi^{*}$, integrating them over $\mathbf{r}_{1},\dots,\mathbf{r}_{N_{\uparrow}},\mathbf{s}_{1},\dots,\mathbf{s}_{N_{\downarrow}}$ for $\mathbf{r}_{\mu}\neq\mathbf{s}_{\mu^{\prime}}$ for all $\mu,\mu^{\prime}$, we get

$$\sum_{\sigma}\int\mathrm{d}^{2}r\left[V(\mathbf{r})n_{\sigma}(\mathbf{r})-\dfrac{\hbar^{2}}{2m}\sum_{i=1}^{2}v_{\sigma,ii}(\mathbf{r})\right]=E.$$

(40)

Summing Eq. (36) over $\sigma$, we find

	$\displaystyle\sum_{\sigma}J_{\sigma}(\beta)=$	$\displaystyle\,\sum_{\nu\sigma}n_{\nu\sigma}e^{-\beta\epsilon_{\nu}}$
	$\displaystyle=$	$\displaystyle\,N-\dfrac{\hbar^{2}\mathcal{I}_{2\mathrm{D}}\beta}{8\pi m}\left(1+\gamma+\ln{\dfrac{ma_{2\mathrm{D}}^{2}}{2\hbar^{2}\beta}}\right)$
		$\displaystyle-\beta\sum_{\sigma}\int\mathrm{d}^{2}r\left[V(\mathbf{r})n_{\sigma}-\dfrac{\hbar^{2}}{2m}\sum_{i=1}^{2}v_{\sigma,ii}\right]$
		$\displaystyle+O(\beta^{2})$
	$\displaystyle=$	$\displaystyle\,N-E\beta-\dfrac{\hbar^{2}\mathcal{I}_{2\mathrm{D}}\beta}{4\pi m}\left(1+\gamma+\ln{\dfrac{ma_{2\mathrm{D}}^{2}}{2\hbar^{2}\beta}}\right)$
		$\displaystyle+O(\beta^{2}).$		(41)

Let

$$\rho(\epsilon)\equiv\sum_{\sigma}\rho_{\sigma}(\epsilon)=\sum_{\nu\sigma}n_{\nu\sigma}\delta(\epsilon-\epsilon_{\nu}).$$

(42)

Equation (41) can be rewritten as

	$\displaystyle\int_{-\infty}^{+\infty}\rho(\epsilon)e^{-\beta\epsilon}\mathrm{d}\epsilon=$	$\displaystyle\,N-\dfrac{\hbar^{2}\mathcal{I}_{2\mathrm{D}}\beta}{4\pi m}\left(1+\gamma+\ln{\dfrac{ma_{2\mathrm{D}}^{2}}{2\hbar^{2}\beta}}\right)$
		$\displaystyle-E\beta+O(\beta^{2}).$		(43)

Setting $\beta=\eta+is$ where $\eta$ is a positive infinitesimal and $s$ is real, we see that the above equation shows the Fourier transform of $\rho(\epsilon)$ at small $s$, and this Fourier transform has a singular term proportional to $s\ln s$. This singular term is caused by a power law tail of the coarse-grained version of $\rho(\epsilon)$ at $\epsilon\to\infty$. Taking the inverse Fourier transform of this singular term, we find the power law tail shown in Eq. (5).

Applying $\frac{\mathrm{d}}{\mathrm{d}\beta}$ to both sides of Eq. (III), we find

$$E=\int_{-\infty}^{\infty}\rho(\epsilon)\epsilon e^{-\beta\epsilon}\mathrm{d}\epsilon-\dfrac{\hbar^{2}\mathcal{I}_{2\mathrm{D}}}{4\pi m}\left(\gamma+\ln{\dfrac{a_{2\mathrm{D}}^{2}m}{2\hbar^{2}\beta}}\right)+O(\beta).$$

(44)

We divide the domain of integration over $\epsilon$ into two regions: one is $(-\infty,\epsilon_{\mathrm{M}})$ and the other is $(\epsilon_{\mathrm{M}},\infty)$, where $\epsilon_{\mathrm{M}}$ is an energy scale such that $\epsilon_{\mathrm{M}}$ is very large but $\epsilon_{\mathrm{M}}\beta\ll 1$. In $(-\infty,\epsilon_{\mathrm{M}})$ we have

$$\int_{-\infty}^{\epsilon_{\mathrm{M}}}\rho(\epsilon)\epsilon e^{-\beta\epsilon}\mathrm{d}\epsilon\approx\int_{-\infty}^{\epsilon_{\mathrm{M}}}\rho(\epsilon)\epsilon\mathrm{d}\epsilon,$$

(45)

while in $(\epsilon_{\mathrm{M}},\infty)$ we use Eq. (5) to do the integral:

	$\displaystyle\int_{\epsilon_{\mathrm{M}}}^{\infty}\rho(\epsilon)\epsilon e^{-\beta\epsilon}\mathrm{d}\epsilon$	$\displaystyle=\dfrac{\hbar^{2}\mathcal{I}_{2\mathrm{D}}}{4\pi m}\int_{\epsilon_{\mathrm{M}}}^{\infty}\epsilon^{-1}e^{-\beta\epsilon}\mathrm{d}\epsilon$
		$\displaystyle=\dfrac{\hbar^{2}\mathcal{I}_{2\mathrm{D}}}{4\pi m}\mathrm{\Gamma}(0,\epsilon_{\mathrm{M}}\beta)$
		$\displaystyle=\dfrac{\hbar^{2}\mathcal{I}_{2\mathrm{D}}}{4\pi m}[-\gamma-\ln{(\epsilon_{\mathrm{M}}\beta)}]+O(\epsilon_{\mathrm{M}}\beta).$		(46)

Thus, taking $\beta\rightarrow 0$, we get

	$\displaystyle E$	$\displaystyle=\lim_{\epsilon_{\mathrm{M}}\rightarrow\infty}\left[\int_{-\infty}^{\epsilon_{\mathrm{M}}}\rho(\epsilon)\epsilon\mathrm{d}\epsilon-\dfrac{\hbar^{2}\mathcal{I}_{2\mathrm{D}}}{4\pi m}\ln{\dfrac{e^{2\gamma}ma_{2\mathrm{D}}^{2}\epsilon_{\mathrm{M}}}{2\hbar^{2}}}\right]$
		$\displaystyle=\lim_{\epsilon_{\mathrm{M}}\rightarrow\infty}\left(\sum_{\epsilon_{\nu}<\epsilon_{\mathrm{M}}}\epsilon_{\nu}n_{\nu}-\dfrac{\hbar^{2}\mathcal{I}_{2\mathrm{D}}}{4\pi m}\ln{\dfrac{e^{2\gamma}ma_{2\mathrm{D}}^{2}\epsilon_{\mathrm{M}}}{2\hbar^{2}}}\right),$		(47)

which is Eq. (3).

According to Eqs. (10), (29), and (30), we have

$$n_{\nu\sigma}=\int\text{d}^{2}r\int\text{d}^{2}b\,\phi_{\nu}(\mathbf{r})\phi_{\nu}^{*}(\mathbf{r}+\mathbf{b})p_{\sigma}(\mathbf{r},\mathbf{r}+\mathbf{b}).$$

(48)

When $\epsilon_{\nu}$ is large, the integrand as a function of $\mathbf{b}$ oscillates rapidly, which implies that the only important contribution is from the singular term in the expansion of $p_{\sigma}(\mathbf{r},\mathbf{r}+\mathbf{b})$ at $\mathbf{b}\rightarrow 0$ [1], and this singular term is $\frac{1}{8\pi}C_{2\mathrm{D}}(\mathbf{r})b^{2}\ln{\frac{b}{a_{2\mathrm{D}}e}}$. Since $\phi_{\nu}(\mathbf{r})$ satisfies the single-particle Schrödinger equation, Eq. (31), we have

$$\phi^{*}_{\nu}(\mathbf{r}+\mathbf{b})\approx(\hbar^{2}/2m\epsilon_{\nu})^{2}\nabla^{4}_{\mathbf{b}}\phi^{*}_{\nu}(\mathbf{r}+\mathbf{b})$$

(49)

with relative error $\sim O(\epsilon_{\nu}^{-1})$ at $b\sim\sqrt{\hbar^{2}/2m\epsilon_{\nu}}$. Substituting Eq. (49) into Eq. (48) and carrying out the integral over $\mathbf{b}$, we find

$$n_{\nu\sigma}=\frac{1}{k_{\nu}^{4}}\int\text{d}^{2}rC_{2\mathrm{D}}(\mathbf{r})|\phi_{\nu}(\mathbf{r})|^{2}+O(\epsilon_{\nu}^{-5/2}),$$

(50)

where $k_{\nu}=\sqrt{2m\epsilon_{\nu}}/\hbar$.

The calculation procedure in 1D is similar to the one in 2D. We define a normalized $N$-body energy eigenstate $\ket{\Phi}$ in 1D,

	$\displaystyle\ket{\Phi}=$	$\displaystyle(N_{\uparrow}!N_{\downarrow}!)^{-1/2}\int\tilde{D}_{1}^{\uparrow}\tilde{D}_{1}^{\downarrow}\Phi(x_{1}\dots x_{N_{\uparrow}}y_{1}\dots y_{N_{\downarrow}})$
		$\displaystyle\times\psi_{\uparrow}^{\dagger}(x_{1})\dots\psi_{\uparrow}^{\dagger}(x_{N_{\uparrow}})\psi_{\downarrow}^{\dagger}(y_{1})\dots\psi_{\downarrow}^{\dagger}(y_{N_{\downarrow}})\ket{0}.$		(51)

where $N=N_{\uparrow}+N_{\downarrow}$, $N_{\uparrow}$ is the number of spin-up fermions, $N_{\downarrow}$ is the number of spin-down fermions, $x_{1},\dots,x_{N_{\uparrow}}$ are the coordinates of the spin-up fermions, $y_{1},\dots,y_{N_{\downarrow}}$ are the coordinates of the spin-down fermions, $\psi_{\uparrow}^{\dagger}(x)$ is the creation operator for a spin-up fermion at position $x$, $\psi_{\downarrow}^{\dagger}(y)$ is the creation operator for a spin-down fermion at position $y$, $\tilde{D}_{i}^{\uparrow}\equiv\prod_{\mu=i}^{N_{\uparrow}}\mathrm{d}x_{\mu}$, $\tilde{D}_{i}^{\downarrow}\equiv\prod_{\mu=i}^{N_{\downarrow}}\mathrm{d}y_{\mu}$, and $\Phi(x_{1}\dots x_{N_{\uparrow}}y_{1}\dots y_{N_{\downarrow}})$ is the $N$-body wave function which is antisymmetric under the interchange of any two spin-up (spin-down) fermions. The 1D Bethe-Peierls boundary condition is

	$\displaystyle\ket{\Phi}$	$\displaystyle=$	$\displaystyle A\left(\frac{x_{1}+y_{1}}{2};x_{2}\dots x_{N_{\uparrow}}y_{2}\dots y_{N_{\downarrow}}\right)$		(52)
			$\displaystyle\times\left(1-\frac{|x_{1}-y_{1}|}{a_{1\mathrm{D}}}\right)+O(|x_{1}-y_{1}|^{2}),$

which is satisfied by the wave function when $x_{1}$ and $y_{1}$ are close. The one-particle density matrix for spin-$\sigma$ fermions in 1D is defined as

$$p_{\sigma}(x,x+b)=\bra{\Phi}\psi_{\sigma}^{\dagger}(x)\psi_{\sigma}(x+b)\ket{\Phi}.$$

(53)

For spin-up fermions, we substitute Eq. (IV) into the above definition and find

	$\displaystyle p_{\uparrow}(x,x+b)=$	$\displaystyle\,N_{\uparrow}\int\tilde{D}_{2}^{\uparrow}\tilde{D}_{1}^{\downarrow}\Phi^{*}(x,x_{2}\dots x_{N_{\uparrow}}y_{1}\dots y_{N_{\downarrow}})$
		$\displaystyle\times\Phi(x+b,x_{2}\dots x_{N_{\uparrow}}y_{1}\dots y_{N_{\downarrow}}).$		(54)

After finishing calculations analogous to those for the 2D one-particle density matrix, we find

	$\displaystyle p_{\uparrow}(x,x+b)=$	$\displaystyle n_{\uparrow}(x)+u_{\uparrow}(x)b+\frac{1}{2}v_{\uparrow}(x)b^{2}$
		$\displaystyle+C_{1\mathrm{D}}(x)\left(-\frac{b^{2}a_{1\mathrm{D}}}{4}+\frac{|b|^{3}}{12}\right)$
		$\displaystyle+w(x)\frac{2b^{3}}{3a_{1\mathrm{D}}}+w^{*}(x)\frac{b^{3}}{6a_{1\mathrm{D}}}+T^{\prime}_{b}+O(b^{4}),$		(55)

where

	$\displaystyle n_{\uparrow}(x)$	$\displaystyle=N_{\uparrow}\int\tilde{D}_{2}^{\uparrow}\tilde{D}_{1}^{\downarrow}|\Phi(x,\mathbf{X})|^{2},$		(56)
	$\displaystyle u_{\uparrow}(x)$	$\displaystyle=N_{\uparrow}\lim_{\eta\to 0}\int_{\mathcal{D}_{\eta}}\tilde{D}_{2}^{\uparrow}\tilde{D}_{1}^{\downarrow}\Phi^{*}(x,\mathbf{X})\frac{\partial}{\partial x}\Phi(x,\mathbf{X}),$		(57)
	$\displaystyle v_{\uparrow}(x)$	$\displaystyle=N_{\uparrow}\lim_{\eta\to 0}\int_{\mathcal{D}_{\eta}}\tilde{D}_{2}^{\uparrow}\tilde{D}_{1}^{\downarrow}\Phi^{*}(x,\mathbf{X})\dfrac{\partial^{2}}{\partial x^{2}}\Phi(x,\mathbf{X}),$		(58)
	$\displaystyle T^{\prime}_{b}$	$\displaystyle=\frac{N_{\uparrow}b^{3}}{3!}\lim_{\eta\to 0}\int_{\mathcal{D}_{\eta}}\tilde{D}_{2}^{\uparrow}\tilde{D}_{1}^{\downarrow}\Phi^{*}(x,\mathbf{X})\dfrac{\partial^{3}}{\partial x^{3}}\Phi(x,\mathbf{X}),$		(59)
	$\displaystyle C_{1\mathrm{D}}(x)$	$\displaystyle\equiv\frac{4N_{\uparrow}N_{\downarrow}}{a_{1\mathrm{D}}^{2}}\int\tilde{D}_{2}^{\uparrow}\tilde{D}_{2}^{\downarrow}\ |A(x;\mathbf{X}^{\prime})|^{2},$		(60)
	$\displaystyle w(x)$	$\displaystyle\equiv N_{\uparrow}N_{\downarrow}\int\tilde{D}_{2}^{\uparrow}\tilde{D}_{2}^{\downarrow}\ A^{*}(x;\mathbf{X}^{\prime})\frac{\partial A(x;\mathbf{X}^{\prime})}{\partial x},$		(61)

and $\mathbf{X}\equiv(x_{2}\dots x_{N_{\uparrow}}y_{1}\dots y_{N_{\downarrow}})$, $\mathbf{X}^{\prime}\equiv(x_{2}\dots x_{N_{\uparrow}}y_{2}\dots y_{N_{\downarrow}})$. $n_{\uparrow}(x)$ is the spatial density of spin-up fermions at position $x$, $C_{1\mathrm{D}}(x)$ is the 1D contact density, and $w(x)$ is related to the center-of-mass motion of small-distance pairs of fermions in different spin states. We can also find a similar expansion for $p_{\downarrow}(x,x+b)$.