diff --git a/tools/sgwb/GW_spectrum.py b/tools/sgwb/GW_spectrum.py new file mode 100644 index 00000000..33c4dbc4 --- /dev/null +++ b/tools/sgwb/GW_spectrum.py @@ -0,0 +1,369 @@ +#!/usr/bin/env python +# coding: utf-8 + +# In[ ]: + + +import json +import os +import shutil +import sys + +try: + import numpy as np + _numpy_available = True +except ImportError: + _numpy_available = False + +try: + from oda_api.json import CustomJSONEncoder +except ImportError: + from json import JSONEncoder as CustomJSONEncoder + +_galaxy_wd = os.getcwd() + + +# In[13]: + + +# we first load necessary packages +import numpy as np +import matplotlib +import matplotlib.pyplot as plt +from matplotlib.colors import LogNorm +get_ipython().run_line_magic('matplotlib', 'inline') +from numpy import log, pi, sqrt +import astropy.constants as const +from astropy.constants import c +import astropy.units as u +import pandas as pd +from random import seed +from random import random +import plotly.graph_objs as go +from scipy.spatial import ConvexHull +from scipy.integrate import quad +from math import * + +from oda_api.data_products import PictureProduct +from oda_api.data_products import ODAAstropyTable + + +# In[14]: + + +# Parameters of the phase transition +# Temperature in GeV +T_star=0.3 #http://odahub.io/ontology#Energy_GeV + +# Numbers of relativistic degrees of freedom +g_star=50 # http://odahub.io/ontology#Integer +#g_star_s=50 # http://odahub.io/ontology#Integer + +# beta/H : rate of the phase transition compared to Hubble rate +beta_H=1. + +# ratio of the energy density deposited in the bubbles to the radiation energy density +alpha=1 # http://odahub.io/ontology#Float + +# fraction of turbulent energy that goes to gw N.B. arXiv:1004.4187 claims that epsilon_turb=0.05, but checks below show that it is rather 0.01 +epsilon_turb=1.0 # http://odahub.io/ontology#Float + +# terminal velocity of bubbles +v_w=0.99 # http://odahub.io/ontology#Float + + +# In[ ]: + + +with open('inputs.json', 'r') as fd: + inp_dic = json.load(fd) +if '_data_product' in inp_dic.keys(): + inp_pdic = inp_dic['_data_product'] +else: + inp_pdic = inp_dic + +for vn, vv in inp_pdic.items(): + if vn != '_selector': + globals()[vn] = type(globals()[vn])(vv) + + +# Theo, please specify the units of $H_0$ in the cell below + +# In[15]: + + +T_star=T_star*u.GeV +T0=2.73 #Temperature of the universe now in Kelvin +conversion=8.617*10**(-9-5) #conversion factor from kelvin to GeV +H0=3.086*10**(22) #current Hubble rate +h=0.7 # dimensionless Hubble Constant +c_s=3**(-0.5) #speed of sound + + +# In[16]: + + +#Eq.20 of arXiv:1512.06239, corrected according to Appendix A of arXiv:1004.4187 +def ka_v(alpha_v,v): + zeta=((2./3.*alpha_v+alpha_v**2)**0.5+(1/3.)**0.5)/(1+alpha_v) + kappa_a=6.9*v**(6./5.)*alpha_v/(1.36-0.037*alpha_v**0.5+alpha_v) + kappa_b=alpha_v**(2./5.)/(0.017+(0.997+alpha_v)**(2./5.)) + kappa_c=alpha_v**0.5/(0.135+(0.98+alpha_v)**0.5) + kappa_d=alpha_v/(0.73+0.083*alpha_v**0.6+alpha_v) + deltak=-0.9*log(alpha_v**0.5/(1+alpha_v**0.5)) + if(vzeta): + return (zeta-1)**3*(zeta/v)**(5./2)*kappa_c*kappa_d/(((zeta-1)**3-(v-1)**3)*zeta**(5./2.)*kappa_c+(v-1)**3*kappa_d) + else: + return kappa_b+(v-c_s)*deltak+(v-c_s)**3/(zeta-c_s)**3*(kappa_c-kappa_b-(zeta-c_s)*deltak) +kappa_v=ka_v(alpha,v_w) +kappa_therm=(1-kappa_v) +print('Bubble wall limiting velocity:',v_w) +print('Fraction of released energy in bulk motion:',kappa_v,', Thermal energy fraction:',kappa_therm) +Omega_star=(kappa_v)*alpha/(1+alpha) +print('Omega in bulk motion:', Omega_star) + + +# In[17]: + + +#Comoving Hubble rate at the phase transition Eq. 11 of arXiv:1512.06239 + +def H_star_bis(T_star): + return 16.5e-6*(T_star/(100.*u.GeV))*(g_star/100)**(1./6.)*u.Hz + +#With factor of 2pi with definition of scale factor and Hubble rate +def H_star(T_star): + scale=T0*conversion*u.GeV/(T_star*(g_star)**(1/3.)) + Hexp=((4.7*10**(-42)*g_star)**(1/2)*(T_star/(conversion*u.GeV))**2)*u.Hz + return scale*Hexp/(2*pi) + + +R_H_star=(const.c/H_star(T_star)).to(u.Mpc) + +print('Reference frequency:',H_star(T_star),', Horizon size:',R_H_star) + +lambda_star=((8*pi)**(1/3)*max([v_w,c_s])/(beta_H*H_star(T_star))*c) # size of bubbles at percolation in units of Hubble scale +#there are different distance scales that we can define + +R_star=lambda_star # size of bubbles at percolation (?) in units of Hubble scale +Delta_w=sqrt((v_w-c_s)**2)/v_w # <1, the width of the bubble wall in units of R_star +R_peak=Delta_w*R_star # the GW spectrum peaks ad the frequency correspoding to the bubble wall width +print('light crossing distance on for the duration of phase transition:',lambda_star, ', bubble size at percolation:',R_star) + + +# In[18]: + + +#Ellis definition of GW_spectrum with broken power law + +def GW_old(f,T_star,alpha,beta_H,v_w): + A=5.1*10**(-2) + a=2.4 + b=2.4 + c=4. + delta=1.0 + scale=T0*conversion*u.GeV/(T_star*(g_star)**(1/3.)) + #Hubble rate from Friedman + Hexp=((4.7*10**(-42)*g_star)**(1/2)*(T_star/(conversion*u.GeV))**2)*u.Hz + fH=Hexp*scale/(2*pi) + Sh=(1+(f/fH)**((-3+a)/delta))**(-delta) + fp=0.7*fH*beta_H + factor=scale**4*(Hexp/(70*10**(3)/(H0)*u.Hz))**2 + return (beta_H)**2*A*(a+b)**c*Sh/(b*(f/fp)**(-a/c)+a*(f/fp)**(b/c))**c*factor + + + +# In[19]: + + +# HL model formula + +def GW_sound2(f,T_star,alpha,beta_H,v_w): + #f=f*u.Hz + R_H_star=(const.c/H_star(T_star)).to(u.Mpc) + kappa_v=ka_v(alpha,v_w) + Omega_star=(kappa_v)*alpha/(1+alpha) + Omega_gw_tilde=10**(-2) #new parameter (that comes from simulations apparently) + lambda_star=((8*pi)**(1/3)*max([v_w,c_s])/(beta_H*H_star(T_star))*c) #characteristic light-crossing distance scale + Delta_w=sqrt((v_w-c_s)**2)/v_w # <1, the width of the bubble wall in units of R_star + s=Delta_w*lambda_star*f/c # the GW spectrum peaks ad the frequency correspoding to the bubble wall width + # frequencies are measured in units of the peak frequency + s1=lambda_star*f/c + r_b=Delta_w # apparently, there are two different notations for the same + m=(9*r_b**4+1)/(r_b**4+1) + M=16*(1+r_b**(-3))**(2/3.)*(r_b*s1)**3/((1+s1**3)**(2./3.)*(3+s**2)**2) + mu=4.78-6.27*r_b+3.34*r_b**2 + B_cursive= Omega_gw_tilde/mu #new parameter + factor=(Omega_star*lambda_star*H_star(T_star)/c)**2/(Omega_star**(1/2.)+lambda_star*H_star(T_star)/c) + old_factor=3*B_cursive*factor*1.64*10**(-5)*(100/g_star)**(1/3.)*v_w/h**2 + return (old_factor*M).cgs + + +# In[20]: + + +# SSM model formula + +def GW_sound1(f,T_star,alpha,beta_H,v_w): + #f=f*u.Hz + R_H_star=(const.c/H_star(T_star)).to(u.Mpc) + kappa_v=ka_v(alpha,v_w) + Omega_star=(kappa_v)*alpha/(1+alpha) + Omega_gw_tilde=10**(-2) #new parameter (that comes from simulations apparently) + lambda_star=((8*pi)**(1/3)*max([v_w,c_s])/(beta_H*H_star(T_star))*c) #characteristic light-crossing distance scale + Delta_w=sqrt((v_w-c_s)**2)/v_w # <1, the width of the bubble wall in units of R_star + s=Delta_w*lambda_star*f/c # the GW spectrum peaks ad the frequency correspoding to the bubble wall width # frequencies are measured in units of the peak frequency + m=(9*Delta_w**4+1)/(Delta_w**4+1) + M=s**9*((Delta_w**4+1)/(Delta_w**4+((s)**4)))**2*(5/(5-m+m*(s)**2))**(5./2.) + mu=4.78-6.27*Delta_w+3.34*Delta_w**2 + B_cursive= Omega_gw_tilde/mu #new parameter + factor=(Omega_star*lambda_star*H_star(T_star)/c)**2/(Omega_star**(1/2.)+lambda_star*H_star(T_star)/c) + old_factor=3*B_cursive*factor*1.64*10**(-5)*(100/g_star)**(1/3.)*v_w/h**2 + return (old_factor*M).cgs + + + + +# In[21]: + + +# Eq 1 of the new Overleaf, from Alberto + +def GW_turb1(f,T_star,alpha,beta_H,v_w,epsilon_turb): + R_H_star=(const.c/H_star(T_star)).to(u.Mpc) + kappa_v=ka_v(alpha,v_w) + Omega_star=(kappa_v)*alpha/(1+alpha) + lambda_star=(((8*pi)**(1/3)*max([v_w,c_s])/(beta_H*H_star(T_star)))*c) #characteristic light-crossing distance scale + + Delta_w=sqrt((v_w-c_s)**2)/v_w # <1, the width of the bubble wall in units of R_star + R_peak=Delta_w*lambda_star # the GW spectrum peaks ad the frequency correspoding to the bubble wall width + u_star=sqrt(0.75*epsilon_turb*Omega_star) # alfven velocity + dt_fin=(2*lambda_star/u_star/c) #eddy processing time + + T_GW=np.log10(1+(H_star(T_star)*dt_fin/(2*pi)))*(f<1/dt_fin)+np.log10(1+(H_star(T_star)/(2*pi*f)))*(f>=1/dt_fin) + T_GW=T_GW**2 + alpha_pi=2.15 + P_pi=(1+((lambda_star*f)/(2.2*c))**alpha_pi)**(-11/(3*alpha_pi)) + M=(lambda_star*f/c)**3*(4*pi**2*T_GW*P_pi)/(0.84*(lambda_star*H_star(T_star)/c)**2) + res=3*(1.75 * 10**(-3))*(Omega_star*epsilon_turb)**2*(lambda_star*H_star(T_star)/c)**2*1.64*10**(-5)*(100/g_star)**(1/3.)/h**2 + return (res*M) + + + +# In[28]: + + +#Range of frequencies used + +f=np.logspace(-10,-2,50)*u.Hz + +lambda_star=((8*pi)**(1/3)*max([v_w,c_s])/(beta_H*H_star(T_star))*c) +kappa_v=ka_v(alpha,v_w) #kappa_therm=(1-kappa_v) +Omega_star=(kappa_v)*alpha/(1+alpha) +u_star=sqrt(0.75*epsilon_turb*Omega_star) # alfven velocity +dt_fin=(2*lambda_star/u_star/c) +Delta_w=sqrt((v_w-c_s)**2)/v_w +delta=1.0 +scale=T0*conversion*u.GeV/(T_star*(g_star**(1/3.))) +Hexp=((4.7*10**(-42)*g_star)**(1/2)*(T_star/(conversion*u.GeV))**2)*u.Hz +fH=Hexp*scale/(2*pi) +fp=0.7*fH*beta_H + + +#Printing the spectrums + +spectrum_turb=GW_turb1(f,T_star,alpha,beta_H,v_w,epsilon_turb) +plt.plot(f,spectrum_turb,label='turbulence',color='blue',linewidth=4,linestyle='dotted') +spectrum_sound=GW_sound2(f,T_star,alpha,beta_H,v_w) + +spectrum_old=GW_old(f,T_star,alpha,beta_H,v_w) +plt.plot(f,spectrum_sound,label='sound waves',color='red',linewidth=4,linestyle='dashed') +plt.plot(f,spectrum_sound+spectrum_turb,color='grey',linewidth=4,label='total') +#plt.plot(f,spectrum_old,color='purple',linewidth=4,label='total') + + +#plt.axvline(1/(dt_fin*u.Hz).cgs,linestyle='dashed',color='black',label=r'$dt_{fin}**(-1)$') +#plt.axvline((c/lambda_star/u.Hz).cgs,linestyle='dotted',color='black',label=r'$c/\lambda_*$') +#plt.axvline((c/(Delta_w*lambda_star)/u.Hz).cgs,linestyle='dashdot',color='black',label=r'$c/(\Delta_w \lambda_*)$') +#plt.axvline(c/(R_peak*u.Hz),linestyle='dashed',color='blue',label='peak') + +plt.legend() + +#Butterfly of compatibility with the NANOgrav data + +butterfly2=np.array([[8.533691739758777e-8,0.0000016780550400704105],[8.361136018023004e-8,1.0740085826829729e-7],[8.499066585218291e-9,3.882820258009008e-9],[1.029481229115365e-9,3.6152599901104956e-11],[1.0353166841868053e-9,2.766052794917924e-10],[7.3944540651210095e-9,6.632930188159042e-9]]) +plt.fill(butterfly2[:,0],butterfly2[:,1]) +plt.xscale('log') +plt.yscale('log') +plt.xlabel('$f$, Hz') +plt.ylabel('$h_0^2\Omega_{gw}(f)$') +plt.savefig('Spectrum.png',format='png',bbox_inches="tight") + + +# In[29]: + + +bin_image = PictureProduct.from_file('Spectrum.png') +from astropy.table import Table +data=[f,spectrum_sound,spectrum_turb] +names=('f[Hz]','Omega_sound_waves[MJD]','Omega_turbulence[MJD]') +spectrum = ODAAstropyTable(Table(data, names = names)) + + +# In[30]: + + +picture = bin_image # http://odahub.io/ontology#ODAPictureProduct +spectrum_astropy_table = spectrum # http://odahub.io/ontology#ODAAstropyTable + + +# In[ ]: + + + + + +# In[ ]: + + +_simple_outs, _oda_outs = [], [] +_galaxy_meta_data = {} +_oda_outs.append(('out_GW_spectrum_picture', 'picture_galaxy.output', picture)) +_oda_outs.append(('out_GW_spectrum_spectrum_astropy_table', 'spectrum_astropy_table_galaxy.output', spectrum_astropy_table)) + +for _outn, _outfn, _outv in _oda_outs: + _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn) + if isinstance(_outv, str) and os.path.isfile(_outv): + shutil.move(_outv, _galaxy_outfile_name) + _galaxy_meta_data[_outn] = {'ext': '_sniff_'} + elif getattr(_outv, "write_fits_file", None): + _outv.write_fits_file(_galaxy_outfile_name) + _galaxy_meta_data[_outn] = {'ext': 'fits'} + elif getattr(_outv, "write_file", None): + _outv.write_file(_galaxy_outfile_name) + _galaxy_meta_data[_outn] = {'ext': '_sniff_'} + else: + with open(_galaxy_outfile_name, 'w') as fd: + json.dump(_outv, fd, cls=CustomJSONEncoder) + _galaxy_meta_data[_outn] = {'ext': 'json'} + +for _outn, _outfn, _outv in _simple_outs: + _galaxy_outfile_name = os.path.join(_galaxy_wd, _outfn) + if isinstance(_outv, str) and os.path.isfile(_outv): + shutil.move(_outv, _galaxy_outfile_name) + _galaxy_meta_data[_outn] = {'ext': '_sniff_'} + elif _numpy_available and isinstance(_outv, np.ndarray): + with open(_galaxy_outfile_name, 'wb') as fd: + np.savez(fd, _outv) + _galaxy_meta_data[_outn] = {'ext': 'npz'} + else: + with open(_galaxy_outfile_name, 'w') as fd: + json.dump(_outv, fd) + _galaxy_meta_data[_outn] = {'ext': 'expression.json'} + +with open(os.path.join(_galaxy_wd, 'galaxy.json'), 'w') as fd: + json.dump(_galaxy_meta_data, fd) +print("*** Job finished successfully ***") + diff --git a/tools/sgwb/sgwb_astro_tool.xml b/tools/sgwb/sgwb_astro_tool.xml new file mode 100644 index 00000000..4dfc4b9c --- /dev/null +++ b/tools/sgwb/sgwb_astro_tool.xml @@ -0,0 +1,63 @@ + + + ipython + numpy + matplotlib + astropy + pandas + plotly + scipy + oda-api + + nbconvert + + ipython '$__tool_directory__/GW_spectrum.py' + + + + + + + + + + + + + + + + + + + + + + + + + + + + + This service provides a calculaiton of the power spectrum of Stochastic +Gravitational Wave Backgorund (SGWB) from a first-order cosmological +phase transition based on the parameterisztions described by `Roper Pol +et al. (2023) <https://arxiv.org/abs/2307.10744>`__. The power spectrum +includes two components: from the sound waves excited by collisions of +bubbles of the new phase and from the turbulence that is induced by +these collisions. + +The cosmological epoch of the phase transition is described by the +temperature, :math:``T_star`` and by the number(s) of relativistic +degrees of freedom, :math:``g_star`` that should be specified as +parameters. + +The phase transition itself is characterised by phenomenological +parameters, :math:``alpha``, :math:``beta_H`` and +:math:``epsilon_turb``, the latent heat, the ratio of the Hubble radius +to the bubble size at percolation and the fraction of the energy otuput +of the phase transition that goes into turbulence. + + \ No newline at end of file