Motivated by recent developments in perturbative calculations of the nonlinear evolution of large-scale structure, we present an iterative algorithm to reconstruct the initial conditions in a given volume starting from the dark matter distribution in real space. In our algorithm, objects are first moved back iteratively along estimated potential gradients, with a progressively reduced smoothing scale, until a nearly uniform catalog is obtained. The linear initial density is then estimated as the divergence of the cumulative displacement, with an optional second-order correction. This algorithm should undo nonlinear effects up to one-loop order, including the higher-order infrared resummation piece. We test the method using dark matter simulations in real space. At redshift $z=0$, we find that after eight iterations the reconstructed density is more than $95\%$ correlated with the initial density at $k\le 0.35\; h\mathrm{Mpc}^{-1}$. The reconstruction also reduces the power in the difference between reconstructed and initial fields by more than 2 orders of magnitude at $k\le 0.2\; h\mathrm{Mpc}^{-1}$, and it extends the range of scales where the full broadband shape of the power spectrum matches linear theory by a factor of 2-3. As a specific application, we consider measurements of the baryonic acoustic oscillation (BAO) scale that can be improved by reducing the degradation effects of large-scale flows. In our idealized dark matter simulations, the method improves the BAO signal-to-noise ratio by a factor of 2.7 at $z=0$ and by a factor of 2.5 at $z=0.6$, improving standard BAO reconstruction by $70\%$ at $z=0$ and $30\%$ at $z=0.6$, and matching the optimal BAO signal and signal-to-noise ratio of the linear density in the same volume. For BAO, the iterative nature of the reconstruction is the most important aspect.

26 pages, 14 figures, published version. Code available at https://github.com/mschmittfull/iterrec