I am reading Gazzola et al's book and there for a space $H^{m}(\Omega)$ they seem to allow boundary conditions that are of degree as high as $2m-1$. How is this reconciled with the fact that technically no derivative higher than $m$ has to exist in this space? The general form of the boundary conditions that they prescribe is
\begin{equation}
\sum_{|\alpha|\leq m_j} b_{j,\alpha}(x)D^{\alpha}u = h_j, \quad x\in\partial\Omega, \,\, j=1\ldots m,
\end{equation}
where $m_j\leq 2m-1$ as noted. Then they require that $h_j\in H^{m-m_j-\frac{1}{2}}(\partial\Omega)$. And if $m-m_j-\frac{1}{2}<0$ they identify $H^{m-m_j-\frac{1}{2}}(\partial\Omega)$ to be the dual space of $H^{-m+m_j+\frac{1}{2}}(\partial\Omega)$:
\begin{equation}
H^{m-m_j-\frac{1}{2}}(\partial\Omega) = [H^{-m+m_j+\frac{1}{2}}(\partial\Omega)]' = \textrm{the dual space of }H^{-m+m_j+\frac{1}{2}}(\partial\Omega).
\end{equation}
I am quite confused regarding prescribing functionals on the boundary instead of functions. Even more so when I am not sure how they define $D^{\alpha}$ when $|\alpha|>m$. Do they somehow interpret $D^{\alpha}$ as a weak derivative on $\partial\Omega$ in this case? Is this why we get functionals there? Do I need to somehow bring these boundary conditions to a weak form? How would I do that?
Edit: To clarify my question I have written a simple example problem, and my main issue with it. Ideally I would like to relate it to the theory in Gazzola et al's book.
An Example
Here is an example to illustrate things. I will pick one of the simplest polyharmonic formulations I can think of – one that can be decomposed in $m$ pure-Neumann Poisson problems.
\begin{alignat}{3}
(-\Delta)^m u(\vec{x}) &=0, &\quad& \vec{x}\in \Omega, \\
\partial_{\vec{n}}(-\Delta)^lu(\vec{x}) &= 0, &\quad& \vec{x}\in\partial\Omega, \, 0\leq l \leq m-1.
\end{alignat}
This can be decomposed as
\begin{alignat}{3}
u_m(\vec{x}) &=0, &\quad& \vec{x}\in\Omega, \\
(-\Delta) u_i &= u_{i+1}, &\quad& \vec{x}\in\Omega, \, 0\leq i \leq m-1, \\
\partial_{\vec{n}}u_i(\vec{x}) &= 0, &\quad& \vec{x}\in\partial\Omega, \, 0\leq i \leq m-1.
\end{alignat}
Of course, here the solution is not unique, but it can be made unique through some additional integral constraints I have omitted for simplicity.
Variational Formulation for $u,v\in H^{2m}(\Omega)$
If I assume $u,v\in H^{2m}(\Omega)$ then the following variational formulation corresponds to the BVP:
\begin{equation}
\min_{u\in\mathcal{U}}E(u), \quad E(u) = B(u,u), \, B(v,u) = \int_{\Omega} v\,(-\Delta)^mu,
\end{equation}
where the space $\mathcal{U}$ is a subspace of $H^{2m}(\Omega)$ with the boundary conditions $\partial_{\vec{n}}(-\Delta)^l u(\partial\Omega) = 0$ for $0\leq l \leq m-1$. Taking the Gateaux derivative of the above energy yields:
\begin{equation}
\frac{dE}{dv}(u) = \lim_{\epsilon\to 0}\frac{E(u+\epsilon v)-E(u)}{\epsilon} = B(u,v) + B(v,u).
\end{equation}
In order to find the functional derivative $\frac{\delta E}{\delta u}$ I need to bring $B(u,v)$ to a form where there is no $(-\Delta)^m$ in front of $v$. To do this I apply the divergence theorem $2m$ times
\begin{equation}
B(u,v) = \sum_{l=0}^{m-1}\int_{\partial\Omega}(\partial_{\vec{n}}(-\Delta)^lu)\,(-\Delta)^{m-1-l}v – \sum_{l=0}^{m-1}\int_{\partial\Omega}(-\Delta)^lu\,(\partial_{\vec{n}}(-\Delta)^{m-1-l}v) + B(v,u).
\end{equation}
From $u\in\mathcal{U}$ the boundary terms in the first sum vanish, and from $v\in\mathcal{V}$ the boundary terms in the second sum vanish. Then I am left with
\begin{equation}
\int_{\Omega} v \, \frac{\delta E}{\delta u} = 2B(v,u) = 2\int_{\Omega} v\, (-\Delta)^m u.
\end{equation}
Then it is clear that $\frac{\delta E}{\delta u} = 2(-\Delta)^m u$ and setting this to zero brings me to the $m$-harmonic boundary value problem $(-\Delta)^m u = 0$.
Weak Formulation for $u,v\in H^m(\Omega)$
The above problem assumes that $u\in H^{2m}(\Omega)$ and I want to bring this to a weak form so that is in $u,v\in H^m(\Omega)$. I can now start from the same formulation and bring it to a weak one by applying the divergence theorem $m$ times.
\begin{equation}
\int_{\Omega} v\, (-\Delta)^m u =
\begin{cases}
\int_{\Omega} \Delta^{2k} v \, \Delta^{2k} u,&\text{for } m = 2k, \, j\in\mathbb{N}_+, \\
\int_{\Omega} \nabla \Delta^{2k} v \cdot \nabla\Delta^{2k} u,&\text{for } m = 2k+1, \, j\in\mathbb{N}_0. \\
\end{cases}
\end{equation}
For simplicity let $m=2k$ then the boundary terms that pop out from applying the divergence theorem $m$ times are
\begin{equation}
\int_{\Omega} v\, (-\Delta)^{2k} u = \sum_{l=0}^{k-1}\int_{\partial\Omega} (\partial_{\vec{n}} \Delta^lv)\, \Delta^{m-1-l} u – \sum_{l=0}^{k-1}\int_{\partial\Omega} \Delta^lv\, (\partial_{\vec{n}}\Delta^{m-1-l} u) + \int_{\Omega} \Delta^k u\, \Delta^k v.
\end{equation}
It is clear that here the highest derivative for $v$ is $\partial_{\vec{n}}\Delta^{k-1}v$ which is of degree $2k-1 = m-1$, so it is within $H^{m-\frac{1}{2}}(\partial\Omega)$. However, the derivatives for $u$ are all higher than degree $m$, so what even are those terms supposed to be? It seems like some abuse of notation to me.
Gazzola et al's Book
I was looking for some examples on how to handle that problem and I stumbled upon Gazzola et al's book. It seems to cover things in a somewhat abstract setting, so likely a lot went over my head. I wanted to find some concrete examples on how $D^{\alpha}u$ are treated for $\alpha\geq m$, i.e. how they compute the Green's adjoint boundary operators mentioned in the book (e.g. (2.39)). Below I have shortly described my current understanding of the part of the book relevant to my problem.
What I understood from Gazzola et al's book is that when I have boundary conditions $\{B_1,\ldots, B_p\}$ in the strong formulation with a highest derivative lower than $m$ then I can keep them as usual and incorporate them in the solution space. They term those stable (see right before (2.34)). The remaining conditions $\{B_{p+1},\ldots,B_m\}$ with highest derivative greater or equal than $m$ are termed natural, and seem to require special treatment since they do not fit in $H^m$.
They assume that $\{B_i\}_{i=1}^m$ form a normal system (2.36), and that $m_i+m_j\ne 2m+1$ (2.37) (which is the case for the ones I picked since $2k_1+1 + 2k_2+1 = 2(k_1+k_2)+2$ is even). Then they argue that there exists a system (not unique) of boundary operators of order $m$
$$\{F_j\}_{j=1}^m = \{B_1,\ldots,B_p,S_{p+1},\ldots,S_{m}\},$$
which they reorder in Dirichlet order (i.e. the highest degree of $F_j$ is $j-1$).
Finally they need the Green's adjoint boundary operators $\{\Phi_i\}_{i=1}^m$, which they choose so that the natural conditions are all part of it, i.e. $\{B_{j}\}_{j=p+1}^{m} \subset \{\Phi_i\}_{i=1}^m$. This naturally puts some restrictions on how $\{S_j\}_{j=p+1}^m$ can be chosen.
Looking at Theorem 2.13, if I wanted to relate this to the weak form from my example I would guess that $\{F_j\}_{j=1}^m$ are the derivative operators showing up in front of $v$, and $\{\Phi_j\}_{j=1}^m$ are the derivative operators showing up in front of $u$ in the boundary integrals here:
\begin{equation}
\int_{\Omega} v\, (-\Delta)^{2k} u = \sum_{l=0}^{k-1}\int_{\partial\Omega} (\partial_{\vec{n}} \Delta^lv)\, \Delta^{m-1-l} u – \sum_{l=0}^{k-1}\int_{\partial\Omega} \Delta^lv\, (\partial_{\vec{n}}\Delta^{m-1-l} u) + \int_{\Omega} \Delta^k u\, \Delta^k v.
\end{equation}
Luckily for me, the $\{F_j\}_{j=1}^m$ form a Dirichlet system once they are ordered properly, and $\{\Phi\}_{j=1}^m$ form a normal system once they are also ordered properly. Is this a correct interpretation of my example in the context of the book? Still, how do the operators $\{\Phi\}_{j=1}^m$ make sense considering $u$ is supposedly from $H^m$ and not $H^{2m}$?
What would I have done if $\{F_j\}_{j=1}^m$ didn't form a Dirichlet system after I was done with the divergence theorem? I can't think of a specific example of this, but ideally I would want to see something concrete in order to understand it better.
Edit: Taking into account humanStampedist's answers, I have arrived at the following:
Notation Mismatch
It seems that there is a mismatch in the notation of Gazzola et al's book. In (2.29) they define $\Psi$ using the weak formulation without the boundary terms, and in (2.32) they supposedly have the strong formulation that arises from (2.29) by integrating by parts resulting in the boundary terms being there. Up until here everything is consistent. However, the $\Psi$ in (2.42) is actually (2.29) with the boundary terms. This becomes clear when one considers their Steklov example.
Here is the Steklov example:
\begin{alignat}{3}
(-\Delta)^2 u(x) &=0, &\quad&\text{for } x\in\Omega, \\
u(x) &= 0, &\quad&\text{for } x\in\partial\Omega, \\
\Delta u(x) -a\frac{\partial u}{\partial \nu}(x) &= 0, &\quad&\text{for } x\in\partial\Omega.
\end{alignat}
Originally, according to (2.29) and (2.32) we ought to have
\begin{alignat}{3}
\Psi(v,u) &= \int_{\Omega} \Delta v\, \Delta u, &\quad& u,v \in H^2(\Omega) \\
\Psi(v,u) &= \int_{\Omega} v\, \Delta^2 u + \int_{\partial\Omega}\partial_{\nu}v\, \Delta u – \int_{\partial\Omega}v\, \partial_{\nu}\Delta u , &\quad& u,v \in H^4(\Omega)
\end{alignat}
Instead the Steklov example from (2.42) reads
\begin{equation}
\Psi(v,u) = \int_{\Omega} \Delta v \, \Delta u – \int_{\partial\Omega}a \partial_{\nu}v\, \partial_{\nu} u.
\end{equation}
To arrive at said example we can subtract the boundary terms from (2.29) and (2.32):
\begin{alignat}{3}
\Psi(v,u) &= \int_{\Omega} \Delta v\, \Delta u – \int_{\partial\Omega}\partial_{\nu}v\, \Delta u + \int_{\partial\Omega}v\, \partial_{\nu}\Delta u, &\quad& u,v \in H^2(\Omega) \\
\Psi(v,u) &= \int_{\Omega} v\, \Delta^2 u, &\quad& u,v \in H^4(\Omega)
\end{alignat}
Of course, the boundary terms in the first expression do not exactly make sense since $u\in H^2(\Omega)$. However by substituting the Steklov boundary conditions $v(\partial\Omega) = 0$ and $\Delta u = a\partial_{\nu} u$ in the first expression, we get the desired form (which now also makes sense for $u\in H^2(\Omega)$:
\begin{equation}
\Psi(v,u) = \int_{\Omega} \Delta v \, \Delta u – \int_{\partial\Omega}a \partial_{\nu}v\, \partial_{\nu} u.
\end{equation}
I understand the motivation behind defining (2.32) as they did, since the boundary terms do not make sense unless we assume $u,v \in H^{2m}(\Omega)$, however the usage of $\Psi$ in (2.42) is simply misleading considering its prior definition and the Steklov example that contradicts said definition. To arrive at the $\Psi$ from (2.42), one actually needs to subtract the boundary terms and substitute them for the boundary conditions to get the weak formulation.
General Formulation
I am still unclear about how one would proceed with a general formulation where the boundary conditions are not as nice. We are given the boundary operators $B_j,\, 1\leq j\leq m$, with highest degree $m_j\leq 2m-1$ (i.e. $\partial^{m_j}_{\nu} u$ is the highest derivative of $\partial_{\nu} u$ in $B_j$) and we require $m_i+m_j\ne 2m-1$ for $i\ne j$ (2.37). We split the operators into stable ($m_j<m$) and natural $(m_j\geq m)$, and we consider the solution space (2.34):
\begin{equation}
V = \{ v \in H^m(\Omega)\,:\, (B_jv)(\partial\Omega) = 0 \text{ for } m_j<m\}.
\end{equation}
That is, we keep only the stable boundary operators in the space.
Now for simplicity let us consider $m=2k, \, k\in\mathbb{N}$, then we have:
\begin{equation}
\int_{\Omega} v\, (-\Delta)^{2k} u = \sum_{l=0}^{k-1}\int_{\partial\Omega} (\partial_{\vec{n}} \Delta^lv)\, \Delta^{m-1-l} u – \sum_{l=0}^{k-1}\int_{\partial\Omega} \Delta^lv\, (\partial_{\vec{n}}\Delta^{m-1-l} u) + \int_{\Omega} \Delta^k u\, \Delta^k v.
\end{equation}
We can now take the above boundary operators in front of $v$ and $u$ and put them in two adjoint vectors, where the highest degrees of the components sum up to $m=2k$:
\begin{align}
\begin{bmatrix} F' \\ \Phi' \end{bmatrix} =
\begin{bmatrix}
id & \partial_{\nu} & \Delta & \partial_{\nu}\Delta & \ldots & \Delta^{k-1} & \partial_{\nu}\Delta^{k-1} \\
\partial_{\nu}\Delta^{2k-1} & \Delta^{2k-1} & \partial_{\nu} \Delta^{2k-2} & \Delta^{2k-2} & \ldots & \partial_{\nu} \Delta^{k} & \Delta^{k} \end{bmatrix}
\end{align}
If the $B_j$ are of the above form, then the stable ones are a subset of $F'$ and the natural ones a subset of $\Phi'$ and we are done. The only thing left to do is to substitute the boundary conditions in the boundary terms in order to eliminate all terms of degree $\geq m$.
I am not so clear what to do when the $B_j$ do not have the above form however. Even simple Dirichlet boundary conditions $id, \partial_{\nu}, \ldots, \partial^{m-1}_{\nu}$ do not match the above form of the boundary operators $F'$. My understanding is that at that point one must decompose $\Delta^l$ into a linear combination of powers of $\partial_{\nu}$. The authors have provided the following example in Remark 2.11:
\begin{equation}
\Delta u = \partial^2_{\nu}u + H\partial_{\nu} u + \Delta_{\tau} u \quad\text{on } \partial\Omega.
\end{equation}
So supposedly one would try to decompose the terms in $F'$, $\Phi'$, and $B_j$ as linear combinations of $\partial_{\nu}$ and then try to match those, from where they would get the adjoints. More precisely, if I assume that the following decomposition is valid (I do not know whether this is the case):
\begin{equation}
\Delta^i = \sum_{j=0}^{2m-2} \alpha_{i,j} \partial^j_{\nu} + G_i, \quad 0\leq i \leq m-1,
\end{equation}
where $G_i$ contains tangential terms such as $\Delta_{\tau}$, then I should be able to rewrite the boundary terms as follows:
\begin{align}
\int_{\Omega} v\, (-\Delta)^{2k} u &= \sum_{l=0}^{k-1}\int_{\partial\Omega} (\partial_{\vec{n}} \Delta^lv)\, \Delta^{m-1-l} u – \sum_{l=0}^{k-1}\int_{\partial\Omega} \Delta^lv\, (\partial_{\vec{n}}\Delta^{m-1-l} u) + \int_{\Omega} \Delta^k u\, \Delta^k v \\
&= \sum_{i=0}^{k-1}\int_{\partial\Omega}\left(\sum_{j=0}^{2m-2}\alpha_{i,j}\partial^{j+1}_{\nu} + \partial_{\nu}G_i \right)v \, \left(\sum_{j=0}^{2m-2}\alpha_{m-1-i,j}\partial^{j}_{\nu} + G_{m-1-i}\right)u \\
&- \sum_{i=0}^{k-1}\int_{\partial\Omega}\left(\sum_{j=0}^{2m-2}\alpha_{i,j}\partial^{j}_{\nu} + G_i \right)v \, \left(\sum_{j=0}^{2m-2}\alpha_{m-1-i,j}\partial^{j+1}_{\nu} + \partial_{\nu}G_{m-1-i}\right)u \\
&+ \int_{\Omega} \Delta^k v \, \Delta^k u \\
&= \int_{\Omega} \Delta^k v \, \Delta^k u + \int_{\partial\Omega} \begin{bmatrix} v & \partial_{\nu}v & \ldots & \partial^{2m-1}_{\nu}v\end{bmatrix} C \begin{bmatrix} u \\ \partial_{\nu} u \\ \vdots \\ \partial^{2m-1}_{\nu} u \end{bmatrix},
\end{align}
where C is some matrix of coefficients. Let's call the stable operators $B_j$ and the natural $D_j$. If I can group the terms by $B_j v$ and $D_j u$ I get corresponding $B'_j u$ and $D'_j v$ terms attached to those. Supposedly $F = \{-B_j\} \cup \{B'_j\}$ and $\Phi = \{-D_j\} \cup \{D'_j\}$ from (2.39). It also looks like I cannot pick arbitrary $B_j$ and $D_j$ since the matrix $C$ seems to have a specific structure, which would not be satisfied for arbitrary boundary conditions.
Best Answer
First let us talk about trace theorems (see e.g. https://en.wikipedia.org/wiki/Trace_operator). If $u\in H^1(\Omega)$ and the boundary of $\Omega$ has some regularity, then there exists a linear continuous operator $\gamma:H^1(\Omega)\rightarrow L^2(\partial\Omega)$, which has the property, that for $u:\overline{\Omega}\rightarrow\mathbb{R}$ smooth we have $$\gamma(u)=u|_{\partial\Omega}.$$ On the flip side, you cannot expect for a generic $u\in H^1(\Omega)$ to have well defined values of $\nabla u$ on the boundary. To see this, remember that weak derivates are defined by partial integration, i.e. for all $\varphi\in C^\infty_0(\Omega)$ we have $$\int_\Omega u\nabla \varphi\, dx =-\int_\Omega \nabla u\varphi\, dx.$$ If you now change $u$ or $\nabla u$ on a set of Lebesgue measure zero, the integral does not see that and the equation is still correct. Given some boundary regularity (which is the usual case in applications), we have that $\partial\Omega$ is a set of measure zero. Therefore the derivative of $u$ cannot be well defined on the boundary. In essence the trace operator $\gamma$ chooses a good representative, such that it is the identity if $u$ is regular enough. This is in general not possible for the derivative. In the context of Gazzola, Grunau and Sweers book, boundary conditions, which can be prescribed by such a trace operator are called stable.
The natural boundary conditions are handled differently. Let us discuss an example: $$\left\{\begin{array}{cc}-\Delta u =f,& \mbox{ on }\Omega\\ \partial_n u=0,& \mbox{ on }\partial \Omega.\end{array}\right.$$ Here $\partial_n u$ is the derivative of $u$ w.r.t. to the inner normal of $\partial \Omega$. A weak formulation now only wants $u\in H^1(\Omega)$, i.e. only one derivative. The discussion above on the other hand yields $\nabla u|_{\partial\Omega}$ not to be well defined. Hence we incorporate this notion into our definition of a weak solution. We say $u\in H^1(\Omega)$ is a weak solution, if and only if $$\int_\Omega\nabla u\nabla \varphi\, dx =\int_\Omega f\varphi\, dx\mbox{ for all }\varphi\in H^1(\Omega).$$ Please not that $\varphi$ is not coming from $H^1_0(\Omega)$, i.e. is allowed to have nontrivial values on the boundary, if it were more regular. That is exactly the idea. Instead of extending values of a solution to the boundary, we allow our space of test functions (test functions are these $\varphi$) to be bigger.
Let us now check, that we indeed have a good weak formulation: One prerequisite for that is, if $u$ is a weak solution and is smooth up to the boundary, then it is also a classical solution. Let us check this: First choose $\varphi\in C^\infty_0(\Omega)\subseteq H^1(\Omega)$. Then by partial integration (no boundary terms, because $\varphi$ is zero on the boundary) we get $$0=\int_\Omega\nabla u\nabla\varphi\, dx - \int_\Omega f\varphi\, dx =\int_\Omega ((-\Delta)u - f)\varphi\, dx.$$ Then the fundamental lemma of variational calculus (see e.g. https://en.wikipedia.org/wiki/Fundamental_lemma_of_calculus_of_variations) yields $$-\Delta u=f\mbox{ on }\Omega.$$ Since we are allowed to choose $\varphi\in C^\infty(\overline{\Omega})\subset H^1(\Omega)$, we get again by partial integration (this time with boundary terms) $$0=\int_\Omega\nabla u\nabla \varphi - f\varphi\, dx =\int_\Omega((-\Delta)u - f)\varphi\, dx + \int_{\partial\Omega}\partial_nu \varphi\, dA = \int_{\partial\Omega}\partial_nu \varphi\, dA.$$ The integral over $\Omega$ vanishes, because the integrand is zero by our calulations above. Again the fundamental lemma of calculus of variations yields $$\partial_n u=0\mbox{ on }\partial\Omega,$$ i.e. exactly our boundary condition.
To reiterate: Stable boundary conditions can be handled by trace theorems and natural boundary conditions by picking a larger class of test functions.
Now in Gazzola, Grunau and Sweers this is formulated for higher order equations. In (2.42) such a weak formulation is given. After Thm 2.15 they give a similar argument to the one I have given you here, to show, that if a weak solution is regular enough, it is indeed also a classical solution.
Now the question what happens if these boundary conditions do not form a Dirichlet system, i.e. have some natural parts in it. In essence you need some so called compatibility conditions. In our example here it would be $\int_\Omega f\, dx =0$. To see this, integrate the differential equation and use divergence theorem: $$0=\int_\Omega(-\Delta) u -f \, dx =\int_{\partial\Omega}u\partial_n u\, dA - \int_\Omega f\, dx =-\int_\Omega f\, dx.$$ Therefore a solution can only exist, if $\int f\, dx =0$. If the boundary conditions are more complicated, these compatibility conditions become more involved. This essentially results in Gazzola, Grunau, Sweers in these Green's adjoint operators.