Symmetry-Preserving Physics-Informed Neural Network Framework for Relativistic Charged-Particle Dynamics in 3+1 Dimensions

1. Introduction

Relativistic charged-particle dynamics in intense electromagnetic fields constitutes one of the central challenges of modern plasma physics, laser–matter interaction physics, and high-energy accelerator science. The rapid development of petawatt-class laser facilities has made it experimentally feasible to accelerate electrons to GeV-scale energies over centimeter-scale distances via laser wakefield acceleration (LWFA), generating ultra-short, high-brightness particle beams unattainable with conventional radio-frequency accelerators [1,2]. In these extreme regimes—where the normalized laser vector potential satisfies $a_0=eA/m{c}^2\gg 1$—the particle dynamics become intrinsically relativistic, and even small numerical errors in the equations of motion accumulate rapidly, leading to qualitatively incorrect physical predictions. Accurate long-term integration of the relativistic Lorentz-force equation is therefore of paramount importance for the predictive simulation of such systems.

The standard workhorse of particle-in-cell (PIC) codes for solving the relativistic equations of motion is the Boris pusher [3], supplemented by the later improvements of Vay [4] and Higuera and Cary [5]. The Boris algorithm is volume preserving and second-order accurate in time, but it is not symplectic in the strict sense: it does not exactly preserve the symplectic two-form $\omega =\mathrm{d}\mathbf{P}\wedge \mathrm{d}\mathbf{q}$ of the underlying Hamiltonian system [6,7]. As a consequence, long-time integrations with the Boris pusher can accumulate a non-physical secular drift of geometric phase-space invariants, even while the gyro-orbit radius itself is preserved with remarkable accuracy [8,9]. The Vay integrator correctly preserves the $\mathbf{E}\times \mathbf{B}$ drift velocity but likewise does not satisfy strict symplecticity, and higher-order Runge–Kutta (RK4) methods, while accurate over short intervals, fail to conserve the geometric structure of phase space for nonlinear Hamiltonian systems over long integration times [10].

To overcome these limitations, structure-preserving numerical methods—symplectic integrators—have been developed, which preserve the symplectic form exactly and hence, through backward error analysis, conserve the corresponding integrals of motion up to exponentially small corrections over arbitrarily long times [10,11]. For non-relativistic charged particles, explicit high-order symplectic integrators of the Yoshida and Forest–Ruth type are well established [11,12]. For relativistic dynamics the situation is more complex, since the relativistic Hamiltonian $\mathcal{H}=c\sqrt{{(\mathbf{P}-e\mathbf{A}/c)}^2+m^2{c}^2}+e\phi$ is non-separable in general electromagnetic fields, precluding a straightforward splitting. Explicit high-order symplectic integrators for charged particles in general fields were constructed by Tao through an extended phase-space formulation [13]; structure-preserving second-order methods for relativistic trajectories were derived by Higuera and Cary [5]; volume-preserving algorithms for charged-particle dynamics were developed by He et al. [14]; and a comprehensive comparison of relativistic particle pushers was performed by Ripperda et al. [15]. Nevertheless, a framework that simultaneously guarantees (i) a Lorentz-invariant Hamiltonian structure, (ii) explicit high-order symplecticity, and (iii) applicability to prescribed, time-dependent fields such as Gaussian laser pulses remains an open problem.

In parallel, the paradigm of physics-informed neural networks (PINNs), introduced by Raissi, Perdikaris, and Karniadakis [16] (see [17] for a comprehensive review), has provided a powerful approach to embedding physical laws directly in neural-network training through residual loss functions. In Hamiltonian mechanics, Hamiltonian Neural Networks [18] and symplectic networks (SympNets) [19] learn the system Hamiltonian and enforce energy conservation by construction; subsequent work has extended these ideas to generalized and constrained Hamiltonian systems [23,24,25]. In particular, unsupervised Hamiltonian neural networks learn the Hamiltonian surrogate ${\mathcal{H}}_{\theta }(\mathbf{q},\mathbf{p})$ directly from the governing equations of motion—requiring no trajectory data [20]—which can then be advanced with a structure-preserving integrator. To date, however, these approaches have been developed for classical low-dimensional systems and have not been combined with an explicit symplectic map for relativistic dynamics governed by Lorentz-covariant Hamiltonian structure in 3+1 dimensions, which is the gap addressed here.

The present work builds directly on the Lorentz-invariant Lagrangian–Hamiltonian framework developed by the authors. In [36] we derived the Lagrangian and Hamiltonian formalisms for relativistic mechanics in terms of Lorentz-invariant evolution parameters, establishing the covariant canonical structure. In [37] we introduced a coupled-parameter hyperbolic parametrization of the Lorentz coordinates in terms of rapidity $\theta$ and Lorentz factor $\gamma$, demonstrating the mutual invertibility $\gamma \left(\theta \right)$ and $\theta \left(\gamma \right)$ as analytic integrals of motion. More recently [38], integrals of motion for a relativistic particle in 1+1 dimensions were obtained in terms of generalized inverse energy–momentum functions, and the dynamics of a relativistic particle in the field of a Gaussian laser pulse were analyzed in [39]. Invariant descriptions of the classical relativistic particle motion in 3+1 dimensions were subsequently derived in [40], and the dynamics of relativistic particles in an ion-cyclotron trap driven by an external ultrashort laser pulse were studied in [41]. The SP-PINN framework proposed here integrates these results into a unified, symmetry-preserving numerical integrator for the 3+1-dimensional problem.

The main contributions of this work are as follows:

• We formulate the SP-PINN integrator for a relativistic charged particle in a 3+1-dimensional prescribed electromagnetic field, including a Gaussian laser pulse, within the Lorentz-invariant Hamiltonian framework of [36,37].
• An unsupervised PINN learns the relativistic Hamiltonian surrogate ${\mathcal{H}}_{\theta }$ from a loss incorporating the covariant equations of motion, the mass-shell constraint $\mathcal{H}=m{c}^2\gamma$, and boundary conditions from the analytic integrals of motion.
• The learned Hamiltonian is advanced with an explicit symplectic map [11,13], preserving the Poincaré–Cartan invariant and the Lorentz-group orbit structure of phase space; promoting time to a canonical coordinate restores exact symplecticity for time-dependent fields, improving conservation of the light-front invariant by about three orders of magnitude (Appendix A.1).
• Performance is benchmarked against RK4 and the Boris pusher in terms of geometric-invariant conservation, symplecticity diagnostics, and trajectory accuracy, with all simulation code released openly [44].

We delimit the novelty explicitly. The Stage-2 geometric integrator is Tao’s explicit extended-phase-space symplectic map [13], used here unchanged, and the two-stage idea of learning a Hamiltonian surrogate and then integrating it symplectically follows the structure-preserving neural-integrator line [18,19,20]. Most directly, the unsupervised-PINN-Hamiltonian-plus-fourth-order-Yoshida combination has very recently been proposed independently by Liang et al. (SPINI [21]) for non-relativistic, separable systems (a nonlinear pendulum), and symplectic neural networks have been applied to charged-particle dynamics in electromagnetic fields by Drimalas et al. [22], who learn the symplectic map directly and report sub-gyroperiod accuracy against the Boris pusher. Relative to these, the present work is specifically relativistic and covariant: the learned object is a Lorentz-covariant Hamiltonian for the non-separable relativistic problem—where the standard separable Yoshida splitting does not apply and Tao’s extended phase space is required—demonstrated in 3+1 dimensions including a focused laser pulse.

The new elements of the present work are: (i) the formulation of the surrogate-learning loss for a relativistic, Lorentz-covariant Hamiltonian with an explicit mass-shell constraint $\mathcal{H}=m{c}^2\gamma$, which ties the learned object to the coadjoint-orbit (mass-hyperboloid) geometry of the Lorentz group; and (ii) the demonstration of the symplectic pipeline on 3+1-dimensional charged-particle dynamics in prescribed fields (including a focused Gaussian laser pulse), both with the analytic Hamiltonian—which isolates the geometric integrator—and, for the laser, with a learned vector-potential surrogate driven end-to-end through the complete two-stage map. We do not claim a new symplectic integrator; we claim a symmetry-aware way of supplying one with a relativistic Hamiltonian, and an assessment of where this currently helps and where it does not. In particular, every comparison figure uses the analytic Hamiltonian in Stage 2; the learned Stage-1 surrogate is assessed separately (Section 6), where a vector-potential light-cone formulation reproduces the demanding time-dependent 3+1D laser trajectory phase-coherently over essentially the whole interaction, leaving only the final closing of the surrogate-accuracy gap (${\epsilon }_{\theta }\lesssim {10}^{-4}$) to future work. We also flag at the outset that, in the coupled limit $\Omega \propto \Delta t^{-1}$ used to drive the binding floor to zero, the realized order of the Stage-2 map on the true trajectory is second, not fourth: although the map is built from a fourth-order Yoshida composition—and at fixed binding constant $\Omega$ integrates the extended Hamiltonian to fourth order down to an $\Omega$-dependent error floor—the residual binding error of the extended phase space then scales linearly with $\Delta t$ and caps the realized convergence order at two, so the benefit of the method is the absence of secular error growth, not a higher formal order (Section 4; Appendix A.2).

The remainder of the paper is organized as follows. Section 2 presents the covariant Lagrangian–Hamiltonian framework and the field configuration. Section 3 discusses the symplectic structure and the Poincaré–Cartan invariant. Section 4 describes the SP-PINN architecture, loss function, and integration scheme. Section 5 presents numerical results, Section 6 discusses their implications, and Section 7 concludes.

2. Relativistic Hamiltonian Formalism in 3+1 Dimensions

2.1. Covariant Lagrangian and Hamiltonian of a Charged Particle

The classical relativistic mechanics of a charged particle in an external electromagnetic field is naturally formulated within the variational Lagrangian–Hamiltonian framework, which makes the Lorentz symmetry of the equations of motion explicit [34,35]. The action must be a Lorentz scalar, so that the equations derived from $\delta S=0$ are covariant under the Poincaré group; this requirement, up to gauge freedom, determines the form of the Lagrangian. In an inertial frame parametrized by coordinate time t, the Lagrangian of a particle of mass m and charge e in a field with four-potential $A^{\mu }=(\phi /c,\mathbf{A})$ reads [33,35]

$$\mathcal{L}(\mathbf{r},\dot{\mathbf{r}},t)=-m{c}^2\sqrt{1-{\displaystyle \frac{{\dot{\mathbf{r}}}^2}{c^2}}}+{\displaystyle \frac{e}{c}}\mathbf{A}·\dot{\mathbf{r}}-e\phi ,$$

(1)

where $\gamma ={(1-{\dot{\mathbf{r}}}^2/c^2)}^{-1/2}$. The canonical momentum conjugate to $\mathbf{r}$ is

$$\mathbf{P}={\displaystyle \frac{\partial \mathcal{L}}{\partial \dot{\mathbf{r}}}}=\gamma m\dot{\mathbf{r}}+{\displaystyle \frac{e}{c}}\mathbf{A}=\mathbf{p}+{\displaystyle \frac{e}{c}}\mathbf{A},$$

(2)

where $\mathbf{p}=\gamma m\dot{\mathbf{r}}$ is the kinetic momentum. The distinction between canonical and kinetic momenta is fundamental for the Hamiltonian formulation and for the identification of conserved quantities under gauge transformations. The Legendre transformation yields

$$\mathcal{H}(\mathbf{r},\mathbf{P},t)=\mathbf{P}·\dot{\mathbf{r}}-\mathcal{L}=c\sqrt{{\left(\mathbf{P}-{\displaystyle \frac{e}{c}}\mathbf{A}\right)}^2+m^2{c}^2}+e\phi ,$$

(3)

the exact relativistic Hamiltonian, from which Hamilton’s canonical equations follow:

$$\dot{\mathbf{r}}={\displaystyle \frac{\partial \mathcal{H}}{\partial \mathbf{P}}}={\displaystyle \frac{c(\mathbf{P}-\frac{e}{c}\mathbf{A})}{\sqrt{{(\mathbf{P}-\frac{e}{c}\mathbf{A})}^2+m^2{c}^2}}},\dot{\mathbf{P}}=-{\displaystyle \frac{\partial \mathcal{H}}{\partial \mathbf{r}}}.$$

(4)

Upon substitution these reproduce the relativistic Lorentz-force equation $\mathrm{d}\mathbf{p}/\mathrm{d}t=e\left(\mathbf{E}+\frac{1}{c}\dot{\mathbf{r}}\times \mathbf{B}\right)$ [33]. Along any physical trajectory $\mathcal{H}=\gamma m{c}^2+e\phi$; in the radiation gauge ($\phi =0$) this reduces to $\mathcal{H}=\gamma m{c}^2$. Following [36], the proper time $\tau$ may be employed as the evolution parameter, rendering the canonical equations manifestly covariant first-order ordinary differential equations on the extended phase space $(\mathbf{r},\mathbf{P},t,-\mathcal{H})$.

2.2. Lorentz Symmetry and Integrals of Motion

The Lorentz-invariant structure of the system admits exact first integrals that are direct manifestations of Poincaré symmetry [36,37]. The mass-shell constraint

$$p^{\mu }{p}_{\mu }={\left({\displaystyle \frac{\mathcal{E}}{c}}\right)}^2-{\left|\mathbf{p}\right|}^2=m^2{c}^2,$$

(5)

with $\mathcal{E}=\gamma m{c}^2$ and $p^{\mu }=(\mathcal{E}/c,\mathbf{p})$, expresses the conservation of the Lorentz norm of the four-momentum and is preserved exactly along every physical trajectory. In terms of the Lorentz factor it reads ${\gamma }^2-{\left|\mathbf{p}\right|}^2/{\left(mc\right)}^2=1$ and furnishes the fundamental constraint enforced by the SP-PINN loss. In [37], a coupled-parameter hyperbolic parametrization was introduced in which the rapidity $\theta$—the additive parameter of Lorentz boosts in the Lie algebra of $\mathrm{SO}(1,3)$—plays the central role:

$$\gamma =cosh\theta ,{\displaystyle \frac{p_{\Vert }}{mc}}=sinh\theta ,$$

(6)

so that the phase-space point $(\gamma ,p_{\Vert }/mc)$ lies on the unit hyperbola ${\gamma }^2-{(p_{\Vert }/mc)}^2=1$ (the ${\mathbf{p}}_{\perp }=0$ restriction of the mass-shell), a coadjoint orbit of the Lorentz group [37,38]. The additivity of $\theta$ under boosts makes it the natural coordinate on the mass hyperboloid. For a field invariant under transverse translations—such as a plane wave or, approximately, a broad laser beam propagating along z—the transverse canonical momenta $P_x$ and $P_y$ are exact integrals of motion by Noether’s theorem, and together with $\theta$ form a complete set of first integrals used below to define the PINN boundary conditions.

2.3. Electromagnetic Field Configuration: Gaussian Laser Pulse

As a physically motivated and numerically demanding test case we consider a focused, linearly polarized Gaussian laser pulse propagating along z, consistent with the configuration analyzed in [39,43]. In the paraxial approximation, the vector potential $\mathbf{A}=A_x\widehat{\mathbf{x}}$ takes the form

$$A_x=A_0{\displaystyle \frac{w_0}{w\left(z\right)}}exp\left[-{\displaystyle \frac{x^2+y^2}{w{\left(z\right)}^2}}-{\displaystyle \frac{{(z-ct)}^2}{c^2{\tau }_L^2}}\right]cos\left[k_0(z-ct)+\psi (x,y,z)\right],$$

(7)

with $A_0=m{c}^2{a}_0/e$, beam waist $w_0$, pulse duration ${\tau }_L$, carrier wavenumber $k_0={\omega }_0/c$, and

$$w\left(z\right)=w_0\sqrt{1+{(z/z_R)}^2},z_R={\displaystyle \frac{k_0{w}_0^2}{2}},\psi ={\displaystyle \frac{k_0(x^2+y^2)}{2R\left(z\right)}}-arctan{\displaystyle \frac{z}{z_R}},$$

(8)

where $z_R$ is the Rayleigh length. The fields follow from $\mathbf{E}=-(1/c){\partial }_t\mathbf{A}$ and $\mathbf{B}=\nabla \times \mathbf{A}$. In the plane-wave limit ($w_0\to \infty$) the field is invariant under transverse translations, implying exact conservation of $P_x$ and $P_y$; for a finite beam this symmetry is weakly broken by diffraction, while a discrete ${\mathbb{Z}}_2$ symmetry under $x\to -x$ is retained. In the regime $a_0\gg 1$ the dynamics become strongly relativistic—the transverse quiver Lorentz factor is of order $a_0$, while an initially stationary electron acquires a ponderomotive energy gain $\langle \gamma \rangle -1\sim a_0^2/4$ (cycle-averaged), with instantaneous peaks $\gamma \approx 12$ for $a_0=5$ (Figure 3)—making long-time structure preservation particularly important. The parameters used in Section 5 are summarized in Table 1.

3. Symplectic Structure and the Poincaré–Cartan Invariant

3.1. Symplectic Form in Relativistic Phase Space

The geometric foundation of Hamiltonian mechanics is the symplectic structure of phase space [30,31,32]. For canonical coordinates $(\mathbf{q},\mathbf{P})\in {\mathbb{R}}^{2n}$ the canonical symplectic two-form $\omega ={\sum }_i\mathrm{d}{P}_i\wedge \mathrm{d}{q}^i$ is closed and non-degenerate, and Hamilton’s equations are equivalent to ${\iota }_{X_{\mathcal{H}}}\omega =\mathrm{d}\mathcal{H}$. The flow generated by any Hamiltonian vector field is a symplectomorphism: it preserves $\omega$ exactly. The Poincaré–Cartan integral invariant is built from the canonical one-form $\Theta =P_{\mu }\mathrm{d}{q}^{\mu }$,

$$I_1={\oint }_{\mathcal{C}\left(t\right)}\Theta ={\oint }_{\mathcal{C}\left(t\right)}\left(\mathbf{P}·\mathrm{d}\mathbf{q}-\mathcal{H}\mathrm{d}t\right),$$

(9)

and is conserved exactly along the flow for any closed co-moving curve $\mathcal{C}\left(t\right)$ [30]. At fixed time it reduces to the first Poincaré invariant $I_1={\oint }_{\mathcal{C}}\mathbf{P}·\mathrm{d}\mathbf{q}=\mathrm{const}$, which provides the central computable diagnostic of symplecticity. Liouville’s theorem—preservation of the phase-space volume ${\omega }^n/n!$—is a corollary of symplecticity, but is strictly weaker: a map can preserve volume without preserving $\omega$ [7,9].

3.2. Violation of Symplecticity in Conventional Solvers

Standard integrators used in plasma physics do not, in general, preserve the symplectic form. The Boris pusher is volume preserving but not symplectic [7]; Qin et al. [7] showed that it conserves phase-space volume while its energy error is bounded for integrable systems but can grow secularly for non-integrable ones. RK4 is neither symplectic nor volume preserving: its per-step deviation from a symplectomorphism is $\mathcal{O}\left(\Delta t^5\right)$, and because it admits no modified Hamiltonian the geometric invariants drift secularly (linearly in the number of steps) rather than remaining bounded [10]. The first Poincaré invariant is monitored numerically by evolving an ensemble of M nearby phase-space points forming a closed loop $\mathcal{C}$ and evaluating

$$I_1^{\left(k\right)}={\displaystyle \sum _{j=1}^M}{\displaystyle \frac{{\mathbf{P}}_j^{\left(k\right)}+{\mathbf{P}}_{j+1}^{\left(k\right)}}{2}}·\left({\mathbf{q}}_{j+1}^{\left(k\right)}-{\mathbf{q}}_j^{\left(k\right)}\right),{\mathcal{C}}_{M+1}\equiv {\mathcal{C}}_1,$$

(10)

at each step k [28,29]. For a truly symplectic integrator $I_1^{\left(k\right)}=I_1^{\left(0\right)}$ up to the discretization error of the loop quadrature. The relative deviation $\delta I_1\left(t\right)=|I_1\left(t\right)-I_1\left(0\right)|/|I_1\left(0\right)|$ provides a rigorous symplecticity diagnostic; for the integrable benchmarks considered below, the bounded versus secular growth of the exact scalar invariants (the Lorentz factor and the Larmor radius) is a robust signature of the same structure-preserving-versus-non-structure-preserving distinction—necessary though, since a volume-preserving non-symplectic map such as the Boris pusher also keeps them bounded, not by itself sufficient to certify symplecticity—and is the quantity we actually report in Section 5, so the reported evidence does not depend on the loop-discretization size M.

4. SP-PINN: Proposed Method

4.1. Architecture Overview

The SP-PINN integrator is a two-stage hybrid algorithm that decouples system identification from geometric integration [20]. In Stage 1 an unsupervised PINN ${\mathcal{H}}_{\theta }(\mathbf{r},\mathbf{P};t)$ is trained on collocation points in the extended phase space, with a loss combining the residual of Hamilton’s equations, the relativistic mass-shell constraint, and boundary conditions from the analytic integrals of motion [16,36,37]. In Stage 2 the trained surrogate is advanced with an explicit symplectic map. Because the symplectic map is a composition of exact canonical sub-flows, the generated discrete map is exactly symplectic with respect to ${\mathcal{H}}_{\theta }$, independently of the training error [11,13]; it is the fidelity of ${\mathcal{H}}_{\theta }$ to the true Hamiltonian—not the symplecticity of the map—that is governed by ${\epsilon }_{\theta }$, which can be reduced by improving the surrogate parametrization and collocation density. The workflow is illustrated in Figure 1; the dashed boundary separates the offline training phase from the online integration phase, the latter requiring only forward evaluations of ${\mathcal{H}}_{\theta }$ and its gradients via automatic differentiation.

4.2. Stage 1—Unsupervised PINN for Hamiltonian Learning

The surrogate ${\mathcal{H}}_{\theta }:{\mathbb{R}}^7\to \mathbb{R}$ is a fully connected feedforward network with tanh activations,

$${\mathcal{H}}_{\theta }(\mathbf{r},\mathbf{P},t)=W^{(L+1)}\sigma \left(\dots \sigma \left(W^{\left(1\right)}[\mathbf{r};\mathbf{P};t]+{\mathbf{b}}^{\left(1\right)}\right)\dots \right)+b^{(L+1)},$$

(11)

with $\sigma =tanh$. The smoothness of tanh guarantees that all derivatives entering the loss are well defined and computable by automatic differentiation [16]. Weights are initialized with the Xavier scheme. Concretely, the static (free-particle, magnetic-field, and non-integrable) surrogates use $L=5$ hidden layers of width 128, whereas the time-dependent vector-potential laser surrogate of Section 6 uses $L=6$ layers of width 256 ($\approx 3.4\times {10}^5$ trainable parameters); these and the remaining hyperparameters are specified in the released code [44]. The total objective is

$${\mathcal{L}}_{\mathrm{total}}={\mathcal{L}}_{\mathrm{eqs}}+{\lambda }_1{\mathcal{L}}_{\mathrm{constraint}}+{\lambda }_2{\mathcal{L}}_{\mathrm{bc}},$$

(12)

with the residual of Hamilton’s equations

$${\mathcal{L}}_{\mathrm{eqs}}={\displaystyle \frac{1}{N_c}}\sum _k\left[{∥{\dot{\mathbf{r}}}_k-{\displaystyle \frac{\partial {\mathcal{H}}_{\theta }}{\partial \mathbf{P}}}{|}_k∥}^2+{∥{\dot{\mathbf{P}}}_k+{\displaystyle \frac{\partial {\mathcal{H}}_{\theta }}{\partial \mathbf{r}}}{|}_k∥}^2\right],$$

(13)

the relativistic mass-shell constraint

$${\mathcal{L}}_{\mathrm{constraint}}={\displaystyle \frac{1}{N_c}}\sum _k{\left({\mathcal{H}}_{\theta }({\mathbf{r}}_k,{\mathbf{P}}_k,t_k)-m{c}^2{\gamma }_k\right)}^2,{\gamma }_k=\sqrt{1+|{\mathbf{p}}_k{|}^2/{\left(mc\right)}^2},$$

(14)

and boundary conditions ${\mathcal{L}}_{\mathrm{bc}}$ enforcing the analytic integrals of motion (initial energy, rapidity, and—for the laser case—the conserved transverse momenta) at the initial collocation points [37,38]. The mass-shell term is the central symmetry-preserving ingredient: enforcing it throughout training biases ${\mathcal{H}}_{\theta }$ toward the hyperbolic mass-hyperboloid geometry of relativistic momentum space, the coadjoint orbit of the Lorentz group. (Here ${\mathbf{p}}_k={\mathbf{P}}_k-\frac{e}{c}\mathbf{A}({\mathbf{r}}_k,t_k)$ is the kinetic momentum, formed from the canonical network input and the known field. This $\phi =0$ form applies to the magnetic-field and laser cases trained here; for a static scalar potential the surrogate target carries the additional $e\phi$ term, ${\mathcal{H}}_{\theta }=\gamma m{c}^2+e\phi$, and the constraint is imposed on the kinetic mass-shell.) Collocation points are drawn quasi-randomly (Sobol sequence) over the physically relevant domain; optimization proceeds with Adam followed by L-BFGS, with ${\lambda }_1=10$, ${\lambda }_2=1$ following the adaptive-weighting practice of [24,26,27]. We quantify a trained surrogate by ${\epsilon }_{\theta }$, the root-mean-square error of ${\mathcal{H}}_{\theta }$ (or, for the vector-potential variant of Section 4.4, of ${\mathbf{A}}_{\theta }$) together with its gradients, evaluated over a held-out test set.

4.3. Stage 2—Explicit Symplectic Integration

Once ${\mathcal{H}}_{\theta }$ has converged it is treated as a fixed, differentiable Hamiltonian and advanced with an explicit symplectic integrator. Because the relativistic Hamiltonian (3) is non-separable, we employ the extended phase-space construction of Tao [13]: two copies $(\mathbf{r},\mathbf{P})$ and $(\mathbf{x},\mathbf{y})$ are introduced with the extended Hamiltonian

$$\overline{\mathcal{H}}={\mathcal{H}}_{\theta }(\mathbf{r},\mathbf{y})+{\mathcal{H}}_{\theta }(\mathbf{x},\mathbf{P})+{\displaystyle \frac{\Omega }{2}}\left({\parallel \mathbf{r}-\mathbf{x}\parallel }^2+{\parallel \mathbf{P}-\mathbf{y}\parallel }^2\right),$$

(15)

whose three parts each admit an explicit symplectic flow; a binding constant $\Omega$ keeps the two copies synchronized (its value must be chosen with care, as the bounded-error floor depends non-monotonically on $\Omega$ and develops sharp resonances; see Appendix A.2). A symmetric Strang composition of these flows yields a second-order symplectic map ${\Phi }_{\Delta t}^{\left(2\right)}$, which is promoted to fourth order by the Yoshida triple-jump

$${\Phi }_{\Delta t}^{\left(4\right)}={\Phi }_{w_1\Delta t}^{\left(2\right)}\circ {\Phi }_{w_0\Delta t}^{\left(2\right)}\circ {\Phi }_{w_1\Delta t}^{\left(2\right)},w_1={\displaystyle \frac{1}{2-2^{1/3}}},w_0=1-2w_1,$$

(16)

with $w_1\approx 1.3512$, $w_0\approx -1.7024$ [11,12]. All gradients of ${\mathcal{H}}_{\theta }$ are obtained by automatic differentiation. For an autonomous Hamiltonian (the free-particle and magnetic-field cases) the resulting map is exactly symplectic by construction, independently of the accuracy of ${\mathcal{H}}_{\theta }$: even with a small approximation error the discrete flow is symplectic with respect to ${\mathcal{H}}_{\theta }$. While the Yoshida triple-jump is a fourth-order composition for the extended (binding) Hamiltonian, the order realized for the true trajectory is limited by the binding error of the extended phase space; we measure it directly in Appendix A.2 and find a realized second-order convergence in the coupled limit $\Omega \propto \Delta t^{-1}$. This is consistent with the central message of the paper: the value of the symplectic map is not a higher formal order—a high-order non-symplectic method is more accurate at short times—but the absence of secular error growth. For an explicitly time-dependent Hamiltonian (the Gaussian-pulse case) the field is held fixed at the step time within each composition, so the present implementation is non-autonomously frozen; the exact symplecticity and fourth-order guarantees then apply rigorously only to the static problems, while the laser results in Section 5.2 should be read as those of a high-order structure-aware scheme. Promoting time to a canonical coordinate—extending the phase space by the pair $(t,-\mathcal{H})$ as in Section 2—restores exact symplecticity for the time-dependent case; we implement this autonomized variant and validate it on a plane wave in Appendix A.1 (Figure A1), where it conserves the light-front invariant $\gamma -p_z$ roughly three orders of magnitude better than the frozen map—a three-order improvement quantified there.

4.4. Surrogate Variant for Time-Dependent Fields: The Vector-Potential Light-Cone Form

For the static (free-particle, magnetic-field, and non-integrable) problems the surrogate ${\mathcal{H}}_{\theta }$ is learned directly as above. The time-dependent 3+1D laser is substantially harder: the carrier $cos\left[k_0(z-ct)\right]$ produces tens of oscillations across the domain—the spectral bias of plain tanh networks—and the quiver is resonantly forced, so a small force error accumulates into a carrier-phase slip (analyzed in Section 6). For this case we therefore adopt a lower-dimensional, structurally exact reformulation, which is the configuration behind the learned-pipeline result of Figure 7. Rather than the seven-dimensional Hamiltonian, the network represents the four-dimensional vector potential in light-cone (light-front) form,

$${\mathbf{A}}_{\theta }(x,y,z,t)={\mathbf{A}}_{\mathrm{pw}}\left(\eta \right)+\mathrm{NN}(x,y,z,t),\eta =z-ct,$$

(17)

where the analytic plane-wave potential ${\mathbf{A}}_{\mathrm{pw}}\left(\eta \right)$ supplies the carrier exactly and the network learns only the slow focusing correction; the Hamiltonian is then reconstructed analytically as ${\mathcal{H}}_{\theta }=c\sqrt{{(\mathbf{P}-\frac{e}{c}{\mathbf{A}}_{\theta })}^2+m^2{c}^2}$, so the mass-shell constraint holds by construction and the Lorentz force is an exact function of $\partial {\mathbf{A}}_{\theta }$. Training (full hyperparameters in the released code [44]) uses collocation points sampled in a tube around a reference trajectory of the prescribed field—used only to concentrate the sampling domain, not as supervised targets, so the loss remains the unsupervised equation-residual form of Stage 1 and no trajectory data is fitted—together with Fourier features of the carrier phase $\eta$ and a stabilized optimizer—learning-rate warmup, gradient clipping, an upweighted gradient-matching term, residual-based adaptive resampling of the collocation points, and an L-BFGS polish from the best checkpoint. The resulting accuracy and its diagnosis are reported in Section 6.

4.5. Lorentz Symmetry Preservation

Proposition 1 

(Mass-shell preservation). Let ${\mathcal{H}}_{\theta }$ satisfy ${\mathcal{H}}_{\theta }=m{c}^2\gamma$ to within ${\epsilon }_{\theta }$, and let the Hamiltonian beautonomous(the free-particle, magnetic-field, and non-integrable cases; the time-dependent laser case is treated by the frozen-in-time map of Section 4 and is excluded here). Then, under the symplectic map ${\Phi }_{\Delta t}$, the invariant $p^{\mu }{p}_{\mu }={({\mathcal{H}}_{\theta }/c)}^2-{\left|\mathbf{p}\right|}^2$ is preserved with aboundederror $\mathcal{O}({\epsilon }_{\theta }+{\delta }_{\Omega }+\Delta t^r)$—where ${\delta }_{\Omega }=\mathcal{O}\left({\Omega }^{-1}\right)$ (away from the resonances of Appendix A.2) is the bounded binding-floor of the extended-phase-space map at fixed Ω, and the realized order $r=2$ is attained in the coupled limit $\Omega \propto \Delta t^{-1}$ that sends ${\delta }_{\Omega }\to 0$ (Appendix A.2)—with no secular growth up to integration times exponentially long in $1/\Delta t$ (the standard backward-error horizon for symplectic maps—a finite, though exponentially large, time, rather than all $t\ge 0$).

Proof 

(Proof sketch). The map preserves $\omega$ exactly, hence is a canonical transformation. Backward error analysis for symplectic integrators guarantees that the numerical trajectory lies on the level set of a modified Hamiltonian $\tilde{\mathcal{H}}={\mathcal{H}}_{\theta }+\mathcal{O}\left(\Delta t^r\right)$, which is conserved with a bounded, non-accumulating error over times exponentially long in $1/\Delta t$ [10]; the realized order r is examined in Appendix A.2. Hence ${\mathcal{H}}_{\theta }$, and the constraint $p^{\mu }{p}_{\mu }$ derived from it, do not drift secularly. The mass-shell constraint is enforced at training to $\mathcal{O}\left({\epsilon }_{\theta }\right)$, so $p^{\mu }{p}_{\mu }=m^2{c}^2+\mathcal{O}\left({\epsilon }_{\theta }\right)$ initially; combining the two estimates gives the stated bound [10,13]. This bounded behavior—in contrast to the linear-in-N drift of non-symplectic schemes—is exactly what Figure 2 demonstrates.    □

Proposition 2 

(Rapidity integral, longitudinal sector). In the purely longitudinal (1+1-dimensional) case ${\mathbf{p}}_{\perp }=0$, where the rapidity $\theta =\mathrm{arcsinh}(p_{\Vert }/mc)$ obeys $\gamma =cosh\theta$, the SP-PINN map satisfies $|{\theta }^n-{\theta }_{\mathrm{exact}}\left(t_n\right)|=\mathcal{O}({\epsilon }_{\theta }+\Delta t^2)$ with bounded, non-secular error.

We emphasize that the relation $\gamma =cosh\theta$ with $\theta =\mathrm{arcsinh}(p_{\Vert }/mc)$ holds only in the longitudinal sector; in the general 3+1-dimensional case with transverse momentum the genuinely Lorentz-invariant statement is the mass-shell constraint $p^{\mu }{p}_{\mu }=m^2{c}^2$ of Proposition 1, of which $\gamma =cosh\theta$ is the ${\mathbf{p}}_{\perp }=0$ restriction. Proposition 2 then follows from Proposition 1 in that sector. Together the two propositions establish that the SP-PINN integrator preserves the Lorentz-group orbit structure—the mass hyperboloid $p^{\mu }{p}_{\mu }=m^2{c}^2$—to within errors controlled by the training tolerance ${\epsilon }_{\theta }$ and the time step, in contrast to the unbounded constraint drift of non-symplectic schemes.

5. Numerical Experiments

Throughout this section we work in code units ($c=m=k_0={\omega }_0=1$; cf. Table 1), whereas the analytical development of Section 2, Section 3 and Section 4 is written in Gaussian (CGS) units. All experiments use Python with NumPy, SciPy, and PyTorch in double precision; the complete code is openly available [44] (see Data Availability). The Boris pusher follows the standard relativistic algorithm [3,6] and RK4 is the classical explicit scheme. Throughout this section the curve labeled “SP-PINN” is the Stage-2 explicit symplectic map driven by the analytic relativistic Hamiltonian—that is, the ${\epsilon }_{\theta }\to 0$ idealization of the learned pipeline—so that the figures isolate the geometric properties of the integrator from neural-network approximation error. The realistic training-error floor ${\epsilon }_{\theta }$ of the actual Stage-1 surrogate is reported in Section 5.4; no figure here uses the learned surrogate, because at the presently attainable ${\epsilon }_{\theta }\sim 4\times {10}^{-2}$ (static magnetic case) it is not competitive (Section 5.4). A DOP853 (eighth-order) solver with tight tolerances provides the reference trajectory where an analytic solution is unavailable (for the laser ensemble spectrum, a fine-step integration is used as reference, as noted in Figure 4). For the laser case the fields $\mathbf{E}=-\frac{1}{c}{\partial }_t\mathbf{A}$ and $\mathbf{B}=\nabla \times \mathbf{A}$ are evaluated by central finite differences of the analytic vector potential, common to all schemes; the resulting field error ($\sim {10}^{-8}$ at the step $h={10}^{-4}$ used) is well below the integrator differences reported in Figure 3Figure 4.

5.1. Test Case 1: Uniform Magnetic Field

The first case is a relativistic particle in a uniform field $\mathbf{B}=B_0\widehat{\mathbf{z}}$, the prototype integrable relativistic Hamiltonian system. The magnetic force does no work, so the Lorentz factor and the Larmor radius $r_L=p_{\perp }/\left(e{B}_0\right)$ are exact constants of motion; the cyclotron frequency is ${\omega }_c=e{B}_0/\left(\gamma mc\right)$. We take ${\gamma }_0=5$ (so $p_{\perp }=\sqrt{24}mc$), $p_{\Vert }=0$, and integrate for $4\times {10}^3$ cyclotron periods at $\Delta t=T_c/100$ (the SP-PINN map uses a binding constant $\Omega =15$). At this fixed$\Omega$ the SP-PINN error saturates at a bounded floor rather than converging with $\Delta t$; the realized second-order convergence of the map (Proposition 1) is a distinct property of the coupled limit $\Omega \propto \Delta t^{-1}$, examined separately in Appendix A.2. Figure 2 reports the relative error of the Larmor radius and of the Lorentz factor as functions of the number of gyrations.

Figure 2 is the textbook contrast between bounded and secular error growth. RK4, being neither symplectic nor volume preserving, accumulates a secular error in both the Larmor radius and the Lorentz factor that grows linearly with the number of gyrations n (the per-step amplitude defect of a non-symplectic Runge–Kutta map applied to rotational motion accumulates coherently, producing the observed $\propto n$ drift). The Boris pusher conserves both invariants to machine precision: for pure gyration it is an exactly volume-preserving rotation, which is its celebrated property and makes it optimal for this special integrable case [7,8]. The symplectic SP-PINN map keeps the error bounded ($\lesssim 3\times {10}^{-5}$) for all time at the level set by the binding tolerance of the extended-phase-space construction; it does not match the machine-precision conservation of Boris on this linear problem, but unlike RK4 it exhibits no secular drift and therefore overtakes RK4 within a few hundred gyrations, ending a few times below it at $n=4\times {10}^3$. We emphasize what this test does and does not show: on a benign, linear, integrable benchmark a high-order method such as RK4 is highly accurate at short times, and the Boris pusher is unbeatable for the magnetic sub-problem; the distinctive value of the symplectic SP-PINN scheme is the absence of secular drift together with its applicability to the general non-separable relativistic Hamiltonian, for which exact volume-preserving rotations are not available [9,13,30]. We return to this point in Section 6.

5.2. Test Case 2: Gaussian Laser Pulse (3+1D)

The most demanding case is the 3+1-dimensional motion of an electron in the focused Gaussian pulse of Table 1 ($a_0=5$; the Stage-2 map uses a binding constant $\Omega =40$, chosen in a low-error window away from the $\Omega =20,50$ resonances of Appendix A.2). The pulse, centered at $z=ct$, overtakes an electron initially at rest at the origin; the interaction window is $\left|t\right|\lesssim {\tau }_L$. As discussed in Section 4, for this time-dependent field the present Stage-2 map is frozen-in-time within each step, so it is here a high-order structure-aware scheme rather than a strictly symplectic one. Consequently one should not expect SP-PINN to achieve the smallest energy error in this time-dependent case: panel (c) of Figure 3 accordingly shows it lying between RK4 and Boris rather than below both, which is the expected behavior of a frozen (non-strictly-symplectic) map and not a deficiency of the integrator. Figure 3 compares the longitudinal coordinate $z\left(t\right)$, the Lorentz factor $\gamma \left(t\right)$, and the post-interaction energy error against the DOP853 reference.

All three integrators reproduce the qualitative dynamics—transverse quiver at the carrier frequency, longitudinal ponderomotive acceleration, and partial post-pulse deceleration—and track the reference closely over the moderate interaction time, with RK4 attaining the smallest absolute energy error and SP-PINN remaining within a few $\times {10}^{-5}$ of the reference. Figure 4 examines an ensemble of electrons and the preservation of the transverse translational symmetry. Panel (a) shows the final kinetic-energy spectrum $\mathrm{d}N/\mathrm{d}\gamma$ for 256 electrons with initial transverse positions uniformly distributed in ${[-w_0,w_0]}^2$; all schemes reproduce the reference spectrum closely (mean-energy error below ${10}^{-4}$). Panel (b) tracks the conserved transverse canonical momentum $P_x=p_x+(e/c)A_x$—the Noether invariant of transverse translations, exactly conserved in the plane-wave limit—for an on-axis electron, plotting its deviation from the reference. Here an apparent separation appears: the Boris pusher exhibits the largest deviation of this Lorentz-related quantity (of order ${10}^{-2}$), whereas the structure-aware schemes—high-order RK4 and the SP-PINN map—track it several orders of magnitude better. We note a caveat: for the finite focused beam used here ($w_0=5{\lambda }_0$), $P_x$ is conserved only approximately—diffraction weakly breaks the transverse translational symmetry—so panel (b) measures the deviation from the reference trajectory, which conflates a small physical symmetry breaking with numerical error; a strict plane-wave control ($w_0\to \infty$), in which $P_x$ is exactly conserved, isolates the purely numerical violation, and we provide exactly this control in Appendix A.1 (Figure A1b), where the Boris pusher again shows by far the largest, purely numerical, deviation—confirming that the ordering seen here is a numerical effect rather than the residual physical symmetry breaking of the finite beam.

5.3. Test Case 3: A Non-Integrable System

The uniform-field test of Section 5.1 is integrable, and there the Boris pusher is unbeatable; the regime in which symplecticity is genuinely decisive is non-integrable dynamics, where no exact volume-preserving rotation exists. To exhibit this we consider an autonomous but non-integrable system: a charged particle in a uniform magnetic field $\mathbf{B}=B_0\widehat{\mathbf{z}}$ plus a static anharmonic electrostatic well

$$\phi \left(\mathbf{r}\right)=\frac{1}{2}{k}_{\perp }(x^2+y^2)+ϵx^2{y}^2.$$

(18)

The coupling $ϵx^2{y}^2$ (a Hénon–Heiles-type term) breaks integrability, so the motion is chaotic at moderate energy (Figure 5a), while the total energy $\mathcal{H}=\gamma m{c}^2+e\phi$ is an exact constant of motion. We take $B_0=1$, $k_{\perp }=1$, $ϵ=0.3$, and integrate for $3\times {10}^5$ steps at $\Delta t=0.05$ (the SP-PINN map uses $\Omega =10$). Figure 5b shows the relative energy error: RK4 drifts secularly, the Boris pusher—volume-preserving but not symplectic, and tuned for pure magnetic rotation rather than the nonlinear, inhomogeneous force of the $ϵx^2{y}^2$ coupling, so that (consistent with [7]) its energy error, bounded for the integrable gyration of Section 5.1, grows secularly here—exhibits the largest error, and the symplectic SP-PINN map keeps the energy error bounded for all time. This is the qualitative advantage promised by backward error analysis (Proposition 1), now realized on a system where, unlike the pure gyration of Section 5.1, no scheme can fall back on an exact volume-preserving rotation.

5.4. Computational Cost

Figure 6 reports the total wall-clock time per integration step as a function of the number of particles integrated in parallel (batched NumPy), measured on the magnetic-field problem. Per single particle (mean ± standard deviation over five trials), Boris is cheapest ($0.108\pm 0.003$ ms), RK4 intermediate ($0.217\pm 0.004$ ms), and the SP-PINN symplectic map most expensive ($0.726\pm 0.004$ ms), because each step evaluates the Hamiltonian gradients several times within the Yoshida composition. The total per-step time grows sublinearly with the batch size, so the per-particle cost decreases under batching—to $\approx 2.2\times {10}^{-3}$ ms for SP-PINN at $N_p={10}^3$—but the SP-PINN/Boris cost ratio stays roughly constant ($\sim 3$–$7\times$) across all batch sizes; the symplectic map is, and remains, several times more expensive per particle than Boris. These figures are single-CPU; on a GPU the ratio is expected to narrow appreciably, because the dominant per-step cost—the several automatic-differentiation evaluations of $\nabla {\mathcal{H}}_{\theta }$ inside the Yoshida composition—maps directly onto the highly parallel batched tensor kernels for which PyTorch is optimized, whereas the Boris rotation is already arithmetically minimal and benefits less from massive parallelism. For the PIC-style workloads that motivate this work—many particles pushed simultaneously—a GPU implementation should therefore reduce the relative overhead of the symplectic map below the single-CPU ratio reported here; we leave a dedicated GPU benchmark to future work. The one-time Stage-1 training of the PINN surrogate is a fixed offline cost. Training the compact network of Section 4 on the magnetic-field Hamiltonian over the domain $\mathbf{r},\mathbf{P}\in {[-6,6]}^3$ for about $3.5$ min per seed on a single CPU (an Intel Core i7-4770; 2500 Adam epochs followed by 80 L-BFGS iterations) reaches, over five random seeds, a root-mean-square test error of only ${\epsilon }_{\theta }=(4.1\pm 0.5)\times {10}^{-2}$ in ${\mathcal{H}}_{\theta }$ and its gradients, with a larger worst-case error ($\sim 0.3$) near the domain boundary. This surrogate accuracy—not the geometric integrator—is the present bottleneck of the fully learned pipeline: at ${\epsilon }_{\theta }\sim 4\times {10}^{-2}$ the learned-Hamiltonian integration does not yet match the conventional pushers, and the headline curves of Figure 2, Figure 3 and Figure 4 accordingly use the analytic Hamiltonian, i.e. the ${\epsilon }_{\theta }\to 0$ limit. Driving ${\epsilon }_{\theta }$ down for the full 3+1D surrogate—primarily by reformulating what the network represents (the vector-potential light-cone form of Section 6), together with adaptive collocation and improved optimization rather than raw network size—is the principal direction of ongoing work. Representative timings are collected in Table 2.

6. Discussion

6.1. Symplecticity Versus Volume Preservation

The central lesson of Section 5 is the distinction between bounded and secular error growth, and the complementary roles of the three integrators. A common simplification in the particle-pushing literature equates energy conservation, volume preservation, and symplecticity [7]; the magnetic-field test makes the distinctions concrete. The Boris pusher is exactly volume preserving for pure gyration and conserves the Larmor radius and Lorentz factor to machine precision, making it the method of choice for the magnetic sub-problem; it is, however, not symplectic in the strict sense, and its favorable behavior does not automatically extend to general non-separable fields [7,9]. RK4 is highly accurate at short times but drifts secularly, because its discrete map is not structure preserving. The symplectic SP-PINN map, through backward error analysis, keeps the error bounded for all time and therefore overtakes RK4 at long integration times. We deliberately avoid overstating the magnitude of the advantage on this benign, integrable benchmark: a high-order non-symplectic method is competitive here, and the symplectic gain is qualitative (bounded versus secular) rather than a dramatic short-time accuracy improvement. The regime in which symplecticity becomes decisive is long-time integration and non-integrable dynamics in complex fields, where no exact volume-preserving rotation is available and where secular drift of non-symplectic schemes dominates the error budget [9,10,28]. The principal practical contribution of SP-PINN is to make explicit, structure-preserving integration available for the general non-separable relativistic Hamiltonian through the learned surrogate, rather than to outperform Boris on the special magnetic case.

6.2. Role of the Mass-Shell Constraint

The mass-shell term ${\mathcal{L}}_{\mathrm{constraint}}={({\mathcal{H}}_{\theta }-m{c}^2\gamma )}^2$, weighted by ${\lambda }_1\gg 1$, is the single most important design choice for Lorentz-symmetry preservation in Stage 1. Enforcing it throughout training biases the learned Hamiltonian toward the mass-hyperboloid geometry; omitting it allows the surrogate to reproduce the equations of motion while letting $p^{\mu }{p}_{\mu }$ drift. This is consistent with the broader experience that high-weight constraint residuals improve the geometric fidelity of constrained-Hamiltonian PINNs [24,26].

6.3. Limitations

Several limitations deserve emphasis. First, and most importantly, the surrogate ${\mathcal{H}}_{\theta }$ has so far been trained only for the static magnetic-field Hamiltonian, where it reaches only ${\epsilon }_{\theta }=(4.1\pm 0.5)\times {10}^{-2}$ (root-mean-square over the training domain, five seeds) with about $3.5$ min of single-CPU training per seed; at this accuracy the learned pipeline is not competitive with conventional pushers, which is why every figure here uses the analytic Hamiltonian in Stage 2. Training a surrogate for the time-dependent 3+1D laser Hamiltonian is markedly harder still, and exposes a sharper obstacle: a plain tanh network over the full phase-space–time box fails to learn it, because the laser carrier $cos\left[k_0(z-ct)\right]$ produces tens of oscillations across the domain that the network cannot represent—the well-known spectral bias of multilayer perceptrons. With a plain network the loss stalls near ${\epsilon }_{\theta }\approx 0.8$ and the resulting Hamiltonian yields trajectories whose amplitude is wrong by a factor of two. This obstacle is substantially mitigated by (i) sampling collocation points in a tube around the integrated trajectory, concentrating capacity where the integrator evaluates ${\mathcal{H}}_{\theta }$, and (ii) augmenting the input with Fourier features of the carrier phase $\eta =z-ct$. On a single GPU these two changes restore healthy training (the loss falls by four orders of magnitude rather than stalling) and reduce the in-tube error to ${\epsilon }_{\theta }\approx 9\times {10}^{-3}$ (root-mean-square; worst case $\sim 0.7$ near the tube boundary), recovering the correct trajectory amplitude (peak $\gamma \approx 12$); the learned trajectory tracks the exact one closely early in the interaction.

However, a further, more demanding obstacle remains, and we have isolated its origin. The laser-driven quiver is a resonantly forced oscillation, so the ∼1% residual error in the learned force accumulates into a carrier-phase slip; the integrated trajectory keeps the right amplitude but loses phase coherence with the exact solution as the interaction proceeds. Crucially, integrating the learned Hamiltonian with the symplectic Stage-2 map gives the same phase drift as a non-symplectic Runge–Kutta integration: the two learned trajectories are indistinguishable from each other, sharing the same peak instantaneous relative error ${max}_t|{\gamma }_{\mathrm{learned}}-{\gamma }_{\mathrm{exact}}|/{\gamma }_{\mathrm{exact}}\approx 12$ (we use this dimensionless quantity throughout as the trajectory-error metric; a value of order unity corresponds to a full carrier-cycle phase slip, so $\approx 12$ signals a multi-cycle slip). The drift is therefore not an integration artifact but a direct consequence of the surrogate error—both integrators faithfully follow the slightly incorrect learned Hamiltonian. Faithful integration of this strongly resonant system thus requires the surrogate force to be far more accurate (we estimate ${\epsilon }_{\theta }\lesssim {10}^{-4}$) than the mass-shell fit alone provides. A substantial improvement comes from a light-cone residual formulation: the analytic plane-wave Hamiltonian—a function of the light-cone phase $\eta =z-ct$—supplies the carrier exactly, and the network learns only the slow focusing correction (a residual that is ∼0.6% of $\mathcal{H}$). This supplies the carrier phase from closed-form physics rather than from the network, and roughly halves the surrogate error (${\epsilon }_{\theta }\approx 5\times {10}^{-3}$) while substantially reducing the trajectory error (to $\approx 7$, from the $\approx 12$ of the tube/Fourier surrogate), with the learned trajectory tracking the exact one further into the interaction. It confirms the light-cone idea as the right direction, but does not by itself fully close the gap. To determine whether the residual drift is a deficiency of the integrator or of the surrogate, we also recast the evolution in the light-cone phase: taking $\eta$ as the independent variable and tracking the slowly varying light-front invariants $(P_x,P_y,\gamma -p_z)$—from which $\gamma$ is reconstructed algebraically together with the analytic field $A_x\left(\eta \right)$—so that the carrier phase is carried by the exact independent variable and cannot accumulate. With the analytic Hamiltonian this scheme reproduces the reference trajectory to $0.25\%$, confirming both the derivation and that the formulation is phase-coherent by construction; with the learned Hamiltonian it nonetheless drifts essentially as before. Together with the symplectic-versus-Runge–Kutta comparison above, this isolates the bottleneck unambiguously: no reformulation of the integration—time-domain, symplectic, or light-cone—can compensate for the ∼1% surrogate force error, which the resonant quiver inevitably amplifies.

Acting directly on this diagnosis, we obtain our best result by reformulating what the network represents: instead of the seven-dimensional Hamiltonian (or its residual), we learn the four-dimensional vector potential in light-cone form (the surrogate variant specified in Section 4.4), $A_{\theta }(x,y,z,t)=A_{\mathrm{pw}}\left(\eta \right)+\mathrm{NN}(x,y,z,t)$, and reconstruct the Hamiltonian analytically as ${\mathcal{H}}_{\theta }=c\sqrt{{(\mathbf{P}-\frac{e}{c}{\mathbf{A}}_{\theta })}^2+m^2{c}^2}$. This formulation is lower-dimensional and structurally exact—the mass-shell constraint holds by construction and the Lorentz force is an exact function of the learned $\partial A_{\theta }$, so the network need only fit a smooth four-dimensional scalar. On a single GPU, with a stabilized optimizer (learning-rate warmup, gradient clipping, an upweighted gradient-matching term, residual-based adaptive resampling of the collocation points, and an L-BFGS polish from the best checkpoint), it reaches ${\epsilon }_{\theta }=(3.0\pm 0.1)\times {10}^{-4}$ (root-mean-square over the in-tube test set, mean ± standard deviation over three seeds), at a training cost of about 20 min per seed on a single NVIDIA Tesla T4 GPU—roughly $2\times {10}^4$ Adam iterations followed by a short L-BFGS polish. The learned-Hamiltonian trajectory then tracks the analytic reference phase-coherently over essentially the whole interaction, with only the final field peak marginally out of phase (Figure 7); the maximum relative $\gamma$ error falls to $\approx 0.9$ ($0.87\pm 0.19$ over three seeds—the mean below unity and the worst seed only marginally above, so the learned and analytic trajectories stay within about one carrier cycle), from $\approx 12$ for the tube/Fourier surrogate and $\approx 7$ for the Hamiltonian residual—and this worst-case figure is dominated by the sharp $\gamma$-troughs, where a sub-cycle phase offset produces a large instantaneous relative error even though the envelope is tracked throughout. That lowering the surrogate force error—by changing what the network represents, with the integrator left untouched—is precisely what restores trajectory fidelity is the clearest confirmation of the diagnosis above; we also find that naively enlarging the network and prolonging training does not help and can destabilize the fit—starting from the same $3.4\times {10}^5$-parameter network at its pre-stabilization accuracy (${\epsilon }_{\theta }\approx 1.2\times {10}^{-3}$, before the learning-rate warmup, gradient-clipping, and adaptive-sampling fixes that yield the headline $(3.0\pm 0.1)\times {10}^{-4}$), doubling the parameter count to $7.5\times {10}^5$ degraded the surrogate error to ${\epsilon }_{\theta }\approx 3.5\times {10}^{-2}$, the Adam optimization diverging mid-training and never recovering the earlier minimum (training log in [44])—so the effective levers are the reformulation and a stabilized optimizer, not raw capacity. Table 3 summarizes the surrogate variants and the corresponding errors.

Finally, driving this best surrogate through the actual Stage-2 Tao symplectic map—rather than the Runge–Kutta integration used to evaluate it—changes the trajectory by at most $6\times {10}^{-4}$ in $\gamma$, three orders of magnitude below the surrogate-limited error: the complete two-stage SP-PINN pipeline thus reproduces the single-integrator result, confirming directly that the residual error lies entirely in the learned surrogate and not in the geometric integrator. Reaching full quantitative reliability (${\epsilon }_{\theta }\lesssim {10}^{-4}$, now within about a factor of three), for which this vector-potential light-cone formulation is the most promising route, is therefore the single remaining prerequisite for the learned pipeline in the laser regime, and the principal direction of future work.

Second, the surrogate is trained on a bounded phase-space domain; outside it neural-network extrapolation is uncontrolled, and the Poincaré-invariant diagnostic $\delta I_1\left(t\right)$ can serve as an online detector of such excursions.

Third, the method carries hyperparameter sensitivity: the constraint weight ${\lambda }_1$, network size, learning-rate schedule, and collocation distribution all affect ${\epsilon }_{\theta }$. A distinct, integrator-level sensitivity is the binding constant $\Omega$ of the extended-phase-space map: the bounded error floor depends non-monotonically on $\Omega$ and develops sharp parametric (Floquet) resonances at $2\Omega \Delta t=k\pi$, i.e. along the comb $\Omega \approx k\pi /\left(2\Delta t\right)$ (near $\Omega \approx 20$ and $\Omega \approx 50$ here)—derived and numerically verified in Appendix A.2 (Figure A2b)—where the error degrades severely, so $\Omega$ must be placed between bands; this delicacy is a genuine practical limitation of the explicit construction and motivates an adaptive or resonance-avoiding choice of $\Omega$ for production use. The residual-based adaptive collocation refinement [27] adopted for the laser surrogate above is precisely what brought the trajectory error below unity, and a principled automated loss weighting [26] is a natural further step for production use. The supporting numerical studies—a multi-seed assessment of ${\epsilon }_{\theta }$ (Section 5.4), the realized second-order convergence of the Stage-2 map, and the dependence of the bounded error floor on the binding constant $\Omega$ (Appendix A.2)—confirm the robustness of the single-trajectory results reported in the main text.

6.4. Symmetry of the Loss and Geometric Preservation

A conceptual thread directly relevant to this journal connects the symmetry of the training objective to the geometric properties of the resulting integrator. The constraint term is a Lorentz scalar (built from the invariant $p^{\mu }{p}_{\mu }-m^2{c}^2$) and the residual term is formulated from covariant canonical equations, so the loss is Lorentz invariant by construction. Enforcing a Lorentz-covariant objective biases the network toward a Lorentz-covariant surrogate, yielding a hierarchical chain: Lorentz symmetry of the loss⇒covariant surrogate Hamiltonian⇒mass-shell conservation⇒ Lorentz-group orbit preservation under the symplectic map. In classical symplectic integration the geometric structure is guaranteed by the algorithm alone; in SP-PINN the physical symmetry of the system is encoded in the loss and transferred, through the trained network, into the geometric structure of the integrator.

6.5. Prospects

Natural extensions include radiation reaction via the Landau–Lifshitz force, which breaks Hamiltonian structure and could be accommodated by a conformal symplectic Stage 2 with an auxiliary dissipation network [42]; application to accelerated (Rindler) frames building on the hyperbolic formalism of [37]; and implementation as a drop-in Boris-pusher replacement in open-source PIC codes, with the Poincaré-invariant diagnostic acting as a self-monitoring symplecticity guard. The preservation of the Poincaré–Cartan invariant at the single-particle level is expected to improve the long-term fidelity of collective laser–plasma simulations [1,2].

7. Conclusions

We have proposed and tested SP-PINN, a two-stage symmetry-preserving framework for the relativistic equations of motion of a charged particle in 3+1-dimensional electromagnetic fields. An unsupervised physics-informed neural network learns a Lorentz-covariant surrogate Hamiltonian from a mass-shell-constrained loss, and an explicit symplectic map—built on Tao’s extended phase space and therefore applicable to the non-separable relativistic Hamiltonian—advances it while preserving the symplectic structure by construction. On a uniform magnetic field the method exhibits the bounded-versus-secular error contrast: RK4 drifts secularly in the Lorentz factor and Larmor radius, the Boris pusher conserves them to machine precision as a volume-preserving gyro-integrator, and the symplectic SP-PINN map keeps the error bounded for all time, overtaking RK4 within a few hundred gyrations. On a focused Gaussian laser pulse the method reproduces the ponderomotive dynamics, the energy spectrum, and the conserved transverse canonical momentum, tracking the high-order reference closely. The Lorentz-invariant mass-shell constraint $p^{\mu }{p}_{\mu }=m^2{c}^2$—whose longitudinal restriction is the integral $\gamma \left(\theta \right)=cosh\theta$ of [37]—is preserved to within the training and time-step tolerances, in agreement with Propositions 1 and 2. We have been explicit that on benign integrable benchmarks a high-order non-symplectic method is competitive and the Boris pusher is optimal for the magnetic sub-problem; the distinctive value of SP-PINN is the absence of secular drift combined with structure-preserving integration of the general non-separable relativistic Hamiltonian. We demonstrated this on a non-integrable magnetic trap, where—no exact volume-preserving rotation being available—the symplectic map alone keeps the energy error bounded while both RK4 and the Boris pusher drift. For the time-dependent laser surrogate, a vector-potential light-cone reformulation brings the learned-Hamiltonian trajectory into phase coherence with the reference over essentially the whole interaction (${\epsilon }_{\theta }=(3.0\pm 0.1)\times {10}^{-4}$ over three seeds), leaving the final approach to ${\epsilon }_{\theta }\lesssim {10}^{-4}$ as the principal open task. Planned extensions include radiation reaction through conformal symplectic integration, non-inertial reference frames, and integration into a full particle-in-cell loop.

Appendix A.1. Time as a Canonical Coordinate; Plane-Wave Symmetry Tests

For a time-dependent field, the frozen-in-time Stage-2 map of Section 4 is no longer exactly symplectic. Exact symplecticity is restored by adjoining time as a canonical coordinate: the extended position is $(\mathbf{r},t)$ and the extended momentum $(\mathbf{P},p_t)$, with the autonomous extended Hamiltonian $\overline{\mathcal{H}}=\mathcal{H}(\mathbf{r},\mathbf{P},t)+p_t$, so that $\mathrm{d}t/\mathrm{d}\tau =1$ and $\mathrm{d}{p}_t/\mathrm{d}\tau =-\partial \mathcal{H}/\partial t$; the Tao map is then applied in the enlarged phase space. We validate this on a linearly polarized plane wave $\mathbf{A}=a_{0}cos\left[k_0(z-ct)\right]\widehat{\mathbf{x}}$ ($a_0=2$) driving an electron from rest over ${10}^3$ wave periods (Figure A1). The plane wave has two exact invariants: the light-front (equivalently light-cone) quantity $K=\gamma -p_z$ (sensitive to the explicit time dependence) and the transverse canonical momentum $P_x$ (the Noether invariant of transverse translations). Panel (a) shows that the frozen map drifts in K to $\sim 2\times {10}^{-2}$, whereas the autonomized (time-as-coordinate) map keeps $|K-1|\lesssim 2\times {10}^{-5}$—about three orders of magnitude better—confirming the construction. Panel (b) addresses the interpretation of Figure 4b: for a true plane wave $P_x$ is exactly conserved, so any deviation is purely numerical, and the Boris pusher again shows by far the largest violation ($\sim 4\times {10}^{-2}$), confirming that the ordering seen for the focused beam is a numerical effect of the rotation splitting, not the residual physical symmetry breaking of the finite beam.

Figure A1. Plane-wave symmetry tests over ${10}^3$ wave periods ($a_0=2$). (a) Error of the light-front invariant $K=\gamma -p_z$: the frozen Tao map (green dotted) drifts to $\sim 2\times {10}^{-2}$, while the autonomized time-as-coordinate map (green solid) stays bounded at $\sim {10}^{-5}$. (b) Error of the transverse canonical momentum $P_x$, exactly conserved for a plane wave: the Boris pusher (orange) shows the largest, purely numerical, violation.

Figure A1. Plane-wave symmetry tests over ${10}^3$ wave periods ($a_0=2$). (a) Error of the light-front invariant $K=\gamma -p_z$: the frozen Tao map (green dotted) drifts to $\sim 2\times {10}^{-2}$, while the autonomized time-as-coordinate map (green solid) stays bounded at $\sim {10}^{-5}$. (b) Error of the transverse canonical momentum $P_x$, exactly conserved for a plane wave: the Boris pusher (orange) shows the largest, purely numerical, violation.

Appendix A.2. Order of Convergence and the Binding Constant

Figure A2a reports the global trajectory error at a fixed short final time $T=4$ (short enough that the error stays in the convergent regime rather than saturating at the orbit scale). RK4 exhibits its expected fourth-order slope. The Tao symplectic map, integrated in the coupled limit $\Omega =5/\Delta t$ so that the binding floor vanishes with $\Delta t$, converges at second order: although the Yoshida triple-jump is a fourth-order composition for the extended Hamiltonian, the residual binding error of the extended phase space scales linearly with $\Delta t$ in this limit and caps the realized order at two. This quantifies the trade-off at the heart of the method—the symplectic map is of lower formal order than RK4, and is accordingly less accurate at short times, but it preserves the geometric structure and therefore avoids the secular drift that dominates the error budget at long times (Figure 2 and Figure 5). Figure A2b shows the bounded long-time Larmor-error floor of the Tao map at fixed $\Delta t=T_c/100$ as a function of the binding constant $\Omega$. The dependence is non-monotonic and exhibits sharp resonances (e.g. near $\Omega =20$ and $\Omega =50$), where the binding rotation resonates parametrically with the integration step. This structure is now understood quantitatively. For the uniform field the motion is linear (the Hamiltonian is quadratic and $\gamma$ is constant along the gyro-orbit), so one step of the extended-phase-space map is a constant symplectic matrix; its monodromy develops a parametric (Floquet) instability whenever the binding rotation—which the construction performs at angular frequency $2\Omega$—is commensurate with the step,

$$2\Omega \Delta t=k\pi ,k=1,2,3,\dots ⟺{\Omega }_k={\displaystyle \frac{k\pi }{2\Delta t}},$$

(A1)

giving narrow instability bands at $\Omega \approx 5k$ for $\Delta t=T_c/100$. While the band positions follow directly from this commensurability, their relative strengths are organized—empirically, from the computed monodromy eigenvalues—by $kmod3$, a three-fold pattern traced to the three-substep Yoshida composition and its negative central weight $w_0$: the $k\equiv 1\left(\mathrm{mod}3\right)$ bands ($\Omega \approx 5,20,35,50$) are the most violent—hence the pronounced spikes at $\Omega \approx 20$ and $\Omega \approx 50$—whereas the $k\equiv 0\left(\mathrm{mod}3\right)$ bands ($\Omega \approx 15,30,45$) are the weakest, so $\Omega \approx 15$ (off the strong family and at a band edge) yields the lowest long-time error, and values of $\Omega$ lying off the comb (e.g. $\Omega \approx 2$) are linearly stable to machine precision. A direct parameter sweep (study_omega_resonance.py in [44]) confirms Eq. (A1): the resonance locations scale as $\Omega \propto \Delta t^{-1}$—halving or doubling $\Delta t$ doubles or halves the band spacing, collapsing all bands onto fixed $2\Omega \Delta t=k\pi$—and are independent of the cyclotron frequency ${\omega }_c$. The binding constant should therefore be kept off the strong $k\equiv 1\left(\mathrm{mod}3\right)$ bands; the value $\Omega =15$ used in Section 5.1 sits in a low-error window between teeth—clear of the strong family, with the nearest tooth being the weakest ($k\equiv 0$) one ($k=3$), whose narrow width and a small cyclotron detuning of the band centers leave the integer value effectively off-resonance, at the floor minimum. For the time-dependent laser and the nonlinear trap the linear (Floquet) analysis provides only the leading guide, so the values $\Omega =40$ and $\Omega =10$ were selected empirically in low-error windows; an adaptive or resonance-avoiding selection of $\Omega$ would be valuable in production use.

Figure A2. (a) Global trajectory error at $T=4$ versus time step $\Delta t$ for the uniform-magnetic-field problem: RK4 is fourth order, while the Tao symplectic map (in the coupled limit $\Omega =5/\Delta t$) is second order, its realized order capped by the binding error; reference slopes 4 (dashed) and 2 (dotted) are shown. (b) Bounded Larmor-error floor of the Tao map at $\Delta t=T_c/100$ versus the binding constant $\Omega$: the dependence is non-monotonic with sharp resonances (near $\Omega =20$, 50) that must be avoided, and a robust low-error minimum near $\Omega \approx 15$.

Figure A2. (a) Global trajectory error at $T=4$ versus time step $\Delta t$ for the uniform-magnetic-field problem: RK4 is fourth order, while the Tao symplectic map (in the coupled limit $\Omega =5/\Delta t$) is second order, its realized order capped by the binding error; reference slopes 4 (dashed) and 2 (dotted) are shown. (b) Bounded Larmor-error floor of the Tao map at $\Delta t=T_c/100$ versus the binding constant $\Omega$: the dependence is non-monotonic with sharp resonances (near $\Omega =20$, 50) that must be avoided, and a robust low-error minimum near $\Omega \approx 15$.

Appendix A.3. Free-Particle Baseline

As a zero-force baseline confirming the absence of implementation bias, consider a free particle ($\mathbf{E}=\mathbf{B}=0$): the canonical force vanishes identically, so the momentum is constant and every consistent integrator reproduces the exact solution. We verified that for ${\gamma }_0=10$ over $N={10}^4$ steps all three schemes conserve $\gamma$ and ${\left|\mathbf{p}\right|}^2$ to machine precision ($\lesssim {10}^{-15}$, Table A1); the discriminating tests are the state-dependent magnetic-field and laser cases of the main text (Section 5), where the force depends on the state.

Table A1. Lorentz-factor error $\Delta \gamma =|{\gamma }^{\left(N\right)}-{\gamma }_0|$ and relative momentum-invariant error after $N={10}^4$ steps for a free relativistic particle (${\gamma }_0=10$). All schemes are exact to machine precision because the force vanishes.

Integrator	$\mathit{\Delta }\mathit{\gamma }$	${\mathit{\Delta }\left|\mathit{p}\right|}^2/{\mathit{p}}_0^2$
RK4	$<{10}^{-15}$	$<{10}^{-15}$
Boris	$<{10}^{-15}$	$<{10}^{-15}$
SP-PINN	$<{10}^{-15}$	$<{10}^{-15}$

Table A1. Lorentz-factor error $\Delta \gamma =|{\gamma }^{\left(N\right)}-{\gamma }_0|$ and relative momentum-invariant error after $N={10}^4$ steps for a free relativistic particle (${\gamma }_0=10$). All schemes are exact to machine precision because the force vanishes.

Integrator	$\mathit{\Delta }\mathit{\gamma }$	${\mathit{\Delta }\left|\mathit{p}\right|}^2/{\mathit{p}}_0^2$
RK4	$<{10}^{-15}$	$<{10}^{-15}$
Boris	$<{10}^{-15}$	$<{10}^{-15}$
SP-PINN	$<{10}^{-15}$	$<{10}^{-15}$