Entanglement entropy dynamics of non-Gaussian states in free boson systems: Random sampling approach

We develop a random sampling method for calculating the time evolution of the R\'{e}nyi entanglement entropy after a quantum quench from an insulating state in free boson systems. Because of the non-Gaussian nature of the initial state, calculating the R\'{e}nyi entanglement entropy calls for the exponential cost of computing a matrix permanent. We numerically demonstrate that a simple random sampling method reduces the computational cost of a permanent; for an $N_{\mathrm{s}}\times N_{\mathrm{s}}$ matrix corresponding to $N_{\mathrm{s}}$ sites at half filling, the sampling cost becomes $\mathcal{O}(2^{\alpha N_{\mathrm{s}}})$ with a constant $\alpha\ll 1$, in contrast to the conventional algorithm with the $\mathcal{O}(2^{N_{\mathrm{s}}})$ number of summations requiring the exponential time cost. Although the computational cost is still exponential, this improvement allows us to obtain the entanglement entropy dynamics in free boson systems for more than $100$ sites. We present several examples of the entanglement entropy dynamics in low-dimensional free boson systems.