The answer would be "close but not exactly".
What you have here is a Poisson process with rate $\lambda$, stopping after $k$ events.
If you did not stop, then the number of events in each interval of length $1$ would have have a Poisson distribution with mean $\lambda$. But you do stop, and so the distributions are censored by this, slightly reducing their means: this should be most obvious in the final interval observed but actually affecting them all.
Another issue is that the stopping time itself (you call this the overloaded sum[k] and I will call it $S$) is a random variable which has a gamma distribution with mean $\frac{k}{\lambda}$ but you round it up to $T= \lceil S \rceil$ so $T$ has a expected value which is slightly higher than $\frac{k}{\lambda}$.
$\frac{k}{S}$ has a inverse gamma distribution with mean ${\lambda}\frac{k}{k-1}$, so $\frac{k}{T}$ has a expected value which is slightly lower than ${\lambda}\frac{k}{k-1}$ but not necessarily $\lambda$. This $\frac{k}{T}$ is your average number of events per observed interval.
Your simulation gives $\frac{4979}{549}\approx 9.0692$ for mean(num_events) or $\frac kT$, close to $λ=9$ so an absolute difference of $0.0692$. You would get a different difference with a different seed, sometimes positive and sometimes negative before taking the absolute value (never exactly zero as $4979$ is not a multiple of $9$). Here is a simulation with one million cases suggesting that the expected value may be marginally below $\lambda=9$:
set.seed(2023)
k <- 4979
lambda <- 9
cases <- 10^6
T <- ceiling(rgamma(n=cases, shape=k, rate=lambda)) # number intervals
average_events_per_interval <- k/T
mean(average_events_per_interval) # not exactly lambda
# 8.99345
plot.ecdf(average_events_per_interval)
