$\renewcommand{\Re}{\operatorname{Re}}$ $\renewcommand{\Im}{\operatorname{Im}}$ $\newcommand{\erf}{\operatorname{erf}}$ $\newcommand{\dag}{\dagger}$ $\newcommand{\const}{\mathrm{const}}$ $\newcommand{\arcsinh}{\operatorname{arcsinh}}$
In the next several lectures we are going to consider Laplace equation in the disk and similar domains and separate variables there but for this purpose we need to express Laplace operator in pthe olar coordinates. Recall that (from Calculus I) polar coordinates are $(r,\theta)$ connected with Cartesian coordinates by $x=r\cos (\theta)$, $y=r\sin(\theta)$ and conversely
\begin{equation*}
\left\{\begin{aligned}
&r=\sqrt{x^2+y^2},\\
&\theta = \arctan \bigl(\frac{y}{x}\bigr);
\end{aligned}\right.
\end{equation*}
surely the second formula is not exactly correct as changing $(x,y)\to (-x,-y)$ does not change the ratio but replaces $\theta $ by $\theta +\pi$ (or $\theta-\pi)$ since $\theta$ is defined modulo $2\pi n$ with $n\in \mathbb{Z}$. It does not really matter as we are interested only in derivatives:
\begin{equation}
r_x=\cos(\theta),\ r_y=\sin(\theta),\
\theta_x = -r^{-1}\sin(\theta),\ \theta_y = r^{-1}\cos(\theta).
\label{eq-6.3.1}
\end{equation}
Exercise 1. Prove (\ref{eq-6.3.1}).
Then by chain rule \begin{equation} \left\{\begin{aligned} &\partial_x = \cos(\theta)\partial_r - r^{-1}\sin(\theta)\partial_\theta,\\ &\partial_y = \sin(\theta)\partial_r + r^{-1}\cos(\theta)\partial_\theta \end{aligned}\right. \label{eq-6.3.2} \end{equation} and therefore \begin{equation*} \Delta = \partial_x^2+ \partial_y^2= \bigl(\cos(\theta)\partial_r - r^{-1}\sin(\theta)\partial_\theta\bigr)^2 + \bigl(\sin(\theta)\partial_r + r^{-1}\cos(\theta)\partial_\theta\bigr)^2 \end{equation*} and after tedious calculations one can get \begin{equation} \Delta = \partial_r^2 +\frac{1}{r}\partial_r +\frac{1}{r^2}\partial_\theta^2. \label{eq-6.3.3} \end{equation}
Exercise 2. Do it.
Instead we want to use a different method requiring much less error prone calculations but more delicate arguments (useful in more complicated cases).
Note first the identity \begin{equation} \iint\Delta u\cdot v\,dA = -\iint \nabla u \cdot \nabla v\,dA, \label{eq-6.3.4} \end{equation} with $dA=dxdy$ and integrals taken over domain $\mathcal{D}$, provided $v=0$ near $\Gamma$ (boundary of $\mathcal{D}$) and integrals are taken over $\mathcal{D}$.
Now let us express the left- and right-hand expression in the polar coordinates. Recall that polar coordinates are orthogonal (i.e. level lines of $r$ (circles) and level lines of $\theta$ (rays from the origin) are orthogonal in the points where they intersect) and the distance $ds$ between two close points can be calculated as \begin{equation} ds^2 =dx^2+dy^2 =dr^2+r^2d\theta^2 \label{eq-6.3.5} \end{equation} and therefore area element is $dA=dx dy= rdrd\theta$.
But what about $\nabla u\cdot \nabla v$? We claim that \begin{equation} \nabla u\cdot \nabla v = u_r v_r +\frac{1}{r^2}u_\theta v_\theta. \label{eq-6.3.6} \end{equation} Indeed, $\nabla u$ is a vector of the different nature than $d\mathbf{s}=(dx,dy)$. They are connected by $du=\nabla u\cdot d\mathbf{s}$ and when we change coordinates $d\mathbf{s}'=Qd\mathbf{s}$ with some matrix $Q$, and since \begin{equation*} du=\nabla u\cdot d\mathbf{s}=\nabla u'\cdot d\mathbf{s}'= \nabla u'\cdot Qd\mathbf{s}= Q^T \nabla u'\cdot d\mathbf{s} \end{equation*} we conclude that $\nabla u'=Q^{T\,-1}\nabla u$ where $^T$ means transposed matrix. Such dual vectors mathematicians call covectors.
Remark 1.
While mathematicians talk about vectors and covectors, physicists often call them covariant and contravariant vectors.
Also there are notions of pseudo-vectors (and pseudo-covectors) and pseudo-scalars which change signs when right-oriented coordinate system is changed to the left-oriented one. F.e. if we restrict ourselves to Cartesian coordinates, vector-product of two vectors is a pseudo-vector, and oriented volume is a pseudo-scalar. Curl of a vector field is a pseudo-vector field. Intensity of magnetic field is a pseudo-vector.
However for more general coordinate systems there are also densities which in addition to usual transformations "reflect" the change of volume.
For us here important is only the difference between vectors and covectors.
Therefore (\ref{eq-6.3.4}) becomes \begin{multline*} \iint \Delta u\cdot v\,r drd\theta = -\iint \bigl(u_r v_r +\frac{1}{r^2}u_\theta v_\theta\bigr)r\,drd\theta=\\ -\iint \bigl(r u_r v_r +\frac{1}{r}u_\theta v_\theta\bigr)\,drd\theta= \iint \bigl(\bigl(r u_r\bigr)_r + \bigl(\frac{1}{r}u_\theta\bigr)_\theta \bigr)v\,drd\theta, \end{multline*} where we integrated by parts. This identity \begin{equation*} \iint r \Delta u\cdot v\,drd\theta = \iint \bigl(\bigl(r u_r\bigr)_r +\bigl(\frac{1}{r}u_\theta\bigr)_\theta \bigr)v\,drd\theta \end{equation*} holds for any $v$, vanishing near $\Gamma$, and therefore we can nix both integration and $v$: \begin{equation*} r \Delta u = \bigl(r u_r\bigr)_r +\bigl(\frac{1}{r}u_\theta\bigr)_\theta . \end{equation*}
Exercise 3. Think about this. Finally we get \begin{equation*} r\Delta u = r^{-1}\bigl(r u_r\bigr)_r +r^{-1}\bigl(r^{-1}u_\theta\bigr)_\theta, \end{equation*} which is exactly (\ref{eq-6.3.3}).
It may look too complicated for polar coordinates but in more general cases this approach is highly beneficial.
Spherical coordinates are $\rho$ (radius), $\phi$ (latitude) and $\theta$ (longitude): \begin{equation*} \left\{\begin{aligned} &x=\rho \sin(\phi)\cos(\theta),\\ &y=\rho \sin(\phi)\sin(\theta)\\ &z=\rho\cos(\phi). \end{aligned}\right. \end{equation*} Conversely \begin{equation*} \left\{\begin{aligned} &\rho=\sqrt{x^2+y^2+z^2},\\ &\phi =\arctan\bigl(\frac{\sqrt{x^2+y^2}}{z}\bigr),\\ &\theta = \arctan \bigl(\frac{y}{x}\bigr); \end{aligned}\right. \end{equation*} and using chain rule and "simple" calculations becomes rather challenging.
Instead we recall that these coordinates are also orthogonal: if we fix $\phi$ and $\theta$ we get rays from origin, which are orthogonal to the spheres which we get if we fix $r$. On the spheres if we fix $\theta$ we get meridians and if we fix $\phi$ we get parallels and those are also orthogonal. Then \begin{equation} ds^2 =dx^2+dy^2+dz^2 =d\rho^2+\rho^2d\phi^2+\rho^2 \sin^2(\phi)d\theta^2, \tag*{(\ref{eq-6.3.5})'}\label{eq-6.3.5'} \end{equation} where $d\rho$, $\rho d\phi$ and $\rho\sin(\phi)d\theta$ are distances along rays, meridians and parallels and therefore the volume element is $dV=dx dy dz= \rho^2\sin(\theta)d\rho d\phi d\theta$.
Therefore \begin{equation} \nabla u\cdot \nabla v = u_\rho v_\rho +\frac{1}{\rho^2}u_\phi v_\phi+ \frac{1}{\rho^2\sin^2(\phi)} u_\theta v_\theta. \tag*{(\ref{eq-6.3.6})'}\label{eq-6.3.6'} \end{equation}
Plugging this into \begin{equation} \iiint\Delta u\cdot v\,dxdydz = -\iiint \nabla u \cdot \nabla v\,dxdydz, \tag*{(\ref{eq-6.3.4})'}\label{eq-6.3.4'} \end{equation} we get \begin{align*} \iiint &\Delta u\cdot v \rho^2\sin(\phi)\,d\rho d\phi d\theta \\ =&-\iiint \Bigl(u_\rho v_\rho +\frac{1}{\rho^2} u_\phi v_\phi + \frac{1}{\rho^2\sin(\phi)} u_\theta v_\theta\Bigr)\rho^2\sin(\phi) \,d\rho d\phi d\theta\\ =&\iiint \Bigr(\bigl(\rho^2\sin(\phi) u_\rho\bigr)_\rho + \bigl(\sin(\phi) u_\phi\bigr)_\phi+ \bigl(\frac{1}{\sin(\phi)}u_\theta\bigr)_\theta \Bigr) v\,d\rho d\phi d\theta. \end{align*} Then we can nix integration and factor $v$: \begin{equation*} \Delta u\cdot \rho^2\sin(\phi)= \bigl(\rho^2\sin(\phi) u_\rho\bigr)_\rho + \bigl(\sin(\phi) u_\phi\bigr)_\phi+ \bigl(\frac{1}{\sin(\phi)}u_\theta\bigr)_\theta \end{equation*} and then \begin{equation*} \Delta u= \frac{1}{\rho^2\sin(\phi)}\Bigl(\bigl(\rho^2\sin(\phi) u_\rho\bigr)_\rho + \bigl(\sin(\phi) u_\phi\bigr)_\phi+ \bigl(\frac{1}{\sin(\phi)}u_\theta\bigr)_\theta \Bigr) \end{equation*}
and finally \begin{equation} \Delta = \partial_\rho^2 + \frac{2}{\rho}\partial_\rho + \frac{1}{\rho^2}\bigl(\partial_{\phi}^2 +\cot(\phi)\partial_\phi\bigr) +\frac{1}{\rho^2\sin^2(\phi)}\partial_{\theta}^2. \label{eq-6.3.7} \end{equation} (compare with (\ref{eq-6.3.3}))
Definition 1. \begin{equation} \Lambda:= \partial_{\phi}^2 +\cot(\phi)\partial_\phi +\frac{1}{\sin^2(\phi)}\partial_{\theta}^2 \label{eq-6.3.8} \end{equation} is a spherical Laplacian (a.k.a. Laplace-Beltrami operator on the sphere).
If the length element is \begin{equation} ds^2 =\sum_{j,k} g_{jk}d q^j dq^k \tag*{(\ref{eq-6.3.5})''}\label{eq-6.3.5''} \end{equation} where $q=(q^1,\ldots,q^n)$ are new coordinates and we prefer to write $dq^j$ rather than $dq_j$ (to half-follow Einstein's notations) and $g_{jk}$ is symmetric matrix ($g^{kj}=g^{jk}$), then \begin{equation} \nabla u\cdot \nabla v = \sum_{j,k} g^{jk}u_{q^j}v_{q^k} \tag*{(\ref{eq-6.3.6})''}\label{eq-6.3.''} \end{equation} where $(g^{jk})$ is an inverse matrix to $(g_{jk})$: $\sum_{k} g^{jk}g_{kl}= \sum_{k} g_{lk}g^{kj}= \delta^j_l$.
Then the volume element is $d V=|\det (g_{jk})|^{\frac{1}{2}}\,dq^1\cdots dq^n$ and \begin{equation} \Delta u = |\det (g_{jk})|^{-\frac{1}{2}}\sum_{j,k} \frac{\partial \ }{\partial q^j} \Bigl( |\det (g_{jk})|^{\frac{1}{2}} g^{jk}\frac{\partial u}{\partial q^k} \Bigr). \label{eq-6.3.9} \end{equation} Indeed, \begin{align*} \iiint \nabla u\cdot \nabla v\,d V =& \iiint \sum_{j,k} g^{jk}u_{q^k}v_{q^j} |\det (g_{jk})|^{\frac{1}{2}} \,dq^1\cdots dq^n\\ =&\iiint \sum_{j,k} \frac{\partial \ }{\partial q^j} \Bigl( |\det (g_{jk})|^{\frac{1}{2}} g^{jk}\frac{\partial u}{\partial q^k} \Bigr)v\,dq^1\cdots dq^n =\\ =&\iiint |\det (g_{jk})|^{-\frac{1}{2}}\sum_{j,k} \frac{\partial \ }{\partial q^j} \Bigl( |\det (g_{jk})|^{\frac{1}{2}} g^{jk}\frac{\partial u}{\partial q^k} \Bigr)v\,dV \end{align*} where transition between the first and second lines is due to integration by parts.
Remark 2.
Elliptic coordinates on $\mathbb{R}^2$ are $(\sigma, \tau)$: \begin{equation} \left\{\begin{aligned} &x=c\cosh (\sigma) \cos(\tau),\\ &y=c\sinh (\sigma) \sin(\tau). \end{aligned}\right. \label{eq-6.3.10} \end{equation} Level lines $\sigma=\const$ are ellipses with foci at $(-c,0)$ and $(c,0)$ and level lines $\tau=\const$ are hyperbolae with the same focai; so we have confocal ellipses and hyperbolae.
These coordinates are not only orthogonal but they are conformal ($ds^2$ is proportional to $d\sigma^2+d\tau^2$) \begin{equation} ds^2= \bigl(\sinh^2(\sigma)+\sin^2(\tau) \bigr)(d\sigma^2+d\tau^2) \label{eq-6.3.11} \end{equation} and therefore \begin{equation} \Delta = \frac{1}{c^2\bigl(\sinh^2(\sigma)+\sin^2(\tau) \bigr)} (\partial_\sigma^2 +\partial_\tau^2 ). \label{eq-6.3.12} \end{equation}
Elliptic cylindrical coordinates in $\mathbb{R}^3$ are obtained by adding $z$ to elliptic coordinates.
Parabolic coordinates on $\mathbb{R}^2$ are $(\sigma, \tau)$: \begin{equation} \left\{\begin{aligned} &x=\sigma \tau,\\ &y=\frac{1}{2}(\sigma^2-\tau^2). \end{aligned}\right. \label{eq-6.3.13} \end{equation} Level lines $\sigma=\const$ and $\tau=\const$ are confocal parabolae.
These coordinates are also conformal \begin{equation} ds^2= (\sigma^2+\tau^2) (d\sigma^2+d\tau^2), \label{eq-6.3.14} \end{equation} and therefore \begin{equation} \Delta = \frac{1}{\sigma^2+\tau^2} (\partial_\sigma^2 +\partial_\tau^2 ). \label{eq-6.3.15} \end{equation}
Remark 3.
In conformal coordinates angles are the same as in Euclidean coordinates.
In 2D conformal coordinates are plentiful and closely connected to "Complex Variables" but in higher dimensions there are very few of them.
See also 3D coordinate systems (which are not conformal):
Three-dimensional parabolic coordinates are obtained by rotating the two-dimensional system about the symmetry axis of the parabolae.
Parabolic cylindrical coordinates in $\mathbb{R}^3$ are obtained by adding $z$ to parabolic coordinates.