  We study the formation of massive Population III binary stars using a newly
developed radiation hydrodynamics code with the adaptive mesh refinement and
adaptive ray-tracing methods. We follow the evolution of a typical primordial
star-forming cloud obtained from a cosmological hydrodynamics simulation.
Several protostars form as a result of disk fragmentation and grow in mass by
the gas accretion, which is finally quenched by the radiation feedback from the
protostars. Our code enables us, for the first time, to consider the feedback
by both the ionizing and dissociating radiation from the multiple protostars,
which is essential for self-consistently determining their final masses. At the
final step of the simulation, we observe a very wide ($\gtrsim
10^4\,\mathrm{au}$) binary stellar system consisting of $60$ and $70\,M_\odot$
stars. One of the member stars also has two smaller mass ($10\,M_\odot$)
companion stars orbiting at $200$ and $800\,\mathrm{au}$, making up a
mini-triplet system. Our results suggest that massive binary or multiple
systems are common among Population III stars.
