  Developed originally for the Holstein polaron, the Davydov D1 ansatz is an
efficient, yet extremely accurate trial state for time-dependent variation of
the spin-boson model [J. Chem. Phys. 138, 084111 (2013)]. In this work, the
Dirac-Frenkel time-dependent variational procedure utilizing the Davydov D1
ansatz is implemented to study entanglement dynamics of two qubits under the
influence of two independent baths. The Ohmic spectral density is used without
the Born-Markov approximation or the rotating-wave approximation. In the strong
coupling regime finite-time disentanglement is always found to exist, while at
the intermediate coupling regime, the entanglement dynamics calculated by
Davydov D1 ansatz displays oscillatory behavior in addition to entanglement
disappearance and revival.
