next up previous contents
Next: Quantum Monte Carlo Method Up: Superfluid density Previous: Superfluid fraction in the   Contents


calculation of the superfluid fraction from GPE

In this section the superfluid fraction $\rho _s/\rho $ will be obtained directly from Gross-Pitaevskii equation in a perturbative manner. This derivation is new and the result coincides with the one obtained from the Bogoliubov model presented in the previous section. The Gross-Pitaevskii equation (1.11) for the condensate wavefunction in the absence of external field takes the form
$\displaystyle i\hbar\frac{\partial }{\partial t} \Psi({\bf r}, t) =
\left(-\frac{\hbar^2}{2m} \triangle + g\vert\Psi({\bf r}, t)\vert^2\right)
\Psi({\bf r}, t)$     (2.46)

Let us add a moving impurity that creates an external field $V_{ext}({\bf r}, t) = g_{imp} \delta({\bf r - V}t)$ and let us treat it as a perturbation to the solution $\Psi_0({\bf r}, t)$ of the equation (2.46), i.e.
$\displaystyle \Psi({\bf r}, t) = \left[\Psi_0({\bf r}) + \delta\Psi({\bf r}, t)\right]
\exp\left(-i \frac{\mu t}{\hbar} \right),$     (2.47)

Then, substitution of (2.47) into (2.46) gives
$\displaystyle i\hbar\frac{\partial }{\partial t} \delta \Psi =
\left(-\frac{\hb...
...elta \Psi
+ g \Psi_0^2 \delta \Psi^\star
+ g_{imp} \delta ({\bf r - V}t) \Psi_0$     (2.48)

The perturbation follows the moving impurity, so $\delta \Psi$ is a function of ${\bf r - V}t$. Let us introduce the new variable ${\bf r' = r - V}t$. It means that the coordinate derivative can be related to time derivative
$\displaystyle \frac{\partial}{\partial t} \delta \Psi({\bf r - V}t)$ $\textstyle =$ $\displaystyle -{\bf V \nabla} \delta \Psi({\bf r - V}t),$ (2.49)
$\displaystyle \triangle {\bf r'}$ $\textstyle =$ $\displaystyle \triangle {\bf r}$ (2.50)

and instead of (2.48) one has
$\displaystyle \left( i\hbar {\bf V \nabla} - \frac{\hbar^2}{2m} \triangle
- \mu...
...
+ g \Psi_0({\bf r})^2 \delta \Psi^\star
+ g_{imp} \delta ({\bf r}) \Psi_0 = 0,$     (2.51)

(from now on, the subscript over ${\bf r}$ will be dropped) Taking the Fourier transform of this equation and treating $\Psi_0$ as a real constant (i.e. the solution for the uniform case $g\Psi_0^2 = \mu$) one obtains
$\displaystyle \left( -\hbar {\bf k V} + \frac{\hbar^2k^2}{2m}
- \mu + 2 \mu \right) \delta \Psi_k
+ \mu\delta (\Psi_{-k})^\star + g_{imp} \Psi_0 = 0,$     (2.52)

where we used the property of Fourier components $\delta (\Psi^\star)_{k} = \delta (\Psi_{-k})^\star$. The substitution of $(-{\bf k})$ in the equation complex conjugate of (2.52) gives
$\displaystyle \left( \hbar {\bf k V} + \frac{\hbar^2k^2}{2m}
- \mu + 2 \mu \right) (\delta \Psi_{-k})^\star
+ \mu\delta \Psi_{k} + g_{imp} \Psi_0 = 0$     (2.53)

The solution for the system of linear equations (2.52 - 2.53) is given by
$\displaystyle \delta \Psi_{k}
= -\frac{g_{imp} \left(\hbar{\bf k V} + \frac{\hb...
...{\hbar^2 k^2}{2m} \left(\frac{\hbar^2 k^2}{2m}+2\mu\right)
-(\hbar{\bf k V})^2}$     (2.54)

The energy $E' = E-\mu N$ has the minimum at fixed $\mu$ for the ground state function $\Psi_0$. It means that $E'$ does not have terms linear in $\delta \Psi_k$ and $\delta \Psi_k^\star$, so $E'= E^{(0)}+E^{(2)}+g_{imp} (\Psi_0^\star\delta\Psi(0)+\Psi_0\delta\Psi^\star(0))$. Here the last term comes from the linear expansion of the energy $\int \vert\Psi({\bf r})\vert^2 g_{imp} \delta({\bf r}){\bf dr}$. The term $E^{(2)}$ being quadratic in $\delta \Psi_k$ and $\delta \Psi_k^\star$ satisfies the Euler identity:
$\displaystyle 2 E^{(2)} =
\int \left[
\delta \Psi({\bf r}) \frac{\delta E^{(2)}...
...f r}) \frac{\delta E^{(2)}}{\delta(\delta
\Psi^\star({\bf r}))}
\right]{\bf dr}$     (2.55)

which using the variational equation
$\displaystyle i \hbar \frac{\partial\delta(\delta \Psi)}{\partial t} =
\frac{\delta E^{(2)}}{\delta(\delta\Psi^\star)}
+ g_{imp}\Psi_0\delta({\bf r})$     (2.56)

can be rewritten as
$\displaystyle E^{(2)}= E^{(2)}_1 + E^{(2)}_2=
\frac{i\hbar}{2} \int \left[
\del...
...f dr}-
\frac{g_{imp} }{2} (\Psi_0^\star\delta\Psi(0)+\Psi_0\delta\Psi^\star(0))$     (2.57)

To start with, let us Fourier transform the first term. Exchanging time derivatives with gradients by the rule (2.49) one obtains
$\displaystyle E^{(2)}_1 = \int \hbar{\bf k V} \vert\delta \Psi_k\vert^2 \frac{\bf dk}{(2\pi)^3}$     (2.58)

The terms of interest are the ones that are quadratic in the velocity ${\bf V}$. It means that the term $(\hbar{\bf kV})^2$ in the denominator of (2.54) can be neglected and $\vert\delta
\Psi_k\vert^2$ turns out to be
$\displaystyle \vert\delta \Psi_{k}\vert^2 = \frac{g_{imp} ^2 \vert\Psi_0\vert^2...
...eft[\frac{\hbar^2 k^2}{2m}
\left( \frac{\hbar^2 k^2}{2m}+2 \mu\right)\right]^2}$     (2.59)

The energy does not have terms linear in ${\bf V}$, because all terms independent of ${\bf V}$ in (2.59) are even in ${\bf k}$, so multiplied by ${\bf k}$ and integrated over momentum space they provide zero contribution to the energy. The only term that is left is the following
$\displaystyle E^{(2)}_1 =
2g_{imp} ^2\vert\Psi_0\vert^2
\int \frac{(\hbar{\bf k...
... k^2}{2m}
\left( \frac{\hbar^2 k^2}{2m}+2 \mu\right)^2}
\frac{\bf dk}{(2\pi)^3}$     (2.60)

For the calculation of $\delta\Psi(0)$ one should consider $\delta \Psi_k$ taking into account that $\hbar {\bf kV} \ll \mu$ and then integrate it over momentum space
$\displaystyle \delta \Psi_{k}
\approx
\frac{(\hbar{\bf kV})^2}{\frac{\hbar^2 k^...
...right)\right]^2}
+\frac{1}{\frac{\hbar^2 k^2}{2m}+2\mu}
\right\}
g_{imp} \Psi_0$     (2.61)

The integral of the second term over momentum space is equal to zero. The third term is diverging and needs the renormalization of $g_{imp}$ (as discussed in sections 1.3.1 and 2.2.2) in order to be calculated correctly. However, its correction does not depend on ${\bf V}$ and will be omitted. The energy is defined by the following integral
$\displaystyle \delta E = E^{(2)}_1 + E^{(2)}_2 +
g_{imp} (\Psi_0^\star\delta\Ps...
...^2}{2m}
\left( \frac{\hbar^2 k^2}{2m}+2 \mu\right)^2}
\frac{d{\bf k}}{(2\pi)^3}$     (2.62)

In the integral $(\bf {kV})^2 d{\bf k}$ can be replaced by $1/3~k^2
V^2 4\pi k^2 dk$ due to the equivalence of different directions. Then the integral can be easily calculated if one recall the following integral identity
$\displaystyle \int \frac{x^2dx}{(x^2+a^2)^2} = -\frac{x}{2(x^2+a^2)^2}
+\frac{1}{2a} \mathop{\rm atan}\nolimits \frac{x}{a}$     (2.63)

The result is the following
$\displaystyle \delta E = \frac{m^{5/2}g_{imp} ^2\vert\Psi_0\vert^2}{3\pi\hbar^3\sqrt{\mu}}\frac{V^2}{2},$     (2.64)

where
$\displaystyle g = \frac{4\pi\hbar^2a}{m},\quad
g_{imp} = \frac{2\pi\hbar^2b}{m},\quad
\vert\Psi_0\vert^2 = n$      

The term in front of $V^2/2$ in (2.64) can be interpreted as an effective mass $m^\star$ of the particles which follow the external perturbation, i.e. the normal (and not superfluid) component of the fluid. Then normal fraction can then be easily obtained which, as anticipated, coincides with result (2.45).
$\displaystyle \frac{\rho_n}{\rho} = \frac{m^\star}{m} \chi =
\frac{2\sqrt{\pi}}{3} (na^3)^{1/2} \chi\left(\frac{b}{a}\right)^2$     (2.65)

Let us compare the results for the superfluid density (2.65) and the condensate fraction (2.22). It is interesting to note that in both formulae the effect of disorder enters as $\sqrt{na^3} R$, where $R = \chi (b/a)^2$ is the universal scaling parameter which already entered the result (2.17) for the energy. This means that systems with different disorder concentration $\chi $ and size of the impurities $b/a$, but same $R$ experience the same effect due to disorder. Another interesting result is that disorder is more efficient (by a factor $4/3$) in depleting the superfluid density than the condensate. Taking into account that even pure systems ($R = 0$) exhibit a nonzero quantum depletion due to particle interactions, one infers that at the critical amount of disorder $R_c =
16/\pi\approx 5.1$ the depletion of the superfluid density becomes larger than the depletion of the condensate fraction. Huang and Meng [16], who first derived results (2.22) and (2.45) even if for a different model of disorder, have used them at $T = 0$ to predict two distinct transitions as a function of the amount of disorder : first a superfluid-insulator transition where $\rho_s = 0$ followed by a Bose-Einstein transition where $N_0 = 0$. These authors also argue that the intermediate phase corresponding to $\rho_s = 0$ and $N_0 \ne 0$ should be identified with a Bose-glass phase. However, in [19] it is stressed that results ([*]) and (2.45) are valid in the weak disorder regime and cannot be applied if the depletion due to disorder is large. The range of validity of results (2.22) and (2.65) will be investigated in detail in section 4.4.2 using Monte Carlo techniques.
next up previous contents
Next: Quantum Monte Carlo Method Up: Superfluid density Previous: Superfluid fraction in the   Contents