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.