Papers

We present release 3.0 of the ALPS (Algorithms and Libraries for Physics Simulations) project, an open-source software project to develop libraries and application programs for the simulation of strongly correlated quantum lattice models such as quantum magnets, lattice bosons, and strongly correlated fermion systems. As in previous releases, development is centered on common data formats, on libraries to simplify and speed up code development, and on full-featured simulation programs that let non-experts carry out serial or parallel numerical simulations using the important algorithms for quantum lattice models: classical and quantum Monte Carlo (QMC) using non-local updates, extended-ensemble simulations, exact and full diagonalization (ED), the density matrix renormalization group (DMRG), and continuous-time QMC solvers for dynamical mean-field theory (DMFT). Major changes in release 3.0 include distribution of the pyalps binary through the Python Package Index (pip install pyalps) and through Spack for HPC systems; migration of development to GitHub with continuous integration and automated testing; relicensing of the package under the permissive MIT license; a rebuilt documentation and tutorial website, including a set of Jupyter-notebook tutorials and localized content; a broad modernization of the C++ codebase (C++17 compliance, Boost and NumPy 2.0 compatibility, and warning and dead-code cleanup) together with a major DMRG update and associated reliability and build-compatibility fixes; the removal of legacy components (the VisTrails provenance integration and the TEBD, MPS, and directed-worm-algorithm application codes); archival of release 3.0.0 with a Zenodo DOI; and a formal governance and sustainability model developed under the US National Science Foundation (NSF) POSE program. The software is available at https://alps.comp-phys.org.

Purpose: This study extends our fluid-structure interaction (FSI) framework for surgical bioprosthetic heart valves in an experimental pulse duplicator to a transcatheter aortic valve replacement (TAVR) device. We assess hemodynamic and leaflet kinematic agreement across pulse rates and the transferability of calibrated boundary-model parameters. Methods: We compared simulated and experimental pressure, volumetric flow rate, projected dynamic valve area (PDVA), and leaflet kinematics at 60, 70, and 80 beats per minute (bpm). Reduced-order boundary-model parameters were calibrated using only data acquired at 70 bpm and were subsequently held fixed. At each pulse rate, a condition-specific pump-pressure waveform was constructed from the measured upstream pressure and flow rate. The 60 and 80 bpm cases therefore evaluate parameter transfer under measurement-informed upstream forcing. Results: Across the three conditions, normalized root-mean-square errors are 2.88-3.08% for aortic flow rate, 7.64-8.07% for upstream ventricular pressure, 3.50-6.21% for downstream aortic pressure, and 3.07-4.70% for PDVA. The simulations reproduce the principal pressure, flow, valve-opening, and valve-closure features. The downstream aortic pressure and PDVA measurements at 60 and 80 bpm were not used to construct the upstream forcing or recalibrate the model parameters and therefore provide the principal out-of-calibration validation evidence. Conclusion: The results provide conditional experimental bench validation evidence for downstream hemodynamics and valve dynamics under measurement-informed upstream forcing and support boundary-parameter transfer across pulse rates without recalibration. This work advances an experimentally supported, high-fidelity FSI framework for evaluating TAVR device performance.

Charge deposition in the particle-in-cell (PIC) method can be viewed as kernel density estimation on a spatial mesh. This allows the smoothing bandwidth to be selected through statistical criteria such as mean integrated squared error or cross-validation, but it also changes the electric field acting on the particles. Therefore, a selected bandwidth can alter the growth or damping of physical modes. We develop a dynamics-aware bandwidth-selection procedure that combines a density-error objective with constraints on specified linear modes. These constraints follow from the plasma equilibrium and the combined transfer of deposition, filtering and field gathering within the model rate-tolerance boundary. We demonstrate the need for equilibrium information through a controlled two-stream experiment. The controlled two-stream scan varies the beam drift while holding the initial spatial density, thermal samples, wavevector and fixed-filter transfer unchanged. Across four resolved equilibria approaching marginal stability, this filter reduces the measured growth rate by $0.7\%$ to $9.1\%$. Equilibrium-dependent interior settings satisfy the prescribed $10\%$ rate tolerance. We then apply the constrained selection rule to 85 Gaussian candidates, that reduces density and grid-field mean-square errors by $7.40\%$ and $3.70\%$ on independently evaluated snapshots. A separate projected-linear calculation at that exact width gives a ratio of mean individual growth rates of $0.9556$, with paired $95\%$ interval $[0.9453,0.9640]$, consistent with the prediction $0.9513$. Analysis of finite-marker loading and ensemble averaging explains how the measured rates relate to the response model. Together, these results demonstrate a filter-design procedure that combines statistical error reduction with subsequent validation of the selected linear dynamics.

GC–EI–MS is an important technique for analyzing volatile and semivolatile compounds in complex samples. However, conventional methods rely heavily on reference spectral library matching, limiting their ability to identify compounds absent from these libraries and to infer complete molecular structures directly from fragmentation information. Here, we present DiffGCMS, a spectrum-conditioned discrete graph diffusion model for de novo structure elucidation from GC–EI–MS, and further develop a framework that integrates DiffGCMS with second-stage reasoning by a large language model (LLM). In the first stage, DiffGCMS generates candidate molecular structures from input spectra; in the second stage, the LLM uses mass spectral information to validate, repair, and rerank the candidates and provides interpretable analysis of fragment-ion peaks. This framework can generate plausible molecular structures for compounds absent from reference spectral libraries and provide traceable evidence supporting its decisions. On a test set comprising 13,696 spectra from NIST 20, the generative model achieved Acc@1 and Acc@10 of 6.01% and 15.76%, respectively. On the test subset containing molecules with no more than 10 heavy atoms, LLM-assisted molecular graph repair and reranking increased Acc@1 from 21.28% to 21.95%, Acc@10 from 46.91% to 47.99%, and candidate validity from 91.04% to 100%. These results demonstrate that spectrum-aware postprocessing can correct errors produced by the generative model while providing auditable and traceable explanations for the final ranking.

Inverse design has been transformative in nanophotonics, providing an automated means of realizing high-performance, non-intuitive devices. However, the primary bottleneck of existing inverse design approaches is their reliance on time-consuming and computationally expensive full-wave electromagnetic simulations. In this work, we extend the Precomputed Numerical Green Function (PNGF) method for the first time to nanophotonic inverse design, accelerating evaluation of the forward problem by multiple orders of magnitude with no loss in accuracy compared to the full-wave solution. After a single, fully parallelizable precomputation step, both the objective function and its gradient can be evaluated through simple calculations whose cost scales linearly with the size of the design region. A low-rank matrix update technique further reduces the cost of objective function evaluation, yielding sub-millisecond computation times per iteration. Design examples using both gradient-based level-set optimization and tile-flipping direct binary search often achieve more than three orders of magnitude speedup in design time, resulting in an ultrafast photonic inverse design platform that can design new devices from scratch in seconds.

The stochastic variational method (SVM) is one of the most powerful methods to solve quantum few-body systems precisely in various fields, such as nuclear and atomic physics. To the best of our knowledge, no SVM code optimized for GPUs has been reported. We developed a parallel SVM code for few-body clusters of $^4$He atoms and optimized it for two GPU architectures, NVIDIA GH200 and AMD Instinct MI300A, with the entire code written by an AI coding agent. On the algorithmic side, we introduce the secular-equation method with the Gu–Eisenstat prescription for solving the generalized eigenvalue problem within the SVM framework. The main part of the code tuning is to batch the many small matrices so that the GPUs are used efficiently, together with an array layout in memory optimized for coalesced addressing. For the five-body system with 4800 SVM basis states, the tuned SVM code runs 16.2 (14.7) times faster on the NVIDIA GH200 (AMD MI300A) than the same tuned code on the Intel Xeon Max CPU system, and 168 (153) times faster than the first working CPU implementation. The achieved performance brings 6-body and larger cluster systems within reach.

We develop a finite element method (FEM) for partial differential equation (PDE) solver based on the entanglement-basis representation introduced in our companion work. By lifting non-linear finite-element equations into an augmented coefficient space, the governing PDE together with boundary, initial, and inter-element constraints can be expressed through a unified quadratic residual minimization. Although this augmented space grows exponentially with the number of elements, its tensor-product structure allows it to be represented efficiently using tensor networks. Using the matrix product state (MPS) as a concrete example, we show that density matrix renormalization group (DMRG) sweeps enable element-by-element optimization without explicitly constructing the full augmented space. For bounded bond dimension, the resulting computational cost scales polynomially with the number of finite elements. We extend the framework to time-dependent problems through implicit temporal discretization and demonstrate convergence under both mesh and polynomial refinement using diffusion equations.

Equivariant machine learning interatomic potentials (MLIPs) have revolutionized atomistic modeling, but accurate treatment of complex materials and molecular systems demands expensive models. This limits simulation length- and time-scales, with tensor products a key computational bottleneck. The recent emergence of foundation-scale MLIPs further exacerbates this challenge. We present Branch Interatomic Potential (BranchIP), a single-model framework for learned adaptive tensor product computation, trained with a novel distillation loss. In our experiments on two systems of physical interest, a heterogeneous catalysis system and a proton-conducting solid acid electrolyte, BranchIP accelerates MLIPs across model sizes by up to $2.4\times$ while reducing memory usage by up to $2.6\times$. This is achieved while maintaining physical fidelity. Furthermore, the learned adaptive computation provides model interpretability by revealing which interactions demand deeper computation and showing how computational depth relates to chemical complexity and dynamics.

Ultrasound computed tomography is emerging as a promising diagnostic imaging tool. 2D geometries suffer from notorious out-of-plane scattering artifacts. Image reconstruction can be achieved with frequency-domain full waveform inversion and it can be further accelerated by randomly phase-encoding the elementary sources. In this manuscript, we extend our previous results for the 2D geometry of a ring-array to the 3D cylindrical geometry of multiple rings. In particular, we consider the cases of an elementary source described by a single array element (quasi-omni-directional transmission), a line source and a focused transmission respectively. With differences in image quality, we prove that a fully randomized frequency-domain inversion in 3D is capable to reconstruct portions of a human breast surrounded by the cylindrical geometry and detect mm-size masses of varying contrast in dense breast, in reasonable computing times, thus opening the concrete possibility to the design of 3D imaging devices with high sensitivity levels. The methods are applicable to multiple tomographic geometries in 3D and, in principle, can be integrated into next generation medical ultrasound scanners.

We present RFBniCS, an open-source finite-element framework implemented in FEniCSx for simulating redox-flow-battery (RFB) half-cells. RFBniCS solves an established macro-homogeneous porous-electrode model accounting for strongly coupled electrolyte flow, multicomponent species transport, ionic and electronic charge conservation, and interfacial Faradaic charge transfer. Different from other open-access tools, RFBniCS can perform transient RFB half-cell simulations in one-, two-, and three-dimensional geometries while explicitly resolving the transport of both redox-active and supporting-electrolyte species. It also describes electrolyte flow through porous electrodes and, in the three-dimensional formulation, through adjacent flow channels. To verify the implementation, RFBniCS's predictions in parameter-limiting regimes are compared against other published implementations, with RFBniCS generally showing superior accuracy and speed. We further demonstrate RFBniCS's capabilities by a simulation of the transient response of the negative half-cell of a vanadium redox flow battery, with a moderately concentrated supporting electrolyte. In the current implementation, RFBniCS provides the computational basis for studying redox-flow-battery half-cells and can be extended to complex flow-field designs, alternative chemistries, full-cell coupling, detailed membrane transport, and additional multiphysics effects.

This paper presents a cell-centred unstructured finite-volume method for time-harmonic two-phase acoustics. A volume-fraction representation supplies density and acoustic compressibility. Face fractions average available neighbouring planar interface cuts independently of velocity direction, retaining the established one-field face-density closure. Cartesian complex-coordinate stretching provides perfectly matched layers (PMLs). Real-imaginary splitting yields a real-valued block system, discretized with corrected non-orthogonal fluxes and consistent boundary contributions. Verification covers homogeneous waves, layered gas-liquid transmission with PML truncation, baffled-piston radiation with Kirchhoff far-field reconstruction, and resolved rigid-sphere radiation forces. Homogeneous-wave pressure converges at approximately second order on orthogonal meshes and on meshes with non-orthogonal interiors and orthogonal boundary cells. Reconstructed velocity converges at approximately second order on orthogonal meshes and with an observed order of about 1.7 on the latter mesh family. Under aligned refinement, the layered case's relative whole-domain complex-pressure $L_2$ error decreases to $2.657\times10^{-4}$. The resolved-sphere force differs from the Gorkov prediction by at most 1.5% over the tested Rayleigh size range. The results quantify the accuracy and current limitations of the heterogeneous Helmholtz-PML formulation.

The paper proposes a generalization of the fully implicit energy- and charge-conserving electromagnetic particle-in-cell method to the case where the lightest type of plasma particle (electrons) is described in the drift-kinetic approximation. This allows us to remove the very strict time step limitation of this method requiring to resolve the gyrorotation of electrons. Since the drift-kinetic model is only applicable to light particles whose contribution to plasma polarization is small, this model can be further simplified by neglecting the electron polarization drift and including in Maxwell's equations, in addition to the current of gyrocenters, only the magnetization current. In order for such a hybrid model to retain conservative properties in finite-difference form, a method of self-consistent interpolation of $\nabla B$ in the mirror force from grid to particle and the magnetization vector from particle to grid is proposed. Unlike existing drift-kinetic models, this model is not limited to considering small perturbations near a given equilibrium, and therefore allows one to study the formation of plasma equilibria in regimes with a finite ratio of plasma to magnetic field pressure. Testing of the parallel code implemented in C++ using the PETSc library confirmed the fulfillment of the finite-difference laws of energy and charge conservation, as well as the ability of the model (in the drift-kinetic version for all types of particles) to correctly reproduce the diamagnetic effect and longitudinal ion-acoustic wave.

This paper reports a computational scheme which predicts the light-matter interaction in multiple mesoscopic plasmonic nanowires (NWs), each of which includes concentric or eccentric constituent wires. Lying at the core of the proposed scheme is the S-Matrix of each individual cylindrical interface. The optical response of the materials inside and outside each interface can be described by three popular mesoscopic models: the Hydrodynamic Drude model (HDM) and its diffusive variant, namely the Generalized Nonlocal Optical Response (GNOR) model, as well as the Surface Response Model (SRM). To take into account the relative displacements between the cylindrical interfaces, the addition theorem is applied to construct the main equation behind the algorithm, where the possibility of longitudinal waves intrinsic in the HDM and GNOR models is also included. The optical responses calculated by the proposed S Matrix Method are compared with an in-house developed boundary element method (BEM) toolbox for a single gold NW and an eccentric dual-material NW, and excellent agreements are seen. The algorithm is further physically checked for a sodium nanolens trimer featuring large field focusing, and the effects of the HDM and the SRM are compared.

We present CAFES, a correlation-aware agentic framework for electronic-structure simulations of strongly correlated systems. CAFES addresses two challenges: the fragmented software landscape for many-body calculations and the difficulty of selecting appropriate methods across diverse correlation regimes. It combines correlation diagnostics, adaptive method selection, and large language model (LLM) assistance for molecular, crystalline, and model-Hamiltonian systems. A study-task architecture separates study-level planning from task-level execution, while a curated scientific knowledge layer provides reusable guidance for method selection, workflow design, and result interpretation. We demonstrate CAFES through three research-level studies: calculating the low-lying electronic states of lutein using DMRG-CASSCF, probing phase competition in the extended honeycomb Hubbard model using DMET, and generating a quantum-chemical dataset with CCSD labels. These calculations demonstrate the potential of agentic workflows for strongly correlated electronic-structure problems, a regime that has received limited attention in existing agentic computational frameworks.

Machine-learned interatomic potentials are commonly fitted by weighted-sum scalarization, which combines energy and force errors in a single loss. A nominal weight, however, identifies a potential only relative to the complete fitting protocol. We therefore treat energy–force balancing as a protocol-dependent problem of physical model selection. For fixed-basis atomic cluster expansion models of molten LiCl, liquid H$_2$O, and Si, the resulting scalarization paths contain dominated states at extreme weights. Their nondominated subsets depend on the solver, and the out-of-distribution Si path is nonmonotone. We replace direct weight selection by minimizing the regularized energy objective subject to an upper bound on the normalized force loss. Projected dual ascent adjusts the Lagrange multiplier from the force-constraint residual. A frozen-multiplier limited-memory quasi-Newton refinement then returns the final model. Across the three systems, the constrained procedure reaches the solver-matched low-error region and gives measured speedups of $6.7$–$8.9$ over completed scalarization scans under the stated timing convention. Physical-property calculations show that first-shell geometry is comparatively insensitive to the selected balance, whereas transport and solid-state observables vary more strongly. These results identify the fitted model–protocol pair as the relevant object of energy–force model selection and support force-constrained fitting as an explicit selection rule.

SIREN (Ship Inviscid Resistance and Elevation eNgine) is a Rankine-source panel solver for the steady wave resistance of a ship, in the tradition of Dawson's method. The free-surface condition is linearised about a base flow, and the linearisation is iterated. In the present formulation the base flow is the solution of a Neumann problem on the actual boundary: sources on the entire hull, wet and dry, and on the free surface panels, with impermeability imposed on both, so that the free surface acts as a rigid lid at the position it currently occupies. This replaces the classical double-body flow, whose image about z=0 is meaningful only for a flat surface and a hull cut at z=0, and it removes the need for desingularised free-surface panels. In the non-linear iteration the free-surface panels are moved rigidly in the vertical to the relaxed wave elevation at their centroid, the hull panels are marked wet or dry against the moved surface, and the base flow is recomputed on the moved geometry; impermeability applies to the whole hull regardless of this marking, which is used only to select which panels enter the force integral, and under-relaxation is confined to the passes where the non-linear correction accumulates without being reset. The radiation condition is imposed by a one-row upstream shift of the collocation points, and the wave resistance is obtained by pressure integration and, independently, by transverse wave-cut analysis.

A majority of measurable dynamic quantities of a physical system is described as an equilibrium correlation function of the form $G_{AB}(t-t') = \left<A(t)B(t')\right>$, where $t'$ represents the impact time and $t$ represents the response time. Conventional quantum algorithms to compute such quantities via Hamiltonian simulation such as the Hadamard test, variational methods, or linear-response-based algorithms share one feature in common: each time point $t-t'$ is calculated by a different quantum circuit, leading to at least $O(N_t)$ total shots, and at least $O(N_t^2)$ total quantum gates. Here, we provide a different approach: utilizing $O(\log N_t)$ ancilla qubits to store time as a quantum variable, and effectively parallelize the computation of an equilibrium correlation function by requiring only a single quantum circuit for all time points. We show that this circuit needs to be run only $O(\log N_t)$ times, which in total requires $O(N_t \log N_t)$ quantum resources in the large system size limits. We prove that this scaling is optimal, and demonstrate our algorithm calculating the Green's function of a one-dimensional spinless Hubbard model.

Directly solving the full electron-nuclear Schrödinger equation remains one of the grand challenges in quantum mechanics. Here, we propose a quantum-chemically motivated orbital-based neural network wavefunction ansatz with explicit many-body electron-electron correlation terms, with parameters depending on the entire electron-nuclear configuration, enabling an accurate description of strongly correlated and/or multi-reference electronic structures by a single Slater-determinant across a broad nuclear configuration space. To train this wavefunction, we further develop an efficient hybrid optimization strategy that combines the variational Monte Carlo technique with local-energy constraints, mitigating the intrinsic sampling bias of the former and considerably reducing the statistical error. These advances are embedded in a unified neural network framework, Schrödinger, enabling the determination of globally accurate potential energy surfaces with a single training for systems ranging from strongly correlated diatomic molecules, triatomic reactions containing conical intersections, to polyatomic multi-reference molecules. This work offers a practical machine-learning tool for solving the electron-nuclear Schrödinger equation and capturing full quantum effects beyond the Born-Oppenheimer approximation.

Beyond a critical Toner–Tu drive, bacterial turbulence exhibits universal statistics, including a parameter-independent energy spectrum $E(k)\sim k^{-3/2}$. We show that the same universal state can be reached not only by increasing activity, but also by strengthening nonlinear advection, weakening nonlinear saturation, or tuning the Swift–Hohenberg instability. Motivated by the competition between these processes, we identify the independent dimensionless groups governing the dynamics and construct a Reynolds-like control parameter from them. A single critical value, calibrated from the known activity-controlled transition, predicts the onset of universality along the other parameter routes and collapses them onto a common threshold. These results identify a dynamical organizing parameter for the onset of universal active turbulence and provide a Reynolds-like measure for the Toner–Tu–Swift–Hohenberg universality class.

Crystal structure governs material properties, making crystal structure prediction (CSP) a fundamental problem in materials science. Generative models are a promising approach for solving this problem, but the prevalence of polymorphism, coupled with large unit cells and complex packing geometry, makes the molecular CSP task challenging for existing models. To address this, we introduce Coarse-Grained Open Materials Generation (CG-OMatG), an equivariant Riemannian flow-based generative model. CG-OMatG predicts molecular crystal structures via a coarse-grained, hierarchical representation. CG-OMatG treats molecules as rigid bodies—performing both inter- and intra-molecular message passing to construct a geometric representation for molecular packings—and learns to reconstruct molecule centroid positions, orientations, and lattice parameters, conditioned on chemical species and conformer geometry. We train the model on subsets of the Open Molecular Crystals (OMC25) and Cambridge Structural Database (CSD) datasets. Further, we fine-tune the model via policy gradient reinforcement learning to steer the model towards generating low-energy candidate structures. We validate the generated structures on the CSP blind test benchmark, assessing agreement with experimentally determined crystals using COMPACK packing-similarity analysis. CG-OMatG exhibits strong performance for generative molecular crystal structure prediction, paving the way for accelerated polymorph screening and organic solid-state materials discovery.

The mixing-plane method plays a pivotal role in steady simulations of multiple turbomachinery components. With advances in computational power, a whole-engine gas turbine simulation is no longer an elusive approach. However, whole-engine simulations require simultaneous coupling of the compressor, turbine and combustor. A key challenge is that the working fluid after the combustor can no longer be treated as a perfect gas, since the combustion introduces composition variations and combustion products. Conventional mixing-plane methods based on a fixed-composition gas cannot guarantee a thermodynamically consistent mixed-out state. To overcome the difficulty, this paper extends the classic mixing-plane approach to handle compressible reacting flows. The nonlinearity in the thermodynamic closure is addressed using a novel nested algorithm combining an outer pressure-root search with an inner enthalpy inversion. Physical state checks and power-law pressure sampling are incorporated to improve numerical robustness at low normal velocities where the pressure root approaches its upper physical bound. The performance of the method is demonstrated in three configurations: a quasi-1D interface test elucidating the thermodynamic stiffness of the proposed mixing-plane formulation for different fuel types and flow conditions; the Darmstadt transonic compressor case that verifies the proposed method reduces to the classic mixing-plane formulation when the mixture fraction reduces to zero; and a whole-engine simulation of the KJ66 micro-turbojet is used to demonstrate the performance of the method for multiphysics turbomachinery simulations. Good agreement with experimental data is observed, and the maximum relative mass-flow error across the mixing planes remains below 0.11%, demonstrating conservative and thermodynamically consistent interface coupling for compressible reacting turbomachinery flows.

Melting temperature calculation is one of the important topics in computational materials science. In high-throughput in silico screening and artificial intelligence assisted design of materials, it usually requires a rapid and autonomous assessment of the melting temperature of the target. Unfortunately, molecular dynamics (MD) simulations of the melting point require many cumbersome and manual operations, making large-scale calculation of the melting point challenging. In this work, we introduce MyTm, a toolkit that employs MD to automatically determine the melting point. The method is fully modularized, and by combining these modules, the program enables fully automated melting calculations by using commonly adopted approaches, including the direct-heating method, the void method, the modified void method, the solid-liquid coexistence method, and the Z method. Moreover, a machine learning (ML) method is proposed and employed to recognize and classify the solid like and liquid like atoms, which effectively resolve the low accuracy issue in conventional classification approaches, thus making the automated high throughput pipeline of melting-point calculation possible. The robustness and efficacy of MyTm have been demonstrated by several well studied systems.

Classical phase-field models typically use the Allen–Cahn evolution law to model interface motion as a consequence of energy-gradient descent. This implies a linear kinetic relation between the velocity of an interface (e.g., a domain wall in a ferroelectric) and its driving force. When nonlinear kinetic relations are to be modeled, as commonly found in ferroelectric ceramics, an alternative evolution law and hence an alternative model is needed. Here, we propose an electromechanically (stress-driven) coupled multiphase-field framework with a general kinetic formulation, which integrates a prescribed kinetic relation directly into the evolution law. We show that this model correctly evolves domain walls with the assigned nonlinear kinetics, e.g., of mixed exponential-power law type, following the so-called Merz–Stadler law. The multiphase formulation further enables distinct kinetic relations and interfacial energies to be assigned to different order-parameter pairs, which reflects the distinguishable properties of the different types of ferroelectric domain walls. The proposed framework is broadly applicable for simulating kinetic relations in materials with applications beyond ferroelectrics.

Explicit lattice Boltzmann simulations of slow transient flows are constrained by the acoustic time step. Dual-time stepping can remove this restriction, but slow inner convergence has limited reported wall-clock gains to roughly four- to tenfold, without a mechanism that strengthens under refinement. We present the macroscopic-shadow-corrected lattice Boltzmann method (MSC-LBM), which applies defect correction to the unsplit kinetic residual. At each Fourier wavenumber, a Stokes-like system is solved exactly for the conserved residual moments while preserving the second-order backward-differentiation formula (BDF2) fixed point and temporal accuracy. Along the two-relaxation-time axis $\Lambda=1/4$, $\omega^+=\omega^-=1$ makes one collision eliminate all non-hydrodynamic perturbations; $\eta(\nu)=|6\nu-1|/(6\nu+1)$ vanishes at $\nu=1/6$, leaving hydrodynamic slow modes whose long-wave shadow is inverted by the corrector. In completed timing campaigns, MSC-LBM achieves 27.95- and 52.71-fold wall-clock speedups at $N=256$, $\mathrm{Re}=1$ under matched-accuracy and $1\%$ gates. The matched speedup rises monotonically to 122.02 at $N=1024$, with gains persisting across $\mathrm{Re}=10^{-4}$–$100$. The three-dimensional D3Q19 extension reaches 27.85 and 8.84 under the same gates at $N=128$, $\mathrm{Re}=1$. A bounce-back-consistent kinetic coarse solver with adaptive relinearisation extends MSC-LBM to fully enclosed cavities, yielding setup-excluding physical-march speedups of 26.86 at $\mathrm{Re}=10$ and 3.84 at $\mathrm{Re}=100$ for $N=128$. An exact per-wavenumber symbol inverse establishes the attainable off-design contraction envelope. The contraction, parameter-sweep, and timing results jointly delimit the demonstrated operating regime: low-to-moderate-Reynolds-number transients for which acoustic stepping need not dictate computational cost.

Earth's lower mantle hosts large low shear-velocity provinces (LLSVPs) beneath Africa and the Pacific, but their origin remains debated. Davemaoite, a major lower-mantle mineral, has been proposed as a contributor to these anomalies; however, its experimentally measured shear modulus is approximately 30% lower than first-principles predictions, a longstanding discrepancy between theory and experiment. Here, using million-atom machine-learning molecular dynamics simulations, we show that nanoscale grain-boundary disorder, invisible to conventional x-ray diffraction, can reconcile this gap. Introducing 5.8 vol% disordered regions reduces the shear and bulk moduli by 37% and 12%, respectively. This shear-selective softening can reproduce LLSVP-like seismic anomalies in basalt-rich assemblages, but only at disorder levels corresponding to grain sizes smaller than those expected in the lower mantle. Combining the disorder-elasticity relationship with seismic constraints, we derive a minimum davemaoite grain size of approximately 100 nm. This limit is consistent with mantle grain growth models and with attenuation-based evidence for coarse-grained LLSVPs.

In fixed-source Monte Carlo particle-transport problems, obtaining an acceptable variance on the sought response is of paramount importance. Currently, the only rigorous tool to analyze the variance of such games is the framework of the moment equations, which is unfortunately unwieldy to use in practice. In this paper, we establish a formula that enables expressing the variance of a Monte Carlo simulation as a sum of 'variance contributions' collected throughout the underlying particletransport process. This formula applies to both direct and adjoint games, and provides a new tool to pinpoint the variance-inducing mechanisms in Monte Carlo simulations, understand common variance-reduction techniques, and even conceive new ones. We showcase the use of the variance-decomposition formula on several applications. In particular, we revisit zero-variance schemes and analyze existing variance-reduction techniques, underlining their strengths and weaknesses. A few relevant numerical examples substantiate our theoretical findings.

Gausslet orbitals on a cubic lattice, introduced in [Steven R. White, J. Chem. Phys. 147, 244102 (2017)], represent a localized, smooth, and systematically refineable basis set featuring a "diagonal" approximability of the electron repulsion integral tensor. The present work develops efficient algorithms for evaluating overlap integrals required for electronic-structure simulations using these Gausslets, specifically the kinetic, nuclear, and electron repulsion integrals (both with and without the diagonal approximation). Computational efficiency improvements rest on the exploitation of translation, permutation, and octahedral symmetries, a reordering and precomputation of nested sums, as well as an early truncation of small coefficients. Our algorithms reduce the number of electron repulsion integrals on a $5 \times 5 \times 5$ grid from $125^4 = 244140625$ to $138187$ due to symmetries, and achieve a wall-clock runtime for evaluating the remaining integrals with a truncation tolerance of $10^{-5}$ in under 1 second on a laptop computer. We apply the developed methodology to compute the ground state of the hydrogen atom and molecule as a demonstration.

Accurate long-time prediction of Cahn-Hilliard dynamics, particularly in the late-stage coarsening regime, is computationally demanding because it requires simulations over extended physical time. Various machine learning approaches have been developed to speed up the simulations through efficient surrogate evaluations. While the learned map can be accurate over a single learned interval, small prediction error will accumulate when the map is repeatedly applied in a long-time prediction, leading to large numerical errors. In this work, we identify the structure function for the coarsening dynamics of the Cahn-Hilliard equation and propose a two-stage structure-preserving Fourier Neural Operator (FNO) learning framework. Stage~I introduces a bound-conforming FNO (bcFNO) with a soft, magnitude-dependent amplitude penalty. Stage~II retains the bound-conforming objective and supplements it with structure-function regularization, yielding a structure-preserving FNO (spFNO) that corrects late-stage coarsening statistics. Numerical experiments on the various spatial resolutions and reference-solver time steps suggest that the bcFNO stabilizes the long-time prediction by suppressing large amplitude excursions, without changing the low online cost of the FNO. Furthermore, the structure-function term consistently reduces discrepancies in the normalized structure function and provides further correction in associated late-stage coarsening quantities, including the $L^3(t)$ growth trend and the scaling collapse. These results indicate that our multi-stage implementation of physical principles in operator learning could be applied to other multiscale systems with statistical scaling, such as the functionalized Cahn–Hilliard equations and turbulence, provided that the state constraints and statistical observables are adapted to the governing dynamics.

Partially miscible ganglia trapped within porous media, spanning one or multiple pores and evolving through diffusive mass transfer, are common in subsurface (e.g., CO$_2$ and H$_2$ storage) and manufacturing (e.g., fuel cells) applications. We present the ganglion network model (GNM), a reduced-order method for simulating how a population of such ganglia evolves inside an arbitrary porous microstructure. GNM operates on a tree graph, the ganglion network, extracted from the pore-scale image of a given sample. Each point on the graph encodes a possible ganglion configuration in the void space, without any loss of geometric or topological complexity. The evolution of a population is modeled by representing each ganglion as a particle on the graph and tracking it according to a set of rules formulated herein. The rules capture capillary events such as pore invasion, retraction, snap-off, fragmentation, and merger. Unlike pore-network models, another graph-based modeling tool at the pore scale, GNM solves no system of equations and its cost scales with ganglion count, not domain size. We validate GNM against an image-based pore-network model in 2.5D and 3D domains with populations undergoing ripening, dissolution, and growth. We find good agreement in ganglion statistics, aggregate properties, and spatial configuration. We further argue that the ganglion network is the statistical space needed for extending kinetic theories of Ostwald ripening from single- to multi-pore ganglia, and provide an outline for how to do this. GNM opens the door to modeling other dynamics of trapped phases in porous media.

Discussions of the sign problem usually focus on its role as the primary limitation of the applicability of quantum Monte Carlo methods for reliably solving quantum many-body models. Its behavior as a function of spatial lattice size, doping, temperature, and interaction strengths has been carefully characterized within the determinant quantum Monte Carlo algorithm. Here, we consider the much less well-explored question of the temperature- and size-dependence of the sign problem in the world-line quantum Monte Carlo method. We review how, in this algorithm, even sampling the partition function of free fermions yields a sign problem, which we connect to non-analytic behavior in the corresponding reference model, hard-core bosons on a lattice. The latter exhibits a finite-temperature Kosterlitz-Thouless phase transition in 2D, and we show how this non-analytic behavior is imprinted on the average sign of the weights for the corresponding free fermion 2D tight-binding model. Our results thus show how phase-transition information can remain encoded in the sign problem even after direct sign sampling becomes impractical.