Let $B$ be a quaternion algebra over $\mathbb{Q}$ that is a division algebra over $\mathbb{R}$. The theorem of Hasse-Schilling tells us that the image of the reduced norm of $B^\times$ is $\mathbb{R}^+$.
Let $M$ be a maximal order in $B$. What is known regarding the image of the reduced norm of $M$? To be more specific: let $S$ be the set of primes for which $B$ ramifies. Can we find for every $p \notin S$ an element of reduced norm $p^k$ for some arbitrary odd $k$?
This is a wonderfully rich question! I'll refer to my book (http://quatalg.org).
I'll write $\mathcal{O}$ for your $M$, a definite quaternion order over $\mathbb{Z}$. Let $\mathrm{nrd} \colon \mathcal{O} \to \mathbb{Z}$ be the reduced norm. Then $\mathrm{nrd}$ is a quadratic form on $\mathcal{O}$ (more generally, see Exercise 10.2).
From $B \hookrightarrow B \otimes_{\mathbb{Q}} \mathbb{R} \simeq \mathbb{H}$, we have $\mathbb{H} \simeq \mathbb{R}^4$ as real vector spaces and $\mathrm{nrd}$ on $\mathbb{H}$ in the basis $1,i,j,ij$ is the sum of squares, so positive definite. Now $\mathcal{O} \subseteq B$ sits as a Euclidean lattice (a $\mathbb{Z}$-basis for $\mathcal{O} \simeq \mathbb{Z}^4$ is a $\mathbb{Q}$-basis for $B$), or put another way, $\mathrm{nrd}$ on $\mathcal{O}$ defines a positive definite quaternary quadratic form. Visually, the elements of a given reduced norm are the lattice points on a ball. (See more generally section 17.5, where this is used to prove finiteness of the class set in a "geometry of numbers" manner like for imaginary quadratic fields.)
For a nonnegative integer $n \in \mathbb{Z}_{\geq 0}$, let $r(n) := \#\{\alpha \in \mathcal{O} : \mathrm{nrd}(\alpha)=n\}$ count the number of elements in $\mathcal{O}$ of reduced norm $n$, a finite count by the preceding paragraph. We keep track of these numbers in a generating series $\Theta(q) := \sum_{n=0}^{\infty} r(n)q^n \in \mathbb{Z}[[q]]$. Your question asks: when do we know that a coefficient of this series is nonzero?
It turns out that, owing to its symmetric description (summing over a lattice, which can be compared to summing over the dual lattice using Poisson summation), taking $q := e^{2\pi i z}$ for $z$ in the upper half-plane, the series $\Theta(q)$ is a classical modular form of weight $2$ (and trivial character)! (See section 41.1.)
One can now leverage the theory of modular forms to prove that every sufficiently large positive integer $n \in \mathbb{Z}_{>0}$ that is relatively prime to the discriminant $\mathrm{disc}(\mathcal{O})=\mathrm{disc} \mathrm{nrd}$ (equivalently, $n$ not divisible by any prime in $S=\mathrm{B}$) has $r(n)>0$. In a nutshell, this follows from writing $\Theta$ as a linear combination of Eisenstein series (positive coefficients) and cusp forms (small coefficients) to conclude that the coefficients coprime to the discriminant are eventually positive. And much more! In general, there are subtle issues for integers not coprime, and there is hard work that goes into making the "sufficiently large" effective! I recommend a recent paper by Rouse (https://arxiv.org/pdf/1802.03437.pdf) for more information.
