With energy harvesting capability,the Internet of things(IoT)devices transmit data de-pending on their available energy,which leads to a more complicated coupling and brings new techni-cal challenges to delay optimization.In this paper,we study the delay-optimal random access(RA)in large-scale energy harvesting IoT networks.We model a two-dimensional Markov decision process(MDP)to address the coupling between the data and energy queues,and adopt the mean field game(MFG)the-ory to reveal the coupling among the devices by uti-lizing the large-scale property.Specifically,to obtain the optimal access strategy for each device,we derive the Hamilton-Jacobi-Bellman(HJB)equation which requires the statistical information of other devices.Moreover,to model the evolution of the states dis-tribution in the system,we derive the Fokker-Planck-Kolmogorov(FPK)equation based on the access strat-egy of devices.By solving the two coupled equations,we obtain the delay-optimal random access solution in an iterative manner with Lax-Friedrichs method.Fi-nally,the simulation results show that the proposed scheme achieves significant performance gain com-pared with the conventional schemes.