Implicit finite-difference (FD) methods are widely used in seismic wave modeling for their accuracy and coarse-sampling advantages, yet balancing temporal accuracy and computational cost remains a challenge. Conventional implicit schemes achieve high spatial accuracy but remain second-order in time, requiring fine time steps to mitigate temporal dispersion. Recently developed spatiotemporal high-order implicit methods improve temporal accuracy but introduce substantial computational overhead due to the inclusion of numerous additional grid points. To address this limitation, we propose an asymmetric factorization of implicit FD operators for time-domain acoustic wave modeling. This novel approach derives FD coefficients through a joint approximation of temporal and spatial derivatives by matching the discrete wave equation’s dispersion relation, rather than approximating spatial derivatives alone. Combined with a variable-substitution-based Taylor-series expansion, our method achieved (2M+2)th-order spatial accuracy and (2N)th-order temporal accuracy while requiring only half as many additional grid points as conventional spatiotemporal implicit schemes (per spatial operator in the 2D case). Dispersion analysis demonstrated that the proposed scheme effectively suppressed numerical dispersion across a wide range of wavenumbers and Courant numbers. Stability analysis revealed that our method permitted larger time steps than conventional implicit schemes. Numerical experiments in both homogeneous and heterogeneous media validated the method’s superior accuracy and computational efficiency. The proposed asymmetric factorization offers a practical solution for seismic imaging and full-waveform inversion applications in 2D and 3D subsurface settings where both accuracy and efficiency are critical.