gw.h

// Copyright (C) 2024 Lisa Biermann, Margarete Mühlleitner, Rui Santos, João
// Viana
//
// SPDX-FileCopyrightText: 2024 Lisa Biermann, Margarete Mühlleitner, Rui
// Santos, João Viana
//
// SPDX-License-Identifier: GPL-3.0-or-later
 
#pragma once
 
#include <BSMPT/bounce_solution/bounce_solution.h> // BounceSolution
#include <BSMPT/models/SMparam.h>
#include <BSMPT/utility/Logger.h>
#include <algorithm>
#include <cmath>
#include <functional>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_odeiv2.h>
#include <iostream>
#include <numeric>
#include <vector>
namespace BSMPT
{
 
struct BPLParameters
{
  std::optional<double> Omega_b;
  std::optional<double> f_b;
  std::optional<double> n1;
  std::optional<double> n2;
  std::optional<double> a1;
 
  inline bool IsDefined() const
  {
    return Omega_b.has_value() && f_b.has_value() && n1.has_value() &&
           n2.has_value() && a1.has_value();
  }
};
 
struct DBPLParameters
{
  std::optional<double> Omega_2;
  std::optional<double> f_1;
  std::optional<double> f_2;
  std::optional<double> n1;
  std::optional<double> n2;
  std::optional<double> n3;
  std::optional<double> a1;
  std::optional<double> a2;
 
  inline bool IsDefined() const
  {
    return Omega_2.has_value() && f_1.has_value() && f_2.has_value() &&
           n1.has_value() && n2.has_value() && n3.has_value() &&
           a1.has_value() && a2.has_value();
  }
};
 
struct GravitationalWaveData
{
  bool swON             = true;  // enable sound wave contribution by default
  bool turbON           = true;  // enable turbulence contribution by default
  bool collisionON      = true;  // enable collision contribution by default
  double transitionTemp = false; // transition temperature
  double reheatingTemp  = false; // reheating temperature
  double PTStrength     = false; // strength of EW phase transition
  double betaH          = false; // inverse time scale beta/H
  double vw             = false; // bubble wall velocity
  double vCJ            = false; // Chapman-Jouguet velocity
  double XiShock        = false; // Shock speed
  double Csound_false   = false; // speed sound false vacuum
  double Csound_true    = false; // speed sound true vacuum
  double kappa_col      = false; // efficiency factor for collision
  double kappa_sw       = false; // efficiency factor for sw
  double K_sw           = false; // kinetic energy fraction for sound waves
  double HR             = false; // time scale x max. velocity for sound waves
  int pnlo_scaling      = false; // pressure scaling at NLO, 1 -> N processes
  double Epsilon_Turb   = false; //  fraction of overall kinetic energy in bulk
                                 //  motion that is converted to MHD
  double gstar = false;          // number of eff. d.o.f.
  double FGW0  = false; // Redshift factor for the fractional energy density
  double Hstar0 =
      false; // Reduced Hubble rate at transition temperature redshift today
 
  BPLParameters CollisionParameter;
  DBPLParameters SoundWaveParameter;
  DBPLParameters TurbulanceParameter;
 
  StatusGW status = StatusGW::NotSet; // gw calculation status
};
 
class GravitationalWave
{
private:
public:
  GravitationalWaveData data;
  GravitationalWave(BounceSolution &BACalc,
                    const TransitionTemperature &which_transition_temp =
                        TransitionTemperature::Percolation);
  ~GravitationalWave();
 
  double CalcEpsTurb(double epsturb_in);
 
  const double AbsErr = 0;
 
  const double RelErr = 1e-6;
 
  double h = 0.674;
 
  void CalcPeakCollision();
 
  double CalculateXiShell();
 
  void CalcPeakSoundWave();
 
  void CalcPeakTurbulence();
 
  double BPL(const double &f, const BPLParameters &par) const;
 
  double DBPL(const double &f, const DBPLParameters &par) const;
 
  double CalcGWAmplitude(double f);
 
  double GetSNR(const double fmin, const double fmax, const double T = 3);
 
  friend double snr_integrand(double freq, void *params);
};
 
double SIfunc(const double f);
 
double Rfunc(const double f);
 
double powspec_density(const double f);
 
double h2OmSens(const double f);
 
struct resultErrorPair
Nintegrate_SNR(GravitationalWave &obj, const double fmin, const double fmax);
 
double GetK_sw(const double &alpha, const double &kappa_sw);
 
double GetHstar0(const double &temp, const double &gstar);
 
double GetKtilde(const double &alpha);
 
double GetHtauSH(const double HR, const double K_sw);
 
double GetHtauSW(const double HR, const double K_sw);
 
double GetYpsilon(const double HR, const double K_sw);
 
namespace kappa
{
// Compute kappa_sw https://arxiv.org/abs/2010.09744
 
enum class ExpansionMode
{
  Deflagration,
  Hybrid,
  Detonation
};
 
double mu(double a, double b);
double getwow(double a, double b);
void custom_error_handler(const char *reason,
                          const char *file,
                          int line,
                          int gsl_errno);
std::pair<double, ExpansionMode> getvm(double al, double vw, double cs2b);
int dfdv(double v, const double y[], double dydv[], void *params);
std::vector<std::vector<double>> solve_ode(double vw, double v0, double cs2);
double integrate(const std::vector<double> &y, const std::vector<double> &x);
std::pair<double, double>
getKandWow(double vw,
           double v0,
           double cs2,
           std::vector<std::vector<double>> &vprofile);
double alN(double al, double wow, double cs2b, double cs2s);
double kappaNuMuModel(double cs2b,
                      double cs2s,
                      double al,
                      double vw,
                      std::vector<std::vector<double>> &vprofile);
double kappaNuMuModel(double cs2b, double cs2s, double al, double vw);
 
double Getkappa_col(const double &Tstar,
                    const int &pnlo_scaling,
                    const double &HR,
                    BounceSolution &BACalc);
} // namespace kappa
} // namespace BSMPT