How to adapt the documentation example for curved universes?

[alexninogh]

Problem

I would like to reproduce the plot of the primordial power spectrum (fully numerically integrated) shown in the first row of Figure 14 in your paper (https://arxiv.org/abs/2205.07374). Is it possible to have an example script that does this?

Suggestion

A valuable addition to the documentation could be an explanation of the recommended functions to use based on the input we wish to provide. Furthermore, it would be helpful to include suggestions on the range of parameters that could aid in achieving convergence of the numerical solution to the differential equation.

What I have tried so far

I've tried following your paper and thesis to understand how to correctly set up the initial conditions using the functions you provide in PrimPy, but I've failed to understand the link between the code and the necessary operations.

I have taken the following parameters from you paper:

Omega_K0h2 = -0.005
h = 0.7
A_s = 2e-9
N_star = 60
f_i = 5

I have placed these parameters in the example code from the documentation:

import numpy as np
import matplotlib.pyplot as plt

from primpy.parameters import K_STAR
import primpy.potentials as pp
from primpy.events import UntilNEvent, InflationEvent, CollapseEvent
from primpy.initialconditions import InflationStartIC, ISIC_NsOk, ISIC_Nt
from primpy.time.inflation import InflationEquationsT as InflationEquations
from primpy.solver import solve
from primpy.oscode_solver import solve_oscode

t_eval = np.logspace(5, 8, 2000)
K = 1
N_star = 60      # number of e-folds of inflation after horizon crossing
N_end = 70       # end time/size after inflation, arbitrary in flat universe
delta_N_reh = 2  # extra e-folds after end of inflation to see reheating oscillations
A_s = 2e-9       # amplitude of primordial power spectrum at pivot scale
Pot = pp.StarobinskyPotential

Omega_K0h2 = -0.005
h = 0.7
Omega_K0 = Omega_K0h2 / h**2
f_i = 5
Omega_Ki = f_i * Omega_K0

Lambda, phi_star, N_star = Pot.sr_As2Lambda(A_s=A_s, N_star=N_star, phi_star=None)
pot = Pot(Lambda=Lambda)
eq = InflationEquations(K=K, potential=pot, track_eta=False)
ev = [UntilNEvent(eq, value=N_end+delta_N_reh),  # decides stopping criterion
      InflationEvent(eq, +1, terminal=False),    # records inflation start
      InflationEvent(eq, -1, terminal=False)]    # records inflation end

ic_fore = ISIC_NsOk(equations=eq, N_star=N_star, Omega_Ki=Omega_Ki, Omega_K0=Omega_K0, h=h, 
                    phi_i_bracket=[phi_star-1, phi_star+1], t_i=t_eval[0])

fward = solve(ic=ic_fore, events=ev, t_eval=t_eval)
print("fward._N_beg = ", fward._N_beg)
print("fward._N_end = ", fward._N_end)

I seem to have an issue with the initial conditions.
I changed the initial conditions input from ISIC_Nt to ISIC_NsOk, because I think this is what is needed to reproduce figure 14, right?
If I understand correctly, _N_beg and _N_end are 'nan' because the solver is unable to solve the equations given the input I'm providing.

[lukashergt]

Hi Alessandro, thanks for opening this issue.

I would like to reproduce the plot of the primordial power spectrum (fully numerically integrated) shown in the first row of Figure 14 in your paper (https://arxiv.org/abs/2205.07374). Is it possible to have an example script that does this?

The code to produce the figures is actually available on zenodo:

It probably won't be possible to use the plotting code therein as is, given that there is quite a gap from primpy version 2.3.6 to 2.9.7, but it should be useful nonetheless.

---

I changed the initial conditions input from ISIC_Nt to ISIC_NsOk, because I think this is what is needed to reproduce figure 14, right?

Correct.

If I understand correctly, _N_beg and _N_end are 'nan' because the solver is unable to solve the equations given the input I'm providing.

Actually, the solver works fine, it is just stopping too early, because you inadvertently told it to... To understand what is going on, I recommend checking plots of the inflationary background, e.g.:

fig, ax = plt.subplots()
ax.plot(fward._N, fward.phi)
ax.axvline(N_end+delta_N_reh, c='k', ls='--')
ax.set_ylim(0, 7)
ax.set_ylabel(r"$\phi$")
ax.set_xlabel(r"$N=\ln(a)$")

The problem, here, is that you changed the example from N_star=50 to N_star=60, while telling the solver to not go beyond _N=72, such that you never reach the end of inflation.

To make things work, just change the ending criterion for the solver. Currently the ending criterion is defined by UntilNEvent(eq, value=N_end+delta_N_reh) with N_end=70 and delta_N_reh=2. For figure 14, you don't need to evolve the solver into reheating, so the simplest is to tell the solver to stop at the end of inflation, by changing ev to:

ev = [InflationEvent(eq, +1, terminal=False),   # records inflation start
      InflationEvent(eq, -1, terminal=True)]    # records inflation end and stops

With that change the solver goes through to the end of inflation:

Now you can follow the rest of the example script to get you to the full numeric primordial power spectrum:

[alexninogh]

Thank you! Everything works perfectly now. Moreover, the Jupyter notebooks on Zenodo were very helpful.