Can we use polar coordinates to evaluate $\int_{x=0}^{1} \int_{y=x}^{1}\left(x^{2}+y^{2}\right) d y d x$

Yeah of course.

$y=1\implies r\sin(\theta)=1$.

So $r=\csc(\theta)$.

So you take the Jacobian and change the limits.

$$\int_{\frac{\pi}{4}}^{\frac{\pi}{2}}\int_{0}^{\csc(\theta)}r^{2}\cdot r\,dr\,d\theta=\frac{1}{3}$$

Yes. If we measure the angle $\theta $ from the positive $y$ axis, the integral takes the form $$I=\int_{\theta=0}^{\pi/4}\int_{r=0}^{\sec\theta}r^2\cdot r\,\mathrm dr\,\mathrm d\theta.$$ The inner integral is $[r^4/4]_0^{\sec\theta}=\frac14\sec^4\theta,$ and so $$I=\int_{\theta=0}^{\pi/4}\frac14\sec^4\theta\,\mathrm d\theta.$$ Now the transformation $t:=\tan\theta$, with $\mathrm dt=\sec^2\theta\,\mathrm d\theta$, yields $$I=\int_{t=0}^1\frac14\sec^2\theta\,\mathrm d t=\int_{t=0}^1\frac14(1+t^2)\,\mathrm d t=\frac14\left[t+\frac{t^3}3\right]_0^1=\frac13.$$