Let $U\subsetneqq\mathbb{R}^{n}$ be open, $\varepsilon>0$ and consider the function $$f_{\varepsilon}(x)=\frac{\pi^{-\frac{n}{2}}}{\varepsilon^{n}}\int_{U}\exp\left\{-\left\|\frac{x-y}{\varepsilon}\right\|^{2}\right\}dy$$
I want to prove that $f_{\varepsilon}$ is analytic on $\mathbb{R}^{n}$. Any tips on how I should proceed? I think the integral (over $\mathbb{R}$, that is) is an error function which has the property that it is entire (i.e. analytic)?
Let's write the function in the form
$$f_{\varepsilon}(x) = \frac{\pi^{-n/2}}{\varepsilon^n} \int_{\mathbb{R}^n} \exp \biggl(-\frac{1}{\varepsilon^2}\sum_{k = 1}^n (x_k - y_k)^2\biggr)\cdot \chi_U(y)\,dy = (\varphi_{\varepsilon} \ast \chi_U)(x),$$
where $\varphi$ is a Gaussian kernel, and $\varphi_{\varepsilon}(x) := \varepsilon^{-n}\varphi(x/\varepsilon)$.
Generalities:
We recall that
This criterion is often the most convenient to show that a convolution involving a real-analytic function is again real-analytic. If $g \colon \mathbb{R}^n \to \mathbb{C}$ is the restriction of a holomorphic function $G$ defined on a uniform neighbourhood $B_{\delta}(\mathbb{R}^n) = \{ z \in \mathbb{C}^n : \lvert\operatorname{Im} z_k\rvert < \delta,\, 1 \leqslant k \leqslant n\}$ of $\mathbb{R}^n$ in $\mathbb{C}^n$, and $h \colon \mathbb{R}^n \to \mathbb{C}$ is measurable, then $K \colon B_{\delta}(\mathbb{R}^n) \times \mathbb{R}^n \to \mathbb{C}$ given by
$$K(z,y) := G(z-y)\cdot h(y)$$
is measurable, and if $\delta$ is chosen small enough, $K(z,\,\cdot\,)$ is often integrable over $\mathbb{R}^n$ for all $z \in B_{\delta}(\mathbb{R}^n)$, such that the function
$$F \colon z \mapsto \int_{\mathbb{R}^n} K(z,y)\,dy$$
is continuous on $B_{\delta}(\mathbb{R}^n)$. Then Morera's theorem shows that $F$ is separately holomorphic if Fubini's theorem can be applied, so that we have
$$\int_{\partial \Delta} F(\zeta_1,\dotsc,\zeta_{k-1},z_k,\zeta_{k+1},\dotsc,\zeta_n)\, dz_k = \int_{\mathbb{R}^n} \int_{\partial \Delta} K(\zeta_1,\dotsc, \zeta_{k-1},z_k,\zeta_{k+1},\dotsc,\zeta_n, y)\,dz_k\,dy$$
for every triangle $\Delta$ contained in the disk $\lvert z_k\rvert < \delta$. Since separately holomorphic functions are holomorphic, it then follows that $g\ast h$ is real-analytic. Note: It would suffice to show that $F$ is separately continuous, since the implication "separately holomorphic $\implies$ holomorphic" doesn't require a continuity assumption. But typically, separate continuity of $F$ is not easier to show than continuity.
An alternative way of showing that $F$ is holomorphic is to use Cauchy's integral formula. For every fixed $y$ and $z \in B_{\rho}(\mathbb{R}^n)$ with $0 < \rho < \delta$ we have
$$K(z,y) = \frac{1}{(2\pi i)^n} \int_{\lvert \zeta_1\rvert = \rho}\dotsc\int_{\lvert \zeta_n\rvert = \rho} \frac{K(\zeta,y)}{(\zeta_1 - z_1)\cdot\dotsc\cdot (\zeta_n-z_n)}\,d\zeta_1\,\dotsc\,d\zeta_n,$$
and if Fubini's theorem can be applied, so that we have
\begin{align} \int_{\mathbb{R}^n}&\Biggl(\frac{1}{(2\pi i)^n} \int_{\lvert \zeta_1\rvert = \rho}\dotsc\int_{\lvert \zeta_n\rvert = \rho} \frac{K(\zeta,y)}{(\zeta_1 - z_1)\cdot\dotsc\cdot (\zeta_n-z_n)}\,d\zeta_1\,\dotsc\,d\zeta_n\Biggr)dy\\ &= \frac{1}{(2\pi i)^n} \int_{\lvert \zeta_1\rvert = \rho}\dotsc\int_{\lvert \zeta_n\rvert = \rho} \frac{1}{(\zeta_1 - z_1)\cdot\dotsc\cdot (\zeta_n-z_n)} \Biggl(\int_{\mathbb{R}^n} K(\zeta,y)\,dy\Biggr) d\zeta_1\,\dotsc\,d\zeta_n, \end{align}
then $F$ is recognised as holomorphic, since it is a Cauchy integral.
A third way to establish the holomorphicity of $F$ is differentiation under the integral. If differentiation under the integral is legitimate, and the partial derivatives are continuous, the holomorphicity follows since the Cauchy-Riemann equations hold for the integrand, and thus also for the integral.
If the real-analytic function $g$ has no extension to a holomorphic function defined on a uniform neighbourhood of $\mathbb{R}^n$, the above methods don't work. Then one can try to derive the necessary bounds for the partial derivatives of $g\ast h$ that guarantee that the error in the Taylor approximations tends to $0$ as the order of the approximation grows. But that is usually a lot of work if at all feasible. When possible, the route through the complex domain is preferable.
Back to our specific case:
The crucial point that guarantees the existence of an integrable dominating function is that for all $c,K \in (0,+\infty)$ the function
$$t \mapsto \exp(K\cdot \lvert t\rvert - ct^2)$$
is integrable over $\mathbb{R}$. This gives us a dominating function for the integrand uniform in $z \in B_R(0) \subset \mathbb{C}^n$ for every $R \in (0,+\infty)$ and thus allows us to deduce that
$$z \mapsto \int_{\mathbb{R}^n} \varphi_{\varepsilon}(z-y)\cdot \chi_U(y)\,dy$$
is an entire function, whence $f_{\varepsilon}$ is real-analytic.
We make it a bit more general, and show that $\varphi_{\varepsilon} \ast h$ is an entire function when $h \colon \mathbb{R}^n \to \mathbb{C}$ is a measurable function of at most exponential growth, that is, there are $A,B\in [0,+\infty)$ such that $\lvert h(y)\rvert \leqslant A\exp(B\lVert y\rVert)$ for (almost) all $y$. Then for $z \in \mathbb{C}^n$ with $\lVert z\rVert \leqslant M$ we have
$$\Biggl\lvert \exp\Biggl(-\varepsilon^{-2}\sum_{k = 1}^n (z_k - y_k)^2\Biggr)h(y)\Biggr\rvert \leqslant A\exp(M^2/\varepsilon^2)\cdot \exp \Biggl( \sum_{k = 1}^n \biggl( \frac{2M}{\varepsilon^2} + B\biggr)\lvert y_k\rvert - \varepsilon^{-2}y_k^2\Biggr),\tag{1}$$
where the right hand side is an integrable function independent of $z$. Since $\varphi_{\varepsilon}(z-y)\cdot h(y)$ depends continuously on $z$, the dominated convergence theorem asserts that
$$F_{\varepsilon}(z) = \int_{\mathbb{R}^n} \varphi_{\varepsilon}(z-y)\cdot h(y)\,dy$$
is continuous on $\{ z \in \mathbb{C}^n : \lVert z\rVert \leqslant M\}$. Since $M$ was arbitrary, $F_{\varepsilon}$ is continuous on $\mathbb{C}^n$. Using $(1)$ and the joint measurability of $(z,y) \mapsto \varphi_{\varepsilon}(z-y)\cdot h(y)$, we see that we can apply Fubini's theorem, and by our preferred method (Morera's theorem, Cauchy integral formula) we deduce that $F_{\varepsilon}$ is entire.
We can also see that differentiation under the integral is legitimate: $t \mapsto \lvert t\rvert^k \exp(K\lvert t\rvert - ct^2)$ is integrable for all $k,c,K \in (0,+\infty)$ - and obtain the result by differentiating under the integral.
A direct power series expansion
The particular properties of the exponential function allow us to directly obtain a power series expansion of the integral. Using multi-index notation - for $\alpha \in \mathbb{N}^n$ and $w \in \mathbb{C}^n$, we let $\lvert\alpha\rvert = \alpha_1 + \dotsc + \alpha_n$, $\alpha! = \alpha_1!\cdot \dotsc\cdot \alpha_n!$, $w^{\alpha} = w_1^{\alpha_1}\cdot \dotsc\cdot w_n^{\alpha_n}$ and $\lvert w\rvert^{\alpha} = \lvert w^{\alpha}\rvert$ - we obtain the power series representation
$$e^{w_1 + \dotsc + w_n} = \sum_{\alpha \in \mathbb{N}^n} \frac{w^{\alpha}}{\alpha!}$$
which converges absolutely (and locally uniformly) on all of $\mathbb{C}^n$. The estimate
$$\Biggl\lvert\sum_{\lvert \alpha\rvert \leqslant m} \frac{w^{\alpha}}{\alpha!}\Biggr\rvert \leqslant \sum_{\lvert\alpha\rvert \leqslant m} \frac{\lvert w\rvert^{\alpha}}{\alpha!} \leqslant e^{\lvert w_1\rvert + \dotsc + \lvert w_n\rvert}$$
allows us to interchange summation and integration in
\begin{align} \int_{\mathbb{R}^n} \exp\Biggl(\varepsilon^{-2}\sum_{k = 1}^n 2x_k y_k - y_k^2\Biggr)h(y)\,dy &= \int{\mathbb{R}^n} \Biggl(\sum_{\alpha \in \mathbb{N}^n} \frac{2^{\lvert\alpha\rvert}\varepsilon^{-2\lvert\alpha\rvert} y^{\alpha} x^{\alpha}}{\alpha!}\Biggr) \exp(-\varepsilon^2 \lVert y\rVert^2)h(y)\,dy \\ &= \sum_{\alpha \in \mathbb{N}^n} \frac{2^{\lvert\alpha\rvert} \varepsilon^{-2\lvert\alpha\rvert}}{\alpha!} \Biggl(\int_{\mathbb{R}^n} y^{\alpha} \exp(-\varepsilon^{-2}\lVert y\rVert^2)h(y)\,dy\Biggr)x^{\alpha} \end{align}
by the dominated convergence theorem and $(1)$. Multiplying this with the power series
$$\exp( - \varepsilon^{-2} \lVert x\rVert^2) = \sum_{\beta \in \mathbb{N}^n} \frac{(-1)^{\lvert\beta\rvert} \varepsilon^{-2\lvert\beta\rvert} x^{2\beta}}{\beta!}$$
gives a power series expansion of $f_{\varepsilon}$ converging on all of $\mathbb{R}^n$. The coefficients in
$$f_{\varepsilon}(x) = \sum_{\gamma \in \mathbb{N}^n} c_{\gamma} x^{\gamma}$$
are given by
$$c_{\gamma} = \sum_{2\beta \leqslant \gamma} \frac{(-1)^{\lvert\beta\rvert}2^{\lvert\gamma-2\beta\rvert}\varepsilon^{-2\lvert\gamma - \beta\rvert}}{\beta!\cdot (\gamma - 2\beta)!}\cdot \frac{1}{(\varepsilon\sqrt{\pi})^n}\int_{\mathbb{R}^n} y^{\gamma - 2\beta} \exp(-\varepsilon^{-2}\lVert y\rVert^2)h(y)\,dy.$$