-
Notifications
You must be signed in to change notification settings - Fork 0
Scalar Field Potential
Kay Lehnert edited this page Feb 18, 2026
·
2 revisions
We present all potentials with proper formatting in the thesis. Here, we would like to motivate our implementation in background.c. The potentials have up to 4 constants (NaN. To avoid evaluating the potential several times, we decided to do it once, store the value, and then just point to the stored value. This is done once at each time step. Furthermore, the third and fourth derivative are only needed for some of the swampland constraints. It is optional to evaluate those (it is done if the long_info is requested). This is implemented as follows in background.c:
/* Scalar field */
if (pba->has_scf == _TRUE_)
{
phi = pvecback_B[pba->index_bi_phi_scf];
phi_prime = pvecback_B[pba->index_bi_phi_prime_scf];
pvecback[pba->index_bg_phi_scf] = phi; // value of the scalar field phi
pvecback[pba->index_bg_phi_prime_scf] = phi_prime; // value of the scalar field phi derivative wrt conformal time
/* Compute V, V', V'' (and optionally V''', V'''') in a single call.
This avoids redundant extraction of c1–c4 and redundant switch dispatch
that the separate V_scf/dV_p_scf/ddV_scf/d3V_scf/d4V_scf would each do. */
double V_phi, dV_p, ddV_val;
if (return_format == long_info)
{
double d3V_val, d4V_val;
V_scf_derivs(pba, phi, &V_phi, &dV_p, &ddV_val, &d3V_val, &d4V_val);
pvecback[pba->index_bg_d3V_scf] = d3V_val;
pvecback[pba->index_bg_d4V_scf] = d4V_val;
}
else
{
V_scf_derivs(pba, phi, &V_phi, &dV_p, &ddV_val, NULL, NULL);
}
double KE = phi_prime * phi_prime / (2 * a2); // kinetic energy term (reuse a2)
pvecback[pba->index_bg_V_scf] = V_phi;
pvecback[pba->index_bg_dV_p_scf] = dV_p; // pure derivative (no coupling) — cached for swampland diagnostics
pvecback[pba->index_bg_dV_scf] = dV_p + coupling_scf(pba, rho_cdm_prime(pba, phi, pvecback), pvecback); // V'_eff = V'_pure + coupling
pvecback[pba->index_bg_ddV_scf] = ddV_val;
/* Check for NaN/Inf in scalar field potential and its derivatives.
This catches extreme parameter combinations (e.g. very large phi * c2
in the DoubleExp potential) before they propagate further. */
class_test(isnan(V_phi) || isinf(V_phi) || isnan(dV_p) || isinf(dV_p) || isnan(ddV_val) || isinf(ddV_val),
pba->error_message,
"NaN or Inf in scalar field potential at phi=%e (V=%e, dV=%e, ddV=%e). "
"The scalar field parameter combination is likely unphysical.",
phi, V_phi, dV_p, ddV_val);
class_test(isnan(pvecback[pba->index_bg_dV_scf]) || isinf(pvecback[pba->index_bg_dV_scf]),
pba->error_message,
"NaN or Inf in effective scalar field derivative dV_eff at phi=%e. "
"The coupling or scalar field parameters may be unphysical.",
phi);
pvecback[pba->index_bg_rho_scf] = (KE + V_phi) / 3.; // energy of the scalar field
pvecback[pba->index_bg_p_scf] = (KE - V_phi) / 3.; // pressure of the scalar field
rho_tot += pvecback[pba->index_bg_rho_scf];
p_tot += pvecback[pba->index_bg_p_scf];
dp_dloga += 0.0; /** <-- This depends on a_prime_over_a, so we cannot add it now! */
// KBL: Do NOT add scalar field to rho_m when it acts as dark energy (quintessence).
// The original CLASS code below splits scf into relativistic and non-relativistic parts,
// but this is incorrect for tracking quintessence models where the scf IS dark energy.
// Otherwise, Omega_scf is added 4 times to Omega_m, which yields roughly Omega_m=0.7*4=2.8 in MCMC as mean value. This is probably because DE is then dominating the entire energy budget.
// rho_r += 3. * pvecback[pba->index_bg_p_scf]; // field pressure contributes radiation
// rho_m += pvecback[pba->index_bg_rho_scf] - 3. * pvecback[pba->index_bg_p_scf]; // the rest contributes matter
// printf(" a= %e, Omega_scf = %f, \n ",a, pvecback[pba->index_bg_rho_scf]/rho_tot );
}The potentials themselves are implemented in a single function, to avoid retrieving the value of
/**
* Compute the scalar field potential and all its derivatives up to
* the requested order in a single function call. This avoids redundant
* extraction of c1–c4, redundant switch dispatch, and allows sharing
* common subexpressions (exp, pow, trig) across derivatives.
*
* @param pba Input: pointer to background structure
* @param phi Input: scalar field value
* @param V Output: V(phi)
* @param dV Output: V'(phi) (pure, without coupling)
* @param ddV Output: V''(phi)
* @param d3V Output: V'''(phi) (only written if d3V != NULL)
* @param d4V Output: V''''(phi) (only written if d4V != NULL)
*/
void V_scf_derivs(
struct background *pba,
double phi,
double *V,
double *dV,
double *ddV,
double *d3V,
double *d4V)
{
double c1 = pba->scf_parameters[0];
double c2 = pba->scf_parameters[1];
double c3 = pba->scf_parameters[2];
double c4 = pba->scf_parameters[3];
switch (pba->scf_potential)
{
case 1: /* power-law: V = c1^(4-c2) * phi^c2 + c3 */
{
double pref = pow(c1, 4. - c2);
double phic2 = pow(phi, c2);
double phic2m1 = phic2 / phi; /* phi^(c2-1) */
double phic2m2 = phic2m1 / phi; /* phi^(c2-2) */
*V = pref * phic2 + c3;
*dV = c2 * pref * phic2m1;
*ddV = c2 * (c2 - 1.) * pref * phic2m2;
if (d3V)
*d3V = (c2 - 2.) * (c2 - 1.) * c2 * pref * phic2m2 / phi;
if (d4V)
*d4V = (c2 - 3.) * (c2 - 2.) * (c2 - 1.) * c2 * pref * phic2m2 / (phi * phi);
break;
}
case 2: /* cosine: V = c1 * cos(c2*phi) */
{
double s = sin(c2 * phi);
double c = cos(c2 * phi);
*V = c1 * c;
*dV = -c1 * c2 * s;
*ddV = -c1 * c2 * c2 * c;
if (d3V)
*d3V = c1 * c2 * c2 * c2 * s;
if (d4V)
*d4V = c1 * c2 * c2 * c2 * c2 * c;
break;
}
case 3: /* hyperbolic: V = c1 * [1 - tanh(c2*phi)] */
{
double t = tanh(c2 * phi);
double ch = cosh(c2 * phi);
double sech2 = 1.0 / (ch * ch);
*V = c1 * (1. - t);
*dV = -c1 * c2 * sech2;
*ddV = 2. * c1 * c2 * c2 * t * sech2;
if (d3V)
*d3V = -2. * c1 * c2 * c2 * c2 * (cosh(2. * c2 * phi) - 2.) * sech2 * sech2;
if (d4V)
{
double ch5 = ch * ch * ch * ch * ch;
*d4V = 2. * c1 * c2 * c2 * c2 * c2 * (sinh(3. * c2 * phi) - 11. * sinh(c2 * phi)) / ch5;
}
break;
}
case 4: /* pNG: V = c1^4 * [1 + cos(phi/c2)] */
{
double c1_4 = pow(c1, 4.);
double invc2 = 1.0 / c2;
double s = sin(phi * invc2);
double c = cos(phi * invc2);
*V = c1_4 * (1. + c);
*dV = -c1_4 * invc2 * s;
*ddV = -c1_4 * invc2 * invc2 * c;
if (d3V)
*d3V = c1_4 * invc2 * invc2 * invc2 * s;
if (d4V)
*d4V = c1_4 * invc2 * invc2 * invc2 * invc2 * c;
break;
}
case 5: /* iPL: V = c1^(4+c2) * phi^(-c2) */
{
double pref = pow(c1, 4. + c2);
double phimc2 = pow(phi, -c2);
double inv_phi = 1.0 / phi;
*V = pref * phimc2;
*dV = -c2 * pref * phimc2 * inv_phi;
*ddV = c2 * (c2 + 1.) * pref * phimc2 * inv_phi * inv_phi;
if (d3V)
*d3V = -c2 * (c2 + 1.) * (c2 + 2.) * pref * phimc2 * inv_phi * inv_phi * inv_phi;
if (d4V)
*d4V = c2 * (c2 + 1.) * (c2 + 2.) * (c2 + 3.) * pref * phimc2 * inv_phi * inv_phi * inv_phi * inv_phi;
break;
}
case 6: /* exponential: V = c1 * exp(-c2*phi) */
{
double e = c1 * exp(-c2 * phi); /* = V */
*V = e;
*dV = -c2 * e;
*ddV = c2 * c2 * e;
if (d3V)
*d3V = -c2 * c2 * c2 * e;
if (d4V)
*d4V = c2 * c2 * c2 * c2 * e;
break;
}
case 7: /* SqE: V = c1^(c2+4) * phi^(-c2) * exp(c1*phi^2) */
{
double pref = pow(c1, c2 + 4.);
double phimc2 = pow(phi, -c2);
double ep = exp(c1 * phi * phi);
double base = pref * phimc2 * ep;
double phi2 = phi * phi;
double inv_phi = 1.0 / phi;
*V = base;
*dV = base * inv_phi * (2. * c1 * phi2 - c2);
double phi4 = phi2 * phi2;
*ddV = base * inv_phi * inv_phi * (4. * c1 * c1 * phi4 + (-4. * c1 * c2 - 2. * c1) * phi2 + c2 * (c2 + 1.));
if (d3V)
{
double phi6 = phi4 * phi2;
*d3V = base * inv_phi * inv_phi * inv_phi * (8. * c1 * c1 * c1 * phi6 + (-12. * c1 * c1 * (c2 - 1) * phi4) + 6. * c1 * c2 * c2 * phi2 - c2 * (c2 + 1.) * (c2 + 2.));
}
if (d4V)
{
double phi4v = phi2 * phi2;
double phi6 = phi4v * phi2;
double phi8 = phi6 * phi2;
*d4V = base * inv_phi * inv_phi * inv_phi * inv_phi * (16. * c1 * c1 * c1 * c1 * phi8 + (16. * c1 * c1 * c1 * phi6 * (3. - 2. * c2)) + (12. * c1 * c1 * phi4v * (2. * (c2 - 1.) * c2 + 1.)) - 4 * c1 * phi2 * c2 * (c2 + 1.) * (2. * c2 + 1.) + c2 * (c2 + 1.) * (c2 + 2.) * (c2 + 3.));
}
break;
}
case 8: /* Bean: V = c1 * [(c4-phi)^2 + c2] * exp(-c3*phi) */
{
double e = exp(-c3 * phi);
double dp = c4 - phi;
double dp2 = dp * dp;
double bracket = dp2 + c2;
*V = c1 * bracket * e;
*dV = -c1 * e * (2. * dp + c3 * bracket);
*ddV = c1 * e * (2. - 4. * dp * c3 + c3 * c3 * bracket);
if (d3V)
*d3V = -c1 * c3 * e * (c3 * (c3 * (c2 + phi * phi) - 6. * phi) + 6.);
if (d4V)
*d4V = c1 * c3 * c3 * e * (c3 * (c3 * (c2 + phi * phi) - 8. * phi) + 12.);
break;
}
case 9: /* DoubleExp: V = c1 * (exp(-c2*phi) + c3 * exp(-c4*phi)) */
{
double e2 = exp(-c2 * phi);
double e4 = exp(-c4 * phi);
*V = c1 * (e2 + c3 * e4);
*dV = -c1 * (c2 * e2 + c3 * c4 * e4);
*ddV = c1 * (c2 * c2 * e2 + c3 * c4 * c4 * e4);
if (d3V)
*d3V = -c1 * (c2 * c2 * c2 * e2 + c3 * c4 * c4 * c4 * e4);
if (d4V)
*d4V = c1 * (c2 * c2 * c2 * c2 * e2 + c3 * c4 * c4 * c4 * c4 * e4);
break;
}
default:
*V = 0.;
*dV = 0.;
*ddV = 0.;
if (d3V)
*d3V = 0.;
if (d4V)
*d4V = 0.;
break;
}
}