Spin-incoherent Luttinger liquid of one-dimensional spin-1 Tonks-Girardeau Bose gas: Spin-dependent properties

TL;DR

This paper investigates the spin-dependent properties of a one-dimensional spin-1 Bose gas in the spin-incoherent Luttinger liquid regime, revealing how spin incoherence affects momentum distributions and correlations.

Contribution

It provides analytical wave functions and numerical results for the spin-dependent momentum distributions in the Tonks-Girardeau limit, highlighting differences from spinless Bose gases.

Findings

Momentum distributions are broadened in the SILL regime.

At large N, spin components' distributions converge.

Momentum distributions follow a 1/p^4 decay with reduced coefficients.

Abstract

Spin-incoherent Luttinger liquid (SILL) is a different universal class from the Luttinger liquid.\ This difference results from the spin incoherence of the system when the thermal energy of the system is higher than the spin excitation energy.\ We consider one-dimensional spin-$1$ Bose gas in the SILL regime and investigate its spin-dependent many-body properties.\ In Tonks-Girardeau limit, we are able to write down the general wave functions in a harmonic trap.\ We numerically calculate the spin-dependent (spin-plus, minus, and $0$) momentum distributions in the sector of zero magnetization which allows to demonstrate the most significant spin-incoherent feature compared to the spinless or spin-polarized case.\ In contrast to the spinless Bose gas, the momentum distributions are broadened and in the large momentum limit follow the same asymptotic $1/p^4$ dependence but with reduced coefficients.\ While the density matrices and momentum distributions differ between different spin components for small $N$, at large $N$ they approach each other.\ We show these by analytic arguments and numerical calculations up to $N$ $=$ $16$.

Equations

$$|\Psi ⟩=\sum _{s_1,s_2,...s_N}{\psi }_{s_1,s_2,...s_N}(\vec{x})|s_1,s_2,...,s_N⟩,$$

$${\psi }_{\vec{n}}^{sym}(\vec{x})$$

$${\varphi }_n(y)$$

$$x<x_2<x_3,⟨(+,s_2,s_3)|\chi ⟩,$$

$$x_2<x<x_3,⟨(s_2,+,s_3)|\chi ⟩=⟨P_{12}(+,s_2,s_3)|\chi ⟩,$$

$$x_2<x_3<x,⟨(s_2,s_3,+)|\chi ⟩=⟨P_{123}(+,s_2,s_3)|\chi ⟩,$$

$${\rho }_{+}(x,x^{\prime })=N\sum _{{\vec{s}}^{\prime }}\int d\bar{x}{\psi }_{+,{\vec{s}}^{\prime }}^{\ast }(x,\bar{x}){\psi }_{+,{\vec{s}}^{\prime }}(x^{\prime },\bar{x}),$$

$${\rho }_{+}(x<x^{\prime })=$$

$$(E,P_{12}{)}_{+}=\sum _{s_2,s_3}⟨E(+,s_2,s_3)|\chi ⟩⟨P_{12}(+,s_2,s_3)|\chi ⟩.$$

$$(P_{\mathrm{12...}j},P_{\mathrm{12...}k}{)}_{+}=\frac{{\sum }_{{\vec{s}}^{\prime }}⟨P_{\mathrm{12...}j}(+,{\vec{s}}^{\prime })|P_{\mathrm{12...}k}(+,{\vec{s}}^{\prime })⟩}{{\mathrm{Tr}}_{\chi }(\mathrm{E})},$$

$${\rho }_{+}(x<x^{\prime })=$$

$$(P_{m-1,m}...P_{23}{P}_{12}{)}^{-1}{P}_{m...m+j-1}(P_{m-1,m}...P_{23}{P}_{12})$$

$$=P_{12}^{-1}{P}_{23}^{-1}...P_{m-1,m}^{-1}{P}_{m...m+j-1}{P}_{m-1,m}...P_{23}{P}_{12},$$

$$=P_{1,m+\mathrm{1...}m+j-1},$$

$$w_N\equiv {\mathrm{Tr}}_{\chi }(E)=\sum _{n=0}^{\frac{N}{2}\mathrm{or}\frac{N-1}{2}}\frac{N!}{(n!{)}^2(N-2n)!},$$

$$|+\underset{n-1}{\underbrace{+...+}}\underset{n}{\underbrace{-...-}}\underset{N-2n}{\underbrace{\mathrm{00...0}}}⟩.$$

$$(E,E{)}_{+}$$

$$|\underset{j}{\underbrace{+...+}}\underset{n-j}{\underbrace{+...+}}\underset{n}{\underbrace{-...-}}\underset{N-2n}{\underbrace{\mathrm{00...0}}}⟩,$$

$$(E,P_{\mathrm{12...}j}{)}_{+}$$

$$(E,E{)}_0$$

$$(E,P_{\mathrm{12...}j}{)}_0$$

$$2(E,E{)}_{+}+(E,E{)}_0=1.$$

$$2(E,P_{\mathrm{12...}j}{)}_{+}+(E,P_{\mathrm{12...}j}{)}_0=w_{jN}/w_N,$$

$$w_{jN}$$

$${\rho }_{\pm (0)}(p)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }dx{\int }_{-\infty }^{\infty }d{x}^{\prime }{e}^{ip(x-x^{\prime })}{\rho }_{\pm (0)}(x,x^{\prime }),$$

$${\rho }_{\pm (0)}(p)\underset{p\to \infty }{=}$$

$$(E,E{)}_{+}$$

$$(E,E{)}_0$$

$$f_k^{(N)}\equiv \sum _{j=0}^{\frac{N-k}{2}}\frac{N!}{(k+j)!j!(N-2j-k)!}.$$

$$\left|\frac{(E,P_{\mathrm{12...}j}{)}_{+}-(E,P_{\mathrm{12...}j}{)}_0}{(E,P_{\mathrm{12...}j}{)}_0}\right|,$$

$$\left|\frac{{\bar{f}}_j^{(N-j)}-{\bar{f}}_0^{(N-j)}}{{\bar{f}}_0^{(N-j)}}\right|,$$

$$\displaystyle\int dx\left|\begin{array}[]{cc}\phi_{n_{i}}^{\prime}(x)&\phi_{n_{j}}^{\prime}(x)\\ \phi_{n_{i}}(x)&\phi_{n_{j}}(x)\end{array}\right|^{2}$$

$$\displaystyle=\frac{1}{L^{2}}\int dx\left|\begin{array}[]{cc}ik_{i}e^{ik_{i}x}&ik_{j}e^{ik_{j}x}\\ e^{ik_{i}x}&e^{ik_{j}x}\end{array}\right|^{2},$$

$$=\frac{(k_i-k_j{)}^2}{L}.$$

Full text

I Introduction

A plethora of studies on one-dimensional (1D) quantum systems Giamarchi2004 ; Haldane1981 of gaseous atoms thrive recently due to the experimental achievements of 1D confined bosons Paredes2004 ; Kinoshita2004 ; Haller2009 . Many studies focus on the ground state properties of spinless bosons Girardeau2001 ; Papenbrock2003 such as spatial and momentum distributions Minguzzi2002 ; Olshanii2003 ; Xu2015 , quantum magnetism in a spinful Bose gas Deuretzbacher2008 ; Deuretzbacher2014 ; Volosniev2014 ; Yang2015 ; Yang2016 ; Deuretzbacher2016 , and low-energy excitations in the Luttinger liquid model. Meanwhile a spinful quantum system in the spin-incoherent regime Fiete2007 also provides a new avenue for studying 1D quantum many-body systems. This regime is termed as spin-incoherent Luttinger liquid (SILL) which forms a different universality class from the Luttinger liquid, where the temperature is high enough that different spin configurations can be regarded as degenerate while low enough that charge excitation is suppressed. For 1D spin-$1$ Bose gas Ho1998 ; Olshanii1998 with $s$-wave scattering lengths satisfying $|a_{0}-a_{2}|$ $\ll$ $a_{0,2}$, there exists a window of temperature for the gas in SILL regime. This happens since the sound velocity is much larger than the spin velocity.

In the crossover regime between Luttinger liquid and SILL, 1D fermions with tunable spins Pagano2014 and their high momentum tails Decamp2016 have been studied, which show an evident broadening in the momentum distributions Cheianov2005 ; Feiguin2010 . Quantum criticality Hazzard2013 and Pomeranchuk effect Zhou2014 in the spin-incoherent regime are also theoretically predicted in the two-dimensional Hubbard model. Here in contrast we investigate 1D spin-1 Bose gas Cazalilla2011 ; Stamperkurn2013 in the SILL regime in a harmonic trap, which is studied only recently Jen2016_spin1 . We shall focus on the Tonks-Girardeau (TG) regime Tonks1936 ; Girardeau1960 where the density is sufficiently low that the effective repulsion between particles can be regarded as infinite. TG spinor Bose gas is a special case of SILL since the exchange energy vanishes in this limit Volosniev2014 ; Yang2015 . Therefore TG gas automatically is in the regime of SILL. In TG gas limit, we can write down the exact spatial wave functions since bosons are fermionized and impenetrable due to effectively infinitely-strong atom-atom interactions. We then numerically calculate the momentum distributions for the three individual components of the spin-1 Bose gas (spin-plus, minus, and [math]). These predictions can be measurable in spin-resolved matter-wave experiments, either the time-of-flight experiment Shvarchuck2002 ; Davis2012 ; Jacqmin2012 ; Fang2016 or Bragg scattering spectroscopy Kozuma1999 ; Stenger1999 ; Steinhauer2002 ; Papp2008 ; Veeravalli2008 ; Pino2011 . This system allows for better demonstrations of SILL physics which is within reach of present experimental conditions. As compared with electronic spin-$1/2$ systems Fiete2007 ; Feiguin2010 , ultracold atom experiments not only provide controllable spatial dimensions but also tunable atom-atom interactions via Feshbach resonances, thus making our investigations testable in the quantum many-body systems.

In Ref. Jen2016_spin1 , we have derived the wave functions and density matrix for 1D spin-$1$ Bose gas in TG limit. We numerically calculate its momentum distributions, summed over spin components, up to six bosons. The momentum distributions are uniformly broadened as the number of bosons $N$ grows. We have also derived the analytical large momentum ($p$) asymptotic in one-body momentum distributions, which shows the universal $1/p^{4}$ dependence. The coefficients of the asymptotic $1/p^{4}$ are also formulated for arbitrary $N$. Here we present the spin-dependent properties of density matrix in 1D spin-$1$ Bose gas in TG limit, and show the spin-dependent momentum distributions up to $N$ $=$ $16$. We also obtain the spin-dependent coefficients of $1/p^{4}$ for asymptotic large $p$. Though the momentum distributions vary between different spin components for small $N$, they approach each other as $N$ increases. We show this from the numerical results accompanied by analytical arguments in the large $N$ limit.

The rest of the paper is organized as follows. In Sec. II we introduce the general wave functions for 1D spin-1 Bose gas. In Sec. III, we derive the general forms of density matrices for each spin components with individual spin function overlaps in SILL regime, and present the numerically calculated results using Monte Carlo integration method implemented with Gaussian unitary ensemble. In Sec. IV. we discuss the analytical derivation of high momentum asymptotic for each component, which we compare with numerically calculated momentum distributions. We also investigate the momentum distributions in large $N$ limit using the method of steepest descent or stationery phase, and compare them with the numerical results. Finally we conclude in Sec. V.

II General wave functions in TG limit

In general we can express the wave function of $N$ bosons as

[TABLE]

where we denote $\vec{x}$ $=$ $(x_{1},x_{2},...,x_{N})$ and $|s_{1},s_{2},...s_{N}\rangle$ $\equiv$ $|\vec{s}\rangle$ as the spatial distributions and the spin configurations respectively. Here we can label $s_{i}$ $=$ $+$, $-$, or [math], respectively for spin-plus, minus, and zero components for the $i$th particle. The total wave function must satisfy the bosonic symmetry, therefore it is sufficient to just consider the ordered region of $x_{1}$ $<$ $x_{2}$ $<$ $...$ $<$ $x_{N}$, and we can obtain all other regions via permutations of this ordered region. In TG gas limit, the atoms become fermionized that their spatial wave functions take the Slater determinant form of noninteracting fermions. For the symmetrized spatial part of the wave function, we denote it as $\psi_{\vec{n}}^{sym}(\vec{x})$ which can be expressed in terms of the eigenfunctions $\phi_{n_{j}}(x_{j})$ of the noninteracting fermions in a harmonic trap,

[TABLE]

We denote the sign function as sgn and the anti-symmetrizer as $\mathbb{A}$ for later convenience. The orbital indices are ($n_{1}$, $n_{2}$, $...n_{N}$), and the prefactor $\sqrt{N!}$ normalizes the wave function. For convenience we use the dimensionless forms of the eigenfunctions $\phi_{n}(y)$,

[TABLE]

where $H_{n}$ are Hermite polynomials. The harmonic oscillator length is $x_{ho}$ $\equiv$ $\sqrt{\hbar/(M\omega)}$ where $\omega$ is the trap frequency and $M$ is the atomic mass.

To eventually evaluate the density matrix for say the ”$+$” component, we need to obtain the wave function amplitude where at least one particle has spin ”$+$”. First we consider some degenerate and normalized spin configuration state $|\chi\rangle$ in some sector of magnetization, and the wave function can be expressed as $|\Psi\rangle$ $=$ $\psi_{\vec{n}}^{sym}(\vec{x})|\chi\rangle$. Take $N$ $=$ $3$ for an example, we obtain the probability amplitude $\psi_{+,s_{2},s_{3}}^{sym}(x,x_{2},x_{3})$ for the first particle having spin ”$+$” when we project $|\Psi\rangle$ to $\langle s_{1}=+,x_{1}=x|$ in the ordered region of $x$ $<$ $x_{2}$ $<$ $x_{3}$. To access the probability amplitudes in the other regions, we use the permutation operators $P_{12}$ and $P_{123}$ on the projected states, obtaining

[TABLE]

where we have suppressed the common $\psi_{\vec{n}}^{sym}(\vec{x})$ factors. Similar construction applies for other $N$’s. In the next section we proceed to calculate the spin-dependent density matrices for spin-1 Bose gas in the SILL regime.

III Density matrices for SILL of spin-1 Bose gas

The spin-dependent single-particle density matrix can be straightforwardly written down from the wave function described above. For example of the spin-plus component, we have

[TABLE]

where $\bar{x}$ $\equiv$ $(x_{2},x_{3},...,x_{N})$ and $\vec{s}^{\prime}$ $\equiv$ $(s_{2},s_{3},...,s_{N})$. A factor of $N$ represents $N$ possible choices of $x$ and $x^{\prime}$.

Again we take $N$ $=$ $3$ as an example, and consider only the region of $x$ $<$ $x^{\prime}$ which is symmetric to $x$ $>$ $x^{\prime}$. The spin-plus single-particle density matrix for $N$ $=$ $3$ then becomes

[TABLE]

where the parentheses $()$ in various integral regions represent the spin function overlaps. Take $(E,P_{12})_{+}$ as an example where $E$ is the identical permutation operator, we define

[TABLE]

Similar forms apply to the other spin function overlaps in Eq. (6). Also the factor of $2$ in Eq. (6) comes from the contribution of the integral region $x_{2}$ $>$ $x_{3}$, where the spin function overlaps are the same as those with $x_{2}$ $<$ $x_{3}$.

In the SILL regime, we average the above individual spin function overlaps by the total number of spin state configurations, which is denoted as Tr${}_{\chi}(E)$ $\equiv$ $\sum_{\chi}$ $\langle\chi|E|\chi\rangle$. It is simply the trace (Tr) of the identical operator over all spin configurations $|\chi\rangle$ since $\langle\chi|E|\chi\rangle$ $=$ $1$. We define the normalized spin function overlap in general as (using the same notation as the non-normalized one for simplicity)

[TABLE]

where $P_{12...j}$ are $j$-particle permutation operators in the symmetric group $S_{N}$. To derive Eq. (8), we have used the identity $\sum_{\chi}$ $|\chi\rangle\langle\chi|$ $=$ $1$. $(P_{12...j},P_{12...k})_{+}$ in general represents the spin function overlap from the integration region where the particle at $x$ permutes to just behind $x_{j}$ while the particle at $x^{\prime}$ permutes to just behind $x_{k}$.

In general for arbitrary $N$, we obtain the spin-plus density matrix as

[TABLE]

where we have used the properties of $(P_{12...j},P_{12...j})_{+}$ $=$ $(E,E)_{+}$ and $(E,P_{12...j})_{+}$ $=$ $(P_{12...m},P_{12...m+j-1})_{+}$ for $j,m$ $\geq$ $2$. The first property can be proved from Eq. (8) by using $P_{12...j}^{-1}P_{12...j}$ $=$ $E^{-1}E$ $=$ $1$. To prove the second property, we can reduce $P_{12...m}^{-1}P_{12...m+j-1}$ to

[TABLE]

such that $(E,P_{1,m+1...m+j-1})_{+}$, using again Eq. (8), is exactly the same as $(E,P_{12...j})_{+}$.

Other spin components of the density matrices, $\rho_{-}(x<x^{\prime})$ and $\rho_{0}(x<x^{\prime})$, can be derived similarly from characterizing the respective normalized spin function overlaps $()_{-,0}$ which we will evaluate below.

From now on we limit ourselves to the specific sector of total $S_{z}$ $\equiv$ $\sum_{i=1}^{N}s_{i}$ $=$ [math]. For $S_{z}$ close to $N$, spin-$1$ Bose gas will behave not much different from the polarized or spinless one. Therefore we choose the sector of zero $S_{z}$, which allows the SILL of spin-$1$ Bose gas to distinguish most significantly from the spinless bosons for $S_{z}$ $\lesssim$ $N$. The spin configurations $|\chi\rangle$ in this sector generally involve $n$ pairs of $(+-)$, that is $|+++---00...0\rangle$ with $n$ $=$ $3$ for example. The total number of states can then be calculated as

[TABLE]

which we obtain by permuting $n$ $(\pm)$’s and $(N-2n)$ $(0)$’s. For the spin-plus component of the single-particle density matrix in Eq. (9), the spin configuration $|00...0\rangle$ with $n$ $=$ [math] never contributes. Therefore we consider only the spin configurations of at least one pair of $(+-)$, and $|\chi\rangle$ can be generally expressed as

[TABLE]

The first $+$ is projected out in $\rho_{+}(x<x^{\prime})$, and thus we have the normalized spin function overlap $(E,E)_{+}$,

[TABLE]

which is averaged by $w_{N}$, the total number of states. We note that all the arguments of the factorials should be equal and larger than zero. $(E,E)_{+}$ is proportional to the number of states obtained by permuting the rest of $(n-1)$ $(+)$’s, $n$ $(-)$’s, and $(N-2n)$ $(0)$’s. For $(E,P_{12...j})_{+}$, it has a contribution only when the first $j$ entries are $(+)$’s,

[TABLE]

such that we have

[TABLE]

which denotes the number of states obtained by permuting the rest of $(n-j)$ $(+)$’s, $n$ $(-)$’s, and $(N-2n)$ $(0)$’s. In this specific sector of zero $S_{z}$, we note that in general $(E,P_{12...j})_{+}$ is nonvanishing only when $j$ $\leq$ $N/2$.

These corresponding spin function overlaps in $\rho_{-}(x<x^{\prime})$, which are $(E,E)_{-}$ and $(E,P_{12...j})_{-}$, should be the same as those in $\rho_{+}(x<x^{\prime})$. While for $\rho_{0}(x<x^{\prime})$, we have the spin function overlaps as

[TABLE]

which respectively denote the number of states contributed from the spin configurations with the first one and the first $j$ entries of $(0)$’s. We note of the identity that

[TABLE]

This also corresponds to the particle number conservation, that is $2N_{+}$ $+$ $N_{0}$ $=$ $N$, where $N_{\pm(0)}$ $\equiv$ $\int dx\rho_{\pm(0)}(x,x)$. Thus the number of particles is proportional to the spin function overlaps, $N_{\pm(0)}$ $=$ $N(E,E)_{\pm(0)}$. Furthermore we note that

[TABLE]

where $w_{jN}$ was defined in Ref. Jen2016_spin1 ,

[TABLE]

III.1 Spatial correlation in SILL of spin-1 Bose gas

The effect of the spin function overlaps in the SILL regime can be seen in Fig. 1 where we compare the spatial correlations of spin-$1$ [$\rho(x,x^{\prime})$] and spinless bosons [$\rho_{\rm spl}(x,x^{\prime})$] at some chosen $x$ in a harmonic trap. Spinless or spin-polarized bosons show a wider spatial distribution than the spin-$1$ case, indicating a sharper momentum distribution. Large $|x-x^{\prime}|$ in the spatial correlation corresponds to the small $p$ region. For spinless bosons, it has been shown in the bulk that $\rho_{\rm spl,b}(x,x^{\prime})$ $\propto$ $|x-x^{\prime}|^{-1/2}$, thus small $p$ behavior is proportional to $|p|^{-1/2}$ Lenard1964 ; Vaidya1979 ; Jimbo1980 ; Forrester2003 . This narrow momentum distribution resembles the one of a Bose-Einstein condensate but not quite since no condensation is allowed Dalfovo1999 due to large quantum fluctuations in 1D system. Therefore no off-diagonal long-range order can be present in the density matrix of 1D Bose gas. However a superfluid phase can exist in 1D quantum systems, possessing a power-law decay in spatial correlations. This power-law decay can be well described in the Luttinger liquid model using the bosonization method Giamarchi2004 ; Haldane1981 .

In a harmonic trap as shown in Fig. 1, the spatial correlations of $\rho_{\rm spl}(x,x^{\prime})$ are similar to the one in a bulk in a moderate region of $|x-x^{\prime}|$ until the correlation decays faster at the edge of the trap ($x^{\prime}$ $\gtrsim$ $4x_{\rm ho}$). In the trap, $\rho_{\rm spl}(p=0)$ is finite. It has also been shown that $\rho_{\rm spl}(p=0)$ $\propto$ $N$ in the large $N$ limit Papenbrock2003 .

In sharp contrast to the spinless bosons in Fig. 1, the spin-$1$ Bose gas in the SILL regime shows an exponential decay in its spatial correlation, which is therefore not condensed. This exponential decay has been predicted in the single-particle Green’s function of quantum wires in the SILL regime Cheianov2004 ; Fiete2004 ; Fiete2007 , distinguishing from the Luttinger liquid with only power-law decays. In the momentum distributions on the other hand, spin-incoherence tends to broaden the distributions, which has been investigated in the $t$-$J$ model Feiguin2010 ; Penc1996 ; Penc1997 or the system of uniform two-component gas Cheianov2005 . Similarly the spin-$1$ bosons in the SILL regime will also have a broadened momentum distribution due to the averaging of the spin function overlaps, which we discuss in more details below. Large $p$ behavior will be discussed later in Sec. IV. A.

III.2 Momentum distribution in SILL of spin-1 Bose gas

We define the spin-dependent momentum distributions as

[TABLE]

where we set $\hbar$ $=$ $1$. We then numerically calculate the momentum distributions of the three components in 1D TG Bose gas based on Eq. (9), $\rho_{-}(x<x^{\prime})$, and $\rho_{0}(x<x^{\prime})$. In Fig. 2, both spin components of $\rho_{+}(p)$ and $\rho_{0}(p)$ are uniformly broadened as $N$ grows, and $\rho_{+}(p)$ $\neq$ $\rho_{0}(p)$ for finite $N$. The effect of spin-incoherence also averages out the oscillatory structure that is present in the momentum distribution for specific spin state of spinor Bose gas Deuretzbacher2008 . Furthermore the peaks of $\rho_{0}(p)$ are larger than $\rho_{+}(p)$ up to $N$ $=$ $16$. This is due to $N_{0}$ $\geq$ $N_{\pm}$ in general and the spin function overlaps $(E,P_{12...j})_{0}$ are always larger than $(E,P_{12...j})_{+}$, which we will show more specifically in Fig. 9 in Appendix. For spinless bosons, the peaks of $\rho_{\rm spl}(p)$ have a scaling of $\rho_{\rm spl}(p=0)$ $\propto$ $N$ Papenbrock2003 . Here the spin-$1$ Bose gas in the SILL regime shows fitted scalings of $\rho_{+}(p=0)$ $\propto$ $N^{0.49}$ and $\rho_{0}(p=0)$ $\propto$ $N^{0.66}$ from Fig. 2. These reduced scalings again show the feature of broadened momentum distributions in the SILL regime

To calculate $\rho_{\pm(0)}(x,x^{\prime})$ we implement Gaussian unitary ensemble (GUE) Papenbrock2003 to speed up the convergence in the Monte Carlo (MC) integration method. The GUE draws a series of $(N-1)$ random numbers in $\bar{x}$, which are repulsively distributed due to the joint probability density of $\Pi_{1\leq i<j\leq N-1}(x_{i}-x_{j})^{2}$. This implementation of GUE thus enables our MC integration to simulate up to $N$ $=$ $16$, which in this case takes about $140$ hours with MC simulations of $M$ $=$ $10^{6}$ sets of random numbers using $200$ parallel CPU cores. All MC simulations in Fig. 2 use $M$ $=$ $10^{7}$ except for $N$ $=$ $16$ with $M$ $=$ $10^{6}$.

In the next section we investigate their asymptotic forms in large momentum limit, which show $1/p^{4}$ decay, and their momentum distributions in large $N$ limit.

IV Momentum distributions in high $\boldsymbol{p}$ and large $\boldsymbol{N}$ limits

IV.1 Asymptotic high $\boldsymbol{p}$ limit

For spinless bosons in the TG limit, relative wave function between two particles in short distance is $\psi_{\rm rel}(x,x^{\prime})$ $\propto$ $|x-x^{\prime}|$, indicating of impenetrable bosons and corresponding to the feature of fermionic repulsion. Again it has been shown Lenard1964 ; Vaidya1979 ; Jimbo1980 ; Forrester2003 in a bulk where $\rho_{\rm spl,b}(x,x^{\prime})$ in short distance is proportional to $[1+...+|x-x^{\prime}|^{3}/(9\pi)+...]$. Thus the non-analytic $|x-x^{\prime}|^{3}$ term in the short-distance correlation gives a universal $1/p^{4}$ asymptotic in large momentum limit. This universal $1/p^{4}$ asymptotic is not unique for a Bose gas with two-body contact interactions Minguzzi2002 ; Olshanii2003 ; Xu2015 . It also shows up in Tan’s relation Tan2008 ; Barth2011 in the two-component Fermi gas Braaten2008-1 ; Braaten2008-2 ; Werner2009 ; Zhang2009 .

For 1D spin-$1$ TG Bose gas on the other hand, the analytical results for a high $p$ asymptotic total momentum distribution $\rho(p)$ have been derived Jen2016_spin1 , showing also a universal $1/p^{4}$ dependence. Similarly for its spin-dependent components, they can be straightforwardly written as

[TABLE]

where $(n_{i},n_{j})$ denotes any possible pairs of $N$ harmonic oscillator eigenfunctions. The asymptotic form depends on the spin function overlap $(E,P_{12})_{\pm(0)}$ because it has significant contributions only from the integral regions of $x$ $<$ $x_{j}$ $<$ $x^{\prime}$ and $x^{\prime}$ $<$ $x_{j}$ $<$ $x$ for all $x_{j}$ $\in$ $\bar{x}$ with $x$ $\approx$ $x^{\prime}$. The asymptotic form for the spinless bosons can be also obtained by replacing $[(E,E)_{\pm(0)}+(E,P_{12})_{\pm(0)}]$ with $2$ in Eq. (21). We note that using Eqs. (21) and (17), we have $2\rho_{+}(p)$ $+$ $\rho_{0}(p)$ $=$ $\rho(p)$ where the last quantity was computed in Ref. Jen2016_spin1 .

The spin-$1$ Bose gas in the SILL regime shows very different properties from the spin-coherent ones in the coefficients of high $p$ asymptotics. The coefficients are always less than the ones in spinless bosons since $(E,E)_{\pm(0)}$, $(E,P_{12})_{\pm(0)}$ $<$ $1$. And for large $N$, $[(E,E)_{\pm(0)}+(E,P_{12})_{\pm(0)}]$ $\rightarrow$ $[1/3+(1/3)^{2}]$ $=$ $4/9$ from Eq. (59), less than $2$ for the case of spinless bosons as well. As an example, in Fig. 3 we compare the numerical and analytical results of $\rho_{+}(p)$ in high $p$ limit. The numerically calculated high $p$ asymptotics approach approximately to the analytical ones. For even larger $px_{ho}$ $\gtrsim$ $7$, the trends either drop and cross the analytical asymptotics, or bounce back and oscillate, indicating the inaccuracy of numerical results in these regions. To reach accurate high $p$ asymptotics is quite demanding in MC integrations and consuming more CPU time for even larger $N$. However, MC simulations have already achieve the accuracy of $10^{-3}$ and $10^{-2}$ of the momentum distributions for $N$ $=$ $2-3$ and $10$ respectively.

We have also evaluated numerically the potential $(\langle V\rangle)$ and kinetic energies $(\langle K\rangle)$. Since our 1D bosonic TG gas has the same density distribution as the one of a Fermi gas, we have $\langle V\rangle$ $=$ $\langle K\rangle$ $=$ $N^{2}\hbar\omega/4$, equivalent to half of the total energy, which complies with the Virial theorem Werner2006 ; Werner2008 . In Ref. Jen2016_spin1 , we concatenate $\rho(p)$ with the asymptotic tails analytically derived to improve the energy calculations. Here we directly use the momentum distributions of $\rho(p)$ calculated by MC simulations implemented with GUE. We find that the numerical results of these energies improve to the relative errors below $7\%$ and $10\%$ for $N$ $=$ $2-7$ and $16$ respectively to the exact values of $\langle V\rangle$ and $\langle K\rangle$. This further shows the advantage of GUE in the convergence and accuracy of our numerical results.

In Fig. 4, we plot the difference of spin-plus and zero momentum distributions numerically, $\rho_{+}(p)$ $-$ $\rho_{0}(p)$, from Fig. 2. The difference goes away gradually as $N$ increases, indicating these two components approach each other in large $N$ limit. The dips at around $p$ $\sim$ [math] demonstrate that the peaks of $\rho_{0}(p)$ are always higher than $\rho_{+}(p)$, which is due to larger spin function overlaps for the spin-[math] component. A special feature of the peaks for the case of $N$ $=$ $2$ shows a wider $\rho_{+}(p)$ than $\rho_{0}(p)$ while this feature is not obvious for larger $N$.

IV.2 Large $\boldsymbol{N}$ limit

Due to the limits of numerical integration, we can only calculate the single-particle density matrix of spin-$1$ Bosons up to $N$ $=$ $16$. For finite $N$, we have demonstrated numerically $\rho_{+}(p)$ $\neq$ $\rho_{0}(p)$ since in general $N_{0}$ $\geq$ $N_{+}$ and spin function overlaps $(E,P_{12...j})_{0}$ are always larger than $(E,P_{12...j})_{+}$. Thus the peaks of $\rho_{0}(0=0)$ are larger than $\rho_{+}(p=0)$. In this subsection we attempt to investigate the momentum distributions in large $N$ limit. The study in this limit can give insight to practical experiments where several hundreds or thousands of atoms are involved.

To investigate the individual components of spin-$1$ momentum distributions in large $N$ limit, we need the asymptotic forms for various spin function overlaps. These spin function overlaps in general can be written as

[TABLE]

where

[TABLE]

$f^{(N)}_{k}$ are just the coefficients of $x^{k}$ in the binomial expansions of $(x+x^{-1}+1)^{N}$. We then further express them in terms of the complex integral as shown in Appendix, and find the asymptotic form of $f^{(N)}_{k}$ in large $N$ limit using the method of steepest descent or stationary phase Bender1999 . In Fig. 9 of Appendix, the asymptotic forms in Eq. (57) are used to compare with the exact ones of Eqs. (11), (12), (13), and (14), and they approach to the exact ones in large $N$ limit for small $j$ in Eq. (22). Therefore we shall use the asymptotic forms to compare the spin-plus with the spin-[math] components of the spin function overlaps for large $N$.

We define their relative deviations as

[TABLE]

which asymptotically approaches to

[TABLE]

in large $N$ limit, where $\bar{f}^{(N-j)}_{j}$ is the asymptotic form of $f^{(N-j)}_{j}$. This asymptotic form allows us to compute the deviations for even larger $N$ than using the exact formulas. In Fig. 5, the relative deviations decay as $N$ increases for small $j$ and become below $10^{-2}$ for $N$ $\sim$ $600$ with $j$ $\leq$ $3$. Note that for a moderate $j$ $=$ $10$, it only reaches $0.1$ for as high as $N$ $=$ $600$. Since the spin function overlaps with smaller $j$ contribute much more to the momentum distributions than $j$ $\lesssim$ $N/2$ [see Figs. 9(a) and (b)], along with much less relative deviations, we expect that the spin-plus momentum distribution approaches the one of spin-[math] as $N$ increases. We have also shown this trend in Fig. 4 for finite $N$ up to sixteen particles.

This can be further confirmed by studying respective contributions from the integral regions in the single-particle density matrix $\rho_{+}(x<x^{\prime})$ of Eq. (9). In Fig. 6 we plot the results of most significant $12$ out of a total $15$ integral regions for the case of $N$ $=$ $5$ with and without the multiplications of spin function overlaps. The order of the integral regions can be seen from Eq. (9). We denote the first five regions as $x$ $<$ $x_{2}$ with $x^{\prime}$ $<$ $x_{2}$ moving sequentially to $x_{5}$ $<$ $x^{\prime}$, the next four regions ($6$th to $10$th) as $x_{2}$ $<$ $x$ $<$ $x_{3}$ with $x_{2}$ $<$ $x^{\prime}$ $<$ $x_{3}$ moving sequentially to $x_{5}$ $<$ $x^{\prime}$, and so on for the rest of the regions. A feature of hard core bosons in 1D harmonic trap is that atoms prefer to distribute evenly in space. The particles at $\bar{x}$ repel each other due to their strongly repulsive interactions in TG limit. Thus more significant peaks of the integral contributions occur roughly at the corresponding integral regions depending on $(x,x^{\prime})$. For example in Fig. 6(a), the spatial correlation at $(x,x^{\prime})$ $=$ $(-1,0)$ has the largest contribution at the seventh integral region, $x_{2}$ $<$ $x$ $<$ $x_{3}$ $<$ $x^{\prime}$ $<$ $x_{4}$ $<$ $x_{5}$. Similarly the spatial correlation at $(x,x^{\prime})$ $=$ $(-1,3)$ has the largest contribution at the ninth integral region, $x_{2}$ $<$ $x$ $<$ $x_{3}$ $<$ $x_{4}$ $<$ $x_{5}$ $<$ $x^{\prime}$. In Fig. 6(b) we multiply these values with the spin function overlaps (inset), and show that significant contributions come from the overlaps with small $j$. For four specific spatial correlations we choose here, significant values dwell in the $6$th, $7$th, $10$th, and $11$th integral regions, which correspond to the spin function overlaps $(E,E)_{+}$, $(E,P_{12})_{+}$, $(E,E)_{+}$, and $(E,P_{12})_{+}$ respectively. Moreover as expected the spatial correlations decay as $|x-x^{\prime}|$ increases, reminiscent of the exponential decay discussed in III. A.

For TG bosons in the bulk with a length $L$ in thermodynamic limit, we can use the analytical expression of Eq. (21) and let $\phi_{n_{i}}(x)$ $=$ $e^{ik_{i}x}/\sqrt{L}$ with various eigenmodes $n_{i}$ $\rightarrow$ $k_{i}$. The spatial integral in Eq. (21) can then be calculated as

[TABLE]

In the continuous limit of $(k_{i},k_{j})$ corresponding to the pairs of $(n_{i},n_{j})$ in the summation, we have $\rho_{\rm spl,b}(p\rightarrow\infty)$ in the bulk

[TABLE]

where $k_{F}$ $=$ $\pi N/L$. Let $k_{i}$ $=$ $k_{F}x$ and $k_{j}$ $=$ $k_{F}y$, the above becomes

[TABLE]

which is the same as Eq. (67) in Ref. Forrester2003 up to a factor of $2\pi$ in our definition of Fourier transform in Eq. (18). Our derivation is parallel to using short-distance expansions in $\rho_{\rm spl,b}(x,x^{\prime})$ Forrester2003 where its non-analytic term $|x-x^{\prime}|^{3}/(9\pi)$ after Fourier transformed gives the same result.

For large $N$ spinless bosons in a harmonic trap, we can use the local-density approximation (LDA) for the local chemical potential in Thomas-Fermi limit, which reads $\mu(x)$ $=$ $\mu$ $-$ $\alpha x^{2}/2$ with $\alpha$ $=$ $M\omega^{2}$. The cutoff momentum can be derived as $k_{F}(x)$ $=$ $\sqrt{2M(\mu-\alpha x^{2}/2)}$. Since the maximum mode in 1D hard core bosons at infinite interactions is approximately $n_{\rm max}$ $\approx$ $N$, the chemical potential can be determined as $\mu$ $=$ $n_{\rm max}\omega$ $\approx$ $N\omega$. We further define the boundary of $x_{\rm max}$ $=$ $\sqrt{2N\omega/\alpha}$ such that we re-express $k_{F}(x)$ as $\sqrt{(M\alpha)(x_{\rm max}^{2}-x^{2})}$. Applying LDA to Eq. (34), we have the density matrix for spinless bosons in a harmonic trap as

[TABLE]

where the scaling of $N^{5/2}$ disagrees with Ref. Olshanii2003 which gave a different scaling of $N^{3/2}$. $N^{5/2}$ scaling has also been reported for the 1D SU($\kappa$) Fermi gas with $\kappa$ $\neq$ $1$ Decamp2016 . The coefficient in Eq. (35) is the same as the TG Fermi gas in the $\kappa$ $\rightarrow$ $\infty$ limit. We note that the coefficient of the scaling depends on the many-body state and is related to the slope of energy ($-dE/dg_{\rm 1D}^{-1}$) Olshanii2003 ; Decamp2016 .

For the spinful case of our spin-$1$ Bose gas in large $N$ limit, the asymptotic coefficient of $[(E,E)_{\pm(0)}+(E,P_{12})_{\pm(0)}]$ becomes $[1/3+(1/3)^{2}]$ $=$ $4/9$ according to Eq. (59), such that the spin-dependent and total momentum distributions respectively in thermodynamic limit become

[TABLE]

and

[TABLE]

which is $2/3$ of Eq. (35).

We denote the momentum distribution of spinless bosons in high $p$ limit as $\rho_{\rm spl}(p\rightarrow\infty)$ which can be derived by replacing $[(E,E)_{\pm(0)}+(E,P_{12})_{\pm(0)}]$ with $2$ in Eq. (21). We then define $c(N)$ $\equiv$ $\rho_{\rm spl}(\infty)(2^{6}\sqrt{2})^{-1}N^{-5/2}(45\pi^{3}p^{4}x_{ho}^{3})$ and $c(\infty)$ $=$ $2$ according to Eq. (35). $c_{\pm(0)}(N)$ can be also defined in the same way for $\rho_{\pm(0)}(\infty)$. These coefficients can be calculated using Eq. (21). In Fig. 7 we plot $c(N)$ and $c_{+(0)}(N)$ respectively to show how they approach the asymptotic large $N$ values. We find that Eqs. (35) and (36) already give good enough estimates for $N$ $\gtrsim$ $20$ and $30$ for spinless and spin-$1$ bosons respectively. The relative deviations $|c(N)-c(\infty)|/c(\infty)$ $=$ $0.14\%$ for $N$ $=$ $40$, and for the cases of $c_{+(0)}(40)$, they reach $1.9\%(1.5\%)$. In Fig. 7(b), $c_{0}$ $>$ $c_{+}$, which again indicates that the spin function overlaps $(E,P_{12...j})_{0}$ $>$ $(E,P_{12...j})_{+}$ and $N_{0}$ $>$ $N_{+}$.

V Conclusion

In conclusion, we have investigated the spin-dependent properties of spin-$1$ Bose gas in the regime of spin-incoherent Luttinger liquid (SILL). Three components (spin-plus, zero, and minus) of the single-particle density matrix for this universal class can be calculated by deriving respective spin function overlaps. These spin function overlaps result from the highly degenerate spin configurations in the SILL regime. In contrast with the spinless bosons with a power-law decay in its spatial correlation, spin-$1$ bosons in TG limit show an exponentially decaying spatial correlation, which indicates a broadened momentum distribution and a different universal class from Luttinger liquid. The universal $1/p^{4}$ dependence in high $p$ limit is also present in the spin-dependent momentum distributions. This asymptotic has a scaling of $N^{5/2}$ with a reduced coefficient than the one of the spinless bosons. The coefficients of the asymptotic are proportional to Tan’s contact and can be observed in experiments as one of the signatures of SILL. We compare these analytical predictions with the numerical results calculated by Monte Carlo (MC) integration with Gaussian unitary ensemble (GUE) up to sixteen bosons. The method of MC integration implemented with GUE converges faster and gives more accurate results such that we are able to calculate higher $p$ regions. The high momentum tails approximately and asymptotically follow the reduced coefficients we analytically derived.

For the $S_{z}$ $=$ [math] sector, we show that the spin-[math] component always has a larger peak than the spin-plus momentum distribution for finite $N$. This can be explained by the spin function overlaps which are larger for the spin-[math] density matrix than the spin-plus case. While they differ for small $N$, they coincide in the large $N$ limit. This indicates that highly incoherent bosons form in this limit with each component occupying exactly one third of the total number of particles. The ultracold spinor Bose gas allows for a potential realization of this universal class of SILL, and our results offer a testable paradigm to study quantum many-body phenomena in 1D strongly interacting bosons.

Acknowledgements

This work is supported by the Ministry of Science and Technology, Taiwan, under Grant No. 104-2112-M-001-006-MY3. The work of SKY was partially supported by a grant from the Simons Foundation, and was performed at the Aspen Center for Physics supported by National Science Foundation Grant No. PHY-1066293.

Appendix A Spin function overlaps in large $\boldsymbol{N}$ limit

Before we derive the spin function overlaps in large $N$ limit, first we express them in terms of an integral function, and then introduce the method of stationary phase or steepest descent to solve for their asymptotic forms Bender1999 . Consider the following function,

[TABLE]

where the second line above can be derived by binomial expansions. When we set $k$ $=$ $j_{1}$ $-$ $j_{2}$ and $j$ $=$ $j_{2}$, we have

[TABLE]

where the upper bound of index $j$ can be derived by solving $j_{2}$ in the equations of $j_{1}$ $+$ $j_{2}$ $=$ $N$ and $j_{1}$ $-$ $j_{2}$ $=$ $k$. From Eq. (39), we define the coefficient of $x^{k}$ in $f^{(N)}(x)$ as $f^{(N)}_{k}$ which is the same as the one of $x^{-k}$. Compare with Eqs. (11), (12), (13), and (14) in the main paper, we find that the spin function overlaps can be expressed as Eq. (22) in the main text. These coefficients can be calculated as follows

[TABLE]

where the coefficient $f^{(N)}_{k}$ (of $z^{-k}$ in this case) is exactly the residue at the pole of $z$ $=$ [math] with a complex number $z$. $C$ in the above denotes a contour integration on a circle in a counterclockwise direction around the origin.

The asymptotic behavior of the integral can be obtained in large $N$ limit. We let $z$ $=$ $e^{i\theta}$, the coefficient becomes

[TABLE]

where $\mathcal{F}(\theta)$ $\equiv$ $N\ln(1+2\cos\theta)+ik\theta$. Let $\theta$ $\rightarrow$ $w$ in the complex plane, using the method of steepest descent Bender1999 for the above highly-oscillating integrals in large $N$ limit, we first find the saddle points which satisfy the first derivative $\mathcal{F}^{\prime}(w)$ $=$ [math]. The integrals are then dominated by the local maxima passing the saddle points along the integration contour.

The saddle points are therefore located at

[TABLE]

which, after replacing trigonometry functions with exponentials, becomes

[TABLE]

If $k$ $\rightarrow$ [math], the above suggests a multiple of roots that satisfy $e^{iw}$ $=$ $\pm 1$, which are $w$ $=$ $2n\pi$ and $\pm(2n+1)\pi$ for integers $n$, indicating multiple saddle points in this integral. We consider an integration path in Fig. 8, where only three saddle points ($n$ $=$ [math]) are involved. Below we demonstrate why the contour is valid and guarantees to follow the valleys of steepest descent between these three saddle points, which can be determined by the sign of $\mathcal{F}^{\prime\prime}(w)$.

First to calculate $\mathcal{F}^{\prime\prime}(w_{0})$, we define $Q$ $\equiv$ $\sqrt{4-3(k/N)^{2}}$, and we have from Eq. (43) with the ”$+$” sign,

[TABLE]

where $w_{0}$ should be purely imaginary in general. We further use the above to reinterpret

[TABLE]

Now the second derivative of $\mathcal{F}(w_{0})$ becomes

[TABLE]

which is always less than zero. For example, $\mathcal{F}^{\prime\prime}(w_{0})$ $=$ $-2N/3$ and $-N\epsilon$ respectively at small and large $k$ limit ($k/N$ $=$ $1$ $-$ $\epsilon$ with $\epsilon$ $\gtrsim$ [math]).

Next for $\mathcal{F}^{\prime\prime}(w_{\pm 1})$ at the other two saddle points which we denote as $w_{\pm 1}$, we have from Eq. (43) with the ”$-$” sign,

[TABLE]

where $w_{\pm 1}$ in general are complex. Setting $w_{\pm 1}$ $=$ $\pm\pi$ $+$ $iy_{1}$ with real $y_{1}$, we have $e^{-y_{1}}$ $\lesssim$ $1$ and $e^{-y_{1}}$ $\approx$ $1/2$ respectively when $k$ $\rightarrow$ [math] and $k$ $\rightarrow$ $N$, suggesting $y_{1}$ $\gtrsim$ [math] at small $k$ limit. Again we can use the above to reinterpret

[TABLE]

Now the second derivative of $\mathcal{F}(w_{\pm 1})$ becomes

[TABLE]

which is again always less than zero, for example, $\mathcal{F}^{\prime\prime}(w_{0})$ $=$ $-2N$ and $-2N/3$ respectively for small and large $k$ limit.

Now we have located three saddle points, which are $w_{0}$ and $w_{\pm 1}$. The asymptotic form of the integral for $f^{(N)}_{k}$ in Eq. (41) can then be calculated using the contour in Fig. 8, which traverses through these three saddle points on the paths tangent to the real axis. Following the integration contour for Eq. (41), we have the asymptotic form $\bar{f}_{k}^{(N)}$ of $f^{(N)}_{k}$ as

[TABLE]

We obtain the contributions of $w_{0}$ in the path $C_{0}$ as

[TABLE]

where its next correction term is at least $\mathcal{O}(N^{-1/2})$ smaller. Since the integration path $C_{0(\pm 1)}$ follows the valley of steepest descent, the integration is dominated near the region of the saddle point $w_{0(\pm 1)}$, where we are also able to allow the boundary of $w$ to $\pm\infty$. Similarly for the paths $C_{\pm 1}$, we have the contributions from $w_{\pm 1}$,

[TABLE]

where $\mathcal{F}(w_{-1})$ $=$ $\mathcal{F}(w_{1})$. Let $w$ $=$ $w_{\pm 1}$ $+$ $x$ in the paths $C_{\pm 1}$ respectively, we have

[TABLE]

where we have used $\mathcal{F}^{\prime\prime}(w_{-1})$ $=$ $\mathcal{F}^{\prime\prime}(w_{1})$. For the paths of $D_{\pm 1}$, though they do not follow the valleys of steepest descent, their contributions cancel with each other since $\mathcal{F}(w)$ $=$ $\mathcal{F}(w+2\pi)$.

Finally we obtain the asymptotic form of $f^{(N)}_{k}$ in large $N$ limit as

[TABLE]

Before we write down the explicit form for the above, it is useful to derive

[TABLE]

which become $3^{N}$ and $1/\epsilon^{N\epsilon}$ respectively for small and large $k$. Also we have

[TABLE]

which become $(-1)^{N+k}$ and $(3/4)^{N}2^{N\epsilon}/(-1)^{N\epsilon}$ respectively for small and large $k$. Inserting the functions of the saddle points from Eqs. (46), (49), (55), and (56), we have

[TABLE]

We note that the main contribution in Eq. (57) comes from the saddle point $w_{0}$. To have some estimates of Eq. (57), we have

[TABLE]

From the above result at small $k$ and according to Eq. (22), we can show that

[TABLE]

which respectively indicates one third decrease for the spin function overlaps when $j$ increases by one, and the populations of $N_{+}$ and $N_{0}$ coincide in large $N$ limit. In the large $N$ limit, despite the constraint $S_{z}$ $=$ [math], the probability of finding a particle in any of the spin states, ”$+$”, ”[math]”, or ”$-$”, is $1/3$ and is irrespective of the spins of the other particles (if $j$ $\ll$ $N$). This decrease in spin function overlaps also reflects on the exponential decay in spatial correlations discussed in Sec. III. A.

To have some estimates of the spin function overlaps and their asymptotic forms in large $N$ limit, in Fig. 9 we compare $(E,P_{12...j})_{+,0}$ of Eqs. (12) and (14) with $\bar{f}^{(N-j)}_{j,0}/\bar{f}^{(N)}_{0}$ from Eq. (57). In Figs. 9(a) and (b), the values of spin function overlaps decrease rather fast in logarithmic scales as $j$ increases, while $(E,P_{12...j})_{0}$ is always larger than $(E,P_{12...j})_{+}$. As a comparison, we define the relative deviations as

[TABLE]

which we show in Figs. 9(c) and (d), indicating of a good asymptotic form from our derivations for small $j$ in large $N$ limit. However a slow decay in the relative deviation of $(E,P_{12...N/2-1})_{+}$ to $\bar{f}^{(N/2+1)}_{N/2-1}/\bar{f}^{(N)}_{0}$ in (c) shows the worst case in the asymptotic form. It is due to a rather small $\mathcal{F}^{\prime\prime}(w_{0})$ $\rightarrow$ $-N\epsilon$ in Eq. (46) when we set $k/N$ $=$ $1$ $-$ $\epsilon$ with a small value of $\epsilon$, which makes the method of steepest descent less accurate unless we go to $N$ $\rightarrow$ $\infty$.