From a mathematical perspective, it is desirable to apply the adjoint operator for back propagating the receiver wavefield for migration of data residuals for inversion. For frequency-domain direct solvers, it is straightforward to apply the adjoint operator, whereas in most real applications of time-domain finite-difference (TDFD) stencils, the forward-propagation kernel function is reused for backward propagation for simplicity. However, when applying the exact adjoint operator, the migration result will be different, especially when strong variations exist in the velocity model. Actually, these three operators (forward, backward, and adjoint of the forward) can be expressed in matrix forms under the Born approximation for acoustic wave equations. These expressions linearly relate model perturbations to recorded seismic data. Every element in the matrices is well-defined, and all involved operations, such as the imaging condition and wavefield sampling, are included in the closed-form matrix expressions. These matrix expressions provide a platform for analyzing seismic modeling and inversion via mature linear algebra methodologies and provide clear strategies for developing computer algorithms. By analyzing the similarity of the matrix expressions, one can find that the time stepping approaches for all three operators are essentially the same. Based on this observation, a new time-marching stencil can be designed to realize the TDFD adjoint operator. Compared with traditional reverse time migration, the new method using the adjoint operator can provide better image quality, especially at sharp velocity boundaries.