2026-09-09 | | Total: 35
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 $λ$-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.
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.
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.
FBPIC (Fourier-Bessel particle-in-cell) is a high-performance simulation code for relativistic plasma and accelerator physics. Its original accelerator backend relies on Numba CUDA, which limits its direct deployment on accelerators using the HIP (Heterogeneous-Compute Interface for Portability) programming environment, such as DCU (Deep Computing Unit) accelerators. In this work, we develop an accelerator backend compatible with HIP that enables FBPIC to run efficiently on DCU platforms while preserving its Python user interface and high level simulation workflow. For the evaluated LWFA (laser-wakefield acceleration) workloads, the proposed backend achieves 1.32-1.54x speedups over the original FBPIC implementation on an NVIDIA V100 GPU and enables efficient execution on the DCU platform. We also summarize the key lessons learned from porting FBPIC to the DCU platform. Multi-DCU experiments achieve a 1.88x strong-scaling speedup on four accelerators and a 2.72x increase in aggregate throughput at approximately 68\% weak-scaling efficiency, with communication analysis identifying inter-node communication and synchronization as the main scalability limitations. Beyond FBPIC, the proposed approach provides a practical reference for porting and optimizing other scientific computing applications developed with Python on heterogeneous accelerator platforms.
A nonlinear EM inversion method regularized by $\ell_1 / \ell_2$-norm of the model gradient is proposed. The $\ell_1 / \ell_2$-norm, defined as the ratio of $\ell_1$-norm to $\ell_2$-norm, exhibits a scale-invariant property that enables it to more accurately characterize sparsity in the solution compared with the conventional $\ell_1$-norm. To address the nonconvexity and nonsmoothness introduced by the quotient structure of $\ell_1 / \ell_2$-norm, the inversion problem is reformulated and solved using the alternating direction method of multipliers, resulting in an efficient optimization algorithm comprising five sub-steps. The Gauss--Newton method is incorporated into the first sub-step to linearize the nonlinear inversion. The proposed method is validated by a series of numerical examples using both synthetic and experimental datasets, and is compared against standard regularization methods such as Tikhonov and total variation. The results demonstrate that the proposed method yields superior reconstructions of permittivity profiles, particularly in preserving sharp edges and suppressing artifacts, without incurring additional computational cost. Furthermore, the method exhibits strong robustness to measurement noise and reduced sensitivity to the choice of the regularization weight.
Reconstructing physical fields from sparse observations is central to system identification, forecasting, and control, yet sparse measurements generally underdetermine the full field. This makes reconstruction an ill-posed inverse problem rather than simple interpolation. Although many deterministic and generative methods have been developed, there is still no clear consensus on when a single point estimate is sufficient and when a distribution of plausible reconstructions is more useful. We conduct a fair comparison of a deterministic U-Net, conditional diffusion, and prior-guided diffusion under matched experimental settings, including 2D Poisson equation, 2D Navier-Stokes flow, and 1D Kuramoto-Sivashinsky dynamics. Through this comparison, we make three observations. First, accuracy is field- and regime-dependent, with no systematic advantage for diffusion under higher complexity or sparser observations. Second, ensemble means improve phase-aligned accuracy, whereas individual samples better preserve variability and can retain high-wavenumber power in selected regimes. Third, conditional diffusion provides more reliable uncertainty estimates at lower cost, while prior-guided diffusion is more robust to mask-distribution shifts but requires substantially higher inference cost and guidance tuning. These results clarify when generative reconstruction is useful and provide guidance for improving uncertainty estimation, fine-scale sample fidelity, robustness, and computational efficiency in sparse field reconstruction.
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 \texttt{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.
This paper proposes a level-set-based physics-driven neural network solver (LSPDNN) for 3-D electromagnetic inverse scattering. To mitigate boundary blurring and reconstruction artifacts in voxel-wise contrast reconstruction, the proposed solver exploits the piecewise homogeneity of practical scatterers by representing unknown targets with multiple coordinate-dependent neural level-set components. Specifically, a soft-union multi-material model is proposed to separately describe the object support and material distribution. The global support is formed by the union of multiple level-set components, while the local contrast is determined by normalized component weights and learnable complex permittivity candidates. In addition, a model-consistent total variation (TV) regularization is imposed on the material-region indicators, rather than directly on the reconstructed contrast, to suppress fragmented material assignments without excessively smoothing material interfaces. An adaptive loss balancing strategy is further introduced to reduce the dependence on manually selected regularization weights. For each measurement instance, the neural level-set parameters and material candidates are optimized by minimizing a physics-consistent objective function. Numerical and experimental results demonstrate that LSPDNN can reconstruct scatterers with clear boundaries, more uniform material regions, and substantially reduced background artifacts. The results highlight the advantage of the neural level-set parameterization in challenging 3-D inverse scattering cases involving irregular shapes, closely spaced objects, multiple materials, and measurement noise.
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 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.
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.
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.
The capacitance of thin-film diodes based on undoped semiconductors are dominated by injected charge carriers rather than doping-induced carriers. However, an analytical framework of the capacitance in these devices has remained elusive. Here, we derive an analytical description of the capacitance, fully accounting for injected charge carrier and electrode charge effects. Based on these findings, a method to extract the built-in voltage in these devices is presented. The theoretical framework is substantiated numerically by drift-diffusion simulations and experimentally on organic solar cells.
Grid-based fluid dynamics solvers routinely struggle with exhaustive meshing demands and ill-posed inverse problems. A practical mesh-free alternative is given by Physics-informed neural networks (PINNs) but, applying them to highly elastic Oldroyd-B fluids results in many training failures. The High Weissenberg Number Problem (HWNP), driven by the exponential stress growth near stagnation points is a major issue, which causes standard PINN optimizers to diverge. To prevent the network from crashing, we apply a Cholesky decomposition to the conformation tensor. This mathematical constraint stabilizes the gradients by guaranteeing a positive-definite stress field. Beyond mathematical stability, the inherent spectral bias of deep learning models can hinder the network from accurately capturing the highly elastic wake structures. Therefore, we used sparse data assimilation to force the model toward the actual physical solution. By anchoring the physics loss with targeted CFD data points and accelerating training via transfer learning, we successfully pushed the network past non-physical local minima. We validated this Cholesky-PINN approach on flow past cylindrical geometries for Reynolds numbers (Re) between 5-25, for the single-cylinder setup. In addition, the relaxation time (λ) is increased from 0.1 to 0.5 to test the stability of the network. Finally, we scale the framework to a complex 3-cylinder array which proved our constructive solid geometry approach completely bypasses the tedious re-meshing steps of traditional CFD. The combined framework accurately captured sharp viscoelastic wakes, providing a stable computational tool for complex rheological modeling.
Many systems in soft and living matter bind through a limited number of bonds per particle: proteins associate via discrete surface patches, nucleic acids form one-to-one contacts, and the phase behaviour of multivalent biomolecules is governed by the number of binding sites they carry. In simulations, valency limits are typically enforced with patchy particles, whose anisotropic potentials require integration of rotational degrees of freedom and combine hard cores with narrow patches, which forces small timesteps and commits the model to a fixed binding-site geometry that is often unknown, flexible or mobile. We introduce the valency-bounding correction (VBC), a many-body modification of generic short-range pairwise potentials that smoothly suppresses attraction once the neighbour count of either interacting particle exceeds a prescribed valency. The correction carries no angular degrees of freedom, applies on top of soft repulsive cores and evaluates in two passes over the neighbour list at the cost of a standard pairwise potential. The VBC drives the coordination number to the prescribed valency with low error, while its cluster statistics depart from Wertheim and Flory-Stockmayer predictions through unrestricted ring formation. A tuned variant exchanges bonded partners through ordinary molecular dynamics, reducing bond lifetimes at high saturation by an order of magnitude. On GPUs the cost of the correction is nearly independent of valency, reaching an almost tenfold advantage over a patchy-particle reference. We illustrate large-scale applications by reproducing the reentrant aggregation of repeat-expanded RNA, and show that the VBC also remedies the Fisher-Ruelle thermodynamic instability of soft-core potentials with attraction. The VBC is available as an open-source GPU plugin for HOOMD-blue.
Physics-informed neural networks (PINNs) struggle on PDEs whose governing physics varies across the domain. We trace this to a structural property of standard coordinate networks: their neural tangent kernel (NTK) is translation-variant and lets training points of large coordinate magnitude disproportionately influence predictions elsewhere, producing long-range coupling and gradient conflict during training. We show analytically and empirically that mixture-of-experts (MoE) architectures with centered, compact-support routers yield a uniformly banded NTK whose kernel-regression weights decay exponentially with distance, localizing the learning. Building on this, we propose \emph{Latent-MoE}, which interleaves domain-aware MoE blocks within a shared backbone. Unlike FB-PINNs or X-PINNs, which rigidly partition both the domain and the parameters so that the parameters on different subdomains are updated independently, Latent-MoE is designed to preserve the localization benefit of domain-aware routing while allowing capacity to flow across regions through the shared backbone. On standard homogeneous-physics benchmarks Latent-MoE is competitive with established baselines; on benchmarks with multi-stage time-variable physics, where global models and rigid domain decompositions both fall into spurious solutions, it improves over them by more than an order of magnitude, with markedly reduced gradient conflict during training.
Classical crystallography determines orientation and unit-cell geometry from a rotation series collected from a single specimen, whereas serial crystallography acquires one still exposure from each of many randomly oriented microcrystals and must reconstruct a complete dataset by pooling partial measurements across the ensemble. In serial experiments, each frame typically contains only a subset of reflections, and frames with few reliable peaks (the small-N regime) are especially difficult to index blindly, because multiple incorrect lattices or orientations can explain sparse observations. The approach implemented in GLINT addresses this small-N regime by treating the dataset, rather than the individual frame, as the unit of inference. Candidate solutions are generated for many frames, weak but recurrent lattice hypotheses are pooled across the ensemble to identify a shared unit cell, and each frame is then registered against that consensus cell through a known-cell orientation search with lower dimensionality than blind indexing. In benchmark tests, GLINT matched the strongest blind indexer in the comparison on indexing yield at roughly 40x lower per-frame blind-indexing latency against XGANDALF's fastest configuration, and batched known-cell registration ran entirely on one GPU, allowing cell discovery on the fly during serial data collection. On experimental serial data, frames indexed blind by GLINT merged to crystallographic quality (CC* = 0.90), while the method rejected non-crystal images rather than forcing unsupported indexing assignments. These results demonstrate that indexing performance can be evaluated not only by indexing yield and throughput, but also by the quality of the resulting merged crystallographic data.
Neural quantum states (NQS) have emerged as a powerful framework for simulating non-equilibrium dynamics in strongly correlated quantum systems, offering scalable variational representations of highly entangled states. Yet, accurate NQS simulations have been found to be surprisingly challenging in some physical regimes of limited complexity. Here, we address paradigmatic quench dynamics of a one-dimensional quantum Ising model as a controlled benchmark. Through supervised state reconstruction we establish substantially tighter empirical upper bounds on the required parameter count than previous estimates, ruling out representational limitations as the key obstruction. Instead, we uncover a geometric inflation of small deviations as a hitherto overlooked challenge for accurate solutions of the infinitesimal time-dependent variational principle (TDVP): the dynamical rotation of the kernel of the quantum geometric tensor (QGT) can suddenly lend physical significance to previously irrelevant parameter deviations. The stability of matrix product state solutions of the same TDVP suggests that the non-linearity of the neural network ansatz is the origin of the sensitivity. These results identify QGT-null-space rotation as a geometric diagnostic of sensitive NQS dynamics and as a concrete target for improving TDVP algorithms
Rare-event search experiments continue to extend their sensitivities to unprecedented levels. This requires increasingly small backgrounds, and shielding schemes that can suppress external backgrounds by several orders of magnitude. Moreover, detailed simulations are required to attain a good understanding of these experimental backgrounds. Studies have suggested, however, that simulations and computing-related tasks contribute approximately 10\% to the average particle physicist's carbon footprint. With backgrounds frequently below 0.01 counts per kg of target per keV of energy, simulations require more computing resources; growing a rare-event researcher's computing-related carbon footprint. Event biasing is often applied in shielding simulations as a remedy, yet there is a lack of systematic guidance on how to maximally benefit from them. An optimisation study for the importance-splitting biasing technique is discussed, focused on balancing statistical precision with simulation CPU-time, and how this can benefit the average researcher's carbon footprint.
Critical universality does not determine the microscopic conditional structure of a probability measure. We study a bicolored cluster recursion constrained to remain critical at every generation, with no equilibrium spin measure or fixed coupling imposed. In both two and three dimensions, the resulting history-dependent sequence develops a common Ising/Fortuin--Kasteleyn (FK) compatibility structure: the second-shell dependence of a one-site conditional law is strongly suppressed, nearest-neighbor effective couplings move progressively toward one another near the critical Ising value, and cluster and interface observables organize around the corresponding FK geometry. An exact cluster-coloring factorization singles out $q=2$ as the point where the residual connectivity weight disappears from the two-color spin marginal. Thus equilibrium-compatible conditional structure can emerge along a trajectory that remains critical throughout.
Data-driven constitutive models for path-dependent plasticity are commonly formulated using nonlinear recurrent neural networks, whose sequential state evolution limits parallel training and whose predictions may depend on the discretization of the applied strain path. We introduce a Constitutive State Space (CSS) model that reformulates structured state-space dynamics as an incremental constitutive operator. The strain increment is decomposed into magnitude and direction: the loading direction drives the latent state-space system, while the increment magnitude enters the zero-order-hold discretization of its continuous-time linear recurrence. This mechanics-tailored construction guarantees stationarity under zero increments, strongly reduces sensitivity to strain-path resolution, and retains the parallel-scan structure of S5 for efficient training on long constitutive histories. The CSS and Minimal State Cell (MSC) architectures are compared for four multiaxial path-dependent material models including isotropic J2 plasticity, pressure-sensitive foam plasticity, and combined isotropic-kinematic hardening. CSS matches or exceeds the prediction accuracy of the MSC, including one order of magnitude lower validation losses for the plastically incompressible materials. Importantly, CSS maintains low errors across large changes in strain-path discretization, whereas the MSC error increases substantially when evaluated at coarser resolutions than used for training. CSS trains substantially faster and requires fewer strain-stress pairs to attain comparable or better accuracy. Analysis of the learned state further reveals latent structure consistent with the dimensionality of the underlying physical constitutive models. These results establish mechanics-tailored structured state-space dynamics as a computational framework for efficient and discretization-robust data-driven constitutive modeling.
We consider a mixed quantum-classical formulation for open nonequilibrium molecular systems. We employ the pseudoparticle nonequilibrium Green's function (PP-NEGF) method as a fully quantum description and a starting point for introducing classical nuclear dynamics for an open system. We use this formulation to derive a density-matrix equation of motion describing the quantum-classical evolution. We compare our results with previous formulations and highlight important differences related to the nonequilibrium and open character of the molecular system. We also analyze the approximations needed to reduce full quantum-classical dynamics to the fewest-switches surface hopping (FSSH) method and show that these approximations are unreasonable when adiabatic surfaces come close to each other.
Physics-Informed Neural Networks (PINNs) embed PDE residuals into neural network training, but their internal representations remain opaque: it is unknown what physical features their hidden layers encode or whether those features have a localized causal role. We present PhysSAE, a mechanistic interpretability framework that trains overcomplete sparse autoencoders (SAEs) on PINN penultimate-layer activations and evaluates dictionary atoms through direct causal intervention in the original frozen hidden state: $h_{\mathrm{cf}} = h - αz_k d_k$, bypassing the SAE decoder entirely. Across six PDE families, with 3 PINN seeds and 3 SAE seeds each---we show that (i) Our discovered SAE atoms align with independently-defined physical observables (max Pearson $|r|=0.951$, always $\gg$ permutation null), (ii) the causal footprint of top-aligned atom ablation is 1.2--4.2$\times$ more spatially concentrated canonical than PCA or ICA interventions, and (iii) top-aligned atoms outperform matched random controls on causal localization for structured physical concepts (ESF$_{80}$ advantage 0.04-0.44). Two-atom bilateral representations improve concept regression R$^2$ by $ΔR^2\!=\!0.05\text{-}0.15$ over single atoms, while random pairs decrease it by up to 0.60. These results demonstrate that PINNs develop sparse, physically structured latent representations that can be identified and causally interrogated post-hoc, opening a path toward interpretability-aware scientific machine learning.