首页 > itarticle > nsr calculate

nsr calculate

admin 11月 13, 2020 0

Free Energy vs. Renormalize Paramaters

p-wave

cal

begin{align}
delta^p(vec{q}, z) =& mathrm{Arg}left[ frac{Mk_{n^2}}{2}frac{1}{R}
left(
frac{1}{4pi}cdot frac{2R}{k_n^2 v}
+ tilde{z}cdotfrac{1}{4pi}
+ frac{2 R}{M k_n^2}Pi_r(vec{q},z)
right)
right] \
=& mathrm{Arg}left[
frac{1}{4pi}cdot frac{2R}{k_n^2 v}
+ tilde{z}cdotfrac{1}{4pi}
+ 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 = 6pi^2n) , (n =
N/V)

begin{align}
frac{2 R}{M k_n^2}Pi_r(vec{q},z) =& (k_n R)cdotPi_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}{2pi^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}{(2pi)^3}int mathrm{d}tilde{k}
left[
k^2 cdot 4pi |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}cdottilde{k}^4left[
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}{(2pi^3)}int mathrm{d}^3vec{q}
cdot int frac{mathrm{d}omega}{pi}cdot frac{1}{e^{betaomega}-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{partialOmega}{partial mu}
end{align}

决定.

以 (varepsilon) 为单位

若以某一能量 (varepsilon) 为单位, 对应的长度单位 (k_{varepsilon} =
sqrt{2Mvarepsilon}) , 密度单位 (n_{varepsilon} =
k_{varepsilon}^3/(6pi^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}
cdotfrac{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{partialOmega}{partialmu}
=- frac{1}{n_{varepsilon}V}
frac{partialOmega/mu}{partialtilde{mu}}\
=& - frac{1}{n_{varepsilon}V}
frac{partial}{partialtilde{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}{partialtilde{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}{Nvarepsilon}left( frac{varepsilon}{E_n} right)
+ tilde{mu} left( frac{varepsilon}{E_n} right) \
=& left( frac{n_{varepsilon}}{n} right)^{5/3}
intmathrm{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}

result

code

计算 (Delta F) 

from matplotlib import pyplot as plt
import numpy as np
from scipy import integrate

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 (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()
微海报分享

近期文章

近期评论

标签

热门

文章归档

分类目录

功能