Where this polarization vector is coming from? When dealing with vector fields in his QFT book, Schwartz writes the classical field in terms of a basis which I don't know how he is getting.
He first introduces the Proca Lagrangian
$$\mathcal{L}=-\dfrac{1}{4}F_{\mu\nu}F^{\mu\nu}+\dfrac{1}{2}m^2A_\mu A^\mu.$$
The equations of motion are then $(\Box + m^2)A_\mu = 0$ and $\partial^\mu A_\mu = 0$.
He then says:

Let us now find explicit solutions to the equations of motion. We start by Fourier transforming our classical fields. Since $(\Box+m^2)A_\mu=0$, we can write any solution as
$$A_\mu(x)=\sum_i \int \dfrac{d^3 p}{(2\pi)^3}\tilde{a}_i(p)\epsilon^i_\mu(p) e^{ipx}, \quad p_0=\omega_p=\sqrt{|\mathbf{p}|^2+m^2}$$
for some basis vectors $\epsilon^i_\mu(p)$. For example, we could trivially take $i=1,2,3,4$ and use four vectors $\epsilon^i_\mu(p)=\delta^i_\mu$ in this decomposition. Instead, we want a basis that forces $A_\mu$ to automatically satisfy also its equation of motion $\partial^\mu A_\mu=0$. This will happen if $p^\mu \epsilon_\mu^i(p) =0$. For any fixed $4$-momentum $p^\mu$ with $p^2=m^2$, there are three independent solutions to this equation given by three $4$-vectors $\epsilon_\mu^i(p)$, necessarily $p^\mu$ dependent, which we call polarization vectors. Thus we only have to sum over $i=1,2,3$. We conventionally normalize the polarizations by $\epsilon_\mu^\ast \epsilon^\mu =-1$.

I don't understand this at all. Taking the Fourier transform of $\Box A_\mu + m^2 A_\mu = 0$ and denoting $\hat{A}_\mu$ the Fourier transform we have
$$A_\mu(x)=\int \dfrac{d^3 p}{(2\pi)^3} (a_\mu(p) e^{-ipx}+a_\mu^\ast(p)e^{ipx}).$$
This is done by taking the three-dimensional Fourier transform, getting the equation $\partial_t^2\hat{A}+\omega_p^2 \hat{A}=0$, realizing that $\hat{A}_\mu= a_\mu(p) e^{-i\omega_p t}+b_\mu(p) e^{i\omega_p t}$ and finally using the condition that $A_\mu$ is real so that $b_\mu(p)=a_\mu^\ast(-p)$ which leads directly to the formula above.
There is no $\epsilon^i_\mu(p)$ anywhere. Nor can I see why there should be. The sum has only two terms anyway.
So I'm really missing something here. How is this result properly derived?
 A: You're almost there. Given
$$
A_\mu(x)=\int \dfrac{d^3 p}{(2\pi)^3} (a_\mu(p) e^{-ipx}+a_\mu^\dagger(p)e^{ipx}).
$$
you may pick any set of four linearly independent vectors $\{\epsilon_\mu^1,\epsilon_\mu^2,\epsilon_\mu^3,\epsilon_\mu^4\}$ and expand the $a_\mu$ in terms of them:
$$
a_\mu=\sum_{i=1}^4a_i \epsilon_\mu^i
$$
for some coefficients $a_i$. The "trivial basis" mentioned by S. is $\epsilon^i_\mu\equiv\delta^i_\mu$, in which case this expression becomes
$$
a_\mu=\sum_{i=1}^4a_i \delta_\mu^i=a_\mu
$$
i.e. the coefficients $a_i$ are just the cartesian components of $a_\mu$. In principle, this is a valid basis, but we can do better. For one thing, in this basis the transversality condition $p\cdot a=0$ is not immediate to implement.
If we choose a basis that changes with $p$, such that
$$
p\cdot\epsilon^1=p\cdot\epsilon^2=p\cdot\epsilon^3=0,\qquad \epsilon^4=p/m
$$
then we may write, as before,
$$
a_\mu=\sum_{i=1}^4a_i \epsilon_\mu^i
$$
but now the condition $p\cdot a$ is equivalent to $a_4=0$, so that in effect
$$
a_\mu=\sum_{i=1}^3a_i \epsilon_\mu^i
$$
With this,
$$
A_\mu(x)=\sum_{i=1}^3\int \dfrac{d^3 p}{(2\pi)^3} (a_i \epsilon_\mu^i e^{-ipx}+a^\dagger_i \epsilon_\mu^{i*} e^{ipx}).
$$
which is the final expression for a free Proca field.
