# Untitled

unknown
plain_text
2 years ago
5.7 kB
2
Indexable
Never
import math
from scipy import optimize as opt
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,300,.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)

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
r = np.linspace(1e-5,2,100)
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]

# Finding lowest valid y
def findY0():
for i in range(1,integer_y0_max):
y0 = i
y, z = get_sol(y0)
print("Trying y =",(y0), end = "... ")
if (isCresc(y)):
prev = y0
else:
print("That's it!\nFound integer", (prev), "\ngoing decimal\n")
break
print("----------")

y0 = prev
for dPlace in range (1,7):
print("New decimal place: currently at", dPlace)
for j in range (1,10):
decy0 = y0 + j / 10 ** dPlace
y, z = get_sol(decy0)
print("trying dec y = " + (str)(decy0), end = "... ")
if (isCresc(y)):
prev = decy0
print("----------")
else:
y0 = prev;
print ("That's it!\ny must be", y0, "with further decimal places\n")
break
print ("Found y0 with", dPlace, "decimal places:", y0)
return y

# Calc if graph is crescent for y
def isCresc(y):
for indice in range(10):
if(y[indice] < y[indice+1]):
print("Still not it")
return True
else:
return False

# IVP for y
def get_sol(y0):
sol = integrate.solve_ivp(f, (1e-5,2),(y0,z0),t_eval=r)
return sol.y

phi = findY0()

plt.plot(r,phi)
plt.grid(True)