  The "Brownian bees" model describes a system of $N$ independent branching
Brownian particles. At each branching event the particle farthest from the
origin is removed, so that the number of particles remains constant at all
times. Berestycki et al. (2020) proved that, at $N\to \infty$, the
coarse-grained spatial density of this particle system lives in a spherically
symmetric domain and is described by the solution of a free boundary problem
for a deterministic reaction-diffusion equation. Further, they showed that, at
long times, this solution approaches a unique spherically symmetric steady
state with compact support: a sphere which radius $\ell_0$ depends on the
spatial dimension $d$. Here we study fluctuations in this system in the limit
of large $N$ due to the stochastic character of the branching Brownian motion,
and we focus on persistent fluctuations of the swarm size. We evaluate the
probability density $\mathcal{P}(\ell,N,T)$ that the maximum distance of a
particle from the origin remains smaller than a specified value $\ell<\ell_0$,
or larger than a specified value $\ell>\ell_0$, on a time interval $0<t<T$,
where $T$ is very large. We argue that $\mathcal{P}(\ell,N,T)$ exhibits the
large-deviation form $-\ln \mathcal{P} \simeq N T R_d(\ell)$. For all $d$ we
obtain asymptotics of the rate function $R_d(\ell)$ in the regimes $\ell \ll
\ell_0$, $\ell\gg \ell_0$ and $|\ell-\ell_0|\ll \ell_0$. For $d=1$ the whole
rate function can be calculated analytically. We obtain these results by
determining the optimal (most probable) density profile of the swarm,
conditioned on the specified $\ell$, and by arguing that this density profile
is spherically symmetric with its center at the origin.
