Integrating $\int_{1}^{\infty}\int_{1}^{\infty}(x+y)^2e^{-(x+y)}dydx$ 
How to evaluate the following double integral
  $$\int_{1}^{\infty}\int_{1}^{\infty}(x+y)^2e^{-(x+y)}dydx\,?$$

Is there an easier way to evaluate this than to brute force it through integration by parts many many times? I don't think polar coordinate transformation would help for this but I could be mistaken.
 A: Let's consider the following integral: $$\int_1^\infty \int_1^\infty e^{-t(x+y)}dxdy=\int_1^\infty e^{-tx}dx\int_1^\infty e^{-ty}dy=\left(\frac{e^{-t}}{t}\right)^2$$
Then the integral in the question appears by taking two derivatives with respect to $t$ and setting it to $1$.
$$\int_1^\infty \int_1^\infty (x+y)^2 e^{-(x+y)}dxdy=\left.\left(\frac{d^2}{dt^2}\frac{e^{-2t}}{t^2}\right)\right|_{t=1}=\frac{18}{e^2}$$
A: @LeBlanc has already given probably the neatest answer, but I want to offer a slight variant that allows us to make sense of what $18$ is doing there. Consider$$\int_1^\infty\int_1^\infty\exp(-tx-uy)dxdy=\frac{\exp(-t-u)}{tu}.$$We can multiply the integrand by $x^2$, $xy$ or $y^2$ by respectively applying $\partial_t^2$, $\partial_t\partial_u$ or $\partial_u^2$. For example,$$\int_1^\infty\int_1^\infty x^2\exp(-tx-uy)dxdy=\left.\frac{\exp(-t-u)}{tu}\left(\left(-1-\frac1t\right)^2+\frac{1}{t^2}\right)\right|_{t=u=1}=\frac{2^2+1}{e^2}=\frac{5}{e^2},$$and similarly with the $y^2$ term. Finally,$$\int_1^\infty\int_1^\infty xy\exp(-tx-uy)dxdy=\left.\frac{\exp(-t-u)}{tu}\left(-1-\frac1t\right)\left(-1-\frac1u\right)\right|_{t=u=1}=\frac{2^2}{e^2}=\frac{4}{e^2}.$$
