First principle of stationary action
Consider a real Klein-Gordon scalar field $\phi$ living in a $D$ dimensional flat spacetime. The field is considered off shell (the on shell condition is defined below). Suppose for simplicity that its action on an arbitrary region of spacetime $\Omega$ is
\begin{equation}\tag{1}
S = \int_{\Omega} \frac{1}{2} \big((\partial_a \, \phi )(\partial^a \, \phi) – m^2 \phi^2 \big) d^D x.
\end{equation}
The on shell field is defined as the one which render the action stationary under an arbitrary compactly supported variation of the field. The variation $\delta \phi$ is an arbitrary smooth function with compact support (it is not necessarily analytic). It is vanishing on the boundary $\partial \, \Omega$, and all its derivatives are also vanishing there ; $\delta \phi = 0$ and $\partial_a \, \delta \phi = 0$ on $\partial \, \Omega$. An arbitrary variation of the field induces a variation of its action :
\begin{align}
\delta S &= \int_{\Omega} \big( (\partial_a \, \phi )(\partial^a \, \delta \phi) – m^2 \phi \, \delta \phi \big) d^D x \\[18pt]
&= \int_{\Omega} \partial^a \big( (\partial_a \phi) \, \delta \phi \big) \, d^D x – \int_{\Omega} \big( \partial^a \, \partial_a \phi + m^2 \phi \big) \, \delta \phi \; d^D x. \tag{2}
\end{align}
The first integral gives a surface term, by virtue of the Gauss theorem. It's vanishing if $\delta \phi = 0$ on $\partial \, \Omega$. Since $\delta \phi$ is arbitrary inside the bulk of $\Omega$, we get the Klein-Gordon equation, which defines the on shell condition :
\begin{equation}
\partial^a \, \partial_a \phi + m^2 \phi = 0. \tag{3}
\end{equation}
This is all fine with the usual variational principle. However, to solve the on shell differential equation (i.e the equation of motion), we need some proper boundary conditions that should be imposed on the scalar field. Obviously, they should be compatible with the equation of motion. Without them, the equation of motion cannot be solved. What is the "law" that defines the boundary conditions to be imposed on the field ?
Second principle of stationary action
(hypothetical method to find the boundary conditions on the field)
Now consider an on shell field $\phi$ with some unknown boundary conditions on $\partial \, \Omega$. An arbitrary small variation of the boundary conditions induces a variation of the field ; $\phi' = \phi + \delta \phi$, which is still on shell. In this case, the variation $\delta \phi$ and its derivatives do not necessarily vanish on the boundary ! ($\delta \phi$ is not anymore of compact support). The change of boundary conditions also produces a change in the action :
\begin{equation}\tag{4}
\delta S = \int_{\Omega} \partial^a \big( (\partial_a \phi) \; \delta \phi \big) \, d^D x – \int_{\Omega} \big( \partial^a \, \partial_a \phi + m^2 \phi \big) \, \delta \phi \; d^D x.
\end{equation}
Since the field is on shell, the equation of motion is satisfied in the bulk and the second integral vanishes. We now get a surface integral :
\begin{equation}\tag{5}
\delta S = \int_{\partial \, \Omega} (\partial_a \phi) \, \delta \phi \; d\sigma^a,
\end{equation}
where $d\sigma^a$ are the components of the outward boundary normal. Lets suppose that the action is still stationary under the variation of the boundary conditions of an on shell field. The condition $\delta S = 0$ then imposes
\begin{equation}\tag{6}
(d\sigma^a \; \partial_a \phi) \, \delta \phi = 0,
\end{equation}
everywhere on the boundary $\partial \, \Omega$ (I'm not sure this is right, since the surface integral is a flux. Maybe it is just the integral which vanishes). This suggest two choices :
\begin{align}\tag{7}
\delta \phi &= 0 \; \text{(Dirichlet conditions),}
&&\text{or}
&d\sigma^a \; \partial_a \phi &= 0 \; \text{(Neumann conditions).}
\end{align}
So to summarize: I use the stationary action principle to get the field equations, and then use the principle again but now together with the field equations in order to see what are the possible boundary conditions.
Now, the question is this :
Do the previous procedure actually make sense ?
How can we make the boundary conditions more precise, in details ?
And more specifically, how should we translate the Dirichlet conditions above ; $\delta \phi = 0$ on the boundary $\partial \, \Omega$ ? I'm unable to make sense of this part.
Take note that the arbitrary region of spacetime $\Omega$ and its boundary $\partial \, \Omega$ are fixed here, and there is no variation on coordinates (which are fixed). The boundary conditions that I'm talking about refer to the field configuration on $\partial \, \Omega$, which is a closed hypersurface in spacetime, enclosing the arbitrary region $\Omega$.
What are your opinion on this hypothetical (unconventional ?) application of the stationary action principle ?
EDIT: Please, use the same variables (i.e. a scalar field) in your answer, to talk about "boundary conditions" on $\partial \Omega$ of a field in spacetime, instead of "initial conditions". To me, there's a huge distinction between "field boundaries" and "initial conditions".
Very important: Take note that I may be using the "Nature" Hamilton-Jacobi action and not the "observer" Euler-Lagrange action (I'm not sure yet), as defined in this paper :
https://arxiv.org/abs/1203.2736
As a reference to this question, see section 2 (page 4) of the following paper from Padmanabhan:
Best Answer
Let us here for simplicity consider point mechanics with generalized positions $q^k(t)$ defined on a time interval $[t_i,t_f]\subseteq \mathbb{R}$. The generalization to field theory with fields $\phi^{\alpha}(x)$ defined on a spacetime region $\Omega\subseteq\mathbb{R}^D$ is straightforward.$^1$
Given an (off-shell) action functional $$ I[q]~=~\int_{t_i}^{t_f} \! dt~L,\tag{A}$$ it seems that OP in the first half of his post mainly confirms that the functional/variational derivative $$ \frac{\delta I}{\delta q^k} \tag{B}$$ (if it exists!) does not depend on the choice of boundary conditions (BCs): It's always given by the Euler-Lagrange (EL) formula $$ EL_k~=~\frac{\partial L}{\partial q^k} -\frac{d}{dt}\frac{\partial L}{\partial \dot{q}^k}+\left(\frac{d}{dt}\right)^2 \frac{\partial L}{\partial \ddot{q}^k} - \ldots, \tag{C}$$ mainly because of the fundamental lemma of calculus of variations, cf. eq. (D) & (E) below.
Ideologically, OP seems to operate with much more restrictive BCs for the infinitesimal variations $\delta q^k$ than the BCs for the configurations $q^k$ themselves. This is all fine and dandy to recover the EL formula (C) but it is artificial/unnatural from a variational point of view: The BCs on the infinitesimal variations $\delta q^k$ should descend directly from the BCs on the configurations $q^k$, nothing else. We will assume this from now on: There is only one set of BCs in a well-posed variational problem.
As one might expect the caveat is that the functional derivative (FD) only exists for some BCs. If the Lagrangian $L(q,\dot{q};t)$ doesn't depend on higher time-derivatives, there are only two types of BCs that make the functional/variational derivative well-defined, because of the need to eliminate boundary terms in the infinitesimal variation $$\begin{align} \delta I ~=~& \int_{t_i}^{t_f} \! dt\left( \underbrace{EL_k ~\delta q^k}_{\text{bulk}}+\frac{d}{dt}\underbrace{(p_k~\delta q^k)}_{\text{boundary}} \right), \cr p_k~:=~&\frac{\partial L}{\partial \dot{q}^k} ,\end{align}\tag{D}$$ in order to comply with the defining property $$\delta I ~=~ \int_{t_i}^{t_f} \! dt~\frac{\delta I}{\delta q^k} ~\delta q^k \tag{E} $$ of a FD. Comparing eqs. (D) & (E), one indeed finds the possible BCs are
Symmetry of the system may further restrict appropriate BC choices.
It seems that OP basically has the same possible BCs as above in the second half of his post. Note however that a natural BC is not necessarily a
Neumann BC: $\quad \dot{q}(t_i)~=~0\quad\text{and}\quad \dot{q}(t_f)~=~0,$
because momentum $p$ does not need to be proportional to velocity $\dot{q}$.
In the second half of OP's post OP seems to only analyse BCs around on-shell configurations. A similar analysis also applies to off-shell configurations. More importantly: BCs should be imposed to all configurations in order to make the variational problem well-posed, not just on-shell configurations.
Let us mention for completeness that in the second half of OP's post OP is quite close to introduce a (Dirichlet) on-shell action function $S(q_f,t_f;q_i,t_i)$, see e.g. this Phys.SE post.
Example: For a free non-relativistic point particle with Lagrangian $L=\frac{m}{2}\dot{q}^2$, the on-shell action is $$\begin{align} S(q_f,t_f;q_i,t_i)~=~&\frac{m(\Delta q)^2}{2\Delta t},\cr \Delta q~:=~&q_f-q_i,\cr \Delta t~:=~&t_f-t_i. \end{align}\tag{F}$$
Finally, let us repeat that a variational problem typically is ill-defined or un-physical if we don't impose (appropriate) BCs. This is true for both off-shell and on-shell actions, $I[q]$ and $S(q_f,t_f;q_i,t_i)$, respectively. This crucial point seems to essentially invalidate the on-shell analysis of the second half of OP's post. To summarize: it is unnecessary and likely inconsistent to use the on-shell action to deduce the BCs.
Example: If we minimize the on-shell action (F), we would find the periodic BC $q(t_i)=q(t_f)$, which may or may not be physical relevant, depending on context.
--
$^1$ Disclaimer: This answer was originally made as a responds to v8 of the question. From v10 and onwards OP requests to consider the field theoretic (rather than the point mechanical) case. The two cases are quite similar. The main new features in field theory (as compared to point mechanics) are:
Higher spacetime derivatives $\partial_{\mu_1}\ldots \partial_{\mu_r}\phi^{\alpha}$ are totally symmetric in the indices $(\mu_1, \ldots,\mu_r)$, and hence not all independent variables in the Lagrangian density.
The natural BC becomes of the form $n_{\mu}\frac{\partial {\cal L}}{\partial (\partial_{\mu}\phi^{\alpha})}=0$, where $n_{\mu}$ is a normal unit vector to the boundary $\partial \Omega$.
One must chose the same type of BC on a connected component of the boundary $\partial \Omega$.