For non-relativistic quantum mechanics, the boundary conditions are rather simple to deal with, they are just
\begin{equation} \langle x_1, t_1 \vert x_2, t_2\rangle = \int_{x_1(t_1)}^{x_2(t_2)} \mathcal{D}x(t)e^{\frac{i}{\hbar}S[x(t)]} \end{equation}
And then we can solve the path integral by hand with the boundary conditions that the position must be $x_1$ at time $t_1$ and $x_2$ at time $t_2$, with whatever method is used to solve the integral (since they usually can be solved exactly in some manner).
For quantum field theory, though, boundary conditions seem to usually ignored, and only the vacuum expectation value is calculated
\begin{equation} \langle \Omega \vert \hat{A}\vert\Omega\rangle = \int \mathcal{D}\varphi(x)A(\varphi(x))e^{\frac{i}{\hbar}S[\varphi(x)]} \end{equation}
What would be the generic method to compute an arbitrary transition amplitude from one state to the other? It is usually written as
\begin{equation} \langle \varphi_1\vert \hat{A}\vert\varphi_2\rangle = \int_{\varphi_1(x_1)}^{\varphi_2(x_2)} \mathcal{D}\varphi(x)A(\varphi(x))e^{\frac{i}{\hbar}S[\varphi(x)]} \end{equation}
But the indications on how to do it are rather scarce. I suppose the method that suggests itself the most naturally would be to just express the wavefunction as a function of the vacuum of the theory, $\vert \varphi_1 \rangle = \hat{B} \vert \Omega \rangle$, but this seems to be both a bit contrary to the spirit of the formalism (ideally something that does not refer to wavefunctions) and also possibly problematic with regard to operator ordering (I know that operator ordering can be dealt with stochastic integrals in the action, but I'm not sure how it applies outside of it).
Is there a way to deal with boundary conditions without solving by hand the path integral for a quantum field for some boundary conditions of that field (which, I see, is extremely rarely done)? What would be a simple example for it, for instance for a boundary condition of a free field such that
\begin{equation} \varphi_i(\vec x, t_i) \propto e^{-\frac{(\vec{x} - \vec{x}_i)^2}{2\sigma^2}} \end{equation}
to compute the transition amplitude between two gaussian wavepackets, ideally without referring to wavefunctions?