#code author: Zhiqi Huang, huangzhq25@mail.sysu.edu.cn import numpy as np from scipy.interpolate import interp1d import matplotlib.pyplot as plt from os import path plot_background = False iso_curv = False plot_item = 'eps' kMpc = np.array([ 0.05 ]) bound_scale = 0.2 class two_field_inflation_model: """The model defined in Eq.(2.1) of 2004.00672""" Mpl = 256. #reduced Planck mass, you can set it to be any positive number check_conserve = "EP" #set None for doing nothing; to "EP" for energy-momentum conservation; or "ZETA" for super-horizon curvature perturbation conservation def __init__(self, As, num_efolds, mu_by_Mpl, mphi_by_Mpl = 5.e-6, mchi_by_Mpl = 0., b1Mpl = 0.): self.Mplsq = self.Mpl**2 self.As = As self.num_efolds = num_efolds self.mu = mu_by_Mpl * self.Mpl self.mphi = mphi_by_Mpl * self.Mpl self.mchi = mchi_by_Mpl * self.Mpl self.b1 = b1Mpl/self.Mpl self.verbose = True def func_b(self, phi): #e^{2b(phi)}/2 * (\partial chi)^2 return self.b1 * phi def dbdphi(self, phi): return self.b1 def d2bdphi2(self, phi): return 0. def func_V(self, phi, chi): return (self.mphi * phi)**2 / 2. + (self.mchi * chi)**2 / 2. def dVdphi(self, phi, chi): return self.mphi**2 * phi def dVdchi(self, phi, chi): return self.mchi**2 * chi def d2Vdphi2(self, phi, chi): return self.mphi**2 def d2Vdchi2(self, phi, chi): return self.mchi**2 def d2Vdphidchi(self, phi, chi): return 0. def epsilon_V(self, phi, chi): return (self.dVdphi(phi, chi)**2 + self.dVdchi(phi, chi)**2*np.exp(-2.*self.func_b(phi)))*(self.Mpl/self.func_V(phi, chi))**2/2. def bg_eqs(self, y): # y = [H, phi, dot phi, chi, dot chi, lna] e2b = np.exp(2.*self.func_b(y[1])) bphi = self.dbdphi(y[1]) return np.array([ -(y[2]**2 + e2b*y[4]**2)/(2.*self.Mplsq), y[2], bphi*e2b*y[4]**2-self.dVdphi(y[1], y[3])-3.*y[0]*y[2], y[4], -self.dVdchi(y[1], y[3])/e2b - (3.*y[0]+2.*bphi*y[2])*y[4], y[0] ]) def bg_evolve(self, y, dt): k1 = self.bg_eqs(y) k2 = self.bg_eqs(y+k1 * (dt/2.)) k3 = self.bg_eqs(y+k2 * (dt/2.)) k4 = self.bg_eqs(y + k3*dt) return y+(dt/6.)*(k1+(k2+k3)*2.+k4) def get_bg(self, phi_ini, chi_ini, epsilon_end = 1., nsteps = 8192): n = nsteps t = np.empty(n) y = np.empty((n, 6)) eps = np.empty(n) V = self.func_V(phi_ini, chi_ini) V_phi = self.dVdphi(phi_ini, chi_ini) V_chi = self.dVdchi(phi_ini, chi_ini) e2b = np.exp(2.*self.func_b(phi_ini)) bphi = self.dbdphi(phi_ini) self.phi_ini = phi_ini self.chi_ini = chi_ini self.dotphi_ini = 0. self.dotchi_ini = 0. for i in range(7): self.H_ini = np.sqrt((V+(self.dotphi_ini**2 + self.dotchi_ini**2 * e2b )/2.)/3.)/self.Mpl self.dotphi_ini = (bphi*e2b*self.dotchi_ini**2-V_phi)/3./self.H_ini self.dotchi_ini = -V_chi/e2b / (3.*self.H_ini + 2.*bphi *self.dotphi_ini) self.H_ini = np.sqrt((V+(self.dotphi_ini**2 + self.dotchi_ini**2 * e2b )/2.)/3.)/self.Mpl y[0, :] = np.array([self.H_ini, phi_ini, self.dotphi_ini, chi_ini, self.dotchi_ini, 0.]) #H, phi, dot phi, chi, dot chi, ln a eps[0] = (self.dotphi_ini ** 2 + self.dotchi_ini **2 *e2b)/2./(self.Mpl*self.H_ini)**2 if(eps[0] > epsilon_end): print('initial condition error: not in an inflationary regime') exit() step = (self.num_efolds + 10.) / n self.t_end = -1.e30 t[0] = 0. n_used = n while(self.t_end < 0.): step *= 1.2 for i in range(1, n): dt = step/y[i-1, 0]/(1.+eps[i-1]*10.) y[i, :] = self.bg_evolve(y[i-1, :], dt) t[i] = t[i-1]+dt eps[i] = (y[i, 2]**2 + y[i, 4]**2*np.exp(2.*self.func_b(y[i, 1])))/2./(self.Mpl*y[i, 0])**2 if(self.t_end < 0. and eps[i-1] < epsilon_end and eps[i] >= epsilon_end): n_used = i+1 acc = abs( ((y[i, 2]**2 + y[i, 4]**2 * np.exp(2.*self.func_b(y[i,1])))/2. + self.func_V(y[i, 1], y[i,3]) )/(3.*(y[i, 0]*self.Mpl)**2)-1.) if(acc > 1.e-4): print("phi = %10.3f M_p, chi = %10.3f M_p, i=%8i" %(phi_ini/self.Mpl, chi_ini/self.Mpl, i)) print('Energy conservation check fails. Relative error: %10.2e' % (acc)) exit() self.t_end = ((eps[i]-epsilon_end)*t[i-1] + (epsilon_end - eps[i-1])*t[i])/(eps[i]-eps[i-1]) self.lna_end = ((eps[i]-epsilon_end)*y[i-1, 5] + (epsilon_end - eps[i-1])*y[i, 5])/(eps[i]-eps[i-1]) break y[0:n_used, 5] += (self.num_efolds - self.lna_end) self.lna_ini = y[0, 5] if(self.lna_ini <= 0.): i = 0 while(y[i, 5] <= 0. ): i += 1 self.t_start = (y[i, 5]*t[i-1] - y[i-1, 5]*t[i])/(y[i, 5]-y[i-1, 5]) else: self.t_start = -1.e30 self.background = interp1d(x = t[0:n_used], y = y[0:n_used, :], axis = 0, kind="cubic", bounds_error = True, copy=False, assume_sorted=True) self.t_of_lnaH = interp1d(x = y[0:n_used, 5]+np.log(y[0:n_used, 0]), y = t[0:n_used], axis = 0, kind="linear", bounds_error = True, copy=False, assume_sorted=False) self.lnaH_min = y[0, 5] + np.log(y[0, 0]) self.lnaH_max = y[n_used-1, 5] + np.log(y[n_used-1, 0]) y_end = self.background(self.t_end) self.H_end = y_end[0] self.phi_end = y_end[1] self.chi_end = y_end[3] self.dotphi_end = y_end[2] self.dotchi_end = y_end[4] self.t_end = t[n_used-1] if(self.t_start >= 0.): y_start = self.background(self.t_start) self.H_start= y_start[0] self.phi_start= y_start[1] self.chi_start = y_start[3] self.dotphi_start = y_start[2] self.dotchi_start = y_start[4] eb = np.exp(self.func_b(self.phi_start)) dsigma_sq = self.dotphi_start ** 2 + (eb*self.dotchi_start)**2 self.epsilon_start = dsigma_sq/2./(self.H_start*self.Mpl)**2 self.etav_start = (self.d2Vdphi2(self.phi_start, self.chi_start)*self.dotphi_start**2 + 2.*self.d2Vdphidchi(self.phi_start, self.chi_start)*self.dotphi_start * self.dotchi_start + self.d2Vdchi2(self.phi_start, self.chi_start)*self.dotchi_start**2)/dsigma_sq*self.Mplsq/2./self.func_V(self.phi_start, self.chi_start) self.r_slowroll = self.epsilon_start * 16. self.As_slowroll = self.H_start**2 / (8.*np.pi**2*self.Mplsq*self.epsilon_start) y_shift = self.background(self.t_start + 0.2/self.H_start) eps_shift = (y_shift[2]**2 + np.exp(2.*self.func_b(y_shift[1]))*y_shift[4]**2)/(y_shift[0]*self.Mpl)**2/2. lnk_shift = y_shift[5]-y_start[5]+np.log(y_shift[0]/y_start[0]) As_shift = y_shift[0]**2/(8.*np.pi**2*self.Mpl**2*eps_shift) self.ns_slowroll = 1.+np.log(As_shift/self.As_slowroll)/lnk_shift self.nt_slowroll = -2.*self.epsilon_start def get_bg_long(self, phi_ini, chi_ini, t_end, nsteps = 16384): n = nsteps n1 = n * 2 // 3 n2 = n - n1 self.t_start = 0. t = np.empty(n) dt1 = (self.t_end - self.t_start)/(n1-1.) dt2 = (t_end - self.t_end)/n2 t[0 : n1] = np.linspace(self.t_start, self.t_end, n1) t[n1 : n] = np.linspace(self.t_end + dt2, t_end, n2) self.t_infend = self.t_end #inflation end self.t_end = t_end #actual end y = np.empty((n, 6)) V = self.func_V(phi_ini, chi_ini) V_phi = self.dVdphi(phi_ini, chi_ini) V_chi = self.dVdchi(phi_ini, chi_ini) e2b = np.exp(2.*self.func_b(phi_ini)) bphi = self.dbdphi(phi_ini) self.phi_ini = phi_ini self.chi_ini = chi_ini self.dotphi_ini = 0. self.dotchi_ini = 0. for i in range(5): self.H_ini = np.sqrt((V+(self.dotphi_ini**2 + self.dotchi_ini**2 * e2b )/2.)/3.)/self.Mpl self.dotphi_ini = (bphi*e2b*self.dotchi_ini**2-V_phi)/3./self.H_ini self.dotchi_ini = -V_chi/e2b / (3.*self.H_ini + 2.*bphi *self.dotphi_ini) self.H_ini = np.sqrt((V+(self.dotphi_ini**2 + self.dotchi_ini**2 * e2b )/2.)/3.)/self.Mpl y[0, :] = np.array([self.H_ini, phi_ini, self.dotphi_ini, chi_ini, self.dotchi_ini, 0.]) for i in range(1, n1): y[i, :] = self.bg_evolve(y[i-1, :], dt1) lna_shift = self.num_efolds - y[i, 5] self.H_infend = y[i, 0] for i in range(n1, n): y[i, :] = self.bg_evolve(y[i-1, :], dt2) y[:, 5] += lna_shift self.lna_end = y[-1, 5] self.H_end = y[-1, 0] self.background = interp1d(x = t, y = y, axis = 0, kind="cubic", bounds_error = True, copy=False, assume_sorted=True) iup = n-1 idown = 0 assert(y[iup, 5] > 0. and y[idown, 5] < 0.) while(iup - idown > 1): imid = (iup+idown) // 2 if(y[imid, 5] > 0.): iup = imid else: idown = imid self.t_pivot = t[iup] - y[iup, 5]/y[iup, 0] self.y_pivot = self.background(self.t_pivot) self.lnaH_pivot = self.y_pivot[5] + np.log(self.y_pivot[0]) def try_mphi(self, mphi, mchi_by_mphi = 0., chi_ini = 0., evolve_long=False, phi_ini = 0.): self.mphi = mphi self.mchi = mchi_by_mphi * mphi if(phi_ini == 0.): phi_ini = self.Mpl * 0.1 while( self.epsilon_V(phi_ini, chi_ini) > 0.002): phi_ini *= 1.001 if(phi_ini > 30.*self.Mpl): print('initial chi may not work: ', chi_ini/self.Mpl) exit() while( self.epsilon_V(phi_ini, chi_ini) < 1.e-8): phi_ini *= 0.9999 if(phi_ini < 0.01*self.Mpl): print('initial chi may not work: ', chi_ini/self.Mpl) exit() self.get_bg(phi_ini = phi_ini, chi_ini = chi_ini, nsteps = 8192) while(self.lna_ini < -50.): print(phi_ini/self.Mpl, self.lna_ini) phi_ini *= 0.999 self.get_bg(phi_ini = phi_ini, chi_ini = chi_ini, nsteps = 8192) while(self.lna_ini > -8.): print(phi_ini/self.Mpl, self.lna_ini) phi_ini *=1.001 self.get_bg(phi_ini = phi_ini, chi_ini = chi_ini, nsteps = 8192) print(phi_ini/self.Mpl, self.lna_ini) else: self.get_bg(phi_ini = phi_ini, chi_ini = chi_ini, nsteps = 8192) if(evolve_long): self.get_bg_long(phi_ini = phi_ini, chi_ini = chi_ini, t_end = self.t_end*1.06, nsteps = 30000) else: self.get_bg(phi_ini = phi_ini, chi_ini = chi_ini, nsteps = 20000) if(self.t_start < 0.): print("Error: try_mphi cannot find initial conditions") exit() def set_mphi(self, mchi_by_mphi, chi_ini = 0., mphi_by_Mpl = None, phi_ini = 0.): if(mphi_by_Mpl is not None): self.try_mphi(mphi = mphi_by_Mpl*self.Mpl, mchi_by_mphi = mchi_by_mphi, chi_ini = chi_ini, evolve_long=True, phi_ini = phi_ini) return mphi_b = 1.e-6*self.Mpl mphi_u = 1.e-4*self.Mpl self.try_mphi(mphi = mphi_b, mchi_by_mphi = mchi_by_mphi, chi_ini = chi_ini, phi_ini = phi_ini) As_b = self.As_slowroll self.try_mphi(mphi = mphi_u, mchi_by_mphi = mchi_by_mphi, chi_ini = chi_ini, phi_ini = phi_ini) As_u = self.As_slowroll if((As_b - self.As)*(As_u - self.As) > 0.): print('set_mphi fails, you may want to set the mass of inflaton manually!') exit() mphi_m = np.exp((np.log(mphi_u)*(np.log(self.As)- np.log(As_b)) + np.log(mphi_b)*(np.log(As_u) - np.log(self.As)))/(np.log(As_u) - np.log(As_b))) self.try_mphi(mphi = mphi_m, mchi_by_mphi = mchi_by_mphi, chi_ini = chi_ini, phi_ini = phi_ini) As_m = self.As_slowroll if(abs(As_m/self.As - 1.) < 1.e-4): return if(abs(As_m/self.As - 1.) < 0.02): mphi_b = mphi_m * 0.8 mphi_u = mphi_m * 1.2 self.try_mphi(mphi = mphi_b, mchi_by_mphi = mchi_by_mphi, chi_ini = chi_ini, phi_ini = phi_ini) As_b = self.As_slowroll self.try_mphi(mphi = mphi_u, mchi_by_mphi = mchi_by_mphi, chi_ini = chi_ini, phi_ini = phi_ini) As_u = self.As_slowroll if((As_b - self.As)*(As_u - self.As) > 0.): print('set_mphi fails, you may want to set the mass of inflaton manually!!') exit() niter = 0 while(mphi_u / mphi_b > 1.001): niter += 1 if(niter < 4): mphi_m = np.exp((np.log(mphi_u)*(np.log(self.As)- np.log(As_b)) + np.log(mphi_b)*(np.log(As_u) - np.log(self.As)))/(np.log(As_u) - np.log(As_b))) else: #for non-smooth models that do not converge quickly mphi_m = np.sqrt(mphi_b * mphi_u) self.try_mphi(mphi = mphi_m, mchi_by_mphi = mchi_by_mphi, chi_ini = chi_ini, phi_ini = phi_ini) As_m = self.As_slowroll if(abs(As_m/self.As - 1.) < 1.e-4): return if( (As_m - self.As) * (As_u - self.As) <= 0.): mphi_b = mphi_m As_b = As_m elif( (As_m - self.As) * (As_b - self.As) <= 0.): mphi_u = mphi_m As_u = As_m else: print('set_mphi fails, you may want to set the mass of inflaton manually!!!') exit() def check_potential(self): #now I check the potential to avoid typos dphi = self.phi_start/200. phi = self.phi_end + (self.phi_start - self.phi_end)*np.random.rand() dchi = dphi*0.789 chi = np.random.normal() * 0.1*self.Mpl dV = (self.func_V(phi + dphi, chi) - self.func_V(phi-dphi, chi))/(2.*dphi) if(abs(dV - self.dVdphi(phi, chi))/(1.e-12*self.Mpl**3 +abs(dV)) > 1.e-4): print('Warning: you may have some typo in dVdphi: %16.7e ~=? %16.7e' % ( dV, self.dVdphi(phi, chi))) dV = (self.func_V(phi , chi+dchi) - self.func_V(phi, chi-dchi))/(2.*dchi) if(abs(dV - self.dVdchi(phi, chi))/(1.e-12*self.Mpl**3 +abs(dV)) > 1.e-4): print('Warning: you may have some typo in dVdchi: %16.7e ~=? %16.7e' % (dV, self.dVdchi(phi, chi))) d2V = (self.dVdphi(phi + dphi, chi) - self.dVdphi(phi-dphi, chi))/(2.*dphi) if(abs(d2V - self.d2Vdphi2(phi, chi))/(1.e-9*self.Mplsq+abs(d2V)) > 1.e-4): print('Warning: you may have some typo in d2Vdphi2: %16.7e ~=? %16.7e' % (d2V, self.d2Vdphi2(phi, chi))) d2V = (self.dVdchi(phi, chi+dchi) - self.dVdchi(phi, chi-dchi))/(2.*dchi) if(abs(d2V - self.d2Vdchi2(phi, chi))/(1.e-9*self.Mplsq+abs(d2V)) > 1.e-4): print('Warning: you may have some typo in d2Vdchi2: %16.7e ~=? %16.7e' % (d2V, self.d2Vdchi2(phi, chi))) d2Va = (self.dVdchi(phi + dphi, chi) - self.dVdchi(phi-dphi, chi))/(2.*dphi) d2Vb = (self.dVdphi(phi , chi+dchi) - self.dVdphi(phi, chi-dchi))/(2.*dchi) d2Vc = self.d2Vdphidchi(phi, chi) if(abs(d2Va - d2Vc)/(1.e-9*self.Mplsq+abs(d2V)) > 1.e-4 or abs(d2Vb - d2Vc)/(1.e-9*self.Mplsq+abs(d2V)) > 1.e-4): print('Warning: you may have some typo in d2Vdphidchi: %16.7e ~=? %16.7e ~=? %16.7e' % (d2Va, d2Vb, d2Vc )) def perturb_eqs(self, t, lnk, z): #z = [ Re (\Phi, \dot\Phi, \delta\phi, \dot{\delta\phi}, \delta\chi, \dot{\delta\chi}), Im (...) y = self.background(t) # y = [H, phi, dot phi, chi, dot chi, lna] yp = self.bg_eqs(y) kbya = np.exp(lnk - y[5]) kbyaH = kbya/y[0] kbya_sq = kbya**2 eb = np.exp(self.func_b(y[1])) e2b = eb**2 bphi = self.dbdphi(y[1]) bphiphi = self.d2bdphi2(y[1]) Vphi = self.dVdphi(y[1], y[3]) Vchi = self.dVdchi(y[1], y[3]) Vphiphi = self.d2Vdphi2(y[1], y[3]) Vphichi = self.d2Vdphidchi(y[1], y[3]) Vchichi = self.d2Vdchi2(y[1], y[3]) zp = np.empty(12) zp[0] = z[1] if( kbyaH > 3. ): #subhorizon mode zp[1] = -7.*y[0]*z[1] - (2.*yp[0]+6.*y[0]**2+kbya_sq)*z[0] - (Vphi*z[2]+Vchi*z[4])/self.Mplsq zp[7] = -7.*y[0]*z[7] - (2.*yp[0]+6.*y[0]**2+kbya_sq)*z[6] - (Vphi*z[8]+Vchi*z[10])/self.Mplsq elif(kbyaH > 0.3): zp[1] = -4.*y[0]*z[1] - (yp[0]+3.*y[0]**2)*z[0] + (y[2]*z[3] + e2b*y[4]*z[5]+(bphi*e2b*y[4]**2-Vphi)*z[2]-Vchi*z[4])/(2.*self.Mplsq) zp[7] = -4.*y[0]*z[7] - (yp[0]+3.*y[0]**2)*z[6] + (y[2]*z[9] + e2b*y[4]*z[11]+(bphi*e2b*y[4]**2-Vphi)*z[8]-Vchi*z[10])/(2.*self.Mplsq) else: zp[1] = -y[0]*z[1] + kbya_sq*z[0] + (y[2]*z[3] + e2b*y[4]*z[5] + bphi*e2b*y[4]**2*z[2])/self.Mplsq zp[7] = -y[0]*z[7] + kbya_sq*z[6] + (y[2]*z[9] + e2b*y[4]*z[11]+ bphi*e2b*y[4]**2*z[8])/self.Mplsq zp[2] = z[3] zp[3] = -3.*y[0]*z[3] - (kbya_sq + Vphiphi - (bphiphi + 2*bphi**2)*y[4]**2*e2b)*z[2] - Vphichi*z[4] + 2*bphi*e2b*y[4]*z[5] + 4.*y[2]*z[1]- 2.*Vphi*z[0] zp[4] = z[5] zp[5] = -(3.*y[0]+2.*bphi*y[2])*z[5]-(kbya_sq + Vchichi/e2b)*z[4] - 2.*bphi*y[4]*z[3]-((Vphichi - 2.*bphi*Vchi)/e2b+2.*bphiphi*y[2]*y[4])*z[2] + 4.*y[4]*z[1]-2.*Vchi*z[0]/e2b zp[6] = z[7] zp[8] = z[9] zp[9] = -3.*y[0]*z[9] - (kbya_sq + Vphiphi - (bphiphi + 2*bphi**2)*y[4]**2*e2b)*z[8] - Vphichi*z[10] + 2*bphi*e2b*y[4]*z[11] + 4.*y[2]*z[7]- 2.*Vphi*z[6] zp[10] = z[11] zp[11] = -(3.*y[0]+2.*bphi*y[2])*z[11]-(kbya_sq + Vchichi/e2b)*z[10] - 2.*bphi*y[4]*z[9]-((Vphichi - 2.*bphi*Vchi)/e2b+2.*bphiphi*y[2]*y[4])*z[8]+ 4.*y[4]*z[7]-2.*Vchi*z[6]/e2b return zp def perturb_evolve(self, t, lnk, z, dt): k1 = self.perturb_eqs(t, lnk, z) k2 = self.perturb_eqs(t+dt/2., lnk, z+k1 * (dt/2.)) k3 = self.perturb_eqs(t+dt/2., lnk, z+k2 * (dt/2.)) k4 = self.perturb_eqs(t+dt, lnk, z + k3*dt) z += (dt/6.)*(k1+(k2+k3)*2.+k4) def get_perturb(self, lnk, iso_curvature = False, nsave = 2048): lnks = lnk - 4.6 #~k/100 t_ini = self.t_pivot + (lnks - self.lnaH_pivot)/self.y_pivot[0] assert(t_ini >= self.t_start and t_ini < self.t_end) for i in range(3): y_ini = self.background(t_ini) t_ini += (lnks - y_ini[5] - np.log(y_ini[0]))/y_ini[0] assert(t_ini >= self.t_start and t_ini < self.t_end) n1 = nsave // 3 n2 = nsave - n1 dt1 = (self.t_infend - t_ini)/n1 dt2 = (self.t_end - self.t_infend)/n2/1.25 #leave some space to t_end because we will do random drift y = self.background(t_ini) yp = self.bg_eqs(y) kbya = np.exp(lnk - y[5]) eb = np.exp(self.func_b(y[1])) e2b = eb**2 bphi = self.dbdphi(y[1]) bphiphi = self.d2bdphi2(y[1]) Vphi = self.dVdphi(y[1], y[3]) Vchi = self.dVdchi(y[1], y[3]) dsigma = np.sqrt(y[2]**2 + y[4]**2*e2b) cos_theta = y[2]/dsigma sin_theta = eb*y[4]/dsigma amp = kbya/(2.*np.pi) z = np.zeros(12) if(iso_curvature): z[2] = -amp*sin_theta #\delta\phi z[4] = amp*cos_theta/eb #\delta\chi else: z[2] = amp*cos_theta #\delta\phi z[4] = amp*sin_theta/eb #\delta\chi z[9] = -z[2]*kbya z[11] = -z[4]*kbya z[0] = -(y[2]*z[3] + e2b*y[4]*z[5] + (bphi*e2b*y[4]**2+Vphi)*z[2] + Vchi*z[4]+3.*y[0]*(y[2]*z[2]+e2b*y[4]*z[4]))/(2.*self.Mplsq)/(kbya**2+yp[0]) #Re(newtonian potential) z[6] = -(y[2]*z[9] + e2b*y[4]*z[11] + (bphi*e2b*y[4]**2+Vphi)*z[8] + Vchi*z[10]+3.*y[0]*(y[2]*z[8]+e2b*y[4]*z[10]))/(2.*self.Mplsq)/(kbya**2+yp[0]) #Im(newtonian potential) z[1] = (y[2]*z[2] + e2b*y[4]*z[4])/(2.*self.Mplsq) - y[0]*z[0] z[7] = (y[2]*z[8] + e2b*y[4]*z[10])/(2.*self.Mplsq) - y[0]*z[6] t = t_ini lnalist = np.empty(nsave) tlist = np.empty(nsave) Hlist = np.empty(nsave) delta_phi = np.empty(nsave) delta_chi = np.empty(nsave) Phi = np.empty(nsave) zeta = np.empty(nsave) epslist = np.empty(nsave) mplist = np.empty(nsave) mplist[0] = np.sqrt((z[2]*y[2]+e2b*z[4]*y[4])**2 + (z[8]*y[2]+e2b*z[10]*y[4])**2)/self.Mpl/y[0]**2 tlist[0] = t epslist[0] = -yp[0]/y[0]**2 Hlist[0] = y[0] lnalist[0] = y[5] Phi[0] = np.sqrt(z[0]**2 + z[6]**2) #newtonian potential delta_phi[0] = np.sqrt(z[2]**2+z[8]**2) #delta phi delta_chi[0] = np.sqrt(z[4]**2+z[10]**2) #delta chi zeta[0] = np.sqrt((z[0]+y[0]*(z[2]*cos_theta+eb*z[4]*sin_theta)/dsigma)**2 + (z[6]+y[0]*(z[8]*cos_theta+eb*z[10]*sin_theta)/dsigma)**2) nsteps = np.array([ 8, 32, 64, 128, 256, 512 ], dtype=np.int32) dts = dt1 / nsteps dsigma_bound = 0. for i in range(1, nsave): if(i == n1): #switch to smaller recording time dts = dt2/nsteps if(lnalist[i-1] < self.num_efolds-0.2): if( lnalist[i-1]+np.log(Hlist[i-1]) < lnk): istep = 1 elif(lnalist[i-1] > self.num_efolds - 0.4): istep = 2 else: istep = 0 else: if(lnalist[i-1] > self.num_efolds + 0.4): istep = 5 elif(lnalist[i-1] > self.num_efolds): istep = 4 else: istep = 3 for j in range(nsteps[istep]): self.perturb_evolve(t, lnk, z, dts[istep]) t += dts[istep] y = self.background(t) #y = [H, phi, dot phi, chi, dot chi, lna] dsigma = np.sqrt(y[2]**2 + y[4]**2*e2b) while(dsigma < dsigma_bound): #where the expression of zeta has a 0/0 uncertainty, we evolve the system a random dt rdt = dts[istep]*np.random.rand() self.perturb_evolve(t, lnk, z, rdt) t += rdt y = self.background(t) #y = [H, phi, dot phi, chi, dot chi, lna] dsigma = np.sqrt(y[2]**2 + y[4]**2*e2b) dsigma_bound = (dsigma*bound_scale + dsigma_bound*99)/100. #dynamically adjust dsigma_bound to bound_scale* dsigma Phi[i] = np.sqrt(z[0]**2 + z[6]**2) #newtonian potential delta_phi[i] = np.sqrt(z[2]**2+z[8]**2) #delta phi delta_chi[i] = np.sqrt(z[4]**2+z[10]**2) #delta chi Hlist[i] = y[0] tlist[i] = t epslist[i] = 0.5*dsigma**2/self.Mplsq/Hlist[i]**2 #\epsilon = -\dot{H}/H^2 mplist[i] = np.sqrt((z[2]*y[2]+e2b*z[4]*y[4])**2 + (z[8]*y[2]+e2b*z[10]*y[4])**2)/self.Mpl/y[0]**2 eb = np.exp(self.func_b(y[1])) e2b = eb**2 cos_theta = y[2]/dsigma sin_theta = eb*y[4]/dsigma zeta[i] = np.sqrt((z[0]+y[0]*(z[2]*cos_theta+eb*z[4]*sin_theta)/dsigma)**2 + (z[6]+y[0]*(z[8]*cos_theta+eb*z[10]*sin_theta)/dsigma)**2) lnalist[i] = y[5] if( i % 10 == 0): mom1 = z[1] + y[0]*z[0] mom2 = dsigma/(2.*self.Mplsq)*(z[2]*cos_theta+eb*z[4]*sin_theta) print("%8.5f %8.3e %12.5e %12.5e %10.4e %10.4e" % (lnalist[i], (mom1-mom2)/(mom1+mom2), mom1, zeta[i], delta_phi[i]/self.Mpl, delta_chi[i]/self.Mpl)) return Hlist, zeta, Phi, delta_phi, delta_chi, lnalist, epslist, tlist, mplist class alpha_attractor_E(two_field_inflation_model): def func_V(self, phi, chi): return ((self.mphi*self.mu)**2/2.)*(1.-np.exp(-phi/self.mu))**2 + (self.mchi*chi)**2/2. def dVdphi(self, phi, chi): ephi = np.exp(-phi/self.mu) return (self.mphi**2 * self.mu ) * (1.- ephi)*ephi def d2Vdphi2(self, phi, chi): ephi = np.exp(-phi/self.mu) return (self.mphi**2) * (2.* ephi - 1.)*ephi class alpha_attractor_T(two_field_inflation_model): def func_V(self, phi, chi): return ((self.mphi*self.mu)**2/2.) * np.tanh(phi/self.mu)**2 + (self.mchi*chi)**2/2. def dVdphi(self, phi, chi): return (self.mphi**2 * self.mu ) * np.tanh(phi/self.mu)*(1.-np.tanh(phi/self.mu)**2) def d2Vdphi2(self, phi, chi): return (self.mphi**2 ) * (1.-np.tanh(phi/self.mu)**2) * (1. - 3.* np.tanh(phi/self.mu)**2) class alpha_attractor_P(two_field_inflation_model): def func_V(self, phi, chi): return ((self.mphi*self.mu)**2/2.) * phi**2/(phi**2+self.mu**2)+ (self.mchi*chi)**2/2. def dVdphi(self, phi, chi): return (self.mphi**2 * self.mu**4) * phi/(phi**2+self.mu**2)**2 def d2Vdphi2(self, phi, chi): return (self.mphi**2 * self.mu**4) /(phi**2+self.mu**2)**2*(1. - 4.* phi**2/(phi**2+self.mu**2)) tfm = alpha_attractor_T(As = 2.1e-9, num_efolds = 50.3, mu_by_Mpl = 0.04, b1Mpl = 25.) tfm.set_mphi(mchi_by_mphi = 0.3, chi_ini = 1.e-6*tfm.Mpl, mphi_by_Mpl=7.216e-6, phi_ini = 0.2652*tfm.Mpl) #print("checking the potential") #tfm.check_potential() #print("initial phi = %10.3f M_p" %(tfm.phi_ini/tfm.Mpl)) #print("m_phi = %10.3e M_p" %(tfm.mphi/tfm.Mpl)) #print("A_s = %10.4e, n_s = %6.3f, r = %10.3e n_t = %10.3e" % ( tfm.As_slowroll, tfm.ns_slowroll, tfm.r_slowroll, tfm.nt_slowroll)) colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf'] lss = ['solid', 'dotted', 'dashed', 'dashdot', 'dotted', 'dashed'] kp = 0.05 kbykps = kMpc/kp if(plot_background): n = 5000 t = np.linspace(tfm.t_start, tfm.t_end, n) lna = np.empty(n) phi = np.empty(n) chi = np.empty(n) for i in range(n): y = tfm.background(t[i]) phi[i] = y[1]/tfm.Mpl chi[i] = y[3]/tfm.Mpl lna[i] = y[5] plt.xlabel(r'$\ln a$') plt.xlim([tfm.num_efolds-1., tfm.lna_end]) plt.plot(lna, phi, label=r'$\phi$') plt.plot(lna, chi*1.e6, ls='dotted', label=r'$10^6\chi$') plt.legend() plt.show() else: plt.xlabel(r'$\ln a$') if(plot_item == 'zeta'): plt.ylabel(r'$\zeta$') plt.yscale('log') elif(plot_item == 'Phi'): plt.ylabel(r'$\Phi$') plt.yscale('log') elif(plot_item == 'delta_phi'): plt.ylabel(r'$|\delta\phi|/M_p$') plt.xlim([tfm.num_efolds-1., tfm.lna_end]) plt.yscale('log') elif(plot_item == 'delta_chi'): plt.xlim([tfm.num_efolds-1., tfm.lna_end]) plt.ylabel(r'$|\delta\chi|/M_p$') plt.yscale('log') elif(plot_item == 'eps'): plt.xlim([0., 1.]) plt.yscale('log') if(iso_curv): datadir = r'data_isocurv/' else: datadir = r'data/' for i in range(len(kbykps)) : kbykp = kbykps[i] if(kbykp*kp < 1.e3): postfix = str(np.round(kbykp*kp, 4)) else: postfix = "{:.1e}".format(kbykp*kp) if(postfix[-1] == "."): postfix = postfix[:-1] print(postfix) if(path.exists(datadir+r'lna_list_'+postfix+r'.npy')): H_list = np.load(datadir+r'H_list_'+postfix+r'.npy') lna_list = np.load(datadir+r'lna_list_'+postfix+r'.npy') zeta = np.load(datadir+r'zeta_list_'+postfix+r'.npy') Phi = np.load(datadir+r'Phi_list_'+postfix+r'.npy') delta_phi = np.load(datadir+r'delta_phi_'+postfix+r'.npy') delta_chi = np.load(datadir+r'delta_chi_'+postfix+r'.npy') if(plot_item == "eps"): eps_list = np.load(datadir+r'eps_list_'+postfix+r'.npy') t_list = np.load(datadir+r't_list_'+postfix+r'.npy') mp_list = np.load(datadir+r'mp_list_'+postfix+r'.npy') else: H_list, zeta, Phi, delta_phi, delta_chi, lna_list, eps_list, t_list, mp_list = tfm.get_perturb( lnk=tfm.lnaH_pivot+np.log(kbykp), iso_curvature = iso_curv, nsave = 4096) np.save(datadir+r't_list_'+postfix+r'.npy', t_list) np.save(datadir+r'lna_list_'+postfix+r'.npy', lna_list) np.save(datadir+r'H_list_'+postfix+r'.npy', H_list) np.save(datadir+r'zeta_list_'+postfix+r'.npy', zeta) np.save(datadir+r'Phi_list_'+postfix+r'.npy', Phi) np.save(datadir+r'delta_phi_'+postfix+r'.npy', delta_phi) np.save(datadir+r'delta_chi_'+postfix+r'.npy', delta_chi) np.save(datadir+r'eps_list_'+postfix+r'.npy', eps_list) np.save(datadir+r'mp_list_'+postfix+r'.npy', mp_list) if(plot_item == 'zeta'): plt.plot(lna_list, zeta, color=colors[i], ls=lss[i], label=r'$k = '+postfix+r'\,\mathrm{Mpc}^{-1}$') elif(plot_item == 'Phi'): plt.plot(lna_list, Phi, color=colors[i], ls=lss[i], label=r'$k = '+postfix+r'\,\mathrm{Mpc}^{-1}$') elif(plot_item == 'delta_phi'): plt.plot(lna_list, abs(delta_phi)/tfm.Mpl*1.e4, color=colors[i], ls=lss[i], label=r'$k = '+postfix+r'\,\mathrm{Mpc}^{-1}$') elif(plot_item == 'delta_chi'): plt.plot(lna_list, abs(delta_chi)/tfm.Mpl, color=colors[i], ls=lss[i], label=r'$k = '+postfix+r'\,\mathrm{Mpc}^{-1}$') elif(plot_item == 'eps'): plt.plot(lna_list-tfm.num_efolds, mp_list*np.sqrt(H_list/tfm.Mpl), color=colors[0], ls=lss[0], label=r'$\left\vert\frac{\dot\phi\delta\phi + e^{2b}\dot\chi\delta\chi}{(HM_p)^{3/2}}\right\vert$') plt.plot(lna_list-tfm.num_efolds, 2.*eps_list, color=colors[1], ls=lss[1], label=r'$\frac{\dot\phi^2+e^{2b}\dot\chi^2}{H^2M_p^2}$') plt.legend() plt.savefig(plot_item + r'_history.png') plt.show()