We present a theory that incorporates the vibrational degrees of freedom in a high-order harmonic generation (HHG) process with ultrashort intense laser pulses. In this model, laser-induced time-dependent transition dipoles for each fixed molecular geometry are added coherently, weighted by the laser-driven time-dependent nuclear wave packet distribution. We show that the nuclear distribution can be strongly modified by the HHG driving laser. High-order harmonic generation (HHG) has attracted a great deal of attention over the past two decades for its application as a tabletop coherent extreme ultraviolet source [1,2] and as a source of attosecond pulses [3,4]. In recent years, it has been shown that HHG signals also encode information about the target [5][6][7][8]. Since driving laser pulses as short as a few femtoseconds are available, HHG spectroscopy has been perceived as a possible tool for probing chemical processes that evolve in fewfemtosecond time scales. Indeed a few pioneering pumpprobe experiments have been performed so far on different targets, e.g., SF 6 [9], N 2 O 4 [10], Br 2 [11,12], and NO 2 [13], where a dynamic system is initiated either by IR light or its second harmonics, and the time evolution of the system is probed by observing the high harmonics generated by another IR light at varying time delays. Unfortunately, with no accurate theories available, these experiments have been interpreted in terms of simple models, thus leaving out many interesting reaction dynamics buried in the measured data.HHG is a nonlinear process. It is not easily amendable to full ab initio calculations, especially for polyatomic molecules. There have been several attempts to include the nuclear degrees of freedom in HHG theory [14][15][16][17][18][19][20]. They are based on simplified models or direct solution of time-dependent Schrödinger equation (TDSE) for simple systems. Most of the models have not been fully calibrated; thus, their validity is not known. While there has been some success, for example, the prediction [14] of isotope effect in HHG and its experimental confirmation [21], on the whole, it is fair to say that there still exists no reliable theoretical tool for calculating HHG spectra from a timeevolving molecular system. In fact, even within the fixednuclei approximation, there are few reliable calculations of HHG from molecules.The goal of this Letter is twofold. First, we develop a general theory for HHG which includes the vibrational degrees of freedom; see Eqs. (1) and (2) below. Second, as an application we combine this theory with the recently developed quantitative rescattering (QRS) theory [7,8] to calculate HHG from vibrating N 2 O 4 , which was the subject of a recent experiment by Li et al. [10]. Developed initially for the case of fixed nuclei, the QRS is computationally efficient and has been well tested [22,23]. Furthermore, with the inclusion of macroscopic propagation effect [24][25][26], the predicted HHG spectra based on QRS have been shown to agree well with experiments for d...