# Untitled

unknown
plain_text
2 years ago
4.5 kB
4
Indexable
Never
from scipy import optimize as opt
import math
import numpy as np
from scipy import integrate
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from scipy.special import zeta, gamma
mw2 = lambda phi: (m_w**2)*(phi**2/v**2)
mz2 = lambda phi: (m_z**2)*(phi**2/v**2)
mt2 = lambda phi: (m_t**2)*(phi**2/v**2)
cb = 223.6
cf = 14.8
pi = np.pi
m_h = 125
m_w = 80.385
m_z = 91.1876
m_t = 173.07
v = 246.22
M = 600
e = 0.01

# Low temperature approximation to the above integral
JLowT = lambda x: -np.sqrt(2)/4.*(x/pi)**1.5*np.exp(-x)*(1 + 15/(8.*x) + 105/(128.*x*x))
# High temperature
c = pi**2*np.exp(1.5-2*np.euler_gamma)
l = np.arange(20)
Cb = np.sqrt(pi)/2.*(-1)**l*zeta(2*l+3, 1)*gamma(l+1.5)/(gamma(l+4)*(2*pi)**(2*l+5))
Cf = (1-2**(2*l+3))*Cb
def JbHT(x,n=5):
return -pi**2/90. + x**2/24. - x**3/(12*pi)\
- x**4*np.log(np.maximum(abs(x**2),1e-100)/(16*c))/(64*pi**2)\
+ np.sum([Cb[i]*x**(2*i+6) for i in range(n)])
def JfHT(x,n=0):
return -7*pi**2/720. + x**2/48.\
+ x**4*np.log(np.maximum(x**2, 1e-100)/c)/(64*pi**2)\
+ np.sum([Cf[i]*x**(2*i+6) for i in range(n)])
def Jb(x):
""" Interpolation for bosonic contribution to thermal potential
(max error ~ 0.054%)
"""
xlim = 3.72439780262
if x.real < 0 or x.imag != 0:
return np.nan
if 0. <= x and x < xlim: return JbHT(x,n=5)
elif x > xlim: return JLowT(x)
elif np.isnan(x): return np.nan

def Jf(x):
""" Interpolation for bosonic contribution to thermal potential
(max error ~ 3.89%)
With this interpolation, the numerical derivative is continuous
and has a max error ~ 5.5% w.r.t. the full integral
"""
if x.real < 0 or x.imag != 0:
return np.nan
if 0. <= x and x < 1.1:
return -(-7*pi**2/720. + x**2/48. + x**4*np.log(np.maximum(x**2, 1e-100)/c)/(64*pi**2))
elif 1.1 <= x and x <= 3.4: return -(0.0855712017*x - 0.6086856492)/\
(6.321477772 - 0.7249924836*x + x**2)
elif x > 3.4:
return np.sqrt(2)/4.*(x/pi)**1.5*np.exp(-x)*(1 + 15/(8.*x))
elif np.isnan(x): return np.nan

def jb(y_t):
return ((2*pi**2)**-1)*quad((lambda x: x**2 * (np.log(1-np.exp(-np.sqrt(x**2+y_t**2))))), 0, np.infty)[0]
def jf(z_t):
return ((2*pi**2)**-1)*quad((lambda x: x**2 * (np.log(1+np.exp(-np.sqrt(x**2+z_t**2))))), 0, np.infty)[0]

def V(phi,T, M):
A = 3*v**4/(4*M**2) - m_h**2/2
B = -3*v**2/(2*M**2) + m_h**2/(2*v**2)
V0 = A*(phi**2)/2 + (B*phi**4)/4 + (1/(8*M**2))*phi**6
MASSA_W = ((mw2(phi)**2)*np.log(mw2(phi)/m_w**2)) -((3*mw2(phi)**2)/2) + (2*mw2(phi))*(m_w**2) - (m_w**4)/2
MASSA_Z = ((mz2(phi)**2)*np.log(mz2(phi)/m_z**2)) -((3*mz2(phi)**2)/2) + (2*mz2(phi))*(m_z**2) - (m_z**4)/2
MASSA_t = ((mt2(phi)**2)*np.log(mt2(phi)/m_t**2)) -((3*mt2(phi)**2)/2) + (2*mt2(phi))*(m_t**2) - (m_t**4)/2
V1 = (1/(64*(pi**2)))*(6*MASSA_W+3*MASSA_Z-12*MASSA_t)
if T!=0:
del_V1 = (T**4*((6*Jb(np.sqrt(mw2(phi))/T) + (3*Jb(np.sqrt(mz2(phi))/T)) - (12*Jf(np.sqrt(mt2(phi))/T)))))
return V0+V1+del_V1
else:
return V0+V1
def v_linha(phi,T):
return V(phi,T, M)-V(1e-15,T, M)

def find_tc(low, high):
def vmin(T):
phimin = minimize(lambda phi: v_linha(phi,T),v).x
return v_linha(phimin, T)
return opt.bisect(vmin, low, high)

# Plot effective potential graph
plt.figure(num=0,dpi=120.0)
philist = np.arange(1e-15,250,.1)
plt.plot(philist, list(map(lambda x: v_linha(x,80), philist)))
plt.xlabel('$\phi$')
plt.ylabel('V($\phi$,T)')
plt.title('Potencial Efetivo')
plt.grid(True)
plt.show()

def minv_T(Tc):
""" Find the minimum value for phi (x.axis) """
minimo_v = minimize(lambda phi: v_linha(phi,Tc),v).x
csi = minimo_v / Tc
return minimo_v, csi

# Sphaleron with boundary conditions: y0;v0

rf = 2
r = np.linspace(1e-5, rf, 1000)
z0 = 0
T = 80

y0_max = minv_T(T)[0]
integer_y0_max =  math.ceil(y0_max)
print(y0_max, integer_y0_max)

def deriv1(phi):
return (V(phi+e,T,M) - V(phi-e,T,M))/(2*e)

def f(r,S):
y,z = S
dy_dr = z
dz_dr = (-2*z)/r + deriv1(y)
return dy_dr, dz_dr

# ivp for y
def getSol(y0):
sol = integrate.solve_ivp(f, (1e-5, rf), (y0, z0), t_eval=r)
return sol.y