Papers

We systematically investigate finite-difference (FD) derivative computation in Physics-Informed Neural Networks (PINNs) as an alternative to automatic differentiation (AD). On three benchmark PDEs we show that, with a properly calibrated step size, FD matches AD in accuracy on every problem while running faster across the full tested batch-size range and using substantially less GPU memory, and that a stochastic variant we propose outperforms AD on a stationary problem. We further show that for neural architectures with inter-sample dependencies (e.g. BatchNorm, self-attention) the standard PyTorch autograd idiom is silently incorrect; the correct per-sample alternative is computationally infeasible at PINN-relevant batch sizes, while FD provides a forward-only approximation that is empirically an order of magnitude closer to the true per-sample derivative.

Stark acceleration and deceleration based techniques for generating particle ensembles with low velocity spread are useful in many experimental applications. For a given velocity distribution of a particle ensemble, these techniques accelerate or decelerate a small subset of the total population, with a low velocity uncertainty. However, narrowing the original velocity distribution by accelerating or decelerating the ensemble particles near the mean velocity is fundamentally limited and not yet explored. We present a numerical study of particle dynamics using neutral cesium atoms as an example. We investigate different interaction regimes, identify key limitations, and propose an interaction regime in which optical Stark deceleration can be used to narrow the velocity distribution of a propagating ensemble about its mean velocity. These findings have potential implications for optical manipulation and control, controlled collisions, and matter interferometry.

Degradation in lead halide perovskite solar cells is analysed by inverse modelling of published measurements of characteristics of a single solar cell at ages 0, 90, 280, 480 minutes. We employ machine learning to deduce distributions of material parameter values and hence the physics linked to measured changes. Bayesian parameter estimation is coupled with drift diffusion simulations using the IonMonger code combined with an optical model. We accurately replicated measured changes in device performance with age through variations in model input parameters. Our key result is that degradation is influenced by correlated changes in the concentrations and diffusion coefficients of mobile ions and by interface recombination at large mobile ion concentrations. This study demonstrates the power of machine learning combined with simulations to reliably interpret experimental results, a task which is problematic if using simulation models with only manual exploration of the input parameter space.

Computational head models are typically applied to isolated impacts, leaving repeated head loading largely unexplored. An Ogden-Roxburgh Mullins damage formulation was implemented in a high-fidelity finite element head model to represent loading-history-dependent softening during cyclic brain-tissue deformation. Repeated-loading histories derived from mixed martial arts head-impact data were applied and compared with damage-free hyperelastic (HE) and linear visco-hyperelastic (LVHE) model variants. Under five identical single-axis cycles, Mullins-type softening progressively increased strain and strain rate metrics relative to the HE model. Mullins-based injury probabilities progressively exceeded strain-based HE predictions and diverged from unchanged kinematics-based predictions, indicating that neglecting prior softening may underestimate injury risk. In a randomized twenty-cycle multiaxial sequence, cycles of similar kinematic intensity produced different deformation and injury-risk estimates depending on prior softening. HE and LVHE models predicted higher injury probabilities initially, whereas the Mullins-based model produced the largest later-cycle estimates and highest probability of at least one injury over the sequence. Regional amplification depended on loading direction and prior softening, with no direction-independent trend among brain substructures. Gyral elements exhibited higher cumulative maximum principal strain than sulcal elements, which showed greater amplification relative to initial responses. These findings demonstrate that short-term damage-mediated softening can substantially amplify tissue deformation and injury-risk estimates beyond damage-free head models under the same loading histories. Further experimental characterization of cyclic brain-tissue softening is needed to improve models of repeated head loading and traumatic brain injury.

High-order shock-capturing schemes often face severe computational bottlenecks due to wide stencils, stagewise nonlinear weights, and expensive local characteristic decompositions. To address this, we propose a compact finite-average-moment (FAM) framework for hyperbolic conservation laws that completely decouples the formal spatial order from the number of evolved local degrees of freedom. By evolving only a $\mathbb{P}^1$ moment state (i.e., cell averages and scaled first-order moments), we achieve up to sixth-order accuracy ($k=3, 4, 5, 6$) via compact, polynomially exact linear moment reconstructions. A key algorithmic innovation is consolidating the nonlinear stabilization into a single, derivative-free oscillation-elimination (OE) procedure applied only after the final Runge–Kutta stage. This OE procedure exactly preserves cell averages while applying an explicit exponential correction to the first-order moments, avoiding repeated stagewise nonlinear weights or limiting. Fourier analysis of the linear backbone demonstrates $k$th-order accuracy, strong spurious mode damping, and $(k+1)$th-order cell-average superconvergence. Extensive numerical experiments on scalar laws and the Euler equations verify the expected smooth accuracy and robust, nonoscillatory shock resolution. Notably, by performing componentwise reconstruction directly in conservative variables and streamlining stabilization, the FAM schemes deliver significantly lower complete-run wall-clock times compared to the multi-resolution weighted essentially non-oscillatory (MR-WENO), unified-stencil Hermite WENO (HWENO-U), and OE-HWENO methods.

A verbal theory does not run: translating it into an executable model requires choices about variables, interventions, and interactions. We tested whether the presentation of otherwise identical theoretical content systematically changes the programs produced by large language models. In a prospective preregistered randomized experiment, 16 independent assignment bits allocated 32 paired renderer slots between a structured intervention contract and connected prose. Two pinned LLM snapshots translated five anonymous theoretical accounts, yielding 320 preauthorized single-shot programs in a frozen sparse quadratic language. A deterministic evaluator measured atomic finite-difference responses (H1) and mixed interaction responses (H2). Two primary endpoints assessed cross-model matched-distance reduction and closed-set same-account identifiability, with exact randomization inference and 27 preregistered support criteria evaluated across signed-linear and magnitude-rank pipelines. Both H1 and H2 returned the registered verdict NOT_SUPPORTED; only 19 of 108 criterion evaluations passed. Same-account identifiability remained near chance (AUC 0.469-0.523, against a registered 0.80 threshold). One H2 matched-distance endpoint moved and survived multiplicity correction in the signed-linear pipeline, but the corresponding magnitude-rank result missed the registered effect-size floor, so it did not satisfy the joint support rule. Thus renderer format did not produce the uniform, family-invariant, classifiable behavioral geometry predicted in advance. The result places a concrete boundary on specification-format effects in LLM theory-to-program translation and provides a fully auditable randomized design, datasets, and software for studying executable formalization.

The finite element immersed boundary method (FE IBM) with a distributed Lagrange multiplier is an attractive framework for fluid structure interaction (FSI) because it couples an Eulerian description of an incompressible fluid to a Lagrangian description of an immersed solid on two independent, nonconforming meshes, avoiding the costly remeshing required by body fitted arbitrary Lagrangian Eulerian methods. When the immersed solid is modelled as a fully incompressible hyperelastic material, however, a direct finite element discretization of the deviatoric solid stress fails to enforce the Lagrangian incompressibility constraint, and the computed structure exhibits spurious volumetric instabilities and locking. In this work we present a mixed formulation of the distributed Lagrange multiplier FEIBM that restores volumetric stability. Following the theory of nearly incompressible hyperelasticity, we augment the solid stress with a volumetric contribution derived from a dilatational strain energy and introduce an additional solid pressure field, enforced weakly, that plays the role of the Lagrange multiplier for the Lagrangian incompressibility constraint J = 1. The resulting formulation is discretized in space by finite elements and in time by an unconditionally stable semi implicit scheme, and is implemented as a reusable Immersed Boundary physics kernel within the GRINS multiphysics framework, built on the libMesh finite element library. The method is verified on three FSI benchmarks, an elliptically displaced thick ring returning to equilibrium, a radially stretched incompressible ring, and a disk falling under gravity in a viscous fluid, for which the mixed formulation removes the volumetric failure observed with the unstabilized formulation and reproduces the analytical terminal velocity of the falling disk to within 1%.

The direct-sum evaluation of Madelung constants is complicated by the conditional convergence of lattice sums, which gives rise to a shape-dependent boundary term. In this work, we present, for the first time, a closed-form analytic expression for this boundary term that is valid for arbitrary crystal geometries. For general triclinic lattices, this boundary term maps exactly onto the electrostatic potential generated by a set of uniformly charged parallelograms. In addition, we demonstrate that the residual finite-size correction for a crystal of characteristic size $p$ decays as $(2p+1)^{-2}$. Building on these results, we develop a robust direct-sum method for the accurate computation of Madelung constants in arbitrary triclinic lattices and validate its effectiveness through explicit calculations on representative Bravais lattices.

Carleman lifting converts nonlinear polynomial dynamics into finite linear systems, but the resulting truncations are generally nonunitary and need not correspond to physical quantum evolution. We prove that a finite linear endpoint admits an autonomous Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) realization on vacuum coherences if and only if it is invertible and power bounded. The construction is explicit and realizes the nonunitary map directly as open-system dynamics, with no endpoint postselection and with one encoding and one decoding over repeated timesteps. We apply the result to the complete D2Q9 multiple-relaxation-time lattice Boltzmann (LB) timestep by compiling collision and periodic streaming into a single Carleman endpoint. The resulting GKSL evolution reproduces the classical Carleman trajectory over multiple timesteps, while the remaining discrepancy from nonlinear LB dynamics is the expected Carleman truncation error. The result establishes a general criterion for autonomous open-quantum realization of finite nonunitary dynamics, with Carleman–LB dynamics as a concrete example.

Graph-theoretical analysis offers a principled framework for quantifying emergent dynamics in neuronal cultures. Here, we present an analytical pipeline for inferring network-level properties of in vitro cortical cultures from multichannel electrophysiological recordings. The approach builds on a recently proposed Reservoir Computing (RC) framework (Auslender et al., 2025), which enables direct extraction of an Intrinsic Connectivity Map (ICM) from neural activity. We interpret the ICM as an effective adjacency matrix and apply graph-theoretic centrality measures to quantify node- and edge-level contributions to the culture's collective dynamics. We systematically evaluate both local and global graph metrics and examine their relationships with experimentally measured activity features, including firing rates and network-level descriptors. To validate the inference procedure, we also simulate the experimental environment, enabling controlled benchmarking of the RC-derived connectivity against a known ground-truth adjacency matrix and assessment of model performance as a function of graph structure. Our results demonstrate statistically robust associations, of varying strength, between graph-theoretic measures and experimentally observed activity patterns. These findings additionally support the validity of the RC-based connectivity inference and establish a scalable, data-driven framework for functional network characterization in neuronal culture systems.

Recently, tensor-based physics-informed neural networks (T-PINNs) have received increasing attention. However, existing T-PINNs still face a fundamental challenge: they mainly rely on pre-specified low-rank tensor decompositions with manually tuned ranks, which limits their ability to capture the underlying structures of multivariate solution functions and hinders their practical deployment. To address this challenge, we propose a hierarchical rank-evolving (abbreviated as HRE) representation for multivariate functions, which endows us to faithfully capture the underlying structure of the targeted multivariate function accompanying with automatic rank determination. Concretely, in the hierarchical design of HRE representation, the target multivariate function is decomposed as a small-scale inner tensor with a set of univariate functions along each mode, where a customized tensor network decomposition can be readily deployed to capture the underlying structure of the small-scale inner tensor. In HRE representation, the crucial hyperparameters, ranks, can be adaptively revealed during the decomposition, freeing us from manual rank tuning and making HRE practically applicable to real-world problems. Besides, we build the HRE-PINNs correspondingly. Extensive numerical experiments, including high-dimensional static problems (Helmholtz equation and Poisson equation), nonlinear time-dependent problems (Klein-Gordon equation), and complex fluid-dynamics problems (flow mixing equation and Navier-Stokes equation), demonstrate that HRE-PINNs consistently outperform existing state-of-the-art approaches in terms of accuracy.

Electromagnetic inverse scattering is a nonlinear and ill-posed computational imaging problem, where accurate reconstruction is challenging due to measurement limitations, noise, and high computational costs, especially for 3-D imaging. Although physics-driven neural networks (PDNNs) reduce the dependence on labeled training data, existing accelerated PDNN frameworks often rely on preliminary reconstruction-based region selection, which may introduce instability when the selected region is inaccurate. In this paper, a coordinate-residual physics-driven neural network (CRPDNN) is proposed for 3-D electromagnetic inverse scattering. CRPDNN represents the unknown complex contrast distribution using normalized spatial coordinates and a residual convolutional network, whose parameters are optimized by enforcing consistency between the measured and model-predicted scattered fields. Unlike existing subregion-accelerated PDNN approaches, CRPDNN does not require a preliminary reconstruction, thereby avoiding dependence on its accuracy. For the reported noise-free 3-D synthetic cases, CRPDNN achieves an average relative error of 2.10%, compared with 7.97% for CSI and 3.99% for $L_{2/3}$-FBE-WCIE, while providing approximately 5.5- and 12.1-fold speedups over the two baselines, respectively. Additional 2-D comparisons further demonstrate its stability and computational efficiency relative to existing PDNN frameworks. CRPDNN also maintains reliable reconstruction performance under noisy measurements, and the 3-D Fresnel experiments further indicate its potential for practical imaging applications.

A discrete entropy inequality is the principal nonlinear stability estimate available for systems of conservation laws, and evaluating it presupposes a physically admissible state. So far, however, the two have been secured separately. Entropy-stable schemes are almost always semi-discrete, are built around one selected entropy pair, and take for granted the positivity of density and pressure that makes the entropy well defined in the first place, whereas bound-preserving limiters keep the solution admissible but deliver no entropy estimate. For the special relativistic Euler equations, the two cannot be separated at all, since the conservative-to-primitive map is implicit, and an inadmissible state therefore has no entropy to correct. Here we construct high-order discontinuous Galerkin and finite volume schemes that, to our knowledge, for the first time, are entropy stable in the fully discrete sense for an arbitrary prescribed finite family of convex entropy pairs, a property we call multi-entropy stability, and are provably admissible wherever an entropy is evaluated. All of this is achieved by a single cellwise projection, and neither conservation nor the design order is lost. The construction rests on relativistic causality, which bounds every characteristic speed by the speed of light. Consequently, the numerical viscosity can be fixed once for all states and all equations of state, and one two-point building block then serves the whole entropy family. Since only the convexity of the admissible set and this speed bound are used, the same route remains open for related systems. Finally, in computations with four equations of state, the schemes retain high-order accuracy close to vacuum, produce no inadmissible state in strong shocks, near-vacuum shock–vortex interaction or jets with Lorentz factor above $70$, and confirm the monotone decay of every enforced discrete entropy.

The electron-phonon interaction underpins many material properties, for example, the conductivity of metals and the optoelectronic response of semiconductors. First-principles calculations of the electron-phonon interaction are a powerful tool to quantitatively describe many of these properties in increasingly complex materials. However, one key assumption of all calculations is that phonons have infinite lifetimes, an approximation that may break down when anharmonic phonon-phonon interactions are strong. In this work, we present a theory for the interaction of electrons with finite-lifetime phonons experiencing dephasing. Using a first-principles implementation of the theory, we find that anharmonic dephasing dramatically enhances electron-phonon scattering rates in metallic MgB2. Microscopically, phonon-phonon interactions create new scattering channels that increase the phase space available for electron-phonon scattering. As a result, anharmonic dephasing strongly suppresses conductivity in MgB2, bringing the calculated values substantially closer to experiment within the Boltzmann transport equation framework. This example establishes the importance of finite phonon lifetimes in the evaluation of electron-phonon scattering, and the microscopic mechanism suggests that anharmonic dephasing could play an important role in the conductivity of many metals. More broadly, our theory and first-principles implementation of anharmonic dephasing in the electron-phonon interaction provides a solid foundation to explore this regime in other materials.

Electron-phonon coupling has been a central topic in condensed matter physics for decades, and firstprinciples methods have demonstrated remarkable success in quantitatively capturing its role in a wide variety of physical phenomena and materials. Conventional calculations of electron-phonon coupling typically assume that phonons have infinite lifetimes, but phonons can exhibit finite lifetimes due to anharmonic phonon-phonon interactions. In this work, we derive an expression for the electron-phonon coupling scattering rates including the effects of anharmonic three-phonon interactions, which lead to phonon dephasing and finite phonon lifetimes. We also describe a first-principles implementation of this anharmonic electron-phonon coupling which can be seamlessly integrated within existing workflows for the evaluation of electron-phonon and phonon-phonon coupling interactions. Finally, we present calculations of electron-phonon scattering rates including phonon dephasing in a range of materials, and discuss the different microscopic mechanisms by which anharmonic phonons influence electron-phonon coupling. This study establishes the importance of finite phonon lifetimes in the evaluation of electron-phonon coupling, and provides a platform to explore these effects in a wide range of materials and phenomena.

Using first-principles calculations, we predict a new perovskite compound SrAsO$_3$ . The undoped cubic phase has pronounced soft-phonon instabilities, which are gradually suppressed upon K doping the Sr site. The cubic phase becomes dynamically stable for K-doping levels above approximately 60%, and the stabilized K-doped phases are metallic with predicted conventional phonon mediated superconductivity. Moreover, the inclusion of nonlocal exchange interactions broadens the electronic bandwidth, enhances the electron-phonon coupling (EPC) strength, and increases the superconducting transition temperature ($T_c$) of these doped compounds. In particular, the HSE06 hybrid exchange-correlation functional corrected EPC constant $\lambda$ reaches 1.41 for Sr$_{0.4}$K$_{0.6}$AsO$_3$, corresponding to a predicted $T_c$ of 44.3 K. These results suggest that SrAsO$_3$ is a BaBiO$_3$-like superconducting perovskite driven by strong electron-phonon coupling.

Laser wakefield acceleration promises compact electron accelerators for applications in medicine, industry, and fundamental science. Yet, despite rapid progress, accurately predicting the electron energy attainable in a given experimental configuration and the acceleration length required to reach it remains an open challenge. Here we use Bayesian optimization combined with advanced particle-in-cell simulation techniques to determine the maximum electron energy that a self-guided laser wakefield accelerator driven by a laser of a given energy and wavelength can produce. By systematically optimizing the accelerator performance across a range of laser energies and wavelengths, we derive energy-optimized scaling laws. These scaling laws yield the highest electron energy over the shortest acceleration length possible, are expressed solely in terms of laser energy and wavelength, and are accompanied by the complete set of laser and plasma parameters required to enable the scaling. The resulting scaling laws provide practical guidance for designing state-of-the-art laser wakefield acceleration experiments operating at their fundamental performance limits.

We propose a general framework for automatic force field refinement guided by phase diagrams, grounded in Lee-Yang phase transition theory. The central idea is to directly use the partition function modulus as a phase-diagram-guided optimization target. Evaluating the modulus at points close to the real axis, where the Lee-Yang zeros are mostly associated with the phase transition, is more physically meaningful and avoids the numerical difficulty of explicitly solving for the zeros. This approach requires no system-specific order parameters or response properties for characterizing phase transition points, making it universal across various discontinuous phase transitions and material systems. We validate the method on refining parameters of a Lennard-Jones potential covering both gas-liquid and solid-liquid transitions, and a Cu embedded-atom method potential based on experimental melting curves. The refined force fields reproduce the target phase diagrams with significant improvement across all systems. For Cu, the refinement simultaneously improves predictions of enthalpy and heat capacity, which are observables beyond the optimization target. These results establish Lee-Yang theory as a practical tool for contemporary force field development.

The vector spherical harmonics (VSH) are crucial ingredients for representation of arbitrary vector fields in curvilinear coordinate systems throughout mathematical physics. Performant, efficient computations of the VSH are essential to high-performance computing applications that utilize spectral analysis, including magnetohydrodynamics (MHD) simulations for geophysics and astrophysics. The free open-source FORTVSH library provides a catalog of numerical routines for evaluating various VSH forms and related mathematical functions. This manuscript establishes the theoretical basis, addresses the computational need for efficient VSH solvers, and presents a comprehensive overview of the FORTVSH routines.

Hypersonic vehicles and atmospheric-entry bodies can experience radio-frequency communication blackout when shock-heated gas surrounding the vehicle becomes sufficiently ionized that plasma cutoff and collisional attenuation restrict electromagnetic transmission. Magnetic-window approaches attempt to reduce this loss by exploiting the anisotropic dispersion of a magnetized plasma, in which selected right-hand or whistler-like modes may propagate along preferred directions. This paper develops a reduced-order scaling model for magnetic-window transmission through a finite-thickness hypersonic plasma sheath when the magnetic source is represented as an onboard axial dipole. The model gives a closed-form estimate of the collisionless angular aperture, compares the underlying projected-cyclotron criterion with full cold-plasma dispersion roots, and extends the aperture estimate to a loss-limited cone using a collisional optical-depth approximation. A simplified neutral-density, speed, and ionization closure is used only to generate qualitative parametric maps and sensitivity trends. The results clarify how dipole-field decay, sheath thickness, radio frequency, vehicle scale, electron density, and collisions jointly constrain the candidate transmission window. The contribution is intended as a screening framework for selecting cases for full-wave electromagnetic simulation, nonequilibrium aerothermochemistry, antenna-coupling analysis, and laboratory validation, rather than as a demonstrated engineering solution to plasma blackout.

The chemical potentials of multi-component mixtures underlie many physical and chemical phenomena, but remain challenging to compute. The S0 method enables the computation of chemical potentials from equilibrium molecular dynamics simulations, by leveraging the thermodynamic relationship between particle number fluctuations and derivatives of chemical potentials, followed by numerical integration along different compositions. Here we generalize the S0 method from two-component mixtures to neutral multi-component mixtures. We first extend the statistical mechanical formalism to high-dimensional compositional space, and then introduce a Gaussian process integration scheme combined with active learning to efficiently integrate chemical potentials and sample diverse compositions. We use this method to compute the mixing free energies of a molten metal alloy, and the solubilities of two paracetamol polymorphs in water-ethanol solvents. The extended S0 method provides a practical and scalable route for computing chemical potentials in neutral bulk multi-component mixtures from atomistic simulations.

Adaptive loss-balancing schemes for physics-informed neural networks rest on a premise that every residual should be driven to zero. For level-set advection with an eikonal regulariser that premise fails: the eikonal term penalises deviation of $\lVert\nabla\phi\rVert$ from unity, a property transport preserves only under rigid motion; where the exact solution departs from a signed-distance function the eikonal residual of the correct answer is nonzero, and driving it to zero moves the network away from that answer. We show that standard gradient-norm balancing fails in exactly this way, its weight remaining near its initial value throughout training on benchmarks where the property is violated, and we introduce SDF-Aware Weighting (SAW), which combines a residual-quantile gate with a gradient-norm ratio so that points exhibiting legitimate departure are excluded before the surviving term is scaled. Across four three-dimensional benchmarks SAW selects an eikonal weight within an order of magnitude of the value located by an eighteen-run manual sweep, spanning four decades from $10^{-1}$ to $10^{-5}$ with a single fixed configuration. On the slotted sphere, where the initial field is non-differentiable at reentrant edges, SAW attains a lower error than any weight in that sweep. Two smooth rigid benchmarks serve as controls: SAW is worse there, as expected when its premise does not hold. An ablation with the gate disabled shows the slot is nearly entirely filled while the relative $L_2$ error reads $1.06\%$, indistinguishable from a field that never represented the slot. We give a feature-restricted measure that separates the two cases.

Although the Alcubierre Warp Drive is theoretically capable of providing faster-than-light travel, it may be difficult to use for this purpose. But is it useful for slower-than-light travel? We begin by observing the that the warp bubble will act to protect the ship from dust particles and other space debris (a potentially serious hazard even at 10% the speed of light). We then explore several modifications of the Alcubierre Warp Drive, e.g. a "deflector shield", with the perspective of keeping a ship safe from collisions with particles, projectiles, rogue planets, and other dangers of space travel.

This paper presents a family of high-resolution weighted essentially non-oscillatory compact least-squares schemes with implicit time integration for the compressible Navier-Stokes equations on curvilinear grids. Compared with the original compact least-squares schemes, the proposed method introduces two main improvements. First, instead of enforcing the accuracy constraints over the entire computational domain, including discontinuous regions, it constructs the reconstruction matrix only along smooth reconstruction lines using an accuracy-preserving weighting strategy. This treatment effectively suppresses the persistent high-frequency oscillations observed in the original compact least-squares schemes and yields sharper profiles near discontinuities. Second, the method simplifies the shock-capturing procedure and improves efficiency by using a single set of polynomials, whereas the original compact least-squares schemes require both unlimited and limited polynomials. Combined with spectral optimization, the proposed method exhibits more favorable spectral properties than conventional weighted essentially non-oscillatory schemes. The smoothness indicators and penalty matrices are constructed through an efficient iterative procedure with modest additional cost. Numerical results for inviscid and viscous one-, two-, and three-dimensional flows demonstrate that the proposed method provides robust shock-capturing capability while maintaining high resolution in smooth regions and across contact discontinuities.

In aircraft engines, functionally graded abradable coatings are used to control blade-tip clearance, but their durability is governed by effects that are often treated separately in existing models, including temperature-dependent viscoelastic softening, progressive damage, deposition-induced microstructural modulation, and geometric tolerances. This study integrates these effects within a unified multiphysics–probabilistic framework. The solved domain is a local through-thickness coating column driven by prescribed strain and temperature histories that include thermal eigenstrain and coating–substrate expansion mismatch. The results show that periodic property modulation increases end-of-cycle damage relative to the monotonic gradient and concentrates the maximum stress within a localized modulation crest. Propagating symmetric geometric tolerances through the coupled solver produces a strongly right-skewed damage distribution and a higher classification-level exceedance probability for the modulated gradient than for the monotonic gradient. This exceedance behavior cannot be obtained from a deterministic analysis performed at nominal geometry. The implementation is verified through spatial and temporal refinement, closed-form unit tests, nested Monte Carlo convergence with Wilson confidence intervals, and a hold-out-validated polynomial response surface. Numerical verification is clearly distinguished from experimental validation, and the simulated case is presented as a generic benchmark rather than an identified material model. The framework links deposition parameters and tolerance bands to damage-exceedance risk, supporting reliability-oriented screening of graded and multilayer coating systems.

Simulating the coupled, nonequilibrium dynamics of electrons and nuclei is a central challenge in chemistry, physics, and materials science, governing phenomena from photocatalysis to quantum information. The primary bottleneck has been the lack of a general, accurate, and efficient method for modeling the complete excited-state landscape: the potential energy surfaces, forces, and non-adiabatic couplings for multiple electronic states. While machine learning has revolutionized ground-state simulations and shown promise for excited states in molecules, a unified framework that solves the complete multi-state problem for general condensed matter systems has remained elusive. Here we introduce on-the-fly N${^2}$AMD (Neural network NAMD), a machine learning framework that makes on-the-fly NAMD in solids a reality. By employing an equivariant neural network to predict the system Hamiltonian, the framework delivers excited-state energies, forces, and non-adiabatic coupling vectors at a fraction of the cost of ab initio calculations. Crucially, it allows simulations with hybrid functional accuracy, a level of approach previously inaccessible for NAMD. We showcase its capabilities with three topical examples: correcting order-of-magnitude errors in carrier dynamics predicted by conventional procedure in a MoS$_2$/WS$_2$ heterostructure, simulating previously inaccessible photoinduced ferroelectric switching, and capturing real-time polaron formation in TiO$_2$ at the hybrid-functional level. On-the-fly N${^2}$AMD moves beyond the limitations of equilibrium theory, establishing a new paradigm for the predictive, first-principles design of materials operating far from equilibrium.

What is "structure," and how can we exploit it in combinatorial optimization? Reasoning tasks such as system design, state estimation, and prediction can be cast as combinatorial optimization problems (COPs), traditionally attacked by dedicated algorithms that fail outside their narrow subclass, while communities working on equivalent problems reinvent each other's techniques. The weighted constraint satisfaction problem (WCSP) is a general framework that subsumes COPs from many communities and represents complex physical and non-physical systems, from spin glasses to social networks. Is such generality a curse? Our answer is to exploit "structure": a general-purpose WCSP algorithm should automatically imitate the specialized algorithm for whatever subclass its input belongs to. A WCSP has macro (graphical) structure, which variables interact, and micro (numerical) structure, how they interact. Separate schools of thought exploit one or the other, never both. The constraint composite graph (CCG), introduced in 2008, unifies them: it is an undirected graph over a WCSP instance's variables plus auxiliary ones, on which minimum weighted vertex cover (MWVC) solves the original instance. It is efficiently constructible, but largely unexploited. This dissertation answers three questions affirmatively. The CCG has theoretical advantages beyond identifying tractable classes (new properties proved here); it is practically useful (efficient implementation and experiments); and it extends promisingly to non-Boolean variables (new encodings). We use the CCG to (a) kernelize a WCSP instance, fixing optimal values of some variables by maxflow before search begins, (b) improve min-sum message passing, (c) exploit integer linear programming solvers, and (d) solve COPs on quantum annealers. Solving general COPs faster classically also yields better baselines for the debated quantum advantage.

We present PASCHEN-1D (Plasma Advanced Solver with Coupled High-fidelity Emission and external Network), a one-dimensional time-dependent fluid plasma solver developed for self-consistent simulation of gas discharges and plasma breakdown with coupled electrode surface emission and flexible external circuit networks. The code solves drift-diffusion continuity equations for electrons and ions together with Poisson's equation. It is dynamically coupled to lumped RLC circuits, which self-consistently treat plasma transport, plasma-surface interaction, dielectric effects, and circuit response within a single framework. The electrode emission module includes ion-induced secondary electron emission, Fowler-Nordheim and Murphy-Good field emission, Richardson-Dushman thermionic emission, and photoemission based on a general, exact quantum mechanical emission theory. A finite-volume formulation with Kurganov-Tadmor fluxes, explicit diffusion, and fourth-order Runge-Kutta time integration is employed to ensure stable transient (sometimes ultrafast) evolution across breakdown and glow regimes. The solver is validated against multiple benchmark cases, including nanosecond pulsed dielectric-barrier discharges, DC breakdown and glow transitions, and Paschen curve construction for argon and nitrogen, with results consistent with published studies. With high-fidelity emission physics and a flexible circuit-coupling framework, PASCHEN-1D provides a versatile and efficient tool for modeling breakdown and transient discharge phenomena.

Sampling high-dimensional probability distributions is a central task in scientific computing, with applications ranging from Bayesian inference to statistical physics and molecular simulation. Despite decades of methodological developments, two major challenges remain: scaling to high dimensions and efficiently exploring multimodal distributions characterized by metastable states. Classical approaches such as Markov chain Monte Carlo, tempering methods, or enhanced sampling based on collective variables have achieved major successes, but they also face intrinsic limitations. This tutorial review explores a new paradigm that has recently emerged at the interface of machine learning and computational statistical physics: the use of generative models as tools for sampling. In this context, models such as normalizing flows and diffusion models are not used in their traditional data-driven setting, but rather as flexible probabilistic models that can assist the sampling of distributions known only up to a normalization constant. This manuscript reviews the early development of this rapidly evolving field and discusses several methodological directions, including exact samplers based on generative models and strategies to train such models in the absence of data. While an exhaustive survey of the literature is not attempted, we present a selection of key ideas and methods, along with a discussion of their strengths and limitations. The review is intended to be an accessible tutorial for both physics and machine learning audiences, and it aims to provide a starting point for researchers interested in exploring this exciting area of research.

Vacancy formation energies govern diffusion, irradiation damage, phase stability, and dynamic failure in high-entropy alloys (HEAs), yet their strong dependence on local chemical environments and mechanical deformation makes atomistic calculations prohibitively expensive for large-scale studies. Here, we develop an atomistically informed permutation invariant machine learning framework for predicting strain dependent vacancy formation energies in FCC HEAs from local atomic environments. The model employs a vacancy-centered representation constructed from objective geometric descriptors together with invariants of the local deformation gradient, enabling the coupled effects of chemical disorder and finite deformation to be learned within a unified framework. Atomistic simulations reveal that volumetric deformation is the dominant factor controlling the average variation in vacancy formation energy, whereas shear deformation has a comparatively minor influence. At the same time, substantial site to site variability persists under identical macroscopic loading, demonstrating that local chemical environments govern the statistical distribution of vacancy energetics beyond species-averaged trends. The proposed framework accurately predicts vacancy formation energies across diverse deformation states while providing orders-of-magnitude faster evaluation than direct atomistic simulations. These results establish an efficient route for incorporating stress-dependent defect energetics into multiscale models of diffusion, irradiation damage, and dynamic failure in chemically complex alloys.