Cornwall-Jackiw-Tomboulis effective field theory to nonuniversal equation of state of an ultracold Bose gas

The equation of state (EOS) of quantum physical systems, arising from quantum fluctuations, occupies a pivotal position at the very core of quantum many-body physics [1, 2, 3]. A prototypical illustration lies within the realm of quantum droplets [4, 5, 6] in the context of ultracold quantum gases. The basic principle underlying the stabilization of quantum droplets hinges on a delicate equilibrium between the attractive mean-field force and the repulsive force emanating from quantum fluctuations as described by the EOS [7, 8, 9, 10, 11, 12]. Consequently, delving into the EOS can provide invaluable insights into the fundamental properties of quantum droplets, as evidenced across diverse systems, spanning from dipolar Bose gases [7, 8, 9] to Bose-Bose mixtures [10, 11, 12].

The interest in the EOS of ultracold quantum gases, particularly underscored by the aforementioned quantum droplets, can be traced back to the seminal works [13, 14] in the 1950s. Particularly, the Lee-Huang-Yang (LHY) correction to the EOS of weakly-interacting bosonic systems was first explored in Ref. [14]. In more details, the chemical potential has been derived to be proportional to $[32/3\sqrt{\pi}]\left(\rho a_{\text{s}}^{3}\right)^{1/2}$ with $a_{\text{s}}$ [15, 16] and $\rho$ being the $s$-wave scattering length and the atomic density, respectively. Furthermore, the next-order correction to the energy density incorporates the universal quantum effect stemming from three-body correlations, yielding a term of $\frac{8}{3}\left(4\pi-3\sqrt{3}\right)\rho a_{\text{s}}^{3}\ln\left(\rho a_{\text{s}}^{3}\right)$, as outlined in Ref. [17]. Up to now, the highest-order term in the EOS has been derived in Ref. [18], which is given by $[-512/3\pi]\left(\rho a_{\text{s}}^{3}\right)+[8192/9\pi^{3/2}]\left(\rho a_{\text{s}}^{3}\right)^{3/2}$. At present, there are the endless and ongoing research interests and efforts in obtaining beyond the mean-field terms of EOS of weakly-interacting bosonic systems motivated by the advent of well-controlled mixtures of quantum gases with tunable interaction strengths [19, 20, 21, 22, 23].

The EOS obtained in Refs. [13, 14, 15, 16, 17, 18] are characterized by the universality, i.e. in more details, a single parameter $a_{\text{s}}$, eloquently encapsulates both the intricacies of the two-body interaction and, by extension, the overarching physics governing the many-body systems [24]. In sharp contrast, the nonuniversal EOS depends on other than the parameter $a_{\text{s}}$. The EOS of Bose gases becomes nonuniversal when the finite-range effects of the interatomic potential [25, 26, 27, 28, 29] is taken into account. In Refs. [26, 30], the analytical expressions of nonuniversal EOS are derived as $[-64\sqrt{\pi}r_{\text{s}}/a_{\text{s}}]\rho\left(\rho a_{\text{s}}^{3}\right)^{3/2}$, featuring a nonuniversal LHY term in the quantum depletion. Given the tunability of the scattering length $a_{\text{s}}$ [31, 32] through magnetic and optical Feshbach resonances [33, 34, 35] in ultracold atomic gases, the nonuniversal consequences stemming from the finite-range parameter $r_{\text{s}}$ are of paramount importance. Consequently, an immediate challenge is referred as to calculating the next-to-LHY-order correction for the nonuniversal EOS.

The second impetus behind this paper stems from the measurement of collective excitation frequencies, which has emerged as an indispensable and highly accurate instrument for delving into the intricate EOS of atomic Bose-Einstein condensates (BECs) with unprecedented precision [15, 36, 19, 20, 21, 22, 23]. This approach not only solidifies the validity of mean-field predictions but also stands as an exceedingly potent technique to explore effects transcending the mean-field paradigm [37, 38]. From a theoretical standpoint, a gaseous BECs system can be aptly modeled by a single macroscopic wave function, thereby facilitating the derivation of clear-cut hydrodynamic formulations that yield analytical or semi-analytical insights into the dynamical characteristics of BECs systems. Thus, a timely and natural question arises as to observe the frequency shifts in the collective excitations induced by the nonuniversal beyond-LHY EOS, although tuning the finite-range interatomic interactions in this case remains experimentally challenging.

In this work, using Cornwall-Jackiw-Tomboulis (CJT) effective field theory [39, 40, 41, 42], we are interested in the nonuniversal EOS [43, 26, 30, 28, 29] of a three-dimensional (3D) Bose gas with finite-range effective interactions at absolute zero temperature. Accordingly, we derive the analytical expressions of the quantum depletion and the chemical potential up to the order of $\left(\rho a_{\text{s}}^{3}\right)^{2}$. Our results not only include Salasnich’s remarkable EOS, but also give rise to nonuniversal beyond-LHY terms. Moreover, we explore the physical consequences of the nonuniversal beyond-LHY terms in EOS. As such, we extend the superfluid hydrodynamic equations by incorporating the aforementioned next-to-LHY corrections to the nonuniversal EOS. Leveraging these refined hydrodynamic equations, we delve into the fractional frequency shifts in the breathing modes, which are observable within the current experimental facilities. Observing this frequency shifts induced by the nonuniversal beyond-LHY EOS paves the way for a deeper understanding of quantum fluctuation of quantum many-body systems.

The paper is structured as follows. In Sec. II, we revisit the key principles of the CJT effective field theory and subsequently derive the corresponding effective potential for the model system. In Sec. III, we utilize the CJT effective potential within a two-loop approximation to deduce analytical expressions for the nonuniversal EOS. Sec. IV presents the fractional shift in the breathing mode frequency, aided by the next-to-LHY corrections to the nonuniversal chemical potential. Finally, in Sec. V, we provide a comprehensive summary of our paper and discuss the potential experimental conditions for realizing our proposed scenario.

In this work, we are interested in a 3D weakly-interacting Bose gas, paying particular attention to the finite-range effects of the interatomic potential [27]. To investigate this system, we employ the path-integral formalism [1, 2], and the corresponding Euclidean partition function of the model system is presented below

$$\mathcal{Z}=\int\mathcal{D}\left[\mathbf{\Phi},\mathbf{\Phi}^{*}\right]\exp\left\{-\frac{S\left[\mathbf{\Phi},\mathbf{\Phi}^{*}\right]}{\hbar}\right\},$$

(1)

with the action functional $S\left[\mathbf{\Phi},\mathbf{\Phi}^{*}\right]=\int_{0}^{\beta\hbar}d\tau\int d^{3}\mathbf{r}\mathcal{L}$ in Eq. (1). Here, the concrete Lagrangian density $\mathcal{L}$ denotes as follows [43, 26, 30, 44]

$\displaystyle\mathcal{L}\left[\mathbf{\Phi},\mathbf{\Phi}^{*}\right]=\mathbf{\Phi}^{*}\left(\mathbf{r},\tau\right)\left[\hbar\partial_{\tau}-\frac{\hbar^{2}\nabla^{2}}{2m}-\mu\right]\mathbf{\Phi}\left(\mathbf{r},\tau\right)+\frac{g_{0}}{2}\left|\mathbf{\Phi}\left(\mathbf{r},\tau\right)\right|^{4}-\frac{g_{2}}{2}\left|\mathbf{\Phi}\left(\mathbf{r},\tau\right)\right|^{2}\nabla^{2}\left|\mathbf{\Phi}\left(\mathbf{r},\tau\right)\right|^{2}.$

(2)

In Equation (2), the complex field $\mathbf{\Phi}\left(\mathbf{r},\tau\right)$ represents the atomic bosons, varying in both space $\mathbf{r}$ and imaginary time $\tau$. Here, $\mu$ signifies the chemical potential, while $\beta=1/k_{\text{B}}T$ defines the inverse of the thermal energy scale, with $k_{\text{B}}$ representing the Boltzmann constant and $T$ denoting the temperature of the BECs. The parameters $g_{0}=4\pi\hbar^{2}a_{\text{s}}/m$ and $g_{2}=2\pi\hbar^{2}a_{\text{s}}^{2}r_{\text{s}}/m$ [43], with $a_{\text{s}}$ and $r_{\text{s}}$ being the $s$-wave scattering length and the finite-range length respectively.

Note that the LHY correction to nonuniversal EOS of Lagrangian density functional (2) has already been derived [26, 30] within one-loop approximation, i.e. rewriting the partition function (1) in the form of $\mathcal{Z}=e^{-\beta\mathcal{V}V_{\text{eff}}[\phi]}$, with $\mathcal{V}$ being the volume of the system. In this context, the effective potential $V_{\text{eff}}[\phi]$ [45] within one-loop approximation solely relies on $\phi(x)$, which is the expected value of the quantum field $\hat{\mathbf{\Phi}}(x)$.

In contrast, the emphasis and value of this paper lie in employing the CJT [40, 42, 41, 46, 45] effective action approach or the two particle-irreducible framework, surpassing the limitations of the one-loop approximation [45]. CJT theory enables us to compute the beyond-LHY corrections to the nonuniversal EOS for the Lagrangian density functional (2). A comprehensive introduction to the CJT effective field theory is provided in Appendix A.

The central step of CJT effective field theory [39, 40, 41, 42] is to obtain the effective potential denoted as $V_{\text{eff}}[\phi,G]$ (see Eq. (31) in Appendix A) meticulously by taking both the background field $\phi(x)$ and the dressed propagator $G(x,y)$ into account. Here, $G(x,y)$ is a potential expectation value of the time-ordered product $T\hat{\mathbf{\Phi}}^{\dagger}(x)\hat{\mathbf{\Phi}}(y)$. Then, physical solutions are ascertained by ensuring that the generalized effective potential satisfies the stationarity conditions: $\delta V_{\text{eff}}/\delta\phi(x)=0$, $\delta V_{\text{eff}}/\delta G(x,y)=0$. Finally, at the core of the CJT effective field theory lies in the self-consistent loop expansion of the effective potential $V_{\text{eff}}$, intricately tied to the full propagator $G$. This expansion offers a powerful tool for systematically exploring higher-order corrections, thereby enhancing the accuracy and predictive capabilities of the system with finite-range effects taken into account. We remark that our work, together with Refs. [26, 30], provides a reasonable description of the nonuniversal EOS of a 3D interacting Bose gas with finite-range effects of the interatomic potential.

The subsequent goal of Sec. II is to derive the effective potential $V_{\text{eff}}$ for functional (2) within the framework of CJT effective field theory. Then, the obtained effective potential $V_{\text{eff}}$ is used to calculate the nonuniversal beyond-LHY EOS.

The starting point of the CJT effective field theory commences with expressing the field $\mathbf{\Phi}$ of functional (2) as a superposition of the condensate field $\phi_{0}$ and real fluctuation fields of $\phi_{1}$ and $\phi_{2}$, i.e. $\mathbf{\Phi}=\left(\phi_{0}+\phi_{1}+i\phi_{2}\right)/\sqrt{2}$. A rigorous methodology to attain $V_{\text{eff}}$ involves executing a double Legendre transformation on the action functional $S\left[\mathbf{\Phi},\mathbf{\Phi}^{*}\right]$, subjecting it to the conditions $\delta V_{\text{eff}}/\delta\phi_{0}=0$ and $\delta V_{\text{eff}}/\delta G=0$, where $G$ represents a relevant variable (such as the two-point function or propagator) that may need to be optimized simultaneously with $\phi_{0}$. This ensures the effective potential accurately captures the dynamics of the system, incorporating both condensate and fluctuation effects.

Adhering closely to the established CJT effective field theory as outlined in Refs. [47, 48, 46, 49, 18], we proceed by performing Fourier transformations on the fluctuation fields $\phi_{1}$ and $\phi_{2}$ to map them into the momentum-frequency domain. Subsequently, we select the Luttinger-Ward functional [45], denoted as $\Phi\left[\phi_{0},G\right]$, as the foundational element. Consequently, the effective potential $V_{\text{eff}}$ corresponding to partition function in Eq. (1) can be analytically formulated as (see the detailed derivation in Eq. (39) in Appendix B)

	$\displaystyle V_{\text{eff}}\left[\phi_{0},G\right]=$	$\displaystyle-$	$\displaystyle\frac{\mu}{2}\phi_{0}^{2}+\frac{g_{0}}{8}\phi_{0}^{4}$		(3)
		$\displaystyle+$	$\displaystyle\frac{1}{2}\int_{\beta}\text{Tr}\left[\ln G^{-1}\left(k\right)+G_{0}^{-1}\left(k\right)G\left(k\right)-\mathbf{1}\right]$
		$\displaystyle+$	$\displaystyle\frac{3g_{0}}{8}\left(P_{11}^{2}+P_{22}^{2}\right)+\frac{g_{0}}{4}P_{11}P_{22}.$

In Equation (3), $k$ denotes the magnitude of the wave vector $\mathbf{k}$. The quantities $G^{-1}_{0}\left(k\right)$ and $G^{-1}\left(k\right)$ represent the inverse propagators within the one-loop and two-loop approximations, respectively. The notation $\int_{\beta}$ encapsulates the integration over momentum space combined with a summation over bosonic Matsubara frequencies, specifically given by $\int_{\beta}=\frac{1}{\beta}\sum_{n=-\infty}^{+\infty}\int\frac{d\mathbf{k}}{\left(2\pi\right)^{3}}$; $\omega_{n}=2\pi n\beta^{-1}$ are the Matsubara frequencies. The functions $P_{aa}=\int_{\beta}G_{aa}\left(k\right)$ are termed momentum integrals, with $a$ indexing the two fields $\left(1,2\right)$.

Minimizing the CJT effective potential $V_{\text{eff}}$ in Eq. (3) with respect to the components of the propagator $G\left(k\right)$, we obtain

$$G^{-1}\left(k\right)=G_{0}^{-1}\left(k\right)+\Sigma.$$

(4)

In Equation (4), the context form of $G_{0}^{-1}$ can be written as (see the detailed derivation in Eq. (38) in Appendix B.2)

$$\!\!G_{0}^{-1}\!\left(k\right)\!=\!\begin{bmatrix}\!\frac{\hbar^{2}k^{2}}{2m}\!-\!\mu\!+\!\frac{3g_{0}}{2}\phi_{0}^{2}\!+g_{2}\phi_{0}^{2}k^{2}&\!\!-\omega_{n}\\ \!\!\omega_{n}&\frac{\hbar^{2}k^{2}}{2m}\!-\!\mu\!+\!\frac{g_{0}}{2}\phi_{0}^{2}\end{bmatrix}.\!\!$$

(5)

Meanwhile, $\Sigma$ in Eq. (4) represents the self-energy matrix in the form of

$$\Sigma=\begin{bmatrix}\Sigma_{1}&0\\ 0&\Sigma_{2}\end{bmatrix},$$

(6)

with the matrix entries $\Sigma_{1}$ and $\Sigma_{2}$ reading

	$\displaystyle\Sigma_{1}$	$\displaystyle=$	$\displaystyle\frac{3g_{0}}{2}P_{11}+\frac{g_{0}}{2}P_{22},$		(7a)
	$\displaystyle\Sigma_{2}$	$\displaystyle=$	$\displaystyle\frac{3g_{0}}{2}P_{22}+\frac{g_{0}}{2}P_{11}.$		(7b)

Before proceeding with further calculations based on Eq. (3), we conduct a crucial verification to ensure that Eq. (3) can indeed reproduce the previous findings presented in Refs. [43, 26, 30]. Specifically, under the conditions where the propagator $G\left(k\right)$ reduces to its bare form $G_{0}\left(k\right)$ and both $P_{11}$ and $P_{22}$ vanish, the effective potential simplifies significantly to align with the one-loop approximation scenario. This simplification yields $V_{\text{eff}}\left[\phi_{0},G_{0}\right]=-\frac{\mu}{2}\phi_{0}^{2}+\frac{g_{0}}{8}\phi_{0}^{4}+\frac{1}{2}\int_{\beta}\text{Tr}\ln\left[G_{0}^{-1}\left(k\right)\right]$ within the context of the path-integral formalism. Then, after minimizing $V_{\text{eff}}$ with respect to order parameter $\phi_{0}$ and applying the thermodynamic relationship $\rho=-\partial V_{\text{eff}}/\partial\phi_{0}$, the nonuniversal EOS [43, 26, 30, 28, 29] can be obtained within one-loop approximation. Note that the CJT effective field theory not only introduces two-particle irreducible (2PI) terms in the last line of Eq. (3) but also modifies the propagator $G$ in a more intricate manner, affecting the second term of $V_{\text{eff}}$ in Eq. (3) in a complex and non-trivial way.

In the preceding Sec. II, we have delineated the framework of the CJT effective field theory. Moving forward, in Sec. III, our objective is to derive the explicit analytical expressions for the next-to-LHY-order correction to the nonuniversal EOS of a 3D Bose gas, utilizing the CJT effective potential of Eq. (3) within two-loop approximation. The starting point for this endeavor is to deduce the $V_{\text{eff}}$ in Eq. (3) as (refer to Appendix B for a more comprehensive derivation),

$$V_{\text{eff}}=V_{0}+\frac{1}{2}\int_{\beta}\text{\text{Tr}}\left[\ln G^{-1}\left(k\right)\right]+\frac{g_{0}}{8}\left(P_{11}^{2}+P_{22}^{2}\right)+\frac{3g_{0}}{4}P_{11}P_{22}+\frac{1}{2}\left(-\mu+\frac{3g_{0}}{2}\phi_{0}^{2}-M^{2}\right)P_{11}+\frac{1}{2}\left(-\mu+\frac{g_{0}}{2}\phi_{0}^{2}\right)P_{22}.$$

(8)

In Equation (8), $V_{0}=-\frac{\mu}{2}\phi_{0}^{2}+\frac{g_{0}}{8}\phi_{0}^{4}$, $P_{11}=\frac{\sqrt{2}m^{*3/2}M^{3}}{3\pi^{2}\hbar^{3}}\sqrt{\frac{m^{*}}{m}}$ and $P_{22}=\frac{-m^{*3/2}M^{3}}{3\sqrt{2}\pi^{2}\hbar^{3}}\sqrt{\frac{m}{m^{*}}}$ (see Eqs. (54a) and (54b) in Appendix B for calculation details) with $m^{*}=m/\left(1+\frac{2mg_{2}\phi_{0}^{2}}{\hbar^{2}}\right)$.

The pivotal parameter $M$ in Eq. (8) holds a crucial role in determining not only $P_{11}$ and $P_{22}$, but also the other constituent terms in Eq. (8). Subsequently, we embark on deriving the precise expression for $M$, as outlined below:

(i). The parameter $M$ in Eq. (8) fulfills the Schwinger-Dyson (SD) equation (as elaborated in detail in Eq. (51) in Appendix B),

$$-\mu+\frac{3g_{0}}{2}\phi_{0}^{2}+\frac{\sqrt{2}g_{0}m^{*3/2}M^{3}}{12\pi^{2}\hbar^{3}}\left(2\sqrt{\frac{m^{*}}{m}}-3\sqrt{\frac{m}{m^{*}}}\right)=M^{2}.$$

(9)

Equation (9) contains three variables, i.e. $M$, $\mu$, and $\phi_{0}$. To ascertain their values, we require two additional equations in conjunction with Eq. (9).

(ii). Next, we proceed to seek for the second equation between $M$, $\mu$, and $\phi_{0}$. Specifically, the chemical potential of $\mu$ satisfies the gap equation by setting $\delta V_{\text{eff}}/\delta\phi_{0}=0$ (for a detailed exposition, refer to Eq. (44) in Appendix B),

$$-\mu+\frac{g_{0}}{2}\phi_{0}^{2}+\frac{\sqrt{2}g_{0}m^{*3/2}M^{3}}{12\pi^{2}\hbar^{3}}\left(6\sqrt{\frac{m^{*}}{m}}-\sqrt{\frac{m}{m^{*}}}\right)=0.$$

(10)

(iii). Finally, the third equation originated by determining the atomic density of $\rho$, can be obtained by taking the first-order derivative of the $V_{\text{eff}}$ in Eq. (8) with respect to $\mu$,

$$\text{{\hbox{\rho}}}=-\frac{\partial V_{\text{eff}}}{\partial\mu}=\frac{\phi_{0}^{2}}{2}+\frac{\sqrt{2}m^{*3/2}M^{3}}{12\pi^{2}\hbar^{3}}\left(2\sqrt{\frac{m^{*}}{m}}-\sqrt{\frac{m}{m^{*}}}\right).$$

(11)

Equations (9), (10) and (11) constitute a closed set of equations. By eliminating the variables $\mu$ and $\phi_{0}$ from Eq. (9) with the aid of Eqs. (10) and (11), we obtain a cubic equation solely in terms of $M$,

$$M^{3}+\frac{\left(1+4\frac{mg_{2}}{\hbar^{2}}\rho\right)^{3/2}}{8\sqrt{2}a_{\text{s}}m^{*1/2}}M^{2}-\frac{3\pi^{2}\hbar^{3}\left(1+4\frac{mg_{2}}{\hbar^{2}}\rho\right)^{1/2}}{\sqrt{2}m^{*3/2}}\rho=0,$$

(12)

with $g_{0}=4\pi\hbar^{2}a_{\text{s}}/m$ and $m^{*}\simeq m/\left(1+4m\frac{g_{2}}{\hbar^{2}}\rho\right)$.

By solving Equation (12) using perturbation theory (details provided in Appendix C) and substituting $g_{2}$ with $2\pi\hbar^{2}a_{\text{s}}^{2}r_{\text{s}}/m$, we derive the analytical expression of $M$ expanded as the gas parameter $\rho a_{\text{s}}^{3}$

$$M=\sqrt{2\rho g_{0}}\left\{1-\frac{16}{3\sqrt{\pi}\left(1+8\pi\frac{r_{\text{s}}}{a_{\text{s}}}\rho a_{\text{s}}^{3}\right)^{2}}\sqrt{\rho a_{\text{s}}^{3}}\left[1-\frac{40}{3\sqrt{\pi}\left(1+8\pi\frac{r_{\text{s}}}{a_{\text{s}}}\rho a_{\text{s}}^{3}\right)^{2}}\sqrt{\rho a_{\text{s}}^{3}}\right]+\mathcal{O}\left[\left(\rho a_{\text{s}}^{3}\right)^{3/2}\right]\right\}.$$

(13)

In Equation (13), the effective mass $M$, emerges as a pivotal parameter characterizing weakly-interacting BECs. Our findings, encapsulated in Eq. (13) accurately reproduce previously established results in specific limiting cases. (i) In the absence of finite-range effects (i.e. $r_{\text{s}}=0$), we double check that our result aligns with Refs. [49, 18] in terms of the gas parameter $\rho a_{\text{s}}^{3}$. Firstly, in the limit where $\rho a_{\text{s}}^{3}\ll 1$, Eq. (13) decouples from the gas parameter $\rho a_{\text{s}}^{3}$, contributing to the LHY term in the universal EOS. Secondly, under the additional constraint $\rho a_{\text{s}}^{3}\ll\sqrt{\rho a_{\text{s}}^{3}}$, Eq. (13) truncates to the first order of $\sqrt{\rho a_{\text{s}}^{3}}$, contributing to beyond-LHY terms in the universal EOS. (ii) When the finite-range effects are considered (i.e. $r_{\text{s}}\neq 0$), our results accurately recover those in Ref. [30] when calculating the EOS of nonuniversal systems using the expression for $M$. This validation underscores the accuracy and applicability of our approach.

Having acquired the knowledge of $M$ in Eq. (13), we are now poised to compute EOS of the model system, leveraging the CJT effective potential as presented in Eq. (8). Specifically, the EOS explored in this paper pertains to two key quantities: the quantum depletion, denoted as $\rho_{\text{ex}}$, and chemical potential, represented by $\mu$.

Using Equation (11), the expression for quantum depletion is given by $\rho_{\text{ex}}=\frac{\sqrt{2}m^{*3/2}M^{3}}{12\pi^{2}\hbar^{3}}\left(2\sqrt{\frac{m^{*}}{m}}-\sqrt{\frac{m}{m^{*}}}\right)$. By substituting the solution for $M$ from Eq. (13) into $\rho_{\text{ex}}$, we derive the analytical expression for the quantum depletion expanded in terms of the gas parameter $\rho a_{\text{s}}^{3}$

	$\displaystyle\rho_{\text{ex}}=$		$\displaystyle\!\!\!\!\!\!\frac{8\rho}{3\sqrt{\pi}}\left(\rho a_{\text{s}}^{3}\right)^{1/2}-\frac{128\rho}{3\pi}\left(\rho a_{\text{s}}^{3}\right)+\frac{2048\rho}{9\pi^{3/2}}\left(\rho a_{\text{s}}^{3}\right)^{3/2}$		(14)
		$\displaystyle-$	$\displaystyle 64\sqrt{\pi}\rho\frac{r_{\text{s}}}{a_{\text{s}}}\left(\rho a_{\text{s}}^{3}\right)^{3/2}+\frac{5120}{3}\rho\frac{r_{\text{s}}}{a_{\text{s}}}\left(\rho a_{\text{s}}^{3}\right)^{2}.$

Next, integrating the gap equation from Eq. (10) with the thermodynamic relation for the density in Eq. (11), we systematically derive the analytical expression for the chemical potential. The expression for the chemical potential takes the form $\mu=g_{0}\rho+g_{0}\frac{\sqrt{2}m^{*3/2}M^{3}}{3\pi^{2}\hbar^{3}}\sqrt{\frac{m^{*}}{m}}$. By plugging the expression for $M$ from Eq. (13) into the formula of chemical potential, we are able to express $\mu$ as an expansion in terms of the gas parameter $\rho a_{\text{s}}^{3}$

	$\displaystyle\mu\!=\!g_{0}\rho\Big{[}1$	$\displaystyle+$	$\displaystyle\frac{32}{3\sqrt{\pi}}\left(\rho a_{\text{s}}^{3}\right)^{1/2}\!-\!\frac{512}{3\pi}\left(\rho a_{\text{s}}^{3}\right)\!+\!\frac{8192}{9\pi^{3/2}}\left(\rho a_{\text{s}}^{3}\right)^{3/2}$		(15)
		$\displaystyle-$	$\displaystyle\frac{512}{3}\sqrt{\pi}\frac{r_{\text{s}}}{a_{\text{s}}}\left(\rho a_{\text{s}}^{3}\right)^{3/2}+\frac{16384}{3}\frac{r_{\text{s}}}{a_{\text{s}}}\left(\rho a_{\text{s}}^{3}\right)^{2}\Big{]}.$

Equations (14) and (15) embody the key analytical expressions for nonuniversal EOS of a weakly-interacting Bose gas, with finite-range effects taken into consideration. Our results of Eqs. (14) and (15), demonstrate their versatility by successfully recovering previously established results under specific limiting conditions. (i) In the absence of the finite-range effects (i.e. $r_{\text{s}}=0$), our results align perfectly with the Refs. [14, 13, 49, 18], as verified through the order of the gas parameter $\rho a_{\text{s}}^{3}$. First, at the leading order $\left(\rho a_{\text{s}}^{3}\right)^{0}$, all particles condense at the mean field $\phi_{0}$, resulting in zero excess density ($\rho_{\text{ex}}=0$) and chemical potential $\mu=g_{0}\rho$. Second, proceeding to the next order $\left(\rho a_{\text{s}}^{3}\right)^{1/2}$, Eq. (15) recovers the universal LHY term [14, 13]. Third, truncating at the order of $\left(\rho a_{\text{s}}^{3}\right)^{3/2}$ for small gas parameter values, our results coincide with those reported in Ref. [18] as they should be, emphasizing the consistency and robustness of our approach. (ii) Then, when the finite-range effects are taken into account (i.e. $r_{\text{s}}\neq 0$), our results maintain their consistency with the Refs. [26, 30, 43], expanded in terms of the gas parameter $\rho a_{\text{s}}^{3}$. Notably, truncating within the order of $\left(\rho a_{\text{s}}^{3}\right)^{3/2}$, the nonuniversal EOS in Eqs. (14) and (15) match precisely with those in Refs. [26, 30, 43], obtained via the functional path-integral method. Crucially, our approach offers a significant enhancement in precision. Not only do we approximate the universal EOS in Eqs. (14) and (15) up to the order of $\left(\rho a_{\text{s}}^{3}\right)^{3/2}$, but also uncover the nonuniversal next-to-LHY EOS, given by $[16384r_{\text{s}}/3a_{\text{s}}]\left(\rho a_{\text{s}}^{3}\right)^{2}$, utilizing the CJT effective field theory. This additional insight underscores the robustness and precision of our analytical framework.

In Equations (14) and (15), the EOS comprises two distinct components: the universal terms, described by the scattering length $a_{\text{s}}$, and the nonuniversal terms, described by both the scattering length $a_{\text{s}}$ and the finite-range length $r_{\text{s}}$ based on Lagrangian density (2). However, in Refs. [17, 25, 50], the EOS has been calculated with both two-body correlations and three-body correlations taken into account, which we have mentioned in the Introduction. In particular, the universal next-to-LHY term of energy density induced by three-body correlations is $\varepsilon=\frac{g_{0}\rho^{2}}{2}\left[\frac{8\left(4\pi-3\sqrt{3}\right)}{3}\rho a_{\text{s}}^{3}\ln\left(\rho a_{\text{s}}^{3}\right)\right]$, as shown in Table 1. The relative hamiltonian density can be written as $\mathcal{H}=\frac{\hbar^{2}}{2m}\nabla\mathbf{\Phi}^{*}\cdot\nabla\mathbf{\Phi}-\mu\left|\mathbf{\Phi}\right|^{2}+\frac{U_{0}}{2}\left|\mathbf{\Phi}\right|^{4}+\frac{W}{6}\left|\mathbf{\Phi}\right|^{6}$. In the tight-bingding limit, the relationship between dimensionless $\overline{W}$ and $\overline{U}_{0}$ is $\overline{W}=\left(3\pi\right)^{-3/2}\ln\left(C\eta^{2}\right)\left(\frac{V_{0}}{E_{\text{r}}}\right)^{3/4}e^{-2\sqrt{\frac{V_{0}}{E_{\text{r}}}}}a_{\text{s}}^{2}k^{2}\overline{U}^{2}_{0}$ [51, 52]. Here, $\eta=\sqrt{\rho a_{\text{s}}^{3}}$, $C$ constant has been given in Ref. [53], $V_{0}$ represents the tunable barrier height of a homogeneous periodic lattice potential, while $k$ denotes the wave vector and $E_{\text{r}}=\hbar^{2}k^{2}/2m$ is the recoil energy. In contemporary experiments, $a_{\text{s}}^{2}k^{2}$ typically spans from $10^{-8}$ to $10^{-2}$ [54]. Within this range, the influence of three-body interactions is markedly negligible in comparison to two-body interactions, implying that the effects stemming from three-body interactions are exceedingly difficult to discern experimentally. Therefore, we neglect the three-body correlations when we calculate the EOS for our boson-boson interacting system.

Meanwhile, the known next-to-LHY term of EOS induced by three-body correlations resides at the order of $\left(\rho a_{\text{s}}^{3}\right)^{1}$, aligning with the first term of universal next-to-LHY EOS though CJT effective field theory, as demonstrated in the third line in Table 1. It can be roughly interpreted that the CJT effective field theory essentially effects a higher-order perturbation of the two-body interacting system in comparison to a direct consideration of three-body interactions. We note that the EOS in Eq. (15) is at the order of $\left(\rho a_{\text{s}}^{3}\right)^{2}$, which is higher than the results in Refs. [25, 50]. In this sense, we need to consider the three-body correlations. Meanwhile, there is no need to use CJT effective field theory when calculating the nonuniversal EOS induced by three-body effect, which is supposed to give the higher order than $\left(\rho a_{\text{s}}^{3}\right)^{2}$. It is enough to obtain a proper EOS through adding the terms induced by three-body correlations in Refs. [25, 50] into our results in Eq. (15) by hand. We conclude our result together with the logarithmic term of the next-to-LHY correction to the universal EOS gives a reasonable description of the EOS of an ultracold Bose gas.

Next, based on the chemical potential $\mu$, the inverse compressibility $\kappa^{-1}=\rho\frac{\partial\mu}{\partial\rho}$ can be readily derived

	$\displaystyle\kappa^{-1}=$		$\displaystyle\!\!\!\!\!\!g_{0}\rho\Big{[}1+\frac{3}{2}\alpha\rho^{1/2}-3\alpha^{2}\rho+\frac{15}{8}\alpha^{3}\rho^{3/2}$		(16)
		$\displaystyle-$	$\displaystyle\frac{45\pi^{2}}{128}\frac{r_{\text{s}}}{a_{\text{s}}}\alpha^{3}\rho^{3/2}+\frac{81\pi^{2}}{64}\frac{r_{\text{s}}}{a_{\text{s}}}\alpha^{4}\rho^{2}\Big{]},$

with $\alpha=\left(32/3\sqrt{\pi}\right)a_{\text{s}}^{3/2}$.

To delve deeper into our analysis, utilizing the Eq. (16), we explore the nonuniversal quantum effect present in the speed of sound, denoted as $c_{\text{s}}=\sqrt{\left(\kappa m\right)^{-1}}$,

	$\displaystyle c_{\text{s}}=\!\!\sqrt{\frac{g_{0}\rho}{m}}\Bigg{\{}$	$\displaystyle 1$	$\displaystyle+\frac{8\sqrt{\rho a_{\text{s}}^{3}}}{\sqrt{\pi}}\Big{[}1-\frac{64\sqrt{\rho a_{\text{s}}^{3}}}{3\sqrt{\pi}}+\frac{1280}{9\pi}\rho a_{\text{s}}^{3}$		(17)
		$\displaystyle-$	$\displaystyle\frac{80\pi}{3}\frac{r_{\text{s}}}{a_{\text{s}}}\rho a_{\text{s}}^{3}+1024\sqrt{\pi}\frac{r_{\text{s}}}{a_{\text{s}}}\sqrt{\rho a_{\text{s}}^{3}}^{3}\Big{]}\Bigg{\}}.$

In Equation (17), the coefficient associated with the term $\sqrt{\rho a_{\text{s}}^{3}}$ is $8/\sqrt{\pi}$, differing from LHY’s prediction of $16/\sqrt{\pi}$. The deviation is readily traceable to the corresponding corrective terms present in the definition of the inverse compressibility $\kappa^{-1}$, given by Eq. (16). As a consequence of the nonuniversal effect, our derived expression for the sound speed $c_{\text{s}}$ incorporates higher-order terms in the gas parameter $\rho a_{\text{s}}^{3}$, distinguishing it from the findings presented in Ref. [49]. This feature underscores the paramount importance and significance of our approach.

In Sec. III, we have rigorously derived the nonuniversal next-to-LHY term in EOS (see Eqs. (14) and (15)) to the order of $\left(\rho a_{\text{s}}^{3}\right)^{2}$. The purpose of Sec. IV is to propose an experimental protocol to observe the nonuniversal beyond-LHY corrections to the EOS by calculating the frequency shifts in the breathing modes. In contemporary experiments, the gas parameter $\rho a_{\text{s}}^{3}$ typically remains below $10^{-4}$, rendering even the initial LHY correction to the mean-field energy, scaling as $\sqrt{\rho a_{\text{s}}^{3}}$, negligible at best, contributing less than $1\%$ of the total energy. Consequently, the higher-order corrections beyond LHY in Eq. (15) are anticipated to be too subtle to discern in density profiles or release energies. While enhancing quantum fluctuations by tuning $a_{\text{s}}$ via Feshbach resonance [31, 32, 55, 33, 35, 34] can amplify their impact, we propose an alternative route: studying the frequency shifts in collective excitations [56, 57, 58, 59] induced by quantum fluctuations [60]. To this end, we augment our model system with a 3D harmonic trap, defined by $V_{\text{ext}}=\frac{1}{2}m\omega^{2}_{0}r^{2}$, and embark on calculating the resulting shifts in the breathing mode frequency, attributed to beyond-LHY terms in Eq. (15). This endeavor is motivated by the potential to experimentally verify and quantify the nonuniversal contributions to the EOS.

Adhering closely to the standard methodology outlined in Ref. [38] for calculating frequency shifts, we formulate the hydrodynamic equation as follows:

$$m\frac{\partial^{2}\delta\rho\left(\mathbf{r},t\right)}{\partial t^{2}}-\nabla\cdot\left[\rho\left(\mathbf{r}\right)\nabla\left(\frac{\partial\mu_{l}}{\partial\rho}\delta\rho\left(\mathbf{r},t\right)\right)\right]=0.$$

(18)

Equation (18) incorporates the density fluctuation, denoted as $\delta\rho\left(\mathbf{r},t\right)$, which arises around the targeted ground state density $\rho\left(\mathbf{r}\right)$ and the local chemical potential $\mu_{l}$ as specified in Eq. (15).

As an initial step, utilizing Eq. (15), we can iteratively derive the expansion of the ground state density, expressed in terms of $\alpha=\left(32/3\sqrt{\pi}\right)a_{\text{s}}^{3/2}$,

	$\displaystyle\rho\left(\mathbf{r}\right)=$		$\displaystyle\!\!\!\!\!\!\rho_{\text{TF}}-\alpha\rho_{\text{TF}}^{3/2}+\frac{3}{2}\alpha^{2}\rho_{\text{TF}}^{2}-\frac{3}{4}\alpha^{3}\rho_{\text{TF}}^{5/2}$		(19)
		$\displaystyle+$	$\displaystyle\frac{9}{64}\pi^{2}\frac{r_{\text{s}}}{a_{\text{s}}}\alpha^{3}\rho_{\text{TF}}^{5/2}-\frac{27}{64}\pi^{2}\frac{r_{\text{s}}}{a_{\text{s}}}\alpha^{4}\rho_{\text{TF}}^{3}.$

In Equation (19), $\rho_{\text{TF}}=\left(g_{0}\rho-V_{\text{ext}}\right)/g_{0}$ [37, 38], is the so-called Thomas-Fermi result for the ground state density. Examining Eq. (19), we observe that its first line comprises three distinct terms: the mean-field term being $\rho_{\text{TF}}$, the universal LHY term represented by $-\alpha\rho_{\text{TF}}^{3/2}$, and the universal beyond-LHY term $3/2\alpha^{2}\rho_{\text{TF}}^{2}-3/4\alpha^{3}\rho_{\text{TF}}^{5/2}$. Notably, the second line of Eq. (19) vanishes solely when nonuniversal effects are disregarded, emphasizing their significance. Consequently, $\frac{9}{64}\pi^{2}\frac{r_{\text{s}}}{a_{\text{s}}}\alpha^{3}\rho_{\text{TF}}^{5/2}$ is designated as the nonuniversal LHY term, while $-\frac{27}{64}\pi^{2}\frac{r_{\text{s}}}{a_{\text{s}}}\alpha^{4}\rho_{\text{TF}}^{3}$ corresponds to the nonuniversal next-to-LHY term.

Subsequently, we derive the expansion of $\left[\rho\partial\mu_{l}/\partial\rho\right]$ as a series in terms of $\rho_{\text{eff}}$ by substituting Eq. (19) into Eq. (16). Following this substitution, we insert both Eqs. (16) and (19) into the Eq. (18). Through meticulous algebraic manipulations, we ultimately arrive at

	$\displaystyle\!\!\!\!\!\!\!\!m\omega^{2}\delta\rho+\nabla\cdot\left[g_{0}\rho_{\text{TF}}\nabla\delta\rho\right]=$	$\displaystyle-$	$\displaystyle\frac{1}{2}\alpha\nabla^{2}\left(g_{0}\rho_{\text{TF}}^{3/2}\delta\rho\right)+\frac{15}{4}\alpha^{2}\nabla^{2}\left(g_{0}\rho_{\text{TF}}^{2}\delta\rho\right)-3\alpha^{2}\nabla\cdot\left(g_{0}\rho_{\text{TF}}\delta\rho\nabla\rho_{\text{TF}}\right)$		(20)
		$\displaystyle+$	$\displaystyle\left(-\frac{45}{4}+\frac{63\pi^{2}}{128}\frac{r_{\text{s}}}{a_{\text{s}}}\right)\alpha^{3}\nabla^{2}\left(g_{0}\rho_{\text{TF}}^{5/2}\delta\rho\right)+\frac{57}{4}\alpha^{3}\nabla\cdot\left(g_{0}\rho_{\text{TF}}^{3/2}\delta\rho\nabla\rho_{\text{TF}}\right)$
		$\displaystyle+$	$\displaystyle\left(\frac{165}{8}-\frac{261\pi^{2}}{128}\frac{r_{\text{s}}}{a_{\text{s}}}\right)\alpha^{4}\nabla^{2}\left(g_{0}\rho_{\text{TF}}^{3}\delta\rho\right)+\left(-\frac{495}{16}+\frac{27\pi^{2}}{16}\frac{r_{\text{s}}}{a_{\text{s}}}\right)\alpha^{4}\nabla\cdot\left(g_{0}\rho_{\text{TF}}^{2}\delta\rho\nabla\rho_{\text{TF}}\right).$

Equation (20) simplifies into the basic hydrodynamic equation $m\omega^{2}\delta\rho+\nabla\cdot\left[g_{0}\rho_{\text{TF}}\nabla\delta\rho\right]=0$ in case of $\alpha=0$. For this simplified scenario, the calculated frequency $\omega$ exhibits analytical solutions of the form $\omega\left(n_{r},l\right)=\omega_{0}\left(2n_{r}^{2}+2n_{r}l+3n_{r}+l\right)^{1/2}$, with $n_{r}$ representing the number of radial nodes and $l$ denoting the angular momentum associated with the excitation. This analytical expression provides a direct link between the frequency of the excitation and its quantum numbers.

Finally, Equation (20) can be routinely tackled by considering its right-hand side as a minor perturbation. By adopting this approach, we can derive the analytical expression for the frequency shifts,

	$\displaystyle\frac{\delta\omega}{\omega}=$	$\displaystyle-$	$\displaystyle\frac{\alpha g_{0}}{4m\omega^{2}}\frac{\int d^{3}{\mathbf{r}}\nabla^{2}\left(\delta\rho^{*}\right)\rho_{\text{TF}}^{3/2}\delta\rho}{\int d^{3}{\mathbf{r}}\delta\rho^{*}\delta\rho}+\frac{15\alpha^{2}g_{0}}{8m\omega^{2}}\frac{\int d^{3}{\mathbf{r}}\nabla^{2}\left(\delta\rho^{*}\right)\rho_{\text{TF}}^{2}\delta\rho}{\int d^{3}{\mathbf{r}}\delta\rho^{*}\delta\rho}+\frac{3\alpha^{2}g_{0}}{2m\omega^{2}}\frac{\int d^{3}{\mathbf{r}}\nabla\left(\delta\rho^{*}\right)\cdot\nabla\left(\rho_{\text{TF}}\right)\rho_{\text{TF}}\delta\rho}{\int d^{3}{\mathbf{r}}\delta\rho^{*}\delta\rho}$		(21)
		$\displaystyle+$	$\displaystyle\frac{\left(-\frac{45}{4}+\frac{63\pi^{2}}{128}\frac{r_{\text{s}}}{a_{\text{s}}}\right)\alpha^{3}g_{0}}{2m\omega^{2}}\frac{\int d^{3}{\mathbf{r}}\nabla^{2}\left(\delta\rho^{*}\right)\rho_{\text{TF}}^{5/2}\delta\rho}{\int d^{3}{\mathbf{r}}\delta\rho^{*}\delta\rho}-\frac{\frac{57}{4}\alpha^{3}g_{0}}{2m\omega^{2}}\frac{\int d^{3}{\mathbf{r}}\nabla\left(\delta\rho^{*}\right)\cdot\nabla\left(\rho_{\text{TF}}\right)\rho_{\text{TF}}^{3/2}\delta\rho}{\int d^{3}{\mathbf{r}}\delta\rho^{*}\delta\rho}$
		$\displaystyle+$	$\displaystyle\frac{\left(\frac{165}{8}-\frac{261\pi^{2}}{128}\frac{r_{\text{s}}}{a_{\text{s}}}\right)\alpha^{4}g_{0}}{2m\omega^{2}}\frac{\int d^{3}{\mathbf{r}}\nabla^{2}\left(\delta\rho^{*}\right)\rho_{\text{TF}}^{3}\delta\rho}{\int d^{3}{\mathbf{r}}\delta\rho^{*}\delta\rho}-\frac{\left(-\frac{495}{16}+\frac{27\pi^{2}}{16}\frac{r_{\text{s}}}{a_{\text{s}}}\right)\alpha^{4}g_{0}}{2m\omega^{2}}\frac{\int d^{3}{\mathbf{r}}\nabla\left(\delta\rho^{*}\right)\cdot\nabla\left(\rho_{\text{TF}}\right)\rho_{\text{TF}}^{2}\delta\rho}{\int d^{3}{\mathbf{r}}\delta\rho^{*}\delta\rho}.$

The integrals in Eq. (21) are confined to the domain where the Thomas-Fermi density remains positive. To delve into effects that transcend the LHY theory, it becomes imperative to concentrate on compressional modes, which are highly responsive to alterations in the EOS. Among these, the breathing mode in a spherical trap, characterized by $\left(n_{r}=1,l=0\right)$, stands as the fundamental excitation. This mode is distinguished by its zeroth-order dispersion, yielding a frequency of $\omega=\sqrt{5}\omega_{0}$ and exhibits density oscillations that adhere to the pattern $\delta\rho\sim\left(r^{2}-3/5R^{2}\right)$. In this specific scenario, Eq. (21) furnishes:

	$\displaystyle\frac{\delta\omega}{\omega}=$		$\displaystyle\!\!\!\!\!\!\frac{63\sqrt{\pi}}{128}\left[\rho\left(0\right)a_{\text{s}}^{3}\right]^{1/2}-\frac{16384}{135\pi}\left[\rho\left(0\right)a_{\text{s}}^{3}\right]$		(22)
		$\displaystyle+$	$\displaystyle\frac{3682-\frac{2205}{16}\pi^{2}\frac{r_{\text{s}}}{a_{\text{s}}}}{3\sqrt{\pi}}\left[\rho\left(0\right)a_{\text{s}}^{3}\right]^{3/2}$
		$\displaystyle+$	$\displaystyle\left(\frac{-71200}{\pi^{2}}+6465\frac{r_{\text{s}}}{a_{\text{s}}}\right)\left[\rho\left(0\right)a_{\text{s}}^{3}\right]^{2},$

revealing the fractional shift in the breathing mode frequency, where $\rho\left(0\right)$ represents the density evaluated at the center of the trap.

Equation (22) constitutes another pivotal finding, encapsulating corrections to the breathing mode frequency that scale with the gas parameter $\rho\left(0\right)a_{\text{s}}^{3}$. The corrections in Eq. (22) stem from diverse origins, encompassing both LHY and beyond-LHY contributions to the EOS. The first term on the right-hand side of Eq. (22) represents the ubiquitous LHY-mediated fractional shift of the breathing mode frequency, which scales as $\left[\rho\left(0\right)a_{\text{s}}^{3}\right]^{1/2}$. All the subsequent terms in Eq. (22) are of higher order than $\left[\rho\left(0\right)a_{\text{s}}^{3}\right]^{1/2}$, offering insights into more intricate effects. Notably, among these beyond-LHY terms, the one scaling with $\left[\rho\left(0\right)a_{\text{s}}^{3}\right]$ is a so-called next-to-LHY universal fractional shift. This term, specifically $\frac{-16384}{135\pi}\left[\rho\left(0\right)a_{\text{s}}^{3}\right]$, underscores the persistence of universal behavior beyond the leading LHY correction. Particularly, terms in Eq. (22) involving $r_{\text{s}}$ represent nonuniversal contributions, arising from the finite-range effective potential, further illuminating the intricacies of the system’s response. Thus, the complexity deepens with the term $\frac{3682-\frac{2205}{16}\pi^{2}\frac{r_{\text{s}}}{a_{\text{s}}}}{3\sqrt{\pi}}\left[\rho\left(0\right)a_{\text{s}}^{3}\right]^{3/2}$ in Eq. (22), which encapsulates both the next-next-to-LHY universal fractional shift $\frac{3682}{3\sqrt{\pi}}\left[\rho\left(0\right)a_{\text{s}}^{3}\right]^{3/2}$ and a LHY nonuniversal contribution $\frac{-\frac{2205}{16}\pi^{2}\frac{r_{\text{s}}}{a_{\text{s}}}}{3\sqrt{\pi}}\left[\rho\left(0\right)a_{\text{s}}^{3}\right]^{3/2}$. Similarly, the last term $\left(\frac{-71200}{\pi^{2}}+6465\frac{r_{\text{s}}}{a_{\text{s}}}\right)\left[\rho\left(0\right)a_{\text{s}}^{3}\right]^{2}$ in Eq. (22) comprises two distinct components: the next-next-next-to-LHY universal fractional shift $\frac{-71200}{\pi^{2}}\left[\rho\left(0\right)a_{\text{s}}^{3}\right]^{2}$ and the next-to-LHY nonuniversal fractional shift $6465\frac{r_{\text{s}}}{a_{\text{s}}}\left[\rho\left(0\right)a_{\text{s}}^{3}\right]^{2}$. We remark that to our best knowledge, the analytical expressions of beyond-LHY-induced fractional shift of the breathing mode frequency in Eq. (22) are obtained for the first time.