How do I correctly evaluate this simple expression with delta-commutator operators? Let $a$ denote a bosonic annihilation operator for the mode $q$. The different-$q$ commutator is
$$ [a(q),a^{\dagger}(q')] = \delta(q-q') $$
Then, how do I evaluate this expression correctly using the commutator:
$$ \left\langle0\right|a(q)a^{\dagger}(q)\left|0\right\rangle $$
which should be equal to $1$? Applying the commutator leads to
$$ \left\langle0\right|a(q)a^{\dagger}(q)\left|0\right\rangle = \delta(0)\left\langle0\middle|0\right\rangle = \delta(0) $$
which clearly is not a well-defined expression.
How do I correctly evaluate this simple expression with delta-commutator operators? What I mean is, what is the rigorous way of / the explanation for properly calculating this expectation value given the operators and their commutator?
 A: Beware of long post with lots of formal atrocities. A formal explanation would require notions of functional analysis and and rigged Hilbert spaces and would probably be better suited for Math.SE. I'm going to give you an intuitive picture for why this happens. Also, regardless of Cosmas' comment reagrding differing normalizations, something like this will always happen. The choice of normalization is just a choice of when it happens.
$a^\dagger_i|\emptyset\rangle$ is, in general, not a normalizable state. This usually shows up in continuous state spaces, but the reasoning is as follows.
Any physical state is never just created by a single mode $i$, but by a linear combination of many of them. A single particle state $|\psi\rangle$ can usually be written as a linear combination involving coefficients $\psi_i$
$$|\psi\rangle=\sum_i\mu_i\psi_ia^\dagger_i|\emptyset\rangle,$$
where $\mu_i$ is some positive function of $i$ that is independent of our $|\psi\rangle$ and is what we call a measure on the set of modes. This measure can be read off the canonical commutation relation, since in more general cases it may read
$$[a_i,a^\dagger_j]=\frac{\delta_{ij}}{\sqrt{\mu_i\mu_j}}.$$
The normalization of the state $|\psi\rangle$ is then imposed,
$$1=\langle\psi|\psi\rangle=\sum_{ij}\mu_i\mu_j\psi_i^*\psi_j\langle\emptyset|a_ia^\dagger_j|\emptyset\rangle=\sum_{ij}\sqrt{\mu_i\mu_j}\psi_i^*\psi_j\delta_{ij}\langle\emptyset|\emptyset\rangle$$
$$\sum_i\mu_i|\psi_i|^2=1,$$
and $p_i=\mu_i|\psi_i|^2$ gives us the probability distribution on the set of modes. For continuous spaces, it's impossible to choose a nonzero $\mu$ that gives a finite value for the "size" of our space,
$$V=\sum_i\mu_i,$$
so we end up taking $\mu\rightarrow0^+$, or informally, $\mu=dx$ and replacing every sum with an integral, and end up with a Dirac delta:
$$\frac{\delta_{xx'}}{\sqrt{dxdx'}}"="\delta(0)\delta_{xx'}=\delta(x-x'),$$
where we have intuitively "identified" $\delta(0)=1/dx$, which matches a lot of the properties of the Delta function.
All in all, your result is correct. It's just that the normalizations involved to create a continuous set of modes stop the pure $a^\dagger(q)|\emptyset\rangle$ from being a normalizable, physically meaningful state. In fact, it's clear from the previous discussion that the closest thing one could get to a state comprised purely of mode $q$ is
$$|q\rangle=\frac{1}{\sqrt{\delta(0)}}a^\dagger(q)|\emptyset\rangle,$$
whatever it could mean. The way to make sense of this is pretty simple. Let's do a couple of manipulations,
$$|q\rangle=\frac{1}{\sqrt{\delta(0)}}\int dq'\delta(q-q')a^\dagger(q')|\emptyset\rangle=\int dq'\frac{\delta(q-q')}{\sqrt{\delta(0)}}a^\dagger(q')|\emptyset\rangle$$
$$|q\rangle=\int dq'\sqrt{\delta(q-q')}a^\dagger(q')|\emptyset\rangle,$$
and now we choose some representation of the Delta as a limit of some other function. Could be the Gaussian, for example:
$$\delta(x)=\lim_{\epsilon\rightarrow0^+}\frac{e^{-x^2/\epsilon}}{\sqrt{\pi\epsilon}},$$
and write the delta like that. We're goint to get, for once, something that actually kind of makes sense (up to commuting limits that don't actually commute):
$$|q\rangle=\lim_{\epsilon\rightarrow0^+}\int dq'\frac{e^{-(q-q')^2/2\epsilon}}{(\pi\epsilon)^{1/4}}a^\dagger(q')|\emptyset\rangle=\lim_{\epsilon\rightarrow0^+}|q_\epsilon\rangle,$$
where the $|q_\epsilon\rangle$ are perfectly well defined states with norm 1, being states centered at $q$ with size $\sqrt\epsilon$. This limit, the physical state $|q\rangle$, does not actually exist, however. Saying it does is equivalent to nonsense such as putting $\sqrt{\delta(0)}$ as a real number in a denominator for something. If you wanted a result for a very specific value of $q$, you would just carry out all the calculations with states of the form $|q_\epsilon\rangle$, and once you've reached your result, take the limit $\epsilon\rightarrow0^+$, which is usually well defined enough to make sense.
In practice other regulators for the delta function are used, such as finite volume for spacetime and momentum space. All of these hinge on the idea that the pure mode $a^\dagger_i|\emptyset\rangle$ is not a physical and normalizable state.
