We apply the method of regions to the massive two-loop integrals appearing in the Higgs pair production cross section at the next-to-leading order, in the high energy limit. For the non-planar integrals, a subtle problem arises because of the indefinite sign of the second Symanzik polynomial. We solve this problem by performing an analytic continuation of the Mandelstam variables such that the second Symanzik polynomial has a definite sign. Furthermore, we formulate the procedure of applying the method of regions systematically. As a by-product of the analytic continuation of the Mandelstam variables, we obtain crossing relations between integrals in a simple and systematic way. In our formulation, a concept of “template integral” is introduced, which represents and controls the contribution of each region. All of the template integrals needed in the computation of the Higgs pair production at the next-to-leading order are given explicitly. We also develop techniques to solve Mellin-Barnes integrals, and show them in detail. Although most of the calculation is shown for the concrete example of the Higgs pair production process, the application to other similar processes is straightforward, and we anticipate that our method can be useful also for other cases
