Neural Hamiltonian Diffusions for Modeling Structured Geometric Dynamics

Dear Reviewer VyhR,

We sincerely appreciate your valuable and insightful feedback, which significantly helps us to improve the manuscript. While fully acknowledging your comments and committing to addressing each of them thoroughly, we respectfully ask you to reconsider our main contribution from the perspective of theoretical novelty and significance.

---

($**General Comment**.$) We humbly emphasize that the core contribution of our paper lies in its theoretical foundation, providing, for the first time, a rigorous and unified mathematical framework for explicitly modeling neural Hamiltonian dynamics across variant differentiable manifolds. Although each manifold individually may be familiar, our theoretical formulation is novel and fundamental: it systematically unifies stochastic Hamiltonian neural modeling in a general geometric setting, supported by rigorous theoretical guarantees (Propositions 2.2 and 3.1–3.2). Without introducing the concept of the horizontal flow, the resulting dynamics would not remain frame-consistent, making it impossible to consistently propagate the learned Hamiltonian dynamics across arbitrary differentiable manifolds.

We would deeply appreciate your reconsideration of the theoretical originality and foundational significance of our approach. Once again, thank you very much for your detailed and constructive comments, which we address carefully and comprehensively in our revised version. Below, we provide point-by-point responses to your specific questions and comments.

---

Q) ($**Originality**$) The paper contains a substantial amount of mathematical detail in Section 2, which mainly summarizes existing geometric analysis tools rather than presenting original contributions...

A) We sincerely thank the reviewer for raising this important point. With all due respect, we wish to clarify that the mathematical framework presented in our manuscript constitutes the entirely new core theoretical contribution of our work, rather than simply enumerating previously established geometric tools. Specifically, while the general stochastic Hamiltonian model in Definition 2.1 has appeared in prior studies, our explicit construction of a learnable Hamiltonian diffusion via the concept of horizontal lifting (Proposition 2.2) is completely new and has never before been proposed or utilized within machine learning. This horizontal lift is not merely technical machinery but an essential and indispensable component for consistently modeling frame-invariant Hamiltonian dynamics across arbitrary differentiable manifolds. To the best of our knowledge, our rigorous approach of explicitly embedding learnable Hamiltonian dynamics into a horizontal lifted frame bundle setting represents a groundbreaking theoretical advancement.

Indeed, we have made a careful and deliberate effort to include only the minimal essential mathematical definitions in Section 2, while moving all detailed derivations, rigorous proofs, and supplementary explanations to the mathematical background section of the appendix to enhance readability. Although we fully acknowledge that our theoretical framework involves advanced geometric concepts potentially unfamiliar to a broader machine learning audience, we emphasize that introducing these geometric constructs is indispensable for developing a rigorous and mathematically coherent framework. Such foundational rigor is essential, especially in scientific domains (e.g., theoretical physics, computational chemistry) that demand precise mathematical guarantees and unerring model consistency to support highly exacting physical experiments and measurements.

In this view, we believe that Eq. (7) is an incremental modification while it adopts essentially the same loss‑function structure as classical Hamiltonian Neural Network (HNN) formulations.

---

Q) ($**Clear Novelty**.$) Although the paper emphasizes its main contribution as integrating stochasticity into Hamiltonian systems, the core content predominantly ...

A) We sincerely thank the reviewer for raising this important issue. Upon careful consideration, we agree that the primary contribution and novelty of our work indeed lies in establishing a rigorous theoretical framework for modeling Hamiltonian dynamics on general manifolds i.e., explicitly addressing geometric issues such as curvature and the horizontal lifting of dynamics onto the frame bundle. While our approach naturally incorporates stochastic modeling, we fully acknowledge that the central novelty and emphasis of our paper is not only stochasticity itself, but rather on extending Hamiltonian dynamics to manifold settings in a geometrically principled manner. In the revised manuscript, we will carefully adjust our narrative and emphasize our primary contribution: $**Hamiltonian dynamics meets manifolds**$. This perspective shift will clarify the core novelty and help avoid potential confusion or overstatement regarding our contribution. We greatly appreciate the reviewer's insightful comment, which helps us significantly improve the clarity and positioning of our work.

---

Q) ($**Explicit Machinery**.$) For example, in Equation (7), it is unclear how Ut​ is obtained explicitly. More explanation is needed on how Ut​ is initialized; the criteria for selecting the numerical schemes for discretization are not discussed...

A) We sincerely appreciate the reviewer's helpful comment and acknowledge that the explicit procedure for sampling $U_t$ was not clearly illustrated in our original manuscript. We now briefly summarize Algorithm 2 below a simplified two-step procedure for explicitly obtaining the trajectory $U_t$:

$\quad$ Step 1) Choose the initial configuration $x_0 = (q_0, p_0) \in \mathcal{M}$ from the physical data (e.g., the observed atom positions), and attach a local orthonormal frame $E_0$ (i.e., random orthonormal matrix) by using Gram-Schmidt orthonormalization. Together, they form the initial state of horizontal flow $U_0 = (x_0, E_0)$.

$\quad$ Step 2) At each step $n$, we update our horizontal state by Euler–Heun Stratonovich discretization, $U_{n+1} \approx U_{n} + \delta U$, where $\delta U = G(U_n)[ v_x ] \epsilon$. Here, $G(U_n)$ is defined in Eq (5) and $v_x \approx \nabla H_\theta(q, p) \approx \\{\mathbf{m}, H_{\theta} \\}$ is gradient of proposed U-Net.

We fully agree with your feedback, and will explicitly integrate this simplified and more intuitive description of the sampling procedure into Algorithm 2, clarifying precisely how $U_t$ is constructed step-by-step. Additionally, we will relocate the Transformer-based neural network architecture clearly into the main body of the manuscript for better readability.

---

Q) ($**Challenges in Practice**.$) It is also not evident whether the terms appearing in the governing functions are straightforward to compute in practice or if they pose significant implementation challenges...

A) We agree that explicitly addressing the ease of implementation is vital for practical applicability. To clarify, all terms in the governing equations, such as the Hamiltonian vector field and the connection tensor (Christoffel symbols), are indeed straightforwardly computable given the manifold's local coordinates and metric tensor. For standard manifolds commonly encountered in scientific applications (e.g., spheres, torii), these terms have known, closed-form analytic expressions and thus pose no significant computational challenge. In the revised manuscript, we will explicitly include a simple and concrete example (such as dynamics on a sphere) that demonstrates precisely how each term is computed and implemented. We believe this will effectively clarify the practicality and accessibility of our methods.

---

Q4) ($**General Manifolds**.$) The paper seems to assume that the manifold is known and well-behaved...

A) Indeed, we agree that our manuscript implicitly assumes that the underlying manifold is known and well-behaved such as sphere and flat-torus. We acknowledge that this assumption, while common in geometric modeling frameworks, can be explicitly mentioned in our manuscript. In our revised manuscript, we will explicitly clarify this assumption by clearly stating the specific geometric conditions required for our theoretical results and algorithms to hold. Additionally, we will briefly discuss potential implications, challenges, and future directions regarding cases where the manifold is not fully known or does not satisfy these regularity conditions.

---

Q) ($**Comparison to the paper [1]**$)

A) The cited reference [1] provided by the reviewer is inherently restricted to flat Euclidean spaces, thus limiting its applicability to general scientific modeling (e.g., blackholes). In contrast, our work explicitly targets breaking through the theoretical limitations of existing Hamiltonian Neural Networks (HNNs) and extending their applicability to curved manifolds. Given this goal, it was inevitably necessary to introduce and briefly summarize certain essential geometric concepts. However, in the main text, we have deliberately omitted detailed mathematical derivations and have presented only the minimal set of definitions required to communicate our core contributions.

---

We would greatly appreciate it if the reviewer could kindly reconsider our manuscript's rating in light of our clarifications and the proposed revisions. We have carefully taken your valuable feedback into account and believe that our revised manuscript clearly highlights the theoretical novelty, mathematical rigor, and practical significance of our contributions. We hope that these detailed clarifications adequately address your initial concerns, and would be sincerely grateful if you could reconsider raising your evaluation accordingly. Thank you again for your thoughtful and constructive feedback

We sincerely thank the reviewer for the excellent questions and constructive suggestions, which have greatly helped us sharpen the manuscript. The detailed responses are provided below.

---

Q) **Points 1, 4, 8, and 9**\

A) We sincerely thank the reviewer for their detailed and constructive feedback, which has substantially aided us in improving both the clarity and scholarly rigor of our manuscript. In response to the suggestions, we will expand related work section to incorporate a critical discussion of relevant non-Euclidean Hamiltonian neural network literature, including key works such as Duong & Atanasov (2021), Li et al. (2022), along with subsequent developments, and provide a concise taxonomy (Euclidean vs. Riemannian vs. frame-bundle formulations) to clearly highlight our distinctive contributions. Additionally, we will move the architectural details of the frame-equivariant Transformer U-Net, currently described in the appendix, into the main manuscrip for enhanced clarity. We will also systematically audit and correct our use of citation commands throughout the manuscript, ensuring proper distinctions between parenthetical citations ($`citep`$) and textual citations ($`citet`$). Finally, we will replace vague references to the appendix with explicit pointers (e.g., "see Appendix D.2 for proof details") and include a comprehensive lookup table that clearly maps main-text propositions, algorithms, and figures to their corresponding appendix sections to facilitate easier navigation for readers.

Q2) (Geodesic Random Walk)

A2) Geodesic random walks (GRWs) implement tangent-space diffusion processes on Riemannian manifolds by applying exponential-map updates directly within the tangent bundle, followed by a retraction onto the manifold. While this approach provides computational simplicity and effectiveness for basic stochastic modeling, it introduces fundamental shortcomings when rigorously addressing the requirements of symplectic dynamics.

First, as a direct exponential-map increment, does not incorporate any structural mechanism to ensure the preservation of the canonical symplectic two-form $\omega$. This absence of a symplectic correction implies an intrinsic incompatibility between GRWs and symplectic geometry, which poses significant challenges in accurately representing Hamiltonian systems on curved manifolds. Second, although classical functional central limit theorems guarantee the convergence of GRWs to Brownian motion, the introduction of a Stratonovich lift which is essential for properly encoding Hamiltonian stochastic flows introduces additional drift terms driven by manifold curvature and torsion. The analytic characterization of this limiting behavior, specifically its weak convergence and preservation of the Poisson bracket structure fundamental to Hamiltonian dynamics, remains unproven in the existing literature. Next, GRW updates inherently depend on arbitrary choices of local coordinates and orthonormal frames. This lack of equivariance introduces potential biases of the stochastic processes modeled on manifolds.

Our proposed method addresses these critical theoretical and practical deficiencies by rigorously enforcing horizontality through a carefully constructed horizontal connection. This ensures that each update step preserves the canonical symplectic two-form $\omega$. By employing a construction that injects stochastic perturbations through frames that are equivariant with respect to the full $O(d)$ action, we achieve intrinsic gauge invariance, thereby removing arbitrary coordinate dependence and ensuring consistency and interpretability.

---

Q3) Implmentation of $\pi$

A3) In our implementation, each state in the frame-bundle is represented as a pair $U = (x, e) \in FM$, where $x \in M \subset \mathbb{R}^{d}$ denotes the embedded base manifold coordinates, and $e \in Mat_{d\times d}(\mathbb{R})$ is an orthonormal frame satisfying $e^{\top}g(x)e = I_d$. The projection $\pi: FM \to M$ is a straightforward tensor slice operation: $`return U[..., :d]`$, directly extracting the embedded coordinates and fully differentiable via PyTorch autograd. The inverse map $\pi^{-1}: M \to FM$ constructs a canonical orthonormal frame $e_{\mathrm{can}}(x)$ via Gram–Schmidt orthonormalization (implemented by Cholesky-based QR decomposition) applied to the metric tensor $g_{ij}(x)$ computed by autograd. To ensure gauge variability, we optionally multiply by a random orthogonal matrix $h \sim \mathrm{Haar}(O(d))$, forming $e(x,h)=e_{\mathrm{can}}(x)\,h$. The final tensor is $`torch.cat([x, e.reshape(-1)])`$. The horizontal lift operator $\mathrm{d}\pi^{-1}$, for a tangent vector $v \in T_xM$, is explicitly computed as: $H(U)[v] = [e v, -\Gamma(x)[ev]]^{\top}$ with the Christoffel symbols given by $\Gamma^k_{ij}(x) = \frac{1}{2} g^{k\ell}\left(\partial_i g_{\ell j} + \partial_j g_{i\ell} - \partial_\ell g_{ij}\right).$ This implementation is detailed in \textsc{Algorithm 2}, fully differentiable, analytically explicit for standard manifolds, and generalizable via autograd. In summary, projection $\pi$ is a tensor slice, and inverse projection $\pi^{-1}$ constructs orthonormal frames via QR decomposition where horizontal lift $\mathrm{d}\pi^{-1}$ is explicitly computed using a closed-form operator, ensuring differentiability and applicability without additional problem-specific hyperparameters.

---

Q5) What is the computational cost of the approach, both at training and testing time?

A5) To address the reviewer's query regarding computational efficiency, we provide detailed measurements of GPU memory usage, average runtime per training epoch, and inference step duration for our proposed method and baseline methods in the table below.

---

Q6) Explanation of Figure 1

A6) The provided figure depicts a conceptual illustration of a horizontally lifted stochastic Hamiltonian diffusion on the orthonormal frame bundle of a manifold. Let’s address the reviewer’s specific questions individually for improved clarity:

1. Does the manifold represent a torus? Yes, the manifold shown here conceptually represents a torus-like geometry. The toroidal shape visually symbolizes a nontrivial, curved configuration space $M = T^*Q$, where $Q$ typically denotes the configuration manifold, and $T^*Q$ is its cotangent bundle. Although the exact geometry of $Q$ might differ depending on the specific application or dataset, a toroidal shape effectively conveys a compact, periodic, and curved structure common in molecular dynamics or spin systems.

2. Is the horizontal space a flat space? Yes, within this illustration, the horizontal space $H_U$ at each frame $U$ (visualized by the tangent plane in the figure) represents a flat, Euclidean-like approximation to the local geometry of the frame bundle. Conceptually, the horizontal space is constructed by decomposing tangent vectors at a frame into horizontal and vertical components using a principal connection. The horizontal component (depicted as a flat plane tangent to the frame bundle) ensures compatibility with the manifold's geometry by enforcing that stochastic dynamics remain consistent with the intrinsic structure. This does not imply global flatness; rather, each horizontal space provides a local linear approximation at each point on the manifold.

3. What are the grey squares representing? The grey squares in the figure visually represent the tangent spaces at different points of the frame bundle. Each tangent space is attached to a specific frame, representing a linear (flat) approximation of the manifold at that particular point. These tangent planes symbolize the local horizontal spaces on which stochastic updates (geodesic random walks or Hamiltonian diffusions) are conducted before projecting back to the manifold. Specifically, the stochastic trajectories (depicted as paths $X_t$ and $U_t$) are initially computed in this tangent space (the flat, grey squares) and then projected back onto the curved base manifold through the exponential map or similar retraction operation.

Q7) How do predicted trajectories look like... A7) We regret that the rebuttal timeline does not allow us to include new visualizations. However, we can report a concise quantitative surrogate that faithfully reflects trajectory quality. As summarised below, our model attains accurate positional predictions and maintains physical fidelity, whereas a representative baseline such as GeoTDM reproduces positions but fails to conserve energy:

In other words, GeoTDM’s purely trajectory-matching loss indeed fits the positional path closely, but the resulting dynamics violate Hamiltonian conservation laws by an order of magnitude. By contrast, our symplectic, frame--consistent formulation yields trajectories that not only track the observed positions but also preserve energy up to ${\sim}3.5\%$ relative error, indicating faithful recovery of the underlying physics even in the absence of additional visual plots.

We sincerely thank the reviewer for their insightful and detailed feedback, which greatly helps us to clarify important aspects of our work. Apart from the major points explicitly addressed below, we will carefully incorporate all minor corrections and suggestions indicated by the reviewer into the revised manuscript.

---

Q1-1) First, how are the time derivatives...

In standard molecular dynamics (MD) simulations, rich trajectory datasets typically include explicit velocity and momentum data obtained directly from the simulation outputs. In our case, specifically for the synthetic data used in the experimental section, the full details of velocity and momentum values at each timestep are explicitly provided. Formally, these values are computed directly via the given Hamiltonian by evaluating its partial derivatives with respect to the momentum and position coordinates at each point in time, as follows: $\dot{q}_t = \frac{\partial H}{\partial p}(q_t, p_t), \dot{p}_t = -\frac{\partial H}{\partial q}(q_t, p_t)$ for known Hamiltonian function.

Furthermore, the core objective of our study is precisely to leverage this physically consistent, rich dynamical information for improved learning. Previous methodologies, specifically those categorized under "brute-force non-physical consistent" methods in our Table 1, typically do not utilize or explicitly incorporate this physically consistent dynamical information. Consequently, these methods exhibit inferior performance compared to our physically grounded approach. Please see Table 1 for a detailed performance comparison. Please see the following table:

---

Q1-2) It shoud depend on it as seen in Eq (4)-(5), ...

A) Indeed, as noted by the reviewer, the simulated state $U_t$ explicitly appears only in the second term of Eq.~(7), the trajectory matching objective, where the projected position $q_t$ obtained via the manifold projection of $U_t$ (using Eq. (4)-(5) and Algorithm 2) is directly compared with the true trajectory data. Specifically, $U_t$ is computed through the neural network as shown explicitly in line 11 of Algorithm 2, which ensures the gradient flows back through all intermediate computations. Therefore, $U_t$ is indeed optimized in an end-to-end manner. As suggested, we will move the table and additional details on how $U_t$ propagates through the computational pipeline, currently in the appendix, into the main manuscript to clarify this point explicitly.

---

Q2) I am wondering if the present approach...

We appreciate the reviewer's insightful question. As suggested, if explicit physical inductive biases or prior knowledge about the true Hamiltonian structure were already available such as separability between kinetic and potential terms, then one could certainly exploit this information to achieve more precise balancing between kinetic and potential energies. However, the core assumption of our current study is solely focused on the geometric structure of the given data manifold. Beyond this geometric structure, no additional system-specific prior information is assumed. Specifically, as formulated in Eq. (8), our model seeks to address the physical inverse problem by parameterizing three commonly used and interpretable terms including kinetic energy, single-particle potential, and pairwise interactions, aiming to recover the underlying physical structure directly from data. Even with this fully data-driven approach, our method remarkably preserves energy conservation properties of the original true Hamiltonian, as demonstrated empirically in Figure 4. This indicates that our learned model indeed effectively approximates the true underlying Hamiltonian structure and system's physical structure.

---

Q3) Top of page vi explains briefly how to enforce frame consistency, but in my view this is a particularly important contribution of this paper and it should be clearly explained in the main text how this is done.

A3) We sincerely thank the reviewer for highlighting this important point. We will explicitly incorporate a clear and detailed explanation of our approach to enforcing frame consistency into the main manuscript.

---

Q4) On the Proposition 3.1 A4) Thank you for this sharp observation. Proposition 3.1 is indeed formulated intentionally as a negative or impossibility result. Formally, it states:

\mathbb{P}\Bigl\[ \cdots \leq\delta\Bigr] \lesssim \exp\bigl(-\Omega,\delta^{1/2}|\Gamma|\_\infty^{3/2}R^{1/2}(\log R)^{1/4}\bigr). Equivalently, we have: \mathbb{P}\Bigl\[\sup\_{t,\theta}\mathcal{W}\bigl(P\_t(\theta),P\_{t,\mathrm{data}}\bigr)>\delta\Bigr] \geq1-\exp(\cdots)

which explicitly indicates that with overwhelming probability, at least some time or parameter choice will deviate by more than $\delta$. In other words, achieving uniformly small deviations across all times and parameter choices is exponentially unlikely.

This contrasts classical PAC-type generalization bounds, which upper-bound large-deviation probabilities. Here, our objective is to rule out the possibility of uniformly good approximation when either the curvature $\|\Gamma\|_\infty$ or the network capacity $R$ grows large. Condition (C1) already assumes an optimistic scenario ($\theta^*$ perfectly fits all training trajectories), yet the proposition demonstrates that even under this ideal scenario, the learned model cannot remain consistently close to the data distribution across all times and parameters. In our revision, we will explicitly restate Proposition 3.1 in this complementary form, clearly label it as an impossibility or general-lower-bound result.

---

Q5) On the Proposition 3.2.

A5) Yes, the rotation $h$ in Proposition 3.2 can indeed be considered time-dependent. To clarify further, the gauge transformation defined by the frame rotation $h \in O(d)$ represents the arbitrary choice of a local orthonormal frame at each point along the manifold. Such a choice can naturally vary with respect to time $t$, meaning that the frame bundle allows time-dependent gauge transformations $h(t)$. Formally, this flexibility is already accommodated in the general construction of horizontally lifted Hamiltonian diffusion (Algorithm 2). Thus, allowing $h$ to vary with time does not violate any assumption or theoretical result in our manuscript. Indeed, considering a time-dependent $h(t)$ aligns with the conceptual and practical motivations for gauge equivariance in the orthonormal frame bundle, enabling consistent and coherent representations of dynamics across temporally evolving frames. We appreciate this important clarification from the reviewer and will explicitly note in the revised manuscript that the gauge transformations $h$ can indeed be time-dependent.

---

Q6) how is stochasticity...

A6) We appreciate the reviewer's insightful question regarding the role of stochasticity in our experimental datasets and the behaviour of our method in the deterministic limit. First, the synthetic $n$-body gravitational dataset contains only extrinsic randomness: each trajectory is generated by sampling initial positions and momenta from an isotropic Gaussian distribution confined to a fixed energy shell, after which the trajectories evolve under the exactly deterministic Newtonian Hamiltonian. No noise is injected during integration. By contrast, the toroidal molecular-dynamics dataset embodies intrinsic randomness. The trajectories are produced with an OpenMM Langevin thermostat at $T = 300~\mathrm{K}$, so the stochastic Stratonovich term that models solvent–particle coupling is retained without alteration. Consequently, the data naturally reflect thermal fluctuations that drive conformational transitions. Regarding recovery of the deterministic case, our model is a neural Stratonovich SDE of the form $dZ_t = X_H(Z_t)\, dt + \sqrt{\beta_t} \Sigma(Z_t) \circ dW_t$. Setting $\beta_t \equiv 0$ eliminates the diffusion term and collapses the system to the ordinary Hamiltonian ODE $\dot{Z}_t = X_H(Z_t)$. All network components, loss functions, and optimization steps remain well-defined under this limit, so the architecture reduces exactly to a classical neural ODE or Hamiltonian neural network without any modification. In the revised appendix we provide an ablation in which we retrain on the deterministic $n$-body dataset with $\beta_t = 0$; the resulting test error coincides with that of a standard HNN, confirming that our framework incurs no loss of fidelity in the absence of stochastic forcing.

Dear Reviewer TByV

We sincerely thank the reviewer for the thoughtful and constructive comments, which significantly help improve the quality and clarity of our manuscript. We fully acknowledge and address your concerns as follows.

---

Q) ($**Physical Consistency**$) We would like to explicitly highlight a fundamental limitation in the trajectory matching methods included as baselines in our experimental comparisons. Although we intentionally adopted an evaluation setting favorable to these trajectory matching methods, to ensure fairness and transparency, it is crucial to emphasize that these methods generally fail to preserve essential physical properties such as energy conservation and symplecticity. $**In other words, their methods have fundamental flaws in real-world physical simulation scenarios.**$

To be more explicit, trajectory matching methods such as GeoTDM and EqMotion inherently fail to maintain fundamental physical consistency, particularly regarding crucial properties such as energy conservation. As clearly demonstrated by the quantitative results summarized in the table below, these baseline trajectory matching approaches suffer from severe violations of fundamental physical laws, reflecting intrinsic methodological shortcomings that significantly limit their reliability in handling physical data.

---

Q) ($**Simulation Complexity**$) This quadratic complexity may become prohibitive for long trajectories...

A) We appreciate the reviewer for raising the important point regarding computational complexity. However, we respectfully clarify that our proposed Gauge Equivariant Transformer U-Net does not process the entire trajectory sequence simultaneously, thus it does not exhibit quadratic $O(N^2)$ complexity. Rather, our model processes each individual timestep independently, resulting in an overall linear complexity $O(N)$. We believe the misunderstanding arose because Transformer architectures, in general, are commonly known for sequential attention mechanisms that yield quadratic complexity in the trajectory length, we fully understand the reviewer’s concern.

However, we respectfully clarify that our transformer operates as a set-based encoder-decoder model (inspired by Set Transformer) at each timestep individually, resulting in a complexity that is linear complexity with respect to number of particles rather than quadratic temporal complexity. The self-attention mechanism is explicitly employed to encode gauge-equivariant relationships among particles at each timestep independently, not temporal interactions across timesteps. The Transformer’s inherent self-attention mechanism naturally encodes symmetry, particularly gauge and permutation invariances, crucial for physically consistent modeling of manifold-based Hamiltonian systems. Such symmetry handling is essential to preserve geometric consistency without additional architectural complexity. In our preliminary experiments, we indeed observed that a gauge-equivariant MLP network, which lacks inherent permutation-equivariance, suffered from a noticeable accuracy drop (up to approximately -10%) compared to our Transformer-based approach.

---

Q) ($**Hybrid Loss**$) The training objective in Equation (7) does not appear to be simulation-free, as the second term ...

A) We acknowledge that this term indeed involves to acquire model trajectory, and thus is not strictly simulation-free. However, we clarify that nearly all existing Hamiltonian Neural Networks (HNNs) in this paper are typically trained using a proposed hybrid loss, consisting of trajectory-level reconstruction (simulation-dependent) and force-field matching terms.

The reason we chose this hybrid training objective in the current manuscript is to simultaneously achieve accurate trajectory matching and strict preservation of physical consistency, both of which are essential for physically-informed Hamiltonian models. As the reviewer correctly pointed out, virtually all physically-informed neural networks in our community proposed thus far have been designed for moderately-sized datasets (MD17) rather than super-large-scale datasets (ImageNet). Thus, incorporating a trajectory-level matching loss remains computationally practical and highly beneficial in such contexts. Furthermore, from a fair comparison perspective, methods in the second group (trajectory matching baselines) inherently optimize a training objective directly aligned with the test-time evaluation metrics. In contrast, conventional Hamiltonian Neural Network (HNN)-type methods do not directly optimize trajectory-level matching, and thus adopting this hybrid objective is essential while it fairly compensate and robustly evaluate the true predictive accuracy of HNN-based methods.

We fully agree that for future research aiming to further reduce training cost especially when dealing with significantly large datasets, it is indeed feasible to simplify the training by completely removing the trajectory matching term and relying exclusively on the force-field matching loss (thus making the training procedure fully simulation-free). Under such a simplified training scenario, we expect a notable improvement in computational efficiency; however, predictive accuracy may slightly degrade due to the absence of explicit trajectory-level guidance. We conducted additional experiments under this simulation-free setup and confirmed that the prediction error slightly increases in AD-3 dataset (but ours still outperforms the second-best), reflecting a modest trade-off between computational efficiency and trajectory accuracy.

Still, our simulation-free approach significantly outperforms trajectory-matching methods such as GeoTDM. However, we emphasize that this simplified objective (force-field only) may not be optimal if the primary goal is to minimize trajectory-based evaluation metrics (ADE/FDE).

---

Q) ($**Comparison with New Baselines**$) While the authors use GeoTDM as a strong baseline, several other relevant works on trajectory generation—such as [1, 2, 3]— ...

A) While we initially chose GeoTDM as a representative strongest baseline due to its state-of-the-art performance, we completely agree that referencing and discussing these additional works would enhance the context and completeness of our evaluation. We note that [1, 2] do not currently provide publicly available code implementations, making additional quantitative comparison challenging. In the new experiment during rebuttal, we checked that reference [3] reports performance metrics on the AD-3 and Spin benchmarks. We will explicitly include these comparative details and further clarify the performance differences in our revised manuscript.

---

Q) However, the proposed method requires access to the complete trajectory ${(q_t, p_t, \dot{q}_t, \dot{p}_t)}$ ...

A) We appreciate the reviewer’s comment regarding the requirement for accessing the complete trajectory ${(q_t, p_t, \dot{q}_t, \dot{p}_t)}$. We respectfully emphasize that such data is generally easy and straightforward to obtain in practice. For instance, molecular dynamics (MD) simulations naturally provide this complete set of dynamical variables. Indeed, comprehensive datasets containing $\textcolor{red}{highly ~ complex ~ large ~ dynamical ~ information}$ can be accessed by using free softwares (e.g, AMBER, LAMMPS) for various materials and scales. Rather, we would like to point out that traditional trajectory matching methods, despite the wide availability of complete dynamical datasets, have overlooked this rich physical information and remain limited by an "position-only-matching"-based approach (the majority of them only use position + connectivity graph information). Consequently, they frequently exhibit poor physical fidelity due to insufficient exploitation of the full dynamical structure inherent in these datasets.

Our future vision focuses on continuously developing current methodology that simultaneously outperform existing position-based approaches in predictive accuracy while strictly ensuring physical fidelity. Rather than prioritizing mere dataset scalability, we aim to incorporate rigorous verification procedures from domain experts such as strict physical conservation laws directly into our modeling frameworks, thereby ensuring models are robust, interpretable, and scientifically trustworthy.

Dear Reviewer TKAP,

Thank you very much for your constructive and detailed feedback, which has significantly helped us improve the manuscript. We carefully acknowledge your comments and commit to fully addressing each of your points. In the revised manuscript, we will explicitly emphasize that trajectory-matching methods fundamentally lack consideration for physical consistency, inherently limiting their ability to effectively handle physical data and resulting in substantially degraded performance.

We will explicitly integrate each of these clarifications and experimental details into a clearly structured appendix section, ensuring readability and comprehensibility for broader ML audiences without sacrificing theoretical rigor. We genuinely believe these concrete steps will fully address your constructive suggestions and significantly enhance the manuscript’s clarity and completeness. We humbly ask you to positively reconsider $**your evaluation**$ in light of these promised revisions.

---

Q) ($**General Audicence**.$) This work pays a lot of attention to theory, moving algorithm details..., One cannot assume ML audiences..

A) We sincerely appreciate your valuable suggestion regarding the readability and accessibility of our manuscript for a general audience. We fully acknowledge that the current manuscript, being heavily theoretical, might pose challenges to readers less familiar with geometric machine learning concepts.

To effectively address your constructive feedback, we will explicitly make the following improvements in the revised manuscript:

1. Providing a Clear and Concise Background Section} (currently Appendix A.1): The background information, currently detailed in Appendix A.1, will be significantly condensed and rewritten in simplified, accessible language and integrated directly into the main manuscript. We will include intuitive explanations, concise definitions, and simple illustrative examples to ensure general ML audiences easily grasp the fundamental concepts such as Hamiltonian dynamics, manifold geometry, and horizontal lifts.

2. Enhancing Writing Clarity and Readability: In addition to the above structural changes, we will systematically revise our manuscript to minimize unnecessary technical jargon, and clearly define essential technical terms at their first use in the main manuscript. We will also provide intuitive flowcharts to aid reader understanding whenever possible, thereby significantly improving accessibility to broader ML and interdisciplinary scientific audiences.

Once again, thank you very much for your thoughtful and valuable feedback that has significantly improved the manuscript's quality and accessibility.

---

Q) ($**Experimental Setups**.$) The author provides experimental setup for the proposed method in appendix. However, ...

A) We fully acknowledge your valid concern that our strong theoretical focus might have limited the accessibility and clarity of experimental details for broader ML audiences. As you rightly suggested, we will explicitly add a dedicated and detailed section in the Appendix of our revised manuscript, clearly summarizing comprehensive experimental settings and comparisons between our proposed method and all baseline methods. Specifically, we plan to carefully address your following points step-by-step:

1. Model Architectures and Hyperparameters (currently partial Appendix B): Provide explicit details for $**both our method and the baselines**$ such as proposed architecture, hidden dimension size, number of layers, attention heads, number of parameters, and precise optimizer and learning rate settings. Similarly, explicitly document hyperparameters for all baseline models to ensure fair, transparent comparisons.

2. Computational Complexity and Resource Requirements: Clearly document the practical GPU memory requirements (training and inference) and runtime measurements for both our method and baseline methods. We will explicitly include the detailed computational complexity analysis you requested, in the form of a concise table clearly summarizing GPU memory consumption and runtime (both training epoch and inference step).

3. Detailed Comparison to Prior Work: Provide an explicit and detailed discussion of the differences between our proposed approach and other relevant trajectory-generation baselines suggested (e.g., GeoTDM, EqMotion).

---

Q) "Do other baselines use the same input?"

A) We sincerely thank the reviewer for raising this insightful point, as it provides us the opportunity to clarify key distinctions between the baseline methods considered in our manuscript.

We have two primary groups of baseline methods:

1. $**Methodology 1 (variant of HNNs)**$ : These methods share the same underlying data input format (q, p) and training objective function as our approach $\textcolor{blue}{without~ connectivity~graphs}$. Thus, they inherently leverage complete physical-state information. However, a fundamental distinction arises from the type of force field employed: traditional HNN variants assume a purely Euclidean force field, whereas our method generalizes this by introducing a Riemannian force field on various differentiable manifolds. As similar to our model, all the methods in this category produce (q,p) by processing the same size of input (q,p)

2. $**Methodology 2 (Trajectory-Matching)**$ : These methods paradoxically utilize overly complex neural architectures, relying heavily on $\textcolor{red}{graph~inputs}$ even when processing inherently physical datasets. They thus fail to explicitly leverage physical information (positions q with momenta p), leading to poor physical consistency (Energy error). Our newly conducted experiments below demonstrate that these trajectory-matching baselines experience substantial performance degradation when explicitly provided with complete physical data and further degrade severely without graph priors. In contrast, our proposed method explicitly focuses solely on intrinsic physical structure, without relying on any additional graph information, thereby achieving both superior accuracy and rigorous physical consistency.

---

Q) "Do other baselines use the same size of model?"

A) We thank the reviewer for raising this practical consideration. Since the baseline models we compare against differ significantly in their architectures and structures, directly comparing the number of parameters would not yield a fair or meaningful evaluation. Above, we provide a practical comparison in terms of GPU memory usage and runtime performance during training and inference, which directly reflect real-world computational costs.

Thanks for your reply and additional results.

These experimental results and setups do not weaken the theorical results. Instead, they support the theory, as well as the "assumptions" in theory.

I raise my score to 5.