For a fixed $y \in \mathbb{R}$, how do I evaluate $$\int_{\mathbb{R}}e^{-(x+iy)^2}dx$$
It seems like a Gaussian integral. The integration is really over a shift of the real line. Let's look at the function $f(z)=e^{-z^2}e^{i2bz}$ for a real number $b$. Without loss of generality assume that $b \gt 0$. If $b \lt 0$, we draw the contour below the $x$-axis. Since $f$ is entire, the integration of $f$ over the contour 
is $0$. Now along the vertical segments $II$ and $IV$, $|f(z)|=\left|e^{-(\pm R+iy)^2}e^{i2b(R+iy)}\right|=e^{-R^2+y^2}e^{-2by} \le e^{-R^2}e^{\tau^2+2b\tau}$. Thus the integrals $\left|\int_{II}f \right|$ and $\left|\int_{IV}f \right|$ are bounded by $\tau e^{-R^2}e^{\tau^2+2b\tau}$. Now with $\tau$ and $b$ fixed, this goes to $0$ as $R \to \infty$.
Thus, letting $R \to \infty$, the integral along the contour becomes $$\int_{-\infty}^{\infty}e^{-x^2}e^{2bix}dx-e^{\tau^2-2b\tau}\int_{-\infty}^{\infty}e^{-x^2}e^{2i(b-\tau)x}dx=0$$
This holds for all $b$ and $\tau$. Letting $b=\tau$, we have $$\int_{-\infty}^{\infty}e^{-x^2}e^{2bix}dx=e^{-b^2}\sqrt{\pi}$$
Using this we get that $$\int_{\mathbb{R}}e^{-(x+iy)^2}dx=\sqrt{\pi}$$
How do I proceed to evaluate this without using Contour Integration? Hints and comments are greatly appreciated.
Due to the factor $e^{-x^2}$ the function $$y\mapsto f(y) = \int_{\mathbb{R}}e^{-(x+iy)^2}dx $$
is smooth (please check it with usual routine). Moreover since $e^{-x^2}\to 0,~as~|x|\to\infty$ we have
$$ f'(y)=-\int_{\mathbb{R}} 2i(x+iy)e^{-(x+iy)^2}dx = -i \int_{\mathbb{R}} \frac{\partial}{\partial x}\left(e^{-(x+iy)^2}\right)dx = \left[e^{-(x+iy)^2}\right]_{x=-\infty}^{x=\infty} =0$$ Hence $f$ is constant and therefore for all $y\in \Bbb R$
$$\int_{\mathbb{R}}e^{-(x+iy)^2}dx=\int_{\mathbb{R}}e^{-x^2}dx=\sqrt{\pi}$$