How to impose boundary condition for 2-dimensional curl-curl weak problem

First about the appropriate boundary condition we would like to impose: The cross product(scalar-$\mathrm{curl}$) in $\mathbb{R}^2$ is done by we embed a vector into $\mathbb{R}^3$: $\boldsymbol{v} = \langle v_1,v_2,0\rangle$, then apply the cross product, what we get would be

\[\nabla \times \boldsymbol{v} = \left\langle 0,0,\left|\begin{matrix}\partial_x&\partial_y\\ v_1&v_2 \end{matrix}\right|\right\rangle\]

Similarly the vector-$\mathbf{curl}$ is we embed a scalar function into the third dimension of $\mathbb{R}^3$: $u = \langle0,0,u\rangle$, and what we get is the rotation of the $\nabla$-operator normally denoted by $\nabla^{\perp}$:

\[\nabla^{\perp} u = \left\langle\frac{\partial u}{\partial y},-\frac{\partial u}{\partial x},0 \right\rangle\]

Reference of how to define 2D-curl operators: Finite element methods for Navier-Stokes equations: theory and algorithms.

Then in $\mathbb{R}^2$ our problem can be written as:

\[\nabla^{\perp}(\mu^{-1}\nabla \times \boldsymbol{A}) = \mu \boldsymbol{J} \tag{1}\]

Suppose we multiply a test vector $\boldsymbol{v}$ vanishing on boundary on both side or the equation and integration by parts:

\[\displaystyle\int_{\Omega}\mu^{-1} \nabla\times \boldsymbol{A} \cdot \nabla \times \boldsymbol{v} - \int_{\partial \Omega} \mu^{-1}(\nabla \times \boldsymbol{A})(\boldsymbol{v} \times \boldsymbol{n}) = \int_{\Omega} \mu\boldsymbol{J}\cdot \boldsymbol{v}\]

where recall in $\mathbb{R}^2$, $B = \mu^{-1}(\nabla \times \boldsymbol{A})$ is really a scalar instead of vector, hence for the tangential Neumann boundary condition, doing the $\mathbb{R}^3$-embedding we have:

\[\boldsymbol{n}\times B = \langle n_2 B, -n_1 B,0 \rangle\]

this is the value of the scalar $B$ along the tangential direction of the domain $\boldsymbol{t} = \langle n_2, -n_1,0 \rangle$, dot product with $\boldsymbol{t}$ itself we could get $B$($\boldsymbol{t}$ is the rotation of a unit outward normal vector), hence the Neumann type boundary condition you really wanna impose here for 2D is:

\[\mu^{-1}(\nabla \times \boldsymbol{A}) = g\]

where $g$ is a scalar function, and our variational form becomes:

\[\int_{\Omega}\mu^{-1} \nabla\times \boldsymbol{A} \cdot \nabla \times \boldsymbol{v} = \int_{\Omega} \mu\boldsymbol{J}\cdot \boldsymbol{v} + \int_{\partial \Omega} g(\boldsymbol{v} \cdot \boldsymbol{t}) \tag{2}\]

The finite element method is we triangulating the domain, construct a finite dimensional Hilbert space $\boldsymbol{V}_h$ (usually a subspace of the solution space $\boldsymbol{V}$ to above variational problem), the boundary integral is treated as:

\[\int_{\partial \Omega} g(\boldsymbol{v} \cdot \boldsymbol{t}) = \sum_{e\subset \partial \Omega} \int_{e} g (\boldsymbol{v}\cdot \boldsymbol{t}_e)\]

where $e$ is a boundary edge, and here is the place quadrature kicks in, that you wanna evaluate a line integral on each edge, it depends on what our finite element test function space $\boldsymbol{V}_h$, nowadays a consensus is to choose Nédélec elements, in lowest order case it is $\boldsymbol{v}_e = \lambda_A\nabla \lambda_B - \lambda_B\nabla \lambda_A$ for $e = \overrightarrow{AB}$.

---

Secondly, what is really problematic with our statement here is: for a smooth vector field $\boldsymbol{A} \in \boldsymbol{C}^{\infty}(\mathbb{R}^2)$ or $\boldsymbol{C}^{\infty}_c(\Omega)$, the equation (1) doesn’t simply reduce to the scalar Laplace equation as you mentioned in OP, a direct observation would be that homogeneous Neumann problem for $-\Delta u = f$ has a constant as null space because a solution plus a constant is still a solution, while $\boldsymbol{curl}$-problem has all the gradient vector fields as null space!

Similar thing applies for variational formulation (2), the bilinear form is not coercive, simply plugging a gradient field $\boldsymbol{A} = \boldsymbol{v} = \nabla \varphi$ into it. What we numerical analyst usually do here is to find a solution in the quotient space $\boldsymbol{H}(\mathrm{curl})/\nabla H^1$ by adding a Lagrange multiplier to let the problem fall into mixed finite element(google MFEM if you would like to know more) framework:

\[\displaystyle\int_{\Omega}\mu^{-1} \nabla\times \boldsymbol{A} \cdot \nabla \times \boldsymbol{v} + \int_{\Omega} \boldsymbol{v}\cdot \nabla \varphi= \int_{\Omega} \mu\boldsymbol{J}\cdot \boldsymbol{v} + \int_{\partial \Omega} g(\boldsymbol{v} \cdot \boldsymbol{t}) \text{ for } \forall \boldsymbol{v}\in \boldsymbol{H}(\mathrm{curl}) \\[4pt] \displaystyle\int_{\Omega} \boldsymbol{A} \cdot \nabla \psi = 0 \text{ for } \forall \psi\in H^1/\mathbb{R}\]

to make sure problem (2) is solvable using finite element, not only we have to use Lagrange multiplier to circumvent the “null space” problem, we also have to guarantee that the right side is in the range of the curl-operator, here plug $\boldsymbol{v} = \nabla \varphi$ into (2):

\[\displaystyle0=\int_{\Omega} \mu\boldsymbol{J}\cdot \nabla \varphi + \int_{\partial \Omega} g(\nabla \varphi \cdot \boldsymbol{t}) = -\int_{\Omega} \nabla \cdot(\mu\boldsymbol{J}) \varphi + \int_{\partial \Omega} (\mu\boldsymbol{J} \cdot \boldsymbol{n}) \varphi + \int_{\partial \Omega} g(\nabla \varphi \cdot \boldsymbol{t})\]

in our case, $g=0$, also by the arbitrariness of $\varphi$, for eg we choose $\nabla \varphi \cdot \boldsymbol{t} = 0$ for $\varphi\in H^1_0$, we would get so called “compatibility condition” in finite element method for Neumann problem:

\[\nabla \cdot(\mu\boldsymbol{J}) = 0\text{ in }\Omega\quad \text{ and }\quad \mu\boldsymbol{J}\cdot \boldsymbol{n} = 0\text{ on }\partial\Omega.\]