Perturbation of coupled spin I am given a system with Hamiltonian (all 1/2 spins)
$$H_0=\alpha(S_1\cdot S_2)$$
I broke it down and found that there were four eigenstates: $|1,[0,\pm1]\rangle$ and $|0,0\rangle$. Each has an associated eigenvalues.
I am then told that a weak magnetic field is applied giving a perturbation
$$H_V=-\mu B_0\cdot(\sigma_1+\sigma_2)$$
I am supposed to use stationary state perturbation theory to find $E_n^{(1)}$.  So I know that the perturbed energies are given as
$$E_n^{(1)}=\langle \psi_n^{(0)}|H_V|\psi_n^{(0)}\rangle$$
I am a bit confused on what to use for the $\psi_n$ and then the $\sigma$ matrices. I know that the eigenstates are those given above, but I do not see how I can form a matrix from those and then I am unsure of what sigma would be in order to calculate the matrix element
 A: It sounds like you're trying to find the shifts in the energy levels caused by $H_V$. To find the energy shift that goes as the first order in your small parameter you compute
$$E_n^{(1)} = \langle \psi_n^{(0)}|H_V|\psi_n^{(0)} \rangle \quad (1)$$
as you noted. In this expression $|\psi_n^{(0)}\rangle$ means "the $n^{\textrm{th}}$ eigenstate of the unperturbed Hamiltonian." That superscript (0) means "to zeroth order in the perturbation" ie. without any perturbation.
So to compute the energy level shifts you have to know the eigenstates of the unperturbed Hamiltonian. You say you already found those to be $|1,\pm 1\rangle$, $|1,0\rangle$ and $|0,0\rangle$. To make contact with the notation in equation (1) let's go ahead and number these
$$|\psi_1^{(0)}\rangle = |0,0 \rangle \quad |\psi_2^{(0)}\rangle = |1,0 \rangle \quad |\psi_3^{(0)}\rangle = |1,-1 \rangle \quad |\psi_4^{(0)}\rangle = |1,+1 \rangle$$
Now you want to compute the energy shifts. Let's just do one of them as an example. Let's compute $E_1^{(1)}$.
$$\begin{eqnarray}
E_1^{(1)} &=& \langle \psi_1^{(0)}|H_V|\psi_1^{(0)}\rangle \\
&=& -\mu B \cdot \langle0,0|\sigma_1 + \sigma_2 | 0,0 \rangle \\
&=& -\mu B \cdot \left( \langle0,0|\sigma_1|0,0\rangle + \langle0,0|\sigma_2|0,0\rangle \right) \end{eqnarray}$$
Now, those two operators $\sigma_1$ and $\sigma_2$ are the spin operators for each individual particle. The state $|0,0\rangle$ is expressed in a basis that is not the single particle basis. $|0,0\rangle$ is a state designed to be simple in the basis of total spin. In order to work out the matrix elements you need for the computation you have to either convert the operators to the total spin basis, or re-express the states in the individual spin basis. I think the latter is easier in this case. Just use
$$|0,0\rangle \equiv \left( \uparrow \downarrow - \downarrow \uparrow \right)/\sqrt{2}.$$
Now it's easy to compute what you need. Here's an example
$$B \cdot \langle 0,0|\sigma_1| 0,0\rangle = \frac{1}{2} \langle\uparrow \downarrow-\downarrow \uparrow|B^x\sigma^x_1 + B^y\sigma^y_1 + B^z\sigma^z_1|\uparrow\downarrow - \downarrow\uparrow\rangle.$$
You asked about what the matrices for $\sigma$ should be. Now that we've got everything expressed in the single particle spin basis, it's easy. For example,
$$\sigma^x_1 = \left(\begin{array}{cc} 0 & 1 \\ 1 & 0 \end{array} \right).$$
The thing to note is that this only acts on the first particle, so you get things like this
$$\begin{eqnarray}
\langle \uparrow \downarrow | \sigma^x_1 | \downarrow \uparrow \rangle &=& \langle \uparrow | \sigma^x_1 | \downarrow \rangle\cdot\langle \downarrow|\uparrow \rangle \\
&=& \langle \uparrow|\uparrow \rangle \cdot \langle \downarrow|\uparrow \rangle \\
&=& 1 \cdot 0 \\
&=& 0.
\end{eqnarray}$$
At this point I'm sure you can figure out what's going on. I think your problem was one of notation and hopefully seeing this written out helps.
Extension:
I've shown how to compute simple terms like
$$\langle \uparrow \downarrow | \sigma^x_1 | \downarrow \uparrow \rangle$$
but the original expression to compute was more complex. In particular we had
$$\langle\uparrow \downarrow-\downarrow \uparrow|B^x\sigma^x_1 + B^y\sigma^y_1 + B^z\sigma^z_1|\uparrow\downarrow - \downarrow\uparrow\rangle.$$
The trick here is that this is all linear so we can break it into little parts like this
$$\langle\uparrow \downarrow-\downarrow \uparrow|B^x\sigma^x_1 + B^y\sigma^y_1 + B^z\sigma^z_1|\uparrow\downarrow - \downarrow\uparrow\rangle$$
$$=B_x\langle\uparrow \downarrow-\downarrow \uparrow|\sigma^x_1|\uparrow\downarrow - \downarrow\uparrow\rangle + B_y\langle\uparrow \downarrow-\downarrow \uparrow|\sigma^y_1|\uparrow\downarrow - \downarrow\uparrow\rangle + B_z\langle\uparrow \downarrow-\downarrow \uparrow|\sigma^z_1|\uparrow\downarrow - \downarrow\uparrow\rangle.$$
Each one of these terms can be broken down further:
$$\langle\uparrow \downarrow-\downarrow \uparrow|\sigma^x_1|\uparrow\downarrow - \downarrow\uparrow\rangle$$
$$\langle \uparrow \downarrow|\sigma^x_1|\uparrow\downarrow\rangle - \langle \uparrow\downarrow | \sigma^x_1 | \downarrow \uparrow\rangle - \langle \downarrow \uparrow | \sigma^x_1 | \uparrow \downarrow \rangle + \langle \downarrow \uparrow | \sigma_1^x | \downarrow \uparrow \rangle$$
From here you can surely see how to to calculate the desired quantity. Note that if you think carefully you can tell that many terms are zero without having to calculate them.
