1. In the abstract the authors state "Our observations support the idea that superconductivity emerges out of spin gapped phases on ladders, driven by a spin-pairing mechanism". This observation is of interest with regards to high-temperature superconductivity in the cuprates, and is not so clearly stated in the conclusions. I think it would be good to state (and possibly discuss) this explicitly there, and it may also be worth noting that similar ideas lie at the heart of the Yang-Rice-Zhang ansatz for the single-electron propagator in the cuprates (see, e.g., the review Rep. Prog. Phys. 75 016502 (2012)) and recent works by Tsvelik on the cuprate problem (see, e.g., Phys. Rev. B 95, 201112 (2017)).
2. As mentioned by the authors in the introduction, the density matrix renormalization group (DMRG) has been used extensively to study $M$-leg ladders. How does the presented VQMC compare to DMRG for, e.g., the value of the ground state energy or for reproducing the phase diagram as a function of filling and U?
3. As shown explicitly in Ref.  by Lin, Balents and Fisher, many details of the $M>2$ phase diagram depend on the boundary conditions along the rungs. Following Eq. (12), the authors state the boundary conditions they use for $M=2$ (open), $M=4$ (antiperiodic) and $M=6$ (periodic), but its unclear why these choices were made. In light of Figs. 9 and 10 of Ref.  for $M=4$, this needs to be discussed and justified. Also, for $M=2$ aren't open and periodic boundary conditions equivalent?
4. It would also be good to mention somewhere about the statistical errors. If the error bars are smaller than the points, this should be stated explicitly.
5. More details on the extrapolation to the $M=\infty$ limit would be welcome - is the "2D limit" shown in Fig. 1 from an extrapolation or fully 2D numerical calculations?