Taylor expansion of a surface

Using dyadic products, we can write the Taylor expansion of a surface as

(1)
\begin{align} \boldsymbol{x}(\theta^1+d\theta^1,\theta^2+d\theta^2)=\boldsymbol{x}(\theta^1,\theta^2)+\boldsymbol{I}\cdot d\boldsymbol{x}+\frac{1}{2}d\boldsymbol{x}\cdot\boldsymbol{H}\cdot d\boldsymbol{x}+\cdots. \end{align}

Here we defined $d\boldsymbol{x},\boldsymbol{I},\boldsymbol{H}$ by

(2)
\begin{align} d\boldsymbol{x}\equiv d\theta^i\boldsymbol{g}_i, \end{align}
(3)
\begin{align} \boldsymbol{I}= \boldsymbol{x}\otimes\nabla=\frac{\partial \boldsymbol{x}}{\partial \theta^i}\otimes\boldsymbol{g}^i=\boldsymbol{g}_i\otimes\boldsymbol{g}^i=\boldsymbol{g}^i\otimes\boldsymbol{g}_i \end{align}

and

(4)
\begin{align} \boldsymbol{H}= \nabla\otimes\left(\boldsymbol{x}\otimes\nabla\right)=\boldsymbol{g}^\alpha\otimes\frac{\partial}{\partial \theta^\alpha}\left(\boldsymbol{g}_i\otimes\boldsymbol{g}^i\right). \end{align}

In this page, we derive those quantities.
Actually, if one uses the second fundamental form and normalized normal vector $\boldsymbol{n}$, this Tayler expansion can also be written in a simpler form as

(5)
\begin{align} =\boldsymbol{x}(\theta^1,\theta^2)+d\theta^i\boldsymbol{g}_i+\frac{1}{2}h_{ij}d\theta^i d\theta^j\boldsymbol{n}+\cdots. \end{align}

Indeed this is more readable representation but we derive this in the second fundamental form and only derive the representation expressed with dyadic product in this page.
Original Taylor expansion
We first try to recall the normal Taylor expansion.
First, we suppose that

(6)
\begin{align} f(x_0+\epsilon)=f(x_0)+a\epsilon+b\epsilon^2+c\epsilon^3+\cdots. \end{align}

Second, we take differentiation of both sides with respect to $\epsilon$ recursively.
(Before that) let us make a convention, such that we write

(7)
\begin{align} f'(x_0+\epsilon)=\frac{df}{d\epsilon} \end{align}

and

(8)
\begin{align} f''(x_0+\epsilon)=\frac{d^2f}{d\epsilon^2}. \end{align}

Then, we can write

(9)
\begin{align} f'=a+2b\epsilon+3c\epsilon^2+\cdots, \end{align}
(10)
\begin{align} f''=2b+6c\epsilon+\cdots, \end{align}

and

(11)
\begin{align} f'''=6c+\cdots. \end{align}

By substituting $\epsilon=0$ to them, we obtain

(12)
\begin{equation} a=f', \end{equation}
(13)
\begin{align} b=\frac{1}{2}f'', \end{align}

and

(14)
\begin{align} c=\frac{1}{6}f'''. \end{align}

When those are substituted to the first equation, we obtain

(15)
\begin{align} f(x_0+\epsilon)=f(x_0)+f'\epsilon+\frac{1}{2}f''\epsilon^2+\frac{1}{6}f'''\epsilon^3+\cdots. \end{align}

This procedure is the original Taylor expansion.
Taylor expansion of a surface
In the following, $\hat{\cdot}$ means that $\cdot$ is fixed.
Basically, a differentiation is defined with respect to one-parameter only. Hence we assume that a "fixed" coordinates

(16)
\begin{align} \left\{\hat{\theta}^1,\hat{\theta}^2\right\} \end{align}

and direction of differentiation

(17)
\begin{align} d\hat{\boldsymbol{x}}(\epsilon)=\epsilon\left(d\hat{\theta^i}\hat{\boldsymbol{g}}_i\right) \end{align}

is given. Note that there is only one parameter, $\epsilon$, in those expressions.
We start with

(18)
\begin{align} \hat{\boldsymbol{x}}(\epsilon)\equiv\boldsymbol{x}(\hat{\theta}^1+\epsilon d\hat{\theta}^1,\hat{\theta}^2+\epsilon d\hat{\theta}^2), \end{align}

where the right-hand side is the position vector that represents the shape of the surface and the left-side hand is a one-parameter vector-valued function that was determined by fixing coordinates and direction of small movement on the surface.
Note that $\hat{\cdot}$ indicates that all the parameters are fixed except $\epsilon$.
Thus, functions became one-parameter functions so we can perform the normal procedure of the Taylor expansion.
By the way, let us suppose

(19)
\begin{align} \hat{\boldsymbol{x}}(\epsilon)=\hat{\boldsymbol{x}}(0)+a\epsilon+b\epsilon^2+c\epsilon^3+\cdots. \end{align}

Here note that

(20)
\begin{align} \hat{\boldsymbol{x}}(0)=\boldsymbol{x}(\hat{\theta}^1,\hat{\theta}^2). \end{align}

By way of trial, we perform the differentiation of the left-hand side with respect to $\epsilon$. Because of the chain-rule, we obtain

(21)
\begin{align} \frac{d\hat{\boldsymbol{x}}}{d\epsilon}=\frac{\partial\boldsymbol{x}}{\partial\theta^1}\frac{d\theta^1}{d\epsilon}+\frac{\partial\boldsymbol{x}}{\partial\theta^2}\frac{d\theta^2}{d\epsilon}=\frac{\partial \boldsymbol{x}}{\partial \theta^1}d\hat{\theta}^1+\frac{\partial \boldsymbol{x}}{\partial \theta^2}d\hat{\theta}^2. \end{align}

Here, we used

(22)
\begin{align} \theta^1(\epsilon)=\hat{\theta}^1+\epsilon d\hat{\theta}^1 \end{align}

and

(23)
\begin{align} \theta^2(\epsilon)=\hat{\theta}^2+\epsilon d\hat{\theta}^2. \end{align}

In addition, because we can write

(24)
\begin{align} \otimes\nabla d\hat{\boldsymbol{x}}=\epsilon\left(\frac{\partial}{\partial \theta^1} d\hat{\theta^1}+\frac{\partial}{\partial \theta^2}d\hat{\theta^2}\right) \end{align}

in genral, in the following calculations we can extensively exploit the following relation:

(25)
\begin{align} \frac{d\hat{\ \ }}{d \epsilon}\epsilon=\otimes\nabla\cdot d\hat{\boldsymbol{x}}. \end{align}

It is not mistaken that $\epsilon$ is not written as $d\epsilon$. This is OK because of the conventions we prepared in the above.
Consequently, we can conclude that differentiation with respect to $\epsilon$ and directional differentiation in $d\hat{\boldsymbol{x}}$ are essentialy identical.
By using those, we can write the first order differentiation as

(26)
\begin{align} \frac{d\hat{\boldsymbol{x}}}{d\epsilon}\epsilon=\boldsymbol{x}\otimes\nabla\cdot d\hat{\boldsymbol{x}}=\boldsymbol{I}\cdot d\hat{\boldsymbol{x}}. \end{align}

The second order differentiation is written as

(27)
\begin{align} \frac{d\hat{\boldsymbol{x}}}{d\epsilon^2}\epsilon^2=d\hat{\boldsymbol{x}}\cdot\nabla\otimes\left(\boldsymbol{I}\cdot d\hat{\boldsymbol{x}}\right). \end{align}

Although it's complicated appearance, because $d\hat{\boldsymbol{x}}$ is only a function of $\epsilon$, we can rewrite

(28)
\begin{align} \frac{d\hat{\boldsymbol{x}}}{d\epsilon^2}\epsilon^2=d\hat{\boldsymbol{x}}\cdot\boldsymbol{H}\cdot d\hat{\boldsymbol{x}}. \end{align}

We are almost in the final stage.
We recursively differentiate Eq. (19).

(29)
\begin{align} \frac{d\hat{\boldsymbol{x}}}{d\epsilon}=a+2b\epsilon+3c\epsilon^2\cdots \end{align}
(30)
\begin{align} \frac{d^2\hat{\boldsymbol{x}}}{d\epsilon^2}=2b+6c\epsilon\cdots \end{align}

Hence, by substituting $\epsilon=0$, we obtain

(31)
\begin{align} a\epsilon=\boldsymbol{I}\cdot d\boldsymbol{x} \end{align}

and

(32)
\begin{align} b\epsilon^2=\frac{1}{2}d\hat{\boldsymbol{x}}\cdot\boldsymbol{H}\cdot d\hat{\boldsymbol{x}}. \end{align}

Finally, we obtain

(33)
\begin{align} \hat{\boldsymbol{x}}(\epsilon)=\hat{\boldsymbol{x}}(0)+\boldsymbol{I}\cdot d\hat{\boldsymbol{x}}(\epsilon)+\frac{1}{2}d\hat{\boldsymbol{x}}(\epsilon)\cdot\boldsymbol{H}\cdot d\hat{\boldsymbol{x}}(\epsilon)+\cdots \end{align}

Without loss of generality, we can replace $d\hat{\boldsymbol{x}}$ with $d\boldsymbol{x}$. This means that, we can generalize the above derivation to arbitrary directions. Hence, we can dettatch $\hat{}$ and obtain

(34)
\begin{align} \boldsymbol{x}(\theta^1+d\theta^1,\theta^2+d\theta^2)=\boldsymbol{x}(\theta^1,\theta^2)+\boldsymbol{I}\cdot d\boldsymbol{x}+\frac{1}{2}d\boldsymbol{x}\cdot\boldsymbol{H}\cdot d\boldsymbol{x}+\cdots. \end{align}

Here, as mentioned previously, we define

(35)
\begin{align} d\boldsymbol{x}= d\theta^i\boldsymbol{g}_i \end{align}
(36)
\begin{align} \boldsymbol{I}= \boldsymbol{x}\otimes\nabla=\frac{\partial \boldsymbol{x}}{\partial \theta^i}\otimes\boldsymbol{g}^i=\boldsymbol{g}_i\otimes\boldsymbol{g}^i=\boldsymbol{g}^i\otimes\boldsymbol{g}_i \end{align}
(37)
\begin{align} \boldsymbol{H}= \nabla\otimes\left(\boldsymbol{x}\otimes\nabla\right)=\boldsymbol{g}^\alpha\otimes\frac{\partial}{\partial \theta^\alpha}\left(\boldsymbol{g}_i\otimes\boldsymbol{g}^i\right). \end{align}

Thus we obtained the Taylor expansion of a surface expressed with dyadic products.