Population III (Pop III) stars, forming as the first generation objects in the history of the Universe, bring an end to the dark age of the Universe and pollute the pristine intergalactic medium with heavy elements. Based on recent theoretical studies, there is general agreement that the typical mass scale of Pop III stars is larger than ordinary stars. However, although most massive stars are found in binary or multiple systems in local observations, we know little about the multiplicity and other related properties of Pop III stars. In this talk, I will present the result of our radiation hydrodynamics simulations in which we observe that Pop III stars generally form as massive multiples. For the first time, we self-consistently follow the multiple Pop III star formation under radiation feedback from cosmological initial conditions. For this purpose, we have developed a new code from an adaptive mesh refienement (AMR) self-gravitational magneto-hydrodynamics code, SFUMATO, by adding modules for ray tracing of radiation from Pop III proto-stars and for primordial-gas chemistry and thermodynamics. I will also discuss important physical processes in multiple Pop III star formation, as well as a possible link to observations such as gravitational wave events due to BH mergers.