In our paper we analyzed nonlocally coupled bosons in ladder lattices, or equivalently two-component bosons in onedimensional ͑1D͒ lattices. We showed that due to the direct Mott-insulator ͑MI͒ to pair-superfluid ͑PSF͒ transition, the MI lobes become severely distorted, leading to a counterintuitive behavior of the MI spatial regions in harmonic traps.In our paper, we characterized the possible phases, i.e., MI, PSF, and a phase of two superfluids ͑2SF͒, by employing as order parameters ͗b i,␣ ͘ and ͗b i,1 b i,2 ͘, where b i,␣ denotes the annihilation operator for the rung i and the wire ␣ =1,2. These order parameters characterize indeed the phases in the grand-canonical ensemble, since for PSF ͗b i,1 b i,2 ͘ 0 and ͗b i,␣ ͘ = 0, for 2SF both are nonzero, and for MI all are zero. Although our matrix-product-state ͑MPS͒ calculations conserve in principle the number of particles ͑and hence work in the canonical ensemble͒, a slight numerical inaccuracy in our simulations allowed for the violation of the number conservation. This effectively allowed for the identification of the MI, PSF, and 2SF phases. However, more accurate simulations do conserve the number of particles, and hence ͗b i,1 b i,2 ͘ = ͗b i,␣ ͘ = 0 in any case. Consequently, we have to rely in number conserving order parameters for the phases. In particular we employ the correlationsWe have calculated these correlation functions as a function of the site distance ⌬. For the MI phase both correlations decay exponentially with ⌬, for the PSF G 1 decays exponentially and G 2 polynomically, whereas for 2SF both decay polynomically. By fitting the correlation functions using exponential functions or power laws ͑see Fig. 1͒, we have identified for different values of the chemical potential and the tunneling rate J the different regions ͑see Fig. 2͒, which coincide almost exactly with those in our paper.