According to Azpeitia 1982, under the conditions that $f'''(x)$ is continuous at $a$ and $f'''(x) \neq 0$, the Lagrange remainder for a second order Taylor Series expansion around $a$ can be expressed as:
$$f(x) = f(a) + f'(a)(x-a) + \frac{1}{2}f''(a^\star)(x-a)^2.$$
Where $a^\star$ lies on the interval between $x$ and $a$ and at the limit $x \to a$ we obtain:
$$\lim_{x \to a} \frac{a^\star - a}{x - a} = {3\choose 2}^{-1} = \frac{1}{3},$$
and therefore
$$\lim_{x \to a} a^\star = \frac{x +2a}{3}.$$
Now I am interested in a similar expression for a multivariate function $f:\mathbb{R}^d \to \mathbb{R}$ where a $d$-dimensional input vector $\mathbf{x}$ is approximated around a $d$-dimensional vector $\mathbf{a}$. In this case we have (I believe based on this answer):
$$f(\mathbf{x}) = f(\mathbf{a}) + \nabla_x f(\mathbf{a})(\mathbf{x}-\mathbf{a}) + \frac{1}{2}(\mathbf{x}-\mathbf{a})^T\mathbf{H}(\mathbf{a}^\star)(\mathbf{x}-\mathbf{a}),$$
where $\mathbf{H}$ represents the Hessian matrix of $f$ and is evaluated at $\mathbf{a}^\star$ which lies on the line between $\mathbf{x}$ and $\mathbf{a}$.
My question is can it be shown that the same limit holds in the multivariate case (e.g. $\lim_{\mathbf{x} \to \mathbf{a}} \mathbf{a}^\star = \frac{1}{3}(\mathbf{x} +2\mathbf{a})$)? If not, why not and is there an alternative solution?
Edit based on comments:
For the single variable case: considering $g(t) = f(a + t(x-a))$, we can take a taylor expansion around $t=0$:
$$g(t) = f(a) + f'(a)(x-a) + \frac{1}{2}f''(a + t^\star(x-a))(x-a)^2,$$
where $t^\star \in [0,t]$. We observe that $g(t=0) = f(x=a)$. My guess is that the next step is to find the $t^\star$ such that the full Taylor expansion around $a$ (i.e. $h(x) = \sum_{n=0}^\infty\frac{1}{n!}f^{(n)}(a)(x-a)^n$) is equal to $g(t)$.
$${\exists \; t^\star}{ : g(t) = h(x)}$$
However, I am not too sure how to find this.
Based on the hints in the comments, I have obtained the following solution. As suggested, I begin by solving for the 1-variable case and then extend to the multivariate case.
1-Variable Case
We wish to solve the following equation for $a^\star = a + t^\star(x - a)$:
$$f(x) = f(a) + f'(a)(x-a) + \frac{1}{2}f''(a^\star)(x-a)^2. \tag{1}$$
Note that, assuming the third derivative exists, we can take the first order taylor expansion of $f''(a^\star)$ around $a$, which (for $\bar{a} = a + \bar{t}(a^\star - a)$) can be written as:
$$f''(a^\star) = f''(a) + f'''(\bar{a})(a^\star - a).$$
Plugging this back into $(1)$ we obtain:
$$f(x) = f(a) + f'(a)(x-a) + \frac{1}{2}(f''(a) + f'''(\bar{a})(a^\star - a))(x-a)^2. \tag{2} $$
We can also use the same method to obtain the third-order Taylor expansion of $f(x)$ directly (where $\hat{a} = a + \hat{t}(x - a)$):
$$f(x) = f(a) + f'(a)(x-a) + \frac{1}{2}f''(a)(x-a)^2 + \frac{1}{3!}f'''(\hat{a})(x-a)^3. \tag{3}$$
Setting $(2) = (3)$ we obtain that:
$$\newcommand\underrel[2]{\mathrel{\mathop{#2}\limits_{#1}}} \frac{1}{2}f'''(\bar{a})(a^\star - a)) = \frac{1}{3!}f'''(\hat{a})(x-a)$$ $$\implies \frac{a^\star - a}{x-a} = \frac{1}{3}\frac{f'''(\hat{a})}{f'''(\bar{a})} \underrel{\lim x\to a}{=} \frac{1}{3}$$
We note that this holds because $\hat{a} \underrel{\lim x\to a}{=} a$ and $\bar{a} \underrel{\lim x\to a}{=} a$.
Multivariate Case
We now extend the same reasoning to the multivariate case. We will use some simplified notation. The tensor or partial derivatives of the Hessian $\mathbf{H}$ with respect the $k$ elements of $\mathbf{x}$ can be expressed as: $\mathbf{M}(\mathbf{a^\star}) = (M_1, ... , M_k) = (\frac{\partial \mathbf{H}}{\partial \mathbf{x}_1}(\mathbf{a^\star}), \cdots, \frac{\partial \mathbf{H}}{\partial \mathbf{x}_k}(\mathbf{a^\star}))$. For a vector $\mathbf{v}$ of lenth $k$ we define $\langle\mathbf{v}, \mathbf{M}(\mathbf{a^\star})\rangle = \sum_{i=1}^k v_i M_i$.
In this case we substitute the multivariate version for $(1)$:
$$f(\mathbf{x}) = f(\mathbf{a}) + \nabla_x f(\mathbf{a})(\mathbf{x}-\mathbf{a}) + \frac{1}{2}(\mathbf{x}-\mathbf{a})^T\mathbf{H}(\mathbf{a}^\star)(\mathbf{x}-\mathbf{a}). \tag{4}$$
Again we take the first order Taylor expansion of $\mathbf{H}(\mathbf{a}^\star)$:
$$\mathbf{H}(\mathbf{a}^\star) = \mathbf{H}(\mathbf{a}) + \langle\mathbf{a^\star} - \mathbf{a}, \mathbf{M}(\mathbf{\bar{a}})\rangle$$
And substituting into $(4)$ we obtain:
$$f(\mathbf{x}) = f(\mathbf{a}) + \nabla_x f(\mathbf{a})(\mathbf{x}-\mathbf{a}) + \frac{1}{2}(\mathbf{x}-\mathbf{a})^T\mathbf{H}(\mathbf{a})(\mathbf{x}-\mathbf{a}) + \frac{1}{2}(\mathbf{x}-\mathbf{a})^T\langle\mathbf{a^\star} - \mathbf{a}, \mathbf{M}(\mathbf{\bar{a}})\rangle(\mathbf{x}-\mathbf{a}) \tag{5}$$
Then as we did in equation $(3)$ we take the third-order expansion of $f(\mathbf{x})$:
$$f(\mathbf{x}) = f(\mathbf{a}) + \nabla_x f(\mathbf{a})(\mathbf{x}-\mathbf{a}) + \frac{1}{2}(\mathbf{x}-\mathbf{a})^T\mathbf{H}(\mathbf{a})(\mathbf{x}-\mathbf{a}) + \frac{1}{3!}(\mathbf{x}-\mathbf{a})^T\langle\mathbf{x} - \mathbf{a}, \mathbf{M}(\mathbf{\hat{a}})\rangle(\mathbf{x}-\mathbf{a}) \tag{6}$$
Comparing $(5)$ and $(6)$ and noting that $(\mathbf{x}-\mathbf{a})^T\langle\mathbf{a^\star} - \mathbf{a}, \mathbf{M}(\mathbf{\bar{a}})\rangle(\mathbf{x}-\mathbf{a}) = (\mathbf{a^\star} - \mathbf{a})^T\langle\mathbf{x}-\mathbf{a}, \mathbf{M}(\mathbf{\bar{a}})\rangle(\mathbf{x}-\mathbf{a})$, we obtain:
$$\frac{1}{2}(\mathbf{a^\star} - \mathbf{a})^T\langle\mathbf{x}-\mathbf{a}, \mathbf{M}(\mathbf{\bar{a}})\rangle(\mathbf{x}-\mathbf{a}) = \frac{1}{3!}(\mathbf{x}-\mathbf{a})^T\langle\mathbf{x} - \mathbf{a}, \mathbf{M}(\mathbf{\hat{a}})\rangle(\mathbf{x}-\mathbf{a}) \tag{7}$$
Again, when we limit $\mathbf{x} \to \mathbf{a}$ we find that $\mathbf{\bar{a}} = \mathbf{a}$ and $\mathbf{\hat{a}} = \mathbf{a}$. By comparing elementwise the LHS and RHS of $(7)$ to find that:
$$\frac{(\mathbf{a^\star} - \mathbf{a})_i}{(\mathbf{x} - \mathbf{a})_i} = \frac{1}{3}.$$