There are many ways to derive the Rarita-Schwinger equation. Below, I will describe the most common two for you.
1) In the representation theory of Lorentz group, Dirac spinor \psi^{n} transforms in the representation space \left(0 , 1/2 \right) \oplus \left(1/2 , 0\right), and a Lorentz vector V_{\mu} transforms in \left(1/2 , 1/2\right). One way of obtaining the spin-3/2 representation is by taking the tensor product V_{\mu} \otimes \psi^{n} \equiv \Psi^{n}_{\mu} \in \left(1/2 , 1/2\right) \otimes \left[ \left(0 , 1/2 \right) \oplus \left(1/2 , 0\right) \right] \ , and reducing it by subtracting all invariant subspaces: Notice that \Psi^{n}_{\mu} is a 16-component object, i.e., twice the number needed for the description of spin-3/2 particle and its anti-particle. Indeed, from the rule of angular momentum addition (in the rest frame) we get \mbox{Spin} (\Psi_{\mu}) : \left( 1 + 0 \right) \otimes \frac{1}{2} = \frac{3}{2} + \frac{1}{2} + \frac{1}{2} \ . Thus, beside the pure spin-3/2 component, the tensor product field \Psi^{n}_{\mu} has two (extra) spin-1/2 components. So, we need to get rid of these extra components. We do this by the following decomposition (from now on I will suppress the Dirac indices) \Psi_{\mu} = \left( \Psi_{\mu} - \frac{1}{4} \gamma_{\mu} \gamma^{\sigma} \Psi_{\sigma} \right) + \frac{1}{4} \gamma_{\mu} \gamma^{\sigma} \Psi_{\sigma} \ . In the last term, the (4-component) factor \gamma^{\rho} \Psi_{\rho} \in [ (1/2 , 0) \oplus (0 , 1/2)] represents one of those extra spin-1/2 components. So, the pure spin-3/2 component is contained in the first term: \psi_{\mu} \equiv \Psi_{\mu} - \frac{1}{4} \gamma_{\mu} \gamma^{\rho} \Psi_{\rho} \in [ (1 , 1/2) \oplus (1/2 , 1)] \ . Now, the above-defined spinor-vector \psi_{\mu} satisfies the following 4 constraints \gamma^{\mu}\psi_{\mu} = \gamma^{\mu} \Psi_{\mu} - \frac{4}{4} \gamma^{\rho} \Psi_{\rho} = 0 \ . So, there are 16 - 4 = 12 components left in \psi_{\mu}, 8 of them describe the irreducible spin-3/2 representation, and the other 4 are carried by a spin-1/2 component still hidden in the spinor-vector \psi_{\mu}. This is, however, not a problem because for each value of the Lorentz index (\mu = 0 , \cdots , 3) the field \psi_{\mu} satisfies the following Dirac equation \left( i \gamma^{\sigma} \partial_{\sigma} - m \right) \psi_{\mu} = 0 \ . \ \ \ \ \ (1) From this it follows that the hidden spin-1/2 component is carried by \partial^{\mu}\psi_{\mu} since it is a Dirac spinor: \left( i \gamma^{\sigma} \partial_{\sigma} - m \right) \partial^{\mu}\psi_{\mu} = 0 \ . However, eq(1) together with the constraint \gamma^{\mu}\psi_{\mu} = 0 imply that \partial^{\mu}\psi_{\mu} = 0. To see this, contract eq(1) from the left with \gamma^{\mu} and then use the relation \gamma^{\mu}\gamma^{\sigma} = 2 \eta^{\mu\sigma} - \gamma^{\sigma} \gamma^{\mu} and the constraint \gamma^{\mu}\psi_{\mu} = 0. Thus, eq(1) together with \gamma^{\mu}\psi_{\mu} = 0 describe the irreducible (massive) spin-3/2 representation of the Lorentz group. This is indeed equivalent to the fancy looking Rarita-Schwinger equation \left( \epsilon^{\mu\nu\rho\tau} \gamma_{5} \gamma_{\nu} \partial_{\rho} - i m \sigma^{\mu\tau}\right) \psi_{\tau} = 0 \ , where \sigma^{\mu\nu} = \frac{i}{2}[\gamma^{\mu} , \gamma^{\nu}]. To see this, we start with i \gamma^{\mu}\partial_{\mu} \psi_{\nu} - m \psi_{\nu} = 0 \ , \ \ \ (1)\partial^{\mu}\psi_{\mu} = 0 \ . \ \ \ \ \ \ \ \ \ \ (2) We will, for a moment, forget about the constraint \gamma^{\mu}\psi_{\mu} = 0. Contracting eq(1) from the left with \gamma^{\nu} and using \gamma^{\nu} \gamma^{\mu} = 2 \eta^{\nu\mu} - \gamma^{\mu} \gamma^{\nu} and eq(2), we obtain i \gamma^{\mu}\gamma^{\nu} \partial_{\mu} \psi_{\nu} + m \gamma^{\nu}\psi_{\nu} = 0 \ . \ \ \ \ (3) Side remark: Eq(3) shows that the field \gamma^{\nu}\psi_{\nu} is not a Dirac spinor because the mass term has a plus sign. Multiply eq(3) from the left with \gamma^{\rho} to obtain i \gamma^{\rho}\gamma^{\mu}\gamma^{\nu} \partial_{\mu} \psi_{\nu} + m \gamma^{\rho}\gamma^{\nu}\psi_{\nu} = 0 \ . \ \ (4) Now, in the first term we use the identity \gamma^{\rho}\gamma^{\mu}\gamma^{\nu} = \eta^{\rho \mu} \gamma^{\nu} - \eta^{\rho \nu} \gamma^{\mu} + \eta^{\mu \nu} \gamma^{\rho} + i \epsilon^{\rho \mu \nu \tau} \gamma_{5} \gamma_{\tau} \ , and, in the second term of eq(4), we use the identity \gamma^{\rho} \gamma^{\nu} = \eta^{\rho \nu} - \frac{i}{2} \epsilon^{\rho \nu \mu \tau} \gamma_{\mu} \gamma_{\tau} \gamma_{5} \ . This leads, after using eq(1) and eq(2), to \epsilon^{\rho \tau \mu \nu} \gamma_{5} \gamma_{\tau} \partial_{\mu} \psi_{\nu} + \frac{im}{2} \epsilon^{\rho \mu \nu \tau} \gamma_{\nu} \gamma_{\tau} \gamma_{5}\psi_{\mu} - i \partial^{\rho}(\gamma^{\nu} \psi_{\nu}) = 0 . Finally, we use the following identity in the second term \sigma^{\rho \mu} = - \frac{1}{2} \epsilon^{\rho \mu \nu \tau} \gamma_{\nu} \gamma_{\tau} \gamma_{5} \ . This gives you \epsilon^{\rho \tau \mu \nu} \gamma_{5} \gamma_{\tau} \partial_{\mu}\psi_{\nu} - i m \sigma^{\rho \mu} \psi_{\mu} - i \partial^{\rho} (\gamma^{\nu} \psi_{\nu}) = 0 \ . Thus, the Rarita-Schwinger equation is obtained when the constraint \gamma^{\nu}\psi_{\nu} = 0 is satisfied. In general, a massive spin (k + 1/2) field with integer k can be described by a spinor-tensor \psi_{\mu_{1} \mu_{2} \cdots \mu_{k}} (x), totally symmetric in the space-time indices, satisfying the generalized Rarita-Schwinger equations \begin{align*}\left( i \gamma \cdot \partial - m\right) \psi_{\mu_{1} \cdots \mu_{k}} (x) & = 0 , \\ \gamma^{\mu_{1}} \psi_{\mu_{1} \cdots \mu_{k}} (x) & = 0 , \\ \partial^{\mu_{1}} \psi_{\mu_{1} \cdots \mu_{k}} (x) & = 0 , \\ \eta^{\mu_{1} \mu_{2}} \psi_{\mu_{1} \mu_{2} \cdots \mu_{k}} (x) & = 0 . \end{align*}
2) The Rarita-Schwinger equation can also be obtained as a special case of the Bargmann-Wigner equations. Recall that Dirac spinor \psi_{n} , \ n = 1 ,2 , \cdots \ , 4 can be regarded as a vector (i.e., one index object or rank-1 tensor) on the complex vector space \mathbb{C}^{4}. It describes relativistic particles of mass m and spin s = 1/2. Notice that the rank of the tensor \psi_{n} \in \mathbb{C}^{4} is 2s. The Bargmann-Wigner equations are a generalization of the Dirac equation in the following sense. A massive field of arbitrary spin s is represented by a totally symmetric tensor \Psi_{n_{1}n_{2} \cdots n_{2s}} \in \mathbb{C}^{4} of rank 2s satisfying Dirac equations in all spin indices n_{i} = 1 , \cdots , 4:
\left( i \gamma \cdot \partial - m \right)_{n_{1}m} \Psi_{m n_{2} \cdots n_{2s}} = 0 \ ,\left( i \gamma \cdot \partial - m \right)_{n_{2}m} \Psi_{n_{1} m \cdots n_{2s}} = 0 \ ,\ \ \ \ \vdots \ \ \ \ \left( i \gamma \cdot \partial - m \right)_{n_{2s}m} \Psi_{n_{1} n_{2} \cdots m} = 0 \ .
The tensors \Psi_{n_{1}n_{2} \cdots n_{2s}} \in \mathbb{C}^{4} are usually referred to as multi-spinors to avoid confusion with spinor-tensors \psi_{n}^{\mu_{1}\mu_{2} \cdots} which are space-time tensors with values in the spin space.
Okay, let us now apply this method to a field with finite mass m and spin s = 3/2. In this case the Bargmann-Wigner multi-spinor is of rank 2s = 3, i.e., it has three indices \Psi_{mnp}(x), and it is completely symmetric. It is assumed that \Psi (x) satisfies the following set of Bargmann-Wigner equations \begin{align*}D_{mq} \Psi_{qnp}(x) & = 0 \ , \\ D_{nq} \Psi_{mqp}(x) & = 0 \ , \\ D_{pq} \Psi_{mnq}(x) & = 0 \ , \end{align*} where D_{mn} are the matrix elements of the Dirac differential operator D_{mn} = \left( i \gamma \cdot \partial - m \right)_{mn} = i \gamma^{\mu}_{mn} \partial_{\mu} - m \delta_{mn} . To make connection with physical objects, we try to express the abstract multi-spinor \Psi_{mnp} field in terms of something more “physical”, i.e., (Clifford algebra)-valued space-time tensors, such as spinor-vector \psi^{\mu}_{n}(x) and spinor-tensor \chi^{\mu\nu}_{n}(x). In 4-dimentional space-time, the generators of the Clifford algebra form a complete set of sixteen 4\times 4 matrices \Gamma^{a} = \{ I , \ i \gamma_{5}, \ \gamma^{\mu}\gamma_{5} , \ \gamma^{\mu} , \ \sigma^{\mu\nu} \} . But \Psi_{mnp} is totally symmetric, so if we are to write the expansion \Psi_{mnp}(x) = (M_{\mu})_{mn}\psi^{\mu}_{p} + (N_{\mu\nu})_{mn} \chi^{\mu\nu}_{p} \ , the set of 4 \times 4 matrices ( M_{\mu} , N_{\mu\nu}) must form a complete set of symmetric Clifford numbers. The Clifford algebra, \Gamma^{a} given above, does not have such a closed subset of symmetric matrices. Now comes the beauty of the charge conjugation matrix C = i\gamma^{2}\gamma^{0} with its many properties. As you can check for yourself, C divides the Clifford algebra into two closed subsets \Gamma^{a}C = \big \{ \ \{ C , \ i \gamma_{5} \ C , \ \gamma^{\mu}\gamma_{5} \ C \} , \ \{ \gamma^{\mu} \ C , \ \sigma^{\mu\nu} \ C \} \ \big \} \ , consisting of 6 antisymmetric matrices \{ C , \ i \gamma_{5} \ C , \ \gamma^{\mu}\gamma_{5} \ C \} \ , and 10 symmetric matrices \{ \gamma^{\mu} \ C , \ \sigma^{\mu\nu} \ C \} . And so, we have the following expansion \Psi_{mnp}(x) = \left( \gamma_{\mu} C \right)_{mn} \psi^{\mu}_{p}(x) + \frac{1}{2} \left( \sigma_{\mu\nu} C \right)_{mn} \chi^{\mu\nu}_{p}(x) \ . \ \ \ (2.1) As you might have guessed by now, this is a lengthy business with a lot of algebra, so I will only describes the main ideas and leave you to fill in the algebraic details.
So, the first thing you need to do is to use the symmetry of the multi-spinor \Psi_{mnp} and see what conditions you get on the fields \psi^{\mu} and \chi^{\mu\nu}: Use the antisymmetric subset of Clifford algebra to obtain the following equations \begin{align*} \Psi_{mnp} C^{-1}_{np} & = 0 \ , \ \ \ \ \ \ (2.2a) \\ \Psi_{mnp} \left( C^{-1} i \gamma_{5}\right)_{np} & = 0 \ , \ \ \ \ \ \ (2.2b) \\ \Psi_{mnp}\left( C^{-1} \gamma_{5} \gamma^{\rho}\right)_{np} & = 0 \ . \ \ \ \ \ \ (2.2c) \end{align*} If you use the properties of C and the \gamma’s, the equations (2.2) should give you the following (I am suppressing the spin indices from now on)
\begin{align*}\gamma_{\mu}\psi^{\mu}(x) + \frac{1}{2} \sigma_{\mu\nu} \chi^{\mu\nu}(x) & = 0 \ , \ \ \ (2.3a) \\ \gamma_{\mu}\gamma_{5} \psi^{\mu}(x) + \frac{1}{2} \sigma_{\mu\nu} \gamma_{5} \chi^{\mu\nu}(x) & = 0 \ , \ \ \ (2.3b) \\ \gamma_{\mu}\gamma_{5} \gamma_{\rho} \psi^{\mu}(x) + \frac{1}{2} \sigma_{\mu\nu} \gamma_{5} \gamma_{\rho} \chi^{\mu\nu}(x) & = 0 \ . \ \ \ (2.3c) \end{align*}
Now, you should check to see if these conditions are independent of each other. Since \{ \gamma_{5} , \gamma^{\mu} \} = 0 and [ \sigma_{\mu\nu} , \gamma_{5}] = 0, the first two equations are equivalent to \gamma_{\mu} \psi^{\mu}(x) = 0 \ , \ \ \ \ \ \ \ (2.4a)\sigma_{\mu\nu} \chi^{\mu\nu}(x) = 0 \ . \ \ \ \ \ \ (2.4b) In eq(2.3c), commute \gamma_{\rho} to the left, multiply by \gamma_{5} from the right and then use [\gamma_{\rho} , \sigma_{\mu\nu}] = 2i ( \eta_{\rho \mu}\gamma_{\nu} - \eta_{\rho \nu} \gamma_{\mu} ) \ , to transform eq(2.3c) into \psi^{\rho}(x) + i \gamma_{\mu} \chi^{\rho \mu}(x) = 0 \ . \ \ \ \ (2.4c) Now you can see that eq(2.4b) can be obtained from eq(2.4a) and eq(2.4c): multiply (2.4c) from the left by \gamma_{\rho}, use (2.4a) and \chi^{\rho \mu} = - \chi^{\mu\rho}. Thus, we see that eq(2.4b) is a redundant condition. So, you are left with the two independent conditions: \gamma_{\mu} \psi^{\mu}(x) = 0 \ , \ \ \ \ \ \ \ (2.5a) \ \ \mbox{4 equations}\psi^{\rho}(x) = i \gamma_{\mu}\chi^{\mu\rho}(x) \ . \ \ \ (2.5b) \ \mbox{16 equations} These are 20 linear equations between the 16 + 24 = 40 components of the fields \psi^{\mu}_{n} and \chi^{\mu\nu}_{n}. So, there are 40 - 20 = 20 independent components in the fields (\psi^{\mu} , \chi^{\mu\nu}) and this is exactly the number of independent components of a totally symmetric rank 3 tensor in a 4-dimensional vector space. In other words, the spinor-vector \psi^{\mu}_{n} and the spinor-tensor \chi^{\mu\nu}_{n} together with the conditions (2.5) contain the same information as the Bargmann-Wigner multi-spinor \Psi_{mnp}. Finally, you need the Bargmann-Wigner equation D_{mq} \Psi_{qnp} (x) = 0. Contracting this with \left(C^{-1} \sigma_{\rho \tau} \right)_{pm}: (C^{-1}\sigma_{\rho \tau})_{pm}(i \gamma_{\sigma}\partial^{\sigma} - m)_{mq} \left( (\gamma_{\mu} C)_{qn} \psi^{\mu}_{p} + \frac{1}{2} (\sigma_{\mu\nu}C)_{qn}\chi^{\mu\nu}_{p}\right) = 0 . So, have fun simplifying this equation. Notice that we are taking the trace, so you need to use the fact that the trace of odd number of gamma matrices is zero. You will also need the trace identities for \mbox{Tr}(\sigma_{\rho\tau} \sigma_{\mu\nu}) and \mbox{Tr}(\sigma_{\rho \tau} \gamma_{\mu} \gamma_{\nu}). And if you don’t make mistakes, you end up with m \chi^{\mu\nu} = \partial^{\mu}\psi^{\nu} - \partial^{\nu}\psi^{\mu} \ . \ \ \ (2.6) Now contract this with i\gamma_{\mu} and use the conditions (2.5) to obtain the Rarita-Schwinger equation (i \gamma \cdot \partial - m )\psi^{\nu}(x) = 0 we considered in the beginning of this post. So, you know how to proceed from this point to derive the on-shell condition \partial_{\mu} \psi^{\mu} = 0 and the rest of the strory.
