Proving this completeness relation for polarization of photon? 
In the Coulomb gauge, I want to prove the completeness relation $$\sum_{s=1,2} \epsilon_i^{(s)} (\epsilon_j^{(s)} )^{*} = \delta_{ij} - \hat{p}_i \hat{p}_j$$ for a photon traveling in the $z$-direction and $\vec{\epsilon}^{(1)} = (1,0,0)$ and $\vec{\epsilon}^{(2)} = (0,1,0)$. The vector $\hat{p}_i$ is a unit vector pointing in the direction of travel. 

I see this relation in every textbook on QED, but nobody proves it.
Since the above equation is a matrix equality, I wanted to write it in terms of matrices. I let $\hat{p}_i = (0,0,1)$. Also, I think I can assume that $i = j$, otherwise the sum is zero because of orthogonality of the polarization vectors (is this correct?). So I wrote the sum as $$\sum_{s=1,2} \epsilon_i^{(s)} (\epsilon_j^{(s)} )^{*} = \epsilon_i^{(1)} (\epsilon_i^{(1)} )^{*} + \epsilon_i^{(2)} (\epsilon_i^{(2)} )^{*} = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \end{pmatrix} + \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix} \begin{pmatrix} 0 & 1 & 0 \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix}. $$ Then the RHS of the above completeness relation is then (if I assume $i = j$) : $$ \mathbb{1}_{3 \times 3} - \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \begin{pmatrix} 0 & 0 & 1 \end{pmatrix} \mathbb{1}_{3 \times 3} = \mathbb{1}_{3 \times 3} - \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}  = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{pmatrix}. $$ 
Is this how you are suppose to prove this relation, and is my reasoning correct? Thanks in advance.
 A: Just an addition to the comment section. 
Instead the matrix-based proof You should try to use simpler approach.
You define the polarization 3-vectors $\epsilon(\mathbf p)$ of the photon with 3-momentum $\mathbf p$ to be orthogonal to the $\mathbf p$. This means that there are two such vectors. Let's denote the number of these vectors as $(a)$. Let's also denote the i-th component of the a-th polarization vector as $\epsilon_{i}^{(a)}$. Finally, let's also norm it to one. 
Because of definitions, You have two relations valid for arbitrary direction of $\mathbf p$:
$$
\tag 1 \epsilon^{(1,2)}_{i}(p)p_{i} = 0, \quad (\mathbf{\epsilon}^{(1,2)}(p))^{2} = 1
$$
valid for arbitrary direction 3-vector $p_{i}$. Let's expand the sum $\sum_{l}\epsilon^{l}_{i}(p)\epsilon^{*l}_{j}(p)$ on possible tensors being $g_{ij}, p_{i}p_{j}$:
$$
\tag 2 \sum_{l}\epsilon^{l}_{i}(p)\epsilon^{*l}_{j}(p) = Ag_{ij} + Bp_{i}p_{j}
$$
Then by using $(1)$ you'll obtain that $A = 1, B = -\frac{1}{\mathbf{p}^{2}}$.
The remark about decompositions
Let's understand how to guess the expression $(2)$.   
First, note that the vectors 
$$
\{\mathbf e^{(k)}, k = 1,2,3\} \equiv \{\mathbf e^{(1)}(\mathbf p), \mathbf e^{(2)}(\mathbf p), \frac{\mathbf p}{|\mathbf p|}\}
$$ 
forms the orthonormal basis in 3-dimensional space of vectors $\mathbf p$. These vectors have obvious transformation law under the rotations generated by the matrix $\hat{R}$:
$$
e^{(k)}_{i} \to \hat{R}_{i}^{\ j}e^{(k)}_{j}
$$
Next, note that the $\epsilon^{a}_{i}\epsilon^{b}_{j}$ is in fact the tensor products of vectors, i.e.,
$$
\epsilon^{(a)}_{i}\epsilon^{(b)}_{j} \equiv (\epsilon^{(a)}\otimes \epsilon^{(b)})_{ij}
$$
Therefore it can be expanded on matrices which have definite transformation laws under rotations with some coefficients $A, B, ...$ (people call this "helicity decomposition", or something like that; in general it is called "the decomposition of the reducible representation of the given group on irreducible representations"). 
Finally, note that the sum $P_{ij} = \sum_{a = 1,2}\epsilon^{(a)}_{i}\epsilon^{(a)}_{j}$ is explicitly invariant under the rotations around $\hat{\mathbf p}$ axis. I general it can also contain the part which is invariant under any rotation.
Therefore $P^{ij}$ is expanded on following matrices: $g_{ij}$, which isn't transformed under the rotations $R$, and $p_{i}p_{j}$, which isn't changed under the rotations along $\hat{\mathbf p}$. We don't include terms like $\mathbf{e}^{(1,2)}\otimes \mathbf{e}^{(1,2)}$ or $\mathbf{e}^{(1,2)}\otimes \hat{\mathbf p}$ since they are not invariant under the rotation around $\hat{\mathbf p}$.
Therefore we obtain $(2)$.
A: The answer given above has already mentioned the necessary ingredients for you to find the relation you are after. Here I will reiterate with allegedly simpler manner:
Coming from D.Tong's QFT Lecture Notes, page-130.
Define the polarization vectors $\vec{\epsilon}_1(\vec{p}), \vec{\epsilon}_2(\vec{p}) $ as
$$
\vec{\epsilon}_r(\vec{p}) := \epsilon^i_r(\vec{p}) \ ;
$$
$$
 r=1,2. \ , \ i=1,2,3.
$$
with the properties:
$$
\vec{\epsilon}_r(\vec{p}) \cdot \vec{p}  =  \epsilon^i_r(\vec{p}) p^i =0 
$$
$$
\vec{\epsilon}_r(\vec{p}) \cdot \vec{\epsilon}_s(\vec{p}) = \epsilon^i_r(\vec{p})\epsilon^i_s(\vec{p})=\delta_{rs}
$$
Above two is satisfied by polarization vectors by definition.
Using these we can show, that
$$
\tag 1
\sum_r \epsilon^i_r(\vec{p})\epsilon^j_r(\vec{p}) = \delta^{ij} - \frac{p^ip^j}{|\vec{p}|^2}
$$
This can be done by using the properties of the polarization vectors in a following way.
When we take a scalar product of the L.H.S. with a polarization vector we expect to get
$$
\epsilon^i_s(\vec{p})\sum_r \epsilon^i_r(\vec{p})\epsilon^j_r(\vec{p}) = \sum_r \delta_{rs}\epsilon^j_r(\vec{p})=\epsilon^j_s(\vec{p})
$$
To achieve this, we just need a delta function on the R.H.S. Hence the first entry $\delta^{ij}$ in our expression $(1)$.
Now let's take a scalar product of L.H.S. with a momentum vector,
$$
p^i\sum_r \epsilon^i_r(\vec{p})\epsilon^j_r(\vec{p}) =0
$$
For the R.H.S. to produce the same result we must add an extra term to $\delta^{ij}$, which will cancel out the momentum vector. However this extra term should also be consistent with the previous case, when we multiplied both sides with a polarization vector.
Little thought will lead you to the expression $(1)$.
Here's a check:
$$
\epsilon^i_s(\vec{p}) \left(\delta^{ij} - \frac{p^ip^j}{|\vec{p}|^2} \right)= \left(\epsilon^j_s(\vec{p})- \frac{\epsilon^i_s(\vec{p})p^ip^j}{|\vec{p}|^2} \right) =\epsilon^j_s(\vec{p})
$$
$$
p^i \left(\delta^{ij} - \frac{p^ip^j}{|\vec{p}|^2} \right)=\left(p^j - \frac{p^ip^ip^j}{|\vec{p}|^2} \right)= 0
$$
