Papers

We expand the correlation-consistent effective core potentials (ccECPs) library by developing semi-local pseudopotentials and matching basis sets by heavy-elements from $5d$ (Hf, Os, Hg) and $6p$ (Tl, Po, At, Rn) blocks. In order to accurately capture scalar relativistic effects, spin-orbit coupling, and electron-electron correlation, we implement a tiered core-valence partitioning strategy across three distinct resolutions. This includes a small 60-core (Hf, Os, Hg) that explicitly correlates subvalence shells, a large 78-core definition for the main-group elements that rigorously accounts for core polarization and relaxation effects in sparse valence environments, and an intermediate 68-core partition for Hg and Tl. This 68-core architecture represents a unique development in the ccECP library, optimizing the balance between accuracy and computational efficiency in a manner unexplored for ligther elements. Optimized against relativistic all-electron CCSD(T) references, the ccECPs deliver outstanding atomic precision, achieving a global average atomic low-lying states deviation of just 0.045 eV. This accuracy translates directly to robust molecular transferability, systematically restricting dissociation energy discrepancies to under 0.03 eV, equilibrium bond lengths to within 0.005 Å. By enforcing a regularized, finite potential at the origin for enhanced numerical stability in stochastic quantum Monte Carlo methods, this library removes a critical methodological bottleneck for predictive many-body simulations of heavy-element systems and materials.

A robust finite-volume framework is presented for the simulation of compressible multiphase flows with surface tension across a wide range of Mach numbers. The method is based on a two-pressure, six-equation diffuse interface model incorporating viscous, gravitational, and capillary effects through the continuum surface force formulation. To consistently account for capillary-induced pressure jumps, an HLLC Riemann solver is developed using generalized Riemann invariant analysis, and the instantaneous pressure relaxation procedure is modified to preserve the Laplace pressure jump during phase equilibration. To overcome the excessive numerical diffusion of conventional approximate Riemann solvers in the low-Mach regime, a robust low-Mach correction is proposed by extending our previous formulation with a modified scaling strategy that remains stable in regions of strong pressure variation. The resulting method retains accuracy from nearly incompressible flows to compressible regimes while preserving the robustness of the six-equation formulation. The numerical framework is validated using a series of benchmark problems involving surface tension, viscosity, gravity, and compressibility. The results demonstrate accurate prediction of interface dynamics, capillary pressure, and low-Mach flow features, while significantly reducing numerical dissipation without compromising stability. The proposed methodology provides an efficient and reliable approach for the simulation of complex multiphase flows spanning a broad range of flow regimes.

This paper presents a unified invariant-domain-preserving (IDP) framework for hybrid discretizations of hyperbolic conservation laws, including active flux and PAMPA methods, in which cell averages are updated conservatively while cell-boundary point values evolve under a possibly non-conservative operator. The main challenge is to preserve admissibility for these two coupled sets of states under a single local stability condition, without adding evolved degrees of freedom or relying on post-update repairs. For the point-value update, we introduce an admissibility transform based on a barrier–Legendre map for convex admissible interiors described by concave constraints, and prove that the inverse map is globally defined and Lipschitz continuous on the relevant sets. For the conservative cell-average update, we establish a structural obstruction theorem: the single-state continuous physical flux built from admissible boundary traces alone cannot provide a conservative IDP guarantee, for any prescribed CFL number, when internal reconstruction values are uncontrolled. This identifies the missing local control that must be supplied by an additional IDP flux mechanism. To realize this mechanism explicitly, we combine cell average decompositions (CAD), geometric quasilinearization, and local a priori scaling to construct admissible, generally discontinuous trace states before flux evaluation. Under an explicit trace-based CFL condition, with constants determined by local CAD weights and trace-state wave-speed bounds, the coupled hybrid update preserves the prescribed invariant domain. Concrete third-order schemes are developed on triangular, Cartesian, convex quadrilateral, and general convex polygonal meshes. Numerical results demonstrate the designed order of accuracy for smooth solutions and the strict preservation of physical admissibility.

Identifying ground-state and low-energy atomic configurations is a central problem in alloy design. The cluster-expansion (CE) method represents configurational energetics on a fixed lattice as an effective Hamiltonian; for binary alloys, higher-order CE models become polynomial Ising Hamiltonians. Using Au-Cu as a model binary alloy, we formulate cubic and quartic cluster-expansion Hamiltonians as penalty-augmented polynomial unconstrained binary optimization (PUBO) problems under fixed-composition constraints. We optimize these PUBO problems using SQBM+, a simulated-bifurcation-based solver that can treat higher-order polynomial binary objectives directly. This direct PUBO treatment avoids the need to construct an explicit quadratic reformulation with auxiliary variables. Composition constraints are imposed through quadratic penalty terms, whose weights are estimated from derivative coefficients of the continuous relaxation of the CE objective. Benchmark calculations for systems up to 2048 atoms show that SQBM+ robustly obtains low-energy feasible configurations for cubic CE models and remains effective for many quartic instances. Formation-energy convex hulls constructed from the optimized configurations recover the CuAu and Cu3Au ordering trends and reveal finite-size effects at off-stoichiometric compositions. These results demonstrate simulated bifurcation as a practical route to constrained higher-order CE optimization for alloy configuration search.

Agitated filter dryers (AFDs) are a crucial unit operation in small molecule manufacturing that enables simultaneous filtration, washing, and drying of active pharmaceutical ingredients. One of the key challenges in AFDs is associated with undesired agglomeration, where the presence of hard agglomerates results in off-spec products, equipment damage, and additional downstream processing. This article presents a novel mechanistic model that describes the formation of soft and hard agglomerates during agitated filter drying. By integrating population balance and multiphysics modeling, the model can accurately predict the evolution of the product temperature, moisture content, and particle size distribution, and hence quantify the extend and impact of undesired agglomeration across various operating conditions. Our proposed model-based framework enables the rational design and operation of AFDs for improving the product quality and process reliability.

A narrowband emitter aligning with the bandgap of the underlying solar cell is essential for improving the spectral efficiency and thermal stability of thermophotovoltaic (TPV) systems. Emitters based on oxide materials present a promising solution to the optical and mechanical performance degradation of traditional emitters, which employ metal-dielectric structures that experience metal oxidation and structural deterioration at high temperatures. Here, we presented a surrogate-assisted inverse-design framework for a narrowband 1D grating emitter comprising ITO and Al2O3 layers on a sapphire substrate. We performed Bayesian optimization over the trained ExtraTrees surrogates on a penalty-augmented objective containing a peak-emission constraint (E_peak > 0.90) while minimizing the full width at half maximum (FWHM) and maximizing the fraction of emission concentrated within the selected peak-centered spectral band to acquire a high, narrowband peak emission. The surrogate model was trained using a dataset of emission spectra obtained from the finite-difference time-domain (FDTD) by systematically varying the layers' thicknesses and the structure's period as the input features, and the narrowband emission's figure of merit (FOM) ((E_peak), wavelength of peak emission (_peak), FWHM, in-band fraction (f_in) for determining the amount of emission outside the peak band, and concentration of peak emission near the peak) was used as the prediction target. The resulting set of predictions for the optimized structure was further validated using the FDTD method.

Directly self-starting implicit methods are attractive for transient analysis because they avoid auxiliary starting procedures while retaining the original first- or second-order governing equations. However, most existing formulations usually fix the last sub-step at the end of each time interval, which restricts the attainable order. This study develops a generalized $s$-sub-step implicit framework by releasing this constraint and treating all sub-step locations as design variables. The resulting methods admit a unified Runge–Kutta representation for both first- and second-order transient systems and preserve identical effective matrices over all sub-steps. Accuracy conditions are derived by simultaneously matching the numerical amplification factor and load operator, thereby accounting for both homogeneous and forced responses. For $s=1,~\cdots,~6$, two complementary families are obtained: $s$th-order members with user-controllable high-frequency numerical dissipation and adjustable sub-step locations, and $(s+1)$th-order members obtained by selecting the sub-step locations, with fixed dissipation. The latter reach up to seventh-order accuracy without increasing the number of sub-steps, although some high-order members are $A(\alpha)$-stable with stability angles extremely close to $90^\circ$. Analytical amplitude and phase errors further reveal parity-dependent superconvergence in undamped systems, and appropriate parameter selections can substantially increase either phase or amplitude accuracy beyond the formal order. Numerical benchmarks confirm the predicted convergence orders and the controllable suppression of spurious high-frequency responses.

A direct, analytical, non-iterative ghost-cell reconstruction framework is developed for Cartesian-grid embedded-boundary methods. Analytical expressions impose Dirichlet and Neumann boundary conditions directly at the embedded boundary, eliminating intermediate image-point reconstruction, matrix inversion, and precomputation or storage of geometry-dependent reconstruction weights. For the Cartesian stencil considered, dependencies among neighboring ghost cells form a directed acyclic graph. A topological ordering partitions ghost cells into dependency levels, permitting level-by-level reconstruction without iterative updates. The formulation is combined with hybrid ghost cells (HGC), whose centers may lie on either side of the embedded boundary. This placement satisfies the linear reconstruction-stability criterion previously derived for scalar advection while retaining a compact nearest-neighbour Cartesian stencil. An extended-stencil classical ghost-cell formulation (CGC_ES) serves as a stability-preserving reference, distinguishing the effects of reconstruction stability and stencil compactness. The same analytical relations provide solution values and spatial gradients directly on the embedded boundary, enabling evaluation of pressure forces, wall stresses, drag, and lift without separate surface reconstruction or filtering. Simulations of flow past a circular cylinder and an airfoil show that the linear criterion remains a useful indicator of reconstruction stability for the nonlinear incompressible Navier-Stokes cases considered. Classical ghost-cell reconstruction develops spurious oscillations when the criterion is violated, whereas CGC_ES and HGC remain stable. HGC additionally satisfies the stability requirement with a compact nearest-neighbour stencil.

Modal methods are particularly well suited to layered photonic structures because electromagnetic propagation within each layer is described analytically. Here, we introduce a simple reformulation of the scattering-matrix (S-matrix) approach in which the S-matrix is defined as an intrinsic property of an individual layer, independently of its neighboring interfaces. This separation between layer propagation and interface coupling preserves the numerical stability of the conventional S-matrix formalism while providing a more modular description of multilayer structures. It also enables direct computation of the scattering coefficients between external plane waves and Bloch modes, as well as between Bloch modes themselves, providing generalized Fresnel coefficients for periodic interfaces. In addition, the reformulation allows the optical response to be evaluated simultaneously for many layer thicknesses or incident wavevectors, with only a modest computational overhead compared with a single calculation. Implemented in the RETICOLO freeware, these capabilities provide a practical framework for modal analysis, parameter sweeps, and the design of layered photonic structures.

The three-dimensional electron gas is a fundamental model in condensed matter physics and quantum chemistry, and the exchange-correlation energy derived from it serves as the starting point of ab initio computations of materials. However, the wave-vector-dependent response, and hence the static local field factor $G(q)$, has remained without accurate ground-state benchmark in the strongly coupled regime for three decades. Using diffusion Monte Carlo, we calculate $G(q)$ and the static density-density response function across the liquid phase and find a pronounced structure in $G(q)$ at intermediate wave vectors, already visible at metallic densities and growing with increasing interaction strength. We identify it as a fingerprint of short-range order by showing that it is required to reproduce the static structure factor. Our parametrization, valid in the entire liquid phase, predicts a low-energy resonant mode inside the particle-hole continuum.

Context. Scattering polarization in the Sr I 4607 Å line is a key diagnostic for small-scale, unresolved magnetic fields in the quiet solar photosphere, inaccessible to Zeeman-based techniques. The complete frequency redistribution (CRD) limit is commonly used to model this line, but its validity and impact on magnetic sensitivity have not been systematically investigated. This is timely given new facilities observing this line with unprecedented accuracy, and future synoptic programs targeting long-term quiet-Sun magnetism. Aims. We assess the CRD limit against the general partial frequency redistribution (PRD) description of scattering for the Sr I 4607 Å polarization, focusing on its magnetic sensitivity via the Hanle effect. Methods. We solved the non-LTE radiative transfer problem for polarized radiation in a semi-empirical 1D plane-parallel static atmosphere, with CRD and PRD scattering, across magnetic field configurations relevant for Hanle diagnostics. Results. For small-scale unresolved fields, where the field acts only via Hanle depolarization of Q/I, CRD and PRD emergent profiles and Hanle sensitivity are essentially identical. For deterministic fields, CRD remains accurate when polarization signals exceed about 0.3%; for weaker signals, PRD effects can appreciably affect the emergent profiles. Conclusions. CRD is adequate for modeling Sr I 4607 Å scattering polarization in the vast majority of observationally relevant cases; PRD effects matter only for weak polarization signals.

This study investigates the nonlinear kinematics of a premixed slit flame using a two-dimensional $G$-equation level-set framework. Results show that combined forcing induces nonlinear saturation in the FDF, characterized by early gain flattening and premature phase drops, which intensify with the transverse forcing amplitude. Kinematic analysis reveals that this geometric nonlinearity manifests as a reduction in the time-averaged flame height, defined as the mean shift. In the quasi-steady limit, this mean shift is analytically quantified via a multivariate asymptotic expansion, where fourth-order terms successfully capture the saturation mechanism at elevated amplitudes. By introducing a scaling parameter to account for transverse dominance, the frequency-dependent decay of the mean shift in the compact limit collapses onto a single master curve, enabling the derivation of a unified theoretical model that integrates this asymptotic response with a second-order low-pass filter. Furthermore, because the mean shift reduces the physical extent of the flame, it alters the wrinkle propagation time. Correcting the Strouhal number using the measured mean shift collapses the dispersed nonlinear FDF curves onto the linear theory prediction. The analysis is further extended to disturbances convected at a finite speed, for which the linear transfer function is derived analytically and the correction with the measured mean shift continues to collapse the nonlinear FDF. These findings establish that the nonlinear FDF behavior under multidimensional forcing is fundamentally governed by the kinematic mean shift, providing a theoretical baseline for decoupling geometric nonlinearities from other thermo-diffusive or hydrodynamic instabilities in turbulent flames.

The Carleman lattice Boltzmann method replaces the quadratic collision of the lattice Boltzmann equation by a linear map on a lifted state. Existing block encodings of this map succeed with probability $10^{-5}$ to $10^{-2}$ per step, so that the probability of a multi-step evolution decays exponentially with the number of steps. This study develops multi-relaxation-time Carleman lattice Boltzmann, a representation split in which the Carleman-lifted collision is performed in a moment basis on the amplitude-encoded state, while streaming remains a permutation in population space. The level-2 hierarchy closes under collision but not under streaming, which generates the work of pressure, advection and stress on the momentum flux. The state therefore remains the full pair array, and the moment representation is an orthogonal change of basis of the population one, with identical dynamics and truncation error. The operator to be block encoded changes with the basis, and with it the success probability of a quantum step. In the weighted Hermite basis with the amplitude encoding $f_i/\sqrt{w_i}$, the relaxation is a diagonal operator of norm one, and the level-2 collision stage becomes this diagonal plus a sparse coupling from products of conserved moments. A structured block encoding of this form raises the single-step success probability on the D2Q9 lattice from $10^{-2}$ for the best population encoding to 0.98. With the pair sector in relative coordinates and scaled, one step requires $O(\log^2 N)$ Toffoli gates and is verified by statevector simulation. The encoding and readout overhead of a $T$-step run is then polynomial in $T$ and $N$ rather than accumulating a factor $p^T$. For resolved flows, the remaining exponential loss in coherent multi-step Carleman lattice Boltzmann evolution is set by the non-equilibrium norm relaxed by dissipation rather than by the block encoding.

Dissimilar Al/Cu joints are essential to battery-pack assemblies; however, their mechanical strength is limited by brittle Al-Cu intermetallic compounds (IMCs) such as Al2Cu and Al4Cu9. Interlayer strategies to suppress these phases remain largely empirical, lacking a predictive framework linking interlayer chemistry to the phases that form and to their intrinsic mechanical character. Here an Ag interlayer is introduced and combines computational thermodynamics, first-principles calculations, microstructural characterization, and mechanical testing into a single self-consistent description of the joint. CALculation of PHAse Diagrams (CALPHAD) equilibrium and Scheil simulations predict the solidification path of the Al-rich Al-Ag fusion zone and explain why Cu incorporation is limited when Ag is present; energy-dispersive X-ray spectroscopy (EDS) and electron backscatter diffraction (EBSD) confirm an FCC Al-Ag solid solution as the dominant constituent. First-principles phonon calculations within the quasiharmonic approximation yield finite-temperature entropy and Gibbs energy, benchmarked against CALPHAD, while elastic constants assess ductility via the Pugh criterion (i.e., the bulk/shear (B/G) modulus ratio). All Al-Ag phases, including the observed solid solution, exceed the Pugh threshold of 1.75, whereas the targeted Al-Cu IMCs do not, giving a mechanistic basis for the interlayer's effectiveness. This microstructural change translates into improved performance: nominal strength rises from 47.9 to 67.4 MPa. Nanoindentation gives a fusion-zone reduced modulus of 82.8 GPa (Young's modulus 82.0 GPa), versus a calculated 0 K Voigt-Reuss-Hill value of 71.4 GPa. The present work establishes a transferable CALPHAD, first-principles, and experiment workflow for rational interlayer selection in dissimilar-metal joining.

N-body simulations are an essential tool for modelling the evolution of gravitating systems. Differentiable N-body simulations allow solving complicated reconstruction problems by connecting observations of evolved systems to simple informative priors on their initial conditions. While so far most reconstruction efforts focus on mildly non-linear large scale dynamics, the internal dynamics of galaxies and clusters should exhibit rich information about their formation. The exploration of reconstruction possibilities in this highly non-linear regime naturally requires highly performant differentiable N-body simulations with small scale forces. Here, we present JZ-FMM, a GPU-native differentiable implementation of the fast multipole method for N-body simulations. We show that the GPU oriented design of tree structure, dual tree traversal and multipole translation operators improves performance over established CPU-native and hybrid N-body codes by notably more than an order of magnitude. Further, we show that gradients can be evaluated efficiently and accurately so that JZ-FMM can be used as an essential building block in future field-level reconstruction efforts. We apply the new code to the reconstruction problem of a tidally stripped satellite to show that deeply non-linear problems can indeed be solved efficiently with simulation gradients, although appropriate care must be taken to navigate the complicated optimization landscape.

In modern density functional theory, effective Hamiltonians are constructed to reproduce densities and forces of electrons at equilibrium. However, their nonlocal potentials leave systematic errors in spectral properties and real-time dynamics. Access to local optimized or inverted effective potentials would remove this limitation, but the numerical fragility in finite orbital bases has long prevented their wide adoption. Here, we introduce multipole splats, a class of trial potentials that carry the correct asymptotic decay required to support the unoccupied spectrum. By connecting the computation of effective potentials to variational and supervised variants of Hamiltonian learning, we recast both problems as stable nonlinear optimization formulated directly in standard orbital basis sets and applicable to any hybrid functional approximation. The resulting solver allows us to resolve spatial profiles of exchange-correlation potential errors during molecular dissociation and accurately reconstruct key excited states without empirical asymptotic corrections. We also show the deviation from the ionization potential theorem for different exchange-correlation approximations on a dataset of molecular systems. Multipole splats extract insights from established approximations and provide capacity for robust dataset generation for downstream processing and learning.

The energetic effect of a force error depends on atomic motion. We establish a directional residual-work coefficient combining directional curvature mismatch with the spatial distribution of the residual response. For conservative potentials force-matched at an anchor, it determines the leading signed work at the first crossing of a small force-error budget. At a 474-atom lithium-electrolyte interface, predictions fixed before future reference evaluations differ from measurements by less than 4.6% of predicted work across 24 prescribed endpoints. Changing only the initial velocity direction at fixed structure and initial total kinetic energy reverses the force-work ranking. At the same admitted time of 0.25 fs, one direction gives an 11.6% larger maximum force residual but 36.1% less work. The reversal recurs at a second structure. The framework connects force tolerances to reference-energy transfer, providing a physical basis for potential assessment and adaptive reference allocation.

In the Path Integral Monte Carlo (PIMC) simulation of fermions in a harmonic trap, with and without pairwise harmonic interactions, the partition functions for any discrete number of imaginary time slices (or beads) and for any choice of the short-time propagator can be analytically obtained from the contracted determinant form of the propagator. This work shows that the resulting recursion relation can be reformulated in the $\lambda$-ring language, yielding a closed-form finite-bead partition function in two dimensions in terms of permutation statistics. This closed-form partition function can be evaluated by a special algorithm in $O(n^2)$ time, providing an exact and numerically stable scheme for reproducing the energies of the original (undoable) fermion PIMC simulation for $n=10^4$ or more fermions. This result provides a concrete framework in which the numerical instability of the sign problem is completely bypassed, and serves as a counterexample to the prevailing view that all truly fermionic PIMC sign problems are NP-hard.

Molecular diffusion in fluctuating amorphous and macromolecular media governs key transport processes across soft-matter physics, energy storage, and biological membranes. Extracting localized trapping states from single-particle tracking trajectories remains a fundamental challenge; because thermal structural breathing continuously reconfigures pore boundaries, conventional geometric algorithms suffer from severe systematic biases, erroneously merging distinct localized states during cyclic molecular returns. Here, we address this deadlock by shifting the paradigm from local geometric recurrence to a structure-informed Bayesian regularization. Leveraging discrete Morse theory, we extract the time-invariant topological skeleton of the fluctuating host matrix to construct robust, gas-specific physical priors that account for individual molecular dimensions. Trajectory steps are sequentially partitioned via a two-stage probabilistic refinement that dynamically adapts to the transport landscape. Benchmarked against a rigorous environment where synthetic particles explore the actual interconnected matrix graph, our approach eliminates systemic biases, restricting macroscopic trapping parameter deviations to just a few percent under optimal linear $O(N)$ computational scaling. Applied to hydrogen and methane transport within a type-I kerogen matrix, serving as a prototype for highly tortuous, flexible macromolecular networks, the method successfully decodes the hidden microscopic mechanisms of confined diffusion. To ensure immediate broad impact, the documented open-source code and data are made publicly available, offering an accessible strategy readily adaptable to a broad spectrum of tracking phenomena, from ion transport in battery polymers to protein trafficking within cellular environments.

Direction-dependent instrumental polarisation is a major systematic limitation in high-fidelity single-dish radio polarimetry, yet a unified characterisation of beam leakage and the conditioning of polarisation recovery is lacking for the Ghana Radio Astronomy Observatory (GRAO) 32-m telescope. We aim to establish a quantitative, direction-dependent polarimetric beam model at 5.0 and 6.7 GHz. High-resolution GRASP electromagnetic simulations of the nominal telescope configuration are used to derive the full complex Jones response and corresponding Mueller matrices, from which beam shape, beam squint, instrumental Stokes leakage, and intrinsic cross-polarisation ratio (IXR) are evaluated across the primary beam. The Stokes $I$ beams have half-power beam widths of $396.4$ and $295.3$ arcsec at 5.0 and 6.7 GHz, respectively, with main-beam efficiencies of $56.4$ and $55.6$ per cent. At 5.0 GHz, the circular-polarisation beam squint is $9.96$ arcsec ($2.5$ per cent of the HPBW), whereas no statistically significant squint is detected at 6.7 GHz. Leakage from Stokes $I$ into linear polarisation remains below $0.1$ per cent within the half-power beam but increases substantially towards the sidelobes. IXR reaches approximately $80$ dB on axis and decreases with angular offset. These results establish the intrinsic electromagnetic polarimetric response of the GRAO 32-m telescope and provide a quantitative baseline for direction-dependent calibration and subsequent observational validation.

Many dielectric targets encountered in 3-D electromagnetic inverse scattering exhibit piecewise-homogeneous material distributions, which provide useful structural information for inversion. To exploit this property, this paper proposes a structured multi-level-set inversion (SMLSI) method within a physics-driven neural network (PDNN) framework. Instead of directly parameterizing the contrast distribution, SMLSI represents it in terms of object support and material distribution. Multiple level-sets are mapped to soft-support functions and combined through a soft-union operation to construct the object support. The normalized contributions of these soft-support functions, together with trainable complex relative permittivities, determine the material distribution within the support. The level-set residuals are parameterized by a coordinate-based neural network, and the trainable parameters are optimized from scattered field data. Structure-aware regularization terms are introduced to promote binary object support and spatially coherent material regions while reducing material mixing and redundant activation of level-set supports. Adaptive regularization weighting is further employed to reduce the dependence on manual parameter tuning. Numerical experiments demonstrate that SMLSI improves reconstruction accuracy over direct contrast parameterization and enables the recovery of irregular 3-D geometries, closely spaced objects, and multiple material regions. Compared with CSI and $L_{2/3}$-FBE-WCIE, SMLSI achieves more accurate reconstructions on representative synthetic targets. Experimental validation using 3-D Fresnel data further demonstrates its applicability to measured scattering data.

Every data-driven equation discovery method assumes (implicitly and without verification) that the target system obeys a single governing law ($K{=}1$). Here we show that this assumption is the primary bottleneck limiting scientific discovery in multi-mechanism systems, and introduce Physical Law Ecology, a framework that makes $K^*$ (the number of coexisting independent mechanisms) itself the first quantity to be determined from data. The framework automatically mines a pool of topologically distinct candidate equations, constructs a continuous dominance weight field across parameter space, and discovers analytic evolution laws governing mechanism succession—with optional monotonicity constraints encoding irreversible physics. Across four unrelated systems (elastomer mechanics, pool boiling, galactic dynamics, and droplet evaporation), BIC consistently identifies $K^*{=}3$ independent governing topologies. Applied to 163 SPARC galaxies (3,269 spatially resolved measurements), the framework autonomously recovers three gravitational laws whose coexistence provides evidence against the single-universal-acceleration hypothesis of MOND ($p<10^{-34}$). In engineering applications, multi-law weighted prediction reduces error by 67-72% over single-equation baselines while retaining full interpretability. By establishing the determination of $K^*$ as the zeroth step of scientific discovery-prior to and independent of equation search—this work opens a direction orthogonal to existing symbolic regression: not finding better equations, but mapping the ecology of mechanisms that govern complex systems.

We develop a nonperturbative functional renormalization group framework within the LPA' to solve the Wetterich flow equation for the $Z_2$-symmetric real singlet extension of the Standard Model at finite temperature, without a low-order polynomial truncation of the loop corrections to the effective potential, using a physics-informed neural network (PINN) representation. In contrast to conventional truncations based on low-order field expansions, our hybrid tree-level-plus-neural-network ansatz yields a continuous, mesh-free description of the effective potential over the full field and scale range. As a proof of concept, we apply the framework to the finite-temperature effective potential relevant to the electroweak phase transition: one-dimensional field-space slices at two benchmark temperatures, and a two-dimensional reconstruction at T=100 GeV yielding a two-step bounce action. The flow is implemented in a multi-domain setup in scale and field space with derivative matching conditions ensuring smoothness and numerical stability. Gauge and Yukawa sectors are incorporated via independently computed perturbative two-loop running couplings, and anomalous dimensions are included in the flow at the LPA' level. We benchmark the network against resummed perturbation theory and against a grid-based relaxation solver of the same equation. Finally, we introduce a soft consistency constraint that keeps the solution close, in sign and order of magnitude, to perturbation theory across the field-space domain. We find this constraint necessary, rather than merely helpful, for selecting a physically sensible solution of the flow equation. The converged result nevertheless retains a residual dependence on this guidance – through hand-tuned weights and an analytic thermal target – which we identify as the central open problem for mesh-free FRG treatments of this kind.

Ice shelves lose mass primarily by calving, a process controlled by the near-front stress field on timescales that span elastic flexure and viscous creep. We formulate a finite-strain Maxwell model for glacier ice in logarithmic strain space. The Hencky strain of a fixed reference configuration is split additively at the level of rates into elastic and viscous parts; the spring is isotropic Hencky elasticity and the dashpot is a Glen-type power law written on the logarithmic strain rate and its work-conjugate stress. At infinitesimal strain the dashpot coincides with Glen's flow law; the elastic strains in the ice-shelf configurations of this paper remain in that regime. The model is integrated with a midpoint evaluation and a backward-Euler correction of the trial dual, and implemented in a finite-element setting. After a viscoelastic column benchmark, the formulation is applied to an idealised ice tongue, including depth-dependent density and moduli, temperature-dependent fluidity, and cliff geometries with a frontal foot or basal undercutting. The resulting stress fields show how viscoelasticity and front morphology control tension near the terminus.

We present an iteration algorithm for vacuum and Einstein scalar-field equations in double-null gauge, which transform the non-linear PDE into systems of ODE. The numerical realization combines characteristic constraint solves, LGL spectral elements, pole-free spherical operators, Galerkin projection, and independent first-order residual and consistency checks.

Energy-based approaches provide a natural and physically consistent framework for a large class of partial differential equations arising in solid and fluid mechanics, where the governing equations follow from variational principles. In contrast to residual-based physics-informed neural networks (PINNs) and their weak-form variants, which enforce the strong or weak form of the equations through loss minimization, the Deep Energy Method (DEM) directly computes the solution as the minimizer of an energy or incremental potential functional. This eliminates the need for residual weighting, avoids high-order derivatives, and enables the direct enforcement of physical constraints through the variational structure.In this work, we systematically revisit the Deep Energy Method, placing it in the broader context of physics-informed learning and variational modeling. We clarify the relationship between DEM, PINNs, and VPINNs, and identify the class of problems for which energy minimization provides intrinsic advantages in terms of stability, robustness and interpretability. Particular emphasis is placed on incremental variational formulations, which allow DEM to be applied to nonlinear, history-dependent and time-dependent problems, including phase-field fracture and dissipative systems. The variational structure underlying DEM further provides a natural foundation for optimization and inverse problems, where the energy functional acts as a physics-based constraint rather than a residual penalty. Through a series of numerical examples, we demonstrate that DEM offers a principled and effective alternative to residual-based methods for variational problems, highlighting its strengths and limitations relative to existing physics-informed approaches.

Machine learning surrogates based on neural operators have shown broad applicability in solving forward PDE problems. However, eigenvalue problems, in which an eigenparameter and one of several valid eigenmodes must be simultaneously solved, remain difficult because standard operator learning formulations assume a unique input-output map. This work demonstrates that Fourier Neural Operators (FNOs), combined with wavelet-based encodings of PDE inputs, can learn and predict multiple eigenmodes of the elastic wave equation, corresponding to deformation modes of acoustic waves propagating through arbitrary metamaterial geometries. We provide a mechanistic explanation and experimental evidence for why wavelet encodings are well matched to the dual spatial-spectral structure of the FNO, enabling deterministic mode selection on both continuous-valued and binary-valued geometries within a single model, and for why prediction accuracy varies with geometric discontinuities. For metamaterial design, the resulting surrogate accelerates the simulation stage of the design cycle by three orders of magnitude relative to finite element analysis on a consumer-grade CPU, while preserving high fidelity. These results also carry broader implications for designing input encodings in other multi-mode PDE solvers based on spectral neural operators.

We tabulate, on the diagonal, the $n$th power of the primitive approximation for the pair Coulomb density matrix on a sphere, with up to $n=5$, and we discuss its behavior at particles contact as $n$ grows. Our results show that even on a curved surface quantum statistical physics is fundamentally different from its classical counterpart. If in the latter one needs to artificially introduce a regularization of the divergent attractive Coulomb potential preventing particles contact, for example with a hard core, in the former this is not necessary.

Encoding classical data into quantum systems is a foundational step in the execution of nearly all quantum algorithms, and a critical bottleneck in realizing practical quantum advantage. This review provides a comprehensive account of the concepts, algorithms, and practical considerations associated with quantum data encoding. We trace the development from its early conceptual foundations to recent advances, considering commonly used access models, such as quantum state preparation, unitary synthesis, QRAM and block encoding. We survey the circuit size, depth, space-time tradeoffs, as well as non-Clifford resources required for fault-tolerant implementation. We also discuss the roles of different access models in quantum algorithms. Special attention is given to structured data, such as sparse data, Boolean functions and data represented by tensor networks. This review bridges theory and applications, serving both as a pedagogical guide for newcomers and as a reference for active researchers. We also highlight the pivotal role of quantum data encoding in quantum computing and provide insights into future directions that will enable quantum advantage.

Runaway-electron avalanches in tokamak plasmas can be limited by kinetic instabilities driven by the non-Maxwellian runaway distribution. We formulate a reduced model for the quasi-steady state in which the total plasma current and bulk electron temperature are prescribed, while the inductive electric field is determined self-consistently from the partition between Ohmic bulk current and runaway-electron current. Because the wave growth time is short compared with the current-decay time, we consider a marginal-stability regime, in which whistler-wave drive by the runaway electrons balances collisional damping. The resulting states separate into three regimes: a subcritical Ohmic regime without an avalanche, an avalanche regime in which runaway growth relaxes the inductive field to the avalanche threshold, and an instability-regulated regime in which self-excited whistler waves enhance momentum-space diffusion and limit the runaway current. In the instability-regulated regime, the whistler wave spectrum forms a narrow ridge, and low-energy runaway electrons carry most of the runaway current.