Sihong Shao

Qiming Fang, Sihong Shao, Yuxuan Wu

Using the concepts of Eulerian-spanning set and coboundary operator, we generalize Hadlock's conversion of the maxcut problem on planar graphs to one on general graphs with non-negative weights. Using our conversion, we can explore algorithms for maxcut beyond the class of planar graphs. We obtain a Fixed-Parameter Tractable algorithm for $k$-contraction apex graphs. Specifically, our algorithm can be applied to graphs with crossing number $k$, giving an $O(2^k(n+k)^{3/2}\log (n+k))$-time algorithm that matches the best known results when restricted to non-negative weights.

Yanjie Li, Sihong Shao

This paper studies the Fourier spectral method (FSM) for the Schrödinger equation with singular potentials $V \in H^{s}$, where $s > \max\{d/2-2,-1\}$ and $d$ denotes the spatial dimension. This setting includes a broad class of singular potentials, such as the 3D Coulomb potential and the 1D Dirac-delta potential. First, we combine the Feshbach-Schur map with a refined perturbation argument to derive sharp convergence orders for FSM, yielding order $2s+2$ for eigenvalues and order $s+1$ for eigenfunctions in the $H^1$ norm. More importantly, the $H^1$ error with respect to the projected eigenfunction converges with a higher order $s+1+b$, where $b=\min\{s+2-d/2-\varepsilon,\; s+1,\; 2\}>0$ for arbitrarily small $\varepsilon>0$, revealing a super-convergence phenomenon. Second, in the presence of potentials with isolated point singularities, we develop an asymptotic-recovery (AR) technique to post-process the FSM solutions. The resulting method, dubbed AR-FSM, fully exploits the super-convergence property and achieves convergence orders $2s+2+2b$ for eigenvalues and $s+1+b$ for eigenfunctions in the $H^1$ norm, while the AR post-processing requires only a computational cost that is linear in the number of FSM degrees of freedom. The analysis introduces a rigorous definition of point singularities and develops a foundational framework for their study. It further establishes an asymptotic expansion of eigenfunctions consisting of a regular component in $H^{s+4}$ together with $d+1$ asymptotic functions associated with each singular point. Numerical experiments confirm the sharpness of these theoretical bounds.

Lin Lin, Sihong Shao, Weinan E

We present for the first time an efficient iterative method to directly solve the four-component Dirac-Kohn-Sham (DKS) density functional theory. Due to the existence of the negative energy continuum in the DKS operator, the existing iterative techniques for solving the Kohn-Sham systems cannot be efficiently applied to solve the DKS systems. The key component of our method is a novel filtering step (F) which acts as a preconditioner in the framework of the locally optimal block preconditioned conjugate gradient (LOBPCG) method. The resulting method, dubbed the LOBPCG-F method, is able to compute the desired eigenvalues and eigenvectors in the positive energy band without computing any state in the negative energy band. The LOBPCG-F method introduces mild extra cost compared to the standard LOBPCG method and can be easily implemented. We demonstrate our method in the pseudopotential framework with a planewave basis set which naturally satisfies the kinetic balance prescription. Numerical results for Pt$_{2}$, Au$_{2}$, TlF, and Bi$_{2}$Se$_{3}$ indicate that the LOBPCG-F method is a robust and efficient method for investigating the relativistic effect in systems containing heavy elements.

Yunfeng Xiong, Zhenzhu Chen, Sihong Shao

As a phase space language for quantum mechanics, the Wigner function approach bears a close analogy to classical mechanics and has been drawing growing attention, especially in simulating quantum many-body systems. However, deterministic numerical solutions have been almost exclusively confined to one-dimensional one-body systems and few results are reported even for one-dimensional two-body problems. This paper serves as the first attempt to solve the time-dependent many-body Wigner equation through a grid-based advective-spectral-mixed method. The main feature of the method is to resolve the linear advection in $(\bm{x},t)$-space by an explicit three-step characteristic scheme coupled with the piecewise cubic spline interpolation, while the Chebyshev spectral element method in $\bm k$-space is adopted for accurate calculation of the nonlocal pseudo-differential term. Not only the time step of the resulting method is not restricted by the usual CFL condition and thus a large time step is allowed, but also the mass conservation can be maintained. In particular, for the system consisting of identical particles, the advective-spectral-mixed method can also rigorously preserve physical symmetry relations. The performance is validated through several typical numerical experiments, like the Gaussian barrier scattering, electron-electron interaction and a Helium-like system, where the third-order accuracy against both grid spacing and time stepping is observed.

Zhenzhu Chen, Yunfeng Xiong, Sihong Shao

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. Using a spectral collocation discretization which 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.

Sihong Shao, Yunfeng Xiong

A branching random walk algorithm for the many-body Wigner equation and its numerical applications for quantum dynamics in phase space are proposed and analyzed. After introducing an auxiliary function, the (truncated) Wigner equation is cast into the integral formulation as well as its adjoint correspondence, both of which can be reformulated into the renewal-type equations and have transparent probabilistic interpretation. We prove that the first moment of a branching random walk happens to be the solution for the adjoint equation. More importantly, we detail that such stochastic model, associated with both importance sampling and resampling, paves the way for a numerically tractable scheme, within which the Wigner quantum dynamics is simulated in a time-marching manner and the complexity can be controlled with the help of an (exact) estimator of the growth rate of particle number. Typical numerical experiments on the Gaussian barrier scattering and a Helium-like system validate our theoretical findings, as well as demonstrate the accuracy, the efficiency and thus the computability of the Wigner branching random walk algorithm.

Sheng Chen, Sihong Shao, Shuai Wu

The high-precision solution of the radial Dirac equation is fundamental to relativistic quantum chemistry, essential for reliable pseudopotential generation and all-electron electronic structure methods. However, standard basis-set approaches struggle to simultaneously capture two distinct physical regimes: the non-polynomial singularities at the origin and the state-dependent, multi-scale asymptotic decay of wavefunctions on semi-infinite domains. In this work, we propose a high-precision adaptive spectral-element framework designed to rigorously resolve these spatial challenges. To capture the diverse exponential decay behavior on $[0, \infty)$ without arbitrary domain truncation, an adaptive generalized Laguerre spectral method is introduced, dynamically optimizing the basis scaling factors. Concurrently, near-origin non-polynomial {$r^s$} singularities are resolved utilizing Log-Orthogonal Functions, a basis that intrinsically approximates complex singular behaviors without requiring prior knowledge of the exact analytical exponent {$s$}. Furthermore, the framework incorporates an inverse operator formulation to guarantee spectral purity and eliminate spurious states. Validated across diverse physical regimes, including Coulomb, finite-nucleus, and screened potentials, the proposed method restores exponential convergence and consistently achieves relative accuracies of $10^{-10}$ {in Hartree atomic units or electron volts}. This work provides a robust, non-pollution computational kernel for atomic structure calculations, establishing a numerical standard for generating high-precision atomic data in complex molecular simulations.

Sihong Shao, Yanli Wang, Jie Wu

We propose an adaptive Hermite spectral method for the Vlasov-Poisson system based on a recently developed frequency indicator that measures the contribution of the high-order expansion coefficients. Precisely, the symmetrically weighted Hermite basis with a scaling factor is utilized to approximate the distribution function to satisfy the increasing resolution requirement, which, for example, is induced by filamentation. To implement the scaling adjustment, a fast conservative projection operator is constructed in two steps. The first step is to formulate the projection as a constrained optimization problem to preserve key invariants, including mass, momentum, energy, and the $L^2$ norm of the distribution function. The second step is an ODE-based approximation developed to compute the updated expansion coefficients with linear complexity. Numerical experiments with 1D1V and 2D2V settings validate the feasibility and efficiency of this proposed adaptive Hermite method.

Andrew Qing He, Wei Cai, Sihong Shao

We extend the Weak Adversarial Neural Pushforward Method to the Wigner transport equation governing the phase-space dynamics of quantum systems. The central contribution is a structural observation: integrating the nonlocal pseudo-differential potential operator against plane-wave test functions produces a Dirac delta that exactly inverts the Fourier transform defining the Wigner potential kernel, reducing the operator to a pointwise finite difference of the potential at two shifted arguments. This holds in arbitrary dimension, requires no truncation of the Moyal series, and treats the potential as a black-box function oracle with no derivative information. To handle the negativity of the Wigner quasi-probability distribution, we introduce a signed pushforward architecture that decomposes the solution into two non-negative phase-space distributions mixed with a learnable weight. The resulting method inherits the mesh-free, Jacobian-free, and scalable properties of the original framework while extending it to the quantum setting.

Jingyang Huang, Zhengyang Lei, Sihong Shao

In order to numerically solve high-dimensional nonlinear PDEs and alleviate the curse of dimensionality, a stochastic particle method (SPM) has been proposed to capture the relevant feature of the solution through the adaptive evolution of particles [J. Comput. Phys. 527 (2025) 113818]. In this paper, we introduce an active birth-death dynamics of particles to improve the efficiency of SPM. The resulting method, dubbed SPM-birth-death, sample new particles according to the nonlinear term and execute the annihilation strategy when the number of particles exceeds a given threshold. A rigorous error estimation for SPM-birth-death is established, elucidating the first-order convergence in time and space, as well as half-order accuracy in the initial sample size with explicit variance estimates. We also extend the analysis framework to SPM and provide theoretical justification for the existing numerical convergence study. Our theoretical results reveal that the introduced active birth-death dynamics of particles results into less frequent resampling and SPM-birth-death is thus able to achieve higher efficiency than SPM. Validating benchmarks are provided. In particular, preliminary numerical experiments on the Allen-Cahn equation demonstrate that SPM-birth-death can achieve smaller errors at the same computational cost compared with the original SPM.

Zhengyang Lei, Sihong Shao

With the aim of generalizing histogram statistics to higher dimensional cases, density estimation via discrepancy based sequential partition (DSP) has been proposed [D. Li, K. Yang, W. Wong, Advances in Neural Information Processing Systems (2016) 1099-1107] to learn an adaptive piecewise constant approximation defined on a binary sequential partition of the underlying domain, where the star discrepancy is adopted to measure the uniformity of particle distribution. However, the calculation of the star discrepancy is NP-hard and it does not satisfy the reflection invariance and rotation invariance either. To this end, we use the mixture discrepancy and the comparison of moments as a replacement of the star discrepancy, leading to the density estimation via mixture discrepancy based sequential partition (DSP-mix) and density estimation via moments based sequential partition (MSP), respectively. Both DSP-mix and MSP are computationally tractable and exhibit the reflection and rotation invariance. Numerical experiments in reconstructing the $d$-D mixture of Gaussians and Betas with $d=2, 3, \dots, 6$ demonstrate that DSP-mix and MSP both run approximately ten times faster than DSP while maintaining the same accuracy.