Publications

A numerical scheme is presented for solving the Helmholtz equation with Dirichlet or Neumann boundary conditions on piecewise smooth open curves, where the curves may have corners and multiple junctions. Existing integral equation methods for smooth open curves rely on analyzing the exact singularities of the density at endpoints for associated integral operators, explicitly extracting these singularities from the densities in the formulation, and using global quadrature to discretize the boundary integral equation. Extending these methods to handle curves with corners and multiple junctions is challenging because the singularity analysis becomes much more complex, and constructing high-order quadrature for discretizing layer potentials with singular and hypersingular kernels and singular densities is nontrivial. The proposed scheme is built upon the following two observations. First, the single-layer potential operator and the normal derivative of the double-layer potential operator serve as effective preconditioners for each other locally. Second, the recursively compressed inverse preconditioning (RCIP) method can be extended to address “implicit” second-kind integral equations. The scheme is high-order, adaptive, and capable of handling corners and multiple junctions without prior knowledge of the density singularity. It is also compatible with fast algorithms, such as the fast multipole method. The performance of the scheme is illustrated with several numerical examples.

We generalize the interpolative separable density fitting (ISDF) method, used for compressing the four-index electron repulsion integral (ERI) tensor, to incorporate adaptive real space grids for potentially highly localized single-particle basis functions. To do so, we employ a fast adaptive algorithm, the recently-introduced dual-space multilevel kernel-splitting method, to solve the Poisson equation for the ISDF auxiliary basis functions. The adaptive grids are generated using a high-order accurate, black-box procedure that satisfies a user-specified error tolerance. Our algorithm relies on the observation, which we prove, that an adaptive grid resolving the pair densities appearing in the ERI tensor can be straightforwardly constructed from one that resolves the single-particle basis functions, with the number of required grid points differing only by a constant factor. We find that the ISDF compression efficiency for the ERI tensor with highly localized basis sets is comparable to that for smoother basis sets compatible with uniform grids. To demonstrate the performance of our procedure, we consider several molecular systems with all-electron basis sets which are intractable using uniform grid-based methods. Our work establishes a pathway for scalable many-body electronic structure simulations with arbitrary smooth basis functions, making simulations of phenomena like core-level excitations feasible on a large scale.

Having an accurate reaction coordinate (RC) is essential for reliable kinetic characterization of molecular processes, but there are few quantitative metrics to evaluate RC quality. In this study, we consider the dimensionless γ metric from the Exponential Average Time-dependent Rate (EATR) method, which represents the fraction of a biasing potential along the RC that contributes to increasing the rate constant. We demonstrate that γ can be used to test whether the utility of a RC for predicting kinetics with a Metadynamics bias improves as the coordinate is iteratively updated to include new data. We evaluate RCs approximated via the iterative State Predictive Information Bottleneck (SPIB) approach, which was previously shown to be accurate across six protein–ligand dissociation systems. For these same systems, we compute γ values and mean accelerated times τ̅accel. After systematically scanning over fitting parameters, the results show that γ increases closer to 1, while τ̅accel decreases, revealing a consistent inverse correlation. These results demonstrate that γ serves as a practical criterion for RC evaluation and offers guidance for selecting SPIB–derived coordinates yielding quantitative kinetic predictions.

The method of fundamental solutions (MFS) is known to be effective for solving 3D Laplace and Stokes Dirichlet boundary value problems in the exterior of a large collection of simple smooth objects. Here, we present new scalable MFS formulations for the corresponding elastance and mobility problems. The elastance problem computes the potentials of conductors with given net charges, while the mobility problem—crucial to rheology and complex fluid applications—computes rigid body velocities given net forces and torques on the particles. The key idea is orthogonal projection of the net charge (or forces and torques) in a rectangular variant of a “completion flow.” The proposal is compatible with one-body preconditioning, resulting in well-conditioned square linear systems amenable to fast multipole accelerated iterative solution, thus a cost linear in the particle number. For large suspensions with moderate lubrication forces, MFS sources on inner proxy-surfaces give accuracy on par with a well-resolved boundary integral formulation. Our several numerical tests include a suspension of 10,000 nearby ellipsoids, using 2.6\times 10^7
total preconditioned degrees of freedom, where GMRES converges to five digits of accuracy in under two hours on one workstation

We reveal a mechanism that may explain the non-metallic out-of-plane resistivity in layered metals. By carefully examining how the Drude-Boltzmann expression for the $c$-axis conductivity emerges out of the Kubo formula, we find, besides the standard metallic term proportional to the carrier lifetime $\tau$, a non-Drude contribution proportional to $1/\tau$. The Drude behavior breaks down when $1/\tau > 2 \eta^*$, the crossover value $\eta^*$ being small (and hence observable) when the $c$-axis velocities vary rapidly with the distance from the Fermi surface. We consider the Hund metal Sr$_2$RuO$_4$ as a test case, which we study within a realistic dynamical mean-field theory approach. The non-Drude behavior observed experimentally in $c$-axis transport is reproduced and explained by our considerations, showing that earlier invoked extrinsic mechanisms that involve either impurities or phonons are unnecessary. We point out that the small value of $\eta^*$ is due to a peculiar accidental cancellation due to destructive interference characteristic of body-centered tetragonal lattices.

We explore a quaternion adjugate matrix-based representation for rotational motion in the Perspective-n-Point (PnP) problem. Leveraging quadratic quaternion terms within a Determinant Ratio Matrix (DRaM) estimation framework, we extend its application to perspective scenarios, providing a robust and efficient initialization for iterative PnP pose estimation. Notably, by solving the orthographic projection least-squares problem, DRaM provides a reliable initialization that enhances the accuracy and stability of iterative PnP solvers. Experiments on synthetic and real data demonstrate its efficiency, accuracy, and robustness, particularly under high noise conditions. Furthermore, our non-minimal formulation ensures numerical stability, making it effective for real-world applications.

This work addresses the question of uniqueness and regularity of the minimizers of a convex but not strictly convex integral functional with linear growth in a two-dimensional setting. The integrand whose precise form derives directly from the theory of perfect plasticity behaves quadratically close to the origin and grows linearly once a speci c threshold is reached. Thus, in contrast with the only existing literature on uniqueness for functionals with linear growth, that is that which pertains to the generalized least gradient, the integrand is not a norm. We make use of hyperbolic conservation laws hidden in the structure of the problem to tackle uniqueness. Our argument strongly relies on the regularity of a vector eld the Cauchy stress in the terminology of perfect plasticity which allows us to de ne characteristic lines, and then to employ the method of characteristics. Using the detailed structure of the characteristic landscape evidenced in our preliminary study [5], we show that this vector eld is actually continuous, save for possibly two points. The di erent behaviors of the energy density at zero and at innity imply an inequality constraint on the Cauchy stress. Under a barrier type convexity assumption on the set where the inequality constraint is saturated, we show that uniqueness holds for pure Dirichlet boundary data devoid of any regularity properties, a stronger result than that of uniqueness for a given trace on the whole boundary since our minimizers can fail to attain the boundary data. We also show a partial regularity result for the minimizer.

As LLMs are increasingly deployed as agents, agentic reasoning - the ability to combine tool use, especially search, and reasoning - becomes a critical skill. However, it is hard to disentangle agentic reasoning when evaluated in complex environments and tasks. Current agent benchmarks often mix agentic reasoning with challenging math reasoning, expert-level knowledge, and other advanced capabilities. To fill this gap, we build a novel benchmark, GSM-Agent, where an LLM agent is required to solve grade-school-level reasoning problems, but is only presented with the question in the prompt without the premises that contain the necessary information to solve the task, and needs to proactively collect that information using tools. Although the original tasks are grade-school math problems, we observe that even frontier models like GPT-5 only achieve 67% accuracy. To understand and analyze the agentic reasoning patterns, we propose the concept of agentic reasoning graph: cluster the environment's document embeddings into nodes, and map each tool call to its nearest node to build a reasoning path. Surprisingly, we identify that the ability to revisit a previously visited node, widely taken as a crucial pattern in static reasoning, is often missing for agentic reasoning for many models. Based on the insight, we propose a tool-augmented test-time scaling method to improve LLM's agentic reasoning performance by adding tools to encourage models to revisit. We expect our benchmark and the agentic reasoning framework to aid future studies of understanding and pushing the boundaries of agentic reasoning.

We analyze and develop numerical methods for time-harmonic wave scattering in metallic waveguide structures of infinite extent. We show that radiation boundary conditions formulated via projectors onto outgoing modes determine the coefficients of propagating modes uniquely, even when the structure supports trapped modes. Building on this, we introduce a fast divide-and-conquer solver that constructs solution operators on subdomains as impedance-to-impedance maps and couples them by enforcing continuity conditions across their interfaces. For Dirichlet waveguides, the computation of impedance-to-impedance maps requires the solution of mixed Dirichlet-Impedance boundary value problems. We construct a second-kind Fredholm integral equation that avoids near-hypersingular operators, requiring only integral operators whose kernels are at most weakly singular. Numerical experiments on large structures with many circuit elements demonstrate substantial efficiency gains: the proposed approach typically outperforms state-of-the-art fast iterative and fast direct solvers by one to two orders of magnitude.

Classical Ewald methods for Coulomb and Stokes interactions rely on ``kernel-splitting," using decompositions based on Gaussians to divide the resulting potential into a near field and a far field component. Here, we show that a more efficient splitting for the scalar biharmonic Green's function can be derived using zeroth-order prolate spheroidal wave functions (PSWFs), which in turn yields new efficient splittings for the Stokeslet, stresslet, and elastic kernels, since these Green's tensors can all be derived from the biharmonic kernel. This benefits all fast summation methods based on kernel splitting, including FFT-based Ewald summation methods, that are suitable for uniform point distributions, and DMK-based methods that allow for nonuniform point distributions. The DMK (dual-space multilevel kernel-splitting) algorithm we develop here is fast, adaptive, and linear-scaling, both in free space and in a periodic cube. We demonstrate its performance with numerical examples in two and three dimensions.