Papers

Snow is a complex geomaterial whose macroscopic response is governed by density, temperature, and the topology of its evolving microstructure. Its mechanical behavior spans elastic, plastic, viscous, and failure dominated regimes, imposing significant challenges for numerical methods, which intends to simulate large deformations, evolving free surfaces, and complex boundary interactions. This work presents the first integration of a Modified Cam-Clay constitutive formulation for snow into a purely meshfree strong-form collocation framework based on the Generalized Finite Difference Method. The main methodological contribution is a numerical coupling that combines a global implicit mixed formulation for pressure and velocity with a constitutive return-mapping algorithm. The hydrostatic pressure contribution is obtained from a Poisson equation and subsequently corrected through the Modified Cam-Clay return-mapping procedure, while the deviatoric response is treated semi-implicitly using a numerical viscosity formulation. This partitioned treatment of the volumetric and deviatoric stress contributions enables stable simulations with comparatively large time steps while producing smooth spatial pressure fields. As a result, forces on complex boundary geometries can be evaluated accurately. Numerical results of this coupling illustrate the algorithmic stability of the framework, the effective imposition of boundary conditions, and the suitability of local spatial refinement. The feasibility of applying the framework to vehicle-snow interaction through rigid-body coupling is also illustrated. The presented formulation provides a robust basis for future simulations of dynamic snow loading on vehicle structures.

The nuclear drip line plays a crucial role in determining the composition of matter under extreme astrophysical conditions. In core-collapse supernovae and neutron-star crusts, matter is driven far from saturation density and nuclear stability; nuclei coexist with a sea of free neutrons, an effect that is present even at zero temperature in neutron-star crusts and becomes more pronounced in the hotter, neutron-rich supernova environment. This makes a careful treatment of drip-line physics essential for a realistic description of the equation of state and composition. In this work, the influence of the nuclear drip line on the baryonic composition of supernova matter is investigated within the framework of nuclear statistical equilibrium (NSE). The composition is evaluated in terms of free nucleons, light clusters, and heavy nuclei at finite temperature and global sub-saturation densities. The results indicate that, at low proton fractions and higher densities, the inclusion of nuclei beyond the drip line enhances the formation of extremely neutron-rich light clusters, leading to a significant reduction in the free-neutron density and the charge fraction of heavy nuclei. These findings demonstrate that drip-line physics has a significant impact on the composition of supernova matter and should be carefully incorporated in supernova modeling and nucleosynthesis studies.

We study the reconstruction of an unknown dynamical system from a single noisy scalar time series. The goal is to recover the underlying dynamics for forecasting. We introduce a method that uses differential embedding coordinates to identify a rational closure of the embedding dynamics directly from data. The closure is identified through a weak-form regression pipeline, which avoids unstable pointwise differentiation of noisy data. When applied to noise-free Lorenz and Rössler systems, the method recovers closures that support long forecasts across a broad ensemble of realizations ($18.1$ and $7.1$ Lyapunov times respectively). Under $15$–$30\%$ additive Gaussian noise, performance becomes system-dependent. For the Lorenz system, forecast horizons remain short even in the best cases, whereas the Rössler system generally performs better in absolute terms, though not once normalized by the Lyapunov time. Our proposed method recovers directly interpretable closure coefficients which we compared against the known analytic closures of the Lorenz and Rössler systems.

Spatially resolved EEDFs/IEDFs provide essential kinetic information about low-temperature plasmas (LTPs) and play a central role in determining transport, chemical reaction rates, and plasma surface interactions. While kinetic simulations directly resolve these distributions, experimental measurements remain challenging and are often invasive, spatially limited, or require assumptions regarding the distribution shape such as a Maxwellian. However, several macroscopic plasma observables can be measured non-invasively using advanced diagnostic techniques, providing spatially resolved information about the plasma state. An important inverse problem is therefore whether readily measurable macroscopic plasma quantities contain sufficient information to reconstruct the underlying kinetic state. In this work, we investigate this problem by learning a nonlinear mapping from spatially resolved macroscopic plasma observables to the corresponding spatially resolved EEDFs/IEDFs using a deep learning framework. Paired datasets comprising 2D macroscopic observables and spatially resolved EDFs are generated using 2D-3V PIC-MCC simulations. Three representative learning paradigms, a U-Net, a FNO, and a MeshGraphNet, are employed in this study to learn this inverse mapping. The predicted EDFs reproduce both bulk plasma and sheath characteristics with good agreement to the PIC-MCC reference data, with the FNO providing the best overall performance. Beyond conventional metrics, physics-based validation demonstrates that the reconstructed EDFs accurately recover the corresponding density and temperature, and rate coefficients. These results demonstrate that macroscopic plasma observables encode sufficient information to infer important kinetic properties in LTPs, providing a potential foundation for surrogate kinetic modeling and next-generation plasma diagnostics.

We construct the natural generalization of stochastic quantization (in the Markovian sense) by considering jump-diffusion processes. This class of stochastic processes exhibits non-continuous paths, so-called Lévy flights. In the presence of jumps, action landscapes with barriers can be efficiently explored, improving and even restoring ergodicity where traditional diffusion approaches become inefficient. We explore different strategies for constructing efficient jump updates, which we deploy to address the benchmark problem of topological freezing in 2d U(1) gauge theory.

We investigate the structural and optoelectronic properties of oligothiophene (nT) molecules encapsulated in boron-nitride (BN) nanotubes using density functional theory, many-body GW and Bethe-Salpeter equation approaches. We show that the binding energy is maximized for tube diameters of approximately 10 Å, decreasing gradually for larger diameters. The nT molecules can slide along the tube with a corrugation potential smaller than room temperature thermal energy, promoting their head-to-tail aggregation. Regarding the electronic properties, we find that hybridization with the tube electronic states has a much smaller effect than that of screening, which can close the nT photoemission gap by as much as an eV. A simple model for the polarization of the BN tube demonstrates how these polarization effects decrease with increasing nT molecule length and BN tube diameter. Compared to the gas phase optical properties, structural relaxation, hybridization, and screening can redshift the absorption onset by up to 250 meV for isolated intercalated nTs. Additionally, the study of a head-to-tail sexithiophene dimer reveals an exciton-exciton interaction that splits the absorption onset into a lowest bright exciton, separated by approximately 75 meV from a dark peak. This suggests possible collective effects upon the formation of nT chains inside BN tubes. Our results confirm and clarify experimental data, mitigating the conclusion that insulating BN tubes just act as a protecting environment for encapsulated molecules.

Bulk materials, as opposed to nanomaterials, require molecular dynamics (MD) simulations on a large spatial scale (~10^9 atoms or more) to adequately capture their atomic-scale physical properties. Previously, the introduction of machine-learning interatomic potentials (MLIPs) has extended MD to this scale, but even single-component bulk systems require tens of thousands of GPUs on high-end supercomputers. However, multi-component bulk MD simulations remain barely achievable, as the HBM footprint of existing MLIPs - already substantial for single-component systems - grows explosively in multi-component scenarios. This paper proposes an MLIP with a small HBM footprint - less than 3% that of existing MLIPs - unlocking multi-component bulk MD using only hundreds of GPUs. This is achieved by first identifying feature vectors and intermediate tensors as the two primary contributors to HBM footprints in existing MLIPs. To address these two sources, the dimensionality of the feature vectors has been reduced by introducing physical and chemical knowledge, and intermediate tensors have been eliminated by aggressively fusing all kernels into a single mega-kernel. In evaluation, the proposed MLIP has used 144 NVIDIA A100 GPUs to perform MD simulations on a 6-component bulk system with 1.14x10^9 atoms, while previously such MD simulation spatial scale has been restricted to unary systems and typically achieved on high-end supercomputers equipped with tens of thousands of GPUs.

Decoupling heat and charge transport is a key challenge in thermoelectrics. Here, we identify a route to spatially separate phonon and carrier transport in quasi-one-dimensional materials through high-throughput screening of the Materials Project database. Representative Sn$_2$S$_3$ and SbTeI exhibit a strong-intrachain–weak-interchain bonding hierarchy that favors phonon propagation along the chains while suppressing transverse lattice heat transport. In contrast, transverse valence-band states provide effective interchain electronic coupling and relatively light hole transport. This mismatch between lattice and electronic transport dimensionalities produces an inverted thermal–electrical anisotropy. Across the screened candidates, interchain lattice thermal conductivity is strongly suppressed, whereas hole transport remains weakly anisotropic or even favors the interchain direction. For SbTeI, this decoupling yields a maximum $zT$ of approximately 2.1 near 900~K. These results establish dimensionality mismatch as a general strategy for decoupling phonon and carrier transport in thermoelectric materials.

Understanding the stochastic evolution in complex networks is a central challenge across physics, biology, engineering, social science, and finance. The most consequential macroscopic events, like cascading failures in communication networks, widespread epidemic outbreaks, and rapid shifts in societal opinions, often emerge from a confluence of rare, localized stochastic processes and need to pass certain bottlenecks. Standard forward-time simulation algorithms like the Gillespie method are inefficient for the investigation of such phenomena due to catastrophic rejection rates. Advanced rare-event techniques like splitting methods and transition-path sampling often suffer from kinetic trapping, path degeneracy, genealogical correlations, or critical slowing down when applied to complex heterogeneous networks. We propose to overcome this challenge by establishing a novel technique called conditional-path Monte Carlo (CPMC), inspired by loop algorithms from equilibrium condensed-matter physics. By employing non-local updates on spacetime clusters without rejections, CPMC generates a Markov chain of trajectories that all strictly respect the targeted macroscopic boundary conditions like the occurrence of a massive network failure. We demonstrate the framework's potential by performing a simple risk factor analysis for rare large-scale epidemic outbreaks in SIS dynamics on kinship networks.

Largest Lyapunov exponents (LLEs) quantify exponential sensitivity, but data-driven estimates are often obtained from finite-precision trajectories. We show that increasing the length of a single reduced-precision chaotic record can eventually degrade a forecast-error LLE estimate. Using the logistic map at r=4, an ESP32 single-precision trajectory is reproduced bit-for-bit by NumPy float32. Across 10,000 random float32 initial conditions, every trajectory reaches an exact recurrence before iteration 7612. For one long float32 record, the estimated LLE changes from 0.6853 at N=15,000 to 0.1827 at N=20,000 and approximately zero at N=30,000 as exact train-test histories saturate. At N=100,000, the long float32 record gives 0.0016, whereas independently restarted length-100 trajectories give 0.6917; matched float64 controls remain near ln(2)=0.6931. The collapse is reproduced for 28 representative initial conditions, and its onset is strongly correlated with the recurrence scale set by transient length and digital period (Pearson r=0.982). Thus, finite-state recurrence can turn additional samples into duplicate futures rather than new dynamical information, while independent restarts substantially delay this saturation.

Neural autoregressive models have rapidly emerged as powerful emulators of high-dimensional chaotic systems, yet their long-term instability and error growth remain poorly understood, leading to ad-hoc solutions. Here, we develop an eigenanalysis framework that reveals the dynamical origin of this error growth. By analyzing the Jacobian of the learned one-step update map with respect to the state, we show how inference-time error growth, and thus model stability, is governed by its spectral radius. Direct-step architectures (models that predict the next state from the previous one) generically admit unstable eigenvalues with magnitudes exceeding one, explaining the rapid divergence of these widely used models. In contrast, integration-constrained models (where the time derivative is estimated and integrated with a higher-order integrator) collapse their eigenspectrum onto the unit circle, yielding neutral stability and a universal linear error-scaling law. The largest eigenvalue of this Jacobian provides an architecture-agnostic, a priori diagnostic of short-term skill, long-term stability, and spectral bias, without requiring an expensive rollout. Leveraging this theory, we introduce a stability-promoting loss that explicitly regularizes Jacobian-driven error amplification, improving both forecast accuracy and dynamical robustness. Demonstrated across $29$ models spanning two architectures, several explicit and implicit integrators, and multiple loss functions on the Kuramoto-Sivashinsky system, our results establish a theoretical foundation for the design and evaluation of neural emulators of chaotic multi-scale dynamics. More broadly, our framework is a step toward the kind of a priori stability analysis that numerical analysis provides for discretizations of differential equations and that scientific machine learning currently lacks.

Adjoint methods provide a powerful route for gradient-based optimization, but their physical implementation is obstructed in generic nonlinear systems because the adjoint dynamics requires backward-time evolution, Jacobian transposition, and terminal-value constraints. Here we show that nonlinear parity-time ($\mathcal{PT}$)-symmetric systems overcome this obstruction. Using a class of nonlinear non-Hermitian resonator networks, we establish symmetry relations that map the formal adjoint dynamics onto experimentally accessible forward-time evolutions supplemented by controlled injections. This construction enables exact in-situ evaluation of adjoint gradients without requiring explicit backward propagation or matrix transposition. We demonstrate the approach in nonlinear $\mathcal{PT}$-symmetric resonator chains, where the resulting optimization protocol autonomously discovers parameter configurations that realize prescribed spatio-temporal functionalities, including uniform energy redistribution and targeted wave transport at predefined time windows. Our results identify $\mathcal{PT}$ symmetry as a resource for implementing computational sensitivities within physical systems and establish a route toward self-optimizing nonlinear machines.

Thermodynamic stability and causality impose fundamental constraints on the equation of state (EoS) of nuclear matter. Verifying these constraints conventionally requires calculating quantities such as the specific heat, baryon-number susceptibility, and speed of sound, which can become computationally expensive when many candidate EoSs must be examined. We investigate whether the normalized pressure surface, $Q(T,\mu_B)=P(T,\mu_B)/T^4$, alone contains sufficient information to determine the physical admissibility of an EoS. We develop a supervised convolutional neural network (CNN) that uses only this pressure representation to classify EoSs as physically admissible or inadmissible. The network is provided with training labels obtained from direct thermodynmaic stability and causality check and its does not get any information about the parameters of the underlying EoS framework. For EoSs generated within an Ising-mapping framework, the model achieves $97.65%$ accuracy on unseen test data. Applied independently to EoSs from a distinct holographic framework, it achieves perfect classification of the test set. These results show that pressure surfaces contain geometric signatures of thermodynamic stability and causality violations that can be learned directly by a CNN. Because the classifier relies only on the pressure surface, it avoids evaluating higher-order thermodynamic observables during inference and is largely independent of the EoS-generation framework. When the pressure surface is supplied as a two-dimensional array, the machine-learning validation is approximately 20 times faster than direct validation. Our results establish a fast, framework-independent approach for identifying physically admissible EoSs directly from their pressure surfaces.

A major challenge in the field of heterogeneous catalysis is selecting an optimal catalyst support structure. Commercially available structures can be easily manufactured at scale, but their stochastic nature makes their chemical and transport properties suboptimal. This is particularly relevant for gas-generation reactions, where non-uniformity of a porous structure leads to bubble trapping. Such trapping impedes the flow of reactants to catalyst sites, leading to conversion inefficiencies. Previous experimental work demonstrated that spinodally-derived architectures, in particular bicontinuous interfacially jammed emulsion gels (bijels), can alleviate these issues and deliver superior performance. However, to the best of our knowledge, numerical studies to optimize the operating conditions for such a morphology have not been performed yet. In this work, we aim to close this gap using color-gradient lattice Boltzmann simulations of reactive flows with a novel central moments collision operator. We develop an analytical model to predict catalyst performance based on our simulation data. Our findings show that this type of morphology can achieve very high conversion efficiencies. Moreover, we demonstrate that its catalyst performance can be optimized using superhydrophilic surface coatings.

Evolving crack fields in structural health monitoring and fracture assessment must often be inferred from sparse mechanical measurements rather than dense full-field observations. We develop CNN2D–ConvGRU, a convolutional-recurrent framework for measurement-conditioned reconstruction of time-dependent phase-field brittle fracture. At each load step, the model maps a fixed-length history of phase and displacement fields and sparse current-step displacement measurements to the current full-field state. New measurements are assimilated during sequential deployment, making the framework a state-inference surrogate rather than an autonomous time integrator. It reproduces crack paths, damage evolution, and bulk displacement response, with the largest errors near propagating crack tips and steep displacement gradients and some drift at late stages. Without retraining, weights learned on a $256 \times 256$ raster are evaluated on a $512 \times 512$ raster of the same physical domain and finite-element discretization using a proportionally refined measurement grid. This empirical raster-and-sensing transfer preserves the principal damage topology and global damage evolution, although fine-scale displacement errors increase near crack tips. Comparisons with alternative spatial and temporal architectures show that CNN2D–ConvGRU offers a favorable balance between reconstruction accuracy and computational cost. Relative to repeated finite-element solutions, sequential reconstruction achieves mean speedups of $175\times$ on CPU and $253\times$ on GPU. These results demonstrate efficient full-field fracture-state reconstruction from sparse observations while retaining the spatial structure and history dependence of phase-field fracture.

Trefftz functions satisfy the differential equation locally and exactly; quasi-Trefftz functions do so approximately to prescribed high order. This paper considers (quasi-)Trefftz approximations for two-dimensional scalar electromagnetic problems. Established discretizations include the Flexible Local Approximation MEthod (FLAME), Trefftz elements ($T$-elements), and Trefftz discontinuous Galerkin (Trefftz-DG) methods. New developments are gradient-enriched FLAME (GEFLAME), conforming Trefftz–FLAME finite elements (TFF), and a full-field Bloch-wavevector-vs-frequency solution with Schur–DtN reduction. Applications cover electrostatics, scattering, singular fields, and Bloch waves in periodic structures. These methods make different compromises between conformity and flexibility. FLAME and GEFLAME incorporate Trefftz functions directly into local difference schemes; TFF uses elementwise FLAME schemes to lift polynomial traces into element interiors; $T$-elements match elementwise Trefftz spaces weakly to such traces; and Trefftz-DG couples broken Trefftz spaces through fluxes and penalties. In the reported wave-scattering examples, GEFLAME gives field and gradient errors several orders of magnitude below those of the quadratic finite-element discretization at comparable algebraic cost. Localized singular-function enrichment removes the dominant reentrant-corner error. Conforming TFF and quasi-conforming $T$-elements admit standard finite-element assembly but expose polynomial edge traces as an accuracy bottleneck. For periodic media, the bilinear Trefftz-DG formulation gives an analytic polynomial wavenumber-vs-frequency eigenproblem without restricting the Bloch multiplier to the unit circle. Schur–DtN reduction yields compact boundary-response problems.

Parameter identification of stochastic dynamical systems driven by mixed noises is challenging due to intractable likelihood functions. We propose PENN-GMD, a parameter estimation neural network that maps partially observed trajectories to a Gaussian mixture distribution (GMD) over the system parameters. Unlike conventional uncertainty estimates, the GMD employs full covariance matrices to explicitly reveal parameter couplings and multi-modal likelihood structures. The network is trained by minimizing the negative log-likelihood via a surjective parameterization that hard-encodes all GMD constraints, thereby approximating the true likelihood. We validate the method on five numerical examples with increasing complexity, including systems driven by fractional Gaussian and Lévy noises, oscillators with colored noise, coupled neurons under different observability, and an aeroelastic airfoil with unidentifiable stochastic disturbances. Results demonstrate that PENN-GMD accurately recovers likelihood distributions, captures parameter couplings, and naturally diagnoses non-identifiability through variance broadening or mode splitting. These capabilities establish PENN-GMD as a practical tool for uncertainty-aware parameter identification in complex stochastic systems where conventional likelihood-based methods are infeasible.

Moiré superlattices reshape lattice dynamics without altering chemical composition, yet how crystal anisotropy modifies this control remains unclear. We combine density-functional-theory (DFT)-calibrated lattice-dynamical calculations with angle-matched untwisted controls to study puckered bilayer SnSe across seven commensurate twist angles ($3.18^\circ$–$8.77^\circ$). At 300 K, twisting suppresses the band-path heat-capacity-weighted mean-square group velocity to 2.6–8.4% of the control values; the suppression spans a broad frequency range rather than a few soft branches. The velocity response crosses over between $4.78^\circ$ and $3.82^\circ$ into a regime where the relaxed stacking textures and frequency-resolved velocity profiles become self-similar, with the normalized mean-square velocity ratio spanning only 11.1% of its mean across the three smallest angles—a finite anisotropic crossover, not a singular-angle condition. Direct DFT–MACE force-constant agreement ($r=0.996$), uniform $4\times4\times1$ stability scans, and acoustic-sum-rule and path-density tests support the trend. The equilibrium trend is defined by six structures after excluding one relaxation-sensitive case. These results extend phonon twistronics to low-symmetry layered materials and identify crystal anisotropy as a key determinant of finite-angle phonon crossover behavior.

Elastoplastic analysis is computationally demanding because its nonlinear, path-dependent constitutive behavior requires incremental loading and repeated iterative solutions. To address this challenge, we propose Plasolver, a physics-informed neural operator framework that combines the efficiency of operator learning with the accuracy and robustness of classical numerical solvers. Plasolver consists of a physics-informed pretraining stage and an optional warm-start stage. During pretraining, the neural operator is trained solely by minimizing the incremental potential energy of elastoplasticity formulated by Simo, without requiring any labeled solution data. It operates directly on unstructured point clouds by encoding spatial coordinates, loading histories, and material properties as unified point-wise prompts. This formulation provides dual invariance to spatial and loading-path discretizations, enabling consistent predictions across different spatial resolutions and different numbers of increments representing the same loading trajectory. The pretrained Plasolver achieves relative errors on the order of 1% while providing approximately two orders of magnitude acceleration over conventional finite element simulations. In the warm-start stage, the pretrained prediction is supplied as the initial solution to a classical iterative solver, preserving its numerical accuracy, robustness, and convergence properties while substantially accelerating convergence. Numerical results show that Plasolver reduces the required number of iterations by approximately 50% compared with conventional zero-initialized solvers and converges to solutions at any prescribed tolerance. Plasolver thus provides an efficient, accurate, and discretization-invariant computational framework for nonlinear, path-dependent elastoplastic problems.

Atom-light interaction noise can limit thermal-vapor sensing. Existing theories often treat internal-state dynamics, finite-mode atomic motion, and stochastic renewal separately, obscuring their coupled contributions to measured noise. We develop a general stochastic Liouville-transport theory, tested against polarization-resolved resonant Cs D$_2$ spectra. Joint experiment-theory analysis identifies atom-light noise below approximately 100 kHz as transit-dominated. Ballistic motion through the finite Gaussian mode modulates both the coupling-weighted effective atom number and trajectory-dependent Rabi coupling, producing predominantly common-mode noise. Boundary renewal introduces atoms with independently sampled ground-state sublevels, generating differential population fluctuations with opposite effects on the circular channels. Under an applied longitudinal magnetic field, experiment and theory show the same qualitative nonmonotonic change in common-mode suppression, supporting Zeeman redistribution of the channel responses. The framework can analyze noise in other thermal-atom sensors, including Rydberg-atom electric-field measurements.

Updated reference data and an analytic parameterisation of elastic electron and X-ray scattering are presented for all 118 neutral atoms. The reference electron densities for the multi-electron elements $Z=2$–$118$ are computed with the relativistic B-spline Dirac–Fock code atomx, while hydrogen is constructed from the exact relativistic one-electron Dirac $1s$ solution; the electron scattering factor $f_e(g)$, X-ray scattering factor $f_x(g)$ and radial moments are derived from these densities. The reported parameterisation extends the fixed-size Lobato–Van Dyck hydrogenic expansion while retaining closed-form expressions for $f_x(g)$, $\rho(r)$, the electrostatic potential $V(r)$ and the projected potential $V(R)$. These extensions are an element-adaptive basis size $n_t(Z)$, a simultaneous real- and reciprocal-space fit, an exact $\langle r^4\rangle$ constraint in place of the non-relativistic Kato cusp, and a charge-carrying Dirac–Padé basis term that adds a polynomial-times-exponential shape channel without replacing the hydrogenic basis by a tabulated Dirac radial function or assigning the term to a physical shell. For the same total parameter count, the Dirac–Padé-enriched basis improves on the parameter-matched non-relativistic basis for 111 of the 118 elements, lowering the mean total cost by 39%. Relative to a controlled fixed five-term refit on the same reference grid and objective, the element-adaptive bases improve the median reciprocal-space deviations by about three to four orders of magnitude and resolve shell structure in $4\pi r^2\rho(r)$ that the fixed five-term basis cannot. The largest changes occur near the nucleus and in the reciprocal-space tail beyond the legacy 12 inverse angstroms range, which is directly relevant to quantitative high-angle scattering and electron-diffraction measurements.

ESPResSo++ is an open-source software package for molecular dynamics (MD) simulations with a particular emphasis on coarse-grained (CG) models of soft matter systems. Written in C++ with a flexible Python interface, it is designed for high-performance computing (HPC) environments and supports massively parallel simulations through MPI. The package enables simulations of polymers, membranes, colloids and complex fluids with a wide range of interaction models and advanced algorithms.

We implemented on GPU a mixed quantum-classical Liouville molecular dynamics simulation based on a momentum-jump-free theory. The trajectory spawning that was previously implemented on CPU for sampling enhancement was eliminated to avoid the overhead of thread divergence and dynamic memory allocation on the GPU. This achieved a speedup of an order of magnitude compared to the CPU computation with spawning, as well as a linear scaling with respect to the number of sampling trajectories.

Effective Hamiltonians connect microscopic interactions to measurable collective behavior in quantum materials, but determining their parameters directly from experiment remains a challenging inverse problem. We introduce a supervised machine-learning framework that infers Hamiltonian parameters from nonlinear terahertz two-dimensional coherent spectra. A calibrated forward model generates spectra from candidate Hamiltonians, a common preprocessing pipeline maps simulated and experimental spectra into the same representation, and a neural network learns the inverse map from spectral fingerprints to microscopic parameters. We demonstrate the approach for rare-earth orthoferrites using a two-sublattice Landau–Lifshitz–Gilbert spin model with exchange, Dzyaloshinskii–Moriya interaction, anisotropies, and damping. Synthetic benchmarks show that nonlinear spectra encode parameters beyond those fixed by the linear response, with inference accuracy tracking the physical spectral sensitivity and robustness against noise improved by using multiple inter-pulse delays. Applied to experimental THz-2DCS data from Sm$_{0.4}$Er$_{0.6}$FeO$_3$, the inferred parameters yield physically reasonable forward simulations, while remaining discrepancies identify limitations of the reduced model. These results establish THz-2DCS as a data-rich platform for effective-Hamiltonian inference and model refinement, enabling experimentally driven identification of microscopic interactions while providing a foundation for understanding, predicting, and ultimately controlling the emergent properties of quantum materials.

Quantifying forced and internal variability in climate data is fundamental to detecting the climate change signal and assessing uncertainty in climate projections. We propose a metric that quantifies the relative magnitude of internal variability in single-model initial-condition large ensembles (SMILEs). We measure forced and internal variability by two 1-Wasserstein distances, a form of optimal transport cost. The proposed metric is the distance for internal variability divided by the sum of the two. Computing this ratio requires only sorting the samples and differencing the resulting quantiles, with no additional parameters. The metric reflects the entire distribution shape and applies to non-Gaussian variables, because the 1-Wasserstein distance quantifies the difference between any two probability distributions. We validate the metric with synthetic climate data from Gaussian, uniform, and lognormal distributions. Unlike the two existing metrics, the proposed metric gives stable estimates irrespective of the distribution shape and the presence of outliers, provided that the ensemble has about 40 members or more. We then apply the proposed and existing metrics to the 2 m air temperature and total precipitation of the Community Earth System Model Large Ensemble (CESM-LE) under two different forcing scenarios. All the metrics indicate that the relative contribution of internal variability decreases as the forcing increases, but the proposed metric shows this response most clearly. The 1-Wasserstein distance thus provides a simple and useful tool for analyzing large-ensemble datasets.

The immersed boundary (IB) method models fluid–structure interaction using the natural Lagrangian and Eulerian formulations of structural mechanics and fluid dynamics, respectively, but explicit time discretization of the IB force imposes a stiffness-dependent upper bound on the time-step size. Treating these forces implicitly removes this restriction, but the resulting coupled linear systems become increasingly difficult to solve as the structural stiffness increases. Algebraically eliminating the Lagrangian degrees of freedom yields a reduced Eulerian velocity–pressure IB system. Here we introduce a coupling-aware Vanka (CAV) smoothing strategy to enable effective multigrid preconditioning of this system. CAV patches are built as unions of standard pressure-centered Vanka patches, with the graph of the Eulerian elasticity matrix determining which patches are combined. Under grid refinement, CAV patch sizes remain bounded, and the computational cost of each multigrid cycle thereby grows linearly with the number of Eulerian degrees of freedom. Tests using target-point, membrane, and beam force laws show that CAV-preconditioned FGMRES reduces the relative residual by ten orders of magnitude in $9--15$ iterations, with little growth under grid refinement. In a nonlinear benchmark modeling flow past a flexible fiber, the average number of FGMRES iterations per Newton solve increases only from $8.6$ to $9.5$ as the Eulerian grid is refined from $32\times32$ to $256\times256$ cells. To our knowledge, CAV provides the first robust multigrid strategy for time-dependent implicit IB formulations.

When a crystal undergoes a second-order structural phase transition, such as in ferroelectrics, Peierls, and charge-density waves, the diverging fluctuations of the order parameter lead to the break- down of the standard phonon quasiparticle picture. Simulating these highly anharmonic regimes is notoriously challenging, as methods such as molecular dynamics suffer from a critical slowdown near the transition point, while the harmonic approximation fails dramatically at saddle points of the energy landscape stabilized by quantum or thermal ionic fluctuations. This work introduces a new approach, based on the variational free-energy principle, to predict critical long-range behavior and dynamical spectra in strongly anharmonic systems, even when quantum ionic fluctuations dominate. The proposed framework builds upon the stochastic self- consistent harmonic approximation but reduces its computational scaling with the number of atoms, N, from O(N^6) to O(N^2) and the memory requirement from O(N^4) to O(N). We benchmark the method on the prototypical lead-free metal-halide perovskite CsSnI3, a promising candidate for photovoltaic engineering, simulating its phase stability and Raman spectrum near the phase transition, where the breakdown of the quasiparticle picture becomes evident. We demonstrate the effectiveness of the method by computing the full free-energy Hessian and the critical temperature in a supercell with 1080 atoms. Such simulations would have required tens of thousands of years with the legacy approach; it is now feasible in a few hours on consumer hardware.

Accurate treatment of long-range interactions in machine learning interatomic potentials (MLIPs) is essential for electrochemical simulations. However, aggregate energy and force errors alone are insufficient to establish an MLIP's physical accuracy since they do not detect qualitative inconsistencies in the model such as the prediction of image-charge attraction, dielectric screening, or charge transfer. We introduce a benchmark suite EPhEct (Electrostatic Phenomena for Electrochemistry) of focused test cases designed to evaluate MLIPs on electrochemically relevant physical phenomena. The tests probe for image-charge attraction at a metal electrode, the splitting between longitudinal and transverse optical phonons as a probe of ionic and electronic screening, the dipole moment of interfacial water, and Fermi-level pinning during ion discharge. These tests establish a qualitative diagnostic routine complementary to aggregate energy-force metrics.

Hybrid functionals are routinely employed self-consistently within the generalized Kohn-Sham framework. The evaluation of the nonlocal Fock exchange operator makes hybrid functional calculations computationally expensive, in particular with plane-wave basis sets. Here, we investigate the advantages of non-self-consistent hybrid functional calculations, focusing on the evaluation of atomic forces. The analytical force terms that arise due to non-self-consistency are computed using density functional perturbation theory (DFPT), as implemented within the Quantum ESPRESSO distribution. A non-self-consistent hybrid force calculation thus consists of self-consistent DFPT calculations with a local or semi-local functional and a single evaluation of the Fock exchange operator. The overall computational cost is, thereby, reduced in general, especially for solids that require a dense Brillouin-zone sampling. Moreover, results for structural parameters are barely affected by self-consistency, and vibrational frequencies are typically agreeing within 0.5%, thus making non-self-consistent calculations an interesting alternative.

We investigate range-separated hybrids (RSHs) and locally range-separated hybrids (LRSHs) for linear-response time-dependent density-functional theory (TDDFT) Sternheimer calculations of the photoionization spectrum of the Ne atom. This system constitutes a stringent test for approximate exchange-correlation treatments because it exhibits both valence and core resonances with very different energy scales. Building on previous work employing a simple one-parameter local range-separation function, we assess here a more flexible two-parameter range-separation function designed to improve the high-density limit. We compare photoionization spectra and resonance parameters obtained with RSHs and LRSHs. We find that the LRSH approach with the two-parameter range-separation function provides an overall satisfactory description of the photoionization spectrum, including energies for both valence and core resonances. The lifetimes of the 2s $\rightarrow$ np valence resonances are also reasonably reproduced since their decay does not involve double excitations. In contrast, the lifetimes of the 1s $\rightarrow$ np core resonances are overestimated by orders of magnitude because their Auger decay channels involve double excitations that are absent in adiabatic, single-determinant TDDFT. Obtaining more accurate resonance widths within linear-response range-separated TDDFT would require using multideterminant schemes and/or adding a frequency-dependent response kernel.