Say $\Omega$ is a domain in $\mathbb R^3$ with a smooth boundary $\Gamma$.
Consider the spaces $$ H_{n,0}=\{v\in H^1(\Omega):n\cdot v \bigr |_{\Gamma} = 0\} $$ and $$ H_{t,0}=\{v\in H^1(\Omega):n\times v \bigr |_{\Gamma} = 0\}, $$ where $H^1$ is the space of square integrable functions with values in $\mathbb R^3$ and with one spatial derivative and where $n$ denotes the normal vector at the boundary.
Using Green's formula (or integration by parts), for a smooth scalar function $\phi$ $$ \int_\Omega \phi \nabla\cdot v ~d\omega = \int_\Gamma \phi v\cdot n ~d\gamma - \int_\Omega \nabla \phi \cdot v ~d\omega, $$ on $H_{n,0}$, one can establish that the Gradient operator $\nabla$ is dual to the divergence operator $\nabla \cdot$. For simplicity, I leave aside all questions related to the right choice of the spaces.
Question: How does one construct or specify the dual operator to the divergence $\nabla \cdot$ on $H_{t,0}$?
Any intuitive idea and any reference is highly appreciated.
It is a good idea to start with the standard approach presented, e.g., in Ch. 13 of Functional Analysis by Walter Rudin for the case of one Hilbert space. The same approach is widely employed in the case of two Hilbert spaces as well. The key idea of this approach is to consider an unbounded operator $$ L={\rm div}\colon\;D_L=H_{t,0}\subset\bigl(L^2(\Omega)\bigr)^3 \to L^2(\Omega)\tag{1} $$ with its dual unbounded operator $$ L^{\ast}=\nabla\colon\;D_{L^{\ast}}\subset L^2(\Omega) \to\bigl(L^2(\Omega)\bigr)^3,\tag{2} $$ the domain $D_{L^{\ast}}$ of which is defined as $$ D_{L^{\ast}}=\{v\in L^2(\Omega)\colon\; |\Lambda_v(u)| \leqslant C_v\|u\|_{L^2}\;\forall\, u\in D_L \}, $$ where $\Lambda_v$ denotes the linear functional $$ \Lambda_v(u)=\int\limits_{\Omega}v\,{\rm div\,}u\,dx\quad\forall\,u\in D_L\,. $$ Since $\bigl(C_0^{\infty}(\Omega)\bigr)^3\subset H_{t,0}\,$, domain $D_L$ is dense in $\bigl(L^2(\Omega)\bigr)^3$, and hence there is a continuos linear extension $\widetilde{\Lambda_v}$ of $\Lambda_v$ from $D_L$ to the whole $\bigl(L^2(\Omega)\bigr)^3$. By the Riesz representation theorem, there is $w\in\bigl(L^2(\Omega)\bigr)^3$ such that a linear continuous functional $\widetilde{\Lambda_v}$ can be represented in the form $$ \widetilde{\Lambda_v}(u)=\int\limits_{\Omega}w\cdot u\,dx \quad\forall\, u\in \bigl(L^2(\Omega)\bigr)^3. $$ Thus we get an integral identity $$ \int\limits_{\Omega}v\,{\rm div\,}u\,dx=\int\limits_{\Omega}w\cdot u\,dx \quad\forall\, u\in H_{t,0}\; \Rightarrow \; \forall\, u\in\bigl(C_0^{\infty}(\Omega)\bigr)^3\tag{3} $$ with $v\in L^2(\Omega)$ and $w\in\bigl(L^2(\Omega)\bigr)^3$. Taking $w=(\varphi,0,0)$ in $(3)$ we find $$ \int\limits_{\Omega}v\,\partial_{x_1}\varphi\,dx=\int\limits_{\Omega}w_1\varphi\,dx \quad\forall\, \varphi\in C_0^{\infty}(\Omega) $$ whence follows the existence of a weak derivative $\partial_{x_1}v=-w_1$. In a similar manner, $(3)$ implies the existence of weak derivatives $\partial_{x_j}v=-w_j\,,\,j=2,3$. Therefore $\varphi\in H^1(\Omega)$ and $w=-\nabla v$, while $$ \int\limits_{\Omega}v\,{\rm div\,}u\,dx=-\int\limits_{\Omega}\nabla v\cdot u\,dx \quad\forall\, u\in H_{t,0} $$ which immediately implies that $$ \int\limits_{\Omega}v\,u\cdot n\,dx=0\;\, \forall\,u\in H_{t,0} \; \Rightarrow \;\int\limits_{\Omega}v\,\psi\,dx=0\;\, \forall\,\psi\in H^1(\Omega)\; \Rightarrow \;v|_{\partial\Omega}=0 $$ whence follows $v\in H^1_0(\Omega)$, i.e., domain $D_{L^{\ast}}=H^1_0(\Omega)$. Such dual operator $L^{\ast}$ is sometimes referred to as a "Hilbert space adjoint" of $L$. Now it is clear that for the operator $L^{\ast\ast}=(L^{\ast})^{\ast}$, its domain $D_{L^{\ast\ast}}$ does contain a Sobolev space $H^1(\Omega)$. More precisely, by definition, $$ D_{L^{\ast\ast}}=\{u\in (L^2\bigl(\Omega)\bigr)^3:\,|\lambda_u(v)| \leqslant c_u\|v\|_{L^2}\;\;\forall\,v\in D_{L^{\ast}}\} $$ where $D_{L^{\ast}}=H^1_0(\Omega)$ while $\lambda_u$ denotes the linear functional $$ \lambda_u(v)=\int\limits_{\Omega}u\cdot\nabla v\,dx \quad\forall\,v\in D_{L^{\ast}}\,. $$ Again, there is a continuos linear extension $\widetilde{\lambda_u}$ of $\lambda_u$ from $D_{L^{\ast}}$ to the whole $L^2(\Omega)$. By the Riesz representation theorem, there is $f\in L^2(\Omega)$ such that a linear continuous functional $\widetilde{\lambda_u}$ can be represented in the form $$ \widetilde{\lambda_u}(v)=\int\limits_{\Omega}f\,v\,dx \quad\forall\, v\in L^2(\Omega). $$ Thus we get an integral identity $$ \int\limits_{\Omega}u\cdot\nabla v\,dx=\int\limits_{\Omega}f\,v\,dx \quad\forall\, v\in H^1_0(\Omega)\; \Rightarrow \; \forall\, u\in C_0^{\infty}(\Omega)\tag{4} $$ with $u\in\bigl(L^2(\Omega)\bigr)^3$ and $f\in L^2(\Omega)$, whence follows the existence of a weak divergence ${\rm div\,}u=-f$, as they say "in distribution sense". The latter implies that $D_{L^{\ast\ast}}=H({\rm div};\Omega)$ where $$ H({\rm div};\Omega)\overset{\rm def}{=} \{u\in \bigl(L^2(\Omega)\bigr)^3\colon\;{\rm div\,}u\in L^2(\Omega)\} $$ with a weak divergence in the sense $(4)$. The double adjoint $L^{\ast\ast}$ is known to coincide with the closure of $L$, i.e., with the minimal closed extension of $L$. Thus we conclude that the unbounded operator $(3)$ is not closed, while its closure is devoid of any prescribed boundary conditions. And of course, this cannot be the case with the operator $$ L={\rm div}\colon\;D_L=H_{n,0}\subset\bigl(L^2(\Omega)\bigr)^3 \to L^2(\Omega), $$ the closure of which $\widetilde{L}$ has a domain $$ D_{\widetilde{L}}=\{u\in H({\rm div};\Omega)\colon\;(u\cdot n)|_{\partial\Omega}=0\} $$ where for every element $u\in H({\rm div};\Omega)$, its normal component $u\cdot n$ does possess a trace $(u\cdot n)|_{\partial\Omega}\in H^{-1/2}(\partial\Omega)$ — for details see, e.g., p. 129 in An Introduction to Navier-Stokes Equation and Oceanography by Luc Tartar.