A solution to the linear Boltzmann equation satisfies an energy bound, which reflects a natural fact: The energy of particles in a finite volume is bounded in time by the energy of particles initially occupying the volume augmented by the energy transported into the volume by particles entering the volume over time. In this paper, we present boundary conditions (BCs) for the spherical harmonic $(P_N)$ approximation, which ensure that this fundamental energy bound is satisfied by the $P_N$ approximation. Our BCs are compatible with the characteristic waves of $P_N$ equations and determine the incoming waves uniquely. Both, energy bound and compatibility, are shown on abstract formulations of $P_N$ equations and BCs to isolate the necessary structures and properties. The BCs are derived from a Marshak type formulation of BC and base on a non-classical even/odd-classification of spherical harmonic functions and a stabilization step, which is similar to the truncation of the series expansion in the $P_N$ method. We show that summation by parts (SBP) finite difference on staggered grids in space and the method of simultaneous approximation terms (SAT) allows to maintain the energy bound also on the semi-discrete level.