I am studying the evolution of a population in time, see figure below. I simulate $N=100$ different trajectories from time $t=0$ to $t=400$. When I compute the average of these trajectories at every time step I get the trajectory in represented in blue, and the point at which it crosses the criticality threshold of $0.5$ is shown in the figure as the vertical dashed black line.
I am especially interested in the time where the trajectories cross the criticality threshold of $0.5$ and want to deal with the uncertainty at which time the trajectory will cross the threshold.
How can I compute a PDF to represent such uncertainty over time? Is there any specific approach that anyone can recommend to do this ?
I want to use the empirical results, e.g. how many trajectories out of $N$ have crossed at times $t=0:500$
I assume that your simulation is based on discrete time steps $t \in \{0,...,T\}$ (in your case $T = 400$).
Lets assume, that you run $n \in \mathbb{N}_{\geq 1}$ simulations, which gives you trajectories $s_1,...,s_n$, where each $$ s_i : \{0,...,T\} \rightarrow \mathbb{R} $$ describes the value of the $i$-th simulation in dependence of time. So $s_i(t)$ would be the value of the $i$-th simulation at time $t$. Each $s_i$ describes exactly one of the grey lines in your plot. Now, for each $s_i$, you can compute a critical time point $$t_{\text{crit}}(s_i) := \min\{ t \in \{0,..,T\} \ \vert \ s_i(t) \geq c_{\text{crit}} \} .$$ In your case, we have $c_{\text{crit}} = 0.5$.
We can see the values given by the expressions $t_{\text{crit}}(s_i)$ as $n$ independent realizations of a discrete random variable $$T_{\text{crit}} : \Omega \rightarrow \mathbb{N}_{\geq 0} .$$ If I understand your question correctly, then you are interested in estimating the probability distribution of $T_{\text{crit}}$.
Given the data you have obtained from computing the $t_{\text{crit}}(s_i)$, you can empirically estimate the probability mass function of $T_{\text{crit}}$ via $$ \Pr(T_{\text{crit}} = k) = \frac{\text{Number of elements in the set } \{ i \in \{1,...,n\} \ \vert \ t_{\text{crit}}(s_i) = k \}}{n} $$ for $k \in \mathbb{N}_{\geq 0}$.
To get a reasonable estimate with this method, take care to choose $n$ as big as possible given your computational resources. Also, you have to choose $T$ big enough, so that all you $s_i$ (or at least almost all of them, you can ignore a few) actually exceed $c_{\text{crit}}$ at some time in $\{0,...,T\}$, because otherwise some of the expressions $t_{\text{crit}}(s_i)$ are the minimum of an empty set and do not exist.
In your plot, clearly your have to increase $T>>400$, in order to follow the method. If you cannot increase $T = 400$, you have to additionally consider the case that the critical value is never hit, but this only requires minor changes to the procedure I just described.