What squeeze theorem says is, if a(t)<=b(t)<=c(t) for all t, and a(t) has the same limit as c(t) does when t goes to some number, then b(t) goes to the same limit. This is quite easy to be seen, just thinking two curves in a plane, being closer and closer as x goes to infinity but never intersect, then every curve between these two curves must be "squeezed" and goes to the same value.
Back to the proof. The first part of the proof shows that I(a) goes to what we want to
compute as a goes to infinity. So all we are going to do is to calculate
(or estimate in some sense) I(a) and apply the limit procedure.
Actually, this proof of Gaussian integral uses Fubini's theorem to make this problem easier, just because it is difficult to calculate the integral of e^(x^2), but it is quite easy when calculating xe^(x^2). But where does this x come from?
When changing Cartesian coordinates into polar coordinates, it is not free because one system is just "uneven" to the other. So we need some coefficients for transformation. Luckily the coefficient is exactly r when changing dxdy into drd(theta). This is actually a motivation of this method.
So, calculating I(a) in polar coordinates is reasonable. The inequality comes from that a circle of radius a is contained in the square with side length 2a, and this square is moreover contained in the circle with radius sqrt(2)a. Also, exp() is always positive. So the integrals over these two circles give us a good estimation of I(a), and they are easy to calculate. Now apply squeeze theorem and everything is done.