How is the number density operator for Dirac fermions defined? In my studies of quantum field theory I stumbled across the question of how to define the number operator for Dirac fermions? I.e. the operator that counts the total number of particle + antiparticles in a given volume (in a certain reference frame). 
Unfortunately I have not found any discussion of this in a textbook (if you happen to know one - please leave me the reference below). 
My current understanding:
My intuition would say that this operator has to commutes with the Hamiltonian for free Dirac fields (after choosing a reference frame), as this would then imply that the total number is conserved for free Dirac fermions. Also, it should not commute with the Hamiltonian of interacting Dirac fermions, due to pair-creation and annihilation. My guess would have been as something proportional to $\int \bar{\psi}\psi  d^3x$, but it turns out this does unfortunately not commute with the free Dirac Hamiltonian $H_D$.
$$H_D = \int \text{d}{\bf{x}} \hspace{0.2cm}\psi^\dagger({\bf{x}}) \gamma^0 [i \gamma^\mu \partial_\mu +m] \psi({\bf{x}}) $$
 A: First, we need to work out the relationship between the field operators $\psi(\mathbf{x})$ and the creation/anninilation operators. That relationship is non-local. Working this out explicitly involves three steps:


*

*Step 1: Diagonalize the Hamiltonian $H$ so that Step 2 is easier.

*Step 2: Choose a Hilbert-space representation of the operator algebra so that the spectrum of $H$ consists of the non-negative real numbers. This makes Step 3 easier.

*Step 3: Determine the number operator by inspection. 
Steps 1 and 2 (and almost step 3) are worked out in section 3.5 of Peskin and Schroeder's book An Introduction to Quantum Field Theory, but that presentation does much more than is needed for answering this question, so I'll outline a more concise version here to help the key ideas shine through more clearly.
In this answer, I'm using the word "particle" to include both particles and antiparticles.
 Step 1: Diagonalize the Hamiltonian 
The Hamiltonian is
$$
 H = \int d^3x\ \psi^\dagger(\mathbf{x})
  \gamma^0[i\gamma^k\nabla_k+m]\psi(\mathbf{x})
 +\text{constant},
\tag{1}
$$
which involves only the space derivatives $\nabla_k$ (not the time derivative $\partial_0$). After expressing the field operators in terms of their Fourier transforms, this becomes
$$
 H = \int \frac{d^3p}{(2\pi)^3}\ 
  \big(\psi(\mathbf{p})\big)^\dagger
  \Omega(\mathbf{p})\psi(\mathbf{p})
 +\text{constant},
\tag{2}
$$
with
$$
 \Omega(\mathbf{p}) := \gamma^0[\gamma^k p_k+m].
\tag{3}
$$
The purpose of using the Fourier transform is that it diagonalizes the Hamiltonian (eliminates the derivative), and the purpose of diagonalizing the Hamiltonian is so that we can determine which combinations of field operators annihilate the vacuum state. This is a prerequisite for determining the number operator, because the number operator should annihilate the vacuum state.
 Step 2: Construct a positive-energy representation 
By definition, the vacuum state is the state of lowest energy. The matrix (3) has both positive and negative eigenvalues. The eigenvalues have magnitude $\omega(\mathbf{p}) := (\mathbf{p}^2+m^2)^{1/2}$, so we can write
$$
 \psi(\mathbf{p}) := \psi^+(\mathbf{p}) + \psi^-(\mathbf{p})
\tag{4}
$$
with $\psi^\pm(\mathbf{p})$ defined by the conditions
$$
 \Omega(\mathbf{p})\psi^\pm(\mathbf{p})
  =\pm\omega(\mathbf{p})\psi^\pm(\mathbf{p}).
\tag{5}
$$
Then, using the anticommutation relations for the field operators (and with the help of a suitable regulator so that everything is well-defined), the Hamiltonian (2) may be written
$$
 H = \int \frac{d^3p}{(2\pi)^3}\ 
  \omega(\mathbf{p})
  \Big(
  \big(\psi^+(\mathbf{p})\big)^\dagger
  \psi^+(\mathbf{p})
  +
  \psi^-(\mathbf{p})
  \big(\psi^-(\mathbf{p})\big)^\dagger
  \Big)
\tag{6}
$$
for a suitable choice of the constant term. Using the field operators, we can generate a Hilbert space representation starting with a state $|0\rangle$ that is annihilated by $\psi^+(\mathbf{p})$ and $\big(\psi^-(\mathbf{p})\big)^\dagger$. These are the annihilation operators for the free Dirac field, and their adjoints are the creation operators. Now we have $H\geq 0$, and the state-vector $|0\rangle$ represents the vacuum state. The purpose of all of the manipulations that led to (6) was to make this last statement obvious.
 Step 3: Determine the number operator by inspection 
The number operator is obtained from (6) simply by omitting the factor $\omega(\mathbf{p})$:
$$
 N = \int \frac{d^3p}{(2\pi)^3}\ 
  \Big(
  \big(\psi^+(\mathbf{p})\big)^\dagger
  \psi^+(\mathbf{p})
  +
  \psi^-(\mathbf{p})
  \big(\psi^-(\mathbf{p})\big)^\dagger
  \Big).
\tag{7}
$$
This clearly annihilates the vacuum state (a prerequisite for any operator that is supposed to count particles) and it clearly commutes with the Hamiltonian (6). Using the anticommutation relations, we can verify that it counts the number of creation operators $\big(\psi^+(\mathbf{p})\big)^\dagger$ and $\psi^-(\mathbf{p})$ that have been applied to the vacuum state. In other words, it counts the number of particles.
 Comments about locality 
The original field operators $\psi(\mathbf{x})$ and $\psi^\dagger(\mathbf{x})$ are local by definition. The creation/annihilation operators (5) are non-local operators: each one involves non-local combination of the original field operators (because of the separation into positive/negative eigenvalues in the Fourier domain). In relativistic QFT, this is unavoidable: a local operator cannot annihilate the vacuum state.
The OP requests an "operator that counts the total number of particle + antiparticles in a given volume", which is not strictly possible in a relativistic QFT: an operator localized in a strictly bounded region cannot annihilate the vacuum state. However, we can construct such an operator for all practical purposes. To do that, simply write the original field operators as
$$
 \psi(\mathbf{x},t) = \psi_C(\mathbf{x},t) + \psi_A(\mathbf{x},t)
\tag{8}
$$
where the subscripts $C$ and $A$ denote the parts that involve creation/annihilation operators, respectively, as defined above. Then the desired operator is
$$
 N(V,t)=\int_V d^3x\ 
  \Big(
  \big(\psi_A(\mathbf{x},t)\big)^\dagger
   \psi_A(\mathbf{x},t)
  +
   \psi_C(\mathbf{x},t)
  \big(\psi_C(\mathbf{x},t)\big)^\dagger
  \Big)
\tag{9}
$$
where the integral is over the desired volume $V$. I added the time argument $t$ because this operator depends on time: particles may enter/leave the given volume as time passes. In other words, the operator (9) doesn't (and shouldn't) commute with the Hamiltonian. The operator $N(V,t)$ is not strictly localized in $V$, because the operators $\psi_{A/C}(\mathbf{x},t)$ are not strictly localized at $\mathbf{x}$: their anticommutators with the original field operators are not zero at non-zero separation (even at equal time). However, the magnitude of the anticommutator falls off exponentially with characteristic length $h/mc\sim 10^{-12}$ meter, so the operator $N(V,t)$ is localized in $V$ for all practical purposes. It counts particles (and antiparticles) in the volume $V$ accurately to the degree that particles can be localized at all in relativistic QFT.
 Comments about interactions 
Beware that this explicit relationship is specific to the free Dirac field. When interactions are included, the operators (5) no longer annihilate the vacuum state (so (7) is no longer the number operator), and we no longer know how to obtain explicit expressions for the operators that do. 
Textbooks focused on the mechanics of perturbation theory sometimes use a misleading language in which the operators constructed above are said to create/annihilate particles even in a model with interactions, but that's sloppy, and it's the root of some widespread confusion about QFT. The definition of "particle" tends to be ambiguous in interacting models, but in flat spacetime, we can at least define it partially by saying that the vacuum state should have zero particles. In a model with interactions, the operators $\psi^+(\mathbf{p})$ and $\big(\psi^-(\mathbf{p})\big)^\dagger$ constructed above don't annihilate the vacuum (lowest-energy) state, so clearly (7) cannot be a particle-counting operator in that case. Equation (7) is the number operator only for the free Dirac field, and we don't have any explicit expression for the number operator in an interacting model like QED. 
