$-
\tau_B$ and $\tau_B$ we match the solutions at times $\mp\tau_B-\epsilon$ and $\mp\tau_B+\epsilon$, neglecting any evolution in the intervening time period. This is similar to what was done in Refs. [@Elisa; @Subodh]. A second justification of this method (in addition to the one given earlier) is that for modes we are interested in, the time interval $2\epsilon$ is so small, thus the mode functions do not have enough time to oscillate inside the Hubble radius. On the other hand, in our first model (with a flat plateau for $a(t)$, model 1) we consider the case that the interval $2\Delta$ is sufficiently long such that some of the modes we are interested in have time to oscillate while the mode is inside the Hubble radius. Very large scale modes, on the other hand, still do not oscillate inside the Hubble radius.

### Between region I and region II

The matching conditions are $$\begin{aligned}
v_1(-\tau_B-\epsilon) = v_2 (-\tau_B+\epsilon), \quad
v'_1(-\tau_B-\epsilon) = v'_2 (-\tau_B+\epsilon) \, .\end{aligned}$$ We write down the equations explicitly in terms of coefficients $c_{ij}$ of the fundamental solutions of the equation of motion. The index $i$ stands for the region, the index $j$ (either $1$ or $2$) running over the 