My question is the following: It is usual in the standard textbooks to firstly choose a gauge (usually the conformal gauge) and then extract the equations of motion from the Polyakov action by variating in the string coordinates $X^M$ and the worldsheet metric $\gamma^{ab}$. I am trying to find the equations of motion of the string prior to choosing a gauge. Let me provide with some details.
Suppose that we have the Polyakov action \begin{equation} S_{\text{Pol}}[\gamma^{ab},X^M]=-\frac{T}{2}\int d^2\xi\sqrt{-\gamma}\hspace{0.1cm}\gamma^{ab}\partial_aX^M\partial_bX^Ng_{MN}\hspace{0.2cm}, \end{equation} where capital indices correspond to the target spacetime and small indices to the 2d worldsheet.For the equations of motion, one has to vary the Polyakov action with respect to changes in both $\gamma^{ab}$ and $X^M$. If $\gamma^{ab}\to \gamma^{ab}+\delta\gamma^{ab}$, then \begin{equation*} \begin{split} \delta S_{\text{Pol}}[\gamma^{ab},X^M]&=-\frac{T}{2}\int d^2\xi\hspace{0.1cm}\delta\left(\sqrt{-\gamma}\hspace{0.1cm}\gamma^{ab}\right)\partial_aX^M\partial_bX^Ng_{MN}\\ &=-\frac{T}{2}\int d^2\xi\sqrt{-\gamma}\hspace{0.1cm}\partial_aX^M\partial_bX^Ng_{MN}\delta \gamma^{ab}\\ &-\frac{T}{2}\int d^2\xi\hspace{0.1cm}\gamma^{ab}\partial_aX^M\partial_bX^Ng_{MN}\delta\left(\sqrt{-\gamma}\right).\\ \end{split} \end{equation*} We also have that \begin{equation}\label{eq:1.2} \delta\left(\sqrt{-\gamma}\right)=\frac{1}{2}\sqrt{-\gamma}\gamma^{ab}\delta\gamma_{ab}=-\frac{1}{2}\sqrt{-\gamma}\gamma_{ab}\delta\gamma^{ab}\hspace{0.2cm}, \end{equation} so one gets that \begin{equation*} \begin{split} \delta S_{\text{Pol}}[\gamma^{ab},X^M]&=-\frac{T}{2}\int d^2\xi\sqrt{-\gamma}\hspace{0.1cm}\partial_aX^M\partial_bX^Ng_{MN}\delta \gamma^{ab}\\ &+\frac{T}{4}\int d^2\xi\sqrt{-\gamma}\hspace{0.1cm}\gamma_{ab}\gamma^{cd}\partial_cX^M\partial_dX^Ng_{MN}\delta\gamma^{ab} \end{split} \end{equation*} leading to the equations of motion for the worldsheet metric field, i.e. the Virasoro constraints: \begin{equation}\label{eq:1.3} T_{ab}\equiv g_{MN}\left(\partial_aX^M\partial_bX^N-\frac{1}{2}\gamma_{ab}\gamma^{cd}\partial_cX^M\partial_dX^N\right)=0\hspace{0.2cm}. \end{equation}
Varying now with respect to $X^M$, we have that if $X^M\to X^M +\delta X^M$ then \begin{equation*} \begin{split} S_{\text{Pol}}[\gamma^{ab},X^M]&\to -\frac{T}{2}\int d^2\xi\sqrt{-\gamma}\hspace{0.1cm}\gamma^{ab}\partial_a(X^M+\delta X^M)\partial_b(X^N+\delta X^N)(g_{MN}+\delta g_{MN})\\ &=S_{\text{Pol}}[\gamma^{ab},X^M]-T\int d^2\xi\sqrt{-\gamma}\hspace{0.1cm}\gamma^{ab}g_{MN}\partial_aX^M\partial_b(\delta X^N)\\ &\qquad-\frac{T}{2}\int d^2\xi\sqrt{-\gamma}\hspace{0.1cm}\gamma^{ab}\partial_aX^M\partial_bX^N\delta g_{MN}+\mathcal{O}(\delta^2), \end{split} \end{equation*} since $g_{MN}=g_{NM}$ and $\gamma_{ab}=\gamma_{ba}$. So, \begin{equation*} \begin{split} \delta S_{\text{Pol}}[\gamma^{ab},X^M]&=-T\int d^2\xi\hspace{0.2cm}\partial_b\left\{\sqrt{-\gamma}\hspace{0.1cm}\gamma^{ab}g_{MN}\partial_aX^M\delta X^N\right\}\\ &\qquad+T\int d^2\xi\hspace{0.2cm}\partial_b\left\{\sqrt{-\gamma}\hspace{0.1cm}\gamma^{ab}g_{MN}\partial_aX^M\right\}\delta X^N\\ &\qquad\qquad-\frac{T}{2}\int d^2\xi\sqrt{-\gamma}\hspace{0.1cm}\gamma^{ab}\partial_aX^M\partial_bX^N\hspace{0.1cm}\frac{\partial g_{MN}}{\partial X^{\Lambda}}\hspace{0.1cm}\delta X^{\Lambda}+\mathcal{O}(\delta^2)\hspace{0.1cm}. \end{split} \end{equation*} Renaming some indices gives \begin{equation*} \begin{split} \delta S_{\text{Pol}}[\gamma^{ab},X^M]&=T\int d^2\xi\hspace{0.1cm}\left\{\partial_b(\sqrt{-\gamma}\hspace{0.1cm}\gamma^{ab}g_{MN}\partial_aX^M)-\frac{1}{2}\sqrt{-\gamma}\hspace{0.1cm}\gamma^{ab}\partial_aX^M\partial_bX^{\Lambda}\hspace{0.1cm}\frac{\partial g_{M\Lambda}}{\partial X^{N}}\right\}\delta X^N\\ &\qquad -T\int d^2\xi\hspace{0.2cm}\partial_b\left\{\sqrt{-\gamma}\hspace{0.1cm}\gamma^{ab}g_{MN}\partial_aX^M\delta X^N\right\}+\mathcal{O}(\delta^2)\hspace{0.1cm}. \end{split} \end{equation*} The last term in the above variation gives rise to four surface terms, namely $$\text{st}_1=T\int_{\sigma_i}^{\sigma_f} d\sigma\left\{\sqrt{-\gamma}\hspace{0.1cm}\gamma^{11}g_{MN}\partial_1X^M\delta X^N\right\}_{\tau_i}^{\tau_f},$$ $$\text{st}_2= T\int_{\tau_i}^{\tau_f} d\tau\left\{\sqrt{-\gamma}\hspace{0.1cm}\gamma^{12}g_{MN}\partial_1X^M\delta X^N\right\}_{\sigma_i}^{\sigma_f},$$ $$\text{st}_3=T\int_{\sigma_i}^{\sigma_f} d\sigma\left\{\sqrt{-\gamma}\hspace{0.1cm}\gamma^{21}g_{MN}\partial_2X^M\delta X^N\right\}_{\tau_i}^{\tau_f},$$ $$\text{st}_4= T\int_{\tau_i}^{\tau_f} d\tau\left\{\sqrt{-\gamma}\hspace{0.1cm}\gamma^{22}g_{MN}\partial_2X^M\delta X^N\right\}_{\sigma_i}^{\sigma_f},$$ which have to vanish so that we get the standard equations of motion $$\partial_b(\sqrt{-\gamma}\hspace{0.1cm}\gamma^{ab}g_{MN}\partial_aX^M)-\frac{1}{2}\sqrt{-\gamma}\hspace{0.1cm}\gamma^{ab}\partial_aX^M\partial_bX^{\Lambda}\hspace{0.1cm}\frac{\partial g_{M\Lambda}}{\partial X^{N}}=0.$$ Now comes the confusing part. If we had already chosen a gauge (say the conformal gauge $\gamma^{ab}=\eta^{ab}$) it would be very easy to vanish the surface terms $\text{st}_i$ by choosing boundary conditions for $X^M$. This procedure would also introduce the notion of closed and open strings.
My question is therefore the following: What are the boundary conditions that make the surface terms $\text{st}_i$ vanish and introduce the notion of closed and open (Neumann & Dirichlet) strings in this "pre-gauge-fixed" Polyakov action?
