Klein Gordon Equation in Quintessence Models

[TimeFall]

Hello! I'm studying various dark energy models, and as a part of the project, I need to be able to numerically solve the Klein-Gordon (KG) equation and the Friedmann Equation (FE) in the context of a canonical scalar field. I wasn't sure whether or not this belonged here or in the computational section, so forgive me if I chose wrong. Anways, the KG equation is:
$$ \ddot{\phi} + 3H\dot{\phi} + V_{,\phi} = 0 $$

And the FE is:
$$ H^2 = \frac{8\pi G}{3} \left [ \frac{1}{2} \dot{\phi}^2 + V(\phi) + \rho_M \right ] $$

Where $H$ is the Hubble parameter, $\phi$ is the scalar field (both functions of time), a dot represents a derivative with respect to time, and $V_{,\phi} = \frac{dV(\phi)}{d\phi}$. $V(\phi)$ is the potential of the scalar field. $G$ is simply Newton's gravitational constant. $\rho_M$ is the energy density of the other matter components in the Universe (both relativistic and non-relativistic). There is also an equation that governs the evolution of $H$:

$$ \dot{H} = -\frac{8\pi G}{2} \left ( \dot{\phi}^2 + \rho_M + P_M \right ) $$

Where $P_M$ is the general relativistic pressure corresponding to $\rho_M$. Additionally, for the energy density I used:

$$ \rho_M = \rho_m + \rho_r = \rho_{m,0}a^{-3} + \rho_{r,0}a^{-4} $$

Where $a$ is the scale factor. The pressure is:

$$ P_M = \frac{\rho_r}{3} = \frac{\rho_{r,0} a^{-4}}{3} $$

What I have tried doing is using the gsl's ode (ordinary differential equations) library to solve this system. I wrote the KG equation as a first order equation by introducing the variable $\theta \equiv \dot{\phi}$ and then changed variables from $t$ to $a$. This gave me three equations in the system that needed to be solved: the KG equation, $\dot{\phi} = \theta$ and the $\dot{H}$ equation. However, when I run my code the gsl appears not to update $H$ despite updating both $\theta$ and $\phi$. $H$ is simply stuck at the initial value that I gave it.

I was just wondering if anyone knew what was wrong with this approach and/or if there is a better way to do it? I've read about the slow roll parameters and how they can be used to make simplifications to the KG and FE equations, but I wasn't sure when in the time evolution this approximation was valid. Thank you very much!

 

[Chalnoth]

Sounds like a coding error. I'd need to see the code to have an idea of what's going on. Note that you can make code look decent in the forum by enclosing it in "code" tags (see the help page).

 

[TimeFall]

Thanks! Here's the function where I actually do the integration (I've left out the #includes and variable declarations in the interest of space)

double y[3] = {Par.theta, Par.phi, Par.H_init};

   // Set up ode system and driver
   gsl_odeiv2_system sys = {system_def, jacobian, 3};
   gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&sys, 
                          gsl_odeiv2_step_rk8pd, Par.init_stepsize, 
                          Par.epsabs, Par.epsrel);

   // Set up starting integration limits
   Nstart = -1.0 * Par.start_logz * log(10.0);
   Nend   = -1.0 * Par.end_logz * log(10.0);

   Nold = Nstart;
   Nnew = Nstart;

   for(i = 1; i <= (STEPS + 1); i++)
   {
      status = gsl_odeiv2_driver_apply(d, &Nold, Nnew, y);

      if(status != GSL_SUCCESS)
      {
         printf("There was an error during integration!\n");
         exit(EXIT_FAILURE);
      }

      // If we're here, we're good so far. Get H, a, and z
      Par.H[i - 1] = y[2];
     
      Par.a[i - 1] = exp(Nnew);
      Par.z[i - 1] = (1.0 / Par.a[i - 1]) - 1.0;

      // Update limits
      Nold = Nnew;
      Nnew = Nstart + (((Nend - Nstart) * (double)i) / (double)STEPS);
   }

   // Clean up
   gsl_odeiv2_driver_free(d);

The variable N (Nstart, Nend, Nold, Nnew) is defined as N = ln(a), where a is the scale factor.
Then, there's the actual system of odes that need to be solved (mentioned in the original post. The first equation is for theta, the second is for phi, and the third is for H), again, not showing variable declarations and such in the interest of brevity:

int system_def(double N, const double y[], double f[], void *parameters)
{
   // System of odes to solve. N is the variable the derivatives are taken
   // with respect to, y[] are the dependent variables (theta, phi, H),
   // f are the derivs (d theta /da, etc.)

      rho_M = (3.0 * pow(H0_gev, 2.0) / (8.0 * M_PI * G * exp(3.0 * N))) * (OMEGAM0 + (OMEGAR0 / exp(N)));

      P_M = pow(H0_gev, 2.0) * OMEGAR0 / (8.0 * M_PI * G * exp(4.0 * N));

      // Theta eq
      f[0] = -1.0 * ((3.0 * y[0]) - ((n * pow(M, 4.0 + n) * pow(y[1], -1.0 * (n + 1.0))) / y[2]));

      // Phi eq
      f[1] = y[0] / y[2];

      // H equation
      f[2] = (-1.0 * 8.0 * M_PI * G / 2.0) * (pow(y[0], 2.0) + rho_M + P_M);

   return GSL_SUCCESS;
}

Where OMEGAM0 and OMEGAR0 are the nonrelativistic and relativistic density parameters, respectively (0.3 and 8e-5). H0_gev is the current value of the Hubble parameter in GeV (this code uses units of c = hbar = k_b = 1). G is also given in GeV. If there's anything else you need to see, just let me know. Thanks!

 

[Chalnoth]

I see you're updating the values of H, a, and z ever iteration of the loop, but I don't see you making use of those values in subsequent iterations of the loop. Is this for debugging purposes, or is it a bug?

Anyway, the first thing I'd do to try and debug this is set a breakpoints in the system_def function and step through it one step at a time, making sure it's doing what I expect. Alternatively, you could just add some output statements.

The other question I have is: do you define the jacobian? Presumably you can just pass a null pointer to the function if the solver you're using doesn't use the jacobian.

 

[TimeFall]

Hello! Thanks for responding. I'm not using H, a, and z in this loop, no. They are saved and used later in the code. I just needed this function to get their values. The way I understood the gsl to work is that it updated the values of y[] and worked with those, and so I was saving their values from each step. I'll try the breakpoints, for sure. I didn't even think to pass a null pointer, to be honest. I just defined it as an empty function that simply returns gsl_success. I've tried this set up with the simple harmonic oscillator system of odes just to make sure I knew how to use it and that it gave correct results, which it did.

 

[Chalnoth]

On the jacobian: I see. It should probably be fine, but passing a null pointer is the safest thing here. That way you don't run the risk of the code expecting the function to do something (because the code will crash if any piece of the code tries to run the function).

Getting this kind of mathematical thing can be a bit tricky. One other way to track down errors in the math is to create test functions which verify each mathematical formula individually. For instance, you could have tests which consider the analytically-solvable cases such as only having matter, radiation, or a cosmological constant as components. You'd expect some degree of error, but it'd be a start.

 

[TimeFall]

Awesome, thanks! I'll give that a whirl and see what happens.