Polyakov From Nambu-Goto Directly, for Strings? The following derivation, for a classical relativistic point particle, of the 'Polyakov' form of the action from the 'Nambu-Goto' form of the action, without any tricks - no equations of motion or Lagrange multipliers just a direct set of equalities, is as follows:
$$S = - m \int ds = - m \int \sqrt{-g_{\mu \nu} \dot{X}^{\mu} \dot{X}^{\nu}} d \tau = - m \int \sqrt{- \dot{X}^2}d \tau = \frac{-m}{2} \int \frac{-2\dot{X}^2}{\sqrt{-\dot{X}^2}}d \tau \\ = \frac{1}{2} \int \frac{\dot{X}^2 + \dot{X}^2}{\sqrt{-\dot{X}^2/m^2}}d \tau = \frac{1}{2} \int \frac{\dot{X}^2 - m^2(-\dot{X}^2/m^2)}{\sqrt{-\dot{X}^2/m^2}}d \tau = \frac{1}{2} \int (e^{-1} \dot{X}^2 - e m^2)d \tau$$
Apart from randomly adding $\frac{m^2}{m^2}$ to only one of the $\dot{X}^2$ terms in the second last equality (can anybody explain this without referring to the EOM or LM's?), this derivation is completely straightforward.
Can a similarly straightforward derivation of the Polyakov string action from the Nambu-Goto string action be given, without knowing the Polyakov action in advance?
The best hope comes from reversing the last line of this Wikipedia calculation:
$$
S = \frac{T}{2}\int\mathrm d^2\sigma\sqrt{-h}h^{ab}G_{ab}=\frac{T}{2}\int\mathrm d^2\sigma\frac{2\sqrt{-G}}{h^{cd}G_{cd}}h^{ab}G_{ab}=T\int\mathrm d^2\sigma\sqrt{-G}
$$
but it's so random, unmotivated and unexplained I cannot see it as obvious to do such a calculation. I can loosely motivate adding $\frac{h^{ab}G_{ab}}{h^{cd}G_{cd}}$ by noting $\sqrt{-G}$ is like the general relativity volume element telling us to add in $1 = $ stuff built from what's under the square root over itself, but that's it, the $2$'s are quite random too...
[This is nice but (maybe I'm wrong) I see it as too distinct from what I'm asking].
 A: I) OP is asking for a direct/forward derivation from the Nambu-Goto (NG) action to the Polyakov (P) action (as opposed to the opposite derivation). This is non-trivial since the Polyakov action contains the world-sheet (WS) metric $h_{\alpha\beta}$ with 3 more variables as compared to the Nambu-Goto action.
Although we currently do not have a natural forward derivation of all 3 new variables, we have for 2 of the 3 variables, see section IV below.
II) Let us first say a few words about the derivation of the relativistic point particle
$$ L~:=~\frac{\dot{x}^2}{2e}-\frac{e m^2}{2}\tag{1} $$
from the square root Lagrangian
$$L_0~:=~-m\sqrt{-\dot{x}^2}.\tag{2} $$
Note that OP's derivation does not explain/illuminate the fact that the einbein/Lagrange multiplier
$$ e~>~0\tag{3}$$
can be taken as an independent variable, and not just a trivial renaming of the quantity $\frac{1}{m}\sqrt{-\dot{x}^2}>0$. It is an important property of the Lagrangian (1) that we can vary the einbein/Lagrange multiplier (3)  independently. OP's request to not use Lagrange multipliers seems misguided, and we will not follow this instruction.
III) It is possible to directly/forwardly/naturally derive the Lagrangian (1) with its Lagrange multiplier $e$ from the square root Lagrangian (2) as follows:

*

*Derive the Hamiltonian version of the square root Lagrangian (2) via a (singular) Legendre transformation. This is a straightforward application of the unique Dirac-Bergmann recipe. This leads to momentum variables $p_{\mu}$ and one constraint with corresponding Lagrange multiplier $e$. The constraint reflects world-line reparametrization invariance of the square root action (1). The Hamiltonian $H$ becomes of the form 'Lagrange multiplier times constraint':
$$H~=~\frac{e}{2}(p^2+m^2).\tag{4}  $$
See also e.g. this & this Phys.SE posts.


*The corresponding Hamiltonian Lagrangian reads
$$\begin{align} L_H~=~&p \cdot \dot{x} - H\cr
~=~&p \cdot \dot{x} - \frac{e}{2}(p^2+m^2).\end{align} \tag{5} $$


*If we integrate out the momentum $p_{\mu}$ again (but keep the Lagrange multiplier $e$), the Hamiltonian Lagrangian density (5) becomes the sought-for Lagrangian (1). $\Box$
IV) The argument for the string is similar.

*

*Start with the NG Lagrangian density
$${\cal L}_{NG}~:=~-T_0\sqrt{{\cal L}_{(1)}}, \tag{6}$$
$$\begin{align} {\cal L}_{(1)}~:=~&-\det\left(\partial_{\alpha} X\cdot \partial_{\beta} X\right)_{\alpha\beta}\cr
~=~&(\dot{X}\cdot X^{\prime})^2-\dot{X}^2(X^{\prime})^2~\geq~ 0.  \end{align}\tag{7}$$


*Derive the Hamiltonian version of the NG string via a (singular) Legendre transformation. This leads to momentum variables $P_{\mu}$ and two constraints with corresponding two Lagrange multipliers, $\lambda^0$ and $\lambda^1$, cf. my Phys.SE answer here. The two constraints reflect WS reparametrization invariance of the NG action (6).


*If we integrate out the momenta $P_{\mu}$ again (but keep the two Lagrange multipliers, $\lambda^0$ and $\lambda^1$), the Hamiltonian Lagrangian density for the NG string becomes
$${\cal L}~=~T_0\frac{\left(\dot{X}-\lambda^0 X^{\prime}\right)^2}{2\lambda^1} -\frac{T_0\lambda^1}{2}(X^{\prime})^2,\tag{8}$$
cf. my Phys.SE answer here.


*[As a check, if we integrate out the two Lagrange multipliers, $\lambda^0$ and $\lambda^1$, with the additional assumption that
$$\lambda^1~>~0\tag{9}$$
to avoid a negative square root branch, we unsurprisingly get back the original NG Lagrangian density (6).]


*Eq. (8) is as far as our forward derivation goes. It can be viewed as the analogue of our derivation for the relativistic point particle in section III.


*Now we will cheat and work backwards from the Polyakov Lagrangian density
$$\begin{align} {\cal L}_P~=~&-\frac{T_0}{2}  \sqrt{-h} h^{\alpha\beta}
\partial_{\alpha}X \cdot\partial_{\beta}X\cr
~=~&\frac{T_0}{2} \left\{\frac{\left(h_{\sigma\sigma}\dot{X}- h_{\tau\sigma}X^{\prime}\right)^2}{\sqrt{-h}h_{\sigma\sigma}} - \frac{ \sqrt{-h}}{h_{\sigma\sigma}}(X^{\prime})^2 \right\} .\end{align} \tag{10}$$


*By classical Weyl symmetry, only 2 out of the 3 degrees of freedom in the WS metric $h_{\alpha\beta}$ enter the Polyakov Lagrangian density (10). If we identify
$$ \lambda^0~=~\frac{h_{\tau\sigma}}{h_{\sigma\sigma}}\quad\text{and} \quad\lambda^1~=~\frac{\sqrt{-h}}{h_{\sigma\sigma}}~>~0, \tag{11} $$
then the Lagrangian (8) becomes the Polyakov Lagrangian density (10). $\Box$
A: One method is to notice that given $$S_{NG} = - T \int d \tau d \sigma \sqrt{- h}$$
where $h = \det (h_{ab})$, $h_{ab} = \partial_a X^{\mu} \partial_b X_{\mu}$ it's variation with respect to $X^{\mu}$ is partially worked out as follows
\begin{align}
\delta S_{NG} &= - T \delta \int d \tau d \sigma \sqrt{-h} \\
         &= - T \int d \tau d \sigma \delta \sqrt{-h} \\
         &= - \frac{T}{2} \int d \tau d \sigma \sqrt{-h}h^{ab} \delta h_{ab} \\
         &= - \frac{T}{2} \int d \tau d \sigma \sqrt{-h}h^{ab} \delta (\partial_a X^{\mu} \partial_b X_{\mu}) 
\end{align}
but the last line is what we'd get as the first line from varying the new action
\begin{align}
S_P = - \frac{T}{2} \int d \tau d \sigma \sqrt{-h}h^{ab} (\partial_a X^{\mu} \partial_b X_{\mu}) 
\end{align}
with respect to $X^{\mu}$ where $h_{ab}$ is just an independent variable (metric).
Another method in line is given in section 3.4.1 of Townsend's string notes using Dirac constrained systems in line with the other answer.
