While computing the convolution of Gaussian functions of type $\varphi_a(x)=e^{-\pi x^2/a}$, $x\in\mathbb{R}^d$, $a\in\mathbb{R}, a\neq 0$, I wrote without thinking:
$\begin{align} \varphi_{a}*\varphi_{b}\left(x\right)&=\int e^{-\pi y^{2}/a}e^{-\pi\left(x-y\right)^{2}/b}\mathrm{d}y=e^{-\pi x^{2}/b}\int e^{2\pi xy/b}e^{-\pi\left(\frac{1}{a}+\frac{1}{b}\right)y^{2}}\mathrm{d}y\\&=e^{-\pi x^{2}/b}\int e^{-2\pi i\left(ix/b\right)y}e^{-\pi cy^{2}}\mathrm{d}y=\varphi_{b}\left(x\right)\mathcal{F}\left[\varphi_{\frac{1}{c}}\right]\left(\frac{ix}{b}\right)\\&=\varphi_{b}\left(x\right)\left(\frac{ab}{a+b}\right)^{d/2}\varphi_{c}\left(\frac{ix}{b}\right), \end{align}$
where $c=\frac{1}{a}+\frac{1}{b}=\frac{a+b}{ab}$ and $\mathcal{F}$ stands for the Fourier transform. In particular, we have for Gaussian functions that $$\mathcal{F}[\varphi_c](\omega)=c^{d/2}\varphi_{1/c}(\omega).$$
Even if this yields the expected correct result, I am not able to justify this fact. Is it legitimate to argue as before and thus evaluate the Fourier transform at a complex point even if it takes $\mathbb{R}^d$ arguments by definition? Perhaps Paley-Wiener theory or Laplace transform are involved, but I cannot unravel their role.
Yes, the Fourier transform can be studied in the complex domain, provided that the integral defining it converges there. A sufficient condition for $$ \hat f(y) = \int_{\mathbb{R}^d} e^{-2\pi i x\cdot y} f(x)\,dx \tag1 $$ to converge for all $y\in \mathbb{C}^d$ is that $$ |f(x)| e^{A|x|}\in L^\infty(\mathbb{R}^d) \quad \forall \ A>0 \tag2 $$ The Gaussian function satisfies this property.
Since the integrand in (1) is a holomorphic function of $y\in\mathbb{C}^d$, the function $\hat f$ is holomorphic as well.
The properties of Fourier transform such as "products become convolutions", etc, remain true. One way to see this is to look at the computations and notice we don't need $y$ to be real other than to justify convergence; and under the condition (2) convergence is guaranteed regardless.
Another way is to recall that an entire function on $\mathbb{C}^d$ is determined by its values on $\mathbb{R}^d$. So if a relation such as "the Fourier transform of a Gaussian is another Gaussian" holds for all real arguments, then it holds for all complex arguments.