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.
Best Answer
I haven't thought about this one before, so here is an approach that will work if you work hard enough at it.
Before I begin banging on, point number 1:
Unquestionably the latter. It is a bipartite system and its state space is the tensor product of the two particle spaces. It simply cannot be anything else.
The basic principle here is conservation of angular momentum, so your basic procedure to solve your problem is:
Work out the matrices for the observables for the three nett angular momentum components (the three nett angular momentum operators);
Find the most general Hamiltonian which commutes for all of these three as commutation with the Hamiltonian is equivalent to invariance with time of all the moments of probability distributions of the measurements.
Part 1: The Three Angular Momentum Operators
The $x$-AM component observable for the spin half particle,
$$\sigma_x=\left(\begin{array}{cc}0&1\\1&0\end{array}\right)$$
has AM eigenvectors:
$$\psi_+=\frac{1}{\sqrt{2}}\left(\begin{array}{c}1\\1\end{array}\right);\quad\psi_-=\frac{1}{\sqrt{2}}\left(\begin{array}{c}1\\-1\end{array}\right)$$
and AM eigenvalues $\lambda_+=+\frac{1}{2}$ and $\lambda_-=-\frac{1}{2}$, respectively.
The $x$-AM component observable for the spin 1 particle,
$$S_x = \left(\begin{array}{ccc}0&0&0\\0&0&i\\0&-i&0\end{array}\right)$$
has AM eigenvectors:
$$\Psi_+=\frac{1}{\sqrt{2}}\left(\begin{array}{c}0\\i\\1\end{array}\right);\quad\Psi_-=\frac{1}{\sqrt{2}}\left(\begin{array}{c}0\\1\\i\end{array}\right);\quad\Psi_0=\left(\begin{array}{c}1\\0\\0\end{array}\right)$$
and AM eigenvalues $\Lambda_+=+1$, $\Lambda_-=-1$ and $\Lambda_0=0$, respectively. So now, for the two particle system, the six $x$-AM eigenstates are:
and so, if we order the eigenstates as above, the eigenvectors as columns are $\mathrm{vec}(\psi_+\otimes\Psi_+),\,\mathrm{vec}(\psi_+\otimes\Psi_0)\cdots$ (see the Wikipedia Vectorization Page) and so at last we get as the total $x$-AM component observable $\Sigma_X = P_X \Lambda_X P_X^\dagger$ where
$$P_X=\left( \begin{array}{cccccc} 0 & \frac{1}{\sqrt{2}} & 0 & 0 & \frac{1}{\sqrt{2}} & 0 \\ 0 & \frac{1}{\sqrt{2}} & 0 & 0 & -\frac{1}{\sqrt{2}} & 0 \\ \frac{i}{2} & 0 & \frac{1}{2} & \frac{i}{2} & 0 & \frac{1}{2} \\ \frac{i}{2} & 0 & \frac{1}{2} & -\frac{i}{2} & 0 & -\frac{1}{2} \\ \frac{1}{2} & 0 & \frac{i}{2} & \frac{1}{2} & 0 & \frac{i}{2} \\ \frac{1}{2} & 0 & \frac{i}{2} & -\frac{1}{2} & 0 & -\frac{i}{2} \\ \end{array} \right)$$
and $\Lambda_X =\mathrm{diag}\left(\frac{3}{2},\,\frac{1}{2},\,\frac{-1}{2},\,\frac{1}{2},\,\frac{-1}{2},\,\frac{3}{2}\right)$. The result is:
$$\Sigma_X=\left( \begin{array}{cccccc} 0 & \frac{1}{2} & 0 & 0 & 0 & 0 \\ \frac{1}{2} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{3}{4} & -\frac{1}{4} & \frac{i}{4} & \frac{3 i}{4} \\ 0 & 0 & -\frac{1}{4} & \frac{3}{4} & \frac{3 i}{4} & \frac{i}{4} \\ 0 & 0 & -\frac{i}{4} & -\frac{3 i}{4} & \frac{3}{4} & -\frac{1}{4} \\ 0 & 0 & -\frac{3 i}{4} & -\frac{i}{4} & -\frac{1}{4} & \frac{3}{4} \\ \end{array} \right)$$
From here on it should be conceptually clear how to go, although tedious. You do the same for the $y$-AM observables:
$$\sigma_y=\left(\begin{array}{cc}0&-i\\i&0\end{array}\right)$$ $$S_y = \left(\begin{array}{ccc}0&0&-i\\0&0&0\\i&0&0\end{array}\right)$$
to find the total system $y$-AM observable $\Sigma_Y$ and for the $z$-AM observables:
$$\sigma_z=\left(\begin{array}{cc}i&0\\0&-i\end{array}\right)$$ $$S_z = \left(\begin{array}{ccc}0&i&0\\-i&0&0\\0&0&0\end{array}\right)$$
to get the total system $z$-AM observable $\Sigma_Z$.
Part 2: Find most general Hamiltonian
Your most general Hamiltonian will be defined by the three commutator relationships expressing conservation of AM:
$$[\hat{H},\,\Sigma_j]=0;\;j=X,\,Y,\,Z$$
You'll need to work out the invariant spaces of the three $\Sigma$s to do this. You'll get a linear space of possible $\hat{H}$s: in the two coupled spin half particles case there is essentially only one possible Hamiltonian that falls out of this approach and that is one proportional to $\sigma_x\otimes\sigma_x+\sigma_y\otimes\sigma_y+\sigma_z\otimes\sigma_z$ (plus a term proportional to the $4\times4$ identity matrix expressing the shift in the ground state energy) but this six dimensional case things will be a bit more complicated. Now as I said, I've never done this before, so I daresay there is a more systematic and less cumbersome way to work this all out. But any method is going to rest on the first principles expressed above.
Magnetic Field
What are the terms for the influence of the magnetic field. Well that's an easy one: in the ordering we have studied above, the uncoupled Hamiltonian will be:
$$\hat{H} = \gamma_{\frac{1}{2}}\left(\sigma_x\,B_x + \sigma_y\,B_y+ \sigma_z\,B_z\right)\otimes 1_{3\times3} + \gamma_1\,I_{2\times2}\otimes\left(S_x\,B_x + S_y\,B_y+ S_z\,B_z\right)$$
where $\gamma_{\frac{1}{2}}$ and $\gamma_1$ are the respective gyromagnetic ratios.
Notes on completing the method. You can also represent a bipartite state $\Phi=\psi\otimes\Psi$ as the literal $2\times 3$ matrix that is the outer product $\Phi=\psi Psi^T$ of the $2\times 1$ and $3\times 1$ column vectors. Then the operator on the first space act on the left and the operators on the second act on the right. So our $x$-component observable would be the linear, homogeneous transformation:
$$\Phi\mapsto \sigma_x\,\Phi\,S_x^T$$
and the vectorization operator (See Vectorization Wiki Page), which reorders our states into a $6\times 1$ column vectors as in my answer, writes this as
$$\mathrm{vec}(\Phi) \mapsto S_x\otimes\sigma_x\,\mathrm{vec}(\Phi)$$
Using the standard formula $\mathrm{vec}(A\,B\,C) = C^T\otimes A \,\mathrm{vec}(B)$. By dint of the formula $(A\otimes B)\, (C\otimes D) = (A\,C)\otimes(B\,D)$, and using the fact that inverse, complex conjugate, Hermitian conjugate and transpose operations distribute over the Kronecker produt, we can diagonalize $S_x\otimes\sigma_x$ inside the Kronecker product and find that the coupled system’s eigenstates are $\Pi_x\otimes \pi_x$, where $P_X,\,p_x$ are the matrices of eigenvectors of the individual multiplicands written as columns. So this will let you calculate $\Sigma_j,\,j=X,\,Y,\,Z$ systematically and fast.
Now to find the most general Hamiltonian, you need to find the invariant space of the group of matrices generated by the three matrices $\exp(i\,\Sigma_j)$ and find the irreducible representation of it: equivalently the smallest vector subspace of $\mathbb{C}^6$ left invariant by the group: by Schur’s lemma, any matrix commuting with all three must be proportional to the identity operator when restricted to this subspace. The scaling factor is possibly nought – i.e. the operator could possibly be the zero endomorphism. This completely characterizes the most general Hamiltionian: it can be any operator which is proportional to the identity when restricted to this irreducible subspace.
You could also find the subspace which is common to all three of the nullspaces of the three $36\times 36$ matrices $1_{36\times 36}\otimes \Sigma_j - \Sigma_j^T\otimes 1_{36\times 36}$ in Mathematica or Matlab, but I suspect there is a much eleganter method grounded on Schur's lemma!