I've recently started writing a library of numerical solvers for elliptic partial differential equations, with particular focus on the Poisson equation. If one considers typical Poisson equation in two dimensions: $$\nabla^2\varphi(x, y)=f(x, y)$$ with corresponding domain and boundary condition definition, there are numerous ways of how to solve this numerically. Its discretization in terms of finite differences becomes rather straight-forward, as it yields: $$\frac{\varphi_{i+1,j}-2\varphi{i,j}+\varphi_{i-1,j}}{\Delta x^2}+\frac{\varphi_{i,j+1}-2\varphi{i,j}+\varphi_{i,j-1}}{\Delta y^2}=f_{i,j}$$
However, there are also several important physical problems where the Laplacian represents the kinetic term in an overall Hamiltonian, and it has the form: $$(\nabla-i\textbf{Q})^2\varphi=f(x, y)$$ where $\textbf{Q}$ is a vector potential. Note here that $f(x, y)$ is just short-hand notation for all the remaining terms in the equation, which are not part of the Laplacian/kinetic term. These kinds of examples can be seen with Schrodinger equation in magnetic field, or in Gross-Pitaevski theory, or in Ginzburg-Landau theory of phase-transitions applied to quasi-classical treatment of superconductivity.
Now, if one wants to solve the latter equation, for instance, numerically via finite differences, how should one proceed with discretizing the kinetic term $(\nabla-i\textbf{Q})^2\varphi$?
My first instinct would be to just discretize directly, without any fancy methodology, as with the typical Poisson/Laplace equation. However, I might not be in the right here as I've recently found in the literature that people use the so-called Wilson's formulation of lattice gauge invariance (see this and this). In particular Eq. (6) of the second paper is where the so-called link variables are introduced. Now, Wilson's original paper was about quarks and the notation within is beyond my knowledge and experience, plus his motivation is highly specific. On the other hand, based on the other literature, it seems that his approach can be generalized, to a certain extent, which also eludes my understanding. For instance, similar approach was implemented in this paper, where authors derived the method to solve Ginzburg-Landau equations for superconductor by using the link variables (Eq. 11) to discretize the time-dependent Ginzburg-Landau equation (Eq. 12). The Eq. 12 if given in steady state and if linearized transforms into a Poisson- (Laplace-) like equation with vector potential.
So, to reiterate, my question is then how should one proceed with discretizing the term $(\nabla-i\textbf{Q})^2\varphi$ in finite difference framework? Does one really need to consider physical properties of the system and based on those approach with/without link-variables? Or are link-variables purely mathematically a formal and mandatory way to deal with the discretization of $(\nabla-i\textbf{Q})^2\varphi$? If so, what is then the mathematical intuition behind the link variables?