How to Taylor expand a scalar quantity about a unit sphere? Consider a vector field $\mathbf{v} (r,\theta)$ expressed in the system of (axisymmetric) spherical coordinates with $r$ denoting the radial distance and $\theta$ the polar angle.
We consider a surface $S$ defined by $r = 1 + \epsilon f(\theta)$ where $\epsilon \ll 1$.
We denote by $\mathbf{n}$ the vector normal to that surface.
Specifically,
$$
\mathbb{n} = \mathbb{e}_r - \epsilon f'(\theta) \, \mathbb{e}_\theta \, , 
$$
wherein $\mathbb{e}_r$ and $\mathbb{e}_\theta$ are the basis unit vectors.
Using perturbation analysis, the field $\mathbb{v}$ can be expressed to leading order as
$$
\mathbb{v} = \mathbb{v}^{(0)} + \epsilon \mathbb{v}^{(1)} \, .
$$
Accordingly,
$$
\mathbb{v} \cdot \mathbb{n} = 
v_r^{(0)} + \epsilon\left( v_r^{(1)} - f'(\theta) v_\theta^{(0)} \right) \, 
$$
wherein $v_r := \mathbb{v} \cdot \mathbb{e}_r$ and $v_\theta := \mathbb{v} \cdot \mathbb{e}_\theta$.

I was wondering how one can Taylor expand the dot product $\mathbb{v} \cdot \mathbb{n}$ about the unit sphere $r=1$.
Specifically,

$$
\bigg. \mathbb{v} \cdot \mathbb{n} \bigg|_{r = 1+\epsilon f(\theta)}
= 
[\bigg. \mathbb{v}^{(0)} \cdot \mathbb{e}_r + \epsilon \left( \text{something} \right) ] \bigg|_{r = 1} 
+ \mathcal{O}(\epsilon^2) \, .
$$
Any help is highly appreciated!
Thank you
 A: We may view $\vec{v} \cdot \hat{n}$ as a scalar function $g=g(r,\theta)$, and then expand this function in the argument $r$ by regular Taylor series
$$
\vec{v} \cdot \hat{n} \bigg\vert_{r=1}=g(r,\theta) \bigg\vert_{r=1} =g(1,\theta)+(r-1)\partial_rg(1,\theta)+\mathcal{O}((r-1)^2)
$$
With (a mess of arguments, to keep track of what is a function of what)
$$
g(r,\theta)=v_r^{(0)}(r,\theta)+\varepsilon \left(v_r^{(1)}(r,\theta)-f'(\theta) v_{\theta}^{(0)}(r,\theta) \right) \\
\partial_rg(1,\theta)= \partial_rv_r^{(0)}(1,\theta)+\varepsilon \left(\partial_r v_r^{(1)}(1,\theta)-f'(\theta) \partial_r v_{\theta}^{(0)}(1,\theta) \right) 
$$
Substituting back
$$
g(r,\theta) \bigg\vert_{r=1} =v_r^{(0)}(1,\theta) + \varepsilon v_r^{(1)}(1,\theta)-\varepsilon f'(\theta)v_{\theta}^{(0)} + (r-1)\left( \partial_r v_r^{(0)}(1,\theta)+\varepsilon \partial_rv_r^{(1)}(1,\theta)-\varepsilon f'(\theta)\partial_r v_\theta ^{(0)}\right)
$$
Collecting powers of $\varepsilon$, and setting $(r-1)=\varepsilon f(\theta)$ we have
$$
g(\varepsilon,\theta) \bigg\vert_{r=1} =v_r^{(0)}(1,\theta)+\varepsilon \left[ v_r^{(1)}(1,\theta)-f'(\theta)v_\theta ^{(0)}(1,\theta)+f(\theta) \partial_r v_r^{(0)}(1,\theta) \right] +\mathcal{O}(\varepsilon^2)
$$
Just to check this is reasonable, I chose some arbitrary functions for $v$ and $f$, and plotted the result above versus the actual functions near $r=1$ across a few values of $\theta$:

We see it does what we want: produce a linear approximant in $\varepsilon$ at $r=1$ for arbitrary $\theta$.
