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?

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.