  Based on the path integral formulation of the reduced density matrix, we
develop a scheme to overcome the exponential growth of computational complexity
in reliably extracting low-lying entanglement spectrum from quantum Monte Carlo
simulations. We test the method on the Heisenberg spin ladder with long
entangled boundary between two chains and the results support the Li and
Haldane's conjecture on entanglement spectrum of topological phase. We then
explain the conjecture via the wormhole effect in the path integral and show
that it can be further generalized for systems beyond gapped topological
phases. Our further simulation results on the bilayer antiferromagnetic
Heisenberg model with 2D entangled boundary across the $(2+1)$D O(3) quantum
phase transition clearly demonstrate the correctness of the wormhole picture.
Finally, we state that since the wormhole effect amplifies the bulk energy gap
by a factor of $\beta$, the relative strength of that with respect to the edge
energy gap will determine the behavior of low-lying entanglement spectrum of
the system.
