The orthogonal space random walk (OSRW) algorithm is generalized to be the orthogonal space tempering (OST) method via the introduction of the orthogonal space sampling temperature. Moreover, a double-integration recursion method is developed to enable practically efficient and robust OST free energy calculations, and the algorithm is augmented by a novel θ-dynamics approach to realize both the uniform sampling of order parameter spaces and rigorous end point constraints. In the present work, the double-integration OST method is employed to perform alchemical free energy simulations, specifically to calculate the free energy difference between benzyl phosphonate and difluorobenzyl phosphonate in aqueous solution, to estimate the solvation free energy of the octanol molecule, and to predict the nontrivial Barnase-Barstar binding affinity change induced by the Barnase N58A mutation. As demonstrated in these model studies, the DI-OST method can robustly enable practically efficient free energy predictions, particularly when strongly coupled slow environmental transitions are involved.