Combining all these expansions, we obtain up to order $\epsilon^2$:

$$
\begin{aligned}
& \lambda_0 + \epsilon^2 \lambda_1 + c \partial_z U_0(z) + \epsilon^2 c \partial_z U_1(z) + m(z) \\
& = \exp \left( U_1(z_0^*) - 2U_1(\bar{z}) + U_1(z) \right) \\
& \quad + \frac{1 + \frac{\epsilon^2}{2} |\partial_z U_1(\bar{z})|^2 - \frac{3\epsilon^2}{4} \partial_z^2 U_1(\bar{z}) - 2\epsilon^2 U_2(\bar{z}) + \epsilon^2 U_2(z)}{1 + \frac{\epsilon^2}{2} |\partial_z U_1(z_0^*)|^2 - \frac{\epsilon^2}{2} \partial_z^2 U_1(z_0^*) - \epsilon^2 U_2(z_0^*)} \\
& = \exp \left( U_1(z_0^*) - 2U_1(\bar{z}) + U_1(z) \right) \\
& \quad \left( 1 + \epsilon^2 \left( \frac{1}{2} |\partial_z U_1(\bar{z})|^2 - \frac{1}{2} |\partial_z U_1(z_0^*)|^2 + \frac{1}{2} \partial_z^2 U_1(z_0^*) \right. \right. \\
& \quad \left. \left. - \frac{3}{4} \partial_z^2 U_1(\bar{z}) + U_2(z_0^*) - 2U_2(\bar{z}) + U_2(z) \right) \right).
\end{aligned}
$$

By identifying contributions of order $\epsilon^2$ on both sides, we deduce the following equation for the next order correction $U_2$,

$$
\begin{aligned}
& U_2(z_0^*) - 2U_2(\bar{z}) + U_2(z) \\
& = \frac{1}{2} |\partial_z U_1(z_0^*)|^2 - \frac{1}{2} |\partial_z U_1(\bar{z})|^2 + \frac{3}{4} \partial_z^2 U_1(\bar{z}) - \frac{1}{2} \partial_z^2 U_1(z_0^*) + \frac{\lambda_1 + c \partial_z U_1(z)}{1 + \mathbf{G}(z)}.
\end{aligned}
$$

By evaluating, and differentiating at $z = z_0^*$, we deduce the following pair of identities,

$$
\left\{
\begin{array}{l}
0 = \frac{1}{4} \partial_z^2 U_1(z_0^*) + \lambda_1 + c \partial_z U_1(z_0^*) \\
0 = -\frac{1}{2} \partial_z^2 U_1(z_0^*) \partial_z U_1(z_0^*) + \frac{3}{8} \partial_z^3 U_1(z_0^*) + c \partial_z^2 U_1(z_0^*)
\end{array}
\right.
\qquad (\text{F.17})
$$

The second identity enables to compute $p^* = \partial_z U_1(z_0^*)$:

$$
p^* = \frac{3\partial_z^3 U_1(z_0^*)}{4\partial_z^2 U_1(z_0^*)} + 2c = \frac{\partial_z^3 m(z_0^*)}{2\partial_z^2 m(z_0^*)} + 2c, \qquad (\text{F.18})
$$

where $\partial_z^2 U_1(z_0^*)$ and $\partial_z^3 U_1(z_0^*)$ are deduced from equation (F.11) after multiple differentiation, or directly from (F.12). This yields the missing part in (F.15).
SI F.2.5. Analytical expressions of the macroscopic corrections terms $\lambda_1$ and $z_1^*$

Description of Malthus rate $\lambda_1$. The first identity in (F.18) provides $\lambda_1 = -\partial_z^2 U_1(z_0^*)/4 - c\partial_z U_1(z_0^*)$. The expression (F.10) differentiated twice and evaluated at $z = z_0^*$, yields $\partial_z^2 U_1(z_0^*) = 2\partial_z^2 m(z_0^*)$. We conclude from the expression of $p^*$ that

$$ \lambda_1 = -2c^2 - c \frac{\partial_z^3 m(z_0^*)}{2\partial_z^2 m(z_0^*)} - \frac{1}{2} \partial_z^2 m(z_0^*) \quad (F.19) $$

Description of the mean relative phenotype correction $z_1^*$. The first order correction of the mean relative phenotype $z_1^*$ is defined such that $z_0^* + \varepsilon z_1^*$ is the critical point of $U_0 + \varepsilon^2 U_1$, that is $\partial_z(U_0 + \varepsilon U_1)(z_0^* + \varepsilon z_1^*) = 0$. Expanding this relation and keeping only the terms of order $\varepsilon^2$, we obtain using the expression of $p^*$,

$$ z_1^* = -\partial_z U_1(z_0^*) = -\frac{\partial_z^3 m(z_0^*)}{2\partial_z^2 m(z_0^*)} - 2c. \quad (F.20) $$

Description of the local shape. The second derivative of $U_0 + \varepsilon^2 U_1$ at the mean relative phenotype $z^*$ is equal to $\partial_z^2(U_0 + \varepsilon^2 U_1)(z_0^* + \varepsilon^2 z_1^*) = \partial_z^2 U_0(z_0^*) + \varepsilon^2 (\partial_z^3 U_0(z_0^*)z_1^* + \partial_z^2 U_1(z_0^*))$, up to the order $\varepsilon^2$. Since $\partial_z^3 U_0$ is equal to 0, we can deduce from the expression of $U_1$ that the local shape around $z^*$ is given by

$$ \partial_z^2(U_0 + \varepsilon^2 U_1)(z_0^* + \varepsilon^2 z_1^*) = 1 + 2\varepsilon^2 \partial_z^2 m(z_0^*) \quad (F.21) $$

SI G. Numerical computation of the equilibrium ($\lambda, F$)

In order to obtain numerical approximations of the pair $(\lambda, F)$, we get back to the time marching dynamics of the density $f(t, z)$ which satisfies the following equation:

$$ \partial_t f(t, z) - c \partial_z f(t, z) = \beta B(f(t, \cdot))(z) - \mu(z) f(t, z) \quad (G.1) $$

The density $f(t, z)$ is expected to behave like $\exp(\lambda t)F(z)$ for large time. It is preferable to introduce the frequency of traits in population: $p(t, z) = f(t, z)/\int f(t, z') dz'$. The equation for $p$ is:

$$ \begin{aligned} \partial_t p(t, z) + (\beta - \bar{\mu}(t)) p(t, z) - c \partial_z p(t, z) \\ = \beta B(p(t, \cdot))(z) - \mu(z) p(t, z), \quad (G.2) \end{aligned} $$
where the additional $\bar{\mu}(t)$ ensures that $\int p$ remains constant:

$$ \bar{\mu}(t) = \int \mu(z')p(t, z') dz'. \qquad (G.3) $$

We expect that the pair $(\beta(1 - \bar{\mu}(t), p)$ does converge to $(\lambda, F)$ as $t \to +\infty$.

Classical numerical methods were used to approximate (G.2)-(G.3) for large time, until some error threshold is reached for $\|\partial_t p(t, \cdot)\|_{\infty}$. The transport term $-c\partial_z p(t, z)$ was handled using an upwind scheme. The convolutions involved in operator $\mathcal{B}$ were handled using the function `conv` in MATLAB software. The grid mesh was adapted to the scales in SI SI B in order to capture the appropriate phenomena at the correct scale.

## SI H. Comparison with an Individual-based model.

In this section we aim to compare our deterministic approximation with the outcome of a stochastic individual based model (IBM model) with a finite population. We first describe briefly the IBM model. Then we compare the equilibrium distribution of the IBM with our approximation distributions described in Fig. 8. Finally, we compare our results on the effect of the speed of environmental change with the outcomes of the IBM model.

### Stochastic Individual Based Model

We consider a stochastic IBM model where each individual is characterized by its trait $X_i$ and its size $1/N$. They reproduce at a rate $\beta$ and die at a rate that depends on their traits $X_i$, the speed of environmental change $c$ and on the number of individuals $N_t$ in the population at time $t$. More precisely, the individuals may die due to their maladaptness in the phenotypic landscape, which happens at a rate $\mu(X_i - ct)$. Or they may die from density dependence at a rate $\kappa N_t/N$, where $N_t/N$ denotes the size of the population (it is similar to $\rho(t)$ in our deterministic model). The density dependence keeps the number of individuals finite, and proportional to $N$. In particular, when $N$ tends to infinity, the number of individuals also tends to infinity and the (renormalized) stochastic model converges to the deterministic model (2) (Champagnat et al., 2006).

In the case of a birth event, the trait of the offspring is drawn according to the operator $\mathcal{B}$. In the asexual model, the offspring trait $X_{\text{offspring}}$ is given by

$$ X_{\text{offspring}} = X_{\text{parent}} + Y. $$
where $Y$ is a random variable with probability distribution $K_{div}$. In the sexual infinitesimal model, the trait of an offspring $X_{offspring}$ with parents traits $X_{parent,1}$ and $X_{parent,2}$ is given by

$$X_{offspring} = \frac{X_{parent,1} + X_{parent,2}}{2} + Y$$

where $Y$ is drawn from a centered normal distribution with variance $V_{div}/2$.

Numerically, this model has a very high computational cost, especially when the number of individuals is large. As a consequence, we performed the simulations using an approximating model, by first fixing $dt$ to a small but deterministic value. Then, for each individual, we draw a time of birth following the exponential law $\mathcal{E}(\beta)$ and a time of death following the exponential law $\mathcal{E}(\mu(X_i - ct) + \kappa N_t/N)$. Then we simply count which individuals led to a reproduction event and which died on the time-window $[t, t + dt]$. This amounts to the supposition that on this time interval, individuals cannot reproduce more than once.

### Deterministic approximation of the phenotypic distribution

We first compare our approximation of the phenotypic distribution with the empirical distribution of the IBM model for the scenarios described in Fig. 8. When the number of individuals is large (of order $N = 10^4$), we see that our second order approximations are accurate and fit with the empirical distribution of the stochastic model (see Fig. S3).

### Deterministic approximation of the effect of the changing speed

Here, we compare our approximation formula described in Table 2, with the outcomes of the stochastic model with a small number of individuals ($N$ is equal to $10^2$ or $10^3$) in the various scenarios described in Fig. 3.

When the speed of change is slow compared to the critical speeds, our approximations seem accurate in the sense that the approximation error usually falls on our confidence intervals (see Fig. S4-S6). In the infinitesimal sexual model, our approximation also does well when the speed is close to the critical threshold. In this model, we know that the population adapt thanks to the bulk of the population, which moves forward. Thus, even if the number of individuals decreases, many individuals remain at the dominant trait. The number of individuals does not have a critical influence on the adaptation response.
However for the asexual model, when the speed increases, our approximations become less accurate. In this model, only the individuals near the optimal trait help the population to adapt. Thus when the speed increases, the proportion of individuals near the optimal trait decreases because the lag increases. Moreover, when the number of individuals decreases, the actual number of individuals at the optimal trait may be zero, which may lead to an additional burden, and possibly the extinction of the population before the critical value $c_c$ is reached (Calvez et al., 2023). In particular, we see in Figures S4-S6 (a) that the mean fitness of the population drops below 0 for fifty percent of the simulations when the speed is close to the critical speed. Thus the effect of the population size is stronger for the asexual model than for the infinitesimal sexual model.
Asexual model

<figure><img src="image_14.png" /><figcaption>(a) Quadratic selection</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) Super-quadratic selection</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(e) Bounded selection</figcaption></figure>

Infinitesimal sexual model

<figure><img src="image_17.png" /><figcaption>(b) Quadratic selection</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(d) Super-quadratic selection</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(f) Bounded selection</figcaption></figure>

Figure S3: Mutation-selection equilibria F in changing environment with three different shapes of selection: (a)-(b) quadratic function $m(z) = z^2/2$ (blue circled marked curves); (c)-(d) super-quadratic function $m(z) = z^2/2 + z^6/64$ (blue star marked curves); (e)-(f) bounded function $m(z) = m_\infty(1-\exp(-z^2/(2m_\infty)))$ (orange diamond marked curves). The speed of environment change is $c = 0.09$ in the asexual model while it is $c = 0.05$ in the infinitesimal sexual model so that it remains below the critical speeds $c_c$ and $c_{\text{tip}}$ and the distribution deviates significantly from the Gaussian distribution approximation. Other parameters are: $\beta = 1$, $V_{\text{sel}} = 1$ and $V_{\text{div}} = 0.01$ and $m_\infty = 0.5$ in the asexual model and $V_{\text{div}} = 0.1$ and $m_\infty = 1$ in the infinitesimal sexual model. We compare our analytical results (second order results plain marked curves) with the histogram of the stochastic model with $N = 10^4$ individuals. For the asexual scenario, we used the Gaussian kernel.
<figure><img src="image_20.png" /><figcaption>Figure S4: Influence of the speed of environmental change c for a population with finite number of individuals ($N = 10^2$ dashed curves and $N = 10^3$ dash-dotted curves) under quadratic selection $m(z) = z^2/2$. Other parameters are: $\beta = 1$, $V_{\text{sel}} = 1$ and $V_{\text{div}} = 0.01$ in the asexual model and $V_{\text{div}} = 0.1$ in the infinitesimal sexual model. In the asexual model, the mutation kernel is Gaussian. The shade region corresponds to the 95% and 5% confidence intervals around the median. The plain curves correspond to the first order approximation in the asexual model and the second order approximation in the sexual infinitesimal model.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Figure S5: Influence of the speed of environmental change c for a population with finite number of individuals ($N = 10^2$ dashed curves and $N = 10^3$ dash-dotted curves) under super-quadratic selection $m(z) = z^2/2 + z^6/64$. Other parameters are: $\beta = 1$, $V_{\text{sel}} = 1$ and $V_{\text{div}} = 0.01$ in the asexual model and $V_{\text{div}} = 0.1$ in the infinitesimal sexual model. In the asexual model, the mutation kernel is Gaussian. The shade region corresponds to the 95% and 5% confidence intervals around the median. The plain curves correspond to the first order approximation in the asexual model and the second order approximation in the sexual infinitesimal model.</figcaption></figure>
Asexual model

<figure><img src="image_22.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(e)</figcaption></figure>

Infinitesimal sexual model

<figure><img src="image_25.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(f)</figcaption></figure>

Figure S6: Influence of the speed of environmental change **c** for a population with finite number of individuals ($N = 10^2$ dashed curves and $N = 10^3$ dash-dotted curves) under bounded selection function $m(z) = m_\infty(1 - \exp(-z^2/(2m_\infty)))$. Other parameters are: $\beta = 1$, $V_{\text{sel}} = 1$ and $V_{\text{div}} = 0.01$ and $m_\infty = 0.5$ in the asexual model and $V_{\text{div}} = 0.1$ and $m_\infty = 1$ in the infinitesimal sexual model. In the asexual model, the mutation kernel is Gaussian. The shade region corresponds to the 95% and 5% confident intervals around the median. The plain curves correspond to the first order approximation in the asexual model and the second order approximation in the sexual infinitesimal model.