Can you represent a diagonal density matrix with creation/annihilation operators? Say I’m work with free fermions. In the computational basis one has two basis elements for each site $|0\rangle$ and $|1\rangle$ with creation $d^\dagger$ and annihilation $d$ operators (i.e. $d|1\rangle =|0\rangle$).
Consider the following
$$
d^\dagger d =Id^\dagger dI \equiv(|0\rangle\langle0|+ |1\rangle\langle1|) d^{\dagger}d(|0\rangle\langle0|+ |1\rangle\langle1|)=|1\rangle\langle1|,
$$
similarly, $dd^\dagger \equiv|0\rangle\langle 0|$. So long as the density’s matrix doesn’t have off-diagonal elements one is tempted to
$$
\rho_{\alpha} = (d^\dagger d)^\alpha(dd^\dagger)^{1-\alpha} \equiv \delta_{\alpha,0}|0\rangle \langle0|+\delta_{\alpha,1}|1\rangle \langle1| \quad \alpha \in \{0,1\}.
$$
If I only work with states of the form $|1,0,0,1…\rangle$, i.e. Fock basis elements can I write the corresponding density matrix
$$
|fock\rangle\langle fock| \equiv \rho_{\vec \alpha} = \prod_j \rho_{\alpha_j} =\prod_j  (d_j^\dagger d_j)^{\alpha_j}(d_jd_j^\dagger)^{1-\alpha_j} \quad ?
$$
it doesn’t have any off-diagonal elements so I can’t see why not?
 A: Absolutely! In fact, you can represent any state if you alow for arbitrary combinations of creation and annihilation operators. You already found $|0\rangle\langle 0|=dd^\dagger$ and $|1\rangle\langle 1|=d^\dagger d$; the only remaining components of a single-fermion density matrix are
$$|0\rangle\langle 1|=d\qquad |1\rangle\langle 0|=d^\dagger.$$ This means that a generic single-fermion density matrix can be written as [$0\leq p\leq 1$, $|a|\leq \sqrt{p(1-p)}$]
$$\rho=p dd^\dagger+(1-p)d^\dagger d+ ad+a^* d^\dagger.$$
Your expression for the Fock basis states is correct. If you want to write the most general density matrices, you'll need sums with a whole bunch of indices, because the most general density matrix can be written in a form like
$$\rho=\sum_i\prod_j (\beta_{i,j} d_j d_j^\dagger+\gamma_{i,j} d_j^\dagger d_j +\zeta_{i,j} d_j+\eta_{i,j} d_j^\dagger).$$

This property is harder to show for bosonic creation and annihilation operators $a$ and $a^\dagger$. One can show, for example, that
$$a^{\dagger\,m}a^m=\sum_{l=m}^\infty \frac{l!}{(l-m)!}|l\rangle\langle l|\qquad a^m a^{\dagger\,m}=\sum_{l=0}^\infty \frac{(l+m)!}{l!}|l\rangle\langle l|,$$ so writing a state like $|0\rangle \langle 0|$ requires an infinite sum of the form
$$|0\rangle\langle 0|=\sum_{m=0}^\infty a^{\dagger\,m} a^m\frac{(-1)^m}{m!}.$$
That the linear combinations require infinite sums of creation and annihilation operators makes this a less useful representation. I might have seen states written as $\rho(a,a^\dagger)$ before but I cannot find a useful reference. What I do know is that, through the Segal-Bargmann function, one can create generic states by taking convex combinations of states of the form $$\rho=\sum_i f_i(a^\dagger)|0\rangle\langle 0| f_i^*(a)$$ for analytic functions $f_i$, so this shows that it is always in principle possible to write bosonic states in terms of only creation and annihilation operators.
