You wouldn't think it, from how easy it is to pose this question, but it is ridiculously nontrivial. As it happens, it is entirely impossible to find the position-basis matrix elements of this propagator.
So far you've done good, and the identification
$$
U\left( t_{2},t_{1}\right)
=e^{\frac {-i} {\hbar }H\left( t_{2}-t_{1}\right) }
=\sum_{n=-\infty}^\infty e^{-i\left( t_{2}-t_{1}\right) E_{n}/\hbar}\lvert n\rangle\langle n\rvert
$$
is correct. Naively, the path is clear, because you know the eigenfunctions,
$$
⟨\theta|n⟩=\frac{1}{\sqrt{2\pi}}e^{in\theta},
$$
so you can in principle just sandwich $U(t_2,t_1)$ between two position kets and you'd be done. This gives you
\begin{align}
\langle \theta_2 | U\left( t_{2},t_{1}\right) | \theta_1\rangle
& =
\sum_{n=-\infty}^\infty e^{-i\left( t_{2}-t_{1}\right) E_{n}/\hbar}\langle \theta_2 | n\rangle\langle n|\theta_1\rangle
\\ & =
\frac{1}{2\pi} \sum_{n=-\infty}^\infty e^{-i\left( t_{2}-t_{1}\right) E_{n}/\hbar}e^{in(\theta_1-\theta_2)},
\end{align}
and everything looks fine, right? You've just got one series to calculate and you'll be on your way. Moreover, you can simplify things here a fair bit: you can set $\theta=\theta_1-\theta_2$ and $\tau=t_2-t_1$ to simplify notation (i.e. the propagator only depends on those differences, which is nice), and you can also calculate the energies,
$$
E_n=\frac{\hbar^2}{2mR^2}\left(n+\frac{e}{h}\phi\right)^2,
$$
so that
$
\frac{(t_2-t_1)E_n}{\hbar}
= \frac{\hbar}{2mR^2}\left(n+\frac{e}{h}\phi\right)^2\tau
=\kappa n^2\tau + \lambda n\phi\tau + \mu\phi^2\tau,
$
and all of this makes for a nicely simplified series:
\begin{align}
\langle \theta_2 | U\left( t_{2},t_{1}\right) | \theta_1\rangle
& =
\frac{e^{-i\mu\phi^2\tau}}{2\pi} \sum_{n=-\infty}^\infty e^{-in^2\kappa \tau }e^{in(\theta-\lambda\phi\tau)},
\end{align}
and this is clean and simple as can be. What could possibly go wrong?
Now, as it happens, this series is exactly summable, and it represents something called a Jacobi $\vartheta_3$ function, and you can read all about it in the DLMF . More precisely, in this case, the relationship comes down to
\begin{align}
\langle \theta_2 | U\left( t_{2},t_{1}\right) | \theta_1\rangle
& =
\frac{1}{\sqrt{4\pi i \kappa\tau}}
e^{-i\mu\phi^2\tau}
e^{\frac{i (\lambda\phi\tau-\theta)^2}{4 \kappa\tau}}
\vartheta_3\mathopen{}\left(
\frac{\pi }{2 \kappa\tau}(\lambda\phi\tau-\theta),
e^{i\frac{ \pi^2}{\kappa\tau}}
\right)\mathclose{}
,
\tag{$\ast$}
\end{align}
and it kinda looks like you're done: you've reduced your series to some special function, and you obviously cannot whittle that down any further, so there's nothing else to do.
Unfortunately, however, you're not actually done yet. If you look more closely at the theta function, it is normally defined as
$$
\vartheta_3(z,q)=1+2\sum_{n=0}^\infty q^{n^2}\cos(2nz),
$$
where $q$ and $z$ are complex parameters, and we will want to set $q=e^{i \pi^2/\kappa\tau}$. It is obvious from the definition that the series will converge if $|q|<1$, and that it will diverge if $|q|>1$. The quantum mechanics result we just got is balanced on the knife edge, at $|q|=1$: there are no convergence guarantees, but the series is oscillatory, so maybe we can hope for some conditional convergence brought about by the oscillations?
Alas, there's no joy there either. The full answer is that $\vartheta_3$ has something called a natural boundary, which means that for any point $q_0$ on the unit circle (so $|q_0|=1$) and any $z\in\mathbb C$, the limit $\lim_{q\to q_0^-} \vartheta_3(z,q)$ diverges to infinity. (I'm using the notation $q\to q_0^-$ to imply that $q$ stays inside the disk during the limit, just to keep things easy.) If you want the big ugly mathematician details, this is because the series is something called a lacunary series, which I asked about here on MathOverflow and here on Maths SE; in particular, here the Fabry gap theorem implies that the theta function has a ring of fire around it, which sort of looks like this:
Mathematica code through Import["https://raw.githubusercontent.com/halirutan/Mathematica-SE-Tools/master/SETools/SEImageExpressionDecode.m"]["http://i.stack.imgur.com/wM9dN.png"]
The problem here is that our $\vartheta_3$ is sitting right on top of the natural boundary: it is an analytical function for $|q|<1$, but there is provably no way to provide an analytical continuation past its circle of convergence (and that is kind of a bummer), and there is provably no way to make any sense of it for $|q|=1$, which is exactly where we are.
So... where does that leave us? Our nice propagator $(\ast)$ looked neat enough, but a closer inspection shows that it is actually mathematical nonsense. What does this mean, and how do we fix it?
Well, there's a couple of ways, but they both essentially mean that you are only allowed to interpret $\langle \theta_2 | U\left( t_{2},t_{1}\right) | \theta_2+\theta \rangle$ as something called a generalized function, also known as a distribution: an object, kind of like the Dirac delta function, that doesn't make sense on its own, but only when integrated in products with other functions.
There's two easy ways to do this:
One is to take the complex-$q$ thing seriously. This is easily achieved by giving the time $\tau=\tau_0-i\delta$ a nonzero imaginary component, which brings $|q|$ down to under $1$, and lets $⟨ \theta_2 | U\left( t_{2},t_{1}\right) | \theta_1 ⟩$ make sense.
If you do this, the key thing to note is that normally the only physically meaningful quantities where you will use the propagator are matrix elements of the form $⟨f|U(t_2, t_1)|g⟩$ for physical states $|f⟩$ and $|g⟩$. Those ultimately boil down to integrals of the form
\begin{align}
⟨f|U(t_2, t_1)|g⟩
& =
\int_0^{2\pi}
⟨0|U(t_2, t_1)|\theta⟩
F(\theta)\mathrm d\theta
\\&=
\int_0^{2\pi}
\frac{
e^{-i\mu\phi^2\tau}
e^{\frac{i (\lambda\phi\tau-\theta)^2}{4 \kappa\tau}}
}{\sqrt{4\pi i \kappa\tau}}
\vartheta_3\mathopen{}\left(
\frac{\pi }{2 \kappa\tau}(\lambda\phi\tau-\theta),
e^{i\frac{ \pi^2}{\kappa\tau}}
\right)\mathclose{}
F(\theta)\mathrm d\theta
,
\end{align}
and here the recipe says that you should read those integrals as the limit of integrals with a nonzero $\delta$ as that $\delta$ goes to zero:
\begin{align}
⟨f|U(t_2, t_1)|g⟩
& =
\lim_{\delta\to 0}
\int_0^{2\pi}
\frac{
e^{-i\mu\phi^2\tau}
e^{\frac{i (\lambda\phi\tau-\theta)^2}{4 \kappa\tau}}
}{\sqrt{4\pi i \kappa\tau}}
\vartheta_3\mathopen{}\left(
\frac{\pi }{2 \kappa\tau}(\lambda\phi\tau-\theta),
e^{i\pi^2/\kappa(\tau-i\delta)}
\right)\mathclose{}
F(\theta)\mathrm d\theta
.
\end{align}
If you do that, what happens is that for each $\delta$ the integrand is well defined, analytical and perfectly well behaved, and it gives you a perfectly normal result for the integral. More importantly, those results will approach a limit as $\delta\to0$ as long as $F(\theta)$ is not absolutely crazy.
The other approach is similar, but it's series-based, i.e. it requires you to take our original series before the evaluation to $\vartheta_3$,
\begin{align}
\langle \theta_2 | U\left( t_{2},t_{1}\right) | \theta_1\rangle
& =
\frac{e^{-i\mu\phi^2\tau}}{2\pi} \sum_{n=-\infty}^\infty e^{-in^2\kappa \tau }e^{in(\theta-\lambda\phi\tau)},
\end{align}
and see it as indicating that you need to take the $|n|\to\infty$ limit after you've done any relevant inner products you need to take. For the integral above, for example, this means the understanding that
\begin{align}
\int_0^{2\pi}
⟨0|U(t_2, t_1)|\theta⟩
F(\theta)\mathrm d\theta
&=
\frac{e^{-i\mu\phi^2\tau}}{2\pi}
\sum_{n=-\infty}^\infty
e^{-in^2\kappa \tau }
\int_0^{2\pi}
e^{in(\theta-\lambda\phi\tau)}
F(\theta)\mathrm d\theta
,
\end{align}
with the series after the integral. If $F(\theta)$ is well-behaved, then its Fourier series coefficients will go down with $|n|$ reasonably fast, and the series will (whew!) converge again - as it needs to do for a physically relevant matrix element.
Alternatively, you can take the series expression we had earlier,
\begin{align}
\langle \theta_2 | U\left( t_{1}+\tau,t_{1}\right) | \theta_2+\theta\rangle
& =
\frac{e^{-i\mu\phi^2\tau}}{2\pi} \sum_{n=-\infty}^\infty e^{-in^2\kappa \tau }e^{in(\theta-\lambda\phi\tau)},
\end{align}
and simply see it as an identity between distributions.
This paradigm requires you to see each mapping $\theta\mapsto e^{-in^2\kappa \tau }e^{in(\theta-\lambda\phi\tau)}$ not as a function but as a distribution, which you are then adding up together as a series. This is a pretty strong thing to do, since it requires you to stop seeing $⟨\theta_2|n⟩⟨n|\theta_1⟩$ as something which has a value if you give it a specific $\theta_1$, $\theta_2$ and $n$, but if you do this you get some nice mathematical advantages.
More specifically, each of the distributions $\theta\stackrel{\mathrm{DS}}{\mapsto} e^{-in^2\kappa \tau }e^{in(\theta-\lambda\phi\tau)}$ is a tempered distribution (i.e. the kind of distribution that plays well with Fourier transforms and quantum mechanics; see these notes, these ones or these ones for more details), and the series is bounded above termwise by the Dirac comb, which is also a tempered distribution, and this is (probably?) plenty to show that the series sums to a tempered distribution.
Thus, in this paradigm, you get that the propagator exists and it is as nice as you'd like it to be (within the bounds that it is a distribution) and that it is indeed the limit of the series on the right hand side, so long as you take it in the distributional sense.
That's sort of it, really, and I think this is as clear as I can make it. I'll close with the two references where I learned this stuff,
L. Schulman. A path integral for spin. Phys. Rev. 176 no. 5, p. 1558. (1968).
and
L. S. Schulman. Techniques and applications of path integration (Dover, New York, 2005), pp. 190-196.
and finally with a quote from Schulman (actually present in both references), which describes the situation here rather accurately:
The degree of pathology exhibited by this Green's function is entertaining, especially in view of the elementary nature of the example.
Hint :
In my opinion you must deal with the inverse problem : that of finding the propagator $\:K(\mathbf{x},t\;;\;\mathbf{x}_{0},t_{0})\:$ from the non-relativistic Schr$\ddot{\rm{o}}$dinger equation, so to understand why it is the transition probability amplitude for a particle present at position $\:\mathbf{x}_{0}\:$ on time $\:t_{0}\:$ to propagate at position $\:\mathbf{x}\:$ on time $\:t\:$ and how to use it for the solutions of the equation.
If by your efforts you could answer your question directly without dealing with the inverse problem then you would gain much on mathematics but lose much on physics.
In the following I sketch the main points of the inverse problem :
1. Let the non-relativistic Schr$\ddot{\rm{o}}$dinger equation
\begin{equation}
i\hbar\;\dfrac{\partial \psi(\mathbf{x},t)}{\partial t}=H(\mathbf{x},t)\;\psi(\mathbf{x},t)
\tag{sk-01}
\end{equation}
where $\:H(\mathbf{x},t)\:$ the Hamiltonian, an hermitian operator.
2. Let $\:H\:$ be dependent only on the position coordinate $\:\mathbf{x}\:$ but not on time, that is $\:H=H(\mathbf{x})$.
3. Now we suppose that $\:H\:$ has (real) eigenvalues $\:E_{\jmath}\:$ with corresponding eigenfunctions $\:u_{\jmath}(\mathbf{x})\:$
\begin{equation}
H(\mathbf{x})u_{\jmath}(\mathbf{x})=E_{\jmath}\;u_{\jmath}(\mathbf{x}) \qquad \jmath=1,2,\cdots
\tag{sk-02}
\end{equation}
where the index $\:\jmath\:$ takes discrete or continuous values. For simplicity we suppose that $\:\jmath\:$ takes discrete values and that the set of eigenfunctions $\:\lbrace u_{\jmath} \rbrace, \jmath=1,2,...\:$ is a complete basis of the Hilbert space of states $\:\psi $. So, any state $\:\psi(\mathbf{x},t)\:$ can be expressed as
\begin{equation}
\psi(\mathbf{x},t)=\sum_{\jmath}a_{\jmath}(t) \; u_{\jmath}(\mathbf{x})
\tag{sk-03}
\end{equation}
4. So
\begin{equation}
\sum_{\jmath}\left[ i\hbar\;\dfrac{da_{\jmath}(t)}{dt} -E_{\jmath}a_{\jmath}(t)\right] u_{\jmath}(\mathbf{x})=\mathbf{0}
\tag{sk-04}
\end{equation}
and from the completeness of the basis $\:\lbrace u_{\jmath}\rbrace\:$
\begin{equation}
a_{\jmath}(t)=a_{\jmath}(t_{0})\; \exp \left[ -i \frac{E_{\jmath}}{\hbar}(t-t_{0})\right] \qquad \jmath=1,2,\cdots
\tag{sk-05}
\end{equation}
5. The coefficients $\:a_{\jmath}(t_{0}) \:$ are determined from the initial state of the system
\begin{equation}
\psi_{0}(\mathbf{x}) \equiv \psi(\mathbf{x},t_{0})=\sum_{\jmath}a_{\jmath}(t_{0}) \; u_{\jmath}(\mathbf{x})
\tag{sk-06}
\end{equation}
6. Under these conditions we have the following general solution of the Schr$\ddot{\rm{o}}$dinger equation (sk-01)
\begin{equation}
\psi(\mathbf{x},t)=\sum_{\jmath}a_{\jmath}(t_{0})\; \exp \left[ -i \frac{E_{\jmath}}{\hbar}(t-t_{0})\right] \; u_{\jmath}(\mathbf{x})
\tag{sk-07}
\end{equation}
7. Now, equation (sk-01) is linear, that is, if to two initial states $\:\psi_{0}^{(1)}(\mathbf{x}) \:$ and $\:\psi_{0}^{(2)}(\mathbf{x}) \:$ at time $\: t_{0}\:$ there correspond the solutions $\:\psi^{(1)}(\mathbf{x},t) \:$ and $\:\psi^{(2)}(\mathbf{x},t) \:$ respectively, then to the initial state $\:\psi_{0}=a_{1}\psi_{0}^{(1)}+a_{2}\psi_{0}^{(2)}\:$, where $\:a_{1},a_{2}\:$ are complex numbers, there corresponds the solution $\:\psi=a_{1}\psi^{(1)}+a_{2}\psi^{(2)}$.
Any initial state $\:\psi_{0}(\mathbf{x})\:$ may be considered as the superposition (linear combination) of impulses expressed via the Dirac function $\:\delta^{3}(\mathbf{x})\:$ as
\begin{equation}
\psi_{0}(\mathbf{x}) = \int \limits_{\mathbf{x}_{0}} \delta^{3}(\mathbf{x}_{0}-\mathbf{x}) \;\psi_{0}(\mathbf{x}_{0}) \; d^{3}\mathbf{x}_{0}
\tag{sk-08}
\end{equation}
8. If our system is a single point particle and this particle is appeared suddenly at point $\:\mathbf{x}_{0}\:$ on time $\: t_{0}\:$, then the expression $\:\delta^{3}(\mathbf{x}_{0}-\mathbf{x})\:$ under the integral represents the state (wave) function of this particle. If to this initial state there corresponds the solution
\begin{align}
K(\mathbf{x},t \; \boldsymbol{;} \;\mathbf{x}_{0},t_{0})\equiv \: & \text{the state function at point $\:\mathbf{x}\:$ on time $\:t\:$}
\nonumber\\
&\text{corresponding to the appearance of the}
\nonumber\\
&\text{particle at point $\:\mathbf{x}_{0}\:$ on time $\:t_{0}\:$}
\tag{sk-09}
\end{align}
then, because of the linearity discussed in the previous paragraph, the solution $\psi(\mathbf{x},t)$ corresponding to the initial state $\psi_{0}(\mathbf{x})$ of equation (sk-08) would be
\begin{equation}
\psi(\mathbf{x},t) = \int \limits_{\mathbf{x}_{0}}K(\mathbf{x},t \; \boldsymbol{;} \;\mathbf{x}_{0},t_{0}) \; \psi_{0}(\mathbf{x}_{0}) \; d^{3}\mathbf{x}_{0}
\tag{sk-10}
\end{equation}
9. In mathematical terms $\:K(\mathbf{x},t \; ; \;\mathbf{x}_{0},t_{0})\:$ is called kernel while in physical terms is called propagator since it is the transition probability amplitude for a particle present at position $\:\mathbf{x}_{0}\:$ on time $\:t_{0}\:$ to propagate at position $\:\mathbf{x}\:$ on time $\:t$.
10. We'll determine $\:K(\mathbf{x},t \; ; \;\mathbf{x}_{0},t_{0})\:$ from equation (sk-07)
\begin{equation}
K(\mathbf{x},t \; \boldsymbol{;} \;\mathbf{x}_{0},t_{0})=\sum_{\jmath}c_{\jmath}\; \exp \left[ -i \frac{E_{\jmath}}{\hbar}(t-t_{0})\right] u_{\jmath}(\mathbf{x})
\tag{sk-11}
\end{equation}
where $\:c_{\jmath}\:$ are the "coordinates" of the initial state $\:\delta^{3}(\mathbf{x}_{0}-\mathbf{x})\:$ relatively to the complete basis $\:\lbrace u_{\jmath} \rbrace, \jmath=1,2,...\:$ according to equation (sk-06)
\begin{equation}
\delta^{3}(\mathbf{x}_{0}-\mathbf{x}) =\sum_{\jmath}c_{\jmath}u_{\jmath}(\mathbf{x})
\tag{sk-12}
\end{equation}
11. Taking the inner product of above equation with $\:u_{k}(\mathbf{x})\:$ we have
\begin{equation}
\left\langle \delta^{3}(\mathbf{x}_{0}-\mathbf{x}),u_{k}(\mathbf{x})\right \rangle =\sum_{\jmath}c_{\jmath}\left\langle u_{\jmath}(\mathbf{x}),u_{k}(\mathbf{x})\right \rangle
\tag{sk-13}
\end{equation}
The complete basis $\:\lbrace u_{\jmath} \rbrace, \jmath=1,2,...\:$ is orthogonal since its members are eigenfunctions of the
hermitian operator $\:H$. Without loss of generality we suppose that these eigenfunctions are normalized (unit norm) so that
\begin{equation}
\left\langle u_{\jmath}(\mathbf{x}),u_{k}(\mathbf{x})\right \rangle = \int \limits_{\mathbf{x}} u_{\jmath}(\mathbf{x}) \;u_{k}^{\boldsymbol{*}}(\mathbf{x})\;d^{3}\mathbf{x}=\delta_{\jmath k}
\tag{sk-14}
\end{equation}
By equations (sk-14), (sk-13)
\begin{equation}
c_{k}=\left\langle \delta^{3}(\mathbf{x}_{0}-\mathbf{x}),u_{k}(\mathbf{x})\right \rangle = \int \limits_{\mathbf{x}} \delta^{3}(\mathbf{x}_{0}-\mathbf{x}) \;u_{k}^{\boldsymbol{*}}(\mathbf{x})\;d^{3}\mathbf{x}=u_{k}^{\boldsymbol{*}}(\mathbf{x}_{0})
\tag{sk-15}
\end{equation}
that is
\begin{equation}
c_{\jmath}=u_{\jmath}^{\boldsymbol{*}}(\mathbf{x}_{0})
\tag{sk-16}
\end{equation}
so equations (sk-12) and (sk-11) yield respectively
\begin{equation}
\delta^{3}(\mathbf{x}_{0}-\mathbf{x}) =\sum_{\jmath}u_{\jmath}^{\boldsymbol{*}}(\mathbf{x}_{0})\;u_{\jmath}(\mathbf{x})
\tag{sk-17}
\end{equation}
and
\begin{equation}
K(\mathbf{x},t \; \boldsymbol{;} \;\mathbf{x}_{0},t_{0})=\sum_{\jmath} u_{\jmath}^{\boldsymbol{*}}(\mathbf{x}_{0})\;\exp \left[ -i \frac{E_{\jmath}}{\hbar}(t-t_{0})\right] u_{\jmath}(\mathbf{x})
\tag{sk-18}
\end{equation}
Best Answer
Yes, you are correct. Actually a few pages later Sakurai derives the same relation: