Can we derive free field expansion formula for the spin-1/2 Dirac field? The Dirac field has the expansion $$\Psi(x)=\int\frac{d^3p}{\sqrt{(2\pi)^32E_p}}\sum\limits_{s=1,2}\Big(b_s(p)u^s(p)e^{-ip\cdot x}+d^\dagger_s(p)v^s(p)e^{+ip\cdot x}\Big)$$ where $b_s$ and $d_s$ are the annihilation operators for the particle and antiparticle respectively with momentum $p$ and spin projection $s$. For a scalar field, such an expansion can be rigorously derived. But I have not seen a derivation of this expansion for $\Psi$; it's written down as if it is very obvious.
Peskin and Schroeder has a derivation but it goes back and forth between Schrodinger and Heisenberg picture while I would like to stick to Heisenberg picture. 
 A: Each component of $\Psi$ satisfies the Klein-Gordon equation, and so we can write (cf. this PSE post)
$$
\Psi_\alpha(x)=\int\frac{d^3p}{\sqrt{(2\pi)^32E_p}}\Big(a_\alpha(p)e^{-ip\cdot x}+b^\dagger_\alpha(p)e^{+ip\cdot x}\Big)
$$
for some operators $a_\alpha,b_\alpha$. If we now require $\Psi$ to satisfy the Dirac equation, we get the algebraic conditions
$$
(\not p+m)a(p)=(\not p-m)b(p)=0
$$
We solve these as follows. Let $u_s(p)\in \mathbb C^4$ with $s=1,2$ be the two linearly linearly independent solutions to $(\not p+m)u(p)=0$, and let $v_s(p)\in \mathbb C^4$ with $s=1,2$ be the two linearly independent solutions to $(\not p-m)v(p)=0$ (there are two and only two solutions because the matrices $\not p\pm m$ have rank 2, as is easily checked). As $u_s,v_s$ are four linearly independent vectors, they are a basis of $\mathbb C^4$, which means we can expand any other vector as linear combinations of them. Thus, we can write
$$
a(p)=\sum_{s=1,2}b_s(p)u_s(p),\qquad b(p)=\sum_{s=1,2} d_s^\dagger v_s(p)
$$
for some scalar operators $b_s,d_s$. Finally, plugging this back into our previous expression, we get
$$
\Psi_\alpha(x)=\int\frac{d^3p}{\sqrt{(2\pi)^32E_p}}\sum_{s=1,2}\Big(b_s(p)u_s(p)e^{-ip\cdot x}+d_s^\dagger v_s(p) e^{+ip\cdot x}\Big)
$$
as required.
For more details see Srednicki §37.
A: Modes decomposition comes from the solution of motion equation. You should start from Dirac equation,
$$[i(\gamma\cdot\partial)-m]\psi=0$$
and consider the following ansatz for $\psi$,
$$\psi=\sum_s\int_{\bf p}\frac{1}{\sqrt{2E_{\bf p}}}\left(b_su_s(p)e^{-ip\cdot x}+d_s^{\dagger}v_s(p)e^{+ip\cdot x}\right).$$
Then you can rewrite your equation in terms of 2$\times$2 block-matrices and consider that $\psi=(\phi,\chi)^T$ (bispinor with two spinors component). Dirac equation gives you system of two equations for $\phi$ and $\chi$. It is convenient to start from $\phi$ and then find $\chi$. Negative frequency solution, $v_s$, can be obtained by charge-congutation matrix $\mathcal{C}=-i\gamma^2\gamma^0$ (it depends on representation of $\gamma$-matrices).
Modes decompoisition arises as you can create/annihilate particle with momentum ${\bf p}$ in different point in space-time. Integration over ${\bf p}$ means that momentum of particle can be arbitrary. Generally, with external fields, spinors $u_s$ and $v_s$ may have more complicated structure.
