Attractor Solutions

Motivation

Quintessence models can exhibit tracking behaviour at early times, where the equation of state satisfies $\omega_\phi = \omega_{\text{background}}$ and the energy density tracks the energy density of the background, which we assume to be radiation dominated at early times in our implementation. If an attractor exists, we can solve for the initial conditions $\phi_{\text{ini}}$ and $\dot{\phi}_{\text{ini}}$.

Since the attractor solution in CLASS provides only the starting value for the scalar field $\phi_i$ and its derivative $\phi_{\text{i}}^\prime$, we broaden the scope of the attractor solutions slightly by allowing to use an approximate solution as well, namely, if a potential can be approximated by a simple exponential potential, we use the attractor solution of the simple exponential potential for parameter combinations for which the actual attractor solution of the potential does not hold.

Derivation

The tracking condition follows from Copeland et al. (1998, arXiv:gr-qc/9711068). For a tracking attractor, the scalar field energy density is related to the background via:

$$\rho_\phi = \rho_{\text{background}} \frac{3\gamma}{\lambda^2(1 - 3\gamma/\lambda^2)},$$

where $\gamma = 1 + \omega$ is the adiabatic index of the background ($\gamma = 4/3$ for radiation), and $\lambda = -V'/V$ is the dimensionless slope of the potential.

From the attractor condition $\dot{\phi}^2 = 2V(1 + \omega)/(1 - \omega)$, we derive:

$$\phi_{\text{ini}}^\prime = a\sqrt{\rho_{\text{tracking}}(1 + \omega)}.$$

For the field value $\phi_{\text{ini}}$, we solve:

$$2V(\phi) = \rho_{\text{tracking}}(1 - \omega)$$

This is potential-dependent and may not always yield closed-form solutions. In our Mathematica notebook Attractor Solutions.nb we derive the large- and small-field attractors, which we implement or approximate in background.c in CLASS.

Implementation

The implementation is controlled by the flag attractor_ic_scf. When enabled, it computes $\phi_{\text{ini}}$ and $\phi_{\text{ini}}^\prime$ based on the potential type.

A potential can have true large-field or small-field attractor solutions, as well as approximate attractors, e.g. if the potential can be approximated by an exponential function. To assess which regime is dominant during an MCMC, each solution has a tag attractor_regime_scf:

• 0: no attractor (using given initial field values)
• 1: true large-field attractor
• 2: approximate large-field attractor
• 3: true small-field attractor
• 4: approximate small-field attractor

If, for example, the MCMC reveals that the 'best-fit' value for attractor_regime_scf is 2, we'd know that the large-field attractor is preferred, though the conditions for the true large-field attractor solution are not met, yet the initial conditions approximating an exponential attractor solution work well. This can then be used to further investigate the tracking condition, e.g. by plotting $\Omega_\phi/\Omega_m$ or $\left(\Omega_m-\Omega_\phi\right)/\Omega_m$ over the redshift, to see when decoupling takes place.

Relation to Shooting

All implemented dark energy potentials have an overall scaling factor $c_1$, which can be used to in a shooting approach to close the energy budget of the universe. If $c_1$ is used as the shooting target, we face an issue with most of the attractor solutions, as $\phi_{ini}$ depends on $c_1$ in such a way that $d\Omega_\phi/dc_1=0$, i.e. the shooting target becomes flat and shooting fails. For example, for the exponential potential $V = c_1\exp\left(-c_2\phi\right)$, the attractor solution is $\phi = \log(f(\text{params})/c_1)/(-c_2)$, such that $V\left(\phi_{ini}\right) = f\left(\text{params}\right)$ is independent of $c_1$. Therefore, either the slope $c_2$ must be used for shooting, or $\phi_\text{ini}$ is to be sampled, instead of using an attractor solution.

Setup

We first retrieve the potential parameters from the user input and set the background to be radiation dominated.

      double c1 = pba->scf_parameters[0];
      double c2 = pba->scf_parameters[1];
      double c3 = pba->scf_parameters[2];
      double c4 = pba->scf_parameters[3];

      double scf_gamma = 4. / 3.;        // radiation domination
      double omega_background = 1. / 3.; // radiation domination

      double rho_background = rho_rad; // assuming only radiation at initial time
      double rho_tracking = 0.;

Potential-Specific Branches

Each potential has their own attractor solution (or no tracking behaviour).

Power-Law Potential

$$V=c^{4 - p} \phi^p$$

This potential does not have an attractor solution, as

$$\Gamma=\frac{VV^{\prime\prime}}{\left(V^\prime\right)^2}=\frac{p-1}{p}<1.$$

Therefore, the initial conditions have to be provided, instead of using an attractor.

case 1: // power-law
        // pvecback_integration[pba->index_bi_phi_prime_scf] = 0;
        class_stop(pba->error_message, "scf_potential=%d does not have an attractor solution", pba->scf_potential);
        break;

Cosine Potential

$$V=c_1\cos\left(c_2 \phi\right)$$

This potential does not have an attractor solution, as

$$\Gamma=\frac{VV^{\prime\prime}}{\left(V^\prime\right)^2}=-\cot\left(c_2\phi\right)^2\leq0.$$

Therefore, the initial conditions have to be provided, instead of using an attractor.

case 2: // cosine
        class_stop(pba->error_message, "scf_potential=%d does not have an attractor solution", pba->scf_potential);
        break;

Hyperbolic Potential

$$V=c_1\left(1-\tanh\left(c_2 \phi\right)\right)$$

The potential has

$$\Gamma=\frac{VV^{\prime\prime}}{\left(V^\prime\right)^2}=2\cosh\left(c_2\phi\right)\sinh\left(c_2\phi\right)\left(1-\tanh\left(c_2\phi\right)\right)$$

as a general attractor condition. There is no small-field attractor, since $\Gamma \to 0$ for $\phi \to 0$. However, for $c_2 &gt; 0$, there is a large-field attractor as $\Gamma \to 1$ for $\phi \to \infty$. In the large field case, we have $\lambda = -c_2$, $\rho_{\text{tracking}} = \rho_{\text{background}} \frac{3\gamma}{c_2^2(1 - 3\gamma/c_2^2)}$, and $\phi_{\text{ini}}^\prime = a\sqrt{\rho_{\text{tracking}}(1 + \omega)}$, implemented as

      case 3: // hyperbolic
        scf_lambda = -c2;
        /** Reject if no tracking attractor exists (lambda^2 <= 3*gamma implies Omega_scf >= 0.5. CLASS would reject it because radiation is not domiant enough.). */
        class_test(1. <= 3. * scf_gamma / scf_lambda / scf_lambda,
                   pba->error_message,
                   "No tracking attractor exists for potential %d: "
                   "|lambda| = %e < sqrt(3*gamma) = %.4f. "
                   "Use non-tracking ICs (attractor_ic_scf = no) or increase |lambda|.",
                   pba->scf_potential, fabs(scf_lambda), sqrt(3. * scf_gamma));
        {
          rho_tracking = rho_background * 3. * scf_gamma / (scf_lambda * scf_lambda) / (1. - 3. * scf_gamma / scf_lambda / scf_lambda);
          /** The tracking fraction Omega_phi = 3*gamma/lambda^2 is an exact
           *  (non-perturbative) solution of the full Friedmann + Klein-Gordon
           *  system. It does not require subdominance. During RD the tracking
           *  field has w = 1/3 and is correctly counted as radiation in
           *  background_functions (conditional rho_r/rho_m decomposition),
           *  so Omega_r ~ 1 at a_ini regardless of the tracking fraction.
           *  The only physical constraint is lambda^2 > 3*gamma (attractor
           *  existence), already checked above. */
          pvecback_integration[pba->index_bi_phi_prime_scf] = a * sqrt(rho_tracking * scf_gamma);          

The class_test is performed to avoid NaN from a complex root.

Solving $2c_1(1 - \tanh(c_2\phi)) = \rho_{\text{tracking}}(1 - \omega)$ for the initial value yields

$$\phi_{\text{ini}} = \frac{1}{c_2}\text{atanh}\left(\frac{8c_1 c_2^2 - 6c_1\gamma - 3\gamma\rho_{\text{tracking}} + 3\gamma\rho_{\text{tracking}}\omega}{2c_1(4c_2^2 - 3\gamma)}\right)$$

under the condition that the denominator is non-zero and the argument of $\text{atanh}$ lies in $(-1, 1)$. This depends on the sign of $c_1$ and the relative magnitudes of $c_2$, $\gamma$, $\omega$, and $\rho_{\text{tracking}}$. We implemented this as follows:

if (c2 > 0.0 && (
                            /* c1 > 0 branches */
                            (
                                c1 > 0.0 &&
                                scf_gamma > (4.0 / 3.0) * c2 * c2 &&
                                rho_tracking > 0.0 &&
                                omega_background > 1.0 &&
                                omega_background < (1.0 / (3.0 * scf_gamma * rho_tracking)) * (-16.0 * c1 * c2 * c2 + 12.0 * c1 * scf_gamma + 3.0 * scf_gamma * rho_tracking)) ||
                            (c1 > 0.0 &&
                             scf_gamma > 0.0 && scf_gamma < (4.0 / 3.0) * c2 * c2 &&
                             rho_tracking > (1.0 / (3.0 * scf_gamma)) * (16.0 * c1 * c2 * c2 - 12.0 * c1 * scf_gamma) &&
                             omega_background > (1.0 / (3.0 * scf_gamma * rho_tracking)) * (-16.0 * c1 * c2 * c2 + 12.0 * c1 * scf_gamma + 3.0 * scf_gamma * rho_tracking) &&
                             omega_background < 1.0) ||
                            (c1 > 0.0 &&
                             scf_gamma > 0.0 && scf_gamma < (4.0 / 3.0) * c2 * c2 &&
                             rho_tracking > 0.0 && rho_tracking <= (1.0 / (3.0 * scf_gamma)) * (16.0 * c1 * c2 * c2 - 12.0 * c1 * scf_gamma) &&
                             omega_background > 0.0 && omega_background < 1.0) ||
                            /* c1 < 0 branches */
                            (
                                c1 < 0.0 &&
                                scf_gamma > (4.0 / 3.0) * c2 * c2 &&
                                rho_tracking > (1.0 / (3.0 * scf_gamma)) * (16.0 * c1 * c2 * c2 - 12.0 * c1 * scf_gamma) &&
                                omega_background > (1.0 / (3.0 * scf_gamma * rho_tracking)) * (-16.0 * c1 * c2 * c2 + 12.0 * c1 * scf_gamma + 3.0 * scf_gamma * rho_tracking) &&
                                omega_background < 1.0) ||
                            (c1 < 0.0 &&
                             scf_gamma > (4.0 / 3.0) * c2 * c2 &&
                             rho_tracking > 0.0 && rho_tracking <= (1.0 / (3.0 * scf_gamma)) * (16.0 * c1 * c2 * c2 - 12.0 * c1 * scf_gamma) &&
                             omega_background > 0.0 && omega_background < 1.0) ||
                            (c1 < 0.0 &&
                             scf_gamma > 0.0 && scf_gamma < (4.0 / 3.0) * c2 * c2 &&
                             rho_tracking > 0.0 &&
                             omega_background > 1.0 &&
                             omega_background < (1.0 / (3.0 * scf_gamma * rho_tracking)) * (-16.0 * c1 * c2 * c2 + 12.0 * c1 * scf_gamma + 3.0 * scf_gamma * rho_tracking))))
        {
          pvecback_integration[pba->index_bi_phi_scf] = atanh(
                                                            (8.0 * c1 * c2 * c2 - 6.0 * c1 * scf_gamma - 3.0 * scf_gamma * rho_tracking + 3.0 * scf_gamma * rho_tracking * omega_background) / (2.0 * c1 * (4.0 * c2 * c2 - 3.0 * scf_gamma))) /
                                                        c2;
        }

If these conditions are not met, the behaviour can still be approximated by the attractor of the single exponential, which is the dominant term: $$\phi_{\text{ini}} = \frac{1}{c_2}\log\left(\frac{-2(c_1 c_2^2 - 3c_1\gamma)}{3\gamma\rho_{\text{tracking}}(\omega - 1)}\right)$$. This solution also exists only under certain condition and is implemented as follows:

else if (
            /* c1 > 0 branches */
            (
                c1 > 0.0 &&
                scf_gamma > (1.0 / 3.0) * c2 * c2 &&
                rho_tracking > 0.0 &&
                omega_background > 1.0) ||
            (c1 > 0.0 &&
             scf_gamma > 0.0 && scf_gamma < (1.0 / 3.0) * c2 * c2 &&
             rho_tracking > 0.0 &&
             omega_background > 0.0 && omega_background < 1.0) ||
            /* c1 < 0 branches */
            (
                c1 < 0.0 &&
                scf_gamma > (1.0 / 3.0) * c2 * c2 &&
                rho_tracking > 0.0 &&
                omega_background > 0.0 && omega_background < 1.0) ||
            (c1 < 0.0 &&
             scf_gamma > 0.0 && scf_gamma < (1.0 / 3.0) * c2 * c2 &&
             rho_tracking > 0.0 &&
             omega_background > 1.0))
        {
          pvecback_integration[pba->index_bi_phi_scf] =
              log(
                  -2.0 * (c1 * c2 * c2 - 3.0 * c1 * scf_gamma) / (3.0 * scf_gamma * rho_tracking * (-1.0 + omega_background))) /
              c2;
        }

If these conditions are also not met, we still need to assign a value to avoid NaN errors. Since it is only a starting point from which the field then evolves, we consider another attractor solution a save bet: the small-field attractor which exists unconditionally for this potential, $$\phi_{\text{ini}} = \frac{2c_1 c_2^2 - 6c_1\gamma - 3\gamma\rho_{\text{tracking}} + 3\gamma\rho_{\text{tracking}}\omega}{2c_1 c_2^3 - 6c_1 c_2 \gamma}.$$ This last resort is implemented as

else
          pvecback_integration[pba->index_bi_phi_scf] = (2.0 * c1 * c2 * c2 - 6.0 * c1 * scf_gamma - 3.0 * scf_gamma * rho_tracking + 3.0 * scf_gamma * rho_tracking * omega_background) /
                                                        (2.0 * c1 * c2 * c2 * c2 - 6.0 * c1 * c2 * scf_gamma);
        break;

This completes the case for this potential.

Pseudo-Nambu–Goldstone Potential

$$V=c_1^4\left(1+\cos\left(\frac{\phi}{c_2}\right)\right)$$

Since

$$\Gamma=\frac{VV^{\prime\prime}}{\left(V^\prime\right)^2}=-\left(1 + \cos\left(\frac{\phi}{c_2}\right)\right)\cot\left(\frac{\phi}{c_2}\right)\csc\left(\frac{\phi}{c_2}\right)$$

is never 1 or bigger, there is no attractor solution for this potential.

case 4: // pNG
        class_stop(pba->error_message, "scf_potential=%d does not have an attractor solution", pba->scf_potential);
        break;

Inverse Power-Law

While the inverse power-law potential $c^{p + 4} \phi^{-p}$ always has an attractor solution, since

$$\Gamma=\frac{VV^{\prime\prime}}{\left(V^\prime\right)^2}=\frac{p+1}{p}>1,$$

there is no analytic, closed-form solution for $\phi_{\text{ini}}$, i.e. we need to use a shooting algorithm to meet the energy budget constraints. We decided that in this case we can as well just use our MCMC run to find the best-fit values for $\phi_{\text{ini}}$.

case 5: // iPL
        // pvecback_integration[pba->index_bi_phi_prime_scf] = 0.;
        // pvecback_integration[pba->index_bi_phi_scf] = 0.;
        class_stop(pba->error_message, "An attractor solution for the inverse power-law potential exists, however, it requires shooting to match the energy budget. Using MCMC, it's equally efficient to just brute-force it by finding phi_ini and setting phi_prime_ini=0.");
        break;

Exponential Potential

$$V=c_1\exp\left(-c_2 \phi\right)$$ $$\Gamma=\frac{VV^{\prime\prime}}{\left(V^\prime\right)^2}=1.$$

There is always an attractor solution. We have $\lambda = -c_2$ (constant for all $\phi$), $\rho_{\text{tracking}} = \rho_{\text{background}} \frac{3\gamma}{c_2^2(1 - 3\gamma/c_2^2)}$, and $\phi_{\text{ini}}^\prime = a\sqrt{\rho_{\text{tracking}}(1 + \omega)}$. If $1 &lt; 3\gamma/c_2^2$, we used the user provided initial parameters for the scalar field instead a tracking solution. From $$2c_1\exp(-c_2\phi) = \rho_{\text{tracking}}(1 - \omega)$$ we find $$\phi_{\text{ini}} = \frac{1}{c_2}\log\left(\frac{-2(c_1 c_2^2 - 3c_1\gamma)}{3\gamma\rho_{\text{tracking}}(\omega - 1)}\right)$$, provided the argument is positive. If not, we use the small-field attractor solution $$\phi_{\text{ini}} = \frac{2c_1 c_2^2 - 6c_1\gamma - 3\gamma\rho_{\text{tracking}} + 3\gamma\rho_{\text{tracking}}\omega}{2c_1 c_2^3 - 6c_1 c_2 \gamma}$$.

As we explained above, $c_1$ cannot be used for shooting together with the attractor solution, as $\phi_\text{ini}$ depends on $c_1$ in such a way that it cancels in $\Omega_\phi$. We implemented a corresponding class_test that stops the computation if $c_1$ is used for shooting and the attractor solution is requested. The second class_test avoids complex roots.

      case 6: // exponential
        /* KBL: For V = c1*exp(-c2*phi), every attractor IC branch computes
           phi = log(f(params)/c1)/(-c2), so V(phi_ini) = f(params) is
           independent of c1.  Shooting on c1 therefore sees a flat target
           function (dOmega_scf/dc1 = 0) and cannot converge. */
        class_test(pba->scf_tuning_index == 0,
                   pba->error_message,
                   "Exponential potential (case 6) with attractor_ic_scf = yes: "
                   "attractor ICs cancel c1 from V(phi_ini), so shooting on c1 "
                   "(scf_tuning_index = 0) cannot converge. "
                   "Either set attractor_ic_scf = no, or shoot on a different "
                   "parameter (scf_tuning_index != 0).");
        scf_lambda = -c2;
        class_test(1. <= 3. * scf_gamma / scf_lambda / scf_lambda,
                   pba->error_message,
                   "No tracking attractor exists for potential %d: "
                   "|lambda| = %e < sqrt(3*gamma) = %.4f. "
                   "Use non-tracking ICs (attractor_ic_scf = no) or increase |lambda|.",
                   pba->scf_potential, fabs(scf_lambda), sqrt(3. * scf_gamma));
        {
          rho_tracking = rho_background * 3. * scf_gamma / (scf_lambda * scf_lambda) / (1. - 3. * scf_gamma / scf_lambda / scf_lambda);
          pvecback_integration[pba->index_bi_phi_prime_scf] = a * sqrt(rho_tracking * (1. + omega_background));
          if (
              /* c1 > 0 branches */
              (
                  c1 > 0.0 &&
                  scf_gamma > (1.0 / 3.0) * c2 * c2 &&
                  rho_tracking > 0.0 &&
                  omega_background > 1.0) ||
              (c1 > 0.0 &&
               scf_gamma > 0.0 && scf_gamma < (1.0 / 3.0) * c2 * c2 &&
               rho_tracking > 0.0 &&
               omega_background > 0.0 && omega_background < 1.0) ||
              /* c1 < 0 branches */
              (
                  c1 < 0.0 &&
                  scf_gamma > (1.0 / 3.0) * c2 * c2 &&
                  rho_tracking > 0.0 &&
                  omega_background > 0.0 && omega_background < 1.0) ||
              (c1 < 0.0 &&
               scf_gamma > 0.0 && scf_gamma < (1.0 / 3.0) * c2 * c2 &&
               rho_tracking > 0.0 &&
               omega_background > 1.0))
          {
            pvecback_integration[pba->index_bi_phi_scf] =
                log(
                    -2.0 * (c1 * c2 * c2 - 3.0 * c1 * scf_gamma) / (3.0 * scf_gamma * rho_tracking * (-1.0 + omega_background))) /
                c2;
          }
          // If these conditions are not met, we still need to assign some value. The small-field attractor solution avoids NaN.
          else
            pvecback_integration[pba->index_bi_phi_scf] = (2.0 * c1 * c2 * c2 - 6.0 * c1 * scf_gamma - 3.0 * scf_gamma * rho_tracking + 3.0 * scf_gamma * rho_tracking * omega_background) /
                                                          (2.0 * c1 * c2 * c2 * c2 - 6.0 * c1 * c2 * scf_gamma);
        }
        break;

Squared Exponential Potential

$$V=c^{4+p}\phi^{-p}\exp\left(c\phi^2\right)$$

From

$$\Gamma=\frac{VV^{\prime\prime}}{\left(V^\prime\right)^2}=\frac{c^{4+p}(c-4-p)e^{c\phi^2}\phi^{-p} - c^{4+p}(-c+4-p)e^{c\phi^2}(-1-p)p\phi^{-2-p} - 4c^{5+p}e^{c\phi^2}p\phi^{-p}}{(-c^{4+p}e^{c\phi^2}p\phi^{-1-p} + 2c^{5+p}e^{c\phi^2}\phi^{1-p})^2}$$

we can deduce that there is always an attractor solution, since the large-field limit is $\Gamma=1$ and the small-field limit is $\Gamma=1+1/p$. However, to get $\phi_{\text{ini}}$ we need a shooting algorithm. We prefer to rely on MCMC and find the best-fit initial conditions through it instead.

case 7: // SqE
        // pvecback_integration[pba->index_bi_phi_prime_scf] = 0.;
        // pvecback_integration[pba->index_bi_phi_scf] = 0.;
        class_stop(pba->error_message, "An attractor solution for the squared exponential potential exists, however, it requires shooting to match the energy budget. Using MCMC, it's equally efficient to just brute-force it by finding phi_ini phi_prime_ini.");
        break;

Bean's Exponential Potential

$$V=c_1\left(\phi^2 + c_2\right)\exp\left(-c_3\phi\right)$$

If $1 &lt; 3 \gamma / \lambda^2$, there is no physical tracking energy density.

        class_test(1. <= 3. * scf_gamma / scf_lambda / scf_lambda,
                   pba->error_message,
                   "No tracking attractor exists for potential %d: "
                   "|lambda| = %e < sqrt(3*gamma) = %.4f. "
                   "Use non-tracking ICs (attractor_ic_scf = no) or increase |lambda|.",
                   pba->scf_potential, fabs(scf_lambda), sqrt(3. * scf_gamma));

Furthermore, as we approximate the large-field attractor by the exponential potential attractor, $\Omega_\phi$ is independent of $c_1$ if the attractor solution is used, i.e. $c_1$ cannot be used for shooting in this case.

        class_test(pba->scf_tuning_index == 0,
                   pba->error_message,
                   "Bean potential (case 8) with attractor_ic_scf = yes: "
                   "attractor ICs make Omega_scf(today) independent of c1 "
                   "(tracking trajectory depends only on c2, c3, c4). "
                   "Shooting on c1 (scf_tuning_index = 0) cannot converge. "
                   "Use non-tracking ICs (attractor_ic_scf = no) instead.");

There are small- and large-field attractors, as $\Gamma_{\phi\rightarrow0}=\frac{2 c1 + c1 c2 c3^2}{c1 c2 c3^2}$ and $\Gamma_{\phi\rightarrow\infty}=1$. However, no analytic closed-form solution could be derived for the cosmological interesting large-field case, and we therefore approximate it by the dominant exponential term. Nevertheless, a small-field solution was obtained and implemented, which differs from the simple exponential. In the small-field case 'emergency', we use $\phi_{\text{ini}} = \frac{2c_1 c_2 c_3^2 - 6c_1 c_2\gamma - 3\gamma\rho_{\text{tracking}} + 3\gamma\rho_{\text{tracking}}\omega}{2c_1 c_2 c_3^3 - 6c_1 c_2 c_3 \gamma}$ to avoid NaN.

          rho_tracking = rho_background * 3. * scf_gamma / (scf_lambda * scf_lambda) / (1. - 3. * scf_gamma / scf_lambda / scf_lambda);
          pvecback_integration[pba->index_bi_phi_prime_scf] = a * sqrt(rho_tracking * (1. + omega_background));
          // There is no analytic solution for the large-field attractor in this potential. It is approximated by the exponential potential attractor.
          if (
              /* c1 > 0 branches */
              (
                  c1 > 0.0 &&
                  scf_gamma > (1.0 / 3.0) * c3 * c3 &&
                  rho_tracking > 0.0 &&
                  omega_background > 1.0) ||
              (c1 > 0.0 &&
               scf_gamma > 0.0 && scf_gamma < (1.0 / 3.0) * c3 * c3 &&
               rho_tracking > 0.0 &&
               omega_background > 0.0 && omega_background < 1.0) ||
              /* c1 < 0 branches */
              (
                  c1 < 0.0 &&
                  scf_gamma > (1.0 / 3.0) * c3 * c3 &&
                  rho_tracking > 0.0 &&
                  omega_background > 0.0 && omega_background < 1.0) ||
              (c1 < 0.0 &&
               scf_gamma > 0.0 && scf_gamma < (1.0 / 3.0) * c3 * c3 &&
               rho_tracking > 0.0 &&
               omega_background > 1.0))
          {
            pvecback_integration[pba->index_bi_phi_scf] =
                log(
                    -2.0 * (c1 * c3 * c3 - 3.0 * c1 * scf_gamma) / (3.0 * scf_gamma * rho_tracking * (-1.0 + omega_background))) /
                c3;
          }
        }
        // If these conditions are not met, we still need to assign some value. The small-field attractor solution---which exists for this potential and is not equal to the small-field attractor of the exponential---avoids NaN. Phi_prime_ini is identical and does not need to be changed.
        else
          pvecback_integration[pba->index_bi_phi_scf] = (-2.0 * c1 * c2 * c2 * c2 * c3 * c3 + 6.0 * c1 * c2 * c2 * c2 * scf_gamma + 3.0 * c2 * c2 * scf_gamma * rho_tracking - 3.0 * c2 * c2 * scf_gamma * rho_tracking * omega_background) /
                                                        (-2.0 * c1 * c2 * c2 * c2 * c3 * c3 * c3 + 6.0 * c1 * c2 * c2 * c2 * c3 * scf_gamma);
        break;

Double Exponential Potential

$$V=c_1\left[\exp\left(-c_2 \phi\right) + c_3\exp\left(-c_4\phi\right)\right]$$

If $1 &lt; 3 \gamma / \lambda^2$, there is no physical tracking energy density.

        class_test(1. <= 3. * scf_gamma / scf_lambda / scf_lambda,
                   pba->error_message,
                   "No tracking attractor exists for potential %d: "
                   "|lambda| = %e < sqrt(3*gamma) = %.4f. "
                   "Use non-tracking ICs (attractor_ic_scf = no) or increase |lambda|.",
                   pba->scf_potential, fabs(scf_lambda), sqrt(3. * scf_gamma));

Otherwise, the double exponential has early-times as well as late-times attractors: $\Gamma_{\phi\rightarrow0}=\frac{(1 + c3) (c2^2 + c3 c4^2)}{(-c2 - c3 c4)^2}$ and $\Gamma_{\phi\rightarrow\infty}=1$ (if $c2&gt;c4&gt;0$). However, we couldn't derive a closed-form large-field solution. We approximate one using the flatter exponential component ($\lambda = -\min(c_2, c_4)$. In that case, $c_1$ cannot be used for shooting, as $\Omega_\phi$ becomes independent of $c_1$, since $\phi_\text{ini}$ depends on $c_1$. We therefore reject this combination:

        class_test(pba->scf_tuning_index == 0,
                   pba->error_message,
                   "DoubleExp potential (case 9) with attractor_ic_scf = yes: "
                   "attractor ICs make Omega_scf(today) independent of c1 "
                   "(tracking trajectory depends only on c2, c3, c4). "
                   "Shooting on c1 (scf_tuning_index = 0) cannot converge. "
                   "Use non-tracking ICs (attractor_ic_scf = no) instead.");

Besides that, we have $$\rho_{\text{tracking}} = \rho_{\text{background}} \frac{3\gamma}{\lambda^2(1 - 3\gamma/\lambda^2)}$$ and $$\phi_{\text{ini}}^\prime = a\sqrt{\rho_{\text{tracking}}(1 + \omega)}.$$ For the large-field attractor, we stick to the simple exponential attractor solution $$\phi_{\text{ini}} = \frac{1}{\lambda}\log\left(\frac{-2(c_1 \lambda^2 - 3c_1\gamma)}{3\gamma\rho_{\text{tracking}}(\omega - 1)}\right).$$ If the conditions for the large-field simple exponential attractor solution are not met, we use the small-field attractor solution of the double exponential potential itself, for which not only $\phi_{\text{ini}}$ differs from the simple exponential, yet also $\phi_{\text{ini}}^\prime&gt;0$:

$$\phi_{\text{ini}}^\prime = \sqrt{-3\frac{(1 + c_3)\gamma\rho_{\text{tracking}}(1 + \omega)}{c_2 + c_3 c_4 + 3\gamma + 3c_3\gamma}}$$ $$\phi_{\text{ini}} = \frac{c_1 c_2^2 + c_1 c_2^2 c_3 + 2c_1 c_2 c_3 c_4 + 2c_1 c_2 c_3^2 c_4 + c_1 c_3^2 c_4^2 + c_1 c_3^3 c_4^2 - 3c_1\gamma - 9c_1 c_3\gamma - 9c_1 c_3^2\gamma - 3 c_1 c_3^3\gamma - 3\gamma\rho_{\text{tracking}} - 6c_3\gamma\rho_{\text{tracking}} - 3c_3^2\gamma\rho_{\text{tracking}}}{c_1 c_2^3 + 3c_1^2 c_2^2 c_3 c_4 + 3c_1 c_2 c_3^2 c_4^2 + c_1 c_3^3 c_4^3 - 3c_1 c_2\gamma - 6c_1 c_2 c_3\gamma - 3c_1 c_2 c_3^2\gamma - 3c_1 c_3 c_4\gamma - 6c_1 c_3^2 c_4\gamma - 3c_1 c_3^3 c_4\gamma}$$ For $\phi_{\text{ini}}^\prime&lt;0$, we fall back on the small field-attractor solution of the simple exponential potential to avoid NaN.

When we first implemented this, we used $V=c_1\exp\left(-c_2 \phi\right) + c_3\exp\left(-c_4\phi\right)$, yet we realised that for shooting to work, we need a shooting parameter that scales the entire potential, and we re-wrote the potential as $V=c_1\left[\exp\left(-c_2 \phi\right) + c_3\exp\left(-c_4\phi\right)\right]$. To facilitate the implementation of the attractor solution, we simply redefine $c_3 = c_1 \cdot c_3$ and keep the 'old' code otherwise.

case 9: // DoubleExp
        // scf_lambda = -(c1 * c2 + c3 * c4) / (c1 + c3);
        // printf("Double exponential potential selected for attractor ICs\n");
        c3 = c1 * c3; // For the shooting to work, we need an overall factor for the entire potential. To achieve this in the double exponential potential, I choose to put c1 in front of the sum of both exponentials. This changes the meaning of c3, which was initially implemented for the potential c1*exp(c2*phi)+c3*exp(c4*phi). Now it is c1*exp(c2*phi)+ (c1*c3)*exp(c4*phi).
        if (c2 > c4)
          scf_lambda = -c4;
        else
          scf_lambda = -c2;
        /* KBL: DoubleExp attractor ICs + shooting on c1 is non-viable,
           exactly as for the exponential (case 6) and Bean (case 8).

           On the tracking attractor with lambda = min(c2, c4), the shallower
           exponential dominates: V ~ c1 * exp(-lambda*phi).  The attractor IC
           formula computes phi_ini = ln(c1 * f(lambda, gamma)) / lambda, so
           c1 * exp(-lambda * phi_ini) = rho_tracking / g is c1-independent.
           The steeper exponential is subdominant and its coefficient c3 is
           already absorbed into c1*c3 above.  No late-time transition exists:
           the field permanently tracks at Omega_phi = 3*gamma/lambda^2,
           making the full trajectory and Omega_scf(today) c1-independent.

           Unlike Bean (where c4 controls a late-time exit from tracking),
           DoubleExp has no structural feature that breaks the c1-degeneracy.
           Shooting on c2 or c4 (the slopes) is technically feasible but
           physically equivalent to the pure exponential — the tracker has
           w = w_bg at all times and never acts as dark energy.
           Use non-tracking ICs (attractor_ic_scf = no) instead. */
        class_test(pba->scf_tuning_index == 0,
                   pba->error_message,
                   "DoubleExp potential (case 9) with attractor_ic_scf = yes: "
                   "attractor ICs make Omega_scf(today) independent of c1 "
                   "(tracking trajectory depends only on c2, c3, c4). "
                   "Shooting on c1 (scf_tuning_index = 0) cannot converge. "
                   "Use non-tracking ICs (attractor_ic_scf = no) instead.");
        /** Reject if no tracking attractor exists (lambda^2 <= 3*gamma implies Omega_scf >= 0.5). */
        class_test(1. <= 3. * scf_gamma / scf_lambda / scf_lambda,
                   pba->error_message,
                   "No tracking attractor exists for potential %d: "
                   "|lambda| = %e < sqrt(3*gamma) = %.4f. "
                   "Use non-tracking ICs (attractor_ic_scf = no) or increase |lambda|.",
                   pba->scf_potential, fabs(scf_lambda), sqrt(3. * scf_gamma));
        {
          rho_tracking = rho_background * 3. * scf_gamma / (scf_lambda * scf_lambda) / (1. - 3. * scf_gamma / scf_lambda / scf_lambda);
          /** See case 3 (hyperbolic) for why no subdominance check is needed. */
          pvecback_integration[pba->index_bi_phi_prime_scf] = a * sqrt(rho_tracking * scf_gamma);
          // There is no analytic solution for the large-field attractor in this potential. It is approximated by the single exponential potential attractor.
          if (
              /* c1 > 0 branches */
              (
                  c1 > 0.0 &&
                  scf_gamma > (1.0 / 3.0) * scf_lambda * scf_lambda &&
                  rho_tracking > 0.0 &&
                  omega_background > 1.0) ||
              (c1 > 0.0 &&
               scf_gamma > 0.0 && scf_gamma < (1.0 / 3.0) * scf_lambda * scf_lambda &&
               rho_tracking > 0.0 &&
               omega_background > 0.0 && omega_background < 1.0) ||
              /* c1 < 0 branches */
              (
                  c1 < 0.0 &&
                  scf_gamma > (1.0 / 3.0) * scf_lambda * scf_lambda &&
                  rho_tracking > 0.0 &&
                  omega_background > 0.0 && omega_background < 1.0) ||
              (c1 < 0.0 &&
               scf_gamma > 0.0 && scf_gamma < (1.0 / 3.0) * scf_lambda * scf_lambda &&
               rho_tracking > 0.0 &&
               omega_background > 1.0))
          {
            pvecback_integration[pba->index_bi_phi_scf] =
                log(
                    -2.0 * (c1 * scf_lambda * scf_lambda - 3.0 * c1 * scf_gamma) / (3.0 * scf_gamma * rho_tracking * (-1.0 + omega_background))) /
                scf_lambda;
            pba->attractor_regime_scf = 2; // large-field attractor (DoubleExp, approx by single exponential)
            if (pba->background_verbose > 0)
            {
              printf("Large-field attractor (approximated by single exponential) is used. We assign phi_ini_scf=%e and phi_prime_ini_scf=%e.\n", pvecback_integration[pba->index_bi_phi_scf], pvecback_integration[pba->index_bi_phi_prime_scf]);
            }
          }
          // If the conditions are not met, we assign the small field attractor solution of the double exponential to avoid NaN.
          else if (((c1 + c3) * scf_gamma * rho_tracking * (1.0 + omega_background)) /
                       (c1 * c2 + c3 * c4 + 3.0 * c1 * scf_gamma + 3.0 * c3 * scf_gamma) <
                   0.0)
          {
            // The small-field attractor solution has a different phi_prime_ini. This was not the case for Bean's exponential.
            pvecback_integration[pba->index_bi_phi_prime_scf] =
                sqrt(
                    -3.0 * (((c1 + c3) * scf_gamma * rho_tracking * (1.0 + omega_background)) /
                            (c1 * c2 + c3 * c4 + 3.0 * c1 * scf_gamma + 3.0 * c3 * scf_gamma)));
            pvecback_integration[pba->index_bi_phi_scf] = (2.0 * c1 * c1 * c1 * c2 * c2 + 2.0 * c1 * c1 * c2 * c2 * c3 + 4.0 * c1 * c1 * c2 * c3 * c4 + 4.0 * c1 * c2 * c3 * c3 * c4 + 2.0 * c1 * c3 * c3 * c4 * c4 + 2.0 * c3 * c3 * c3 * c4 * c4 - 6.0 * c1 * c1 * c1 * scf_gamma - 18.0 * c1 * c1 * c3 * scf_gamma - 18.0 * c1 * c3 * c3 * scf_gamma - 6.0 * c3 * c3 * c3 * scf_gamma - 3.0 * c1 * c1 * scf_gamma * rho_tracking - 6.0 * c1 * c3 * scf_gamma * rho_tracking - 3.0 * c3 * c3 * scf_gamma * rho_tracking + 3.0 * c1 * c1 * scf_gamma * rho_tracking * omega_background + 6.0 * c1 * c3 * scf_gamma * rho_tracking * omega_background + 3.0 * c3 * c3 * scf_gamma * rho_tracking * omega_background) /
                                                          (2.0 * c1 * c1 * c1 * c2 * c2 * c2 + 6.0 * c1 * c1 * c2 * c2 * c3 * c4 + 6.0 * c1 * c2 * c3 * c3 * c4 * c4 + 2.0 * c3 * c3 * c3 * c4 * c4 * c4 - 6.0 * c1 * c1 * c1 * c2 * scf_gamma - 12.0 * c1 * c1 * c2 * c3 * scf_gamma - 6.0 * c1 * c2 * c3 * c3 * scf_gamma - 6.0 * c1 * c1 * c3 * c4 * scf_gamma - 12.0 * c1 * c3 * c3 * c4 * scf_gamma - 6.0 * c3 * c3 * c3 * c4 * scf_gamma);
            pba->attractor_regime_scf = 3; // approximate small-field attractor (DoubleExp-specific)
            if (pba->background_verbose > 0)
            {
              printf("Small-field attractor (DoubleExp-specific) is used. We assign phi_ini_scf=%e and phi_prime_ini_scf=%e.\n", pvecback_integration[pba->index_bi_phi_scf], pvecback_integration[pba->index_bi_phi_prime_scf]);
            }
          }
          // If these conditions are not met, we still need to assign some value. The small-field attractor solution of the simple exponential avoids NaN. It has again the 'old' phi_prime_ini, which does not need to be changed.
          else
          {
            pvecback_integration[pba->index_bi_phi_scf] = (2.0 * c1 * c2 * c2 - 6.0 * c1 * scf_gamma - 3.0 * scf_gamma * rho_tracking + 3.0 * scf_gamma * rho_tracking * omega_background) /
                                                          (2.0 * c1 * c2 * c2 * c2 - 6.0 * c1 * c2 * scf_gamma);
            pba->attractor_regime_scf = 4; // small-field attractor (simple exponential)
            if (pba->background_verbose > 0)
            {
              printf("Approximate small-field attractor (simple exponential) is used. We assign phi_ini_scf=%e and phi_prime_ini_scf=%e.\n", pvecback_integration[pba->index_bi_phi_scf], pvecback_integration[pba->index_bi_phi_prime_scf]);
            }
          }
        }
        break;

Error Handling

In case we run into numerical errors with the initial conditions, we abort the computation.

    class_test(!isfinite(pvecback_integration[pba->index_bi_phi_scf]) ||
                   !isfinite(pvecback_integration[pba->index_bi_phi_prime_scf]),
               pba->error_message,
               "initial phi = %e phi_prime = %e -> check initial conditions",
               pvecback_integration[pba->index_bi_phi_scf],
               pvecback_integration[pba->index_bi_phi_prime_scf]);