In rarefied gas flows, the spatial grid size could vary by several orders of magnitude in a single flow configuration (e.g., inside the Knudsen layer it is at the order of mean free path of gas molecules, while in the bulk region it is at a much larger hydrodynamic scale). Therefore, efficient implicit numerical method is urgently needed for time-dependent problems. However, the integro-differential nature of gas kinetic equations poses a grand challenge, as the gain part of the collision operator is non-invertible. Hence an iterative solver is required in each time step, which usually takes a lot of iterations in the (near) continuum flow regime where the Knudsen number is small; worse still, the solution does not asymptotically preserve the fluid dynamic limit when the spatial cell size is not refined enough. Based on the general synthetic iteration scheme for steady-state solution of the Boltzmann equation, we propose two numerical schemes to push the multiscale simulation of unsteady rarefied gas flows to a new boundary, that is, the numerical solution not only converges within dozens of iterations in each time step, but also asymptotically preserves the Navier-Stokes-Fourier limit in the continuum flow regime, when the spatial grid is coarse, and the time step is large (e.g., in simulating the extreme slow decay of two-dimensional Taylor vortex, the time step is even at the order of vortex decay time). The properties of fast convergence and asymptotic preserving of the proposed schemes are not only rigorously proven by the Fourier stability analysis for simplified gas kinetic models, but also demonstrated by several numerical examples for the gas kinetic models and the Boltzmann equation.