Exponentiating a differential operator - "Pseudo-Schrodinger" equation I am struggling to follow a calculation presented in the paper Statistical Mechanics of one-dimensional Ginzburg-Landau fields. An analogous calculation is presented inthe thesis A Study of the Transfer Matrix Method for the Classical Statistical Mechanics of One Dimensional Systems , openly accessible, please refer to equation 19.
A solution to the following integral eigenvalue equation is to be found (I use the notation of the first paper, removing constants I believe are inconsequential):
\begin{gather}
   \int \mathrm{d}x_i \exp{ [ -\beta \Delta x f(x_{i+1},x_i) ] } \Psi_n (x_i) 
  = \exp{[-\beta \Delta x \epsilon_n]} \Psi_n (x_{i+1})
\end{gather}
where
$$ f(x_{i+1}, x_i) = a |x_{i+1}|^2 + b |x_{i+1}|^4 + c \Big| \frac{x_{i+1} - x_i}{\Delta x} \Big| ^2 $$
The left-hand side is re-written by performing a Taylor expansion
\begin{gather}
   \int \mathrm{d}x_i \exp{ [ -\beta \Delta x f(x_{i+1},x_i) ] } \big[ \Psi_n (x_{i+1} ) + 
   (x_{i} -x_{i+1}  ) {\Psi} ^{\prime}(x_{i+1}) + \\
   \frac{1}{2} (x_{i} -x_{i+1}  )^2
   {\Psi}^{\prime \prime} (x_{i+1}) +\dots
   \big] \\
  = \exp{ [-\beta \Delta x \big( a |x_{i+1}|^2 + b |x_{i+1}|^4 \big) 
  \\
  \times (1+\frac{1}{4}\frac{\Delta x}{\beta} \frac{\partial ^2}{\partial x^2_{i=1}} ) \Psi_n (x_{i+1}) ]}
\label{taylor}
\end{gather}
Incidentally, I cannot even reproduce this result as I cannot understand where the factor $\sqrt {\pi}$ from the integral
\begin{equation}
\int_{-\infty}^{\infty} \exp \big(-\frac{x^2}{a}\big) x^2  \mathrm{d}x = \frac{1}{2}  \sqrt {\pi} a^{3/2}
\end{equation}
ends up, but that is not my main problem right now.
Going back to the Taylor expansion of the LHS of the integral eigenvalue equation, now a puzzling step is made.
The authors state "formally, the derivative term can be exponentiated", getting to
$$ \exp{[-\beta \Delta x H] } \Psi_n  = \exp{[-\beta \Delta x \epsilon_n]} \Psi_n $$
where
$$ H = -\frac{1}{4} \frac{1}{\beta^2} \frac{\partial ^2}{\partial x^2_{i+1}} +  a |x_{i+1}|^2 + b |x_{i+1}|^4 $$
I do not understand at all. In the second reference I linked, the step is given for granted.
I have read about exponentiating the derivative operator, as in the  definition
\begin{gather}
\exp {[D]} = \sum_{i=0}^{\infty} \frac{D^i}{i!}
\end{gather}
but how that applies to the calculation above, I am not so sure I follow. Ok I see that only even terms contribute to the integral, and $ \exp {[D^2]}$ would exactly pick those. Yet it seems a bit of a liberty to take, I would appreciate it if anybody would clarify why the calculation is rigorously possible.
On a third source The Frenkel Kontorova model  I found an equally puzzling calculation, again involving exponentiting operators.
Given the integral eigenvalue equation
$$ \int_{-\infty}^{\infty} K(u,u') \Psi_n (u') \mathrm{d}u' = \lambda_n \Psi_n (u) $$
the kernel $K$ so defined
$$ K(u,u') = \exp{ \Big\{ -\frac{1}{2} \beta [ V(u) + V(u') + g(u-u')^2 ]  } \Big\}$$
it is claimed that using the operator identity
$$ \int_{-\infty}^{\infty} \mathrm{d}y \exp [-b(x-y)^2] f(y) = (\frac{\pi}{b})^{\frac{1}{2}} \exp \Big( \frac{1}{4b} \frac{\mathrm{d}^2}{\mathrm{d}x^2} \Big) f(x)
$$
the following holds, (why?)
$$\exp{\big[ -\frac{\beta}{2} V(x) \big]} \exp {(\frac{1}{2 \beta g} \frac{\mathrm{d}^2}{\mathrm{d} x^2})} \exp{\big[ - \frac{\beta}{2} V(x) \big]} \Psi_n (x) = \lambda_n \Psi(x) $$
and "combining three exponentials ... into a single one"(??)
$$ \exp{ \Big( \frac{1}{2 \beta g} \frac{\mathrm{d}^2}{\mathrm{d} x^2} -\beta V(x) - \beta W    \Big)} \Psi_n (x) = \lambda_n \Psi(x)
$$
where $W$ is to be defined by Taylor expanding the three exponents (also unclear, if it is allowed to use the property of exponentials for an operator too, where would $W$ come from?).
Would be grateful if anybody shared a hint on both these calculations.
EDIT - CHECKING ON A SIMPLER CASE
I thought I would check on a simpler case, maybe it could help clarifying what I am missing.
I will try to solve the integral equation
\begin{gather}
   \int \mathrm{d}x_i \exp{ [ -\beta  g(x_{i+1},x_i) ] } \Psi_n (x_i) 
  = \exp{[-\beta  \epsilon_n]} \Psi_n (x_{i+1})
\end{gather}
where
$$ g(x_{i+1}, x_i) =   ( x_{i+1} - x_i) ^2 $$
Following the method sketched above and detailed in the linked sources, the integral eigenvalue problem is equivalent to the differential equation
$$ \Psi ^{\prime \prime} = \epsilon_n \Psi$$
with solutions, given the boundary conditions considered in the references, $\Psi (0) = \Psi (\pi) = 0$
$$ \Psi (x) = \sin (kx) $$
with $k=0,1,2, \dots, $ and $\epsilon_n = k_n^2$
is this correct?
Actually I could check that
\begin{gather}
   \int \mathrm{d}x \exp{ [ -\beta  (y-x)^2 ] } \Im [{\exp(ikx)}]
  = \sqrt \frac{\pi}{\beta} \Im [{\exp(ikx)}]
\end{gather}
so $\sin(kx)$ seems to be indeed an eigenfunction of the integral operator, as arrived at by the method and the ordinary differential equation. I have no doubt made some mess with constants.
I am not so sure the eigenvalues are correct though.
On top of which, the question would anyhow remain, what is the basis of such "operator exponentiation" calculation.
 A: answering " clarify why the calculation is rigorously possible."
A bit formally, but here's an argument:
We want to show that $\int da e^{-\beta a^2}F(a)\sim e^{\nabla^2 \over \beta}F(0) $
lets start by writing the first integral as a Fourier integral:
$\int da e^{-\beta a^2}F(a)=\int da e^{-\beta a^2}\int dk e^{ika}\hat F_k$
 (up to factors of $\pi$ etc)
Exchanging limits and completing the square:
$\int dk \int da e^{-\beta (a - i{k\over 2\beta})^2-{k^2\over 4\beta}}\hat F_k$
The $da$ integral gives us a normalization factor that we can ignore ( $ \int da e^{-\beta (a - i{k\over 2\beta})^2 }= \int da e^{-\beta a^2}={Const\over \sqrt\beta}$), and we are left with:
$\int da e^{-\beta a^2}F(a)\sim \int dk e^{-k^2\over 4\beta}\hat F_k$
Now, lets write $e^{\nabla^2 \over \beta}F(x)=e^{\nabla^2 \over \beta}\int dk e^{ikx}\hat F_k=\int dk e^{\nabla^2 \over \beta} e^{ikx}\hat F_k= \int dk e^{-k^2 \over \beta}e^{ikx}\hat F_k$. We pushed the derivative into the integral since we are differenting under the integral sign.
Setting $x=0$ yields 
$\int da e^{-\beta a^2}F(a)\sim \int dke^{-k^2 \over \beta}  \hat F_k\sim e^{\nabla^2 \over \beta}F(x)|_{x=0}$
A: I continued thinking about my question, and I think I found some interesting information which I would like to mention for future reference, should anybody ever find it useful. 
A very basic yet rigorous explanantion of exponentiated operators is given here, https://projecteuclid.org/download/pdf_1/euclid.pcma/1416323532.
A useful example  is given considering the Fourier transform the PDE
$$ \frac{\partial f}{\partial t}(x) - \frac{\partial ^2 f}{\partial x^2}(x) = 0 $$
to
$$ \frac{\partial \hat{f}}{\partial t}(p) +p^2 \hat{f}(p) = 0 $$
whose solution is
$$ \hat{f}(p) = \exp {(-tp^2)} \hat{f}_0(p)$$
and inverting $f$ could be expressed as
$$ f = U^{-1} \exp{(-tM)} \, U f_0$$
where $U$ stands for the Fourier transform and $M$ the multiplication by $p^2$, from which the relationship to the Green's function is apparent.
The reference made me also understand the relationship to propagators, as hinted out in user619894 's reference to Schulman's book. The following answer is also very interesting, Meaning of Exponential map.
With regards to the answers I received, there are at best fragmentary and partial, as only one question is covered at best. Yet, even if I were not the gentleman I am, I think it is just fair to award the bounty, thank you very much for your help. 
