The linear Boltzmann equation describes neutron transport in nuclear systems. We consider discretization methods for the time-independent mono-energetic transport equation, and focus on mixed-hybrid primal and dual formulations obtained through an even- and odd-parity flux decomposition. A finite element technique discretizes the spatial variable, and a PN spherical harmonic technique discretizes the angular variable. Mixed-hybrid methods combine attractive features of both mixed and hybrid methods, namely the simultaneous approximation of even- and odd-parity fluxes (thus of flux and current) and the use of Lagrange multipliers to enforce interface regularity constraints. While their study provides insight into purely mixed and purely hybrid methods, mixed-hybrid methods can also be used as such. Mixed and mixed-hybrid methods, so far restricted to diffusion theory, are here generalized to higher order angular approximations. We first adapt existing second-order elliptic mixed-hybrid theory to the lowest-order spherical harmonic approximation, i.e., the P1 approximation. Then, we introduce a mathematical setting and provide well-posedness proofs for the general PN spherical harmonic approximation. Well-posedness theory in the transport case has thus far been restricted to the first-order (integro-differential) form of the transport equation. Proceeding from P1 to PN , the primal/dual distinction related to the spatial variable is supplemented by an even-/odd-order PN distinction in the expansion of the angular variable. The spatial rank condition is supplemented by an angular rank condition satisfied using interface angular expansions corresponding to the Rumyantsev conditions, for which we establish a new derivation using the Wigner coefficients. Demonstration of the practical use of even-order PN approximations is in itself a significant achievement. Our numerical results exhibit an interesting enclosing property when both even- and odd-order PN approximations are employed.