What's the Propagator in the Free Particle Case? (Path Integrals with Source Term) If I take the Lagrangian to be,
$$L(t)=\frac{1}{2}m \dot q(t)^2$$
The Euclidean Path Integral is supposed to be,
$$K=\int D[q(t)] \ e^{-\int L(\dot q) d \tau}$$
If I add a source term $J(\tau)$ we obtain,
$$K[J]=\int D[q(t)] \ e^{-\int L(\dot q) -J(\tau) \cdot q(\tau) d\tau}$$
All the books say that the integral is equal to,
$$K[J]=e^{\int \int \ J(\tau) G(\tau,\tau') J(\tau') \ d\tau d\tau'}$$
Where $G(\tau,\tau')$ is the propagator, Green's Function, for the Lagrangian. However, I can't figure out a sensible propagator here. If I plug the Lagrangian into the Euler-Lagrange formula, and add the inhomogeneous term, I get,
$$\frac{m}{2} \partial_t^2 G(\tau-\tau')=\delta(\tau-\tau')$$
Which seems reasonable. However, I wish to calculate $\langle \int q(t)^2 \rangle$, so I need to be able to take the second variational derivative, with respect to $J$. This yields, according to some version of Wick's Theorem,
$$\int d\tau \left( \cfrac{\delta}{\delta J} \right)^2 K[J]=\left \langle \int q(t)^2 \right \rangle=\int 2 \cdot G(\tau-\tau) \ d \tau$$
However, $G(0)$ could be anything, according to the differential for the Green's Function. How do I generally pick boundary conditions for the Green's Function so I can get proper results? I want to be able to try this later with $L(t)=\frac{1}{2}m \dot q(t)^2+\lambda \cdot q(t)$ and after that with a quadratic. So hopefully, the answer can help guide me towards getting correct answers.
For the bounty: I want to see the derivation for the integral with a source term. I also want to see an example of taking the second functional derivative of this derived functional integral.
 A: Jean Zinn-Justin has a great way of teaching path integral techniques starting with finite dimensional random variables (sometimes called "0-dimensional fields"). Here, you should think of these as discrete lattice approximations to continuous fields. In the spirit of Zinn-Justin's approach, I'll describe how this is done for the 1D free particle system you describe above.
Assuming that $q(t)$ vanishes in the infinite past and infinite future, your Lagrangian $L[q]=\int dt\frac{1}{2}m\dot q^2(t)$ can be rewritten as $\int dt\big[-\frac{1}{2}m q \ddot q\big]$, using integration by parts. In this notation, you can think of this as a bilinear $\int dt\big[\frac{1}{2}mq(t)(-\partial_t^2)q(t)\big]=\frac{1}{2}m\langle q,-\partial_t^2q\rangle=\frac{1}{2}q^T A q$. The last expression uses notation from finite dimensional linear algebra in order to temporarily demystify some of the steps. Here, $A$ is thought of as some sort of matrix realization of the linear operator $-m\partial_t^2$. We can also assume that $A=A^T$.
The way to derive the explicit expression for the generating functional $Z[J]$ is by completing the square in the path integral, performing a linear change of variables. To do this, we need to use some notion of the inverse of $A$:
\begin{align*}
\int \mathcal D q \exp\bigg[-\frac{1}{2}q^TAq - J\cdot q\bigg]=&\int\mathcal D q\exp\bigg[-\frac{1}{2}(q+A^{-1}J)^TA(q+A^{-1}J)+\frac{1}{2}J^TA^{-1}J\bigg]\\
=&\int\mathcal Dq\exp\bigg[-\frac{1}{2}q^T A q\bigg]\exp\bigg[\frac{1}{2}J^TA^{-1}J\bigg]\\
=& Z_0\exp\bigg[\frac{1}{2}J^TA^{-1}J^T\bigg].
\end{align*}
To complete the derivation, we need to solve (and pick a convention) for $A^{-1}$, which is called the Green function. This is tricky in general, because of the existence of harmonic functions $f$ that satisfy $Af=0$ (i.e. $A=-m\partial_t^2$ has a non-trivial kernel, consisting of constant functions and a drift term). Physically, you can interpret these harmonic functions more generally as radiation from sources located in the distant past (or less physically, in the far future), and they lead to the distinction between advanced and retarded Green functions. In general, it is conventional to choose the harmonic part of the Green function so that boundary conditions are easily satisfied with linearly independent combinations of the Green function. For example, the retarded Green function is determined by the causality condition that it should vanish in the distant past.
Naturally, the space of allowed $J$'s is restricted to the domain on which $A^{-1}$ makes sense. In practice, the $J$'s satisfy conservation laws, and must decay sufficiently quickly at the boundaries.  
At a purely formal level, computing functional derivatives with respect to $J$ of $Z[J]$ is exactly analogous to computing partial derivatives with respect to $\vec J$ of some discrete finite dimensional approximation to the partition function: in economical finite dimensional notation, $\frac{1}{Z_0}\frac{\delta^2}{\delta J_1 \delta J_2}Z[ J]\Big|_{J=0}:= \frac{\delta^2}{\delta J_1 \delta J_2}\exp\Big[\frac{1}{2}J_i A_{ij}^{-1}J_j\Big]\Big|_{J=0}=\frac{\delta}{\delta J_1}\frac{1}{2}(J_iA_{i2}^{-1}+A_{2i}^{-1}J_i)\exp(\frac{1}{2}J_iA_{ij}^{-1}J_j)\Big|_{J=0}=\frac{1}{2}(A_{12}^{-1}+A_{21}^{-1})\exp(\frac{1}{2}J_iA_{ij}J_j)\Big|_{J=0}+\mathcal O (J^2)\Big|_{J=0}=A_{12}^{-1}=:G(\tau_1,\tau_2)$. 
The Green function can be viewed as an elementary solution to the inhomogeneous version of a linear differential equation. The 2-point correlation function in your example happens to be a Green function for the 1D diffusion equation, but there is a certain amount of freedom in deciding on the Green function depending on the number of independent harmonic modes. Even the correlation function depends on how the space of fields that is integrated over is defined. One common way to choose a Green function is to impose vanishing boundary conditions. For Brownian motion, the choice of boundary conditions can be interpreted as a choice of reference frame.
Now for an explicit example. Consider a fluctuating elastic 'string' with $q(0)=q(L)=0$. To find the correlation function, we first choose an appropriate basis for describing fluctuations: here, we can expand configurations $q(s)$ in Fourier modes $S_n(s)=\sqrt{\frac{2}{L}}\sin(\frac{n\pi s}{L})$, where $n\in\mathbb N$. In this basis, the operator $-\partial_s^2$ acts as $-\partial_s^2 S_n(s)=\left(\frac{n\pi}{L}\right)^2S_n(s)$. Hence, it can be inverted by defining $(-\partial_s^2)^{-1} S_n(s)=\frac{L^2}{\pi^2n^2}S_n(s)$, extending by linearity. To find the integral kernel representation of $(-\partial_s^2)^{-1}$ in position space, we need to evaluate
\begin{align*}
\sum_{n\geq 1} \frac{L^2}{\pi^2n^2}S_n(s')S_n(s)=-\frac{L}{\pi^2}\sum_{n\geq 1} \frac{1}{n^2}\bigg[\cos(\frac{\pi n}{L}(s+s'))-\cos(\frac{\pi n}{L}(s-s'))\bigg]
\end{align*}
In the limit of large $L$, the sum over $n$ can be replaced by an integral over $\omega\equiv\frac{\pi n}{L}$:
\begin{align*}
G(s,s')= -\int_{\frac{\pi}{L}}^\infty\frac{d\omega}{\pi}\frac{1}{\omega^2}\bigg[\cos(\omega(s+s'))-\cos(\omega(s-s'))\bigg]
\end{align*}
This integral can be performed in general by taking appropriate limits of sums of contour integrals (where the sum of integrands, restricted to the real line, must approach the original function $\frac{1-\cos(x)}{x^2}$). Here, a valid sequence of approximations is 
\begin{align*}
G_\epsilon(s,s')=-\int_{-\infty}^\infty \frac{d\omega}{2\pi}\frac{1}{\omega^2+\epsilon}\bigg[\cos(\omega(s+s'))-\cos(\omega(s-s'))\bigg],\quad \epsilon\rightarrow 0^+.
\end{align*}
To perform the contour integrals, the above function can be broken into positive and negative frequency parts, which have poles at $\omega=\pm i\epsilon$. The result of the limiting integration when $s \geq s'$ is $G(s,s')=s'$, or in general $G(s,s')=\min(s,s')$ from the symmetry between $s$ and $s'$ in this case (and also since $s$ and $s'$ are positive). [Note that this result is valid in the limit as $L$ approaches infinity. Near the opposite boundary, a similar formula for $G(s,s')$ holds, except with $s\mapsto L-s$.]
