The infamous numerical sign problem poses a fundamental obstacle to particle- based stochastic Wigner simulations in high-dimensional phase space. Although the existing particle annihilation (PA) via uniform mesh significantly alleviates the sign problem when dimensionality D <= 4, the mesh size grows dramatically when D >= 6 due to the curse of dimensionality and consequently makes the annihilation very inefficient. In this paper, we propose an adaptive PA algorithm, termed sequential-clustering particle annihilation via discrepancy estimation (SPADE), to overcome the sign problem. SPADE follows a divide-and-conquer strategy: adaptive clustering of particles via controlling their number-theoretic discrepancies and independent random matching in each cluster. The target is to alleviate the oversampling problem induced by the overpartitioning of phase space and to capture the nonclassicality of the Wigner function simultaneously. Combining SPADE with the variance reduction technique based on the stationary phase approximation, we attempt to simulate the proton-electron couplings in six- and 12-dimensional phase space. A thorough performance benchmark of SPADE is provided with the reference solutions in six-dimensional phase space produced by a characteristic-spectral-mixed scheme under a 733*803 uniform grid, which fully explores the limit of grid-based deterministic Wigner solvers.

A novel adaptive spectral method has been recently developed to numerically solve partial differential equations (PDEs) in unbounded domains. To achieve accuracy and improve efficiency, the method relies on the dynamic adjustment of three key tunable parameters: the scaling factor, a displacement of the basis functions, and the spectral expansion order. In this paper, we perform the first numerical analysis of the adaptive spectral method using generalized Hermite functions in both one- and multi-dimensional problems. Our analysis reveals why adaptive spectral methods work well when a “frequency indicator” of the numerical solution is controlled. We then investigate how the implementation of the adaptive spectral methods affects numerical results, thereby providing guidelines for the proper tuning of parameters. Finally, we further improve performance by extending the adaptive methods to allow bidirectional basis function translation, and the prospect of carrying out similar numerical analysis to solving PDEs arising from realistic difficult-to- solve unbounded models with adaptive spectral methods is also briefly discussed.

Numerical resolution for 6-D Wigner dynamics under the Coulomb potential is faced with the combined challenges of high dimensionality, nonlocality, oscillation, and singularity. In particular, the extremely huge memory storage of 6-D grids hinders the usage of all existing deterministic numerical schemes, which is well known as the curse of dimensionality. To surmount these difficulties, we propose a massively parallel solver, termed the characteristic-spectral-mixed (CHASM) scheme, by fully exploiting two distinct features of the Wigner equation: locality of spatial advection and nonlocality of quantum interaction. Our scheme utilizes the local cubic B-spline basis to interpolate the local spatial advection. The key is to use a perfectly matched boundary condition to give a closure of spline coefficients, so that distributed pieces can recover the global one as accurately as possible owing to the rapid decay of wavelet basis in the dual space, and communication costs are significantly reduced. To resolve the nonlocal pseudodifferential operator with a singular symbol, CHASM further adopts the truncated kernel method to attain a highly efficient approximation. Several typical experiments including the quantum harmonic oscillator and the 1s state of hydrogen demonstrate the accuracy and efficiency of CHASM. Nonequilibrium electron-proton couplings are also clearly displayed and illuminate the uncertainty principle and quantum tunneling in phase space. Finally, the scalability of CHASM up to 16000 cores is presented.

An approach to generalize any kind of collinear functional in density functional theory to noncollinear functionals is proposed. This approach satisfies the correct collinear limit for any kind of functional, guaranteeing that the exact collinear functional after generalization is still exact for collinear spins. Besides, it has well-defined and numerically stable functional derivatives, a desired feature for noncollinear and spin-flip time-dependent density functional theory. Furthermore, it provides local torque, hinting at its applications in spin dynamics.

Nonlocal modeling has drawn more and more attention and becomes steadily more powerful in scientific computing. In this paper, we demonstrate the superiority of a first-principle nonlocal model—Wigner function—in treating singular potentials which are often used to model the interaction between point charges in quantum science. The nonlocal nature of the Wigner equation is fully exploited to convert the singular potential into the Wigner kernel with weak or even no singularity, and thus highly accurate numerical approximations are achievable, which are hardly designed when the singular potential is taken into account in the local Schrödinger equation. The Dirac delta function, the logarithmic, and the inverse power potentials are considered. Numerically converged Wigner functions under all these singular potentials are obtained with a fourth-order accurate operator splitting spectral method, and display many interesting quantum behaviors as well.

An accurate description of 2-D quantum transport in a double-gate metal oxide semiconductor field effect transistor (dgMOSFET) requires a high-resolution solver for a coupled system of the 4-D Wigner equation and 2-D Poisson equation. In this paper, we propose an operator-splitting spectral method to evolve such Wigner–Poisson (WP) system in 4-D phase space with high accuracy. After the operator splitting of the Wigner equation, the resulting two sub-equations can be solved analytically with spectral approximation in phase space. Meanwhile, we adopt a Chebyshev spectral method to solve the Poisson equation. Spectral convergence in phase space and a fourth-order accuracy in time are both numerically verified. Finally, we apply the proposed solver to the simulation of a dgMOSFET, develop the steady states via long-time simulations and obtain numerically converged current–voltage (I–V) curves.

When using Laguerre and Hermite spectral methods to numerically solve PDEs in unbounded domains, the number of collocation points assigned inside the region of interest is often insufficient, particularly when the region is expanded or translated in order to safely capture the unknown solution. Simply increasing the number of collocation points cannot ensure a fast convergence to spectral accuracy. In this paper, we propose a scaling technique and a moving technique to adaptively cluster enough collocation points in a region of interest in order to achieve fast spectral convergence. Our scaling algorithm employs an indicator in the frequency domain that both is used to determine when scaling is needed and informs the tuning of a scaling factor to redistribute collocation points in order to adapt to the diffusive behavior of the solution. Our moving technique adopts an exterior-error indicator and moves the collocation points to capture the translation. Both frequency and exterior-error indicators are defined using only the numerical solutions. We apply our methods to a number of different models, including diffusive and moving Fermi--Dirac distributions and nonlinear Dirac solitary waves, and demonstrate recovery of spectral convergence for time-dependent simulations. A performance comparison in solving a linear parabolic problem shows that our frequency scaling algorithm outperforms the existing scaling approaches. We also show our frequency scaling technique is able to track the blowup of average cell sizes in a model for cell proliferation. In addition to the Laguerre and Hermite basis functions with exponential decay at infinity, we also successfully apply the frequency-dependent scaling technique into rational basis functions with algebraic decay at infinity.

A set-pair Lovász extension is established to construct equivalent continuous optimization problems for graph k-cut problems.

Existing critical point theories including metric and topological critical point theories are difficult to be applied directly to some concrete problems in particular polyhedral settings, because the notions of critical sets could be either very vague or too large. To overcome these difficulties, we develop critical point theory for nonsmooth but Lipschitzian functions defined on convex polyhedrons. This yields natural extensions of classical results in critical point theory, such as the Liusternik-Schnirelmann multiplicity theorem. More importantly, eigenvectors for some eigenvalue problems involving graph 1-Laplacian coincide with critical points of the corresponding functions on polytopes, which indicates that the critical point theory proposed in the present paper can be applied to study the nonlinear spectral graph theory.

When using spectral methods, a consistent method for tuning the expansion order is often required, especially for time-dependent problems in which oscillations emerge in the solution. In this paper, we propose a frequency-dependent p-adaptive technique that adaptively adjusts the expansion order based on a frequency indicator. Using this p-adaptive technique, combined with recently proposed scaling and moving techniques, we are able to devise an adaptive spectral method in unbounded domains that can capture and handle diffusion, advection, and oscillations. As an application, we use this adaptive spectral method to numerically solve Schrödinger’s equation in an unbounded domain and successfully capture the solution’s oscillatory behavior at infinity.

The stochastic solutions to the Wigner equation, which explain the nonlocal oscillatory integral operator $\Theta_V$ with an anti-symmetric kernel as the generator of two branches of jump processes,  are analyzed. All existing branching random walk solutions are formulated based on the Hahn-Jordan decomposition $\Theta_V=\Theta^+_V-\Theta^-_V$, i.e., treating $\Theta_V$ as the difference of two positive operators $\Theta^\pm_V$, each of which characterizes the transition of states for one branch of particles. Despite the fact that the first moments of such models solve the Wigner equation, we prove that the bounds of corresponding variances grow exponentially in time with the rate depending on the upper bound of $\Theta^\pm_V$, instead of $\Theta_V$. In other words, the decay of high-frequency components is totally ignored, resulting in a severe numerical sign problem. To fully utilize such decay property, we have recourse to the stationary phase approximation for $\Theta_V$, which captures essential contributions from the stationary phase points as well as the near-cancelation of positive and negative weights. The resulting branching random walk solutions are then proved to asymptotically solve the Wigner equation, but gain a substantial reduction in variances, thereby ameliorating the sign problem. Numerical experiments in 4-D phase space validate our theoretical findings. 

A branching random walk algorithm for many-body Wigner equations and its numerical applications for quantum dynamics in phase space are proposed and analyzed in this paper. Using an auxiliary function, the truncated Wigner equation and its adjoint form are cast into integral formulations, which can be then reformulated into renewal-type equations with probabilistic interpretations. We prove that the first moment of a branching random walk is the solution for the adjoint equation. With the help of the additional degree of freedom offered by the auxiliary function, we are able to produce a weighted-particle implementation of the branching random walk. In contrast to existing signed-particle implementations, this weighted-particle one shows a key capacity of variance reduction by increasing the constant auxiliary function and has no time discretization errors. Several canonical numerical experiments on the 2D Gaussian barrier scattering and a 4D Helium-like system validate our theoretical findings, and demonstrate the accuracy, the efficiency, and thus the computability of the proposed weighted-particle Wigner branching random walk algorithm.

We propose a high order numerical method for computing time dependent 4-D Wigner equation with unbounded potential and study a canonical quantum double-slit interference problem. To address the difficulties of 4-D phase space computations and higher derivatives from the Moyal expansion of nonlocal pseudo-differential operator for unbounded potentials, an operator splitting technique is adopted to decompose the 4-D Wigner equation into two sub-equations, which can be computed analytically or numerically with high efficiency. The first sub-equation contains only linear convection term in $(\bm x, t)$-space and can be solved with an advective method, while the second involves the pseudo-differential term and can be approximated by a plane wave expansion in $\bm k$-space. By exploiting properties of Fourier transformation,  the expansion coefficients for the second sub-equation have explicit forms and the resulting scheme is shown to be unconditionally stable for any higher derivatives of the Moyal expansion, ensuring the feasibility of the 4-D Wigner numerical simulations for quantum double-slit interferences. Numerical experiments demonstrate the spectral convergence in $(\bm x, \bm k)$-space and provide highly accurate information on the number, position, and intensity of the interference fringes for different types of slits, quantum particle masses, and initial states (pure and mixed).

The externally driven nonlinear Dirac (NLD) equation with scalar-scalar  self-interaction studied in [J. Phys. A: Math. Theor. 49, 065402 (2016)] is revisited. By using a variational method and an ansatz with five collective coordinates, the dynamics of the NLD solitons  is well described. It is shown that this new ansatz possesses certain advantages, namely the canonical momentum agrees with the field momentum, the energy associated to the collective coordinate equations agrees with the energy of the NLD soliton, whereas the ansatz with either three or four collective coordinates does not. Thus the study of the  whole phase space of the system is enhanced. It is also shown that this approach is equivalent to the method of moments: the time variation of the charge,  the momentum, the energy, and the first moment  of the charge. The advantages of the new ansatz are illustrated by means of numerical simulations. 

Unbounded potentials are always utilized to strictly confine quantum dynamics and generate bound or stationary states due to the existence of quantum tunneling. However, the existed accurate Wigner solvers are often designed for either localized potentials or those of the polynomial type. This paper attempts to solve the time-dependent Wigner equation in the presence of a general class of unbounded potentials by exploiting two equivalent forms of the pseudo-differential operator: integral form and series form (i.e., the Moyal expansion). The unbounded parts at infinities are approximated or modeled by polynomials and then a remaining localized potential dominates the central area. The fact that the Moyal expansion reduces to a finite series for polynomial potentials is fully utilized. In order to accurately resolve both the pseudo-differential operator and the linear differential operator,  a spectral collocation scheme for the phase space and an explicit fourth-order Runge-Kutta time discretization are adopted. We are able to prove that the resulting full discrete spectral scheme conserves both mass and energy. Several typical quantum systems are simulated with a high accuracy and reliable estimation of macroscopically measurable quantities is thus obtained.  

To implement the Wigner branching random walk, the particle carrying a signed weight, either $-1$ or $+1$, is more friendly to data storage and arithmetic manipulations than that taking a real-valued weight continuously from $-1$ to $+1$. The former is called a signed particle and the latter a weighted particle. In this paper, we propose two efficient strategies to realize the signed-particle implementation. One is to interpret the multiplicative functional as the probability to generate pairs of particles instead of the incremental weight, and the other is to utilize a bootstrap filter to adjust the skewness of particle weights. Performance evaluations on the Gaussian barrier scattering (2D) and a Helium-like system (4D) demonstrate the feasibility of both strategies and the variance reduction property of the second approach. We provide an improvement of the first signed-particle implementation that partially alleviates the restriction on the time step and perform a thorough theoretical and numerical comparison among all the existing signed-particle implementations. Details on implementing the importance sampling according to the quasi-probability density and an efficient resampling or particle reduction are also provided.

The eigenvectors for graph 1-Laplacian possess some sort of localization property: On one hand, the characteristic function on any nodal domain of an eigenvector is again an eigenvector with the same eigenvalue; on the other hand, one can pack up an eigenvector for a new graph by several fundamental eigencomponents and modules with the same eigenvalue via few special techniques. The Courant nodal domain theorem for graphs is extended to graph 1-Laplacian for strong nodal domains, but for weak nodal domains it is false. The notion of algebraic multiplicity is introduced in order to provide a more precise estimate of the number of independent eigenvectors. A positive answer is given to a question raised in Chang (2016) [3], to confirm that the critical values obtained by the minimax principle may not cover all eigenvalues of graph 1-Laplacian.

We consider the nonlinear Dirac equation in 1  +  1 dimension with scalar–scalar self-interaction in the presence of external forces as well as damping of the form ${{\gamma}^{0}}f(x,t)-\text{i}\mu {{\gamma}^{0}} \Psi $ , where both $f,\left\{\,{{f}_{j}}={{r}_{j}}{{\text{e}}^{\text{i}{{K}_{j}}x}}\right\}$ and $ \Psi $  are two-component spinors. We develop an approximate variational approach using collective coordinates for studying the time dependent response of the solitary waves to these external forces. In our previous paper we assumed Kj  =  K, j  =  1, 2 which allowed a transformation to a simplifying coordinate system, and we also assumed the 'small' component of the external force was zero. Here we include the effects of the small component and also the case ${{K}_{1}}\ne {{K}_{2}}$  which dramatically modifies the behavior of the solitary wave in the presence of these external forces.

It is shown that relativistic many-body Hamiltonians and wave functions can beexpressed systematically with Tracy-Singh products for partitioned matrices. The latter gives rise to the usual notion for a relativistic $N$-electron wave function: A column vector composed of $2^N$ blocks, each of which consists of $2^N$ components formed by the Kronecker products of $N$ one-electron 2-spinors. Yet, the noncommutativity of the Tracy-Singh product dictates that the chosen serial ordering of electronic coordinates cannot be altered when antisymmetrizing a Tracy-Singh product of 4-spinors. It is further shown that such algebraic representation uncovers readily the internal symmetries of the relativistic Hamiltonians and wave functions, which are crucial for deriving the electron-electron coalescence conditions.

The electron-electron coalescence conditions for the wave functions of the Dirac-Coulomb (DC), Dirac-Coulomb-Gaunt (DCG), and Dirac-Coulomb-Breit (DCB) Hamiltonians are analyzed by making use of the internal symmetries of the reduced two-electron systems. The results show that, at the coalescence point of two electrons, the wave functions of the DCG Hamiltonian are regular, while those of the DC and DCB Hamiltonians have $r_{12}^{\nu}$ type weak singularities, with $\nu$ being negative and of $\mathcal{O}(\alpha^2)$. Yet, such asymptotic expansions of the relativistic wave functions are only valid within an extremely small convergence radius $R_c$ of $\mathcal{O}(\alpha^2)$. Beyond this radius, the behaviors of the relativistic wave functions are still dominated by the nonrelativistic limit.