A Hydrodynamical Description of Bose Liquid with Fractional Exclusion Statistics

Previously, it is discovered that for hamiltonian of the form as $\hat{H}=-\frac{1}{2}\sum_{i}\partial_{i}^{2}+V(x_{1},\cdots,x_{N})+\frac{1}{2}\sum_{i}\omega_{0}^{2}x_{i}^{2}$, the system has $SU(1,1)$ dynamical symmetry as long as $V(x_{1},\cdots,x_{N})$ satisfying $V(\lambda x_{1},\cdots,\lambda x_{N})=\lambda^{-2}V(x_{1},\cdots,x_{N})$. A simple review of this symmetry can be checked by the commutative relation among $\hat{H}_{0}=-\frac{1}{2}\sum_{i}\partial_{i}^{2}+V(x_{1},\cdots,x_{N})$, $\hat{I}=\frac{1}{2}\sum_{i}\omega^{2}_{0}x_{i}^{2}$, and $\hat{Q}=-\frac{i}{2}\sum_{i}(\partial_{i}x_{i}+x_{i}\partial_{i})$, i. e.,

$\displaystyle[\hat{Q},\hat{H}_{0}]=2i\hat{H}_{0},\hskip 17.22217pt[\hat{Q},\hat{I}]=-2i\hat{I},$

$\displaystyle[\hat{I},\hat{H}_{0}]=i\omega_{0}^{2}\hat{Q}.$

By introducing $L_{1}=(\hat{H}_{0}-\hat{I})/2\omega_{0}$, $L_{3}=\hat{H}/2\omega_{0}$, and $L_{2}=\hat{Q}/2$, we have

$\displaystyle[L_{1},L_{2}]=-iL_{3},\hskip 12.91663pt[L_{2},L_{3}]=iL_{1},\hskip 12.91663pt[L_{3},L_{1}]=iL_{2}.$

Clearly, these three operators form $\mathfrak{su}(1,1)$ algebra. We can check that both of the CS model and the trapped $|\Phi|^{6}$ theory have this scale property. For the $|\Phi|^{6}$ theory, its microscopic first quantization interaction term is a three-body contact interaction

$\displaystyle V(x_{1},\cdots,x_{N})=g_{0}^{2}\pi^{2}\sum_{i<j<k}\delta(x_{i}-x_{j})\delta(x_{j}-x_{k}).$

(3)

The dynamical symmetry is easily checked by noticing $\delta(\lambda x)=\delta(x)/|\lambda|$, which is also inherited by the trapped $|\Phi|^{6}$ theory.

The direct result for this dynamical symmetry is a $2\omega_{0}$ frequency breathing mode in the trap pitaev1997 . We also stress that the symmetry property is correct for any particle number $N$. But the symmetry alone can not ensure the correspondence between the two systems. Actually we can see the deviation between these two theories at large $g_{0}=g_{\rm cs}$ with a finite $N$, where the long-range interaction in the CS model induce a Wigner crystal wigner while this density order is clearly absent in the $|\Phi|^{6}$ theory.

In this section, we first verify that the ground state density profile in the CS model and trapped $|\Phi|^{6}$ theory are identical in the thermodynamic limit. Then, for a finite $N$, we present the numerical results of the ground state density $\rho^{\rm CS}_{\rm g}(x)=|\Psi_{\rm g}(x)|^{2}$, $\rho_{\rm g}(x)$ for the CS model and the $|\Phi|^{6}$ theory respectively.

In the $|\Phi|^{6}$ theory, the ground state density profile can be obtained by an equivalent hydrodynamical description, where the density and velocity field obey the following equations,

$\displaystyle\frac{\partial}{\partial t}\rho+\frac{\partial}{\partial x}(\rho v)=0,$

(4)

$\displaystyle\frac{\partial v}{\partial t}\!+\!\frac{\partial}{\partial x}\left(\frac{1}{2}v^{2}\!+\!\frac{1}{2}\omega_{0}^{2}x^{2}\!+\!\frac{g_{0}^{2}\pi^{2}\rho^{2}}{2}\!+\!P_{Q}\!\right)\!=0.$

(5)

Here, $\rho=\rho(x,t)\equiv|\Phi(x,t)|^{2}$ is the density field, and $v=v(x,t)\equiv\left(\Phi^{*}\partial_{x}\Phi-\partial_{x}\Phi^{*}\Phi\right)/2i\rho(x,t)$ is the velocity field. $P_{Q}=-\sqrt{\rho(x,t)}\partial^{2}_{x}\sqrt{\rho(x,t)}/2$ is the so-called quantum pressure. When the trap varies slowly on the scale of the inter-particle spacing, the density profile $\rho$ becomes smooth and the quantum pressure $P_{Q}$ can be neglected, which is called the Thomas-Fermi (TF) approximation. Because in ground state, the velocity field is zero, $\partial_{t}\rho=0$, then it requires

$\displaystyle\frac{\partial}{\partial x}\left(\omega_{0}^{2}x^{2}+g_{0}^{2}\pi^{2}\rho^{2}\right)=0.$

(6)

Therefore, the ground state density profile is

$\displaystyle\rho_{\rm g}(x)=\frac{\omega_{0}}{\pi g_{0}}\sqrt{\frac{2g_{0}N}{\omega_{0}}-x^{2}}\theta\left(\sqrt{\frac{2g_{0}N}{\omega_{0}}}-|x|\right),$

(7)

where $N=\int{\mathrm{d}x}\rho(x)$ is the total particle number, $\theta(x)$ denotes the Heaviside step function. When $g_{0}=1$, such an equation is the ground state profile of spinless ideal fermions in a harmonic trap, which has been studied in Ref. eugene2000 Although this result is obtained under TF approximation, we see the semicircle law is quickly reached for finite $N$ with our numerical simulations.

Next, we solve the ground state density profile for the CS model. The expression of the ground state of the CS model is given as

$\displaystyle\Psi_{\rm g}(\{x_{j}\})=\sqrt{C}\prod_{1\leq i<j\leq N}|x_{i}-x_{j}|^{g_{\rm cs}}e^{-\frac{\omega_{\rm cs}}{2}\sum_{j}x_{j}^{2}},$

(8)

with the ground state energy $E_{g}=\omega_{\rm cs}[N/2+g_{\rm cs}N(N-1)/2]$, where $x_{j}$ are the positions of particles, and $C$ is the normalization factor, whose inverse is

$\displaystyle C^{-1}=\left(\frac{\omega_{\rm cs}}{2}\right)^{-E_{g}/\omega_{\rm cs}}\left(\frac{\pi}{2}\right)^{N/2}\prod_{j=1}^{N}\frac{\Gamma(1+jg_{\rm cs})}{\Gamma(1+g_{\rm cs})}.$

(9)

The density profile is obtained by $\rho^{\rm CS}_{\rm g}(x)=\int\prod_{j=2}^{N}\mathrm{d}x_{j}|\Psi_{g}(x,x_{2},\cdots,x_{N})|^{2}$. As the wave function itself is just the joint probability of the eigenvalues for random Gaussian matrices ensemble, the density function can be viewed as the charge distribution of log-potential Coulomb gas forrester2010 . In the thermodynamic limit, the solution of the density profile is obtained

$\displaystyle\rho^{\rm CS}_{\rm g}(x)=\frac{\omega_{\rm cs}}{\pi g_{\rm cs}}\sqrt{\frac{2g_{\rm cs}N}{\omega_{\rm cs}}-x^{2}}\theta\left(\sqrt{\frac{2g_{\rm cs}N}{\omega_{\rm cs}}}-|x|\right).$

(10)

By comparing Eq. (7) and Eq. (10), we find that these two density profiles are exactly the same when we match the frequencies and the interaction strength. It is worth mentioning that $\rho_{\rm g}(x)$ also obtained from the zero-temperature distribution of the occupation number $n(\epsilon)$ for anyons satisfying fractional exclusion statistics murthy1994 :

$$\displaystyle\begin{split}n(\epsilon)=g_{0}^{-1}\theta(g_{0}N\omega_{0}-\epsilon)\end{split}$$

(11)

with a semiclassical single-particle energy $\epsilon=p^{2}/2+\omega_{0}^{2}x^{2}/2$. Then, $\rho_{\rm g}(x)$ is given by $\rho_{\rm g}(x)=(2\pi)^{-1}\int\mathrm{d}pn(\epsilon)$.

In the following, we are going to show the numerical results for ground state density profile for the CS model and for the $|\Phi|^{6}$ theory. Here and in the following figures, we set $\omega_{\rm cs}=\omega_{0}=1$. For the CS model, we use the Monte-Carlo integration. Here we keep $x_{1}=x$, and we reparametrize $x_{i=2}$ to $x_{i=N}$ on a $(N-1)$-sphere $\sum_{i=2}^{N}x_{i}^{2}=r^{2}$. The measure can be then written as $\prod_{j=2}^{N}\mathrm{d}x_{i}=r^{N-2}\mathrm{d}r\mathrm{d}\Omega_{N-1}$ with $\Omega_{N-1}$ denoting the solid angle of the $(N-1)$-sphere. The ground state density profile is then expressed as

	$\displaystyle\rho_{\rm g}^{\rm CS}(x)$	$\displaystyle=$	$\displaystyle\int r^{N-2+g_{\rm cs}(N-1)(N-2)}e^{-\omega_{\rm cs}r^{2}}\mathrm{d}r\int\mathrm{d}\Omega_{N-1}$
			$\displaystyle\prod_{j=2}^{N}|x\!-\!x_{j}(r,\Theta)|^{2g_{\rm cs}}\!\!\!\!\!\!\!\prod_{2\leq i<j\leq N}\!\!\!\!\!\!\!\left|\frac{x_{i}(r,\Theta)\!-\!x_{j}(r,\Theta)}{r}\right|^{2g_{\rm cs}}.$

Here $\Theta$ represent all the angle parameters of the $(N-1)$-sphere. From the expression, the integrant is semi-positive definite, thus the multi-dimensional integral is very suitable for Monte-Carlo integration. Meanwhile, because of the measure and the Gaussian distribution in wave-function, the major contribution of above integration comes from very narrow region of $r=0$. The concentration of measure in angular part due to Lévy’s Lemma Ledoux can further reduce the complexity in integration. With all these helps, the ground state density profile of CS model for $N=3$, $N=6$ can be obtained.

On the other hand, we carry out a numerical calculation of $|\Phi|^{6}$ theory as a non-linear Schrödinger equation. A cutoff in real space is taken as $X_{\rm max}$, and a boundary condition is taken as $\Phi(x=X_{\rm max})=\Phi(x=-X_{\rm max})$. To solve the ground state of the $|\Phi|^{6}$ theory, we turn on an evolution in imaginary time, where the ground state is the fixed point of this flow.

The numerical results are summarized in Fig. 1. In Fig. 1(a), (b) and (c), we have shown the ground state density profile of CS model, the $|\Phi|^{6}$ theory and the semicircle law with different interaction strength $g_{\rm cs}=g_{0}=0.5$, $1$ and $1.5$. For $g_{\rm cs}>1$, the density wave order becomes prominent. This is the tendency for Wigner crystal because of the long-range nature of $r^{-2}$ interaction. For $g_{\rm cs}=g_{0}<1$, the ground state density profile of CS model and the $|\Phi|^{6}$ model is very close to each other even for $N=6$. At the same time, in Fig. 1(d), (e) and (f), the ground state density profile of the $|\Phi|^{6}$ theory approaches the semicircle distribution as $N$ increases. Actually, the semicircle law is already a very good approximation for $N=6$.

In the following, we are going to connect the density-wave excitation of the $|\Phi|^{6}$ theory and the CS model in the thermodynamic limit. By comparison, we find that the density-wave excitations of the two models share the same expression. Then, we check the finite-$N$ effect by full numerical calculations of the time evolution of the excitations. We analytically prove the correspondence between the density-wave excitations for any excited modes where $g_{\rm cs}=g_{0}=1$, or for the second excited mode where $g_{\rm cs}=g_{0}$ is arbitrary. A guess for the correspondence between the density-wave excitations of the two models for any excited modes and for arbitrary $g_{\rm cs}=g_{0}$ is made and is checked with finite $N$. From the numerical results we obtained, the guess seems to be hold in the thermodynamic limit.

First of all, let us focus on the $|\Phi|^{6}$ theory. Here our starting point is the hydrodynamical equation Eq. (4) and Eq. (5). Following the procedure in Refs. string1996 ; menotti2002 , we linearize the density $\rho=\rho_{0}+\delta\rho$ (correspondingly $v=\delta v$) with a perturbation $\delta\rho$ obeying the linear expansion of Eqs. (4), (5)

$\displaystyle\frac{\partial^{2}\delta\rho}{\partial t^{2}}-\omega_{0}^{2}\frac{\partial}{\partial x}\left[\sqrt{R^{2}-x^{2}}\frac{\partial}{\partial x}(\sqrt{R^{2}-x^{2}}\delta\rho)\right]=0.$

(13)

where $R=\sqrt{2g_{0}N/\omega_{0}}$ is the radius of the semicircle distribution in the thermodynamic limit. The boundary conditions come from the absence of the particle flow at $x=\pm R$ imam2006

$\displaystyle\rho_{0}\delta v|_{x=\pm R}=0.$

(14)

Then, by performing the transformation $\vartheta=\arccos(x/R)$ ($0\leq\vartheta\leq\pi$), $\delta\tilde{\rho}=\sqrt{R^{2}-x^{2}}\delta\rho,\delta\tilde{v}=\rho_{0}\delta v$, we obtain a standard wave function

$\displaystyle\frac{\partial^{2}\delta\tilde{\rho}}{\partial t^{2}}-\omega^{2}\frac{\partial^{2}\delta\tilde{\rho}}{\partial\vartheta^{2}}=0,$

(15)

with the Neumann boundary condition

$\displaystyle\left.\frac{\partial\delta\tilde{\rho}}{\partial\vartheta}\right|_{\vartheta=0,\pi}=0.$

(16)

Thus, the solution of $\delta\tilde{\rho}$ and $\delta\tilde{v}$ consists of Fourier series

$\displaystyle\delta\tilde{\rho}_{n}(\vartheta,t)$

$\displaystyle=$

$\displaystyle A_{0}\!+\!\!\sum_{n=1}^{\infty}\!A_{n}\cos(\omega_{n}t+\varphi_{n})\cos(n\vartheta),$

(17)

$\displaystyle\delta\tilde{v}(\vartheta,t)=\delta v_{0}\!-\!\omega_{0}\!\!\!\sum_{n=1}^{\infty}\!A_{n}\sin(\omega_{n}t+\varphi_{n})\sin(n\vartheta),$

(18)

where $\omega_{n}=n\omega_{0}$ denotes the quantized excitation energy. This equidistant excitation energy is also consistent with the $\mathrm{SU}(1,1)$ dynamic symmetry. The constant $A_{0}$ as a perturbation actually implies a small derivation from the chemical potential.The constant $\delta v_{0}$ is set to be 0 with a Galilean transformation.

In the above calculation, we assumed that the spatial derivative part of Gross-Pitaevskii (GP) equation can be neglected. The validity of this assumption has to be checked. By GP equation, we mean $\Phi(x,t)$ satisfy a non-linear Schrödinger equation as follows,

$\displaystyle i\partial_{t}\Phi=\!\left(\!-\frac{\partial_{x}^{2}}{2}\!-\!\mu\!+\!\frac{\omega_{0}^{2}x^{2}}{2}\!+\!\frac{g_{0}^{2}\pi^{2}}{2}|\Phi|^{4}\!\right)\Phi,$

(19)

where $\mu$ denotes the chemical potential, $\Phi$ is short for $\Phi(x,t)$. A direct numerical calculation for the excited states from the GP equation is difficult. Here we carry out the check in a reversed way. We set the initial wave function as a superposition of the ground state and the $n$-th excited mode. Then, we check the time evolution of the density wave with this initial condition, where the information of the excited modes can be extracted. Approximately, we take the initial state as

$\displaystyle\Phi_{n}(x,0)=\sqrt{\rho_{g}(x)+0.1\cos[n\arccos(x/R)]}.$

(20)

The dynamics of $\rho_{n}(x,t)\equiv|\Phi_{n}(x,t)|^{2}$ can be calculated by the GP equation. We find the dynamics of $\rho_{n}(x,t)$ can be factorized into several harmonic oscillations,

	$\displaystyle\delta\rho_{n}(x,t)$	$\displaystyle\equiv$	$\displaystyle\rho_{n}(x,t)-\overline{\rho_{n}(x,t)}$		(21)
		$\displaystyle=$	$\displaystyle\sum_{m}A_{m}(x)\cos[\omega_{m}(x)t+\varphi_{m}(x)],$

where $\overline{\rho_{n}(x,t)}=\int_{0}^{T_{\text{max}}}\mathrm{d}t\rho_{n}(x,t)/T_{\text{max}}$ is the average value of $\rho_{n}(x,t)$. $T_{\text{max}}$ is the time cutoff, which is taken as $50\pi$ here. These data are shown in Fig. 2. By a nonlinear fitting of the dynamics, we obtain the values of $A_{m}(x)$, $\omega_{m}(x)$ and $\varphi_{m}(x)$. These parameters are $x$ dependent. If our result matches the theoretical predictions, $\omega_{m}(x)=\omega_{m}$ is $x$ independent. Furthermore, we check that $\omega_{m}$ is interaction $g_{0}$ independent. This property is exact in $N\rightarrow\infty$ limit, but it is approximate for finite $N$. Also, we check $\omega_{m}=m\omega_{1}$. If this equation fails to be true, there are soft modes in the system. The existence of these soft modes will lead to decoherence for the long-time evolution. As we find $\tilde{\omega}_{m}(x)=\omega_{m}(x)/m$ are all very close to $1$, therefore we will give the mean value of $\omega(x)$ as $\overline{\omega}_{m}=\int_{-X_{\rm max}}^{X_{\rm max}}\mathrm{d}x\tilde{\omega}_{m}(x)/(2X_{\rm max})$ and its variance $\Delta\omega_{n}^{2}$ as a function of $g_{0}$ and $N$ in Tables 1, 2. In Fig. 2, one can observe that the dynamics in (b) is more stable than in (a), this is manifested in the frequency difference between $\bar{\omega}_{2}$ and $\bar{\omega}_{4}$. The difference between $\bar{\omega}_{2}$ and $\bar{\omega}_{4}$ shrinks from $0.003$ to $0.0002$ when the particle increase from $N=40$ to $N=125$.

We also study the $x$ dependence in $A_{m}(x)$, which gives us the information about the eigenfunctions of the excited modes,

$\displaystyle A_{m}(x)\propto\frac{\delta\tilde{\rho}_{m}(x)}{\sqrt{1-(x/R)^{2}}}.$

(22)

We find these relations are satisfied both for $N=40$ and $N=125$. In Fig. 3(a) and (b), we show two typical $\delta\rho_{2}(x,t)$ at two different positions $x=0$ and $x=3.711$ ($x/R=0.3$). In Fig. 3 (c), we show $A_{2}(x)$, $A_{4}(x)$ and $A_{6}(x)$ for $g_{0}=1.4$ and $g_{0}=0.6$ with the reference of theoretical prediction of $\delta\tilde{\rho}_{2}(x)/\sqrt{1-(x/R)^{2}}$ and $\delta\tilde{\rho}_{4}(x)/\sqrt{1-(x/R)^{2}}$ ($A_{m}(x)$ are even function of $x$). One finds precise agreement between the numerical results and the theoretical predictions when $|x|<R$. When $|x|\approx R$, the derivation indicates the finite-$N$ effects. The eigenfunctions show $g_{0}$ independence up to a rescaling of $R$ ($R$ is in general $g_{0}$ dependent).

Secondly, we study the density-wave excitation of the CS model. Because the CS model can be mapped identically to a set of free harmonic oscillators with $E_{g}$ as the ground-state energy gurappa1999 , the excited energy is also equidistant. We consider the first excited state in the irreducible subspace involving the ground state, that is csm1969 ; gambar1975

$\displaystyle\Psi_{1}(\{x_{i}\})=\left(\frac{E_{g}}{\omega_{\rm cs}}\right)^{-1/2}\left(\frac{E_{g}}{\omega_{\rm cs}}-\omega_{\rm cs}\sum_{i=1}^{N}x_{i}^{2}\right)\Psi_{\rm g}.$

(23)

More generally, we have

$$\displaystyle\begin{split}\Psi_{m}(\{x_{i}\})=&\!\sqrt{mB\left(\frac{E_{g}}{\omega_{\rm cs}},m\right)}L_{m}^{\frac{E_{g}}{\omega_{\rm cs}}\!-\!1}\!\left(\omega_{\rm cs}\sum_{i=1}^{N}x_{i}^{2}\right)\\ &\times\Psi_{\rm g}(\{x_{i}\}),\end{split}$$

(24)

where the index $m$ in $\Psi_{m}$ refers to the $m$-th excited state in the subspace. Here $B(p,q)$ denotes the Beta-function, and $L_{m}^{\alpha}(x)$ denotes the generalized Laguerre polynomials. Notice that the $j$-th excited state’s energy is $E_{g}+2j\omega_{\rm cs}$.

Let $\Psi(t=0)=c_{0}\Psi_{\rm g}+c_{1}\Psi_{1}$ with $|c_{0}|^{2}+|c_{1}|^{2}=1$. We extract the information from the time evolution of the density $\rho_{1}^{\rm CS}(x,t)=\int\prod_{j=2}^{N}\mathrm{d}x_{j}|\Psi(t)|^{2}=\int\prod_{j=2}^{N}\mathrm{d}x_{j}[|c_{0}|^{2}\Psi_{\rm g}^{2}+|c_{1}|^{2}\Psi_{1}^{2}+(c_{0}c_{1}^{*}e^{2i\omega_{\rm cs}t}+c_{0}^{*}c_{1}e^{-2i\omega_{\rm cs}t})\Psi_{\rm g}\Psi_{1}]$. Since $\delta\rho$ is time dependent, we have $\delta\rho_{1}^{\rm CS}(x,t)=(c_{0}c_{1}^{*}e^{2i\omega_{\rm cs}t}+c_{0}^{*}c_{1}e^{-2i\omega_{\rm cs}t})\delta\rho_{1}^{\rm CS}(x)$, where

			$\displaystyle\delta\rho_{1}^{\rm CS}(x)=\int_{-\infty}^{\infty}\prod_{j=2}^{N}\mathrm{d}x_{j}\Psi_{\rm g}(x,\{x_{j}\})\Psi_{1}(x,\{x_{j}\})$		(25)
		$\displaystyle=$	$\displaystyle\int_{-\infty}^{\infty}\prod_{j=2}^{N}\mathrm{d}x_{j}C\left(\frac{E_{g}}{\omega_{\rm cs}}-\omega_{\rm cs}\sum_{i=1}^{N}x_{i}^{2}\right)C^{-1}\Psi_{\rm g}(x,\{x_{j}\})^{2}$
		$\displaystyle=$	$\displaystyle\int_{-\infty}^{\infty}\prod_{j=2}^{N}\mathrm{d}x_{j}C\left(\frac{E_{g}}{\omega_{\rm cs}}+\omega_{\rm cs}\partial_{\omega_{\rm cs}}\right)\left[C^{-1}\Psi_{\rm g}(x,\{x_{j}\})^{2}\right]$
		$\displaystyle=$	$\displaystyle C\left(\frac{E_{g}}{\omega_{\rm cs}}-\omega_{\rm cs}\frac{\partial}{\partial\omega_{\rm cs}}\right)\left[C^{-1}\rho_{\rm g}^{\rm CS}(x)\right]$
		$\displaystyle\propto$	$\displaystyle\frac{\cos[2\arccos(x/R)]}{\sqrt{1-(x/R)^{2}}},$

which is consistent with Eq. (17) for $n=2$. In the derivation, we have assumed that the partial derivative $\partial_{\omega_{\rm cs}}$ and the thermodynamic limit are interchangeable. However, for higher excited states which involves higher partial derivatives, this assumption is invalid due to the divergent derivative of $\rho_{\rm g}^{\rm CS}(x)$ at the boundary $x=\pm R$. Although we can not prove the general form of $\delta\rho_{m}(x)$, we can guess its form is consistent with Eq. (17), i. e.,

	$\displaystyle\delta\rho_{m}^{\rm CS}(x)$	$\displaystyle=$	$\displaystyle\int_{-\infty}^{\infty}\prod_{j=2}^{N}\mathrm{d}x_{j}\Psi_{\rm g}(\{x_{j}\})\Psi_{m}(\{x_{j}\})$		(26)
		$\displaystyle\propto$	$\displaystyle\frac{\cos[2m\arccos(x/R)]}{\sqrt{1-(x/R)^{2}}}$

For an arbitrary $g_{\rm cs}$, it is too difficult to solve the overlap integral for an any excited state, but Eq. (26) can be proved when $g_{\rm cs}=1$. Here, one of the excited states can be represented by the single-particle energy eigenstates of the harmonic oscillator $\{f_{n}(x)\}_{n=0,1,2,\cdots.}$. Let $|0,1,\cdots,N-1\rangle$ denote the occupied single-particle energy eigenstates in the ground state of $N$ fermions. Then, one of the $m$-th excited states $|\psi_{m}\rangle=\hat{A}\hat{b}^{\dagger}_{N-l+m}\hat{b}_{N-l}|0,1,\cdots,N-1\rangle$, for some $1\leq l\leq m$, where $\hat{A}$ denotes the unit antisymmetric function girardeau1960 , $\hat{b}^{\dagger}_{l}$ ($\hat{b}_{l}$) denotes the creation (annihilation) operator of $l$ single-particle energy eigenstates. The corresponding perturbation of the density is given by (notice $\hat{A}^{2}=1$)

	$\displaystyle\delta\rho_{m}^{\rm CS}(x)$	$\displaystyle=$	$\displaystyle\int_{-\infty}^{\infty}\prod_{j=2}^{N}\mathrm{d}x_{i}\psi_{0}(x,\{x_{j}\})\psi_{m}(x,\{x_{j}\})$		(27)
		$\displaystyle=$	$\displaystyle f_{N-l}(x)f_{N-l+m}(x).$

Using the asymptotic expansion of the product in the thermodynamic limit (appendix A):

$$\displaystyle\begin{split}&f_{N-l}\left(\sqrt{\frac{2N}{\omega_{\rm cs}}}x\right)f_{N-l+m}\left(\sqrt{\frac{2N}{\omega_{\rm cs}}}x\right)\\ &\sim\sqrt{\frac{\omega_{\rm cs}}{2\pi^{2}N}}\frac{\cos(m\arccos x)}{\sqrt{1-x^{2}}}.\end{split}$$

(28)

Comparing this result with Eq. (17), we prove that the density-wave excitations in the trapped $|\Phi|^{6}$ theory and the CS model are the same when $g_{\rm cs}=1$.

Now, we are going to check Eq. (26) for finite-$N$ systems by numerics. The overlap integral is numerically done by using Monte-Carlo method because the integrant is semi-positive definite. Here we introduce a normalization condition $\int\mathrm{d}x\delta\rho_{m}(x)^{2}\rho_{\rm g}^{\rm CS}(x)=\pi/2$ for the comparison between different results, where $\delta\rho_{m}(x)$ and $\rho_{\rm g}^{\rm CS}(x)$ are numerically obtained with a finite $N$. In Fig. 4, we present $\delta\rho_{m=1,2,3}(x)$ for $g_{\rm cs}=0.5$. There are datas for two sizes: $N=3$ and $N=6$. One can see the tendency for $\delta\rho_{m}(x)$ approaching the guessed function in the thermodynamic limit. The results for $g_{\rm cs}=1.5$ are given in (d), (e) and (f). There are clear deviations from the guessed results due to the high-frequency oscillation. The coarse-graining line is taken by replacing $\delta\rho_{m}(x)$ as $\int_{x-\Delta x}^{x+\Delta x}\mathrm{d}x^{\prime}\delta\rho_{m}(x^{\prime})/(2\Delta x)$. Here $\Delta x$ is taken as $0.25$. The purple dashed lines are results after two times of coarse graining, which can be viewed as the average behavior of $\delta\rho_{m}(x)$. One can see their resemblance with the guessed result. From the numerical result we find when $g_{\rm cs}\gg 1$, then the deviation from the theoretical prediction is more serious. In the CS model, when $g_{\rm cs}$ is large, the particles are inclined to be aligned in equal space in one-dimension due to the long-range interaction. On the other hand, when $g_{\rm cs}<1$, $\vartheta=\arccos(x/R)$ characterizes the oscillation of the density wave. But the oscillation in $\vartheta$ coordinate is not consistent with the equal spacing in $x$ coordinate. Thus these two requirements are conflicting with each other.For this reason, when $g_{\rm cs}/N\rightarrow 0$ is a suitable requirement for our guess.

The wave function Eq. (15) implies the canonical quantization of the density-wave excitation of the trapped $|\Phi|^{6}$ theory. Following the procedure of the phenomenological bosonization in Ref. giamarchi2003 , we introduce a scalar field $\phi$ and its conjugate momentum $\Pi$ through $\delta\tilde{\rho}=-\partial\phi/\partial\vartheta,\delta\tilde{v}=\omega(g\pi)^{-1}\Pi$. Then, using the solution (Eqs. 17, 18) and to the second order of $\delta\rho$, the energy functional (Eq. 2) is rewritten in terms of $\phi$ and $\Pi$ as

$\displaystyle H=\frac{1}{2}N^{2}g_{0}\omega_{0}+\frac{1}{2}\int_{0}^{\pi}\mathrm{d}\vartheta\left[\frac{\omega_{0}}{g_{0}\pi}\Pi^{2}+\pi g_{0}\omega_{0}\left(\frac{\partial\phi}{\partial\vartheta}\right)^{2}\right],$

(29)

which describes a Luttinger-liquid system. In the derivation, we neglect the term $\int\mathrm{d}x(\partial\rho_{0}/\partial x)^{2}$ within the TF approximation.

We quantize the scalar field with the commutation relation $[\hat{\phi}(t,y),\hat{\Pi}(t,y^{\prime})]=i\delta(y-y^{\prime})$. Thus, the density wave (Eq. 17, 18) is of bosonic-excitation type. It follows from Eqs. (17), (18) that

$\displaystyle\hat{H}=\frac{1}{2}g_{0}N^{2}\omega_{0}+\sum_{n=1}^{\infty}n\omega_{0}\hat{a}_{n}^{\dagger}\hat{a}_{n},$

(30)

where $\partial\hat{\phi}/\partial\vartheta=\sqrt{2/\pi}\sum_{n}(2g_{0}\pi n)^{-1/2}(\hat{a}_{n}^{\dagger}+\hat{a}_{n})\cos(n\vartheta)$, $\hat{\Pi}=\sum_{n}i(g_{0}\pi n/2)^{1/2}(\hat{a}_{n}^{\dagger}-\hat{a}_{n})\sqrt{2/\pi}\sin(n\vartheta)$ with $\{\hat{a}_{n}^{\dagger},\hat{a}_{n}\}$ are the creation and annihilation operators of the density-wave modes. Compared with the energy spectrum of the CS model gurappa1999 , Eq. (30) is also a set of free harmonic oscillators with $\frac{1}{2}g_{0}N^{2}\omega_{0}$ as the ground-state energy. Because, the only difference between the two models is just a negligible chemical potential deviation $(1-g_{0})\omega_{0}/2$ for $\omega_{0}\to 0$ in the thermodynamic limit, we conclude that the trapped $|\Phi|^{6}$ theory has the CS model-like energy spectra and thus obeys the fractional exclusion statistics with $g_{0}$ as the statistical parameter. About the connection between the Luttinger-liquid and the fractional exclusion statistics, we introduce Refs. Pham1998 ; Wu2001 for further reading .

The linearized density-wave excitations capture the low-energy excitations of the $|\Phi|^{6}$ theory. When the temperature $k_{\text{B}}T$ ($k_{\text{B}}$ denotes the Boltzmann constant) of the field grows up to $N\hbar\omega$, the thermal-excitation energy in Eq. (30) is comparable to the ground-state energy. At this moment, the linearized theory fails and we use the semiclassical approximation to study the finite-temperature property of the $|\Phi|^{6}$ theory in the following.

Using the Hatree-Fock approximation and the WKB approximation huse1982 ; giorgini1997 ; bhaduri2000 , the density $\rho(x)$ is determined self-consistently. First there should no condensate part in $\rho(x)$ at a high temperature. The system consists of thermal particles and the density $\rho(x)$ at finite temperature $\beta=1/k_{\text{B}}T$ reads

$\displaystyle\rho(x)=\int_{-\infty}^{\infty}\frac{\mathrm{d}p/(2\pi\hbar)}{e^{\beta\varepsilon(p,x)}-1},$

(31)

where

$\displaystyle\varepsilon(p,x)=\frac{p^{2}}{2}+\frac{1}{2}\omega_{0}^{2}x^{2}+\frac{g_{0}^{2}\pi^{2}\hbar^{2}}{2}\rho(x)^{2}-\mu$

(32)

denotes the effective single-particle spectrum, and $\mu$ denotes the chemical potential. After performing the $p$-integration of Eq. (31), we find the self-consistent equation of $\rho(x)$:

$\displaystyle u(x)=\mathrm{Li}_{1/2}[e^{\beta\mu-\beta\omega_{0}^{2}x^{2}/2-\pi g_{0}^{2}u(x)^{2}/4}],$

(33)

where $u(x)=\sqrt{2\pi\beta}\rho(x)$ denotes the dimensionless density, and $\mathrm{Li}_{s}(z)=\sum_{n=1}^{\infty}z^{n}/n^{s}$ denotes the polylogarithm function. Also, the constraint is given by

$\displaystyle N\sqrt{2\pi\beta}=\int_{-\infty}^{\infty}\mathrm{d}xu(x).$

(34)

Specially, we find the solution of $\rho(x)$ can always be obtained self-consistently for any temperature, which indicates that there are no condensate in $\rho(x)$ at any temperature.

In Fig. 5, we show the density $\rho(x)$ numerically solved from Eq. (33). Fig. 5a indicates three typical regions: (1) semicircle distribution: $|x|<R$; (2) Gaussian tail: $|x|>R$; (3) near the boundary $|x|\approx R$. Actually, the region 1 and 3 respectively corresponds to the two limits of the polylogarithm function: $\mathrm{Li}_{1/2}(z)\to\sqrt{\pi/(1-z)}$ for $z\to 1^{-}$, and $\mathrm{Li}_{1/2}(z)\to z$ for $z\to 0$. Accordingly, the density in the two regions are

$\displaystyle u^{1}(x)=\sqrt{\frac{(2\mu-\omega_{0}^{2}x^{2})}{\pi g_{0}^{2}k_{\text{B}}T}}\sqrt{1+\sqrt{1+\left(\frac{2\pi g_{0}k_{\text{B}}T}{2\mu-\omega_{0}^{2}x^{2}}\right)^{2}}},$

(35)

$\displaystyle u^{3}(x)=e^{\beta(\mu-\omega_{0}^{2}x^{2}/2)}.$

(36)

Substituting Eq. (35), or Eq. (36) into Eq. (32), we have

$\displaystyle\varepsilon^{1}(p,x)=\frac{p^{2}}{2}+\frac{(\pi g_{0}k_{\text{B}}T)^{2}}{4\mu-2\omega_{0}^{2}x^{2}}$

(37)

which results in the semicircle distribution of the density and implies the existence of the Fermi surface, or

$\displaystyle\varepsilon^{3}(p,x)=\frac{p^{2}}{2}+\frac{1}{2}\omega_{0}^{2}x^{2}-\mu$

(38)

which describes the ideal Bosons, but now the chemical potential contains the effects of interaction. In addition, near the boundary $|x|\approx R$ with a length scale $Rk_{\text{B}}T/\mu$, the analytical solutions (Eqs. 35, 36) are invalid. Fig. 5b compares the density of the $|\Phi|^{6}$ theory ($g_{0}=1$) with it of a free-fermion system for fixed $N$ and $T$. The narrower distribution of the density from the $|\Phi|^{6}$ theory indicates that more particles occupy the low-energy states, which lead to a lower chemical potential.