# Gauge-Invariant Formulation of Time-Dependent Configuration Interaction Singles Method
## Theory
### Evaluation of the Time Derivative of an Observable  
Let us next consider how to compute expectation value of a one-electron operator \(\langle\hat{\mathcal{O}}\rangle(t)=\langle\Psi(t)|\hat{\mathcal{O}}|\Psi(t)\rangle\), and its time derivative \(d\langle\hat{\mathcal{O}}\rangle/dt\). For exact solution of TDSE, \(\partial_{t}|\Psi\rangle=-i\hat{H}|\Psi\rangle\), the time derivative is given by  
\[\frac{d}{dt}\langle\Psi|\hat{\mathcal{O}}|\Psi\rangle = \langle\Psi|\hat{\mathcal{O}}|\partial_{t}\Psi\rangle+\langle \partial_{t}\Psi|\hat{\mathcal{O}}|\Psi\rangle \tag{38a}\] \[= -i\langle\Psi|[\hat{\mathcal{O}},\hat{H}]|\Psi\rangle, \tag{38b}\]  
known as the Ehrenfest expression (where for simplicity a trivial, explicit time-dependence of the operator has been dropped). For an approximate method, however, the Ehrenfest theorem, Equation (38b), generally does not hold, and one should explicitly evaluate the time derivative as Equation (38a). Important exceptions include those theories using time-dependent orbitals that have evolved to satisfy the time-dependent variational principle, such as time-dependent Hartree-Fock (TDHF), TD-MCSCF, and time-dependent density functional theory. See Ref. [41] for more details.  
The TDCIS expectation value of a one-electron operator \(\hat{O}\) is given [2] by  
\[\langle\Psi_{\text{L}}|\hat{\mathcal{O}}|\Psi_{\text{L}}\rangle=2\sum_{j} \langle\phi_{j}|\hat{\mathcal{O}}|\phi_{j}\rangle+\sum_{j}\langle x_{j}|\hat{ \mathcal{O}}|x_{j}\rangle+2\sqrt{2}\operatorname{Re}\left[C_{0}^{*}\sum_{j} \langle\phi_{j}|\hat{\mathcal{O}}|x_{j}\rangle\right]-\sum_{ij}\langle x_{i}| x_{j}\rangle\langle\phi_{j}|\hat{\mathcal{O}}|\phi_{i}\rangle, \tag{39}\]  
in the LG case. That for the VG is given by replacing \(C_{0}\) with \(D_{0}\) in the above equation, and for the rVG by replacing \(\{\phi_{j},x_{j}\}\) with \(\{\phi_{j}^{\prime},x_{j}^{\prime}\}\). The expression for the time derivative, in the LG case, is derived by using Equation (28) in Equation (38a) as  
\[\frac{d\langle\Psi_{\text{L}}|\hat{\mathcal{O}}|\Psi_{\text{L}}\rangle}{dt}=2 \operatorname{Re}\left[\sum_{j}\langle x_{j}|\hat{\mathcal{O}}|\hat{x}_{j} \rangle+\sqrt{2}(\hat{C}_{0}^{*}\langle\phi_{j}|\hat{\mathcal{O}}|x_{j}\rangle +C_{0}^{*}\langle\phi_{j}|\hat{\mathcal{O}}|\hat{x}_{j}\rangle)-\sum_{ij} \langle x_{i}|\hat{x}_{j}\rangle\langle\phi_{j}|\hat{\mathcal{O}}|\phi_{i} \rangle\right], \tag{40}\]  
where \(\hat{C}_{0}\equiv\partial_{t}C_{0}\) and \(|\hat{x}_{j}\rangle\equiv\partial_{t}|x_{j}\rangle\). The VG expression is also given by the above equation with \(C_{0}\) replaced with \(D_{0}\), and that for the rVG is  
\[\frac{d\langle\Psi_{\text{V}}^{\prime}|\hat{\mathcal{O}}|\Psi_{ \text{V}}^{\prime}\rangle}{dt}=2\operatorname{Re}\left[\sum_{j}\langle x_{j}^{ \prime}|\hat{\mathcal{O}}|x_{j}^{\prime}\rangle+\sqrt{2}(\hat{C}_{0}^{*} \langle\phi_{j}^{\prime}|\hat{\mathcal{O}}|x_{j}^{\prime}\rangle+C_{0}^{*} \langle\phi_{j}^{\prime}|\hat{\mathcal{O}}|x_{j}^{\prime}\rangle)-\sum_{ij} \langle x_{i}^{\prime}|\hat{x}_{j}^{\prime}\rangle\langle\phi_{j}^{\prime}| \mathcal{O}|\phi_{i}^{\prime}\rangle\right] \tag{41}\] \[+\sqrt{2}\operatorname{Im}\left[2\mathbf{E}\cdot\sum_{j}C_{0}^{*} \langle\phi_{j}^{\prime}|\hat{\mathcal{O}}|x_{j}^{\prime}\rangle+|\mathbf{A}|^{2} \sum_{j}C_{0}^{*}\langle\phi_{j}^{\prime}|\hat{\mathcal{O}}|x_{j}^{\prime} \rangle\right]-i\mathbf{E}\cdot\sum_{ij}(2\delta_{ij}-\langle x_{i}^{\prime}|x_{ j}^{\prime}\rangle)\langle\phi_{j}^{\prime}|[\hat{\mathcal{P}},\hat{\mathcal{O}} |\phi_{i}^{\prime}\rangle.\]  
Although Equations (40) and (41) look rather complicated, their evaluations are straightforward given the time derivatives of working variables \(C_{0}\), \(\{\chi_{i}\}\), etc, which are necessary, in any case, to propagate the EOMs.