Motivated by applications to the study of ultracold atomic gases near the unitarity limit, we investigate the structure of the operator product expansion (OPE) in non-relativistic conformal field theories (NRCFTs). The main tool used in our analysis is the representation theory of charged (i.e. non-zero particle number) operators in the NRCFT, in particular the mapping between operators and states in a non-relativistic radial quantization Hilbert space. Our results include: a determination of the OPE coefficients of descendant operators in terms of those of the underlying primary state, a demonstration of convergence of the (imaginary time) OPE in certain kinematic limits, and an estimate of the decay rate of the OPE tail inside matrix elements which, as in relativistic CFTs, depends exponentially on operator dimensions. To illustrate our results we consider several examples, including a strongly interacting field theory of bosons tuned to the unitarity limit, as well as a class of holographic models. Given the similarity with known statements about the OPE in SO(2,d) invariant field theories, our results suggest the existence of a bootstrap approach to constraining NRCFTs, with applications to bound state spectra and interactions. We briefly comment on a possible implementation of this non-relativistic conformal bootstrap program.