How do I apply the spin operator to a wavefunction? I think this might be a dumb question so I apologize in advance. I have a model for a particle subject to a radial potential. Specifically the spatial part of the wavefunction takes the form
$$\psi_{\ell,m}(r,\theta,\phi)=R_{\ell}(r)Y_{\ell,m}(\theta,\phi)$$
Here $Y$ denote the Spherical Harmonics, and $R$ is some arbitrary function describing the radial dependence. I'd like to be able to compute this particle's total angular momentum. (In particular the end goal is to compute the eigenstates of the operators $\widehat{\mathbf J_z}$ and $\widehat{J^2}$). From what I've seen on Wikipedia, this is done by using the total angular momentum operator,
$$\hat{\mathbf J}=\hat{\mathbf L}+\hat{\mathbf S}$$
Here the L is orbital angular momentum, and it has the form
$$\hat{\mathbf L}=-i\hbar (\hat{\mathbf r}\times \nabla)$$
Where r is the position operator. This makes sense - I can apply this operator to a scalar valued wavefunction $\psi$ and get a vector out on the other side. However, I can't see any similar forms for the spin. I've seen the spin operator written in other places as
$$\hat{\mathbf S}=\frac{\hbar}{2}(\boldsymbol{\sigma}_x,\boldsymbol{\sigma}_y,\boldsymbol{\sigma}_z)$$
Where the sigmas are the Pauli spin matrices. To be frank I'm not really sure what the tuple of matrices is supposed to mean, and I don't know how to apply it to my wavefunction $\psi$, especially since it is scalar valued.
Can someone help me figure this out? I am an applied mathematician, not a physicist, so I'm having a hard time comprehending some of the literature on quantum mechanics.
 A: The wavefunction you wrote down is, as you say, the spatial part of the vector.  The full Hilbert space is $L^2(\mathbb R^3)\otimes \mathbb C^2$, an element of which has the form
$$|\psi\rangle = \sum_{i=1}^2\int \mathrm d^3 x \ \psi_i(\mathbf x) \bigg(|\mathbf x\rangle  \otimes |\alpha_i\rangle\bigg)$$
where $\{|\alpha_i\rangle\}$ is a basis for $\mathbb C^2$.
If the Hamiltonian doesn't feature any spin-splitting, then it can be expressed in the form $H = H_\mathrm{space} \otimes \mathbb I$, which means that its (at least twice-degenerate) eigenvectors can be expressed as eigenvectors of $H_\mathrm{space}$ tensored with arbitrary elements of $\mathbb C^2$.
In contrast, the spin operators to which you refer are of the form $S_i = \mathbb I \otimes \bigg(\frac{\hbar}{2} \sigma_i\bigg)$, with $\sigma_i$ the Pauli matrices, so they leave the spatial part of the vectors alone while acting on the spin part. The orbital angular momentum operators are of the form $L_i \otimes \mathbb I$, as they act exclusively on the spatial part of the vectors while leaving the spin alone.  Finally, the total angular momentum operators are $J_i = L_i \otimes \mathbb I + \mathbb I \otimes S_i$, which is colloquially written as $J_i = L_i + S_i$ with the tensor product structure swept under the rug.

As an explicit example, the vector corresponding to a spin-up electron with spatial wavefunction $\psi_{n\ell m}(\mathbf x)$ is
$$ |\psi \rangle = \int \mathbf d^3 x \ \psi_{n\ell m}(\mathbf x ) |\mathbf x\rangle \otimes |\alpha_\uparrow\rangle $$
whose wavefunction can be written
$$\langle \mathbf x, \hat z|\psi \rangle = \psi_{n\ell m}(\mathbf x) \pmatrix{1 \\ 0} = \pmatrix{\psi_{n\ell m}(\mathbf x) \\ 0 }$$
The action of the spin operator $\mathbb I \otimes S_x$ on this vector is
$$\mathbb I \otimes S_x |\psi \rangle = \int \mathrm d^3 x \ \psi_{n\ell m}(x) \mathbb I |\mathbf x\rangle \otimes S_x  |\alpha_\uparrow\rangle$$
whose wavefunction is
$$\langle \mathbf x,\hat z|\mathbb I \otimes S_x|\psi\rangle = \psi_{n\ell m}(\mathbf x) \pmatrix{0\\ \frac{\hbar}{2}}$$
A: The expression
$$\psi_{\ell,m}(r,\theta,\phi)=R_{\ell}(r)Y_{\ell,m}(\theta,\phi)$$
is incomplete because it contains only the spatial part, but neglects the spin part.
Actually the wave function $\psi$ is not just a single position-dependent
function, but it is a "vector" of two such functions, together called
a spinor wave function.
$$\psi(r,\theta,\phi)
 =\begin{pmatrix}\psi_+(r,\theta,\phi)\\\psi_-(r,\theta,\phi)\end{pmatrix}$$
Now, the orbital angular momentum operator $\hat{\mathbf L}$
(and also  other spatial operators like $\hat{\mathbf r}$ or $\nabla$)
act on the spatial dependency only.
And they act equally on the 2 spinor components.
For example $L_z=-i\hbar\frac{\partial}{\partial\phi}$,
the $z$-component of the orbital angular momentum,
acts on the spinor wave function like this:
$$\begin{align}
L_z\psi(r,\theta,\phi)
 &=-i\hbar\frac{\partial}{\partial\phi}
   \begin{pmatrix}
   \psi_+(r,\theta,\phi) \\
   \psi_-(r,\theta,\phi)
   \end{pmatrix} \\
 &=-i\hbar\begin{pmatrix}
   \frac{\partial\psi_+(r,\theta,\phi)}{\partial\phi} \\
   \frac{\partial\psi_-(r,\theta,\phi)}{\partial\phi}
   \end{pmatrix}
\end{align}$$
On the other hand, the spin angular momentum operator $\mathbf S$
(and also its components $S_x$, $S_y$, $S_z$ and the Pauli matrices
$\sigma_x$, $\sigma_y$, $\sigma_z$) act differently on the two spinor components.
And they act independent of the position ($r,$, $\theta$, $\phi$).
So these operators are constant $2\times 2$ matrices, i.e. when applied to
a $2$-component spinor they produce another $2$-component spinor.
For example $S_z=\frac{\hbar}{2}\sigma_z$, the $z$-component
of the spin angular momentum operator, acts like this:
$$\begin{align}
S_z \psi(r,\theta,\phi)
 &=\frac{\hbar}{2}\begin{pmatrix}1&0\\0&-1\end{pmatrix}
   \begin{pmatrix}\psi_+(r,\theta,\phi)\\\psi_-(r,\theta,\phi)\end{pmatrix} \\
 &=\frac{\hbar}{2}
   \begin{pmatrix}\psi_+(r,\theta,\phi)\\-\psi_-(r,\theta,\phi)\end{pmatrix}
\end{align}$$
From the two examples above you can derive how $J_z$, the $z$-component
of the total angular momentum, acts on the spinor function.
$$\begin{align}
J_z\psi(r,\theta,\phi)
 &=(L_z+S_z)\psi(r,\theta,\phi) \\
 &=\left(-i\hbar\frac{\partial}{\partial\phi}
  +\frac{\hbar}{2}\begin{pmatrix}1&0\\0&-1\end{pmatrix}\right)
  \begin{pmatrix}\psi_+(r,\theta,\phi)\\ \psi_-(r,\theta,\phi)\end{pmatrix} \\
 &=\hbar\begin{pmatrix}
   -i\frac{\partial\psi_+(r,\theta,\phi)}{\partial\phi}+\frac{1}{2}\psi_+(r,\theta,\phi) \\
   -i\frac{\partial\psi_-(r,\theta,\phi)}{\partial\phi}-\frac{1}{2}\psi_-(r,\theta,\phi)
  \end{pmatrix}
\end{align}$$
