  We propose that loop integrals with internal heavy particles can be evaluated
by expanding in the limit of small external masses. This provides a
systematically improvable approximation to the integrals in the entire phase
space, and works particularly well for the high energy tails of kinematic
distributions (where the usual $1/M$ expansions cease to be valid). We
demonstrate our method using Higgs boson pair production as an example. We find
that at both one-loop and two-loop, our method provides good approximations to
the integrals appearing in the scattering amplitudes. Comparing to existing
expansion methods, our method are not restricted to a special phase space
region. Combining our efficient method to compute the two-loop amplitude with
an infrared subtraction method for the real emission corrections, we expect to
have a fast and reliable tool to calculate the differential cross sections for
Higgs boson pair production. This will be useful for phenomenological studies
and for the extraction of the Higgs self-coupling from future experimental
data. Our method can also be applied to other processes, such as the associated
production of the Higgs boson with a jet or a $Z$ boson.
