PRELIMINARIES
Q1:
Primarily, the result (3.6.18) is correct (as are all eqns (3.6.14)-(3.6.24)), but it is good that you are sceptical. Showing that the Weyl variations of both sides are the same only guarantees that the two expressions are the same up to an ''integration constant'', which in the current context means they are equal up to an additive Weyl-invariant term. To determine the latter one should write down all local operators that can be constructed out of the fields present, keeping in mind that the scaling dimension (under $\sigma\rightarrow \lambda \sigma$) of all terms should match. You will realise that every local operator you write down satisfying these properties will fail to be Weyl-invariant (remember we need offshell Weyl invariance, so using $k^2=0$ is not allowed), and so this indeed proves that if the Weyl variations of both sides of (3.6.18) are equal then the equality in (3.6.18) is justified.
Q2: I'm not sure I see any problem here. The dimensions on both sides match, and there are three dimensionful quantities to play with: $X^\mu$, $\alpha'$ and $k^\mu$. In fact, the result is exact in $\alpha'$ (since the flat background around which we are doing perturbation theory is exact, modulo vacuum instability). So in fact we could in principle have, e.g., arbitrary powers of $\alpha'k^2$ appearing and there would still be no contradiction.
Eqn (3.6.18) is an offshell statement, we are not allowed to enforce $k^2=0$.
Q3: I will present a (perhaps less simple) proof momentarily; but it is not the same as Polchinski's proof: I will work in a different renormalisation scheme (where in his notation $\gamma=-1$ rather than $\gamma=-2/3$). To add some context, Polchinski subtracts self contractions using geodesic distance,
$$
\Delta(\sigma,\sigma')=\frac{\alpha'}{2}\ln d(\sigma,\sigma')^2,
$$
as indicated in eqn (3.6.5) and (3.6.6). Although this is a quick route to the massless (or low mass) vertex operator Weyl invariance conditions, this normal ordering is impractical and clumsy for generic vertex operators (when you are interested in obtaining something globally well-defined), because at higher mass levels one needs successively higher derivatives of $\Delta(\sigma,\sigma')$ (before taking the limit $\sigma'\rightarrow \sigma$), but these are very cumbersome to work out in a covariant fashion. In particular, the standard way to compute geodesic distance is to work in Riemann normal coordinates, but computing geodesic distance covariantly (in terms of Riemann curvature) is very non-trivial if you need an arbitrarily large number of terms; but you do for generic vertex operators since you need an arbitrarily large number of derivatives of $\Delta(\sigma,\sigma')$, and you can't remove curvature dependence by a choice of coordinates (unless you'd rather store curvature in global information, but I won't elaborate on that unless you happen not to know this and ask for more detail). Also, in geodesic normal ordering the operator-state correspondence is not immediate. So I will discuss a different choice that is much simpler, completely general, covariant (globally defined up to an immaterial phase) in that curvature is stored in local information, and the operator-state correspondence will work in precisely the same way as it does using the conformal normal ordering defined in (2.2.7) in Joe's book (and it also works in the BRST formalism, etc.).
SUMMARY
The $\gamma=-1$ scheme that will be discussed below was invented by Polchinski in his 1988 paper 'Factorization of Bosonic String Amplitudes'; it is called Weyl normal ordering (WNO). Before discussing this in detail I should offer a word of caution: in Weyl normal ordering (and also in dimensional regularisation (DR) but not in geodesic normal ordering) the equation of motion can (at least in the absence of contact terms) be used as an operator equation inside normal ordering, so in fact:
$$
\boxed{:\!\nabla^2X^\mu e^{ik\cdot X}\!:\,\,=0\,}\qquad \textrm{(in WNO & DR)},
$$
which is of course not true in geodesic normal ordering where (3.6.18) holds. What I will rather show is that in Weyl normal ordering:
$$\boxed{
\nabla_a:\!\nabla^aX^\mu e^{ik\cdot X}\!:\,=\,:\!\nabla^aX^\mu \,\nabla_ae^{ik\cdot X}\!:\,+\,\frac{i\alpha'\gamma}{4}k^\mu R_{(2)}:\!e^{ik\cdot X}\!:\,\quad ({\rm with}\,\,\gamma=-1)\,}\quad \textrm{(in WNO)}
$$
So derivatives do not commute with normal ordering in WNO. But you see the sense in which (3.6.18) is encoded here: it is as if the product rule for differentiation has been used. But it hasn't: in WNO you need to keep track of whether you are considering derivatives of normal-ordered operators, or whether you are considering normal-ordered derivatives of operators (e.g., when integrating by parts, etc.) - differentiation and normal ordering do not commute in general.
Two more comments before delving into details:
- It is extremely useful to use a renormalisation scheme (like WNO) for which the equation of motion can be used inside normal ordering (as above), because this is the first indication that one might be able to use powerful conformal field theory techniques (despite the fact that curvature is not assumed to vanish); this is realised in WNO (but not in DR).
WNO is the same as conformal normal ordering (CNO), as defined in (2.2.7) in Polchinski vol.1, but with the additional complication (advantage) that in WNO one also keeps track of Ricci curvature as the base point of the holomorphic frame is translated across the surface. This is accomplished by defining local operators and the subtractions in (2.2.7) using holomorphic normal coordinates (that I will denote by $z_{\sigma_1}$ here and define momentarily) instead of general holomorphic coordinates, such as $z_1,z_2$ in (2.2.7). HNC are discussed in a (hopefully) clear and pedagogical manner in Sec. 2.4 in:
DETAILS
Let $\sigma$ denote the point at which a holomorphic coordinate, $z_{\sigma_1}=z_{\sigma_1}(\sigma)$, is evaluated, and let the subscript, $\sigma_1$, denote that value of $\sigma$ at which the chart is based, i.e. $z_{\sigma_1}(\sigma_1)=0$. We can expose the dependence of these coordinates, $z_{\sigma_1}$, on the base point, $\sigma_1$, by identifying them with holomorphic normal coordinates (HNC), which we define as follows.
Consider a local patch on a surface, and go to conformal gauge where the metric and Ricci scalar at a generic point, $\sigma$, read,
$$
ds^2=\rho(z_{\sigma_1},\bar{z}_{\sigma_1})dz_{\sigma_1}d\bar{z}_{\sigma_1},\qquad R_{(2)}=-4\rho^{-1}\partial_{z_{\sigma_1}}\partial_{\bar{z}_{\sigma_1}}\ln\rho (z_{\sigma_1},\bar{z}_{\sigma_1}).\tag{1}\label{1}
$$
We can always choose $z_{\sigma_1}(\sigma)$ such that at $\sigma=\sigma_1$ where the frame is based the metric is ''as flat as possible'', in particular:
$$
\boxed{
\partial_{z_{\sigma_1}}^n\rho(z_{\sigma_1},\bar{z}_{\sigma_1})\Big|_{\sigma=\sigma_1}=\left\{
\begin{array}
11\quad {\rm for}\quad n=0\\
0\quad {\rm for}\quad n>0.
\end{array}\right.}\tag{2}\label{2}
$$
so that all holomorphic derivatives vanish at $\sigma_1$, but mixed derivatives need not since the Ricci scalar, $R_{(2)}$, cannot be made to vanish by a coordinate choice. Notice that we are only requiring (\ref{2}) to hold at a single point, $\sigma=\sigma_1$. You define Weyl normal ordering by replacing $z_1,z_2$ in the subtractions in (2.2.7) in Joe's book by $z_{\sigma_1}(\sigma),z_{\sigma_1}(\sigma')$ for an operator based at $\sigma_1$. The mode expansions, etc., work out precisely as you would expect, e.g., $i\partial_{z_{\sigma_1}}X(\sigma)=\sum_{n\in\mathbf{Z}}\alpha_n^{(z_{\sigma_1})}(\sigma_1)/z_{\sigma_1}(\sigma)^{n+1}$, etc. Of course, as usual the modes, $\alpha_n^{(z_{\sigma_1})}(\sigma_1)$, depend on the frame (which is why I included the superscript $(z_{\sigma_1})$) as well as the base point $\sigma_1$, even though most people (including myself occasionally) omit this from the notation.
Suppose now that we construct another holomorphic coordinate, $z_{\sigma_1'}$, which is based at $\sigma_1'\equiv \sigma_1+\delta\sigma_1$ (with $\delta\sigma_1$ small). (Note that according to our definitions, $z_{\sigma_1'}(\sigma_1')=0$.) It is a defining property of complex manifolds that if $\sigma$ is a point where two holomorphic charts (associated to coordinates, $z_{\sigma_1},z_{\sigma_1'}$) overlap then they can always be related by a holomorphic transformation at $\sigma$,
\begin{equation}
\begin{aligned}
z_{\sigma_1'}(\sigma)&=f_{\sigma_1'\sigma_1}(z_{\sigma_1}(\sigma))\\
&\simeq z_{\sigma_1}(\sigma)+\delta z_{\sigma_1}(\sigma)
\end{aligned}
\tag{3}\label{3}
\end{equation}
Clearly, since the transition function, $f_{\sigma_1'\sigma_1}$, is holomorphic in $z_{\sigma_1}$ so will $\delta z_{\sigma_1}(z_{\sigma_1}(\sigma))$ be holomorphic in $z_{\sigma_1}$. Let us derive $\delta z_{\sigma_1}(\sigma)$, subject to the requirement that shifts $\sigma_1\rightarrow \sigma_1'$ preserve the HNC gauge slice (\ref{2}).
Such a holomorphic transformation (\ref{3}) induces a change in metric at a generic point $\sigma$ of the form,
$$
\delta\ln\rho(\sigma) = \big(\nabla_{z_{\sigma_1}}\delta z_{\sigma_1}+ \nabla_{\bar{z}_{\sigma_1}}\delta \bar{z}_{\sigma_1}\big)(\sigma)\tag{4}\label{4}
$$
Taking the $(n-1)^{\rm th}$ (with $n\geq2$) holomorphic derivative of both sides and evaluating at $\sigma=\sigma_1$, taking into account that $\delta z_{\sigma_1}$ is holomorphic (and $\delta \bar{z}_{\sigma_1}$ anti-holomorphic) and also (\ref{1}) and (\ref{2}) yields,
\begin{equation}
\begin{aligned}
\partial_{z_{\sigma_1}}^{n-1}\delta\ln\rho(\sigma)\big|_{\sigma=\sigma_1} &= \partial_{z_{\sigma_1}}^{n-1}\big(\nabla_{z_{\sigma_1}}\delta z_{\sigma_1}+ \nabla_{\bar{z}_{\sigma_1}}\delta \bar{z}_{\sigma_1}\big)(\sigma)\big|_{\sigma=\sigma_1}\\
&=\partial_{z_{\sigma_1}}^{n-1}\Big[\big(\partial_{z_{\sigma_1}}\delta z_{\sigma_1}+\delta z_{\sigma_1}\partial_{z_{\sigma_1}}\ln\rho\big)(\sigma)+\big( \partial_{\bar{z}_{\sigma_1}}\delta \bar{z}_{\sigma_1}+\delta \bar{z}_{\sigma_1}\partial_{\bar{z}_{\sigma_1}}\ln\rho\big)(\sigma)\Big]\big|_{\sigma=\sigma_1}\\
&=\partial_{z_{\sigma_1}}^n\delta z_{\sigma_1}(\sigma_1)-\frac{1}{4}\delta \bar{z}_{\sigma_1}(\sigma_1)\partial_{z_{\sigma_1}}^{n-2}\big(-4\rho^{-1}\partial_{z_{\sigma_1}}\partial_{\bar{z}_{\sigma_1}}\ln\rho\big)(\sigma_1)\\
&=\partial_{z_{\sigma_1}}^n\delta z_{\sigma_1}(\sigma_1)-\frac{1}{4}\delta \bar{z}_{\sigma_1}(\sigma_1)\nabla_{z_{\sigma_1}}^{n-2}R_{(2)}(\sigma_1)\\
\end{aligned}\tag{5}\label{5}
\end{equation}
Requiring that the shift in base point, $\sigma_1\rightarrow \sigma_1'$, preserves the gauge slice (\ref{2}) amounts to requiring the left-hand side in (\ref{5}) vanishes. That is,
$$
\partial_{z_{\sigma_1}}^n\delta z_{\sigma_1}(\sigma_1)=\frac{1}{4}\delta \bar{z}_{\sigma_1}(\sigma_1)\nabla_{z_{\sigma_1}}^{n-2}R_{(2)}(\sigma_1)
$$
Multiplying left- and right-hand sides by $z_{\sigma_1}(\sigma)^{n}/n!$ and summing over $n=2,3,\dots$, furthermore implies that,
$$
\sum_{n=2}^{\infty}\frac{1}{n!}\big(\partial_{z_{\sigma_1}}^n\delta z_{\sigma_1}(\sigma_1)\big)z_{\sigma_1}(\sigma)^{n}=\delta \bar{z}_{\sigma_1}(\sigma_1)\frac{1}{4}\sum_{n=2}^{\infty}\frac{1}{n!}\big(\nabla_{z_{\sigma_1}}^{n-2}R_{(2)}(\sigma_1)\big)z_{\sigma_1}(\sigma)^{n}.\tag{6}\label{6}
$$
The left-hand side also equals,
\begin{equation}
\begin{aligned}
\sum_{n=2}^{\infty}\frac{1}{n!}\big(\partial_{z_{\sigma_1}}^n\delta z_{\sigma_1}(\sigma_1)\big)z_{\sigma_1}(\sigma)^{n}&=\sum_{n=0}^{\infty}\frac{1}{n!}\big(\partial_{z_{\sigma_1}}^n\delta z_{\sigma_1}(\sigma_1)\big)z_{\sigma_1}(\sigma)^{n}-\delta z_{\sigma_1}(\sigma_1)-\big(\partial_{z_{\sigma_1}}\delta z_{\sigma_1}(\sigma_1)\big)z_{\sigma_1}(\sigma)\\
&=\delta z_{\sigma_1}(\sigma)-\delta z_{\sigma_1}(\sigma_1)-\big(\partial_{z_{\sigma_1}}\delta z_{\sigma_1}(\sigma_1)\big)z_{\sigma_1}(\sigma),
\end{aligned}\tag{7}\label{7}
\end{equation}
where in going from the first to the second equality we recognised that the sum is just the Taylor expansion of $\delta z_{\sigma_1}(\sigma)$. Substituting (\ref{7}) into the left-hand side of (\ref{6}), rearranging and shifting the summation variable in the curvature term yields,
$$
\delta z_{\sigma_1}(\sigma) =\delta z_{\sigma_1}(\sigma_1)+\big(\partial_{z_{\sigma_1}}\delta z_{\sigma_1}(\sigma_1)\big)z_{\sigma_1}(\sigma)+ \delta \bar{z}_{\sigma_1}(\sigma_1)\frac{1}{4}\sum_{n=1}^{\infty}\frac{1}{(n+1)!}\big(\nabla_{z_{\sigma_1}}^{n-1}R_{(2)}(\sigma_1)\big)z_{\sigma_1}(\sigma)^{n+1}\tag{8}\label{8}
$$
From (\ref{4}) and $\delta\ln\rho(\sigma_1)=0$ we have also that ${\rm Re}(\partial_{z_{\sigma_1}}\delta z_{\sigma_1}(\sigma_1))=0$, so that to leading order in the variation, after adding $z_{\sigma_1}(\sigma)$ to both sides of (\ref{8}) and making use of (\ref{4}) we find that under an infinitesimal shift of base point, $\sigma_1\rightarrow \sigma_1'$, the new holomorphic coordinate, $z_{\sigma_1'}(\sigma)$ is:
$$
z_{\sigma_1'}(\sigma) = e^{i{\rm Im}\partial_{z_{\sigma_1}}\delta z_{\sigma_1}(\sigma_1)}\Big(z_{\sigma_1}(\sigma)+\delta z_{\sigma_1}(\sigma_1)+ \delta \bar{z}_{\sigma_1}(\sigma_1)\frac{1}{4}\sum_{n=1}^{\infty}\frac{1}{(n+1)!}\big(\nabla_{z_{\sigma_1}}^{n-1}R_{(2)}(\sigma_1)\big)z_{\sigma_1}(\sigma)^{n+1}\Big)
\tag{9}\label{9}
$$
The overall phase, $e^{i{\rm Im}\partial_{z_{\sigma_1}}\delta z_{\sigma_1}(\sigma_1)}$, is not determined by the gauge slice, it is not globally defined (the obstruction being the Euler number), in fact it can be ignored provided we always work mod $U(1)$. In particular therefore,
$$
\boxed{
z_{\sigma_1'}(\sigma) = z_{\sigma_1}(\sigma)+\delta z_{\sigma_1}(\sigma_1)+ \delta \bar{z}_{\sigma_1}(\sigma_1)\frac{1}{4}\sum_{n=1}^{\infty}\frac{1}{(n+1)!}\big(\nabla_{z_{\sigma_1}}^{n-1}R_{(2)}(\sigma_1)\big)z_{\sigma_1}(\sigma)^{n+1}\,
}\tag{10}\label{10}
$$
This is the holomorphic transition function $f_{\sigma_1'\sigma_1}(z_{\sigma_1})$ we have been aiming for. It is holomorphic in $\sigma$ (where the chart coordinates, $z_{\sigma_1},z_{\sigma_1'}$, are evaluated), but it is not holomorphic with respect to shifts in the base point, $\sigma_1$.
A derivative with respect to the base point at fixed $\sigma$ can now be rewritten as follows using the chain rule, equation \eqref{10} and its complex conjugate:
\begin{equation}
\begin{aligned}
\frac{\partial}{\partial z_{\sigma_1}(\sigma_1)}\Big|_{\sigma}&=\frac{\partial z_{\sigma_1}(\sigma)}{\partial z_{\sigma_1}(\sigma_1)}\Big|_{\sigma}\frac{\partial}{\partial z_{\sigma_1}(\sigma)}+\frac{\partial \bar{z}_{\sigma_1}(\sigma)}{\partial z_{\sigma_1}(\sigma_1)}\Big|_{\sigma}\frac{\partial}{\partial \bar{z}_{\sigma_1}(\sigma)}\\
&=-\underbrace{\Big(-\frac{\partial}{\partial z_{\sigma_1}(\sigma)}\Big)}_{L_{-1}^{(z_{\sigma_1})}}-\frac{1}{4}\sum_{n=1}^{\infty}\frac{1}{(n+1)!}\big(\nabla_{\bar{z}_{\sigma_1}}^{n-1}R_{(2)}(\sigma_1)\big)\underbrace{\Big(-\bar{z}_{\sigma_1}(\sigma)^{n+1}\frac{\partial}{\partial \bar{z}_{\sigma_1}(\sigma)}\Big)}_{\tilde{L}_n^{(z_{\sigma_1})}}
\end{aligned}
\end{equation}
Perhaps you will recognise the quantities in the parentheses as representations of Virasoro generators (when the central charge vanishes, which is always the case in critical string theory). (You can verify, e.g., that $[L_n,L_m]=(n-m)L_{m+n}$ in case these representations are not familiar.) So in fact we have shown that,
$$
\frac{\partial}{\partial z_{\sigma_1}(\sigma_1)}\Big|_{\sigma}=-L_{-1}^{(z_{\sigma_1})}-\frac{1}{4}\sum_{n=1}^{\infty}\frac{1}{(n+1)!}\big(\nabla_{\bar{z}_{\sigma_1}}^{n-1}R_{(2)}(\sigma_1)\big)\tilde{L}_n^{(z_{\sigma_1})}\tag{11}\label{11}
$$
Tracing back the definitions you will notice that this is the derivative with respect to the variation $\delta z_{\sigma_1}(\sigma_1)=(z_{\sigma_1'}(\sigma)-z_{\sigma_1}(\sigma))|_{\sigma=\sigma_1}$. So this is a ''passive variation'', i.e. we are shifting the frame (parametrised by $\sigma_1\rightarrow \sigma_1'$) keeping the coordinate $\sigma$ fixed. We can instead shift the coordinate keeping the frame fixed. This just introduces a minus sign because as you can verify using the above (using that $\delta \sigma_1=\sigma_1'-\sigma_1$ is infinitesimal and Taylor expanding):
$$
\delta z_{\sigma_1}(\sigma_1)=(z_{\sigma_1'}(\sigma)-z_{\sigma_1}(\sigma))|_{\sigma=\sigma_1}=-(z_{\sigma_1}(\sigma_1')-z_{\sigma_1}(\sigma_1))\equiv -\delta z
$$
Writing the corresponding derivative with respect to the location of a local operator, $:\!\mathscr{A}(\sigma_1)\!:_z$, inserted at $\sigma_1$ as $\partial_z$ (which is the usual notation), and normal ordered in the frame $z_{\sigma_1}$ (that I'm writing now as $z$), what we have shown, according to (\ref{11}), is that:
$$
\boxed{
\,\,\partial_z:\!\mathscr{A}(\sigma_1)\!:_z\,\,=\,\,:\!\partial_z\mathscr{A}(\sigma_1)\!:_z+\frac{1}{4}\sum_{n=1}^{\infty}\frac{1}{(n+1)!}\big(\nabla_{\bar{z}}^{n-1}R_{(2)}(\sigma_1)\big)\tilde{L}_n:\!\mathscr{A}(\sigma_1)\!:_z\,\,
}\tag{12}\label{12}
$$
where, having exposed the detailed quantities on which everything implicitly depends we have now lightened the notation. Notice that:
$$
\boxed{\,\,:\!\partial_z\mathscr{A}(\sigma_1)\!:_z\,=L_{-1}:\!\mathscr{A}(\sigma_1)\!:_z\,}\label{13}\tag{13}
$$
which, in case it is not familiar, I suggest you verify with explicit examples. The Virasoro generators are given by:
$$
L_n=\oint \frac{dz}{2\pi iz}z^{n+2}T(z),
$$
and one uses standard OPEs of the energy-momentum tensor $T(z)$ with the operator $:\mathscr{A}(\sigma_1):_z$. (I'm calling $T(z)$ a tensor because I am thinking of it as the total matter plus ghost energy-momentum tensor. This is necessary because in the above we assumed the total central charge vanishes.)
All this work to derive (\ref{12})! Notice that these results are exact, and they work for arbitrary backgrounds (any matter CFT plus ghost), and they are offshell statements (and so can be also used to cut open and translate handles on Riemann surfaces, see the above references for further details in case of interest).
EXPLICIT CALCULATION
Consider the operator:
$$
\nabla_a:\!\nabla^aX^\mu e^{ik\cdot X}(\sigma_1)\!:_z,
$$
Switching to holomorphic normal coordinates, according to the defining property (\ref{2}) covariant derivatives at $\sigma_1$ reduce to ordinary derivatives, while $\rho(\sigma_1)=1$. That is, keeping the above detailed notation understood but implicit, the above displayed operator is also equal to:
\begin{equation}
\begin{aligned}
\nabla_a:\!\nabla^aX^\mu e^{ik\cdot X}(\sigma_1)\!:_z&=\,\partial_{z}:\!4\partial_{\bar{z}}X^\mu e^{ik\cdot X}(\sigma_1)\!:_z\\
&=\,:\!4\partial_{z}\big(\partial_{\bar{z}}X^\mu e^{ik\cdot X}\big)(\sigma_1)\!:_z+\frac{1}{2}R_{(2)}(\sigma_1)\tilde{L}_1:\!\big(\partial_{\bar{z}}X^\mu e^{ik\cdot X}\big)(\sigma_1)\!:_z\\
&=\,:\!4\partial_{z}\big(\partial_{\bar{z}}X^\mu e^{ik\cdot X}\big)(\sigma_1)\!:_z+\frac{1}{2}R_{(2)}(\sigma_1):\!\big(-\frac{i\alpha'}{2}k^\mu e^{ik\cdot X}\big)(\sigma_1)\!:_z\\
&=\,:\!4\partial_{z}\big(\partial_{\bar{z}}X^\mu e^{ik\cdot X}\big)(\sigma_1)\!:_z-\frac{i\alpha'}{4}k^\mu R_{(2)}(\sigma_1):\!e^{ik\cdot X}(\sigma_1)\!:_z\\
\end{aligned}\tag{14}\label{14}
\end{equation}
where we made use of the above result (\ref{12}) and evaluated the OPE involving $\tilde{L}_n$ explicitly (only $n=1$ is non-vanishing). Consider the first term on the right-hand side. According to (\ref{13}),
\begin{equation}
\begin{aligned}
:\!4\partial_{z}\big(\partial_{\bar{z}}X^\mu e^{ik\cdot X}\big)(\sigma_1)\!:_z\,&=\,L_{-1}:\!4\partial_{\bar{z}}X^\mu e^{ik\cdot X}(\sigma_1)\!:_z\\
&=\,:\!4\big(\partial_{\bar{z}}X^\mu\big) \partial_ze^{ik\cdot X}(\sigma_1)\!:_z\\
&=\,:\!\nabla_aX^\mu \nabla^a e^{ik\cdot X}(\sigma_1)\!:_z\\
\end{aligned}\tag{15}\label{15}
\end{equation}
which follows from direct computation (the $T(z)\cdot \partial_{\bar{z}}X(\sigma_1)$ OPE is non-singular along the contour along which the mode $L_{-1}$ is defined). In the last equality we switched back to real coordinates. So since the chain rule can be used inside normal ordering without obstruction, (\ref{15}) is precisely the statement that the equation of motion can be used as an operator equation inside Weyl normal ordering:
$$
\boxed{\,:\!4\big(\partial_{z}\partial_{\bar{z}}X^\mu\big) e^{ik\cdot X}(\sigma_1)\!:_z\,=\,:\!\big(\nabla^2X^\mu\big) e^{ik\cdot X}(\sigma_1)\!:_z\,=0\,}\tag{16}\label{16}
$$
Substituting (\ref{15}) back into (\ref{14}) and rearranging we learn that:
$$
\boxed{\nabla_a:\!\nabla^aX^\mu e^{ik\cdot X}(\sigma_1)\!:_z\,=\,:\!\nabla_aX^\mu\,\nabla^a e^{ik\cdot X}(\sigma_1)\!:_z+\frac{i\alpha'\gamma}{4}k^\mu R_{(2)}(\sigma_1):\!e^{ik\cdot X}(\sigma_1)\!:_z\quad (\gamma=-1)\,}\tag{17}\label{17}
$$
which is what we set out to show. So you see the sense in which this relation (3.6.18) in Polchinski can be understood, despite the fact that in Weyl normal ordering (\ref{16}) holds.
To complete the story we need to discuss Weyl transformations in WNO.
WEYL TRANSFORMATIONS
In a general CFT, local operators transform as follows,
$$
\boxed{
\,\,:\!\mathscr{A}^{(w_{\sigma_1})}(\sigma_1)\!:_{w_{\sigma_1}}\,=\,\,:\!\mathscr{A}^{(z_{\sigma_1})}(\sigma_1)\!:_{z_{\sigma_1}}-\sum_{n=0}^\infty \big(\varepsilon_nL_n^{z_{(\sigma_1})}+\bar{\varepsilon}_n\tilde{L}_n^{(z_{\sigma_1})}\big):\!\mathscr{A}^{(z_{\sigma_1})}(\sigma_1)\!:_{z_{\sigma_1}}\,\,
}\tag{18}\label{18}
$$
under a holomorphic change of frame,
$$
\boxed{z_{\sigma_1}(\sigma)\rightarrow w_{\sigma_1}(\sigma) = z_{\sigma_1}(\sigma)+\sum_{n=0}^\infty \varepsilon_n z_{\sigma_1}(\sigma)^{n+1}}
\tag{19}\label{19}
$$
To determine $\varepsilon_n$ associated to Weyl transformations in particular, note that in HNC at a generic point $\sigma$ and at $\sigma=\sigma_1$ respectively:
$$
ds^2=\rho(z_{\sigma_1},\bar{z}_{\sigma_1})dz_{\sigma_1}d\bar{z}_{\sigma_1},\qquad \partial_{z_{\sigma_1}}^n\rho(z_{\sigma_1},\bar{z}_{\sigma_1})\big|_{\sigma=\sigma_1}=\delta_{n,0}
$$
Under a Weyl transformation,
$$
ds^2\rightarrow d\hat{s}^2=e^{\delta\phi(\sigma)}\rho(z_{\sigma_1},\bar{z}_{\sigma_1})dz_{\sigma_1}d\bar{z}_{\sigma_1}
$$
where $\delta \phi(\sigma)=\delta \phi(z_{\sigma_1},\bar{z}_{\sigma_1})$.
We want to associate this Weyl transformation to a holomorphic change of frame, $z_{\sigma_1}\rightarrow w_{\sigma_1}(z_{\sigma_1})$, in particular:
$$
\boxed{
z_{\sigma_1}(\sigma)\rightarrow w_{\sigma_1}(\sigma)=z_{\sigma_1}(\sigma)+\delta z_{\sigma_1}(\sigma)
}\tag{20}\label{20}
$$
keeping the metric fixed, so that we search for new holomorphic coordinates, $w_{\sigma_1}$, satisfying:
$$
e^{\delta\phi(\sigma)}\rho(z_{\sigma_1},\bar{z}_{\sigma_1})dz_{\sigma_1}d\bar{z}_{\sigma_1}=\rho(w_{\sigma_1},\bar{w}_{\sigma_1})dw_{\sigma_1}d\bar{w}_{\sigma_1},
$$
i.e.,
$$
\boxed{\delta \phi(\sigma) = \big(\nabla_{z_{\sigma_1}}\delta z_{\sigma_1}+\nabla_{\bar{z}_{\sigma_1}}\delta \bar{z}_{\sigma_1}\big)(\sigma),\qquad \delta z_{\sigma_1}(\sigma_1)=0}\tag{21}\label{21}
$$
where we keep leading order terms in the variation and (at a generic point $\sigma$), $\nabla_z \delta z=\partial_z\delta z+\delta z\,\partial_z\ln\rho$. The second relation in (\ref{21}) specifies that Weyl transformations be transverse to rigid shifts, i.e. $w_{z_{\sigma_1}}(\sigma_1)\equiv z_{\sigma_1}(\sigma_1)\equiv 0$ (even though $\partial_{z_{\sigma_1}}^n\delta z_{\sigma_1}(\sigma)|_{\sigma=\sigma_1}\neq0$).
We wish to compute $\delta z_{\sigma_1}(\sigma)$ in terms of $\delta\phi(\sigma)$. As above, we construct $\delta z(\sigma)$ by Taylor series. Take the $(n-1)^{\rm th}$ derivative of (\ref{21}) for $n\geq2$ and evaluate at $\sigma=\sigma_1$; multiply both sides of the resulting equation by $z_{\sigma_1}(\sigma)^n/n!$ and sum over $n=2,3,\dots$:
\begin{equation}
\begin{aligned}
\sum_{n=2}^\infty\frac{1}{n!}\big(\partial_{z_{\sigma_1}}^{n-1}\delta\phi(\sigma)\big)\big|_{\sigma=\sigma_1}z_{\sigma_1}(\sigma)^n&=\sum_{n=2}^\infty\frac{1}{n!}\Big[\partial_{z_{\sigma_1}}^{n-1}\big(\nabla_{z_{\sigma_1}}\delta z_{\sigma_1}+\nabla_{\bar{z}_{\sigma_1}}\delta \bar{z}_{\sigma_1}\big)(\sigma)\Big]\big|_{\sigma=\sigma_1}z_{\sigma_1}(\sigma)^n\\
&=\sum_{n=2}^\infty\frac{1}{n!}\Big[\partial_{z_{\sigma_1}}^{n-1}\big(\partial_{z_{\sigma_1}}\delta z_{\sigma_1}+\delta z_{\sigma_1}\partial_{z_{\sigma_1}}\ln\rho+\partial_{\bar{z}_{\sigma_1}}\delta \bar{z}_{\sigma_1}+\delta \bar{z}_{\sigma_1}\partial_{\bar{z}_{\sigma_1}}\ln\rho\big)(\sigma)\Big]\big|_{\sigma=\sigma_1}z_{\sigma_1}(\sigma)^n\\
&=\sum_{n=2}^\infty\frac{1}{n!}\Big[\partial_{z_{\sigma_1}}^n\delta z_{\sigma_1}+\partial_{z_{\sigma_1}}^{n-1}\big(\delta z_{\sigma_1}\partial_{z_{\sigma_1}}\ln\rho\big)+\delta \bar{z}_{\sigma_1}\partial_{z_{\sigma_1}}^{n-1}\partial_{\bar{z}_{\sigma_1}}\ln\rho\Big]\big|_{\sigma=\sigma_1}z_{\sigma_1}(\sigma)^n\\
\end{aligned}
\end{equation}
where we expanded the covariant derivatives and took into account that $\delta \bar{z}_{\sigma_1}(\sigma)$ is anti-holomorphic.
Enforcing the HNC gauge slice constraint (\ref{2}) implies the second term vanishes at $\sigma=\sigma_1$. According to (\ref{1}) the last term is proportional to the Ricci scalar, but according to the second relation in (\ref{21}) this term also vanishes. So,
\begin{equation}
\begin{aligned}
\sum_{n=2}^\infty\frac{1}{n!}\big(\partial_{z_{\sigma_1}}^{n-1}\delta\phi(\sigma)\big)\big|_{\sigma=\sigma_1}z_{\sigma_1}(\sigma)^n&=
\sum_{n=2}^\infty\frac{1}{n!}\big(\partial_{z_{\sigma_1}}^n\delta z_{\sigma_1}(\sigma)\big)\big|_{\sigma=\sigma_1}z_{\sigma_1}(\sigma)^n\\
&=
\sum_{n=0}^\infty\frac{1}{n!}\big(\partial_{z_{\sigma_1}}^n\delta z_{\sigma_1}(\sigma)\big)\big|_{\sigma=\sigma_1}z_{\sigma_1}(\sigma)^n-
\delta z_{\sigma_1}(\sigma_1)-
\big(\partial_{z_{\sigma_1}}\delta z_{\sigma_1}(\sigma)\big)\big|_{\sigma=\sigma_1}z_{\sigma_1}(\sigma)\\
&=
\delta z_{\sigma_1}(\sigma)-
\partial_{z_{\sigma_1}}\delta z_{\sigma_1}(\sigma_1)z_{\sigma_1}(\sigma)
\end{aligned}\tag{22}\label{22}
\end{equation}
where we used the definition of Taylor expansion (noting that $z_{\sigma_1}(\sigma_1)=0$) to reconstruct $\delta z_{\sigma_1}(\sigma)$, and that (according to (\ref{21})) $\delta z_{\sigma_1}(\sigma_1)=0$. Using (\ref{21}) again in (\ref{22}), rearranging and substituting the resulting relation back into (\ref{20}) yields (to leading order in the variation):
$$
\boxed{\,\,
w_{\sigma_1}(\sigma)=e^{\frac{1}{2}\delta \phi(\sigma_1)}\big(z_{\sigma_1}(\sigma)+
\sum_{n=1}^\infty\frac{1}{(n+1)!}\big(\nabla_{z_{\sigma_1}}^{n}\delta\phi(\sigma)\big)\big|_{\sigma=\sigma_1}z_{\sigma_1}(\sigma)^{n+1}\big)
\,\,}\tag{23}\label{23}
$$
This relation generates the holomorphic change of frame, $z_{\sigma_1}\rightarrow w_{\sigma_1}(z_{\sigma_1})$, induced by a Weyl transformation, $\rho\rightarrow e^{\delta\phi}\rho$.
We dropped an overall ill-defined phase, $e^{i{\rm Im}\partial_{z_{\sigma_1}}\delta z_{\sigma_1}(\sigma_1)}$, that is not determined by the gauge slice, and used that covariant and ordinary derivatives are equal at $\sigma_1$.
To make contact with Polchinski, in WNO we use,
$$
\boxed{
\Delta(\sigma',\sigma)=\frac{\alpha'}{2}\ln \big|z_{\sigma_1}(\sigma')-z_{\sigma_1}(\sigma)\big|^2
}\tag{24}\label{24}
$$
(instead of (3.6.6)) in (3.6.5) (for local operators $\mathscr{F}(\sigma_1)$), whereas the explicit Weyl variation, $\delta_{\rm W}\Delta(\sigma',\sigma)$, in (3.6.7) should be interpreted as:
$$
\boxed{\delta_{\rm W}\Delta(\sigma',\sigma)=\frac{\alpha'}{2}\ln\Big|\frac{w_{\sigma_1}(\sigma')-w_{\sigma_1}(\sigma)}{z_{\sigma_1}(\sigma')-z_{\sigma_1}(\sigma)}\Big|^2
}\tag{25}\label{25}
$$
with $w_{\sigma_1}(\sigma)$ given in (\ref{23}) (one keeps only linear terms in $\delta\phi$). For general CFT's (3.6.7) is replaced by (\ref{18}) with $\varepsilon_n$ read off from (\ref{19}) and (\ref{23}).
Exercise 1: Show that:
\begin{equation}
\begin{aligned}
\delta_{\rm W}\Delta(\sigma',\sigma)\big|_{\sigma'=\sigma=\sigma_1}&=\frac{\alpha'}{2}\delta\phi(\sigma_1)\\
\partial_{z'}\delta_{\rm W}\Delta(\sigma',\sigma)\big|_{\sigma'=\sigma=\sigma_1}&=\frac{\alpha'}{4}\partial_{z}\delta\phi(\sigma_1)\\
\partial_{z'}\partial_{\bar{z}}\delta_{\rm W}\Delta(\sigma',\sigma)\big|_{\sigma'=\sigma=\sigma_1}&=0
\end{aligned}\tag{26}\label{26}
\end{equation}
which replace (3.6.11) and (3.6.15). (Note that $\delta \phi(\sigma) \equiv 2\delta\omega(\sigma)$.)
Exercise 2: Derive (3.6.16) using that in WNO (3.6.14) equals:
\begin{equation}
\begin{aligned}
V_1&= \frac{g_c}{\alpha'}\int d^2\sigma_1\sqrt{g}\Big\{\big(g^{ab}s_{\mu\nu}+i\epsilon^{ab}a_{\mu\nu}\big):\!\partial_a X^\mu \partial_b X^\nu e^{ik\cdot X}\!:_z+\alpha'\big(\phi-\frac{1}{4}s^\mu_{\phantom{a}\mu}\big)R_{(2)}:\!e^{ik\cdot X}\!:_z\Big\}(\sigma_1)\\
&= \frac{2g_c}{\alpha'}\int d^2z\Big\{\big(s_{\mu\nu}+a_{\mu\nu}\big):\!\partial_zX^\mu \partial_{\bar{z}}X^\nu e^{ik\cdot X}\!:_z+\frac{\alpha'}{4}\big(\phi-\frac{1}{4}s^\mu_{\phantom{a}\mu}\big)R_{(2)}:\!e^{ik\cdot X}\!:_z\Big\}(\sigma_1)
\end{aligned}\tag{27}\label{27}
\end{equation}
Hint: use Exercise 1, and (\ref{12}) when integrating by parts (which includes (\ref{17})). (As noted, WNO corresponds to $\gamma=-1$ in Polchinski's classification (p.105). Also, $\phi$ in (\ref{27}) is the same as in (3.6.14), not to be confused with $\delta \phi(\sigma)$ here.)