Reviewer's Comments:
Reviewer: 1

Comments to the Author
The authors present a differentiable symplectic N-body algorithm based on a previously existing algorithm. The authors' two major contributions are to improve the numerical stability of the previous integrator by analytically (as opposed to numerically, as in the previous algorithm) computing the Kepler-then-drift and drift-then-Kepler steps of their integrator, and to calculate the derivatives of the coordinates and velocities with respect to the initial conditions. These derivatives are needed for Hamiltonian MCMC methods (as well as by many nonlinear optimizers), which can be used to speed convergence for complex and/or high-dimensional posteriors. The authors make a Julia implementation of this code freely available, and compare their implementation with other existing open-source packages for N-body integration.

I am not an expert in this subfield (I have often used N-body integrators but have not worked on optimizing their performance), but despite that I found this paper readable and informative. I learned a lot from it, and recommend publication once my (minor) comments are addressed. 

Organizational comments:

- I believe the paper would benefit from moving some of the math to one or more appendices. In particular, much of the detail in the analytical calculation of the Kepler-then-drift and drift-then-Kepler corrections could be moved. I understand that the authors wish to present this derivation as a "main result," and so might be hesitant to make this change, but I think it is possible to make this part of the paper more understandable to novices by keeping the main points in the body of the paper and moving the details. The same applies to the calculation of derivates section.
- I recommend that the authors spend more time in the introduction motivating why calculating derivatives wrt the initial conditions is useful. They spend a lot of time discussing why a general N-body integrator is necessary, but less time motivating the need to compute derivatives. The authors mention HMC methods, but could also spend time discussing nonlinear optimizers that require gradients, as an example. As a side note, they might also mention that photodynamical codes enable computation of covariances between transit parameters and orbital parameters.
- The authors spend some time in section 2 reviewing symplectic integrators, but do not  complete the logic by showing how the Hamiltonian becomes the next position/momentum step. The DH paper has a good review of this if the authors just want to refer novice readers to this section. 
- I recently read an astronomy paper that had published "unit tests" for equations to check their accuracy (https://ui.adsabs.harvard.edu/abs/2021arXiv210306275L/abstract). I think a few tests like this would be very helpful for this paper. Of course, this is a new idea, so this is just a suggestion.

Minor comments:
- Page 3, first line of right column, "V is the total kinetic energy" -> "potential energy"
- Beginning of section 2: "the advantage of this approach" -> "the advantage of this approach relative to [WH method]"
- Section 2, second to last paragraph seems out of place here. Perhaps this could be moved to the end of the section?
- I spent some time puzzling over figure 1. Could the bottom-right-most line (beginning with -T_ij(deltax_kd# be repositioned so that it's clearer which part of the diagram that line is associated with? 

I have done some reading to try to make up for my lack of expertise in this sub-field, but I would not be offended at all if the authors and/or editor wish to request an additional referee for more detailed comments.
