and subsequent Chowla-Selberg lattice summation technique [@Milton2010]. See– Appendix-\[A\] for the use of Chowla-Selberg lattice summation in this subsection too. Following this technique we get $N^{(4)}=0.00914223/R^4$, $D^{(4)}=-0.000263472/R^4$, $N^{(6)}=-0.0185855/R^6$, $D^{(6)}=-0.000389195/R^6$, $N^{(8)}=0.0745135/R^8$, $D^{(8)}=0.00400674/R^8$, and $N^{(1)}=D^{(1)}=N^{(3)}=N^{(D)}=N^{(5)}=D^{(5)}=0$. Substituting these values of $N^{(-2s)}$ and $D^{(-2s)}$ in Eq. , we have $$\begin{aligned}
\label{eqn:70}
{\mathcal E}_0(R,\bar{n})&=&-\frac{0.013561343~\hbar v(\bar{n})L}{R^2}\bigg(1\nonumber\\
&-&\frac{0.008878758}{0.013561343\times4\tilde{\rho}(R,g\bar{n})} \\
&-&\frac{0.018974695}{0.013561343\times32\tilde{\rho}^2(R,g\bar{n})}+...\bigg),~ \nonumber \end{aligned}$$ where $\tilde{\rho}(R,g\bar{n})=\frac{g\bar{n}2MR^2}{\hbar^2}$ is the dimensionless density. The quantum (vacuum) fluctuation part of the Casimir force, $f_R^{\text{(qf)}}=-\frac{\partial}{\partial R}{\mathcal E}_0(R,\bar{n})$, is now easily obtained from Eq. , as $$\begin{aligned}
\label{eqn:71}
f_R^{(\text{qf})}=-0.019\frac{\hbar^2L}{MR^4}\tilde{\rho}^{1/2}\bigg(1-\frac{0.655}{\tilde{\rho}}-\frac{0.262}{\tilde{\rho}^2}+...\bigg).~~\end{aligned}$$ Note that the first term of the Eq.  is the most dominating term for the quantum phonon fluctuations part of the Casimir force, and it is identical to Eq. , which is just based on an intuitive similarity to the case of Casimir force arising the TE modes of the vacuum fluctuations of the electromagnetic field.

Mean field part of the Casimir force
------------------------------------

In this subsection we calculate the contribution to the Casimir force arising from the first term in Eq. , [*i.e. *]{} the mean field part ($\Omega_0$) of the grand canonical Hamiltonian (${\mathcal H}$). Mean-field part of the Casimir force arises due to inhomogeneity of the BEC order parameter \[$m(\textbf{r})=\sqrt{N_0}\phi_0(\textbf{r})$\] in $\Omega_0$. If 