My original answer was going to be much longer, but everything started to be a bit too long and complex, so here's some rough elements to give you how to derive the idea :
Formally, the Hadamard function (which is related to the Green's function by some differing contour integral) is defined on curved spacetimes by
$$G^{(1)}(x, y) = \sum_{\gamma} \frac{\Delta_\gamma^{\frac{1}{2}}(x,y)}{4\pi^2} \left[ \frac{1}{\sigma_{\gamma}(x,y)} + v_\gamma(x,y) \ln(|\sigma_\gamma(x,y)|) + \varpi_\gamma(x,y) \right]$$
Where we sum over every possible geodesic connecting two points. $\Delta$ is the van Vleck determinant, $\sigma$ the geodetic distance between $x$ and $y$, and $v$ and $\varpi$ are non-singular functions defined by the theory.
The Feynmann propagator is defined similarly, although the functions within that sum differ. So the solution is to take the sum of those geodesics.
We also have that the relation between the Feynman propagator and the relativistic point particle is
\begin{eqnarray}
G_F(x, y) &=& \langle T[ \phi(x), \phi(y) ] \rangle\\
&=& \int \mathcal{D}\phi\ \phi(x) \phi(y) e^{i S[\phi]}\\
&=& \int_0^\infty \frac{dT}{T} \int_{x(0) = x}^{x(T) = y} \mathcal{D}x \exp\left[\int_0^T d\tau (g(\dot{x}, \dot{x}) + m^2) \right]
\end{eqnarray}
There's a bunch of additional things to do for this to work properly, ie adding the einbein, ghost fields and so on, but this is roughly the relation we have between the two, and both evaluate to the appropriate propagator. The trick to do this, and I suspect, for the non-unique geodesic case, is to use the decomposition of the action into the extremal solution and variations, ie
\begin{eqnarray}
x(\tau) = x_C(\tau) + y(\tau)
\end{eqnarray}
So that $x_C(\tau)$ is the classical solution that extremize the action, and $y$ is any function on $[0,T]$ with $y(0) = y(T) = 0$. In the case of a cut locus, there are multiple such solutions.
Since we are dealing with massless fields here, for the flat space case we just have
\begin{eqnarray}
G_F(x,y) = \frac{i}{8\pi^2 \sigma(x,y)} - \frac{1}{8\pi} \delta(\sigma(x, y))
\end{eqnarray}
or, for the Hadamard function,
\begin{eqnarray}
G^{(1)}(x,y) = -\frac{1}{4\pi^2 \sigma(x,y)}
\end{eqnarray}
Now let's consider the case of the cylinder. The Minkowski cylinder is the manifold $\mathbb{R} \times S$, with metric
\begin{equation}
ds^2 = -dt^2 + d\theta^2
\end{equation}
On the cylinder, the propagator evaluates to
\begin{equation}
G_C(x,y) = \sum_{k \in \mathbb{Z}} G(x, y + (0, kL))
\end{equation}
This changes the geodetic interval by
\begin{eqnarray}
\sigma(x, y + (0, kL)) &=& |x - y - (0, kL)|^2\\
&=& -(x_t - y_t)^2 + (x_\theta - y_\theta - kL)^2
\end{eqnarray}
This is very much equivalent to the geodetic interval of the geodesic connecting $x$ to $y$ with winding number $k$, which corresponds indeed to every geodesics connecting the two. On the cylinder, two points are connected by the helix
\begin{eqnarray}
\gamma(\lambda) &=& (a \lambda + t_0, b \lambda + \theta_0)
\end{eqnarray}
at $\gamma(0) = x$, $t_0 = x_t$ and $\theta_0 = x_\theta$, and at $\gamma(1) = y$, $a + x_t = y_t$ and $b + x_\theta = y_\theta + k L$, so that the geodesic with winding number $k$ has the geodesic interval
\begin{eqnarray}
\sigma_{\gamma_k}(x, y) &=& \int_0^1 -(x_t - y_t)^2 + (x_\theta - y_\theta - kL)^2 d\lambda\\
&=& -(x_t - y_t)^2 + (x_\theta - y_\theta - kL)^2
\end{eqnarray}
I attempted to do the case of the sphere, but I'm afraid this is getting a bit long for a short answer. This should be roughly equivalent, considering that on a sphere, the geodesics are parametrized by 1) the winding number 2) if the points are antipodal, the azimuthal angle. You can check (in Birrell & Davies, for instance) that in the case of the Einstein static universe (roughly a static sphere spacetime) that similarly to the cylinder, the propagator depends on an infinite sum over the winding number. I assume that in addition, there is an integral over $\varphi$ in the antipodal case that is already worked out here.
