Differentiable physics engines: Warp and Newton

The idea: simulation you can backpropagate through

A classical simulator is a function you can only call. A differentiable simulator is a function you can also differentiate: for any scalar outcome L of a rollout, it returns ∂L/∂θ for the parameters θ that produced it.

Concretely, a simulation step advances state s and parameters θ: s(t+1) = f(s(t), u(t), θ), where u is the control input. Chain many steps into a rollout, compute a loss on the trajectory, and reverse-mode automatic differentiation propagates sensitivities back through every step, exactly the way training a neural network does, except the "network" here is the physics. Three problem classes fall out of this one capability:

• System identification. θ are physical parameters: masses, inertias, spring rates, damping coefficients, friction, tire model coefficients. Fit them to measured trajectories by gradient descent instead of derivative-free search. For a health-management engineer this has a second reading: a drifting θ estimated online is a degradation signal.
• Control and trajectory optimization. θ are the controls u(t) or policy weights. Gradients through the rollout give first-order policy improvement (the SHAC family of methods) and enable differentiable MPC, often with far fewer samples than zeroth-order RL.
• Design optimization. θ are geometric or structural parameters. The simulator becomes an adjoint solver you did not have to derive by hand, the same pattern as the DoMINO design-sensitivities workflow from Part I, but with the real physics instead of a surrogate.

The price of admission is twofold: memory (reverse mode stores or recomputes intermediate states across the rollout, hence gradient checkpointing) and gradient quality through discontinuous events like contact, which gets its own honest section later. Neither price is hidden in Warp or Newton; both are engineered around explicitly.

Feel the gradient: a live identification demo

Before any framework, the idea itself, running in this page. A quarter-car with hidden spring rate and damping drives over a bump; noisy displacement measurements are all you get. Press run and watch gradient descent recover the hidden parameters by backpropagating through the simulation.

Everything below is plain JavaScript in this page's source: a semi-implicit Euler forward pass, a hand-written adjoint (the "tape") of about twenty lines, and Adam on log-parameters. It is the same mathematics Warp generates automatically for every kernel, shrunk to a size you can read in one sitting.

Three things worth noticing as it runs. The stiffness converges before the damping, because spring force scales with displacement while damper force scales with velocity, which the noise punishes harder. Higher noise settings slow and roughen convergence but do not break it; the estimator is averaging over 1,500 timesteps of evidence. And the entire optimization is deterministic given the noise seed: gradients, not guesses.

NVIDIA Warp: the substrate

Warp is not a physics engine. It is the layer that makes writing one in Python viable: kernel functions JIT-compiled to native CPU and CUDA code, with adjoints generated automatically.

You decorate a Python function with @wp.kernel, and Warp compiles it through its own type-checked subset of Python into efficient parallel code, caching the result. The same source runs on CPU or any CUDA device. Because Warp generates both the forward kernel and its adjoint, any program composed of kernels is differentiable end to end via wp.Tape. The published industrial numbers are striking: NVIDIA reports speedups up to 669x over CPU baselines and large gains over JAX in specific workloads, with zero-copy interoperability into PyTorch and JAX tensors so the physics slots directly into ML training loops.

A sixty-second taste

import warp as wp

@wp.kernel
def quarter_car(x: wp.array(dtype=float), v: wp.array(dtype=float),
                f: wp.array(dtype=float), k: float, c: float):
    i = wp.tid()                          # one thread per environment
    f[i] = -k * x[i] - c * v[i]           # spring-damper force

# thousands of parallel suspension corners in one launch
wp.launch(quarter_car, dim=n_envs, inputs=[x, v, f, k, c])

And the gradient, for free

k = wp.array([28000.0], dtype=float, requires_grad=True)   # unknown spring rate

tape = wp.Tape()
with tape:
    for t in range(steps):
        wp.launch(step_dynamics, dim=n, inputs=[state, k, dt])
    wp.launch(trajectory_loss, dim=1, inputs=[state, measured, loss])

tape.backward(loss)
print(k.grad)   # dLoss/dk: descend on this to identify the spring rate

# Sketch of the idiom; consult the Warp docs for current signatures.

The module map

Autodiff under the hood

Warp's differentiability is not a wrapper around PyTorch; it is generated at kernel compile time. Knowing the mechanism tells you exactly what will and will not differentiate, and how to debug it when it goes quiet.

How the adjoint gets made

When Warp compiles a @wp.kernel, it emits two functions: the forward kernel you wrote, and an adjoint kernel derived from it by reversing the dataflow, each primitive operation replaced by its derivative rule. Arrays created with requires_grad=True carry a companion .grad buffer. A wp.Tape then simply records every launch inside its context; calling tape.backward(loss) replays the recorded launches in reverse order using the adjoint kernels, accumulating into the .grad buffers. Reverse-mode over the whole program falls out of reverse order over its launches.

The working rules

• Gradients accumulate. Adjoint buffers are added into, never overwritten, exactly like PyTorch. Between optimizer steps, call tape.zero() or you will be descending on the sum of every gradient you have ever computed. This is the single most common silent bug.
• The tape stores what the backward pass needs. For long rollouts that means intermediate state, and memory grows with horizon length. The remedies, in order: reuse preallocated state arrays where the math permits, checkpoint segments (store every k-th state, recompute the rest during backward), and truncate horizons where the problem allows it.
• Custom gradients are first-class. For operations whose true derivative is useless or undefined, a hard contact response, a lookup table, a clamp at its corner, Warp lets you register a replacement adjoint for a @wp.func. This is the sanctioned mechanism for the smoothed contact gradients discussed later: keep the hard physics forward, supply a soft derivative backward, and document that you did.
• Not everything differentiates. Integer indexing logic, branches on data, and atomic scatter patterns have the derivative semantics you would derive by hand, which is sometimes zero. If a parameter's gradient is exactly zero, suspect a branch or cast severing the path before suspecting the physics.

Debugging the invisible

A wrong gradient looks like a right gradient until it ruins an optimization, so verification is part of the workflow, not an afterthought. The tools, cheapest first: compare against central finite differences on a shrunken problem (the case study in Section 12 does exactly this and agrees to 1.4e-8); inspect the tape's recorded launch graph to confirm the operations you think are being differentiated actually were; enable Warp's floating-point verification mode during development to trap NaN and inf the moment they appear rather than three kernels later; and bisect long rollouts, checking gradient norms per segment, to localize where sensitivity explodes or dies. Gradient norms that grow exponentially with horizon length are not a bug in the tape; they are Lyapunov instability in your dynamics, and the fix is shorter horizons, not different code.

Newton: the engine

Newton is what happens when Warp's simulation module grows up, moves out, and acquires institutional parents: an open-source, GPU-accelerated, multi-solver physics engine created by NVIDIA, Google DeepMind, and Disney Research, governed by the Linux Foundation under Apache 2.0.

Newton explicitly supersedes warp.sim with a more general architecture. Its design commitments: end-to-end GPU execution with CUDA graphs so thousands of environments run without CPU bottlenecks, OpenUSD as the native scene and asset representation, differentiability through inputs, controls, and parameters, and a pluggable solver architecture rather than one monolithic integrator. It reached 1.0 in March 2026, integrated as a first-class backend of Isaac Lab 3.0.

What the solver menu buys you

Two adjacent pieces complete the picture. Hydroelastic contact and an SDF collision library address the contact-quality failure modes that have long plagued sim-to-real transfer. And NeRD, Neural Robot Dynamics, explores the hybrid frontier: learned dynamics models inside the engine that can be fine-tuned with real-world data, which is the engine-side mirror of everything Part I covered.

Anatomy of the Newton API

Like the recipe anatomy in Part I, this is the skeleton every Newton program shares. Learn the lifecycle once and each example reads as a variation.

# The lifecycle, end to end
ModelBuilder  →  finalize()  →  Model  →  State / Control  →  Solver.step  →  Viewer / USD
   build             │        GPU arrays     per-timestep        physics         output
 or import           │       (parameters)      buffers          advance         artifact

• ModelBuilder is the construction phase: procedural calls to add bodies, joints, and shapes, or importers that ingest whole MJCF, URDF, and USD descriptions. Nothing touches the GPU yet; you are describing, not simulating.
• finalize() is the phase change. The description is compiled into a Model: flat, GPU-resident arrays of masses, inertias, joint limits, shape parameters, and contact properties. This flatness is the point, and the opportunity: physical parameters live in plain arrays you can read, write, and, crucially, differentiate with respect to.
• State holds everything that changes per timestep (positions, velocities, in generalized or maximal coordinates depending on solver); Control holds actuation. You allocate them from the model and ping-pong two states through the loop.
• Solver.step advances physics. Solvers are interchangeable behind the same call because they consume the same Model/State contract, which is what makes "swap MuJoCo Warp for XPBD and compare" a one-line experiment.
• Viewers and USD export sit outside the loop: attach a window for interactive work, or write the rollout to a USD file for replay, review, and diffing.

The two idioms that unlock scale

Environment replication. To simulate 4,096 vehicles, you do not create 4,096 models; you replicate one articulation into a single Model whose arrays are 4,096 environments long, then step them in one launch. Batch dimension, not object count.

Parameter randomization by array writes. Because parameters are flat arrays indexed by environment, domain randomization and per-vehicle system identification are the same operation: write different values into each environment's slice of the mass, stiffness, or friction arrays. No scene rebuilding, no per-env objects, and the writes themselves can sit inside the differentiated program.

A short OpenUSD primer

Newton's choice of OpenUSD as its native scene representation looks like a rendering detail and is actually an engineering-data decision.

Universal Scene Description, from Pixar, is a scene graph with a superpower: composition. A USD stage is assembled from layers that override and extend one another non-destructively, so a robot asset, an environment, and a set of simulation-specific tweaks live in separate files composed at load time. Physics enters through the UsdPhysics schemas, which attach rigid-body, collider, joint, and material semantics to scene prims as typed attributes. The practical consequences for a simulation team:

• One source of truth. The same asset flows through CAD-adjacent tools, Omniverse, Isaac Sim, and Newton without lossy format hops; simulation overrides live in a thin layer on top rather than in a forked copy.
• Versionable and reviewable. The ASCII form (.usda) diffs in git like code. A changed inertia or joint limit shows up in review as a changed line, not a changed binary blob.
• Rollouts as artifacts. Newton can write simulations out as USD, which turns every run into a portable, replayable, inspectable file, the difference between "it looked wrong yesterday" and attaching the offending trajectory to the bug report.

Getting started

Both projects are a pip install away. Versions below verified on PyPI at the time of writing.

# Warp (v1.15.0 at time of writing)
pip install warp-lang

# Newton (v1.3.0; the old name newton-physics now redirects here)
pip install newton

# MuJoCo Warp, usable standalone or as Newton's rigid-body backend
pip install mujoco-warp

First contact with Newton

# Discover the shipped examples
python -m newton.examples --list

# Run one with the OpenGL viewer
python -m newton.examples basic_viewer --viewer gl --num-frames 500

# Load a robot description and pick your device
python -m newton.examples basic_urdf --device cuda:0

# Render to an OpenUSD file instead of a window
python -m newton.examples basic_viewer --viewer usd --output-path my_run.usd

The example runner accepts a shared set of flags: viewer backends (gl, usd, rtx, rerun, viser, null), device selection, and frame counts. The usd path matters more than it looks: writing simulations to OpenUSD makes every run a portable, inspectable asset rather than a fleeting window.

The programming idiom

import newton

# 1. Build: describe bodies, joints, shapes (or import MJCF/URDF/USD)
builder = newton.ModelBuilder()
builder.add_ground_plane()
# ... add articulation, e.g. a double-wishbone corner ...
model = builder.finalize()

# 2. Pick a solver and allocate states
solver = newton.solvers.SolverMuJoCo(model)
state_in, state_out = model.state(), model.state()

# 3. Step: the loop is yours, which is the whole point
for t in range(horizon):
    solver.step(state_in, state_out, control, contacts, dt)
    state_in, state_out = state_out, state_in

# API sketch to convey the shape; the library is young and moving,
# so treat the current docs as the source of truth for signatures.

Performance engineering

The throughput headlines are real, but they are earned by a handful of disciplines, and lost by their absence. Most "Warp is slow" reports decompose into one of these rows.

The differentiable engine landscape

Newton and Warp sit in a lively neighborhood. Choosing among the neighbors is mostly a question of ecosystem, autodiff framework, and physics coverage.

Vehicle applications gallery

Eight concrete places a vehicle engineer can point this stack, from published workflows to well-scoped internal projects.

A worked case study: half-car identification and damper-drift detection

The Section 02 demo, grown up: a pitch-plane half car, four unknown suspension parameters, realistic sensor noise, a hand-written adjoint verified against finite differences, and a deliberately degraded damper to detect. Every number below was produced by the companion script shipped alongside this page.

Protocol

The plant is a pitch-plane half car (m = 1400 kg, I = 2100 kg·m², axle distances a = 1.2 m, b = 1.5 m) with front and rear spring-damper corners. Ground truth: kf = 32,000 N/m, cf = 1,500 N·s/m, kr = 28,000 N/m, cr = 1,300 N·s/m. The excitation is a 5 cm half-sine bump traversed at 10 m/s, hitting the rear axle 0.27 s after the front. Sensors are heave and pitch displacement at 1 kHz over a 4 s horizon (4,000 steps, semi-implicit Euler at 1 ms). Identification is Adam on log-parameters, 500 iterations, from a deliberately poor initial guess (18,000 / 800 / 18,000 / 800). Gradients come from a hand-written discrete adjoint, the same reverse pass Warp generates automatically, small enough here to read in full.

# The adjoint core (backward over one step), verbatim from the script
for t in range(N - 1, -1, -1):
    gz  += 2.0 * rz[t] / (SZ * N)      # loss injects adjoint at t+1
    gth += 2.0 * rt[t] / (STH * N)
    gthd += DT * gth;  gthdd = DT * gthd
    gzd  += DT * gz;   gzdd  = DT * gzd
    gFf = gzdd / M - A_ * gthdd / I    # accelerations → forces
    gFr = gzdd / M + B_ * gthdd / I
    gkf += -zf * gFf;  gcf += -zfd * gFf     # parameter gradients
    gkr += -zr * gFr;  gcr += -zrd * gFr
    gzf = -kf * gFf;  gzfd = -cf * gFf       # forces → state at t
    gzr = -kr * gFr;  gzrd = -cr * gFr
    gz  += gzf + gzr;   gth  += -A_ * gzf + B_ * gzr
    gzd += gzfd + gzrd; gthd += -A_ * gzfd + B_ * gzrd

Result 1 · The adjoint is correct

Before trusting any optimization, the gradient itself was checked against central finite differences over all four parameters: maximum relative difference 1.4 × 10⁻⁸. This verification step costs seconds and should be a fixed ritual of any differentiable-simulation project; it is rung R4 of the ladder for a reason.

Result 2 · Identification accuracy vs sensor noise

Two patterns worth internalizing. Noise-free identification is exact to display precision, confirming the estimator and the adjoint jointly. And damping is consistently the hardest parameter: damper force scales with velocity, whose signal-to-noise is worst, so cf carries the largest error at every noise level, the same effect visible in the live demo. Cost: each 500-iteration run over 4,000-step rollouts took about 10 seconds single-threaded on CPU; the identical computation in Warp runs thousands of vehicles as a batch dimension.

Result 3 · Drift detection: the VHM payoff

A second dataset was generated from a vehicle whose front damper had degraded by 25% (cf: 1,500 → 1,125 N·s/m), with everything else unchanged, and re-identified starting from the healthy estimates. The fit localized the fault cleanly:

This is the argument of the whole guide in one table. A threshold on ride acceleration would eventually notice something; the identified parameter vector says which component, by how much, with the healthy parameters standing still as the control group. Run continuously against telemetry, the fitted cf is a physically interpretable degradation trend, prognosable like any other health indicator, and the residual after fitting bounds what the model cannot explain.

Gradients through contact: the fine print

Differentiable simulation's hardest honesty problem: contact is discontinuous, and gradients through discontinuities can be biased, exploding, or subtly wrong while looking perfectly plausible.

Choosing a gradient estimator

"Differentiable" names a capability, not a decision. For any given problem you still choose how to estimate the gradient, and the contact caveats of Section 13 are really a claim that this choice should be deliberate.

A serviceable decision procedure: start with the analytic adjoint and verify it against finite differences on a shrunken problem. If verification fails near contact events or gradients are noisy-looking despite deterministic code, switch the contact-adjacent portion to randomized smoothing or supply a custom smoothed adjoint via Warp's custom-gradient mechanism. If horizons are long enough that gradient norms grow exponentially, truncate and go hybrid. Reserve pure zeroth-order methods for the problems that have defeated everything above; they are the most robust and the most expensive.

Sim-to-real: closing the loop with the physical vehicle

A simulator earns its keep only if conclusions transfer. The gap between simulation and reality is not one thing, and each of its components has a different remedy.

A taxonomy of the gap

Five distinguishable error sources, in roughly increasing order of stubbornness: parameter error (masses, stiffnesses, friction coefficients are wrong), contact model error (the functional form of contact is wrong, not just its coefficients), actuator error (real actuators have torque curves, delays, saturation, backlash, and temperature dependence that idealized joints lack), sensor error (bias, noise, quantization, and latency between the world and the measurement), and unmodeled dynamics (flexibility, slop, and physics your model class cannot express at all).

Identify first, then randomize

Domain randomization, training against a distribution of simulated worlds so the real world looks like one more sample, is the standard robustness tool, but randomizing blindly trades performance for robustness you may not need. The disciplined sequence uses the machinery of this guide: identify what can be identified from real trajectories (Sections 02 and 12), which collapses parameter error to a posterior; then randomize around that posterior, wide where identification was uncertain, narrow where it was confident. In Newton this is mechanically trivial, because per-environment randomization is just writes into the model's parameter arrays; every environment in the batch is a different plausible vehicle.

Model the boring parts

Transfer failures are attributed to contact more often than they are caused by it. Actuator and sensor models are the unglamorous half of sim-to-real: a first-order lag and delay on actuation, torque-speed limits, IMU bias random walks, and measurement latency cost little to implement and routinely dominate the residual gap. They also belong inside the differentiated program, so their parameters are identified alongside the plant's.

Residuals for the last mile, and for health

What remains after identification and honest actuator modeling is structural: the model class itself falls short. This is where learned residual dynamics, the NeRD pattern, earns its place: keep the analytic core, let a network absorb the systematic remainder, fine-tune it on real trajectories. The VHM reading is the same as everywhere in this guide: a residual model whose magnitude grows over a vehicle's life is telling you the physical system has moved, and the direction of its growth says something about where.

A validation protocol worth writing down

Hold out maneuvers, not just trajectories: validate on excitation types the identification never saw. Report reality-gap metrics per channel (position, velocity, load) rather than one aggregate. Track identified parameters over repeated sessions; their scatter is your honest uncertainty. And predeclare the acceptance threshold before running the comparison, the same discipline a DFMEA review would demand of any other verification activity.

How this stacks with PhysicsNeMo

Part I and Part II are not rivals; they are two halves of one workflow. Surrogates learn from simulation; differentiable engines make simulation itself trainable infrastructure.

Four composition patterns cover most of the practical territory:

• Data generation. Newton and Warp-based solvers generate the training corpora that PhysicsNeMo recipes consume, at GPU throughput and with programmatic parameter sweeps, feeding Curator ETL exactly like any other solver output.
• Physics-in-the-loop training. Because Warp interoperates zero-copy with PyTorch, a differentiable solver can sit inside a training loop as a layer: the network proposes, the physics evaluates, and gradients flow through both. This is the solver-in-the-loop pattern that pure surrogate training cannot express.
• Residual physics, both directions. NeRD adds a learned residual to an analytic engine; PhysicsNeMo's predictor-corrector fine-tuning adds a learned corrector to a pretrained surrogate. Same idea, opposite anchors. A mature pipeline may use both, with the residual magnitudes doubling as model-trust diagnostics.
• Cross-validation. A differentiable analytic model is the cheapest independent check on a learned surrogate, and vice versa. Disagreement between them, localized in state space, is a map of where neither should be trusted alone, which is exactly the kind of evidence a V&V case for a digital twin needs.

One sentence to carry away: PhysicsNeMo makes physics learnable; Warp and Newton make it differentiable; the interesting engineering lives where the two meet.

A starter project ladder

Five rungs, each one afternoon to one week, each producing something you keep.

Glossary

The working vocabulary of this guide, in one place.

Adjoint method

Computing sensitivities of a scalar output by solving the reverse (adjoint) problem once, at cost independent of the number of parameters; reverse-mode autodiff is its algorithmic form.

Tape

The record of operations (in Warp, kernel launches) captured during a forward pass, replayed in reverse with adjoint kernels to accumulate gradients.

BPTT

Backpropagation through time: applying reverse-mode differentiation across the steps of a rollout. Memory grows with horizon; gradients can grow with the dynamics' Lyapunov exponents.

Reverse vs forward mode

Reverse mode: one backward pass yields gradients of one output w.r.t. all inputs (ideal for many parameters, one loss). Forward mode: one pass per input direction (ideal for few inputs, many outputs).

Gradient checkpointing

Storing only periodic states during the forward pass and recomputing the segments between them during backward, trading compute for tape memory.

Semi-implicit (symplectic) Euler

Integrator that updates velocity first, then position with the new velocity; the workhorse of real-time physics for its energy behavior at large timesteps.

Reduced vs maximal coordinates

Reduced: state is joint coordinates, constraints exact by construction (Featherstone). Maximal: each body carries full 6-DoF state, joints enforced as constraints (XPBD-style). Different trade-offs in speed, drift, and gradient conditioning.

XPBD

Extended position-based dynamics: compliant constraints solved on positions, stable and fast; the modern form of the position-based family.

VBD

Vertex block descent: a per-vertex optimization scheme for deformables (cloth, cables, volumetric solids) well suited to GPU parallelism.

MPM

Material point method: continuum material carried on particles, computed on a background grid; the tool of choice for granular media, snow, and soil.

SDF

Signed distance field: an implicit surface representation returning distance-to-surface (negative inside), enabling fast, smooth collision queries.

Hydroelastic contact

A compliant contact model producing distributed pressure fields over patch areas instead of point forces, better-behaved physically and for gradients than rigid point contact.

CUDA graph

A recorded sequence of GPU operations replayed as a unit, eliminating per-launch overhead; the key to stepping thousands of small physics kernels efficiently.

Zero-copy interop

Sharing GPU memory between frameworks (Warp ↔ PyTorch/JAX) via DLPack without copying, letting physics sit inside a training loop at no transfer cost.

MJCF / URDF / USD

The three asset dialects Newton ingests: MuJoCo's simulation-native XML, ROS's ubiquitous robot description, and OpenUSD's composable scene description.

System identification

Estimating a model's physical parameters from measured input-output data; with a differentiable simulator, gradient descent on simulation-vs-measurement error.

Domain randomization

Training against a distribution of simulated worlds so the real world resembles a sample from it; best centered on an identified posterior rather than guessed bounds.

Sim-to-real gap

The composite discrepancy between simulated and physical behavior: parameter, contact, actuator, sensor, and unmodeled-dynamics error, each with its own remedy.

Learning resources

References

1. Newton: an open-source, GPU-accelerated physics engine for robotics, built on NVIDIA Warp. Linux Foundation project, Apache 2.0. github.com/newton-physics/newton
2. Macklin, M. "Warp: A High-Performance Python Framework for GPU Simulation and Graphics." NVIDIA GTC, 2022. github.com/NVIDIA/warp
3. NVIDIA Technical Blog. "Build Accelerated, Differentiable Computational Physics Code for AI with NVIDIA Warp." 2026. developer.nvidia.com/blog
4. NVIDIA Developer. "Newton Physics Engine." Overview page with collaborator and design statements. developer.nvidia.com/newton-physics
5. Google DeepMind. MuJoCo Warp (MJWarp): GPU-optimized MuJoCo on NVIDIA Warp. github.com/google-deepmind/mujoco_warp
6. Isaac Lab: a GPU-accelerated simulation framework for multi-modal robot learning, with Newton engine integration. 2025. arXiv:2511.04831

7. Hu, Y., et al. "DiffTaichi: Differentiable Programming for Physical Simulation." ICLR 2020. arXiv:1910.00935
8. Xu, J., et al. "Accelerated Policy Learning with Parallel Differentiable Simulation" (SHAC). ICLR 2022. arXiv:2204.07137
9. Suh, H. J. T., Simchowitz, M., Zhang, K., Tedrake, R. "Do Differentiable Simulators Give Better Policy Gradients?" ICML 2022. arXiv:2202.00817
10. "Differentiable Physics Simulations with Contacts: Do They Have Correct Gradients w.r.t. Position, Velocity and Control?" 2022. arXiv:2207.05060
11. Macklin, M., Müller, M., Chentanez, N. "XPBD: Position-Based Simulation of Compliant Constrained Dynamics." Motion in Games, 2016.
12. Howell, T., Le Cleac'h, S., Kolter, Z., Schwager, M., Manchester, Z. "Dojo: A Differentiable Physics Engine for Robotics." 2022. arXiv:2203.00806
13. Werling, K., Omens, D., Lee, J., Exarchos, I., Liu, C. K. "Fast and Feature-Complete Differentiable Physics for Articulated Rigid Bodies with Contact" (Nimble). RSS 2021. arXiv:2103.16021

14. Todorov, E., Erez, T., Tassa, Y. "MuJoCo: A Physics Engine for Model-Based Control." IROS 2012. MJX (JAX) ships within the MuJoCo repository. github.com/google-deepmind/mujoco
15. Freeman, C. D., et al. "Brax: A Differentiable Physics Engine for Large Scale Rigid Body Simulation." 2021. arXiv:2106.13281
16. Genesis: a generative and universal physics platform for robotics and embodied AI. github.com/Genesis-Embodied-AI/Genesis
17. Ataei, M., Salehipour, H. "XLB: A Differentiable Massively Parallel Lattice Boltzmann Library in Python." Computer Physics Communications, 2024. github.com/Autodesk/XLB
18. Drake: model-based design and verification for robotics. Toyota Research Institute. drake.mit.edu
19. Project Chrono: multi-physics simulation with the Chrono::Vehicle module. projectchrono.org