  We compute the two-loop helicity amplitudes for the production of three
photons at hadron colliders in QCD at leading-color. Using the two-loop
numerical unitarity method coupled with analytic reconstruction techniques, we
obtain the decomposition of the two-loop amplitudes in terms of master
integrals in analytic form. These expressions are valid to all orders in the
dimensional regulator. We use them to compute the two-loop finite remainders,
which are given in a form that can be efficiently evaluated across the whole
physical phase space. We further package these results in a public code which
assembles the helicity-summed squared two-loop remainders, whose numerical
stability across phase-space is demonstrated. This is the first time that a
five-point two-loop process is publicly available for immediate
phenomenological applications.
