The algorithm based on integration over Lefschetz thimbles is a promising method to resolve the sign problem for complex actions. However, this algorithm often meets a difficulty in actual Monte Carlo calculations because the configuration space is not easily explored due to the infinitely high potential barriers between different thimbles. In this paper, we propose to use the flow time of the antiholomorphic gradient flow as an auxiliary variable for the highly multimodal distribution. To illustrate this, we implement the parallel tempering method by taking the flow time as a tempering parameter. In this algorithm, we can take the maximum flow time to be sufficiently large such that the sign problem disappears there, and two separate modes are connected through configurations at small flow times. To exemplify that this algorithm does work, we investigate the (0 + 1)-dimensional massive Thirring model at finite density and show that our algorithm correctly reproduces the analytic results for large flow times such as T = 2.