1 Introduction

The numerical simulation of quantum mechanical systems under the influence of electromagnetic fields plays a pivotal role in modern physics, with applications ranging from condensed matter to quantum computing. The evolution and stationary states of a non-relativistic particle in a magnetic field are governed by the Schrödinger equation, where the magnetic effect is introduced via a vector potential ${\boldsymbol{A}}$ . A fundamental property of this formulation is its gauge invariance: transforming the vector potential as ${\boldsymbol{A}}\to{\boldsymbol{A}}+\nabla\chi$ alongside a local phase shift of the wavefunction $\psi\to\psi e^{i\chi}$ leaves the physical observables unchanged.

From a numerical perspective, preserving this continuous symmetry at the discrete level is a non-trivial task. Standard discretization schemes often break gauge invariance, introducing artificial dependencies on the chosen gauge and potentially leading to spurious eigenvalues or unphysical long-time dynamics. The severe consequences of such discretizations were explicitly highlighted by Governale and Ungarelli [Governale1998], who pioneered the use of Wilson’s lattice gauge theory on uniform finite difference grids to restore gauge invariance. This lattice gauge formalism has since been generalized to broader grid-based and pseudospectral methods by Halvorsen and Kvaal [Halvorsen2009], and recently adapted for the non-perturbative solution of diatomic molecules in strong magnetic fields by Yenugu et al. [Yenugu2025]. For time-dependent problems, other grid-based approaches have also been explored, such as the Finite Difference Time Domain (FDTD) method for complex wavefunctions [Sudiarta2007], splitting methods utilizing rotated potentials to accommodate magnetic rotations [Gradinaru2020], and open-source split-step Fourier solvers [QMsolve].

In the realm of finite element methods (FEM), Alouges and Bonnaillie-Noël [Alouges2006] proposed bypassing standard gauge issues by decoupling the modulus and phase of the wavefunction to study eigenstates in domains with corners. From a preconditioning perspective, Ovall and Zhu [Ovall2025] recently demonstrated that computing a canonical gauge via a Helmholtz-Hodge decomposition can significantly reduce eigenvector oscillations, thereby improving standard FEM efficiency. A major breakthrough in structure-preserving FEM was achieved by Christiansen and Halvorsen [Christiansen2011, Christiansen2015], who successfully merged lattice gauge concepts with finite element exterior calculus to yield gauge-invariant discretizations on simplicial grids, achieving up to second-order convergence. Furthermore, in the context of Discontinuous Galerkin (DG) methods, Yi, Huang, and Liu [Yi2020] developed conservative schemes that preserve mass and energy for nonlinear magnetic Schrödinger equations while achieving optimal $L^{2}$ error estimates. Despite these notable advancements across various discretization paradigms, extending exact discrete gauge invariance to arbitrarily high-order methods on complex unstructured polyhedral meshes, using a purely local formulation remains, to the best of the author’s knowledge, an open problem.

In recent years, the Hybrid High-Order (HHO) method [DiPietro2014] has emerged as a highly flexible framework for the discretization of partial differential equations. By employing fully discontinuous polynomial spaces on mesh elements and faces, HHO methods seamlessly handle arbitrary polyhedral grids and provide dimension-independent constructions.

The primary contribution of this paper is the design and analysis of a gauge-invariant HHO method for the Schrödinger equation. We achieve this by introducing a discrete covariant gradient that intrinsically respects the coupling between the spatial derivative and the magnetic potential. Our approach offers the following main features: (i) Local construction:The discrete covariant gradient is computed element-wise, on a general $d$ -dimensional element, with a controlled accuracy. (ii) Asymptotically discrete gauge covariance:We prove that under a discrete equivalent of the gauge transformation, the local reconstructions and the global bilinear form behave covariantly, ensuring that the physical observables are invariant asymptotically. (iii) Discrete Gårding inequality:Due to the magnetic potential, the standard coercivity of the operator is lost. We bypass this issue by proving a gauge-invariant discrete Gårding inequality, which guarantees that the discrete system possesses a stable ground state and a spectrum bounded from below. The rest of the paper is organized as follows. In Section 2, we recall the continuous problem, establish the continuous Gårding inequality, and review the principles of gauge invariance. Section 3 introduces the HHO discretization, details the construction of the discrete covariant gradient, and presents the main theoretical proofs regarding gauge covariance and stability. Finally, Section 4 is devoted to numerical experiments on 2D and 3D unstructured meshes, validating the theoretical optimal convergence rates and the exactness of the discrete gauge invariance on the computation of the Fock-Darwin spectrum and the long-time behavior of the system.

2 The Covariant Derivative and the Schrödinger Equations

Let ${\boldsymbol{A}}$ denote a vector potential. The central focus of this paper is the covariant derivative operator $(i\boldsymbol{\nabla}-{\boldsymbol{A}})$ . This operator naturally emerges in quantum mechanics when considering charged particles in a magnetic field $\boldsymbol{B}$ , where $B=\boldsymbol{\nabla}{\times}{\boldsymbol{A}}$ , as dictated by the Maxwell-Thomson law. We shall explore two pertinent problems involving this operator: the stationary and unsteady Schrödinger equations. For the sake of simplicity, we set $\hbar=m=q=1$ (the physical constants can be reintroduced without difficulty).

Let $\Omega\subset\mathbb{R}^{d}$ (where $d=2,3$ ) be a Lipschitz bounded domain, and consider a quantum particle $P$ characterized by its wave function $\psi:[0,T]\times\Omega\to\mathbb{C}$ . The probability of finding the particle within an infinitesimal box $dxdy$ at time $t$ is given by the quantity $|\psi(t,x,y)|^{2}dxdy$ , such that $\int_{\Omega}|\psi(t,x,y)|^{2}dxdy=1$ for any $t\in[0,T)$ .

The stationary Schrödinger Equation

Let ${\boldsymbol{A}}\in W^{1,\infty}(\Omega;\mathbb{R}^{d})$ be a vector potential and $V\in L^{\infty}(\Omega;\mathbb{R})$ a scalar potential defined over $\Omega$ . The stationary Schrödinger equation is expressed as an eigenvalue problem: find eigenpairs $\psi$ and $\lambda$ , s.t.

(1)

\begin{cases}\bigl(-i\boldsymbol{\nabla}-{\boldsymbol{A}}\bigr)^{2}\psi+V\psi=\lambda\psi&\text{in }\Omega,\\ \psi=0&\text{on }\partial\Omega.\end{cases}

Here, $\bigl(-i\boldsymbol{\nabla}-{\boldsymbol{A}}\bigr)^{2}$ is interpreted in the sense of $\mathcal{G}_{{\boldsymbol{A}}}^{*}\mathcal{G}_{{\boldsymbol{A}}}$ , where $\mathcal{G}_{{\boldsymbol{A}}}$ and $\mathcal{G}_{{\boldsymbol{A}}}^{*}$ are defined as

\mathcal{G}_{{\boldsymbol{A}}}u:=(-i\boldsymbol{\nabla}u-{\boldsymbol{A}}u\bigr),\qquad\mathcal{G}_{{\boldsymbol{A}}}^{*}\boldsymbol{\tau}=-i\boldsymbol{\nabla}{\cdot}\boldsymbol{\tau}-{\boldsymbol{A}}\cdot\boldsymbol{\tau},

for any $u\in H^{1}_{0}(\Omega;\mathbb{C})$ and $\boldsymbol{\tau}$ in $H^{1}_{0}(\Omega;\mathbb{C})^{d}$ . The eigenvalues are guaranteed to be real, as $\mathcal{G}_{{\boldsymbol{A}}}^{*}\mathcal{G}_{{\boldsymbol{A}}}$ is self-adjoint; however, they may be nonpositive, with the lowest eigenvalue $\lambda_{0}$ being particularly significant in applications. The variational formulation entails identifying $\lambda\in\mathbb{R}$ and $\psi\in H_{0}^{1}(\Omega;\mathbb{C})$ such that

(2)

a(\psi,\phi;{\boldsymbol{A}})=\lambda(\psi,\phi)_{L^{2}(\Omega)}\quad\forall\phi\in H_{0}^{1}(\Omega;\mathbb{C}),\qquad\|\psi\|_{L^{2}(\Omega)}=1,

where the sesquilinear form $a(\cdot,\cdot;{\boldsymbol{A}})$ is defined as

(3)

a(\psi,\phi;{\boldsymbol{A}})\mathrel{\mathop{:}}=\int_{\Omega}\mathcal{G}_{{\boldsymbol{A}}}u\cdot\overline{\mathcal{G}_{{\boldsymbol{A}}}v}\,dx+\int_{\Omega}Vu\overline{v}\,dx.

We will simply write $a(\psi,\phi)$ to refer to $a(\psi,\phi;{\boldsymbol{A}})$ when the context is unambiguous.

The unsteady Schrödinger Equation

Let $\psi_{0}$ represent an initial quantum state and $f\in L^{2}(\Omega;\mathbb{C})$ denote a forcing term. The unsteady (or time-dependent) Schrödinger equation is formulated as follows: find $\psi$ s.t.

(4)

\begin{cases}i\partial_{t}\psi+\bigl(-i\boldsymbol{\nabla}-{\boldsymbol{A}}\bigr)^{2}\psi+V\psi=f&\text{in }\Omega,\\ \psi(0)=\psi_{0}&\text{on }\partial\Omega.\end{cases}

2.1 Gårding Inequality and Spectrum

Before discussing the gauge invariance, we establish a fundamental property of the sesquilinear form (3). Although the form $a(\cdot,\cdot)$ is not strictly coercive due to the presence of the vector potential $\mathbf{A}$ , it satisfies a Gårding inequality [Garding1953]. This property ensures that the quantum system has a stable ground state energy.

Lemma 2.1 (Gårding Inequality).

There exist constants $\alpha>0$ and $\gamma\geq 0$ , depending on $\|\mathbf{A}\|_{L^{\infty}(\Omega)}$ and $\|V\|_{L^{\infty}(\Omega)}$ , such that for all $u\in H_{0}^{1}(\Omega;\mathbb{C})$ ,

a(u,u)\geq\alpha\|\nabla u\|_{L^{2}(\Omega)}^{2}-\gamma\|u\|_{L^{2}(\Omega)}^{2}.

Proof 2.2.

By definition of the sesquilinear form, and taking the real part (which is equal to the form itself since it is Hermitian), we have:

a(u,u)=\int_{\Omega}|(-i\nabla-\mathbf{A})u|^{2}dx+\int_{\Omega}V|u|^{2}dx.

Using the elementary inequality $|a-b|^{2}\geq\tfrac{1}{2}|a|^{2}-|b|^{2}$ for vectors in $\mathbb{C}^{d}$ , we can lower bound the kinetic energy term:

|(-i\nabla-\mathbf{A})u|^{2}\geq\frac{1}{2}|\nabla u|^{2}-|\mathbf{A}u|^{2}.

Integrating this inequality over $\Omega$ and bounding the potential terms by their $L^{\infty}$ norms yields:

a(u,u)\geq\tfrac{1}{2}\|\nabla u\|_{L^{2}(\Omega)}^{2}-\|\mathbf{A}\|_{L^{\infty}(\Omega)}^{2}\|u\|_{L^{2}(\Omega)}^{2}-\|V\|_{L^{\infty}(\Omega)}\|u\|_{L^{2}(\Omega)}^{2}.

The result follows by choosing $\alpha=\frac{1}{2}$ and $\gamma=\|\mathbf{A}\|_{L^{\infty}(\Omega)}^{2}+\|V\|_{L^{\infty}(\Omega)}$ .

Remark 2.3 (Boundedness of the spectrum).

This inequality directly implies that the spectrum of the continuous problem is bounded from below. Indeed, if $(\lambda,\psi)$ is a solution to the eigenvalue problem with $\|\psi\|_{L^{2}(\Omega)}=1$ , then $\lambda=a(\psi,\psi)$ . Applying the Gårding inequality, we obtain:

\lambda\geq\frac{1}{2}\|\nabla\psi\|_{L^{2}(\Omega)}^{2}-\gamma\|\psi\|_{L^{2}(\Omega)}^{2}\geq-\gamma.

The minimal energy of the system is therefore bounded by below by $-\gamma$ . The discrete covariant gradient operator presented in Section 3.4 is specifically designed to preserve this property.

2.2 Gauge Invariance

Definition 2.4.

Let $\chi\in W^{2,\infty}(\Omega;\mathbb{R})$ be an arbitrary gauge function. We define the gauge transformation for any smooth function $\psi$ and vector field ${\boldsymbol{A}}$ as

(5)

\psi_{\chi}:=\psi\,e^{i\chi},\qquad{\boldsymbol{A}_{\chi}}:={\boldsymbol{A}}+\boldsymbol{\nabla}\chi.

Lemma 2.5 (Continuous Gauge Covariance).

We have the following compatibilities with respect to gauge transformation,

\mathcal{G}_{{{\boldsymbol{A}_{\chi}}}}\psi_{\chi}=e^{i\chi}\mathcal{G}_{{\boldsymbol{A}}}\psi,\quad e^{i\chi}\mathcal{G}_{{\boldsymbol{A}}}^{*}(e^{-i\chi}\boldsymbol{\tau})=\mathcal{G}_{{{\boldsymbol{A}_{\chi}}}}^{*}\boldsymbol{\tau},

holds for any smooth functions $\psi$ and $\boldsymbol{\tau}$ , almost everywhere in $\Omega$ . Consequently, the bilinear form $a(\cdot,\cdot)$ , the probability density $|\psi|^{2}$ , and the spectrum remain invariant under this transformation.

Proof 2.6.

A direct computation gives $-i\boldsymbol{\nabla}(\psi e^{i\chi})=e^{i\chi}(-i\boldsymbol{\nabla}\psi+\psi\boldsymbol{\nabla}\chi)$ . Thus,

(-i\boldsymbol{\nabla}-{\boldsymbol{A}}-\boldsymbol{\nabla}\chi)(\psi e^{i\chi})=e^{i\chi}(-i\boldsymbol{\nabla}-{\boldsymbol{A}})\psi\Longleftrightarrow\mathcal{G}_{{{\boldsymbol{A}_{\chi}}}}\psi_{\chi}=e^{i\chi}\mathcal{G}_{{\boldsymbol{A}}}\psi.

The adjoint covariant operator satisfies the exact phase-shifting identity:

	$\displaystyle e^{i\chi}\mathcal{G}_{{\boldsymbol{A}}}^{*}(e^{-i\chi}\boldsymbol{\tau})$	$\displaystyle=e^{i\chi}\Big(-i\boldsymbol{\nabla}\cdot(e^{-i\chi}\boldsymbol{\tau})-{\boldsymbol{A}}\cdot(e^{-i\chi}\boldsymbol{\tau})\Big)$
		$\displaystyle=e^{i\chi}\Big(-\boldsymbol{\nabla}\chi\cdot(e^{-i\chi}\boldsymbol{\tau})-ie^{-i\chi}\boldsymbol{\nabla}\cdot\boldsymbol{\tau}-{\boldsymbol{A}}\cdot(e^{-i\chi}\boldsymbol{\tau})\Big)$
		$\displaystyle=-i\boldsymbol{\nabla}\cdot\boldsymbol{\tau}-({\boldsymbol{A}}+\boldsymbol{\nabla}\chi)\cdot\boldsymbol{\tau}=\mathcal{G}_{{{\boldsymbol{A}_{\chi}}}}^{*}\boldsymbol{\tau}.$

From the first equality, we can readily infer $|\psi_{\chi}|^{2}=|\psi|^{2}$ and

	$\displaystyle a(\psi_{\chi},\phi_{\chi};{\boldsymbol{A}_{\chi}})$	$\displaystyle=\int_{\Omega}\mathcal{G}_{{{\boldsymbol{A}_{\chi}}}}\psi_{\chi}\cdot\overline{\mathcal{G}_{{{\boldsymbol{A}_{\chi}}}}\phi_{\chi}}dx+\int_{\Omega}V\psi_{\chi}\overline{\phi_{\chi}}dx$
		$\displaystyle=\int_{\Omega}e^{i\chi}\mathcal{G}_{{\boldsymbol{A}}}\psi\cdot e^{-i\chi}\overline{\mathcal{G}_{{\boldsymbol{A}}}\phi}dx+\int_{\Omega}Ve^{i\chi}\psi e^{-i\chi}\overline{\phi}dx=a(\psi,\phi;{\boldsymbol{A}}).$

3 The Hybrid High-Order Discretization

In the following discussion, we give a succinct overview of the Hybrid High-Order (HHO) method and expand its discrete gradient to include the discrete counterpart of the covariant gradient $\mathcal{G}_{{\boldsymbol{A}}}$ . To ensure consistency with the notation commonly used in the HHO literature, we will denote the solution to the primary problem as $u$ , reserving the symbol $\psi$ for contexts where the focus is predominantly on the physical aspects.

3.1 Admissible meshes and local polynomial spaces

Denote by ${\mathcal{H}}\subset\mathbb{R}_{*}^{+}$ a countable set of meshsizes having $0$ as its unique accumulation point. We consider $h$ -refined spatial mesh sequences $(\mathcal{T}_{h})_{h\in\mathcal{H}}$ where, for all $h\in\mathcal{H}$ , $\mathcal{T}_{h}=\{T\}$ is a polytopal tessellation of $\Omega$ such that $h=\max_{T\in\mathcal{T}_{h}}h_{T}$ with $h_{T}$ standing for the diameter of the tile $T$ . We assume that mesh regularity holds in the sense of [DiPietro2020, Definition 1.9]; cf. [DiPietro2020, Chapter 1] for a collection of useful geometric and functional inequalities that hold under this assumption.

We define a mesh face $F$ as a hyperplanar closed connected subset of $\bar{\Omega}$ with positive $(d{-}1)$ -dimensional Hausdorff measure and such that (i) either there exist $T_{1},T_{2}\in\mathcal{T}_{h}$ such that $F\subset\partial T_{1}\cap\partial T_{2}$ and $F$ is called an interface or (ii) there exists $T\in\mathcal{T}_{h}$ such that $F\subset\partial T\cap\partial\Omega$ and $F$ is called a boundary face. Interfaces are collected in the set $\mathcal{F}_{h}^{{\rm i}}$ , boundary faces in $\mathcal{F}_{h}^{{\rm b}}$ , and we let $\mathcal{F}_{h}\mathrel{\mathop{:}}=\mathcal{F}_{h}^{{\rm i}}\cup\mathcal{F}_{h}^{{\rm b}}$ . For all $T\in\mathcal{T}_{h}$ , $\mathcal{F}_{T}\mathrel{\mathop{:}}=\{F\in\mathcal{F}_{h}\;|\;F\subset\partial T\}$ denotes the set of faces contained in $\partial T$ and, for all $F\in\mathcal{F}_{T}$ , $\boldsymbol{n}_{TF}$ is the unit normal to $F$ pointing out of $T$ .

We assume throughout the rest of this work that the mesh sequence $(\mathcal{T}_{h})_{\mathcal{H}}$ is admissible in the sense of [DiPietro2020, Chapter 1].

Definition 3.1 (Admissible mesh sequence).

For all $h\in\mathcal{H}$ , $\mathcal{T}_{h}$ admits a matching simplicial submesh $\mathfrak{T}_{h}$ and the following properties hold for all $h\in\mathcal{H}$ with mesh regularity parameter $\varrho>0$ independent of $h$ : (i) for all simplex $S\in\mathfrak{T}_{h}$ of diameter $h_{S}$ and inradius $r_{S}$ , $\varrho h_{S}\leq r_{S}$ ; (ii) for all $T\in\mathcal{T}_{h}$ , and all $S\in\mathfrak{T}_{T}\mathrel{\mathop{:}}=\{S\in\mathfrak{T}_{h}\,|\,S\subset T\}$ , $\varrho h_{T}\leq h_{S}$ .

For an admissible mesh sequence, it is known from [DiPietro2020, Lemma 1.12] that the number of faces of one tile can be bounded uniformly in $h$ , i.e., it holds that

(6)

\forall h\in\mathcal{H},\qquad\max_{T\in\mathcal{T}_{h}}\operatorname{card}(\mathcal{F}_{T})\leq N_{\partial},

for an integer $(d+1)\leq N_{\partial}<+\infty$ depending on $\varrho$ but independent of $h$ . Furthermore, for all $h\in\mathcal{H}$ , all $T\in\mathcal{T}_{h}$ and all $F\in\mathcal{F}_{T}$ , $h_{F}$ is uniformly comparable to $h_{T}$ in the following sense:

\rho^{2}h_{T}\leq h_{F}\leq h_{T}.

Let $X$ be a subset of $\mathbb{R}^{N}$ , $N\geq 1$ , $H_{X}$ the affine space spanned by $X$ , $d_{X}$ its dimension, and assume that $X$ has a non-empty interior in $H_{X}$ (in what follows, the cases $X=T$ and $X=F$ are relevant). For a given integer $k\geq 0$ , we denote by $\mathbb{P}^{k}(X;\mathbb{C})$ the space spanned by $d_{X}$ -variate complex-valued polynomials on $H_{X}$ of total degree $\leq l$ and by $\pi_{X}^{k}:L^{1}(X)\to\mathbb{P}^{k}(X;\mathbb{C})$ the $L^{2}$ -orthogonal projector on $\mathbb{P}^{k}(X;\mathbb{C})$ . We define the global interpolation operator $I_{h}^{k}:H^{1}(\Omega)\to U_{h}^{k}$ such that its local restriction is $I_{T}^{k}u=(\pi_{T}^{k}u,\pi_{\partial T}^{k}u)$ , where $\pi_{\partial T}^{k}$ is the piecewise $L^{2}$ -orthogonal projector on the faces of $T$ . The usual $L^{2}(X)$ -product and norm are denoted by $(u,v)_{X}:=\int_{X}u\overline{v}dx$ and $\|{\cdot}\|_{X}$ , respectively, and we omit the index when $X=\Omega$ . We also recall the following local trace and inverse inequalities (cf. [DiPietro2020, Section 1.2.5]): For all $T\in\mathcal{T}_{h}$ and all $v\in\mathbb{P}^{k}(T)$ ,

(7)

\|v\|_{F}\lesssim h_{F}^{-\frac{1}{2}}\|v\|_{T}\mbox{ for all $F\in\mathcal{F}_{T}$ and }\|\boldsymbol{\nabla}v\|_{T}\lesssim h_{T}^{-1}\|v\|_{T}.

3.2 The degrees of freedom spaces

We follow the standard Hybrid High-Order formulation for the Poisson problem. For a polynomial degree $k\geq 0$ , the local space of degrees of freedom (DOF) on each cell $T$ is

U_{T}^{k}:=\mathbb{P}^{k}(T;\mathbb{C})\times\prod_{F\subset\partial T}\mathbb{P}^{k}(F;\mathbb{C}),

The global DOF space is denoted $U_{h}^{k}$ , with homogeneous Dirichlet boundary conditions stongly enforced:

U_{h}^{k}:=\left\{\underline{u}_{h}=(\underline{u}_{TF})_{T\in\mathcal{T}_{h}},\quad\underline{u}_{TF}\in U_{T}^{k},\quad u_{F}=0,\forall F\in\mathcal{F}_{h}^{{\rm b}}\right\}.

We define the following norm $\|\cdot\|_{1,h}$ over $U_{h}^{k}$ as as

(8)

\|\underline{u}_{h}\|_{1,h}^{2}:=\sum_{T\in\mathcal{T}_{h}}\|\underline{u}_{TF}\|_{1,T}^{2},\qquad\|\underline{u}_{TF}\|_{1,T}^{2}:=\|\boldsymbol{\nabla}u_{T}\|^{2}+\sum_{F\in\mathcal{F}_{T}}h_{F}^{-1}\|u_{F}-u_{T}\|^{2}.

3.3 Potential Reconstruction

For $\underline{u}_{TF}\in U_{T}^{k}$ , the reconstructed potential $p^{k+1}_{T}\underline{u}_{TF}\in\mathbb{P}^{k+1}(T;\mathbb{C})$ is defined such that for all $w\in\mathbb{P}^{k+1}(T;\mathbb{C})$ ,

(9)

(\boldsymbol{\nabla}p^{k+1}_{T}\underline{u}_{TF},\boldsymbol{\nabla}w)_{T}=-(u_{T},\Delta w)_{T}+\sum_{F\subset\partial T}(u_{F},\boldsymbol{\nabla}w\cdot\boldsymbol{n}_{TF})_{F},

or equivalently

(10)

(\boldsymbol{\nabla}p^{k+1}_{T}\underline{u}_{TF},\boldsymbol{\nabla}w)_{T}=(\boldsymbol{\nabla}u_{T},\boldsymbol{\nabla}w)_{T}+\sum_{F\subset\partial T}(u_{F}-u_{T},\boldsymbol{\nabla}w\cdot\boldsymbol{n}_{TF})_{F},

closed with the normalization $(p^{k+1}_{T}\underline{u}_{TF}-u_{T},1)_{T}=0$ . We extend this operator globally to $U_{h}^{k}$ by applying it element-wise, namely we define $p^{k+1}_{h}:U_{h}^{k}\to\mathbb{P}^{k+1}(\Omega;\mathbb{C})$ such that for any $T\in\mathcal{T}_{h}$ ,

(p^{k+1}_{h}\underline{u}_{h})_{|T}:=p^{k+1}_{T}\underline{u}_{TF}.

Proposition 3.2.

There exists $C>0$ independent from $h$ , but from the number of faces and the trace constant $C_{tr}$ such that

(11)

\|\boldsymbol{\nabla}p^{k+1}_{T}\underline{u}_{TF}\|^{2}+s_{T}(\underline{u}_{TF},\underline{u}_{TF})\geq C\|\boldsymbol{\nabla}u_{T}\|^{2}.

Proof 3.3.

If $u_{T}$ is constant over $T$ , the proof is valid as the LHS is always positive. We now assume $u_{T}$ is not constant over $T$ . Using (10), we have

\|\boldsymbol{\nabla}u_{T}\|^{2}=(\boldsymbol{\nabla}p^{k+1}_{T}\underline{u}_{TF},\,\boldsymbol{\nabla}u_{T})_{T}-\sum_{F\subset\partial T}(u_{F}-u_{T},\,\boldsymbol{\nabla}u_{T}\cdot\boldsymbol{n}_{TF})_{F}=\mathfrak{T}_{1}+\mathfrak{T}_{2}

where the first term $\mathfrak{T}_{1}$ can be readily bounded using Cauchy-Schwarz

(\boldsymbol{\nabla}p^{k+1}_{T}\underline{u}_{TF},\,\boldsymbol{\nabla}u_{T})_{T}\leq\|\boldsymbol{\nabla}p^{k+1}_{T}\underline{u}_{TF}\|_{T}\|\boldsymbol{\nabla}u\|_{T},

and the second term $\mathfrak{T}_{2}$ can be bounded as

	$\displaystyle\mathfrak{T}_{2}$	$\displaystyle\leq\sum_{F}\\|u_{F}-u_{T}\\|_{F}\\|\boldsymbol{\nabla}u\cdot\boldsymbol{n}_{TF}\\|_{F}\leq\left(\sum_{F}h_{F}^{-1}\\|u_{F}-u_{T}\\|_{F}^{2}\right)^{1/2}\!\left(\sum_{F}h_{F}\\|\boldsymbol{\nabla}u_{T}\\|_{F}^{2}\right)^{1/2}$
		$\displaystyle\leq s_{T}(\underline{u}_{TF},\underline{u}_{TF})^{1/2}\left(\sum_{F}h_{F}^{1}h_{F}^{-1}C_{tr}^{2}\\|\boldsymbol{\nabla}u_{T}\\|_{T}^{2}\right)^{1/2}\leq s_{T}(\underline{u}_{TF},\underline{u}_{TF})^{1/2}\,C_{tr}\\|\boldsymbol{\nabla}u_{T}\\|_{T}N_{\partial}^{1/2}$
		$\displaystyle\leq C_{tr}N_{\partial}^{1/2}\,s_{T}(\underline{u}_{TF},\underline{u}_{TF})^{1/2}\\|\boldsymbol{\nabla}u_{T}\\|_{T}.$

Hence, dividing by $\|\boldsymbol{\nabla}u_{T}\|\neq 0$ , we have $\|\boldsymbol{\nabla}u_{T}\|_{T}\lesssim s_{T}(\underline{u}_{TF},\underline{u}_{TF})^{1/2}$ . Putting squares on the above inequality and using $(a+b)^{2}\leq 2a^{2}+2b^{2}$ and adding the stabilization terms, we finally obtain

C\|\underline{u}_{TF}\|_{1,T}^{2}\leq\|\boldsymbol{\nabla}p^{k+1}_{T}\underline{u}_{TF}\|_{T}^{2}+s_{T}(\underline{u}_{TF},\underline{u}_{TF}),

with $C:=\dfrac{1}{\max(2,\,2C_{tr}N_{\partial}^{1/2}+1)}$ , independent of $h$ .

3.4 Discrete Covariant Gradient and Bilinear Form

Let ${\boldsymbol{A}}\in L^{\infty}(\Omega)^{d}$ any smooth vector potential and ${\boldsymbol{A}}_{\chi}\in L^{\infty}(\Omega)^{d}$ its transformation with $\chi\in W^{1,\infty}(\overline{\Omega})$ a real-valued gauge phase. Denoting $\underline{v}_{h}=(\underline{u}_{TF})_{T\in\mathcal{T}_{h}}\in U_{h}^{k}$ an arbitrary DOF, $\underline{v}_{\chi,h}$ and ${\boldsymbol{A}_{T,\chi}}$ the DOFs, and discrete vector potential, after a discrete gauge transformation transformation as

\underline{v}_{\chi,TF}:=\left(\pi_{T}^{k}(e^{i\chi}v_{T}),\pi_{F}^{k}(e^{i\chi}v_{F})\right),\quad{\boldsymbol{A}_{\chi,T}}:=\pi_{T}^{k}{\boldsymbol{A}_{\chi}}.

We define the discrete covariant gradient $\mathrm{G}_{T,{\boldsymbol{A}}}^{k}$ , for any vector potential ${\boldsymbol{A}}$ , as

(12)

(\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{u}_{TF},\boldsymbol{\tau})_{T}=(u_{T},\mathcal{G}_{{\boldsymbol{A}}}^{*}\boldsymbol{\tau})_{T}-i\sum_{F\in\mathcal{F}_{T}}(v_{F},\boldsymbol{\tau}\cdot\boldsymbol{n}_{TF})_{F},\qquad\forall\boldsymbol{\tau}\in\mathbb{P}^{k}(T)^{d},

or equivalently, for any $\forall\boldsymbol{\tau}\in\mathbb{P}^{k}(T)^{d}$ ,

(13)

(\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{u}_{TF},\boldsymbol{\tau})_{T}=(\mathcal{G}_{{\boldsymbol{A}}}u_{T},\boldsymbol{\tau})_{T}-i\sum_{F\in\mathcal{F}_{T}}(v_{T}-v_{F},\boldsymbol{\tau}\cdot\boldsymbol{n}_{TF})_{F}.

The local discrete bilinear form $a_{T}$ , for any cell element $T\in\mathcal{T}_{h}$ , is defined using the $L^{2}$ projection of ${\boldsymbol{A}}$ , as

(14)

a_{T}(\underline{u}_{TF},\underline{v}_{TF};{\boldsymbol{A}}_{T})\mathrel{\mathop{:}}=\bigl(\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{u}_{TF},\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF}\bigr)_{L^{2}(T)}+s_{T}(\underline{u}_{TF},\underline{v}_{TF}),

where $s_{T}(\cdot,\cdot)$ denotes a suitable HHO stabilization, in the sense of [DiPietro2020, Assumption 2.4, p. 49], for instance the choice from [DiPietro2020, Eq. (2.22)]:

s_{T}(\underline{u}_{TF},\underline{v}_{TF}):=\sum_{F\in\mathcal{F}_{T}}h_{F}^{-1}((\delta_{TF}^{k}-\delta_{T}^{k})\underline{u}_{TF},(\delta_{TF}^{k}-\delta_{T}^{k})\underline{v}_{TF})_{F}.

where $\delta_{T}^{k}\underline{u}_{TF}:=\pi_{T}^{k}(p^{k+1}_{T}\underline{v}_{TF}-v_{T})$ and $\delta_{TF}^{k}\underline{u}_{TF}:=\pi_{F}^{k}(p^{k+1}_{T}\underline{v}_{TF}-v_{F})$ for all $F\in\mathcal{F}_{T}$ . The global bilinear form over $\mathcal{T}_{h}$ is defined as the sum of all the local contributions (14), i.e.

(15)

\displaystyle a_{h}(\underline{u}_{h},\underline{v}_{h};{\boldsymbol{A}_{h}})

\displaystyle:=\sum_{T\in\mathcal{T}_{h}}a_{T}(\underline{u}_{TF},\underline{v}_{TF};{\boldsymbol{A}}_{T}),

where ${\boldsymbol{A}_{h}}:=\pi_{h}^{k}({\boldsymbol{A}})=(\pi_{T}^{k}{\boldsymbol{A}})_{T\in\mathcal{T}_{h}}$ is the $L^{2}$ projection of ${\boldsymbol{A}}$ over all $T\in\mathcal{T}_{h}$ .

Theorem 3.4 (Discrete Gauge Covariance).

Let $\chi\in C^{\infty}(\overline{\Omega})$ be a smooth real-valued gauge phase, and $\mathrm{G}_{T,\chi}^{k}:=\mathrm{G}_{T,{\boldsymbol{A}}_{\chi,T}}^{k}$ . Then, for any $\underline{v}_{h}\in U_{h}^{k}$ , we have

(\mathrm{G}_{T,\chi}^{k}(\underline{v}_{\chi,TF}),\boldsymbol{\tau})_{T}=(\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF},\pi_{T}^{k}(e^{-i\chi}\boldsymbol{\tau}))_{T}+\mathcal{R}_{\mathrm{gauge}}(\underline{v}_{TF},\boldsymbol{\tau})\qquad\forall\boldsymbol{\tau}\in\mathbb{P}^{k}(T)^{d},

where the remainder $\mathcal{R}_{\mathrm{gauge}}$ gathers all projection commutators and is bounded by:

|\mathcal{R}_{\mathrm{gauge}}(\underline{v}_{TF},\boldsymbol{\tau})|\leq Ch_{T}^{k+1}\|\underline{v}_{TF}\|_{1,T}\|e^{-i\chi}\boldsymbol{\tau}\|_{H^{k+1}(T)}.

Proof 3.5.

Let $\boldsymbol{\tau}\in\mathbb{P}^{k}(T)^{d}$ . By definition of the discrete gradient $\mathrm{G}_{T,\chi}^{k}$ with the projected potential ${\boldsymbol{A}}_{\chi,T}$ , we have for the first term:

(16)		$\displaystyle(\mathrm{G}_{T,\chi}^{k}\underline{v}_{\chi,TF},\boldsymbol{\tau})_{T}$	$\displaystyle=(\pi_{T}^{k}(e^{i\chi}v_{T}),\mathcal{G}_{{{\boldsymbol{A}}_{\chi,T}}}^{*}\boldsymbol{\tau})_{T}-i\sum_{F\in\mathcal{F}_{T}}(\pi_{F}^{k}(e^{i\chi}v_{F}),\boldsymbol{\tau}\cdot\boldsymbol{n}_{TF})_{F}$
(17)			$\displaystyle=(\pi_{T}^{k}(e^{i\chi}v_{T}),\mathcal{G}_{{{\boldsymbol{A}}_{\chi,T}}}^{*}\boldsymbol{\tau})_{T}-i\sum_{F\in\mathcal{F}_{T}}(e^{i\chi}v_{F},\boldsymbol{\tau}\cdot\boldsymbol{n}_{TF})_{F}$

where we have used the fact that $\boldsymbol{\tau}\cdot\boldsymbol{n}_{TF}\in\mathbb{P}^{k}(F)$ , so that the face projector $\pi_{F}^{k}$ can be exactly dropped by orthogonality. For the volumetric term, we add and subtract the continuous potential ${\boldsymbol{A}_{\chi}}$ to reconstruct the continuous adjoint $\mathcal{G}_{{{\boldsymbol{A}_{\chi}}}}^{*}\boldsymbol{\tau}=-i\boldsymbol{\nabla}{\cdot}\boldsymbol{\tau}-{\boldsymbol{A}_{\chi}}\cdot\boldsymbol{\tau}$ :

\mathcal{G}_{{{\boldsymbol{A}}_{\chi,T}}}^{*}\boldsymbol{\tau}=-i\boldsymbol{\nabla}{\cdot}\boldsymbol{\tau}-{\boldsymbol{A}}_{\chi,T}\cdot\boldsymbol{\tau}=\mathcal{G}_{{{\boldsymbol{A}_{\chi}}}}^{*}\boldsymbol{\tau}+({\boldsymbol{A}_{\chi}}-{\boldsymbol{A}}_{\chi,T})\cdot\boldsymbol{\tau}.

Substituting the above equality in (17) yields the decomposition:

(18)

(\mathrm{G}_{T,\chi}^{k}\underline{v}_{\chi,TF},\boldsymbol{\tau})_{T}=(\pi_{T}^{k}(e^{i\chi}v_{T}),\mathcal{G}_{{{\boldsymbol{A}_{\chi}}}}^{*}\boldsymbol{\tau})_{T}+(\pi_{T}^{k}(e^{i\chi}v_{T}),({\boldsymbol{A}_{\chi}}-{\boldsymbol{A}}_{\chi,T})\cdot\boldsymbol{\tau})_{T}\\ -i\sum_{F\in\mathcal{F}_{T}}(e^{i\chi}v_{F},\boldsymbol{\tau}\cdot\boldsymbol{n}_{TF})_{F}.

Adding and removing $e^{i\chi}v_{T}$ in the first volumetric term in the RHS gives

	$\displaystyle(\mathrm{G}_{T,\chi}^{k}\underline{v}_{\chi,TF},\boldsymbol{\tau})_{T}=\;$	$\displaystyle(e^{i\chi}v_{T},\mathcal{G}_{{{\boldsymbol{A}_{\chi}}}}^{*}\boldsymbol{\tau})_{T}-i\sum_{F\in\mathcal{F}_{T}}(e^{i\chi}v_{F},\boldsymbol{\tau}\cdot\boldsymbol{n}_{TF})_{F}$
		$\displaystyle+(\pi_{T}^{k}(e^{i\chi}v_{T}),({\boldsymbol{A}_{\chi}}-{\boldsymbol{A}}_{\chi,T})\cdot\boldsymbol{\tau})_{T}$
		$\displaystyle+(\pi_{T}^{k}(e^{i\chi}v_{T})-e^{i\chi}v_{T},\mathcal{G}_{{{\boldsymbol{A}_{\chi}}}}^{*}\boldsymbol{\tau})_{T}=:\;\mathfrak{T}_{1}+\mathfrak{T}_{2}+\mathfrak{T}_{3}.$

We now show that $\mathfrak{T}_{1}$ simplifies. From Lemma 2.5, we have $\mathcal{G}_{{{\boldsymbol{A}_{\chi}}}}^{*}\boldsymbol{\tau}=e^{i\chi}\mathcal{G}_{{\boldsymbol{A}}}^{*}(e^{-i\chi}\boldsymbol{\tau})$ , the global phase cancels perfectly in the complex inner product:

	$\displaystyle\mathfrak{T}_{1}$	$\displaystyle=(e^{i\chi}v_{T},e^{i\chi}\mathcal{G}_{{\boldsymbol{A}}}^{*}(e^{-i\chi}\boldsymbol{\tau}))_{T}-i\sum_{F\in\mathcal{F}_{T}}(e^{i\chi}v_{F},e^{i\chi}(e^{-i\chi}\boldsymbol{\tau})\cdot\boldsymbol{n}_{TF})_{F}$
		$\displaystyle=(v_{T},\mathcal{G}_{{\boldsymbol{A}}}^{*}(e^{-i\chi}\boldsymbol{\tau}))_{T}-i\sum_{F\in\mathcal{F}_{T}}(v_{F},(e^{-i\chi}\boldsymbol{\tau})\cdot\boldsymbol{n}_{TF})_{F}.$

Let $\boldsymbol{E}:=e^{-i\chi}\boldsymbol{\tau}-\pi_{T}^{k}(e^{-i\chi}\boldsymbol{\tau})$ be the projection error of the gauged test function. Substituting $e^{-i\chi}\boldsymbol{\tau}=\pi_{T}^{k}(e^{-i\chi}\boldsymbol{\tau})+\boldsymbol{E}$ into $\mathfrak{T}_{1}$ , we separate the polynomial part from the remainder:

	$\displaystyle\mathfrak{T}_{1}=\;$	$\displaystyle(v_{T},\mathcal{G}_{{\boldsymbol{A}}}^{*}\pi_{T}^{k}(e^{-i\chi}\boldsymbol{\tau}))_{T}-i\sum_{F\in\mathcal{F}_{T}}(v_{F},\pi_{T}^{k}(e^{-i\chi}\boldsymbol{\tau})\cdot\boldsymbol{n}_{TF})_{F}$
		$\displaystyle+(v_{T},\mathcal{G}_{{\boldsymbol{A}}}^{*}\boldsymbol{E})_{T}-i\sum_{F\in\mathcal{F}_{T}}(v_{F},\boldsymbol{E}\cdot\boldsymbol{n}_{TF})_{F}.$

By the exact definition of the discrete gradient $\mathrm{G}_{T,{\boldsymbol{A}}}^{k}$ given in (12), the first line is exactly equal to $(\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF},\pi_{T}^{k}(e^{-i\chi}\boldsymbol{\tau}))_{T}$ . Thus, we have:

\mathfrak{T}_{1}=(\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF},\pi_{T}^{k}(e^{-i\chi}\boldsymbol{\tau}))_{T}+\mathfrak{T}_{\boldsymbol{E}},

where $\mathfrak{T}_{\boldsymbol{E}}$ encapsulates the action of the continuous operator on the projection error $\boldsymbol{E}$ :

\mathfrak{T}_{\boldsymbol{E}}:=(v_{T},\mathcal{G}_{{\boldsymbol{A}}}^{*}\boldsymbol{E})_{T}-i\sum_{F\in\mathcal{F}_{T}}(v_{F},\boldsymbol{E}\cdot\boldsymbol{n}_{TF})_{F}.

Denoting the total remainder $\mathcal{R}_{\mathrm{gauge}}(\underline{v}_{TF},\boldsymbol{\tau}):=\mathfrak{T}_{2}+\mathfrak{T}_{3}+\mathfrak{T}_{\boldsymbol{E}}$ , it remains to bound these terms. For $\mathfrak{T}_{2}$ and $\mathfrak{T}_{3}$ , standard approximation properties of the $L^{2}$ -projector yield the optimal bound $\mathcal{O}(h_{T}^{k+1})\|v_{T}\|_{T}\|\boldsymbol{\tau}\|_{T}$ . The term $\mathfrak{T}_{\boldsymbol{E}}$ is treated by applying integration by parts to the continuous adjoint $\mathcal{G}_{{\boldsymbol{A}}}^{*}\boldsymbol{E}=-i\boldsymbol{\nabla}{\cdot}\boldsymbol{E}-{\boldsymbol{A}}\cdot\boldsymbol{E}$ :

	$\displaystyle\mathfrak{T}_{\boldsymbol{E}}$	$\displaystyle=(v_{T},-i\boldsymbol{\nabla}{\cdot}\boldsymbol{E})_{T}-(v_{T}{\boldsymbol{A}},\boldsymbol{E})_{T}-i\sum_{F\in\mathcal{F}_{T}}(v_{F},\boldsymbol{E}\cdot\boldsymbol{n}_{TF})_{F}$
		$\displaystyle=(-i\boldsymbol{\nabla}v_{T},\boldsymbol{E})_{T}-(v_{T}{\boldsymbol{A}},\boldsymbol{E})_{T}+i\sum_{F\in\mathcal{F}_{T}}\langle v_{T}-v_{F},\boldsymbol{E}\cdot\boldsymbol{n}_{TF}\rangle_{F}.$

Because $v_{T}\in\mathbb{P}^{k}(T)$ , its gradient $\boldsymbol{\nabla}v_{T}$ is a polynomial of degree $k-1$ . Since $\boldsymbol{E}$ is the error of the $L^{2}$ -orthogonal projection onto $\mathbb{P}^{k}(T)^{d}$ , the term $(-i\boldsymbol{\nabla}v_{T},\boldsymbol{E})_{T}$ vanishes identically by orthogonality. We are left with purely $L^{2}$ terms which can be estimated:

	$\displaystyle\|\mathfrak{T}_{\boldsymbol{E}}\|$	$\displaystyle\leq\\|{\boldsymbol{A}}\\|_{L^{\infty}(T)}\\|v_{T}\\|_{T}\\|\boldsymbol{E}\\|_{T}+\sum_{F\in\mathcal{F}_{T}}\\|v_{T}-v_{F}\\|_{F}\\|\boldsymbol{E}\\|_{F}$
		$\displaystyle\leq C\\|{\boldsymbol{A}}\\|_{L^{\infty}(T)}\\|v_{T}\\|_{T}\left(h_{T}^{k+1}\|e^{-i\chi}\boldsymbol{\tau}\|_{H^{k+1}(T)}\right)+\sum_{F\in\mathcal{F}_{T}}\|v_{h}\|_{1,T}h_{F}^{\frac{1}{2}}\left(h_{T}^{k+\frac{1}{2}}\|e^{-i\chi}\boldsymbol{\tau}\|_{H^{k+1}(T)}\right)$
		$\displaystyle\leq Ch_{T}^{k+1}\|v_{h}\|_{1,T}\\|e^{-i\chi}\boldsymbol{\tau}\\|_{H^{k+1}(T)},$

which concludes the proof.

Remark 3.6.

By employing the adjoint operator $(\pi_{T}^{k})^{*}$ and defining ${\boldsymbol{A}}_{\chi,T}:=\pi_{T}^{k}({\boldsymbol{A}_{\chi}})$ , the preceding identity can be recast as

\mathrm{G}_{T,{\boldsymbol{A}}_{\chi,T}}^{k}\,\underline{v}_{\chi,TF}=e^{i\chi}\,(\pi_{T}^{k})^{*}\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\,\underline{v}_{TF}+\mathcal{O}(h^{k+1}),

which parallels the continuous counterpart $\mathcal{G}_{{\boldsymbol{A}_{\chi}}}\,\psi_{\chi}=e^{i\chi}\,\mathcal{G}_{{\boldsymbol{A}}}\psi$ .

Remark 3.7.

Setting $\boldsymbol{\tau}_{h}:=\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{h}\in\mathbb{P}^{k}(T)^{d}$ and evaluating the $(k+1)$ -th derivatives of the gauged test function via the Leibniz rule yields:

|e^{-i\chi}\boldsymbol{\tau}_{h}|_{H^{k+1}(T)}\leq C_{\chi}\sum_{j=1}^{k+1}\|\boldsymbol{\nabla}^{k+1-j}\boldsymbol{\tau}_{h}\|_{T}.

Since $\boldsymbol{\tau}_{h}$ is a polynomial of degree $k$ , its derivatives up to order $k$ are non-zero. Using the inverse inequality yields $\|\boldsymbol{\nabla}^{m}\boldsymbol{\tau}_{h}\|_{T}\leq C_{inv}h_{T}^{-m}\|\boldsymbol{\tau}_{h}\|_{T}$ . The dominant term corresponds to $j=1$ (the $k$ -th derivative of $\boldsymbol{\tau}_{h}$ ), which introduces a factor $h_{T}^{-k}$ . Consequently, the general algebraic covariance bound degrades to first-order:

|\mathcal{R}_{\mathrm{gauge}}(\underline{v}_{h},\boldsymbol{\tau}_{h})|\lesssim h_{T}^{k+1}\left(h_{T}^{-k}\|\boldsymbol{\tau}_{h}\|_{T}\right)\|\underline{v}_{h}\|_{1,T}=\mathcal{O}(h)\|\underline{v}_{h}\|_{1,T}^{2}.

However, the practical manifestation of this $\mathcal{O}(h)$ limit depends heavily on the physical smoothness of the targeted eigenstate. For the fundamental state and the first few excited states, the physical wavefunctions oscillate on a macroscopic scale $\lambda_{\ell}$ that is much larger than the mesh size ( $\lambda_{\ell}\gg h$ ). The discrete gradient $\boldsymbol{\tau}_{h}$ approximates a highly smooth continuous field. The physical spatial derivatives scale as $||\boldsymbol{\tau}_{h}||_{H^{k+1}(T)}\sim(\lambda_{\ell})^{-m}\|\boldsymbol{\tau}_{h}\|_{T}$ , which is practically bounded independently of $h$ , and the scheme exhibits apparent optimal super-convergence $\mathcal{O}(h^{2k+2})$ on the eigenvalues, effectively bypassing the local $\mathcal{O}(h)$ aliasing limit. For highly excited states, the physical wavelength approaches the grid resolution ( $\lambda_{\ell}\approx h$ ). The eigenfunctions oscillate at the mesh scale, meaning the inverse inequality bound is physically saturated: the spatial derivatives genuinely scale as $h_{T}^{-m}$ . In this regime, the polynomial space is unable to adequately resolve the combined oscillation of the phase $e^{-i\chi}$ and the state. The discrete gauge error strictly degrades to $\mathcal{O}(h)$ , requiring a fine mesh to recover asymptotic coherence. This highlights that strong gauge invariance is lost at the algebraic level due to polynomial aliasing, but remains asymptotically optimal for macroscopic low-energy physical observables.

Corollary 3.8 (Consistency with continuous potential).

Let $\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}$ be the discrete gradient constructed with the projected potential ${\boldsymbol{A}}_{T}=\pi_{T}^{k}({\boldsymbol{A}})$ . For any $\underline{v}_{TF}\in U_{T}^{k}$ , we have

(\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF},\boldsymbol{\tau})_{T}=(\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF},\boldsymbol{\tau})_{T}+\mathcal{O}(h_{T}^{k+1})\|\underline{v}_{TF}\|_{1,T}\|\boldsymbol{\tau}\|_{T},\qquad\forall\boldsymbol{\tau}\in\mathbb{P}^{k}(T)^{d}.

Proof 3.9.

Applying Theorem 3.4 with the trivial gauge phase $\chi=0$ gives $\underline{v}_{0,TF}=\underline{v}_{TF}$ and ${\boldsymbol{A}}_{0,T}={\boldsymbol{A}}_{T}$ . Furthermore, since $e^{-i0}\boldsymbol{\tau}=\boldsymbol{\tau}\in\mathbb{P}^{k}(T)^{d}$ , its $L^{2}$ -projection is exact, meaning the projection error vanishes identically ( $\boldsymbol{E}=0$ ). The total remainder $\mathcal{R}_{\mathrm{gauge}}(\underline{v}_{TF},\boldsymbol{\tau})$ reduces solely to the potential approximation term $\mathfrak{T}_{2}=(v_{T},({\boldsymbol{A}}-{\boldsymbol{A}}_{T})\cdot\boldsymbol{\tau})_{T}$ . By standard polynomial approximation, this is bounded by $Ch_{T}^{k+1}|{\boldsymbol{A}}|_{W^{k+1,\infty}(T)}\|v_{T}\|_{T}\|\boldsymbol{\tau}\|_{T}$ without requiring any inverse inequality on the test function, yielding the optimal $L^{2}$ bound.

3.5 Discrete Gårding Inequality

Before addressing the discrete Gårding Inequality, we need the following Lemma.

Lemma 3.10.

There exists a constant $C>0$ , independent of $h$ , such that the magnetic-free discrete covariante gradient $\mathrm{G}_{T,\boldsymbol{0}}^{k}$ satisfies

||\mathrm{G}_{T,\boldsymbol{0}}^{k}\underline{u}_{TF}||_{T}^{2}+s_{T}(\underline{u}_{TF},\underline{u}_{TF})\geq C||\underline{u}_{TF}||_{1,T}^{2},

for any $\underline{u}_{TF}\in U_{T}^{k}$ .

Proof 3.11.

taking $\boldsymbol{\tau}=\boldsymbol{\nabla}w$ with $w\in\mathbb{P}^{k+1}(T)$ and ${\boldsymbol{A}}=\boldsymbol{0}$ in the definition of $\mathrm{G}_{T,{\boldsymbol{A}}}^{k}$ (13) yields

	$\displaystyle(\mathrm{G}_{T,\boldsymbol{0}}^{k}\underline{u}_{TF},\boldsymbol{\nabla}w)_{T}$	$\displaystyle=-i\Big((\boldsymbol{\nabla}u,\boldsymbol{\nabla}w)_{T}+\sum_{F\in\mathcal{F}_{T}}(u_{T}-u_{F},\boldsymbol{\nabla}w\cdot\boldsymbol{n}_{TF})_{F}\Big)$
		$\displaystyle=(-i\boldsymbol{\nabla}p^{k+1}_{T}\underline{u}_{TF},\boldsymbol{\nabla}w)_{T}.$

This indicates that $-i\boldsymbol{\nabla}p^{k+1}_{T}\underline{u}_{TF}$ is the $L^{2}$ projection of $\mathrm{G}_{T,\boldsymbol{0}}^{k}\underline{u}_{TF}$ on the subspace $\boldsymbol{\nabla}\mathbb{P}^{k+1}(T)\subset\mathbb{P}^{k}(T)^{d}$ . Hence,

||\mathrm{G}_{T,\boldsymbol{0}}^{k}\underline{u}_{TF}||^{2}_{T}\geq||-i\boldsymbol{\nabla}p^{k+1}_{T}\underline{u}_{TF}||^{2}\geq||\boldsymbol{\nabla}p^{k+1}_{T}\underline{u}_{TF}||.

Adding $s_{T}(\underline{u}_{TF},\underline{u}_{TF})$ on both sides and using Proposition 3.2 yields

	$\displaystyle\|\|\mathrm{G}_{T,\boldsymbol{0}}^{k}\underline{u}_{TF}\|\|^{2}_{T}+s_{T}(\underline{u}_{TF},\underline{u}_{TF})$	$\displaystyle\geq C\|\|\boldsymbol{\nabla}u_{T}\|\|_{T}^{2}+s_{T}(\underline{u}_{TF},\underline{u}_{TF}),$
		$\displaystyle\geq\min(1,C)\|\|\underline{u}_{TF}\|\|_{1,T}^{2}.$

Theorem 3.12 (Discrete Gauge-Invariant Gårding Inequality).

There exist constants $\alpha>0$ and $\gamma>0$ , independent of $h$ , such that for all $\underline{v}_{h}\in U_{h}^{k}$ :

a_{h}(\underline{v}_{h},\underline{v}_{h};{\boldsymbol{A}})\geq\alpha\|\underline{v}_{h}\|_{1,h}^{2}-\gamma\left(\|{\boldsymbol{A}}\|_{L^{\infty}(\Omega)}^{2}+\|V\|_{L^{\infty}(\Omega)}\right)\|\underline{v}_{h}\|_{L^{2}(\Omega)}^{2},

where $\|\underline{v}_{h}\|_{L^{2}(\Omega)}^{2}:=\sum_{T\in\mathcal{T}_{h}}\|v_{T}\|_{T}^{2}$ .

Proof 3.13.

By the definition of the full covariant gradient $\mathrm{G}_{T,{\boldsymbol{A}}}^{k}$ :

	$\displaystyle(\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF},\boldsymbol{\tau})_{T}$	$\displaystyle=((-i\boldsymbol{\nabla}-{\boldsymbol{A}})v_{T},\boldsymbol{\tau})_{T}-i\sum_{F\in\mathcal{F}_{T}}\langle v_{T}-v_{F},\boldsymbol{\tau}\cdot\boldsymbol{n}_{TF}\rangle_{F}$
		$\displaystyle=(\mathrm{G}_{T,0}^{k}\underline{v}_{TF},\boldsymbol{\tau})_{T}-({\boldsymbol{A}}v_{T},\boldsymbol{\tau})_{T}.$

Since this equality holds for any test polynomial $\boldsymbol{\tau}\in\mathbb{P}^{k}(T)^{d}$ , it yields the exact algebraic identity:

\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF}=\mathrm{G}_{T,\boldsymbol{0}}^{k}\underline{v}_{TF}-\pi_{T}^{k}({\boldsymbol{A}}v_{T}).

Taking the squared $L^{2}$ -norm of $\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF}$ and using the elementary inequality $\|a-b\|^{2}\geq\frac{1}{2}\|a\|^{2}-\|b\|^{2}$ , we get:

\|\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF}\|_{T}^{2}\geq\frac{1}{2}\|\mathrm{G}_{T,\boldsymbol{0}}^{k}\underline{v}_{TF}\|_{T}^{2}-\|\pi_{T}^{k}({\boldsymbol{A}}v_{T})\|_{T}^{2}.

Since the $L^{2}$ -orthogonal projector $\pi_{T}^{k}$ is a contraction (i.e., $\|\pi_{T}^{k}\boldsymbol{\tau}\|_{T}\leq\|\boldsymbol{\tau}\|_{T}$ ), we can bound the magnetic perturbation as:

\|\pi_{T}^{k}({\boldsymbol{A}}v_{T})\|_{T}^{2}\leq\|{\boldsymbol{A}}v_{T}\|_{T}^{2}\leq\|{\boldsymbol{A}}\|_{L^{\infty}(T)}^{2}\|v_{T}\|_{T}^{2}.

Adding the stabilization term $s_{T}(\underline{v}_{TF},\underline{v}_{TF})$ to both sides (and noting that $\frac{1}{2}s_{T}\leq s_{T}$ ), we obtain:

\|\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF}\|_{T}^{2}+s_{T}(\underline{v}_{TF},\underline{v}_{TF})\geq\frac{1}{2}\left(\|\mathrm{G}_{T,\boldsymbol{0}}^{k}\underline{v}_{TF}\|_{T}^{2}+s_{T}(\underline{v}_{TF},\underline{v}_{TF})\right)-\|{\boldsymbol{A}}\|_{L^{\infty}(T)}^{2}\|v_{T}\|_{T}^{2}.

Using the upper bound in Lemma 3.10 yields:

\|\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF}\|_{T}^{2}+s_{T}(\underline{v}_{TF},\underline{v}_{TF})\geq\frac{C}{2}\|\underline{v}_{TF}\|_{1,T}^{2}-\|{\boldsymbol{A}}\|_{L^{\infty}(T)}^{2}\|v_{T}\|_{T}^{2}.

The local bilinear form is given by $a_{T}(\underline{v}_{TF},\underline{v}_{TF})=\|\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF}\|_{T}^{2}+s_{T}(\underline{v}_{TF},\underline{v}_{TF})+(V\pi_{T}^{k}v_{T},v_{T})_{T}$ . For the scalar potential term, the contraction property of $\pi_{T}^{k}$ yields:

(V\pi_{T}^{k}v_{T},v_{T})_{T}\geq-\|V\|_{L^{\infty}(T)}\|\pi_{T}^{k}v_{T}\|_{T}\|v_{T}\|_{T}\geq-\|V\|_{L^{\infty}(T)}\|v_{T}\|_{T}^{2}.

Summing over all elements $T\in\mathcal{T}_{h}$ , we finally obtain:

	$\displaystyle a_{h}(\underline{v}_{TF},\underline{v}_{TF})$	$\displaystyle=\sum_{T\in\mathcal{T}_{h}}a_{T}(\underline{v}_{TF},\underline{v}_{TF})\geq\sum_{T\in\mathcal{T}_{h}}\left(\frac{C}{2}\\|\underline{v}_{TF}\\|_{1,T}^{2}-\left(\\|{\boldsymbol{A}}\\|_{L^{\infty}(T)}^{2}+\\|V\\|_{L^{\infty}(T)}\right)\\|v_{T}\\|_{T}^{2}\right)$
		$\displaystyle\geq\frac{C}{2}\\|\underline{v}_{TF}\\|_{1,h}^{2}-\left(\\|{\boldsymbol{A}}\\|_{L^{\infty}(\Omega)}^{2}+\\|V\\|_{L^{\infty}(\Omega)}\right)\\|\underline{v}_{h}\\|_{L^{2}(\Omega)}^{2}.$

Setting $\alpha=\frac{C}{2}$ and $\gamma=1$ concludes the proof.

3.6 A Priori Error Estimate for the Stationary Schröedinger problem

In this section, we focus on the stationary Schrödinger problem: find $u\in H^{1}_{0}(\Omega;\mathbb{C})$ , such that

(19)

\mathcal{G}_{{\boldsymbol{A}}}^{*}\mathcal{G}_{{\boldsymbol{A}}}u+Vu=f,\qquad\text{ in }\Omega.

In particular, we establish a priori error estimates for the discrete equivalent of the above problem : find $\underline{u}_{h}\in U_{h}^{k}$ such that

(20)

\displaystyle a_{h}(\underline{u}_{h},\underline{v}_{h};{\boldsymbol{A}_{h}})

\displaystyle=\ell(\underline{v}_{h}),\qquad\forall\underline{v}_{h}\in U_{h}^{k},

where $a_{h}(\cdot,\cdot;{\boldsymbol{A}_{h}})$ is defined using $a_{T}$ defined in (14) with the potential $V$ , i.e.

a_{h}(\underline{u}_{h},\underline{v}_{h};{\boldsymbol{A}_{h}}):=\sum_{T\in\mathcal{T}_{h}}\left\{a_{T}(\underline{u}_{TF},\underline{v}_{TF})+(Vu_{T},v_{T})_{T}\right\}.

We also assume that the mesh is sufficiently fine so that this problem is well-posed.

Theorem 3.14 (Optimal a priori error estimate).

Let $u\in H_{0}^{1}(\Omega;\mathbb{C})$ be the exact solution of the continuous problem (19) and $\underline{u}_{h}\in U_{h}^{k}$ be the discrete HHO solution from (20). Assuming sufficient regularity $u\in H^{k+2}(\Omega;\mathbb{C})$ , ${\boldsymbol{A}}\in W^{k+1,\infty}(\Omega;\mathbb{R}^{d})$ , and $V\in L^{\infty}(\Omega;\mathbb{R})$ , there exists a constant $C>0$ , independent of $h$ , ${\boldsymbol{A}}$ , and $V$ , such that for $h$ sufficiently small:

(21)

\|\underline{u}_{h}-I^{k}_{h}u\|_{1,h}\leq Ch^{k+1}\left(\|u\|_{H^{k+2}(\Omega)}+\big(\|{\boldsymbol{A}}\|_{W^{k+1,\infty}(\Omega)}+\|V\|_{L^{\infty}(\Omega)}\big)\|u\|_{H^{k+1}(\Omega)}\right).

Moreover, an $L^{2}$ -error bound for the cell unknowns can be given as

(22)

\|\underline{u}_{h}-\pi_{h}^{k}u\|_{L^{2}(\Omega)}\leq Ch^{k+2}\left(\|u\|_{H^{k+2}(\Omega)}+\big(\|{\boldsymbol{A}}\|_{W^{k+1,\infty}(\Omega)}+\|V\|_{L^{\infty}(\Omega)}\big)\|u\|_{H^{k+1}(\Omega)}\right).

Proof 3.15 (Proof of Theorem 3.14).

First, we observe that there exists a constant $\beta>0$ , independent of $h$ , such that for all $\underline{w}_{h}\in U_{h}^{k}$ :

\beta\|\underline{w}_{h}\|_{1,h}\leq\sup_{\underline{v}_{h}\in U_{h}^{k}\setminus\{0\}}\frac{|a_{h}(\underline{w}_{h},\underline{v}_{h})|}{\|\underline{v}_{h}\|_{1,h}}.

Indeed, the above discrete inf-sup condition follows from the discrete Gårding inequality via a standard compactness argument in HHO spaces (see [DiPietro2020, Theorem 6.41] with $p=2$ ). Let $\underline{e}_{h}:=\underline{u}_{h}-I^{k}_{h}u\in U_{h}^{k}$ be the discrete error. Applying the inf-sup condition to $\underline{e}_{h}$ yields:

\beta\|\underline{e}_{h}\|_{1,h}\leq\sup_{\underline{v}_{h}\in U_{h}^{k}\setminus\{0\}}\frac{|a_{h}(\underline{u}_{h},\underline{v}_{h})-a_{h}(I^{k}_{h}u,\underline{v}_{h})|}{\|\underline{v}_{h}\|_{1,h}}.

Since $\underline{u}_{h}$ is the discrete solution, it satisfies $a_{h}(\underline{u}_{h},\underline{v}_{h})=\ell(\underline{v}_{h})$ . We define the consistency error functional as $\mathcal{E}_{h}(\underline{v}_{h}):=\ell(\underline{v}_{h})-a_{h}(I^{k}_{h}u,\underline{v}_{h})$ , for any $\underline{v}_{h}\in U_{h}^{k}$ . Since $u$ solves the exact continuous equation (19), we have the following splitting

\mathcal{E}_{h}(\underline{v}_{h})=(f,v_{h})-a_{h}(I^{k}_{h}u,\underline{v}_{h})=:\mathfrak{T}_{grad}+\mathfrak{T}_{stab}+\mathfrak{T}_{pot},

where $\mathfrak{T}_{pot}:=-\sum_{T\in\mathcal{T}_{h}}(V\pi_{T}^{k}u,v_{T})_{T}$ , $\mathfrak{T}_{stab}:=-\sum_{T\in\mathcal{T}_{h}}s_{T}(I^{k}_{T}u,\underline{v}_{TF})$ and

\mathfrak{T}_{grad}:=\sum_{T\in\mathcal{T}_{h}}(\mathcal{G}_{{\boldsymbol{A}}}^{*}\mathcal{G}_{{\boldsymbol{A}}}u+Vu,v_{T})_{T}-(\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}I^{k}_{T}u,\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF})_{T}.

Step 1. We first bound the term $\mathfrak{T}_{pot}$ , by Cauchy-Schwarz, the definition of $h$ , the local approximation property of $\pi_{T}^{k}$ in (30), the global discrete Poincaré inequality in Lemma A.5 for $\underline{v}_{TF}\in U_{T}^{k}$ and the boundedness of the potential $V$ yields

	$\displaystyle\|\mathfrak{T}_{pot}\|$	$\displaystyle=\left\|\sum_{T\in\mathcal{T}_{h}}(Vu-V\pi_{T}^{k}u,v_{T})_{T}\right\|\leq\sum_{T\in\mathcal{T}_{h}}\|\|V\|\|_{L^{\infty}(\Omega)}\\|u-\pi_{T}^{k}u\\|_{T}\\|v_{T}\\|_{T}$
		$\displaystyle\leq\|\|V\|\|_{L^{\infty}(\Omega)}\sum_{T}Ch_{T}^{k+1}\|u\|_{H^{k+2}(T)}\\|v_{T}\\|_{T}$
		$\displaystyle=\|\|V\|\|_{L^{\infty}(\Omega)}Ch^{k+1}\|u\|_{H^{k+2}(\Omega)}\\|v_{h}\\|_{L^{2}(\Omega)}$
		$\displaystyle\lesssim\|\|V\|\|_{L^{\infty}(\Omega)}h^{k+1}\|u\|_{H^{k+2}(\Omega)}\\|\underline{v}_{h}\\|_{1,h}.$

Step 2. Then, for $\mathfrak{T}_{stab}$ , we use Cauchy-Schwarz and Proposition A.1

	$\displaystyle\|\mathfrak{T}_{stab}\|$	$\displaystyle=\left\|\sum_{T\in\mathcal{T}_{h}}s_{T}(I^{k}_{T}u,\underline{v}_{TF})\right\|=\left\|\sum_{T\in\mathcal{T}_{h}}\sum_{F\in\mathcal{F}_{T}}h_{F}^{-1}(\delta_{TF}^{k}(u)-\delta_{T}^{k}(u),v_{F}-v_{T})_{F}\right\|$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}\sum_{F\in\mathcal{F}_{T}}h_{F}^{-\tfrac{1}{2}}\\|\delta_{TF}^{k}(u)-\delta_{T}^{k}(u)\\|_{F}h_{F}^{-\tfrac{1}{2}}\\|v_{F}-v_{T}\\|_{F}$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}\left(\sum_{F\in\mathcal{F}_{T}}h_{F}^{-1}\\|\delta_{TF}^{k}(u)-\delta_{T}^{k}(u)\\|_{F}^{2}\right)^{\tfrac{1}{2}}\left(\sum_{F\in\mathcal{F}_{T}}h_{F}^{-1}\\|v_{F}-v_{T}\\|_{F}^{2}\right)^{\tfrac{1}{2}}$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}s_{T}(I^{k}_{T}u,I^{k}_{T}u)^{\tfrac{1}{2}}s_{T}(\underline{v}_{TF},\underline{v}_{TF})^{\tfrac{1}{2}}$
		$\displaystyle\leq Ch^{k+1}\sum_{T\in\mathcal{T}_{h}}\\|u\\|_{H^{k+2}(T)}\\|\underline{v}_{TF}\\|_{1,T}\lesssim h^{k+1}\|u\|_{H^{k+2}(\Omega)}\\|\underline{v}_{h}\\|_{1,h}$

Step 3. The last term $\mathfrak{T}_{grad}=\sum_{T\in\mathcal{T}_{h}}M_{T}$ reads as a sum of three terms

\small M_{T}=(\mathcal{G}_{{\boldsymbol{A}}}u,\mathcal{G}_{{\boldsymbol{A}}}v_{T})_{T}-i\sum_{F\in\mathcal{F}_{T}}(\mathcal{G}_{{\boldsymbol{A}}}u\cdot\boldsymbol{n}_{TF},v_{T}-v_{F})_{F}-(\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}I^{k}_{T}u,\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF})_{T}.

We first remark that, since $u$ is a regular function, we have $\pi_{F}^{k}(u)=\pi_{T}^{k}(u)_{|F}$ for any $F\in\mathcal{F}_{h}$ , hence

(\mathrm{G}_{T,{\boldsymbol{A}}}^{k}I^{k}_{T}u,\boldsymbol{\tau})_{T}=(\mathcal{G}_{{\boldsymbol{A}}}u,\boldsymbol{\tau})_{T}-i\sum_{F\in\mathcal{F}_{T}}(\pi_{F}^{k}(u)-\pi_{T}^{k}(u),\boldsymbol{\tau}\cdot\boldsymbol{n}_{TF})_{F}=(\mathcal{G}_{{\boldsymbol{A}}}u,\boldsymbol{\tau})_{T},

for any $\boldsymbol{\tau}\in\mathbb{P}^{k}(T)^{d}$ , in other terms

(23)

\mathrm{G}_{T,{\boldsymbol{A}}}^{k}I^{k}_{T}u=\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),\quad\forall{\boldsymbol{A}}\in W^{k+1,\infty}(\Omega).

Replacing $\mathcal{G}_{{\boldsymbol{A}}}u$ with $\mathcal{G}_{{\boldsymbol{A}}}u-\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u)+\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u)$ in the volumetric term and the face terms, yields

	$\displaystyle M_{T}$	$\displaystyle=(\mathcal{G}_{{\boldsymbol{A}}}u,\mathcal{G}_{{\boldsymbol{A}}}v_{T})_{T}-i\sum_{F\in\mathcal{F}_{T}}(\mathcal{G}_{{\boldsymbol{A}}}u\cdot\boldsymbol{n}_{TF},v_{T}-v_{F})_{F}-(\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}I^{k}_{T}u,\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF})_{T}$
		$\displaystyle=(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),\mathcal{G}_{{\boldsymbol{A}}}v_{T})_{T}+\hbox{\pagecolor{gray!20}$\displaystyle(\mathcal{G}_{{\boldsymbol{A}}}u-\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),\mathcal{G}_{{\boldsymbol{A}}}v_{T})_{T}$}-i\sum_{F\in\mathcal{F}_{T}}(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u)\cdot\boldsymbol{n}_{TF},v_{T}-v_{F})_{F}$
		$\displaystyle\hskip 18.49988pt-(\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}I^{k}_{T}u,\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF})_{T}\hbox{\pagecolor{gray!20}$\displaystyle-i\sum_{F\in\mathcal{F}_{T}}((\mathcal{G}_{{\boldsymbol{A}}}u-\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u))\cdot\boldsymbol{n}_{TF},v_{T}-v_{F})_{F}$}$
		$\displaystyle=(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),\mathcal{G}_{{\boldsymbol{A}}}v_{T})_{T}-i\sum_{F\in\mathcal{F}_{T}}(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u)\cdot\boldsymbol{n}_{TF},v_{T}-v_{F})_{F}+\hbox{\pagecolor{gray!20}$\displaystyle\mathfrak{T}_{2}$}$
		$\displaystyle\hskip 18.49988pt-(\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}I^{k}_{T}u,\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF})_{T}$

where we have defined $\displaystyle\mathfrak{T}_{2}$ as the sum of the two terms in grey in the second equality. We notice that, by definition of $\mathrm{G}_{T,{\boldsymbol{A}}}^{k}$ (13) with $\boldsymbol{\tau}=\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u)$ as a test function, we have

(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF})_{T}=(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),\mathcal{G}_{{\boldsymbol{A}}}v_{T})_{T}-i\sum_{F\in\mathcal{F}_{T}}(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u)\cdot\boldsymbol{n}_{TF},v_{T}-v_{F})_{F},

which corresponds exactly to the first two terms, i.e. $M_{T}$ rewrites as

M_{T}=(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF})_{T}-(\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}I^{k}_{T}u,\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF})_{T}+\mathfrak{T}_{2}.

Using (23) with ${\boldsymbol{A}}={\boldsymbol{A}}_{T}$ , we have $\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}I^{k}_{T}u=\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}_{T}}u)$ , therefore $M_{T}$ rewrites as

	$\displaystyle M_{T}$	$\displaystyle=(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF})_{T}-(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}_{T}}u),\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF})_{T}+\mathfrak{T}_{2}$
		$\displaystyle=(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF})_{T}\pm(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),\mathrm{G}_{T,{\boldsymbol{A}}_{T}}^{k}\underline{v}_{TF})-(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}_{T}}u),\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF})_{T}+\mathfrak{T}_{2}$
		$\displaystyle=(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF}-\mathrm{G}_{T,{\boldsymbol{A}}_{T}}^{k}\underline{v}_{TF})_{T}+(\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u)-\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}_{T}}u),\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF})_{T}+\mathfrak{T}_{2}$
		$\displaystyle=:\mathfrak{T}_{1}+\mathfrak{T}_{2}.$

It remains to bound $\mathfrak{T}_{1}$ and $\mathfrak{T}_{2}$ . For the first term, we have

	$\displaystyle\|\mathfrak{T}_{1}\|$	$\displaystyle\leq\|\|\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u)\|\|_{T}\|\|\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF}-\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF}\|\|_{T}+\|\|\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u-\mathcal{G}_{{\boldsymbol{A}}_{T}}u)\|\|_{T}\|\|\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF}\|\|_{T}$
		$\displaystyle\leq\|\|\mathcal{G}_{{\boldsymbol{A}}}u\|\|_{T}Ch_{T}^{k+1}+\|\|({\boldsymbol{A}}_{T}-{\boldsymbol{A}})u\|\|_{T}\|\|\underline{v}_{TF}\|\|_{1,T}$
		$\displaystyle\leq C\|\|\mathcal{G}_{{\boldsymbol{A}}}u\|\|_{T}h_{T}^{k+1}+\|\|{\boldsymbol{A}}-{\boldsymbol{A}}_{T}\|\|_{L^{\infty}}\|\|u\|\|_{T}\|\|\underline{v}_{TF}\|\|_{1,T}$
		$\displaystyle\leq C\|\|\mathcal{G}_{{\boldsymbol{A}}}u\|\|_{T}h_{T}^{k+1}+\widetilde{C}h^{k+1}\|{\boldsymbol{A}}\|_{W^{k+1,\infty(T)}}\|\|u\|\|_{T}\|\|\underline{v}_{TF}\|\|_{1,T}$
		$\displaystyle\lesssim h^{k+1}\|{\boldsymbol{A}}\|_{W^{k+1,\infty}}\|u\|_{H^{k+1}(T)}\|\|\underline{v}_{TF}\|\|_{1,T}$

where we have used Corollary 3.8, the definition of $\mathcal{G}_{{\boldsymbol{A}}}$ together with the linearity of $\pi_{T}^{k}$ , and the boundedness of $\pi_{T}^{k}$ and $\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}$ . For the second term $\mathfrak{T}_{2}$ , we remark that

	$\displaystyle\mathfrak{T}_{2}$	$\displaystyle=(\mathcal{G}_{{\boldsymbol{A}}}u-\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),\mathcal{G}_{{\boldsymbol{A}}}v_{T})_{T}-i\sum_{F\in\mathcal{F}_{T}}((\mathcal{G}_{{\boldsymbol{A}}}u-\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u))\cdot\boldsymbol{n}_{TF},v_{T}-v_{F})_{F}$
		$\displaystyle=(\mathcal{G}_{{\boldsymbol{A}}}u-\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),-i\boldsymbol{\nabla}v_{T}-{\boldsymbol{A}}v_{T})_{T}-i\sum_{F\in\mathcal{F}_{T}}((\mathcal{G}_{{\boldsymbol{A}}}u-\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u))\cdot\boldsymbol{n}_{TF},v_{T}-v_{F})_{F}$
		$\displaystyle=(\mathcal{G}_{{\boldsymbol{A}}}u-\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u),-{\boldsymbol{A}}v_{T})_{T}-i\sum_{F\in\mathcal{F}_{T}}((\mathcal{G}_{{\boldsymbol{A}}}u-\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u))\cdot\boldsymbol{n}_{TF},v_{T}-v_{F})_{F}$

since $\boldsymbol{\nabla}v_{T}\in\mathbb{P}^{k-1}(T)^{d}$ , and $\mathcal{G}_{{\boldsymbol{A}}}u-\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u)$ cancels against any polynomial in $\mathbb{P}^{k}(T)^{d}$ . We can now estimate $\mathfrak{T}_{2}$ as follows:

	$\displaystyle\|\mathfrak{T}_{2}\|$	$\displaystyle\leq\|\|\mathcal{G}_{{\boldsymbol{A}}}u-\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u)\|\|_{T}\|\|{\boldsymbol{A}}\|\|_{L^{\infty}(T)}\|\|v_{T}\|\|_{T}+\sum_{F}\|\|\mathcal{G}_{{\boldsymbol{A}}}u-\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u)\|\|_{F}\|\|v_{T}-v_{F}\|\|_{F}$
		$\displaystyle\leq C\|\|{\boldsymbol{A}}\|\|_{L^{\infty}(T)}h_{T}^{k+1}\|\mathcal{G}_{{\boldsymbol{A}}}u\|_{H^{k+1}(T)}\|\|v_{T}\|\|_{T}+\sum_{F}Ch_{T}^{k+1}\|\mathcal{G}_{{\boldsymbol{A}}}u\|_{H^{k+1}(T)}h^{-\frac{1}{2}}\|\|v_{T}-v_{F}\|\|_{F}$
		$\displaystyle\leq C\|\|{\boldsymbol{A}}\|\|_{L^{\infty}(T)}h_{T}^{k+1}\|\mathcal{G}_{{\boldsymbol{A}}}u\|_{H^{k+1}(T)}\|\|v_{T}\|\|_{T}+N_{\partial}Ch_{T}^{k+1}\|\mathcal{G}_{{\boldsymbol{A}}}u\|_{H^{k+1}(T)}\|\|\underline{v}_{TF}\|\|_{1,T}.$

where we have used (30) for the first term and (31) for the trace terms, and the definition of the norms (8). Summing $\mathfrak{T}_{1}$ and $\mathfrak{T}_{2}$ over all elements $T\in\mathcal{T}_{h}$ , using

|\mathcal{G}_{{\boldsymbol{A}}}u|_{H^{k+1}(T)}\leq C\left(|u|_{H^{k+2}(T)}+\|\mathbf{A}\|_{W^{k+1,\infty}(T)}\|u\|_{H^{k+1}(T)}\right),

and applying the discrete Cauchy-Schwarz and Poincaré inequalities (Lemma A.5) yields the global optimal bound for $|\mathfrak{T}_{grad}|$ , which concludes the proof for the energy norm. The $L^{2}$ -error bound relies on the classic Aubin-Nitsche duality argument.

4 Numerical Results

Refer to caption — Figure 1: Meshes used in the numerical experiments in Section 4 – From first to third: the 8x8x8 cubes, hex and the Voronoi meshes; the fourth is a 3D simplicial mesh utilized for the Aharonov-Bohm experimental test case ( $\approx$ 88k cells).

In this section, we validate our results through two distinct test cases. First, we address the eigenvalue problem in Eq. Eq. 1 and confirm that the discrepancies between the eigenvalues are effectively controlled by the upper bound prediction in Theorem 3.4. Second, we investigate the unsteady Schrödinger equation in Eq. (4) to replicate the Aharonov-Bohm experiment, wherein the presence of a magnetic field induces a phase shift in the propagating photons. Both simulations were conducted using the HArDCore3D library¹¹1https://github.com/jdroniou/HArDCore3D-release, a three-dimensional implementation of the Hybrid High-Order (HHO) schemes [DiPietro2020], in conjunction with the Spectra library [spectralib] for eigenvalue computations.

4.1 The Eigenvalue Problem: Fock-Darwin Spectrum and Gauge Invariance

To assess the accuracy of our HHO scheme for the eigenvalue problem (1) and numerically validate the discrete gauge covariance established in Theorem 3.4, we consider a two-dimensional quantum harmonic oscillator subjected to a uniform perpendicular magnetic field $\boldsymbol{B}=B\boldsymbol{e}_{z}$ .

We operate on a sufficiently large truncated domain $\Omega=[-L,L]^{2}$ with homogeneous Dirichlet boundary conditions. Following our initial non-dimensionalization ( $\hbar=m=q=1$ ), the scalar confinement potential is given by:

(24)

V(x,y)=\frac{1}{2}\omega_{0}^{2}(x^{2}+y^{2}),

where $\omega_{0}>0$ is the oscillator frequency. A constant magnetic field $\boldsymbol{B}$ can be represented by infinitely many vector potentials ${\boldsymbol{A}}$ related by gauge transformations (2.4). To highlight the exact discrete gauge invariance of our scheme under discrete gauge transformations, we solve the generalized eigenvalue problem using three distinct gauges, the symmetric gauge ${\boldsymbol{A}}_{\text{sym}}$ , the Landau ${\boldsymbol{A}}_{\text{Lan}}$ and a manufactured gauge ${\boldsymbol{A}}_{\text{smooth}}$ defined as

{\boldsymbol{A}}_{\text{sym}}=\frac{B}{2}\begin{pmatrix}-y\\ x\end{pmatrix},\qquad{\boldsymbol{A}}_{\text{Lan}}=\begin{pmatrix}-By\\ 0\end{pmatrix},\quad{\boldsymbol{A}}_{\text{smooth}}=\begin{pmatrix}-\tfrac{1}{2}By+0.1\\ \tfrac{1}{2}Bx+0.1\end{pmatrix}.

Note that these quantities are related by ${\boldsymbol{A}}_{\text{Lan}}={\boldsymbol{A}}_{\text{sym}}+\boldsymbol{\nabla}\left(-\frac{B}{2}xy\right)$ and ${\boldsymbol{A}}_{\text{smooth}}={\boldsymbol{A}}_{\text{sym}}+\boldsymbol{\nabla}\left(0.1(x+y)\right)$ . The exact analytical eigenvalues of (1) in three dimensions are derived from the 2D Fock-Darwin spectrum [Governale1998, Eq. (23)], extended by a 1D particle-in-a-box spectrum along the $z$ -axis due to the homogeneous Dirichlet boundary conditions imposed on $z\in[-L,L]$ . They are independent of the gauge ${\boldsymbol{A}}$ and are given, for any quantum numbers $n\in\mathbb{N}$ , $m\in\mathbb{Z}$ , and $n_{z}\in\mathbb{N}^{*}$ , by:

(25)

E_{n,m,n_{z}}=\sqrt{B^{2}+2\omega_{0}^{2}}(2n+|m|+1)-mB+\frac{n_{z}^{2}\pi^{2}}{(2L)^{2}}.

The fundamental energy level, known as the Fock-Darwin ground state, corresponds to $n=m=0$ and $n_{z}=1$ , yielding:

(26)

E_{0,0,1}=\sqrt{B^{2}+2\omega_{0}^{2}}+\frac{\pi^{2}}{(2L)^{2}}.

In what follows, we assume that $w_{0}=1$ and $B=1$ . The discrete eigenvalue problem reads : find $u_{h}\in U_{h}^{k}$ and $\lambda_{h}\in\mathbb{R}$ such that $a_{h}(\underline{u}_{h},\underline{v}_{h};{\boldsymbol{A}}_{T})=\lambda_{h}(\underline{u}_{h},\underline{v}_{h})_{L^{2}(\Omega)}$ . In Figure 2, we report the maximum deviation absolute error

\mathtt{dev\_error}(\mathrm{g}):=\max_{0\leq j\leq N-1}|\lambda_{h,j}^{(\text{sym})}-\lambda_{h,j}^{(\text{gauge})}|,

with $\mathrm{g}\in\{\mathrm{smooth},\mathrm{Landau}\}$ , for the first $N=5$ eigenvalues across polynomial degrees from $k=0$ to $k=3$ , on progressively refined meshes. We observe a superconvergence of one order for $k=0$ (in blue), but an optimal decay rate for $k=2$ (in black), suggesting that the prediction rate $\mathcal{O}(h^{k+1})$ is optimal and confirming that the discrete spectrum behaves covariantly consistent with the discretization error.

Figure 2: Gauge transformation errors : (Left)

\mathtt{dev\_error}(\text{smooth})

; (Right)

\mathtt{dev\_error}(\text{Landau})

, on three different meshes : cubes (Top) and hex (Bottom).

Figure 3: Convergence decay of the relative error (27) on the cubes (left) and Voronoi (right) meshes.

In our context, the fundamental energy (26) is given as $E_{0,0,1}=\sqrt{3}+\frac{\pi^{2}}{64}\approx 1.8862633763358985.$ In Figure 3, we plot the relative errors of the fundamental energy level

(27)

\mathtt{rel\_err}:=|\lambda_{h,0}-E_{0,0,1}|/E_{0,0,1},

against the mesh size $h$ on logarithmic scales, for two families of meshes : cubes and Voronoi meshes, see Fig. 1. We observe that the discrete eigenvalues converge towards the exact Fock-Darwin levels, the convergence is nearly at the optimal rate of $\mathcal{O}(h^{2k+2})$ for the cubes, which confirm the high-order accuracy of the proposed discrete covariant gradient for spectral computations. On the other hand, suboptimal error rates are observed for the Voronoi mesh family, which might be due to boundary conditions together with higher value of $h$ .

4.2 The Aharonov-Bohm Effect

The Aharonov-Bohm (AB) effect is a quantum mechanical phenomenon in which a charged particle is affected by the electromagnetic vector potential ${\boldsymbol{A}}$ even in regions where the magnetic field $\boldsymbol{B}=\boldsymbol{\nabla}{\times}{\boldsymbol{A}}$ vanishes identically [Aharonov1959].

Setup

We simulate a quantum particle traveling through a 3D rectangular domain $\Omega=[-5,5]\times[-2,2]\times[-2,2]$ from which a cylindrical solenoid of radius $r_{s}=0.5$ is excluded, see Figure 4 for a 2D representation and Figure 1 for a 3D representation with a simplicial mesh. The particle state is encoded in the quantum state $\underline{\psi}_{h}$ , solution of the time-dependant discrete magnetic Schrödinger equation expressed with the discrete covariant gradient: find $\underline{\psi}_{h}\in U_{h}^{k}$ such that

	$\displaystyle(\partial_{t}\psi_{h},\phi_{h})_{L^{2}(\Omega)}+a_{h}(\underline{\psi}_{h},\underline{\phi}_{h})$	$\displaystyle=\ell(\underline{\phi}_{h}),\qquad\text{ in }\Omega,$
	$\displaystyle\underline{\psi}_{h}(t=0)$	$\displaystyle=I^{k}_{h}\psi_{0},\qquad\text{ in }\Gamma.$

We use a Crank-Nicholson time scheme with homogeneous Dirichlet boundary conditions on all boundaries (including outer walls and solenoid surface). The magnetic vector potential outside the solenoid is the classical irrotational field,

(28)

{\boldsymbol{A}}(x,y,z)\mathrel{\mathop{:}}=\frac{\Phi}{2\pi}\frac{1}{x^{2}+y^{2}}\begin{pmatrix}-y\\ x\\ 0\end{pmatrix},\qquad(x,y,z)\in\Omega,

where $\Phi$ is the total magnetic flux confined inside the solenoid. Note that $\boldsymbol{\nabla}{\times}{\boldsymbol{A}}=\mathbf{0}$ everywhere in $\Omega$ by construction, so the particle experiences no local magnetic force in the domain; the entire effect is encoded in the topology of $\Omega$ and the circulation of ${\boldsymbol{A}}$ around the solenoid.

The initial wavepacket is a Gaussian centered at $x_{0}=-3.5$ propagating in the $+x$ direction,

(29)

\psi_{0}(x,y,z)=\exp\!\left(-\frac{(x-x_{0})^{2}+y^{2}+z^{2}}{2\sigma^{2}}\right)e^{ik_{0}x},\qquad\sigma=0.8,\quad k_{0}=3.

As the wavepacket encounters the solenoid, it naturally splits into two partial waves propagating along both sides (top and bottom, see Fig. 4), which then recombine downstream to form an interference pattern.

Theoretical Interference Pattern.

The presence of the vector potential ${\boldsymbol{A}}\neq\boldsymbol{0}$ induces a relative phase shift $\Delta\phi$ between the top and bottom partial waves. This shift is given by the line integral of ${\boldsymbol{A}}$ along the closed loop $\Gamma$ formed by the two paths. By Stokes’ theorem, this is exactly the enclosed magnetic flux:

\Delta\phi=\oint_{\gamma_{\mathrm{top}}}{\boldsymbol{A}}\cdot d\boldsymbol{l}-\oint_{\gamma_{\mathrm{bottom}}}{\boldsymbol{A}}\cdot d\boldsymbol{l}=\oint_{\Gamma}{\boldsymbol{A}}\cdot d\boldsymbol{l}=\iint_{\Omega_{S}}\boldsymbol{\nabla}{\times}{\boldsymbol{A}}\cdot d\boldsymbol{S}=\iint_{\Omega_{S}}\boldsymbol{B}\cdot d\boldsymbol{S}=\Phi.

Let us analyze the expected intensity $I(y)$ on a vertical detection screen downstream, at $x=3.0$ . At the exact center of the screen ( $y=0$ ), the geometric paths are perfectly symmetric, see Fig. 4. Let $\psi_{s}$ denote the local amplitude acquired via each path in the absence of a magnetic field. When turning on the vector potential, the total wavefunction $\psi_{end}$ at the center is the coherent superposition of the two paths endowed with their respective AB phases: $\psi_{end}(0)=\psi_{s}e^{i\int_{\gamma_{\text{top}}}{\boldsymbol{A}}\cdot d\boldsymbol{l}}+\psi_{s}e^{i\int_{\gamma_{\text{bottom}}}{\boldsymbol{A}}\cdot d\boldsymbol{l}}$ . Factoring out the global phase of the bottom path, which vanishes when taking the squared modulus, the resulting local intensity on the screen is given by $I(0)=|\psi_{s}e^{i\Phi}+\psi_{s}|^{2}$ . We compare two cases: When $\Phi=0$ , i.e., no flux. The phase shift is $\Delta\phi=0$ , leading to $I(0)=4|\psi_{s}|^{2}$ , a constructive interference, generating a global maximum (a central peak) at $y=0$ . When $\Phi=\pi$ , the phase shift is $\Delta\phi=\pi$ , leading to $I(0)=0$ , a destructive interference. The probability of finding the particle at the exact center is strictly zero, and the displaced probability creates two symmetric peaks bounding this central node.

Furthermore, the local amplitude $\psi_{s}$ arriving at the screen is expected to be drastically smaller than the initial amplitude ( $|\psi_{0}|=1$ ). This amplitude drop is physically consistent and stems from three combined effects: the natural quantum dispersion of the Gaussian wavepacket over time, the strong backward scattering (reflection) of the wave upon hitting the impenetrable Dirichlet cylinder, and the 3D geometric spreading of the wave in the domain.

Results.

Figure 6 displays the computed probability density $|\psi_{h}|^{2}(x_{\mathrm{screen}},0,0)$ at time $t=1.49$ . At the bottom, we display the 1D profiles extracted along the detection line, which match the theoretical predictions. Indeed, for the zero-flux case ( $\Phi=0$ ), the profile exhibits an isolated central peak at $y=0$ , confirming the constructive interference. Conversely, for the case $\Phi=\pi$ , the central peak completely vanishes and is replaced by a deep central minimum bounded by two symmetric peaks, characterizing the $\pi$ topological phase shift. Additionally, the maximum intensity on the screen correctly reflects the expected amplitude drop, peaking around $\approx 0.45$ .

Remark 4.1.

We also tested the cheaper variant

\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{u}_{TF}:=-i\boldsymbol{\nabla}p^{k+1}_{T}\underline{u}_{TF}-{\boldsymbol{A}}_{T}\,p^{k+1}_{T}\underline{u}_{TF},

which reuses the already computed potential reconstruction $p^{k+1}_{T}$ and avoids solving the local system (13). This variant produces similar results for the AB experiment, with similar eigenvalue errors and gauge deviations on all meshes and degrees $k=0,\dots,2$ , at the only theoretical cost of a suboptimal $O(h)$ energy-norm bound in Theorem 3.14.

5 Conclusion and Perspectives

We have proposed a minimal and natural extension of Hybrid High-Order methods to the Schrödinger equation with vector potential. The scheme relies solely on the existing potential reconstruction operator, allowing for concise derivation of discrete gauge invariance and coercivity. Future work includes extending to the time-dependent equation, the Pauli model for spin- $\tfrac{1}{2}$ particles, the Dirac equation, as well as a posteriori error analysis and adaptivity on complex geometries.

Acknowledgments

The author acknowledges the use of Gemini 3.1 Pro GenAI to: (i) polish the English translation, (ii) verify proofs of Section 3, and (iii) generate C++ boilerplate and Python post-processing codes for the numerical experiments.

	$\displaystyle\mathfrak{T}_{2}$	$\displaystyle\leq\sum_{F}\\|u_{F}-u_{T}\\|_{F}\\|\boldsymbol{\nabla}u\cdot\boldsymbol{n}_{TF}\\|_{F}\leq\left(\sum_{F}h_{F}^{-1}\\|u_{F}-u_{T}\\|_{F}^{2}\right)^{1/2}\!\left(\sum_{F}h_{F}\\|\boldsymbol{\nabla}u_{T}\\|_{F}^{2}\right)^{1/2}$
		$\displaystyle\leq s_{T}(\underline{u}_{TF},\underline{u}_{TF})^{1/2}\left(\sum_{F}h_{F}^{1}h_{F}^{-1}C_{tr}^{2}\\|\boldsymbol{\nabla}u_{T}\\|_{T}^{2}\right)^{1/2}\leq s_{T}(\underline{u}_{TF},\underline{u}_{TF})^{1/2}\,C_{tr}\\|\boldsymbol{\nabla}u_{T}\\|_{T}N_{\partial}^{1/2}$
		$\displaystyle\leq C_{tr}N_{\partial}^{1/2}\,s_{T}(\underline{u}_{TF},\underline{u}_{TF})^{1/2}\\|\boldsymbol{\nabla}u_{T}\\|_{T}.$

	$\displaystyle\|\mathfrak{T}_{\boldsymbol{E}}\|$	$\displaystyle\leq\\|{\boldsymbol{A}}\\|_{L^{\infty}(T)}\\|v_{T}\\|_{T}\\|\boldsymbol{E}\\|_{T}+\sum_{F\in\mathcal{F}_{T}}\\|v_{T}-v_{F}\\|_{F}\\|\boldsymbol{E}\\|_{F}$
		$\displaystyle\leq C\\|{\boldsymbol{A}}\\|_{L^{\infty}(T)}\\|v_{T}\\|_{T}\left(h_{T}^{k+1}\|e^{-i\chi}\boldsymbol{\tau}\|_{H^{k+1}(T)}\right)+\sum_{F\in\mathcal{F}_{T}}\|v_{h}\|_{1,T}h_{F}^{\frac{1}{2}}\left(h_{T}^{k+\frac{1}{2}}\|e^{-i\chi}\boldsymbol{\tau}\|_{H^{k+1}(T)}\right)$
		$\displaystyle\leq Ch_{T}^{k+1}\|v_{h}\|_{1,T}\\|e^{-i\chi}\boldsymbol{\tau}\\|_{H^{k+1}(T)},$

	$\displaystyle\|\mathfrak{T}_{pot}\|$	$\displaystyle=\left\|\sum_{T\in\mathcal{T}_{h}}(Vu-V\pi_{T}^{k}u,v_{T})_{T}\right\|\leq\sum_{T\in\mathcal{T}_{h}}\|\|V\|\|_{L^{\infty}(\Omega)}\\|u-\pi_{T}^{k}u\\|_{T}\\|v_{T}\\|_{T}$
		$\displaystyle\leq\|\|V\|\|_{L^{\infty}(\Omega)}\sum_{T}Ch_{T}^{k+1}\|u\|_{H^{k+2}(T)}\\|v_{T}\\|_{T}$
		$\displaystyle=\|\|V\|\|_{L^{\infty}(\Omega)}Ch^{k+1}\|u\|_{H^{k+2}(\Omega)}\\|v_{h}\\|_{L^{2}(\Omega)}$
		$\displaystyle\lesssim\|\|V\|\|_{L^{\infty}(\Omega)}h^{k+1}\|u\|_{H^{k+2}(\Omega)}\\|\underline{v}_{h}\\|_{1,h}.$

	$\displaystyle\|\mathfrak{T}_{stab}\|$	$\displaystyle=\left\|\sum_{T\in\mathcal{T}_{h}}s_{T}(I^{k}_{T}u,\underline{v}_{TF})\right\|=\left\|\sum_{T\in\mathcal{T}_{h}}\sum_{F\in\mathcal{F}_{T}}h_{F}^{-1}(\delta_{TF}^{k}(u)-\delta_{T}^{k}(u),v_{F}-v_{T})_{F}\right\|$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}\sum_{F\in\mathcal{F}_{T}}h_{F}^{-\tfrac{1}{2}}\\|\delta_{TF}^{k}(u)-\delta_{T}^{k}(u)\\|_{F}h_{F}^{-\tfrac{1}{2}}\\|v_{F}-v_{T}\\|_{F}$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}\left(\sum_{F\in\mathcal{F}_{T}}h_{F}^{-1}\\|\delta_{TF}^{k}(u)-\delta_{T}^{k}(u)\\|_{F}^{2}\right)^{\tfrac{1}{2}}\left(\sum_{F\in\mathcal{F}_{T}}h_{F}^{-1}\\|v_{F}-v_{T}\\|_{F}^{2}\right)^{\tfrac{1}{2}}$
		$\displaystyle\leq\sum_{T\in\mathcal{T}_{h}}s_{T}(I^{k}_{T}u,I^{k}_{T}u)^{\tfrac{1}{2}}s_{T}(\underline{v}_{TF},\underline{v}_{TF})^{\tfrac{1}{2}}$
		$\displaystyle\leq Ch^{k+1}\sum_{T\in\mathcal{T}_{h}}\\|u\\|_{H^{k+2}(T)}\\|\underline{v}_{TF}\\|_{1,T}\lesssim h^{k+1}\|u\|_{H^{k+2}(\Omega)}\\|\underline{v}_{h}\\|_{1,h}$

	$\displaystyle\|\mathfrak{T}_{1}\|$	$\displaystyle\leq\|\|\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u)\|\|_{T}\|\|\mathrm{G}_{T,{\boldsymbol{A}}}^{k}\underline{v}_{TF}-\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF}\|\|_{T}+\|\|\pi_{T}^{k}(\mathcal{G}_{{\boldsymbol{A}}}u-\mathcal{G}_{{\boldsymbol{A}}_{T}}u)\|\|_{T}\|\|\mathrm{G}_{T,{\boldsymbol{A}_{T}}}^{k}\underline{v}_{TF}\|\|_{T}$
		$\displaystyle\leq\|\|\mathcal{G}_{{\boldsymbol{A}}}u\|\|_{T}Ch_{T}^{k+1}+\|\|({\boldsymbol{A}}_{T}-{\boldsymbol{A}})u\|\|_{T}\|\|\underline{v}_{TF}\|\|_{1,T}$
		$\displaystyle\leq C\|\|\mathcal{G}_{{\boldsymbol{A}}}u\|\|_{T}h_{T}^{k+1}+\|\|{\boldsymbol{A}}-{\boldsymbol{A}}_{T}\|\|_{L^{\infty}}\|\|u\|\|_{T}\|\|\underline{v}_{TF}\|\|_{1,T}$
		$\displaystyle\leq C\|\|\mathcal{G}_{{\boldsymbol{A}}}u\|\|_{T}h_{T}^{k+1}+\widetilde{C}h^{k+1}\|{\boldsymbol{A}}\|_{W^{k+1,\infty(T)}}\|\|u\|\|_{T}\|\|\underline{v}_{TF}\|\|_{1,T}$
		$\displaystyle\lesssim h^{k+1}\|{\boldsymbol{A}}\|_{W^{k+1,\infty}}\|u\|_{H^{k+1}(T)}\|\|\underline{v}_{TF}\|\|_{1,T}$

Asymptotic gauge-invariant Hybrid High-Order method for magnetic Schrödinger equations

Abstract

keywords: