The governing equations of multibody system (MBS) dynamics are differential algebraic equations (DAEs), and those equations must be solved by DAE integrators. In this work, the details of two families of DAE integrators, the backward difference formula (BDF) family and the generalized $\alpha$ method family, are amply described and discussed. At least two formulations are provided in the description of each family of integrators; the numerical schemes of each formulation for index-3, index-2, and index-1 DAEs are depicted in details; and the accuracy, efficiency, and robustness of all those integrators are discussed. Adaptive step size techniques are implemented in all the integrators, based on error estimations, and adaptive order technique is also implemented in each BDF integrator. The resultant spikes phenomena in calculations of the accelerations and Lagrange multipliers, which are due to index-3 formulations, are automatically resolved in the corresponding index-2 and index-1 schemes; while the index-3 integrators are usually more efficient than the other schemes in the simulations. Critical computational procedures, such as the initial condition analyses and reuse of Jacobian matrices, are described. Moreover, for BDF integrators, an error filtering method is introduced to improve the error estimations for velocity variables in the index-3 integrators, and several stabilization techniques are introduced to solve the issues caused by high order ($\geq3$) BDF schemes, which are not absolutely stable. Furthermore, explicit DAE integrators are briefly introduced, which do not require iterations and might be significant in some specific applications. Two benchmark examples are calculated using those integrators, and the pros and cons of each integrator are depicted and discussed. Typical integrator algorithms are provided in details in the appendix, which can be directly adopted to practical problems.