NSR Calculate

1. Free Energy vs. Renormalize Paramaters
- 1.1. p-wave
  - 1.1.1. cal
- 1.2. 以 \(\varepsilon\) 为单位
2. result
3. code

1. Free Energy vs. Renormalize Paramaters

1.1. p-wave

1.1.1. cal

\begin{align} \delta^p(\vec{q}, z) =& \mathrm{Arg}\left[ \frac{Mk_{n^2}}{2}\frac{1}{R} \left( \frac{1}{4\pi}\cdot \frac{2R}{k_n^2 v} + \tilde{z}\cdot\frac{1}{4\pi} + \frac{2 R}{M k_n^2}\Pi_r(\vec{q},z) \right) \right] \\ =& \mathrm{Arg}\left[ \frac{1}{4\pi}\cdot \frac{2R}{k_n^2 v} + \tilde{z}\cdot\frac{1}{4\pi} + \frac{2 R}{M k_n^2}\Pi_r(\vec{q},z + \mathrm{i}0^+) \right] \end{align}

其中 \(\tilde{z}=z/E_n\) , \(E_n = k_n^2/(2M)\) , \(k_n^3 = 6\pi^2n\) , \(n = N/V\)

\begin{align} \frac{2 R}{M k_n^2}\Pi_r(\vec{q},z) =& (k_n R)\cdot\Pi_r \cdot \frac{2}{Mk_n^3}\\ =& \tilde{R}\left[ \frac{2}{Mk_n^3}\left( -\frac{M}{V} \right)\sum_{\vec{k}}1 - \tilde{z}E_n \frac{M^2}{V}\frac{2}{Mk_n^3}\sum_{\vec{k}}\frac{1}{k^2} + \frac{2}{Mk_n^3}\Pi^{l=1}(\vec{q},z) \right] \\ =& \tilde{R}\left[ -\frac{1}{\pi^2}\int \mathrm{d}\tilde{k}\cdot \tilde{k}^2 -\tilde{z} \frac{1}{2\pi^2}\int \mathrm{d}\tilde{k} +\tilde{\Pi}^{l=1} \right] \end{align}

其中 \(\tilde{R} = k_nR\) , \(\tilde{k} = k/k_n\)

\begin{align} \tilde{\Pi}^{l=1} = &\frac{2}{Mk_n^3}\Pi^{l=1}(\vec{q},\omega) \\ =& \frac{2}{Mk_n^3}\frac{1}{V}\frac{V}{(2\pi)^3}\int \mathrm{d}\tilde{k} \left[ k^2 \cdot 4\pi |Y_{lm}(\hat{k})|^2 \frac{1+n(\xi_{\vec{k}+\vec{q}/2}) + n(\xi_{-\vec{k}+\vec{q}/2})} {\xi_{\vec{k}+\vec{q}/2} + \xi_{-\vec{k}+\vec{q}/2} - \omega} \right] \\ =& \frac{2}{\pi^2}\int \mathrm{d}\tilde{k}\cdot\tilde{k}^4\left[ \frac{1+n(\xi_{\vec{k}+\vec{q}/2}) + n(\xi_{-\vec{k}+\vec{q}/2})} {\tilde{\xi}_{\vec{k}+\vec{q}/2} + \tilde{\xi}_{-\vec{k}+\vec{q}/2} - \tilde{\omega}} \right] \end{align}

其中 \(\tilde{\xi} = \xi/E_n\) , \(\tilde{\omega} = \omega/E_n\) , \(n(\xi) = \frac{1}{e^{\beta \xi}-1}\)

最终

\begin{align} \frac{\tilde{\Omega}}{N E_n} =& \frac{1}{N E_n} \frac{V}{(2\pi^3)}\int \mathrm{d}^3\vec{q} \cdot \int \frac{\mathrm{d}\omega}{\pi}\cdot \frac{1}{e^{\beta\omega}-1} \delta^p \\ =& \frac{3}{\pi} \int \mathrm{d}\tilde{q}\cdot \tilde{q}^2 \int_{-\infty}^{+\infty}\mathrm{d}\tilde{\omega} \cdot \frac{1}{e^{\tilde{\beta}\tilde{\omega}}-1} \tilde{\delta}^p(\vec{q},z) \end{align}

其中 \(\tilde{\beta} = \beta E_n\) . 得自由能

\begin{align} \frac{F}{NE_n} = \frac{\tilde{\Omega}}{N E_n} -\frac{\mu}{E_n} \end{align}

即

\begin{align} f(\tilde{\mu}, \tilde{R}) = \tilde{\Omega}'(\tilde{\mu}, \tilde{R})-\tilde{\mu} \end{align}

其中 \(\tilde{\mu} = \mu/E_n\).

\(\mu\) 由

\begin{align} N = - \frac{\partial\Omega}{\partial \mu} \end{align}

决定.

1.2. 以 \(\varepsilon\) 为单位

若以某一能量 \(\varepsilon\) 为单位, 对应的长度单位 \(k_{\varepsilon} = \sqrt{2M\varepsilon}\) , 密度单位 \(n_{\varepsilon} = k_{\varepsilon}^3/(6\pi^2)\) , 那么

\begin{align} \frac{\Omega}{N \varepsilon} = & \frac{n_{\varepsilon}}{n}\int \mathrm{d}\tilde{q}\cdot \tilde{q}^2 \int_{-\infty}^{+\infty}\mathrm{d}\tilde{\omega} \cdot\frac{3}{\pi}\cdot \frac{1}{e^{\tilde{\beta}\tilde{\omega}}-1} \tilde{\delta}^p(\vec{q},z) \\ = & \frac{n_{\varepsilon}}{n}\int \mathrm{d}\tilde{q}\cdot \tilde{q}^2 \int_{-\infty}^{+\infty}\mathrm{d}\tilde{\omega} \cdot f(\tilde{q}, \tilde{\omega}, \tilde{\mu}, \tilde{\beta}) \end{align}

其中

\begin{align} f(\tilde{q}, \tilde{\omega}, \tilde{\mu}, \tilde{\beta}) = \frac{3}{\pi}\cdot \frac{1}{e^{\tilde{\beta}\tilde{\omega}}-1} \tilde{\delta}^p(\vec{q},z) \end{align}

而

\begin{align} \frac{n}{n_{\varepsilon}} =& - \frac{1}{n_{\varepsilon}V} \frac{\partial\Omega}{\partial\mu} =- \frac{1}{n_{\varepsilon}V} \frac{\partial\Omega/\mu}{\partial\tilde{\mu}}\\ =& - \frac{1}{n_{\varepsilon}V} \frac{\partial}{\partial\tilde{\mu}}\left[ V n_{\varepsilon} \int \mathrm{d}\tilde{q}\cdot \tilde{q}^2 \int_{-\infty}^{+\infty}\mathrm{d}\tilde{\omega} \cdot f(\tilde{q}, \tilde{\omega}, \tilde{\mu}, \tilde{\beta}) \right] \\ =& - \frac{\partial}{\partial\tilde{\mu}}\left[ \int \mathrm{d}\tilde{q}\cdot \tilde{q}^2 \int_{-\infty}^{+\infty}\mathrm{d}\tilde{\omega} \cdot f(\tilde{q}, \tilde{\omega}, \tilde{\mu}, \tilde{\beta}) \right] \end{align}

所以最终要求的为

\begin{align} \frac{\Delta F}{NE_n} =& \frac{\Omega}{NE_n} + \frac{\mu}{E_n} \\ =&\frac{\Omega}{N\varepsilon}\left( \frac{\varepsilon}{E_n} \right) + \tilde{\mu} \left( \frac{\varepsilon}{E_n} \right) \\ =& \left( \frac{n_{\varepsilon}}{n} \right)^{5/3} \int\mathrm{d}\tilde{q}\cdot \tilde{q}^2 \int_{-\infty}^{+\infty}\mathrm{d}\tilde{\omega} \cdot f(\tilde{q}, \tilde{\omega}, \tilde{\mu}, \tilde{\beta}) + \tilde{\mu} \left( \frac{n_{\varepsilon}}{n} \right)^{2/3} \\ \end{align}

横坐标为

\begin{align} \frac{2R}{k_n^2v} = \frac{2R}{k_{\varepsilon v}}\cdot \left( \frac{n_{\varepsilon}}{n} \right)^{2/3} \end{align}

2. result

3. code

计算 \(\Delta F\)

from matplotlib import pyplot as plt
import numpy as np
from scipy import integrate
#from scipy.misc import derivative
from scipy.integrate import fixed_quad
import time

start = time.process_time()

nn = 10

beta = 1
er = 1e-6
R = 1/30
epsabs = 1e-1

def xi(k, mu):
    return k**2 - mu
def n(k, mu):
    x = xi(k,mu)
#    print(x)
    n =  1 / (np.exp(beta*x) - 1)
    return n

def z(omega, q, mu):
    return omega - q**2/2 + 2*mu

def pi(omega, q, k, mu):
    pi = 1 + n(k+q/2, mu) + n(-k+q/2, mu)
    pi = pi / (xi(k+q/2, mu) + xi(-k+q/2, mu) -omega)
    pi = pi * k**4
    pi = pi -k**2/2 - z(omega, q, mu)/4
    pi = pi*2 / np.pi**2
    return pi

def PI(omega, q, mu):
    zz = z(omega, q, mu)
    if zz<0:
	PI, err = fixed_quad(lambda x: pi(omega, q, x, mu), er, 10,
			     n=nn) 
    else:
	a = np.sqrt(zz/2)
	PI1, err = fixed_quad(lambda x: pi(omega, q, x, mu), er, a-er,
			      n=nn) 
	PI2, err = fixed_quad(lambda x: pi(omega, q, x, mu), a+er,
			      10, n=nn)
	PI = PI1 + PI2
    PI = PI * R
    return PI

def delta(omega, q, rkv, mu):
    zz = z(omega, q, mu)
    if zz<0:
	img = 0
    else:
	k = np.sqrt(zz/2)
	img = 1 + n(k+q/2, mu) + n(-k+q/2, mu)
	img = img * R/(2*np.pi)
	img = img * k**3
    rel = PI(omega, q, mu)
    rel = rel + rkv/(4*np.pi)
    rel = rel +zz/(4*np.pi)
    delta = np.angle(rel + 1j*img) - np.pi
    return delta
def f(omega, q, rkv, mu):
    f = 1 / (np.exp(beta*omega) - 1)
    f = f * delta(omega, q, rkv, mu)
    f = 3 * f /np.pi
    return f

def F(rkv, mu):
    ff = lambda y, x: f(y, x, rkv, mu)
    F, err = integrate.dblquad(ff, er, 3, lambda x:er, lambda x:10, epsabs
			       = epsabs)
    return F


M = 1000
N = 10
x = np.linspace(0, 2, M)
y = np.zeros(M*N)
y.shape = (M, N)

mu = np.linspace(-2.1, -1.2, N)

for j in range(N):
    for i in range(M):
	y[i, j] = F(x[i], mu[j])
	print('mu_', j, 'y_', i, '=', y[i, j])

np.savetxt('y.txt', y)
print(y)

density = np.zeros(M*(N-2))
density.shape = (M, N-2)
dd = mu[1] - mu[0]

for j in range(N-2):
    for i in range(M):
	density[i, j] = y[i, j+2] - y[i, j]
	density[i, j] = - density[i, j] / (2*dd)
	print('mu_', j, 'density_', i, '=', density[i, j])

np.savetxt('density.txt', density)
print(density)

for i in range(N):
    plt.plot(x, y[:, i], label=r'$\mu/\epsilon$=%.2f' %mu[i])
plt.legend()

end = time.process_time()
print('time=', end-start, 'seconds')
plt.show()

计算 \(T_{C}\)

N = 100
rkv = np.linspace(-10, -.1, N)

f0 = np.zeros(N)
f1 = np.zeros(N)
density = np.zeros(N)

for i in range(N):
    print('rkv_', i, '=', rkv[i])
    f0[i] = F(rkv[i], -1e-3)
    print('f0_', i, '=', f0[i])
    f1[i] = F(rkv[i], -.1)
    print('f1_', i, '=', f1[i])
    density[i] = - (f0[i] - f1[i]) / .1
    print('density_', i, '=', density[i])

np.savetxt('f0.txt', f0)
np.savetxt('f1.txt', f1)
np.savetxt('density.txt', density)

x = np.zeros(N)
y = np.zeros(N)
for i in range(N):
    x[i] = rkv[i]/(density[i]**(2/3))
    y[i] = 1/(density[i]**(2/3))
plt.plot(x, y)

end = time.process_time()
print('time=', end-start, 'seconds')
plt.xlabel(r'$2R/(k_n^2 v)$')
plt.ylabel(r'$k_BT_C/E_n$')
plt.show()

def tm(omega, q, rkv, mu):
    zz = z(omega, q, mu)
    rel = PI(omega, q, mu)
    rel = rel + rkv/(4*np.pi)
    rel = rel +zz/(4*np.pi)
    return rel
N = 1000
M = 20
mu = np.linspace(-80, -1e-2, N)
y = np.zeros(N)
rkv = np.linspace(2, 100, M)
muRoot = np.zeros(M)
for j in range(M):
    c = 0
    for i in range(N):
	y[i] = tm(0, 0, rkv[j], mu[i])
	if np.abs(y[i])<np.abs(y[c]):
	    c = i
    print('y_', c, '=', y[c])
    muRoot[j] = mu[c]

    plt.plot(mu, y, label=r'$2R/(k_{\epsilon}^2v)=%.1f$'%rkv[j])
plt.legend()
plt.xlabel(r'$\mu/\epsilon$')
plt.ylabel(r'$T^{-1}$')
print(muRoot)
plt.show()

Tc = np.zeros(M)
Rn = np.zeros(M)

dens = np.zeros(M)
for i in range(M):
    print(i)
    dd = np.abs(muRoot[i]) * .1
    print('dd=', dd)
    f1 = F(rkv[i], muRoot[i]+dd)
    print('f1=', f1)
    f2 = F(rkv[i], muRoot[i]-dd)
    print('f2=', f2)
    nnn = - (f1 - f2) / (2*dd)
    print('nnn=', nnn)
    Tc[i] = 1 / (nnn**2/3)
    print('Tc=', Tc[i])
    Rn[i] = rkv[i] / (nnn**2/3)
    print('Rn=', Rn[i])
    dens[i] = nnn
    print('dens=', dens[i])
end = time.process_time()
print('time is', end-start, 'secends')
np.savetxt('density.txt', dens)
plt.plot(Rn, Tc)
plt.xlabel(r'$2R/(k_n^2 v)$')
plt.ylabel(r'$k_BT_C/E_n$')

plt.show()