Deriving the unitary operator $U(R)$ associated with a rotation $R$ using Wigner's theorem A rotation $R(\hat{\textbf{n}},\phi)$ about an arbitrary axis $\hat{\textbf{n}}$ through an angle $\phi$ in the three-dimensional physical space is given by $$R(\hat{\textbf{n}},\phi)=e^{-i(\textbf{j}\cdot\hat{n})\phi}\tag{1}$$ where $\textbf{j}=(j_1,j_2,j_3)$ is dimensionless and $j_i=-j_i^T$ ($i=1,2,3$) owing to the condition $R^TR=\mathbb{1}$. 
In quantum mechanics, by Wigner's theorem, the rotation will be represented by an unitary operator $U(R)$. The operator $U(R)$ is usually quoted to have the form (see Sakurai's Modern Quantum Mechanics, for example) $$U(R(\hat{\textbf{n}},\phi))=e^{-i(\textbf{J}\cdot\hat{\textbf{n}})\phi/\hbar}\tag{2}$$ where $\textbf{J}^\dagger=\textbf{J}$ (follows from $U^\dagger U=\mathbb{1}.$) 
Is it possible to derive equation (2) from Eq. (1)? In other words, given $R$, how do we construct the map U(R) that acts on the Hilbert space? 
I used two different symbols $\textbf{j}$ and $\textbf{J}$ in (1) and (2) respectively because unless we actually construct that map $U: R\to U(R)$ it's not clear what $\textbf{j}$ in relation (1) has to do with $\textbf{J}$ in relation (2).

Attempt Since rotations form a Lie group all representations have the exponential form. Therefore, without any loss of generality, the unitary representation in the Hilbert space must have the form $$U(R)=e^{-i\textbf{J}(\textbf{j})\cdot\hat{\textbf{n}}\phi}$$ where $\textbf{J}(\textbf{j})$ is a function of $\textbf{j}$, and unitarity implies $\textbf{J}(\textbf{j})^\dagger=\textbf{J}(\textbf{j}).$ Now it remains to find what $\textbf{J}(\textbf{j})$ is. One option might be to expand $U(R(\hat{\textbf{n}},\delta\phi))=U\Big(\mathbb{1}-i(\textbf{j}\cdot\hat{\textbf{n}})\delta\phi\Big)$ in a Taylor series? But I'm not sure how to carry it out.
 A: Let an infinitesimal vector with direction $\:\mathbf{n}=\left(n_{1},n_{2},n_{3}\right)\:,\Vert\mathbf{n}\Vert=1$ 
\begin{equation}
\delta \boldsymbol{\theta}=\delta \theta\mathbf{n}
\tag{001}
\end{equation}
An infinitesimal rotation $\:R\left(\delta \boldsymbol{\theta}\right)\:$ around $\:\mathbf{n}\:$ through the angle $\:\delta\theta\:$  in the real space $\:\mathbb{R}^3\:$ is
\begin{equation}
\mathbf{r'}\approx \mathbf{r}+\delta \boldsymbol{\theta}\boldsymbol{\times}\mathbf{r}
\tag{002}
\end{equation}
and in matrix form
\begin{equation}
R\left(\delta \boldsymbol{\theta}\right)\approx
\begin{bmatrix}
  \hphantom{-}1&-\delta\theta n_3&\hphantom{-}\delta\theta n_2\\
  \hphantom{-}\delta\theta n_3&\hphantom{-}1&-\delta\theta n_1\\
  -\delta\theta n_2&\hphantom{-}\delta\theta n_1&\hphantom{-}1
 \end{bmatrix}
\tag{003}
\end{equation}
Suppose now that we have a quantum system in a state with wave function $\:\psi\left(\mathbf{r}\right)$. The infinitesimal rotation $\:R\left(\delta \boldsymbol{\theta}\right)\:$ induces a transformation $\:U\,\left(\delta \boldsymbol{\theta}\right)\:$ of old wave function $\:\psi\left(\mathbf{r}\right)$ to a new one $\:\psi'\left(\mathbf{r'}\right)\:$ such that
\begin{equation}
\psi'\left(\mathbf{r'}\right)=U\left(\delta \boldsymbol{\theta}\right)\psi\left(\mathbf{r'}\right)=\psi\left(\mathbf{r}\right)\,,\quad \mathbf{r'}=R\left(\delta \boldsymbol{\theta}\right)\mathbf{r}
\tag{004}
\end{equation}
that is
\begin{align}
U\left(\delta\boldsymbol{\theta}\right)\psi\left(\mathbf{r}\right) & =\psi\left[R^{\boldsymbol{-}1}\left(\delta \boldsymbol{\theta}\right)\mathbf{r}\right]
\approx \psi\left[ \mathbf{r}-\delta \boldsymbol{\theta}\boldsymbol{\times}\mathbf{r}\right]
\nonumber\\
& \approx \psi\left(\mathbf{r}\right)-\left(\delta \boldsymbol{\theta}\boldsymbol{\times}\mathbf{r}\right)\boldsymbol{\cdot}\boldsymbol{\nabla}\psi\left(\mathbf{r}\right)=\psi\left(\mathbf{r}\right)-\dfrac{i}{\hbar}\left(\delta \boldsymbol{\theta}\boldsymbol{\times}\mathbf{r}\right)\boldsymbol{\cdot}\mathbf{p}\psi\left(\mathbf{r}\right)
\nonumber\\
& =\psi\left(\mathbf{r}\right)-\dfrac{i}{\hbar}\delta \boldsymbol{\theta}\boldsymbol{\cdot}\left(\mathbf{r}\boldsymbol{\times}\mathbf{p}\right)\psi\left(\mathbf{r}\right)=\psi\left(\mathbf{r}\right)-\dfrac{i}{\hbar}\delta \boldsymbol{\theta}\boldsymbol{\cdot}\mathbf{L}\psi\left(\mathbf{r}\right)
\tag{005}
\end{align}
so we could write
\begin{equation}
U\left(\delta\boldsymbol{\theta}\right)\approx I-\dfrac{i}{\hbar}\delta \boldsymbol{\theta}\boldsymbol{\cdot}\mathbf{L}
\tag{006}
\end{equation}
where 
\begin{equation}
\mathbf{L}=\mathbf{r}\boldsymbol{\times}\mathbf{p}
\tag{007}
\end{equation}
the  orbital angular momentum.
Suppose now that we want to find $\:U\left(\boldsymbol{\theta}\right)=U\left(\theta\mathbf{n}\right)\:$  for a rotation through a finite angle $\:\theta\:$. If we apply on this an infinitesimal rotation through
$\:\mathrm d\boldsymbol{\theta}=\mathrm d\theta\mathbf{n}\:$ then according to (006)
\begin{equation}
U\left(\boldsymbol{\theta}+\mathrm d\boldsymbol{\theta}\right)=\left(I-\dfrac{i}{\hbar}\mathrm d \boldsymbol{\theta}\boldsymbol{\cdot}\mathbf{L}\right)U\left(\boldsymbol{\theta}\right)
\tag{008}
\end{equation}
or
\begin{equation}
\dfrac{U\bigl[\left(\theta+\mathrm d\theta\right)\mathbf{n}\bigr]-U\bigl[\left(\theta\right)\mathbf{n}\bigr]}{\mathrm d\theta}=\left(-\dfrac{i}{\hbar}\mathbf{n}\boldsymbol{\cdot}\mathbf{L}\right)U\bigl[\left(\theta\right)\mathbf{n}\bigr]
\tag{009}
\end{equation}
so $\:U(\boldsymbol{\theta})=U(\theta\mathbf{n})\:$ satisfies the differential equation 
\begin{equation}
\dfrac{\mathrm d U\left(\theta\mathbf{n}\right)}{\mathrm d\theta}=\left(-\dfrac{i}{\hbar}\mathbf{n}\boldsymbol{\cdot}\mathbf{L}\right)U\left(\theta\mathbf{n}\right)
\tag{010}
\end{equation}
integrated to
\begin{equation}
U\left(\boldsymbol{\theta}\right)=U\left(\theta\mathbf{n}\right)=\exp\left(-\dfrac{i}{\hbar}\theta\mathbf{n}\boldsymbol{\cdot}\mathbf{L}\right)=\exp\left(-\dfrac{i}{\hbar}\boldsymbol{\theta}\boldsymbol{\cdot}\mathbf{L}\right)
\tag{011}
\end{equation}
since $\:U\left(\boldsymbol{0}\right)=U\left(0\cdot\mathbf{n}\right)=I$.
A: The scheme is bassically to check how infinitesimal rotations work. Let me work with rotations around the z-axis, without loss of generality.
Step 1. The rotation around OZ in $\mathbb{R}^3$ is well known:
$$ \left( \begin{array}{ccc}
\cos(\theta) & -\sin(\theta) & 0 \\
\sin(\theta) & \cos(\theta) & 0 \\
0 & 0 & 1 \end{array} \right)$$
If we rotate about an infinitesimal angle, we can say (Taylor to 1st term):
$$ R_z(\varepsilon)\simeq \left( \begin{array}{ccc}
1 & -\varepsilon & 0 \\
\varepsilon & 1 & 0 \\
0 & 0 & 1 \end{array} \right)$$
Which can be expressed as
$$ R_z(\varepsilon)=\mathbb{I}-i\varepsilon G_z; \qquad \text{ with } G_z=\left( \begin{array}{ccc}
0 & -i & 0 \\
i & 0 & 0 \\
0 & 0 & 1 \end{array} \right)  $$
In the same way, you can build $G_x$ and $G_y$, and you see they commute like the angular momentum:
$[G_x, G_y]=iG_z $ + cyclic permutations.
Note that rotation matrices are orthogonal, while $G$'s are Hermitian.

Step 2.
Rotations about the same axis commute (not being true for different axis). Plus, additivity holds this way:
$$R_z (\theta+\alpha)=R_z (\theta)R_z (\alpha)$$
In particular, for an infinitesimal rotation,
$$R_z (\theta+\varepsilon)=R_z (\theta)R_z(\epsilon)=R_z(\theta)\cdot(\mathbb{I}-i\varepsilon G_z )$$
Solving the parenthesis...
$$R_z (\theta+\varepsilon)=R_z(\theta)-i\varepsilon G_z R_z(\theta)$$
$$R_z (\theta+\varepsilon)-R_z(\theta)=i\varepsilon G_z R_z(\theta)$$
$$\dfrac{R_z (\theta+\varepsilon)-R_z(\theta)}{\varepsilon}=-iG_z R_z(\theta)$$
And if you take the limit $\varepsilon\rightarrow 0$,
$$ \frac{d}{d\theta} R_z(\theta)=-iG_z R_z(\theta)$$
Whose solution is the exponential:
$$R_z(\theta)=e^{-i\theta G_z}$$
In particular, you can solve this series using that $G^{even}=\mathbb{I}_2$, and $G_z^{odd}=G_z$
Thus the series can be written as $\cos(\theta) \mathbb{I}_2 + \sin(\theta) G_z$, which is actually $R_z(\theta)$.
That's why $G_z$ is the infinitesimal generator of rotations around z-axis. Its exponential gives the finite rotation.

step 3
So this is in $\mathbb{R}^3$$. What about the Hilbert space?
The key is that a rotation must keep everything invariant if you are rotating everything at the same time. So
$$\varphi'(\vec{x}_f)=\varphi(\vec{x}_0)$$
In other words, "The new wavefunction in the new point, must have the same value as the "old function" in the previous point." That is, "after the rotation, everything has the same value as the old function in the point before rotation".
We can write $\varphi'(\vec{x})=\varphi(R^{-1} \vec{x})$.
And what's that?
$$ R^{-1}(\varepsilon)\cdot\vec{x}=\left( \begin{array}{ccc}
1 & -\varepsilon & 0 \\
\varepsilon & 1 & 0 \\
0 & 0 & 1 \end{array} \right) \left( \begin{array}{c} x \\ y \\ z \end{array} \right) = \left( \begin{array}{c} x+\varepsilon y \\ y-\varepsilon x \\ z \end{array} \right) $$
So $\varphi'(\vec{x})=\varphi(x+\varepsilon y, y-\varepsilon x, z)$.
Since $\varepsilon$ is meant to be very Small, a Taylor development gives:
$$ \simeq \varphi(x,y,z) + \frac{\partial \varphi}{\partial x} \cdot(\varepsilon y ) + \frac{\partial \varphi}{\partial y} \cdot(-\varepsilon x ) $$
And this is
$$ \simeq \varphi(x,y,z) -i \varepsilon \left(y \frac{\partial \varphi}{\partial x} - x\frac{\partial \varphi}{\partial y} \right) = \left (\mathbb{I}-i\varepsilon \frac{L_z}{\hbar} \right )$$
So we get exactly the same form as infinitesimal generators in $\mathbb{R^3}$. The equation has the same form, and hence the same solution.
Conclusion: $L_z$ is the infinitesimal generator of rotations in Hilbert space. When an infinitesimal rotation is performed in the ordinary space, the wavefunction changes in the same way but instead using $L_z / \hbar$.
Replace $L_z$ by $L_x n_x + Ly n_y + L_z n_z$ to get any general rotation about any axis.
A: In theory, you find the matrix representation of $J_x,J_y$ and $J_z$, construct $\hat n\cdot\vec J\phi/\hbar$ and then exponentiate.  Of course this is a mess unless the rotation is particularly simple and the size of the matrices are small.
In practice, it is possible to go from $e^{-i\hat n\cdot \vec j}$ to the factorized form $e^{-i\omega_z j_z}e^{-i\omega_y j_y} e^{-i \Omega_z j_z}$ where $\omega_z$ etc are complicated functions of $\hat n$ and $\phi$.  This relationship, which does not depend on the size of the matrices, can be found in specialized text, such as Varshalovich, D.A., Moskalev, A.N. and Khersonskii, V.K.M., 1988. Quantum theory of angular momentum.  Other sources include textbooks by Rose, by Edmonds etc.  (See also this wiki page for some additional results related to this factorization.)
With this, the job has been transformed to that of finding the various exponentials, but thankfully there are extensive tables of matrix elements of $e^{-i\omega_y J_y/\hbar}$ (they are the Wigner little-$d$ functions $d_{mm’}(\omega_y)$.). It may or may not be worthwhile to invest time in constructing these functions.  (Wolfram Mathematica has these built-in.)
Faced with the alternative of actually exponentiating $e^{-i\hat n\cdot\vec J /\phi/\hbar}$ one can either give it to a computer (again, the Mathematica function MatrixExp will do that) or - if the size of the matrix is small - you do it by hand by writing
$$
-i\hat n\cdot \vec J\phi/\hbar = T\cdot \Lambda \cdot T^{-1}
$$
where $\Lambda$ is diagonal and contains the eigenvalues of $-i\hat n\vec J\phi/\hbar$. These eigenvalues are nothing but $-im\phi$, with $-J\le m\le J$.  You can then exponentiate directly to get
$$
T\cdot e^{\Lambda}\cdot T^{-1}\, .
$$
Finding the transformation matrix $T$, which is constructed from the eigenvectors of $\hat n\cdot \vec J$, is not so simple when the matrices are large (plus, you need to find all the eigenvectors).
A: None of the answers seem to have noticed that equation (1) in the question is wrong. A rotation takes a real 3-d position vector to another real 3-d position vector so that there cannot be an imaginary $i$ in the argument of equation (1). If $J$ is an element of the Lie algebra of the rotation group then a finite rotation about an angle $\theta$ is,
$$
R=e^{\theta J}
$$
This is the corrected form of equation (1) in the question. Let's check that this makes sense. If we make a rotation about the z axis, the Lie algebra element is,
$$
J=\left[\begin{array}{ccc}
0 & 1 & 0 \\
-1 & 0 & 0 \\
0 & 0 & 0
\end{array}\right]
$$
The $J$ matrix only acts on the $x-y$ plane. Restricting the space to the $x-y$ plane, $JJ=-1$ and then it's easy to work out the exponential using a Taylor expansion. The result is (acting only on the $x-y$ plane)
$$
R=1\cos{\theta}+J\sin{\theta}
$$
If we put the z axis back in, the last formula recovers the usual rotation matrix for a rotation about the z-axis.
$$
R=\left[\begin{array}{ccc}
\cos{\theta} & \sin{\theta} & 0 \\
-\sin{\theta} & \cos{\theta} & 0 \\
0 & 0 & 1
\end{array}\right]
$$
Now go over to quantum mechanics, we make the rotation into a unitary operator with Dirac's prescription (units $\hbar=1$),
$$
J\rightarrow -i\hat{J}
$$ 
The unitary operator for a rotation is now,
$$
\hat{U}(R)=e^{-i\theta\hat{J}}
$$
