Model-based control for robots has increasingly depended on optimization-based methods, such as differential dynamic programming (DDP) and iterative LQR (iLQR). These methods can form the basis of model-predictive control, which is commonly used for controlling legged robots. Computing the partial derivatives of the robot dynamics is often the most expensive part of these algorithms, regardless of whether analytical methods, finite difference, automatic differentiation (AD), or chain-rule accumulation is used. Since the second-order derivatives of the robot dynamics result in tensor computations, they are often ignored, leading to the use of iLQR, instead of the full second-order DDP method. In this article, we present analytical methods to compute the second-order derivatives of inverse and forward dynamics for open-chain rigid-body systems with multi-DoF joints and fixed/floating bases. An extensive comparison of accuracy and run-time performance with AD and other methods is provided, including the consideration of code-generation techniques in C/C++ to speed up the computations. For the 36 DoF ATLAS humanoid, the second-order inverse and forward dynamics derivatives take <inline-formula><tex-math notation="LaTeX">\approx 200 \,\mu \text{s}</tex-math></inline-formula>, and <inline-formula><tex-math notation="LaTeX">\approx \text{2.1}\,\text{ms}</tex-math></inline-formula>, respectively, on a 12th Gen Intel i5-12400 processor with 2.5 GHz clock-speed, resulting in a <inline-formula><tex-math notation="LaTeX">\approx 3.2 \times</tex-math></inline-formula> and <inline-formula><tex-math notation="LaTeX">\approx 3.8 \times</tex-math></inline-formula> speedup, respectively, over the AD approach.