The polar equation of a slightly deformed sphere is $r = r(\theta, \phi)$ so that the position of a point on its surface is given by
\begin{equation}\tag{1}\label{e1}
\vec{r} = r(\theta, \phi)\hat{e}_r.
\end{equation}
We follow the physicists' convention of $\theta$ being the co-latitude and $\phi$ being the longitude of a point. A small area element on the surface of such a sphere can be described
by a parallelogram formed of vectors $d\vec{r}_1$ and $d\vec{r}_2$ where
\begin{eqnarray}
d\vec{r}_1 &=& \frac{\partial\vec{r}_1}{\partial\theta}d\theta \tag{2} \\
d\vec{r}_2 &=& \frac{\partial\vec{r}_2}{\partial\phi}d\phi \tag{3},
\end{eqnarray}
the area being $d\vec{A} = d\vec{r}_1 \times d\vec{r}_2$. Now,
\begin{eqnarray}
\frac{\partial\vec{r}_1}{\partial\theta} &=& \frac{\partial r}{\partial\theta}\hat{e}_r + r\frac{\partial\hat{e}_r}{\partial\theta} = \frac{\partial r}{\partial\theta}\hat{e}_r + r\hat{e}_\theta \tag{4} \\
\frac{\partial\vec{r}_2}{\partial\phi} &=& \frac{\partial r}{\partial\phi}\hat{e}_r + r\frac{\partial\hat{e}_r}{\partial\phi} = \frac{\partial r}{\partial\phi}\hat{e}_r + r\sin\theta\hat{e}_\phi \tag{5}
\end{eqnarray}
so that
\begin{equation}\tag{6}\label{e6}
d\vec{A} = \left(r\sin\theta\frac{\partial r}{\partial\theta}\hat{e}_\theta - r\frac{\partial r}{\partial\phi}\hat{e}_\phi - r^2\sin\theta\hat{e}_r\right)d\theta d\phi.
\end{equation}
The surface area of the deformed sphere is
\begin{equation}\tag{7}\label{e7}
A = \int_0^\pi\int_0^{2\pi}\left(1 + \frac{1}{r^2}\left(\frac{\partial r}{\partial\theta}\right)^2 + \frac{1}{r^2\sin^2\theta}\left(\frac{\partial r}{\partial\phi}\right)^2\right)^{1/2}r^2\sin\theta d\theta d\phi.
\end{equation}
Simplifying the integrand using the binomial theorem and retaining only the lowest ordered derivatives,
\begin{equation}\tag{8}\label{e8}
A = \int_0^\pi\int_0^{2\pi}\left(1 + \frac{1}{2r^2}\left(\frac{\partial r}{\partial\theta}\right)^2 + \frac{1}{2r^2\sin^2\theta}\left(\frac{\partial r}{\partial\phi}\right)^2\right)r^2\sin\theta d\theta d\phi.
\end{equation}
We now express $r$ as a Laplace series
\begin{equation}\tag{9}\label{e9}
r(\theta, \phi) = a_0 + \sum_{n \ge 1}a_n Y_n(\theta, \phi).
\end{equation}
where $Y_n$ is the Laplace function. It is a solution of the partial differential equation
\begin{equation}\tag{10}\label{e10}
\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\left\{\sin\theta\frac{\partial Y_n}{\partial\theta}\right\} + \frac{1}{\sin^2\theta}\frac{\partial^2Y_n}{\partial\phi^2} + n(n+1)Y_n = 0.
\end{equation}
We further assume that $a_n \ll a_0$ in equation \eqref{e9}.
From equation \eqref{e9}
$$
\frac{\partial r}{\partial\theta} = \sum_{n \ge 1}a_n\frac{\partial Y_n}{\partial\theta}
$$
so that
\begin{equation}\tag{11}\label{e11}
\left(\frac{\partial r}{\partial\theta}\right)^2 = \sum_{m,n \ge 1}a_m a_n \frac{\partial Y_m}{\partial\theta}\frac{\partial Y_n}{\partial\theta}.
\end{equation}
Similarly
\begin{equation}\tag{12}\label{e12}
\left(\frac{\partial r}{\partial\phi}\right)^2 = \sum_{m,n \ge 1}a_m a_n \frac{\partial Y_m}{\partial\phi}\frac{\partial Y_n}{\partial\phi}.
\end{equation}
We now express the surface area of the deformed sphere as
\begin{equation}\tag{13}\label{e13}
A = A_1 + A_2 + A_3,
\end{equation}
where
\begin{eqnarray}
A_1 &=& \iint_S r^2\sin\theta d\theta d\phi \\
A_2 &=& \frac{1}{2}\iint_S \left(\frac{\partial r}{\partial\theta}\right)^2 \sin\theta d\theta d\phi \\
A_3 &=& \frac{1}{2}\iint_S \frac{1}{\sin\theta}\left(\frac{\partial r}{\partial\phi}\right)^2 d\theta d\phi,
\end{eqnarray}
where $S$ denotes the surface of the deformed sphere. From equation \eqref{e9}
$$
r^2 = a_0^2 + 2a_0\sum_{n \ge 1}a_nY_n + \sum_{m, n \ge 1} a_m a_n Y_m Y_n
$$
so that
$$
A_1 = 4\pi a_0^2 + 2a_0\sum_{n \ge 1}a_n \iint_S Y_n\sin\theta d\theta d\phi + \sum_{m, n \ge 1}a_m a_n \iint_S Y_m Y_n \sin\theta d\theta d\phi.
$$
From the last equation in Article 170 of Todhunter's book "An Elementary Treatise on Laplace's Functions, Lamé's Functions, and Bessel's Functions",
\begin{equation}\tag{17}\label{e17}
\int_0^{2\pi}Y_n d\phi = 2\pi P_n(\cos\theta)P_n\cos(\theta^\prime)
\end{equation}
so that
$$
\iint_S Y_n\sin\theta d\theta d\phi = 2\pi P_n(\cos\theta^\prime)\int_0^\pi P_n(\cos\theta)\sin\theta d\theta
$$
Using the usual substitution $x = \cos\theta$, we get
\begin{equation}\tag{18}\label{e18}
\iint_S Y_n\sin\theta d\theta d\phi = 2\pi P_n(\cos\theta^\prime)\int_{-1}^1P_n(x)dx = 0,
\end{equation}
where the last equality is a consequence of the property of Legendre functions. The last equation of article 187 and the first equation of article 196 of Todhunter's book
can be combined to get
\begin{equation}\tag{19}\label{e19}
\iint_S Y_n Y_m \sin\theta d\theta d\phi = \frac{4\pi}{2n+1}\delta_{mn}.
\end{equation}
Using equations \eqref{e18} and \eqref{e19} in the expression to get $A_1$ we get
\begin{equation}\tag{20}\label{e20}
A_1 = 4\pi a_0^2 + \sum_{n \ge 1}a_n^2 \iint_S Y_n^2\sin\theta d\theta d\phi.
\end{equation}
We now evaluate the sum
\begin{equation}\tag{22}\label{e21}
A_2 + A_3 =
\frac{1}{2}\sum_{m, n \ge 1} a_ma_n
\iint_S\left(\sin\theta \frac{\partial Y_m}{\partial\theta}\frac{\partial Y_n}{\partial\theta} + \frac{1}{\sin\theta}\frac{\partial Y_m}{\partial\phi}\frac{\partial Y_m}{\partial\phi}\right) d\theta d\phi.
\end{equation}
The integral in this equation can be readily evaluated using the partial differential equation \eqref{e10}. We multiply it by $Y_m$ and integrate the result over $S$. Thus, we have
\begin{equation}\tag{22}\label{e22}
\iint_S Y_m\frac{\partial}{\partial\theta}\left(\sin\theta\frac{\partial Y_n}{\partial\theta}\right)d\theta d\phi + \iint_S\frac{Y_m}{\sin\theta}\frac{\partial^2 Y_n}{\partial\phi^2}d\theta d\phi =
-n(n+1)\iint_S Y_m Y_n \sin\theta d\theta d\phi.
\end{equation}
The first integral in the above equation can be expressed as
$$
\int_0^{2\pi}\left(Y_m \sin\theta\frac{\partial Y_n}{d\theta}\Big|_0^{\pi} - \int_0^\pi \sin\theta \frac{\partial Y_n}{\partial\theta}\frac{\partial Y_m}{\partial\theta}d\theta\right)d\phi =
-\iint_S \frac{\partial Y_m}{\partial\theta}\frac{\partial Y_n}{\partial\theta}\sin\theta d\theta d\phi.
$$
Similarly, the second integral can be simplified to
$$
\int_0^\pi\frac{1}{\sin\theta}\left(Y_m\frac{\partial Y_n}{\partial\phi}\Big|_0^{2\pi} - \int_0^{2\pi}\frac{\partial Y_n}{\partial\phi}\frac{\partial Y_m}{\partial\phi}d\phi\right)d\theta =
-\iint_S\frac{1}{\sin\theta}\frac{\partial Y_m}{\partial\phi}\frac{\partial Y_n}{\partial\phi}d\theta d\phi.A
$$
Equation \eqref{e22} thus becomes
\begin{equation}\tag{23}\label{e23}
\iint_S\left(\sin\theta \frac{\partial Y_m}{\partial\theta}\frac{\partial Y_n}{\partial\theta} + \frac{1}{\sin\theta}\frac{\partial Y_m}{\partial\phi}\frac{\partial Y_n}{\partial\phi}\right) d\theta d\phi =
n(n+1)\iint_S Y_m Y_n \sin\theta d\theta d\phi.
\end{equation}
From equations \eqref{e19}, \eqref{e20} and \eqref{e23}, we get
\begin{equation}\tag{24}\label{e24}
A_2 + A_3 = \frac{1}{2}\sum_{n \ge 1}a_n^2 n(n+1)\iint_S Y_n^2 \sin\theta d\theta d\phi.
\end{equation}
From equations \eqref{e13}, \eqref{e20} and \eqref{e24},
\begin{equation}\tag{25}\label{e25}
A = 4\pi a_0^2\left(1 + \sum_{n \ge 1} \frac{a_n^2}{a_0^2} \left(1 + \frac{n(n+1)}{2}\right) \iint_S \frac{Y_n^2}{4\pi} \sin\theta d\theta d\phi\right).
\end{equation}
The volume element of a deformed sphere is $d\tau = d\vec{r}_3 \cdot (d\vec{r}_1 \times d\vec{r}_3)$, where $d\vec{r}_3 = \partial\vec{r}/\partial r = \hat{e}_r dr$. From equation \eqref{e6}, it follows that
\begin{equation}\tag{26}\label{e26}
d\tau = r^2\sin\theta drd\theta d\phi.
\end{equation}
The volume of the deformed sphere is
\begin{equation}\tag{27}\label{e27}
\tau = \iint_S \frac{r^3}{3}\sin\theta d\theta d\phi.
\end{equation}
From equation \eqref{e9}
$$
r^3 = a_0^3 + 3a_0^2\sum_{n \ge 1} a_n Y_n + 3 a_0 \sum_{m, n \ge 1}a_m a_n Y_m Y_n + O(a_n^3).
$$
We ignore terms of $O(a_n^3)$ so that
$$
r^3 = a_0^3 + 3a_0^2\sum_{n \ge 1} a_n Y_n + 3 a_0 \sum_{m, n \ge 1}a_m a_n Y_m Y_n.
$$
Substituting this in \eqref{e27} and using equations \eqref{e18} and \eqref{e19} we get
$$
\tau = \frac{4\pi a_0^3}{3} + a_0 \sum_{n \ge 1} a_n^2 \iint_S Y_n^2 \sin\theta d\theta d\phi.
$$
or
\begin{equation}\tag{28}\label{e28}
\tau = \frac{4\pi a_0^3}{3}\left(1 + 3 \sum_{n \ge 1} \frac{a_n^2}{a_0^2}\iint_S \frac{Y_n^2}{4\pi}\sin\theta d\theta d\phi\right).
\end{equation}
We now find the potential at the origin due to a charge $Q$ on the surface of a slightly deformed spherical drop of fluid. Let the charge be spread uniformly with a density $\sigma$. An element of
area $dA$ will have a charge $\sigma dA$ and will produce a potential $dV$ at the origin given by
$$
dV = \sigma\frac{dA}{r}
$$
From equation \eqref{e6},
$$
dV = \frac{\sigma}{r}\left(1 + \frac{1}{2r^2}\left(\frac{\partial r}{\partial\theta}\right)^2 + \frac{1}{2r^2\sin^2\theta}\left(\frac{\partial r}{\partial\phi}\right)^2\right)r^2\sin\theta d\theta d\phi
$$
or
$$
dV = \sigma\left(1 + \frac{1}{2r^2}\left(\frac{\partial r}{\partial\theta}\right)^2 + \frac{1}{2r^2\sin^2\theta}\left(\frac{\partial r}{\partial\phi}\right)^2\right)r\sin\theta d\theta d\phi
$$
We now express the potential as a sum
\begin{equation}\tag{29}\label{e29}
V = V_1 + V_2 + V_3,
\end{equation}
where
\begin{eqnarray}
V_1 &=& \sigma\iint_S r\sin\theta d\theta d\phi \\
V_2 &=& \frac{\sigma}{2}\iint_S \frac{1}{r}\left(\frac{\partial r}{\partial\theta}\right)^2 \sin\theta d\theta d\phi \\
V_3 &=& \frac{\sigma}{2}\iint_S \frac{1}{r\sin\theta}\left(\frac{\partial r}{\partial\phi}\right)^2 d\theta d\phi.
\end{eqnarray}
From equations \eqref{e9} and \eqref{e18} it is clear that
\begin{equation}\tag{33}\label{e33}
V_1 = 4\pi\sigma a_0.
\end{equation}
Before considering the other two terms, we note that
$$
\frac{1}{r}\left(\frac{\partial r}{\partial\theta}\right)^2 = \frac{1}{a_0}\left(1 - \sum_{n \ge 1}\frac{a_n}{a_0}Y_n\right)\sum_{l,m \ge 1}a_la_m\frac{\partial Y_l}{\partial\theta}\frac{\partial Y_m}{\partial\theta}.
$$
If we ignore the terms in third order of $a_n$,
\begin{equation}\tag{34}\label{e34}
\frac{1}{r}\left(\frac{\partial r}{\partial\theta}\right)^2 = \frac{1}{a_0} \sum_{l,m \ge 1}a_la_m\frac{\partial Y_l}{\partial\theta}\frac{\partial Y_m}{\partial\theta}.
\end{equation}
Similarly,
\begin{equation}\tag{35}\label{e35}
\frac{1}{r}\left(\frac{\partial r}{\partial\phi}\right)^2 = \frac{1}{a_0\sin\theta} \sum_{l,m \ge 1}a_la_m\frac{\partial Y_l}{\partial\phi}\frac{\partial Y_m}{\partial\phi}.
\end{equation}
Using a procedure similar to the one used to get equation \eqref{e24}, we get
\begin{equation}\tag{36}\label{e36}
V_2 + V_3 = \frac{\sigma}{2a_0}\sum_{n \ge 1}a_n^2 n(n+1)\iint_S Y_n^2 \sin\theta d\theta d\phi.
\end{equation}
Therefore,
\begin{equation}\tag{37}\label{e37}
V = 4\pi\sigma a_0\left(1 + \frac{1}{2}\sum_{n \ge 1}\frac{a_n^2}{a_0^2}n(n+1)\iint_S\frac{Y_n^2}{4\pi} \sin\theta d\theta d\phi \right).
\end{equation}
Substituting $Q/A$ for $\sigma$, where $A$ is given by equation \eqref{e25} we get
\begin{equation}\tag{38}\label{e38}
V = \frac{Q}{a_0}\left(1 - \sum_{n \ge 1}\frac{a_n^2}{a_0^2}\iint_S\frac{Y_n^2}{4\pi} \sin\theta d\theta d\phi + O(a_n^4)\right).
\end{equation}
To get this equation we have used the expression
$$
\frac{1}{A} = \frac{1}{4\pi a_0^2}\left(1 - \sum_{n \ge 1}\frac{a_n^2}{a_0^2}\left(1 + \frac{n(n + 1)}{2}\right)\iint_S\frac{Y_n^2}{4\pi}\sin\theta d\theta d\phi\right).
$$
If we write equation \eqref{e9} as
\begin{equation}\tag{39}\label{e39}
r = a_0\left(1 + \sum_{n \ge 1}F_n\right),
\end{equation}
where
\begin{equation}\tag{40}\label{e40}
F_n = \frac{a_n}{a_0}\frac{Y_n}{\sqrt{n-1}}
\end{equation}
then equation \eqref{e38} becomes
\begin{equation}\tag{41}\label{e41}
V = \frac{Q}{a_0}\left(1 - \sum_{n \ge 1}(n - 1)\iint_S\frac{F_n^2}{4\pi} \sin\theta d\theta d\phi\right).
\end{equation}
