Clusterization of weakly-interacting bosons in one dimension: an analytic study at zero temperature

TL;DR

This paper analytically investigates the clustering transition of weakly-interacting bosons in one dimension at zero temperature, comparing GP and variational approaches, and analyzing superfluid properties.

Contribution

It provides an analytic comparison of GP and variational methods for boson clusterization and characterizes the phase transition and superfluid fraction.

Findings

Clusterization occurs at a specific density with phase transition predicted by all models.

GP theory links clusterization with roton softening.

Analytic expressions for high-density behavior and superfluid fraction are derived.

Abstract

We study a system of penetrable bosons on a line, focusing on the high-density/weak-interaction regime, where the ground state is, to a good approximation, a condensate. Under compression, the system clusterizes at zero temperature, i.e., particles gather together in separate, equally populated bunches. We compare predictions from the Gross-Pitaevskii (GP) equation with those of two distinct variational approximations of the single-particle state, written as either a sum of Gaussians or the square root of it. Not only the wave functions in the three theories are similar, but also the phase-transition density is the same for all. In particular, clusterization occurs together with the softening of roton excitations in GP theory. Compared to the latter theory, Gaussian variational theory has the advantage that the mean-field energy functional is written in (almost) closed form, which enables us to extract the phase-transition and high-density behaviors in fully analytic terms. We also compute the superfluid fraction of the clustered system, uncovering its exact behavior close, as well as very far away from, the transition.

Equations

$$\Psi (x_1,\dots ,x_N)=\prod _{i=1}^N\psi (x_i).$$

$$-\frac{{\hslash }^2}{2m}{\psi }^{\prime \prime }(x)+N{\int }_V\mathrm{d}y|\psi (y){|}^{2}u(x-y)\psi (x)=\mu \psi (x),$$

$${\int }_V\mathrm{d}x|\psi (x){|}^2=1.$$

$$\psi (x)=\frac{1}{\sqrt{V}}\sum _G{c}_G{e}^{iGx},$$

$$\left(\frac{{\hslash }^2{K}^2}{2m}+\rho \tilde{u}(0)\right)c_K+\rho \sum _{G\ne 0}\tilde{u}(G)S_G{c}_{K+G}=\mu c_K,$$

$$\mathcal{E}$$

$$\psi (x)=\frac{1}{\sqrt{V}}\eta (x)e^{i\theta (x)}$$

$$\mathcal{E}=\frac{{\hslash }^2}{8mV}{\int }_V\mathrm{d}x\left(\frac{{\eta }^{\prime 2}(x)}{\eta (x)}+4\eta (x){\theta }^{\prime 2}(x)\right)+\frac{\rho }{2V}{\int }_V\mathrm{d}x{\int }_V\mathrm{d}{x}^{\prime }\eta (x^{\prime })u(x-x^{\prime })\eta (x),$$

$$\mathrm{VT}1:\psi (x)=C_{\alpha }\frac{1}{\sqrt{V}}\sum _{n=-\infty }^{+\infty }{e}^{-\alpha {\left(\frac{x}{a}-n\right)}^2},$$

$$C_{\alpha }={\left(\frac{2\alpha }{\pi I(\alpha {)}^2}\right)}^{1/4},\mathrm{with}I(\alpha )=\sum _{n=-\infty }^{+\infty }{e}^{-\frac{\alpha }{2}{n}^2}.$$

$$I(\alpha )={\vartheta }_3(0,e^{-\alpha /2}).$$

$$\psi (x)=\frac{1}{\sqrt{V}}\sum _G{\psi }_G{e}^{iGx},$$

$${\psi }_G=C_{\alpha }^{\prime }{e}^{-\frac{G^2{a}^2}{4\alpha }}\mathrm{and}{C}_{\alpha }^{\prime }={\left(\frac{2\pi }{\alpha I(\alpha {)}^2}\right)}^{1/4}.$$

$$I(\alpha )={\left(\frac{2\pi }{\alpha }\right)}^{1/2}\sum _G{e}^{-\frac{G^2{a}^2}{2\alpha }}.$$

$${\mathcal{E}}_{\mathrm{kin}}=e_0\frac{\alpha {\sigma }^2}{2a^2}\left(1+2\alpha \frac{I^{\prime }(\alpha )}{I(\alpha )}\right).$$

$${\mathcal{E}}_{\mathrm{pot}}$$

$$J(\alpha )=\sqrt{\frac{2\pi }{\alpha }}\sum _{n=-\infty }^{+\infty }{e}^{-\frac{2{\pi }^2}{\alpha }{\left(n+\frac{1}{2}\right)}^2}.$$

$$\mathrm{VT}2:\eta (x)=\sqrt{\frac{2\alpha }{\pi }}\sum _{n=-\infty }^{+\infty }{e}^{-2\alpha {\left(\frac{x}{a}-n\right)}^2}$$

$$\eta (x)=\sum _G{\eta }_G{e}^{iGx}.$$

$${\eta }_G=\frac{1}{a}{\int }_{\mathcal{C}}\mathrm{d}x{e}^{-iGx}\eta (x)=\frac{1}{V}\sqrt{\frac{2\alpha }{\pi }}\sum _R\underset{1}{\underbrace{e^{-iGR}}}{\int }_V\mathrm{d}x{e}^{-iG(x-R)}{e}^{-2\alpha {\left(\frac{x-R}{a}\right)}^2}=e^{-\frac{G^2{a}^2}{8\alpha }}.$$

$$\eta (x)=\sum _G{e}^{-\frac{G^2{a}^2}{8\alpha }}{e}^{iGx}=\sum _G{e}^{-\frac{G^2{a}^2}{8\alpha }}\cos (Gx).$$

$${\mathcal{E}}_{\mathrm{pot}}$$

$$\tilde{u}(k)=2ϵ\frac{\sin (k\sigma )}{k}\mathrm{and}{\mathcal{E}}_{\mathrm{pot}}=\rho \sigma ϵ+\frac{\rho a}{\pi }ϵ\sum _{n=1}^{\infty }\frac{1}{n}{e}^{-\frac{{\pi }^2}{\alpha }{n}^2}\sin \left(\frac{2\pi }{a}n\sigma \right).$$

$${\mathcal{E}}_{\mathrm{kin}}=\frac{{\hslash }^2}{8mV}{\int }_V\mathrm{d}x\frac{{\eta }^{\prime 2}(x)}{\eta (x)}=\frac{{\hslash }^2}{8ma}{\int }_{\mathcal{C}}\mathrm{d}x\frac{{\eta }^{\prime 2}(x)}{\eta (x)}.$$

$$\eta (x)={\vartheta }_3\left(\frac{\pi }{a}x,e^{-\frac{{\pi }^2}{2\alpha }}\right).$$

$$\frac{{\eta }^{\prime 2}(x)}{\eta (x)}={\eta }^{\prime }(x)\frac{{\eta }^{\prime }(x)}{\eta (x)}=-16\frac{{\pi }^2}{a^2}\sum _{n_1,n_2=1}^{\infty }(-1{)}^{n_1}{n}_2\frac{e^{-\frac{{\pi }^2}{2\alpha }(n_1+n_2^2)}}{1-e^{-\frac{{\pi }^2}{\alpha }{n}_1}}\sin \left(\frac{2\pi }{a}{n}_{1}x\right)\sin \left(\frac{2\pi }{a}{n}_{2}x\right).$$

$${\int }_{-a/2}^{a/2}\mathrm{d}x\sin \left(\frac{2\pi }{a}{n}_{1}x\right)\sin \left(\frac{2\pi }{a}{n}_{2}x\right)=\frac{a}{2}{\delta }_{n_1,n_2},$$

$${\mathcal{E}}_{\mathrm{kin}}=\frac{{\pi }^2{\sigma }^2}{a^2}{e}_0\sum _{n=1}^{\infty }(-1{)}^{n-1}\frac{n{e}^{-\frac{{\pi }^2}{2\alpha }(n+n^2)}}{1-e^{-\frac{{\pi }^2}{\alpha }n}}.$$

$$I(\alpha )=\sqrt{\frac{2\pi }{\alpha }}\sum _{n=-\infty }^{+\infty }{e}^{-\frac{2{\pi }^2}{\alpha }{n}^2}=\sqrt{\frac{2\pi }{\alpha }}{\vartheta }_3(0,e^{-\frac{2{\pi }^2}{\alpha }});$$

$$J(\alpha )=\sqrt{\frac{2\pi }{\alpha }}\sum _{n=-\infty }^{+\infty }{e}^{-\frac{2{\pi }^2}{\alpha }{\left(n+\frac{1}{2}\right)}^2}=\sqrt{\frac{2\pi }{\alpha }}{\vartheta }_2(0,e^{-\frac{2{\pi }^2}{\alpha }}).$$

$${\left(\frac{J(\alpha )}{I(\alpha )}\right)}^2$$

Full text

I Introduction

The quantum mechanics of one-dimensional (1D) many-body systems is, since many decades by now, a well-respected area of research at the boundary between physics and mathematics, with lots of exact results Cazalilla . However, until recent times 1D systems were deemed to be of little interest for real-world physics. Now, their study is experiencing a revival thanks to the availability of ultracold gases of atoms and dipolar molecules confined in highly anisotropic magnetic traps and optical lattices (owing to strong confinement in the transverse directions, only the lowest-energy transverse quantum state needs to be considered, and the three-dimensional problem becomes effectively one-dimensional). Furthermore, by tuning the transverse trap frequency it is possible to change the short-range part of the atom-atom interaction to a certain extent (see, e.g., Ref. Bergeman ), opening the way to a direct comparison between experiment and theory.

Early experimental and simulation works have focused on reproducing the physics of the exactly solvable Lieb-Liniger model Lieb1 ; Lieb2 , describing a 1D system of identical spinless bosons of mass $m$ interacting through a contact repulsion of strength $g$. For all values of the single parameter $\gamma=mg/(\hbar^{2}\rho)$ characterizing the model (with $\rho$ denoting the number density), the spectrum is phonon-like at small momenta, in accordance with Luttinger-liquid (LL) behavior Haldane . The system becomes more weakly-interacting with increasing $\rho$; indeed, it is when $\gamma$ is low that, according to a standard argument (see, e.g., Petrov ), the coherence (or healing) length $\hbar/\sqrt{mg\rho}$ ($m$ is the atom mass) is much larger that the mean interparticle separation $\rho^{-1}$. For small values of $\rho/g$ (large $\gamma$), the 1D fluid acquires fermionic properties as the ground-state wave function strongly decreases at short interparticle distances (Tonks-Girardeau limit) Girardeau . At zero temperature ($T=0$), the one-body density matrix exhibits a power-law decay at large distances for any $g$ (although with a $\rho$-dependent exponent becoming smaller and smaller as $\rho/g$ increases); therefore, strictly speaking there is no Bose-Einstein condensation. It was Hohenberg Hohenberg the first to rigorously prove that no form of long-range order (including off-diagonal long-range order) can exist in one- and two-dimensional quantum systems with continuous group symmetries for non-zero temperatures; later, Pitaevskii and Stringari Pitaevskii2 have extended the proof for 1D bosonic systems also to the zero-temperature case. However, if we add a trapping potential in the axial direction both the amplitude and phase fluctuations are suppressed at low $T$ and one has a true condensate in 1D, provided the number $N$ of particles is large enough Petrov .

Kinoshita and coworkers Kinoshita have prepared a gas of ultracold 87Rb atoms in a 1D, cigar-shaped region. Acting on the harmonic confinement in the transverse directions, they have been able to tune the coupling strength $\gamma$ in the axial direction, making atoms to resemble either a Bose-Einstein condensate ($\gamma\ll 1$) or a Tonks-Girardeau gas of impenetrable bosons ($\gamma\gg 1$). In a wide range of $\gamma$, the experimental data of Kinoshita et al. fit the exact solution for the ground state of the Lieb-Liniger model. Astrakharchik and Giorgini Astrakharchik have studied the same model by Monte Carlo simulation at $T=0$, again observing the crossover from the low-density/free-fermion/Tonks-Girardeau limit to the high-density/weak-interaction/Gross-Pitaevskii (GP) regime Gross1 ; Pitaevskii1 ; Gross2 .

In 1D systems without pinning or trapping potentials, the Tonks-Girardeau behavior is typically confined to low density; for dense systems other ground states may emerge, like in dipolar bosons, which form a crystal at very high density Arkhipov . A further example is a system of softly-repulsive bosons. In a classical context, it has been recognized that bounded potentials favor clustering at low temperature, i.e., crystals with multiply-occupied sites, even when the potential is purely repulsive Likos ; Mladek . Similarly, quantum cluster crystals have been predicted in dense 2D systems of penetrable bosons Pomeau ; Josserand ; Henkel ; Saccani ; Cinti ; Macri1 . In 1D, thermal fluctuations destabilize cluster-crystalline order in favor of cluster-dominated liquid phases with different average occupation Speranza ; Prestipino1 ; Prestipino2 . Quantum mechanics induces coherent delocalization even at $T=0$ Pitaevskii2 , transforming a cluster crystal into a cluster LL (CLL) Rossotti ; Teruzzi1 ; Teruzzi2 . On decompression, the latter phase undergoes a quantum phase transition into a LL without clusters. Rossotti et al. Rossotti have studied in detail the phase transition from LL to “dimer” CLL at $T=0$, which turns out to fall into the Ising universality class.

We here reconsider the $T=0$ phase diagram of weakly-interacting penetrable bosons in 1D, but now using the GP theory for the ground state, namely assuming a pure condensate from the outset. While this cannot be valid in general, the mean-field (MF) approximation is reasonable in the high-density region, which comprises the transition from the fluid to a high-occupancy CLL. As we did for the same system in two and three dimensions Prestipino3 , we describe the CLL phase using a variational wave function written as a sum of evenly spaced Gaussians Tarazona . Even though no truly long-range order can exist in an infinite 1D system at $T=0$ Pitaevskii2 , modeling the CLL phase as a crystal entails an error on the system energy and short-range spatial correlations that is small in the whole range of applicability of MF theory. Otherwise, we may think of a large system of particles arranged in a circle (i.e., a finite system with periodic conditions at its boundary); in this case, clusterization and crystallization become one and the same transition. We show that the optimal Gaussian condensate is an effective approximation to the exact GP ground state. Gaussian variational theory has a major advantage over GP theory, since it allows an analytic study of both the transition region and the ultra-high-density limit, further enabling us to compute the superfluid fraction of the cluster phase.

Before going any further, it is worth discussing the relevance of MF approximation for the effectively 1D systems that can be realized experimentally. For atoms in a cylindrical trap, with transverse level spacing $\hbar\omega_{\perp}$ greatly exceeding the MF energy $g\rho$, the radial extension of the wave function is $a_{\perp}=\sqrt{\hbar/(m\omega_{\perp})}$. As discussed e.g. in Pitaevskii3 , radial motion is frozen only when the product of the 1D density $\rho$ times the 3D scattering length $a_{\rm 3D}$ is very small ($\rho a_{\rm 3D}\ll 1$). On the other hand, the MF regime in 1D corresponds to $\rho|a_{\rm 1D}|\gg 1$ ($a_{\rm 1D}=-2\hbar^{2}/(mg)$ being the 1D scattering length), which is consistent with $\rho a_{\rm 3D}\ll 1$ only provided that $a_{\rm 3D}\ll a_{\perp}$ (since $a_{\rm 1D}=-a_{\perp}^{2}/a_{\rm 3D}$ for values of $a_{\rm 3D}$ far below the threshold $a_{\perp}/1.0326$ of confinement-induced resonance Olshanii ). The smaller the ratio $a_{\rm 3D}/a_{\perp}$ is, the wider the range of densities where the MF approximation holds. In principle, this range can be expanded by reducing the soft-core interaction strength — which, however, might be hard to achieve in practice.

The outline of the paper is the following. In Sec. II we introduce the model and the theory employed to study its phase behavior at $T=0$. In Sec. III we present our results. Concluding remarks are postponed to Sec. IV.

II Model and theory

We investigate a system of $N$ one-dimensional (1D) spinless bosons of mass $m$, interacting through a bounded potential $u(x)$, even function of its argument. The range $\sigma$ and strength $\epsilon$ of the potential set the units of length and energy, respectively. A paradigmatic example of bounded repulsion is the penetrable-sphere model (PSM) potential, $u(x)=\epsilon\Theta(\sigma-|x|)$, $\Theta$ being the Heaviside step function. In the MF approximation, the ground state of the system is represented as a pure condensate:

[TABLE]

The best choice of $\psi$ is that minimizing the expectation value of the Hamiltonian in the state $\Psi$, which corresponds to a single-particle wave function obeying the (time-independent) Gross-Pitaevskii (GP) equation (see, e.g., Ref. Rogel-Salazar ):

[TABLE]

where $\mu$ has to be adjusted so that $\psi$ is normalized:

[TABLE]

Hereafter, the 1D volume $V$ is considered as macroscopic, with $\rho=N/V$ finite (wherever appropriate, an integral over $V$ can be extended to the whole real axis).

In the aim to describe clusterization of the system at $T=0$, we make the ansatz

[TABLE]

where $G=(2\pi/a)n$ are reciprocal-lattice vectors and $\sum_{G}|c_{G}|^{2}=1$. The state function in Eq. (2.4) describes a 1D crystal of spacing $a$. As discussed in the Introduction, no long-range order can actually occur in 1D, even for $T=0$; pretending the opposite is true is clearly an approximation, which predicts the wrong decay of correlation functions at large distances but only barely affects the location of the clusterization transition in the $\epsilon\rightarrow 0$ limit.

Plugging Eq. (2.4) in the GP equation, we obtain Kunimi ; Prestipino3

[TABLE]

with $S_{G}=\sum_{G^{\prime}}c_{G^{\prime}}c_{G^{\prime}+G}^{*}$. The function $\widetilde{u}(k)$ is the real-valued Fourier transform of $u(x)$, satisfying $\widetilde{u}(k)=\widetilde{u}(-k)$. The fluid phase, corresponding to $c_{G}=\delta_{G,0}$, is a special solution to Eq. (2.5), with $\mu=\rho\widetilde{u}(0)$. For the $\psi(x)$ in Eq. (2.4), the MF energy per particle is given by Kunimi ; Prestipino3

[TABLE]

The fluid energy is $\rho\widetilde{u}(0)/2$. We see from Eq. (2.6) that, denoting $e_{0}$ the characteristic energy $\hbar^{2}/(m\sigma^{2})$, the system ground state is only controlled by the dimensionless quantity $\rho\sigma\epsilon/e_{0}$ (which in the following is referred to as the “density”) or, equivalently, by $\rho\widetilde{u}(0)/e_{0}$. A yet different expression of ${\cal E}$ is obtained in the Madelung representation, where the single-particle wave function is written as

[TABLE]

(the amplitude $\eta(x)$ and phase $\theta(x)$ of $\psi$ are real and periodic). One readily obtains Pomeau :

[TABLE]

which makes it clear that the ground-state wave function is necessarily real.

It is useful to discuss the range of applicability of MF theory as a function of system dimensionality $d$. As mentioned in the Introduction, MF theory is expected to hold when the healing length $\hbar/\sqrt{mg\rho}$ (with $g=\widetilde{u}(0)\approx\sigma^{d}\epsilon$), fixing the length scale above which collective physics dominates over single-particle physics Pitaevskii3 , is much larger than the average interparticle separation $\rho^{-1/d}$ note . In 3D, this leads to $\rho\sigma^{3}\epsilon/e_{0}\ll(e_{0}/\epsilon)^{2}$, which corresponds to a density range that is wider the weaker the interaction strength. However, in 1D the MF regime is rather $\rho\sigma\epsilon/e_{0}\gg(\epsilon/e_{0})^{2}$, and the approximation improves with increasing density.

The way to solve Eqs. (2.5) for fixed values of $\rho$ and $a$ is by iteration Kunimi : at each step of the procedure, $S_{G}$ is first estimated from the $c_{G}$ coefficients computed at the previous step; the resulting linear system is then solved, determining eigenvalues $\mu_{n}$ and normalized eigenvectors. Next, the coefficients are updated to the eigenvector with minimum energy. The iterative process comes to an end when self-consistency is attained. The final task to accomplish is the optimization of the lattice parameter $a$, which is stopped when its value is determined to six decimal places. Once the specific energy $e$ has been computed as a function of $\rho$, the identification of the stable ground state at pressure $P$ proceeds via the minimization of the generalized enthalpy $\widetilde{h}(\rho;T=0,P)=e(\rho)+P/\rho$, which contextually determines the equilibrium density as $\rho_{\rm eq}(P)={\rm argmin}\,\widetilde{h}(\rho)$.

Kunimi and Kato have solved Eqs. (2.5) for PSM bosons in two dimensions (2D) Kunimi , showing that the high-density ground state is a triangular crystal. Macrì et al. Macri2 have tested MF results by Monte Carlo simulation, finding that the condensate is indeed only weakly depleted in the fluid region and that the exact freezing point is close to the theoretical estimate. In Ref. Prestipino3 we have extended the ground-state calculation to other 2D and 3D lattices, by employing an accurate variational form of $\psi$ — written as a sum of Gaussians centered at the lattice sites — that reproduces MF data to a high degree of accuracy. By this method, we have shown that the $T=0$ phase diagram of selected 3D potentials can be reconstructed with modest computational effort. Here, we make the same ansatz on the shape of the 1D wave function. While performing well in comparison with unconstrained MF theory, Gaussian variational theory has the distinct virtue of allowing a number of analytic shortcuts that considerably simplify extracting physical predictions from MF theory.

Using the variational method, we represent the single-particle state by the real-valued wave function

[TABLE]

where $C_{\alpha}$ is a suitable normalization constant (observe that the fluid phase, where $\psi=1/\sqrt{V}$, is recovered as a special case of (2.9), for $\alpha\rightarrow 0$). Two variational parameters are present in Eq. (2.9), i.e., $\alpha$ and $a$, related to the width and periodicity of the Gaussians, respectively (we stress that $a$ is an adjustable quantity as well, so as to ensure the possibility of cluster-crystal states). The best parameters ($\overline{\alpha}$ and $\overline{a}$) are those minimizing the restriction ${\cal E}(\alpha,a;\rho)$ of functional (2.6) to the set of functions (2.9); once $\overline{\alpha}$ and $\overline{a}$ have been computed for each density, the energy per particle is given by $e(\rho)={\cal E}(\overline{\alpha}(\rho),\overline{a}(\rho);\rho)$ (there is an energy branch for the crystal and another, $e(\rho)=(\widetilde{u}(0)/2)\rho$, for the fluid). Denoting $N_{c}$ the number of crystal cells, the average number of particles in a cell is $N/N_{c}=(N/V)(V/N_{c})=\rho a$. Hence, in a MF setting a cluster crystal is a crystalline state with $\rho a>1$. We denote VT1 the variational theory based on Eq. (2.9); a different variational approximation, denoted VT2, will be considered below.

In explicit terms,

[TABLE]

Although it does not admit an expression in terms of elementary functions, $I(\alpha)$ is related to a Jacobi theta function [cf. Eq. (A.16) in the Appendix]:

[TABLE]

The periodic function $\psi(x)$ can also be written as a Fourier series:

[TABLE]

with coefficients Prestipino3

[TABLE]

From the normalization condition $\sum_{G}\psi_{G}^{2}=1$ we derive another expression for $I(\alpha)$, namely

[TABLE]

Equation (2.14) proves useful to develop a low-$\alpha$ expansion of the energy functional (see Sec. III).

The advantage of the Gaussian series (2.9) over the more general function (2.4) is to allow analytic manipulations that considerably simplify the energy functional ${\cal E}(\alpha,a;\rho)$. Repeating the same steps followed in Ref. Prestipino3 , we first obtain a closed-form expression for the zero-point kinetic energy:

[TABLE]

As for the potential-energy functional, it simplifies to Prestipino3 :

[TABLE]

with

[TABLE]

As $\alpha$ varies from 0 to infinity, the ratio $J(\alpha)/I(\alpha)$ grows monotonically from 0 to 1.

Instead of the ansatz (2.9), we can directly express the square of $\psi$, that is $\eta(x)=V\psi^{2}(x)$ [see Eq. (2.7)], as a (normalized) sum of Gaussians centered on the lattice positions:

[TABLE]

(notice the different form of the exponent with respect to (2.9), made in order to provide a more meaningful comparison between VT1 and VT2 in the high-density limit). Since $\eta(x)$ is periodic, it can be expanded as a Fourier series:

[TABLE]

Denoting ${\cal C}$ any crystalline cell and $R=na$, the Fourier coefficient $\eta_{G}$ is given by

[TABLE]

Hence,

[TABLE]

Normalization is clearly satisfied, since $(1/V)\int_{V}{\rm d}x\,\eta(x)=\sum_{G}\eta_{G}\delta_{G,0}=\eta_{0}=1$.

We now compute the specific potential energy:

[TABLE]

For the PSM,

[TABLE]

For any $u(x)$, ${\cal E}_{\rm pot}$ depends on $\alpha$ only through the quantity $e^{-\pi^{2}/\alpha}$, which increases monotonically from 0 to 1 when $\alpha$ runs from 0 to infinity.

The kinetic energy per particle is given by the first term in Eq. (2.8):

[TABLE]

To proceed further, one notices that $\eta(x)$ is intimately related to a Jacobi theta function [see Eqs. (2.21) and (A.17)]:

[TABLE]

A graph of this function is plotted in Fig. 1. It is quite remarkable that a simple expression exists for the logarithmic derivative of $\vartheta_{3}$, see Eq. (A.19). Using this formula we can write:

[TABLE]

Since

[TABLE]

we finally obtain:

[TABLE]

Like ${\cal E}_{\rm pot}$, also ${\cal E}_{\rm kin}$ depends on $\alpha$ through $e^{-\pi^{2}/\alpha}$, while its $a$-dependence is simply ${\cal E}_{\rm kin}\propto a^{-2}$.

In the next Section we show that VT1 and VT2 give very similar energies; moreover, VT1 and VT2 also share with GP theory the same transition point.

III Results

Here, we provide results obtained for a few instances of 1D softly-repulsive bosons at $T=0$ using three MF theories, namely GP theory and the variational theories introduced in Sec. II. All theories agree in predicting a continuous quantum transition from a fluid phase to a cluster-crystal phase. While the GP approximation provides by construction the best condensate wave function, i.e., the one with the lowest energy, we shall see that VT1 is indeed as accurate in describing the transition behavior as GP theory.

1 Assessment of the variational approximations

For any given value of $\rho$, we solve the system of equations (2.5) for fixed $a$ by assuming that $c_{K}=0$ for $K=(2\pi/a)n$ and $|n|>8$ (nothing changes if this threshold were rather 12). We cyclically perform the diagonalization of the resulting $17\times 17$ Hermitian matrix of coefficients within the iterative procedure described in Sec. II, until self-consistency is reached. In the end, $a$ is optimized until its value is determined to five decimal places. Then, we solve VT1 and VT2, looking for the minimum of the energy functional on a grid of $(\alpha,a)$ values covering the region where the absolute minimum of ${\cal E}$ lies. The spacing of the grid is progressively reduced around the minimum, until its location is determined to $10^{-6}$ precision.

We show results for PSM bosons on a line in the left panel of Fig. 2, where the excess energy $\Delta{\cal E}={\cal E}-\rho\sigma\epsilon$ is plotted as a function of $\rho$. Below $\rho_{c}\lesssim 10.53$ (units of $e_{0}\epsilon^{-1}\sigma^{-1}$), the minimum $\Delta{\cal E}$ is invariably zero for all theories (fluid phase); above $\rho_{c}$, the minimum excess energy is a negative number, and the system phase is a crystal (we better discuss the nature of this crystal in the following Sec. III.2). We observe that the shape of the best single-particle state is nearly identical for GP theory and VT1 — see the right panel of Fig. 2, where the condensate wave functions for all theories are plotted side by side for $\rho=11$ and 13.

We have verified that the same degree of similarity between GP theory and VT1 also holds for the softened van der Waals (SVDW) repulsion, $u(r)=\epsilon/[1+(r/\sigma)^{6}]$, which is the same interaction investigated in Ref. Rossotti . Again, the transition threshold turns out to be the same in both theories ($\rho_{c}\lesssim 20.65$).

2 Analysis of the transition region

The results of Sec. III.1 indicate that the phase-transition threshold at $T=0$ of 1D penetrable bosons is by all evidence identical in GP theory and VT1 (we provide an explanation of this fact at the end of this Section). To better inquire into the system behavior near the transition, as well as to uncover similarities and differences between the various approaches, we derive below a small-$\alpha$ expansion of the VT1 and VT2 energy functionals that suffices for all purposes.

Let us first consider VT1. To calculate the energy [Eqs. (2.15) and (LABEL:eq2-16)], we need the perturbative expansions of $I(\alpha)$ [Eq. (2.14)] and $J(\alpha)$ [Eq. (2.17)] around $\alpha=0$. In this respect, it comes useful to express these functions in terms of Jacobi theta functions:

[TABLE]

Combining Eqs. (3.2.1) and (A.18) we obtain:

[TABLE]

whence the following small-$\alpha$ expansion of the difference in energy between the crystalline and fluid solutions:

[TABLE]

For the sake of clarity, let us consider the case of PSM bosons. For small $\alpha$, an approximation sufficient for the analysis of the transition behavior is:

[TABLE]

with $X=e^{-\pi^{2}/\alpha}$ and

[TABLE]

Notice that $r$ and $w$ are explicit functions of $a$ and $\rho$. The extremal points of $\Delta{\cal E}(X)$ are $\overline{X}=0$ and (if $r<0$) the non-zero root of $\Delta{\cal E}^{\prime}(X)=0$, that is $\overline{X}=\sqrt{-r/(2w)}$, with specific energies $\Delta{\cal E}_{\rm min}(a,\rho)=0$ and $-r^{2}/(4w)$, respectively (it turns out that $w>0$ in the relevant range of $a$ and $\rho$ values). The non-trivial solution exists providing that $r<0$, namely $\rho>\rho_{0}$, with

[TABLE]

Since $r=4\widetilde{u}(2\pi/a)(\rho-\rho_{0})$, it follows that

[TABLE]

The equilibrium lattice constant, $\overline{a}(\rho)$, is the one providing the minimum value of $\Delta{\cal E}_{\rm min}(a)$ for the given $\rho$.

In order that $\rho_{0}>0$, it is sufficient that $1<a/\sigma<2$; in this interval, $w>0$ as well (in fact, there are infinite other intervals where $\sin(2\pi\sigma/a)<0$, namely $1/2<a/\sigma<2/3,1/3<a/\sigma<2/5,\ldots$, but these other $a$ provide much larger values of $\rho_{0}$ [see Eq. (3.2.6)] and, above this density, also energy minima higher than those in the interval $1<a/\sigma<2$). The smallest density above which the energy is negative is the minimum of $\rho_{0}(a)$, i.e., $\rho_{0}(a_{c})\equiv\rho_{c}$. For $y=-x^{3}\sin(2\pi/x)$, the derivative $y^{\prime}\geq 0$ for

[TABLE]

(the other roots of Eq. (3.2.8) fall outside the range from 1 to 2). Hence, the transition from fluid to crystal occurs for $\rho=\rho_{c}=10.524990629\ldots$; at this density, $\overline{X}$ switches continuously from 0 (stable fluid) to $\sqrt{-r/(2w)}\propto\sqrt{(\rho-\rho_{c})/\rho_{c}}$ (stable crystal). For densities larger than $\rho_{c}$, the minimum of energy occurs at a $\overline{a}$ smaller than $a_{c}$, see Fig. 3. Slightly above $\rho_{c}$, the behavior of $\overline{\alpha}$ and of the excess energy are, up to a ${\cal O}(1)$ factor:

[TABLE]

In particular, the order parameter $\overline{\alpha}$ is continuous for $\rho=\rho_{c}$. At the transition, the crystal pressure equals that of the fluid, as it follows from the general relation $P(\rho)=\rho^{2}e^{\prime}(\rho)$ and the second of Eqs. (3.2.9). The transition pressure is $P_{c}=\sigma\epsilon\rho_{c}^{2}=110.7754\ldots$ (units of $e_{0}^{2}\epsilon^{-1}\sigma^{-1}$). The phase transition is continuous since the energy and its derivative are continuous at $\rho=\rho_{c}$ (the phase with the minimum enthalpy is the one with the minimum energy). At the transition, the average number of particles per site is $\rho_{c}a_{c}=16.2157\ldots$ (hence, the dense phase is a cluster crystal).

Similar results hold for VT2. It follows from Eqs. (2.23) and (2.28) that

[TABLE]

As before, for small values of $X=\exp\{-\pi^{2}/(2\alpha)\}$ (which is different from the definition of $X$ given before) the same approximation (3.2.4) holds, with $w>0$ independent of the density and a $r$ value which is one fourth of that for VT1. This implies that $\rho_{c}$ and $a_{c}$ are identical for VT1 and VT2. In particular, near $\rho=\rho_{c}$ the behavior of $\overline{\alpha}$ and of the excess energy are the following:

[TABLE]

We find an unexpected outcome when computing the isothermal compressibility at $T=0$,

[TABLE]

For $\rho=\rho_{c}$, $K_{T}^{-1}$ has the following values in the two phases:

[TABLE]

Hence, $K_{T}$ undergoes a jump at the transition and, right at the transition point, the crystal is more compressible than the fluid.

Summarizing, all three theories predict the same transition density $\rho_{c}$, whereas small differences arise for $\rho>\rho_{c}$ (see Fig. 3). We can appreciate from this figure that the VT1 approximation is superior to VT2, since it gives a smaller energy; moreover, VT1 has practically the same accuracy of GP theory.

To explain why VT1 and VT2 give exactly the same transition density as GP theory, we consider the nature of fluid excitations in the latter theory. As well known (see, e.g., Kunimi ; Prestipino3 ; Zloshchastiev ), the spectrum of these excitations is Bogoliubov-like:

[TABLE]

According to this dispersion law, if $\widetilde{u}(k)$ is negative in some range of $k$ (as occurs for any bounded interaction that is “fatter” than Gaussian) then the system is superfluid by Landau’s argument. In this case, above a certain density a maxon peak develops in $\omega(k)$, followed at larger $k$ by a roton minimum; increasing the density further, the roton eventually softens and the fluid becomes unstable towards the formation of a density wave (this happens when the quantity within parentheses in Eq. (3.2.14) vanishes). Putting $k=2\pi/a$, roton softening first occurs for a density and a value of the wavelength $a$ such that

[TABLE]

where the rationale behind the second equation is that the density $\rho_{0}(a)$ solving the first equation be as low as possible. These are exactly the same conditions for the occurrence of clusterization in variational theory. In particular, for the SVDW potential the Fourier transform reads

[TABLE]

(this is obtained by evaluating with the residue theorem the integral of $e^{ikz}/(1+z^{6})$ over a large semi-circular contour inscribed in the ${\rm Re}\,z\geq 0$ half-plane). By numerically solving Eqs. (3.2.15) we confirm that clusterization of SVDW bosons in MF theory occurs at $\rho_{c}=20.64654\ldots$, which is the same condition for roton softening in the $\epsilon\rightarrow 0$ limit quoted by Rossotti and coworkers. It is worth observing that the MF locus for roton softening gives a good approximation to the exact transition line even far away from the MF limit (see Fig. 1 in Ref. Rossotti ).

3 PSM bosons: High-density limit

After highlighting the features of clusterization at $T=0$ of 1D soft-core bosons through an analysis of the small-$\alpha$ limit of the energy functional, we now consider the opposite limit of very large $\alpha$ values, which corresponds to high densities, focusing on the case of PSM bosons. It turns out that a duality property of Jacobi theta functions makes this limit accessible analytically.

Starting with VT1, we first approximate the kinetic energy (2.15) by noting that, for $\alpha\gg 1$,

[TABLE]

We immediately obtain:

[TABLE]

Then, we estimate the ratio $[J(\alpha)/I(\alpha)]^{2}$ appearing in the potential-energy formula, Eq. (LABEL:eq2-16). From the second of Eqs. (3.2.1), using a duality property of theta functions [Eq. (A.20)] we obtain:

[TABLE]

so that

[TABLE]

Next, we evaluate for PSM bosons the sum

[TABLE]

with $z=2\pi\sigma/a$. In order to estimate the large-$\alpha$ limit of Eq. (3.3.5), we note that:

[TABLE]

By another duality formula [Eq. (A.21)] we get:

[TABLE]

For $1\leq a/\sigma\leq 2$, the value of $z$ is between $\pi$ and $2\pi$. In this range, the leading terms in the expansion (3.3.7) are (in equal measure) the first and the second one (while the 1 can be ignored), thus obtaining:

[TABLE]

Noting that $f(\pi)=0$, we have:

[TABLE]

The integral returns error functions, whose limiting behavior for large values of the argument is:

[TABLE]

After obvious steps we eventually find:

[TABLE]

Then, we examine the asymptotic behavior of

[TABLE]

with $z=2\pi\sigma/a$. By the same above considerations, we obtain:

[TABLE]

Integrating (3.3.13) from $\pi$ and $z$, and then plugging the result in (3.3.12), we arrive at:

[TABLE]

Putting Eqs. (3.3.2), (3.3.4), (3.3.11), and (3.3.14) together, we obtain the sought-for high-density approximation of the VT1 energy functional.

The treatment is simpler for VT2. As far as kinetic energy is concerned, we start from $\eta(x)=\vartheta_{3}[\pi x/a,\exp\{-\pi^{2}/(2\alpha)\}]$. Observing that

[TABLE]

for $-a/2\leq x\leq a/2$ we obtain

[TABLE]

whence we find [cf. Eq. (2.24)]

[TABLE]

As for the potential energy, for PSM bosons it equals [see Eq. (2.23)]

[TABLE]

with $z=\pi\sigma/a$ and

[TABLE]

On the other hand,

[TABLE]

For $1\leq a/\sigma\leq 2$, $z$ falls between $\pi/2$ and $\pi$. In this interval, the leading term for $\alpha\gg 1$ is the first one in the expansion (3.3.20), while the 1 can be ignored, thus arriving at:

[TABLE]

Integrating from $\pi/2$ (where $f$ vanishes) and $z$, we thus obtain:

[TABLE]

and

[TABLE]

The sum of Eqs. (3.3.17) and (3.3.23) is the asymptotic expression of the VT2 energy functional.

Using the simplified energy functionals, we obtain the data plotted in Fig. 4. It turns out that, for high densities, VT1 and VT2 give practically the same optimal values of the variational parameters. For all densities above $\approx 80$, these values are hardly distinguishable from those extracted from the original functionals. In particular, VT1 energy is imperceptibly smaller than VT2 energy.

4 PSM bosons: Supersolidity of the cluster phase

Consider a quantum solid of cylindrical shape, set in uniform rotation around the axis. Leggett Leggett1 ; Leggett2 has proposed to call superfluid fraction of the system the quantity

[TABLE]

where $I$ is the moment of inertia around the cylinder axis and $I_{0}$ its classical value. It turns out that, like a superfluid, also a quantum solid may exhibit an anomalous response to axial rotations: for low rotation speed, part of the solid may stand still, with the result that $I<I_{0}$ and therefore $f_{s}>0$. In this case, the system is called a supersolid. The specificity of 1D is that, strictly speaking, an infinite crystal does not exist at $T=0$; however, the notion of supersolidity may still be considered for a finite, but large, 1D system in the CLL phase, e.g., for clusters arranged in a ring that rotates around its axis.

Leggett has derived an exact formula for the superfluid fraction of a 1D solid of identical bosons at $T=0$ Leggett2 :

[TABLE]

where $a$ is the volume of the crystalline cell ${\cal C}$ and $\eta(x)=V\psi^{2}(x)$. The same result (3.4.2) has been obtained by Sepulveda et al. within an approximate theory of the supersolid phase Sepulveda .

For the sake of clarity, let us take PSM bosons in 1D. We first consider VT2, for which $\eta(x)=\vartheta_{3}(\pi x/a,\exp\{-\pi^{2}/(2\alpha)\})$. At low density, we can write:

[TABLE]

Hence, we find:

[TABLE]

and finally:

[TABLE]

This indicates that, close to the transition point, the superfluid fraction of the crystal varies as [cf. the first of Eqs. (3.2.11)]:

[TABLE]

Hence, the crystalline solid is a supersolid whose superfluid fraction is exactly one at the transition, then reducing progressively on compression. It is natural to ask whether $f_{s}$ eventually vanishes at a certain large value of the density, and the system then becomes a normal solid. We shall see that the answer is in the negative, at least within variational theory.

In the limit of high densities, Eq. (3.3.15) allows us to write $1/\eta(x)$ as the ratio between two large quantities. However, we were not able to put this ratio in the form of a rapidly-convergent series. It is much simpler to obtain a positive lower limit for $f_{s}$, and this way conclude that the crystalline solid is, like in higher dimensions Prestipino3 , a supersolid also for very large density. Called $\eta_{\rm min}=\eta(a/2)$ and $\eta_{\rm max}=\eta(0)$ the minimum and maximum value of $\eta(x)$ in the cell, we can write:

[TABLE]

We derive from Eq. (3.3.16) that:

[TABLE]

whence it follows:

[TABLE]

as anticipated.

In VT1, where $\sqrt{V}\psi(x)=C_{\alpha}^{\prime}\vartheta_{3}[\pi x/a,\exp(-\pi^{2}/\alpha)]$ and $\eta(x)=V\psi^{2}(x)$, near the transition point the following expansion holds:

[TABLE]

from which we obtain:

[TABLE]

Using Eqs. (2.13), (3.2.1), and (A.18), we estimate:

[TABLE]

and finally:

[TABLE]

where the first term is linear in $\rho-\rho_{c}$ and the second is quadratic.

In the opposite limit of high densities,

[TABLE]

and then

[TABLE]

In conclusion, according to both VT1 and VT2 the crystal is supersolid at all densities.

IV Conclusions

Pressure-driven clusterization of a fluid of soft-core bosons at $T=0$, i.e., the emergence of clumps of overlapping particles (clusters) under compression, is among the simplest examples of a quantum transition. At variance with two or three dimensions, where the formation of clusters is always accompanied by the appearance of crystalline order Pomeau ; Saccani ; Prestipino3 , in 1D a no-go theorem by Pitaevskii and Stringari Pitaevskii2 excludes the possibility of long-range order in the thermodynamic limit, implying loss of crystalline and phase coherence at large distances. However, (truncated) crystalline order is recovered in a large, but finite, 1D system confined in an elongated trap or placed in a narrow torus.

In this paper we apply three different MF theories to the study of clusterization in 1D, representing the ground state of the system as a pure condensate. At variance with 3D, where MF theory holds for small $\rho\epsilon$ values (with $\epsilon$ denoting the interaction strength), a 1D system approaches the weak-coupling regime for decreasing $\epsilon/\rho$. We further assume crystalline ordering of the dense phase. Besides GP theory, which is equivalent to selecting the best MF state, we consider two variational approximations for the single-particle wave function $\psi$: in one case, $\psi$ is written as a sum of Gaussians (a two-parameter ansatz); in the other case, by a similar sum we represent $\psi^{2}$. The virtue of variational theory is that the energy functional is written in almost closed form, which allows us to derive a number of analytic predictions.

In one dimension, the freezing transition turns out to be continuous and occurs at the highest density at which the fluid is still superfluid. As a rule, the crystalline ground state is a cluster crystal, meaning that the average site occupancy is larger than one. Moreover, the crystal is supersolid, meaning that the moment of inertia is smaller than for a classical solid of same mass and size. In more physical terms, supersolidity can be ascribed to the delocalization of the condensate wave function over the whole crystal and, particularly, to a non-zero probability of observing a particle in the interstitial region.

Appendix A Some useful formulas

In this Appendix we collect a few formulas relative to Jacobi theta functions. These are special functions with a relation to elliptic functions. For a complex variable $z$ and a complex number $q$ of modulus less than 1, theta functions are defined as jacobitheta :

[TABLE]

Clearly, a less symmetric form exists for each function, e.g.,

[TABLE]

For $z=0$, the $\vartheta_{2}$ and $\vartheta_{3}$ functions have the following obvious expansions around $q=0$:

[TABLE]

Interestingly, there is a way to express the logarithmic derivative of $\vartheta_{3}$ in the form of a series jacobitheta :

[TABLE]

A similar formula exists for each theta function.

A duality property holds for Jacobi theta functions nistdigitallibrary , by which a particular theta function for $q\lesssim 1$ is related to another theta function for $q\gtrsim 0$. For example,

[TABLE]

where, e.g., $\vartheta_{2}(z|\tau)$ stands for $\vartheta_{2}(z,q)$, with $q=\exp(i\pi\tau)$ and ${\rm Im}\,\tau>0$, whereas $\tau^{\prime}=-1/\tau$. Another useful formula is:

[TABLE]