I was reading this article which derives the value for integral of a negative exponent.
I follow the derivation, from circular symmetry, to integration by substitution, and this is the resulting equation: $$(l(a))^2=\pi\int_0^\infty e^{(-ax)} \, dx=\frac \pi a\tag1$$ (where $l(a)$ is some function on $a$.)
And we take the $\sqrt{}$ of equation (1), which produces: $$l(a)=\int_{-\infty}^\infty e^{(-ax^2)} \, dx=\sqrt{\frac \pi a} \text{ for } \operatorname{Re}(a)>0\tag2$$
I understand that the bound became negative infinity because taking the square root gives negative solutions.
However, I don't understand why the exponential factor became $x^2$, and I have no idea where the $\pi$ from equation (1) went. Why did the integral eat the $\pi$?
I think you're conflating two different things in that article.
First, they discuss the integral $\int_0^\infty e^{-a x}\,dx$ and show that it equals $1/a$. This follows by recognizing $\frac{1}{a}e^{-a x}$ as an antiderivative.
Second, they then want to separately compute the integral $I(a)=\int_{-\infty}^\infty e^{-a x^2}\,dx$. To proceed on this second integral, they write $$I(a)^2 = \int_{-\infty}^\infty e^{-a x^2}\,dx\int_{-\infty}^\infty e^{-a y^2}\,dy=\int_{-\infty}^\infty\int_{-\infty}^\infty e^{-a (x^2+y^2)}\,dxdy.$$ To proceed further, they note that the rotational symmetry permits the integration over all $x,y$ to be converted into an integral over $r=\sqrt{x^2+y^2}$ with weight $2\pi r\,dr$. Hence $$~~I(a)^2 = \int_0^\infty 2\pi r e^{-a r^2}\,dr=2\pi\underbrace{\int_0^\infty\frac12e^{-au}du}_{\large u=r^2,\ du=2rdr}=2\pi\left[\frac1{2a} e^{-ar^2}\right]_0^\infty = \frac{\pi}{a}~~$$ which upon taking the square root gives the desired result.