The combination between non-orthogonal multiple access (NOMA) and hybrid automatic repeat request (HARQ) is capable of realizing ultra-reliability, high throughput and many concurrent connections particularly for emerging communication systems. This paper focuses on characterizing the asymptotic scaling law of the outage probability of HARQ-aided NOMA systems with respect to the transmit power, i.e., diversity order. The analysis of diversity order is carried out for three basic types of HARQ-aided downlink NOMA systems, including Type I HARQ, HARQ with chase combining (HARQ-CC) and HARQ with incremental redundancy (HARQ-IR). The diversity orders of three HARQ-aided downlink NOMA systems are derived in closed-form, where an integration domain partition trick is developed to obtain the bounds of the outage probability specially for HARQ-CC and HARQ-IR-aided NOMA systems. The analytical results show that the diversity order is a decreasing step function of transmission rate, and full time diversity can only be achieved under a sufficiently low transmission rate. It is also revealed that HARQ-IR-aided NOMA systems have the largest diversity order, followed by HARQ-CC-aided and then Type I HARQ-aided NOMA systems. Additionally, the users' diversity orders follow a descending order according to their respective average channel gains. Furthermore, we expand discussions on the cases of power-efficient transmissions and imperfect channel state information (CSI). Monte Carlo simulations finally confirm our analysis.