# Gauge-Invariant Formulation of Time-Dependent Configuration Interaction Singles Method
## Theory
### Review of CI Coefficient-Based TDCIS with Fixed Orbitals  
We write the second-quantized version of Equation (1), for the LG case, as  
\[|\Psi_{\mathrm{L}}(t)\rangle = |\Phi\rangle C_{0}(t)+\sum_{i}^{occ\ \underline{vir}}\sum_{a}^{ \underline{vir}}|\Phi_{ia}\rangle C_{ia}(t), \tag{19}\]  
where \(|\Phi_{ia}\rangle=\sum_{\
u}^{\uparrow\downarrow}\hat{c}_{a\
u}^{\dagger} \hat{c}_{i\
u}|\Phi\rangle/\sqrt{2}\). The equations of motion for the CI coefficients have been derived [2] by inserting Equation (19) into the LG TDSE, the first of Equation (14), and closing from the left with the reference and singly-excited CSFs,  
\[\langle\Phi|(\hat{H}_{\mathrm{L}}-i\partial_{t})\{|\Phi\rangle C_{0 }+\sum_{jb}|\Phi_{jb}\rangle C_{jb}\} = 0, \tag{20a}\] \[\langle\Phi_{ia}|(\hat{H}_{\mathrm{L}}-i\partial_{t})\{|\Phi \rangle C_{0}+\sum_{jb}|\Phi_{jb}\rangle C_{jb}\} = 0. \tag{20b}\]  
Conceptually more proper derivation of Equation (20) is based on the Dirac-Frenkel variational principle, which considers the Lagrangian  
\[L_{\mathrm{L}}(t)=\langle\Psi_{\mathrm{L}}|(\hat{H}_{\mathrm{L}}-i\partial_{t })|\Psi_{\mathrm{L}}\rangle, \tag{21}\]  
and requires \(\partial L_{\mathrm{L}}/\partial C_{0}^{*}=\partial L_{\mathrm{L}}/\partial C _{ia}^{*}=0\). Substituting \(\hat{H}_{\mathrm{L}}\) of Equation (10) into Equation (20), using the Slater-Condon rule for the Hamiltonian matrix elements, and noting the canonical condition \(f_{pq}=\epsilon_{p}\delta_{pq}\), the EOMs for the length gauge are derived as [2]  
\[i\partial_{t}C_{0} = \sqrt{2}E\cdot\sum_{jb}\langle\phi_{j}|\hat{\mathbf{r}}|\phi_{b} \rangle C_{jb}, \tag{22a}\] \[i\partial_{t}C_{ia} = \langle\phi_{a}|\{\sum_{b}(\hat{F}_{i}+\mathbf{E}\cdot\hat{\mathbf{r}})| \phi_{b}\rangle C_{ib}+\sqrt{2}\mathbf{E}\cdot\hat{\mathbf{r}}|\phi_{i}\rangle C_{0} \}-\mathbf{E}\sum_{j}C_{ja}\cdot\langle\phi_{j}|\hat{\mathbf{r}}|\phi_{i}\rangle. \tag{22b}\]  
where the action of the operator \(\hat{F}_{i}\) on a given orbital \(\phi\) is defined as  
\[\hat{F}_{i}|\phi\rangle=(\hat{f}-\epsilon_{i})|\phi\rangle+\sum_{j}(2\hat{W} _{\phi}^{\phi_{j}}|\phi_{i}\rangle-\hat{W}_{\phi_{i}}^{\phi_{j}}|\phi\rangle). \tag{23}\]  
References [17; 18; 19; 20; 21; 22; 23; 24] have used the same expansion in terms of fixed CSFs also in the VG case,  
\[|\Psi_{\mathrm{V}}(t)\rangle = |\Phi\rangle D_{0}(t)+\sum_{i}^{occ\ \underline{vir}}\sum_{a}^{ \underline{vir}}|\Phi_{ia}\rangle D_{ia}(t), \tag{24}\]  
and required Equation (20) to hold, with \(\hat{H}_{\mathrm{L}}\), \(C_{0}\), and \(C_{ia}\) replaced with \(\hat{H}_{\mathrm{V}}\), \(D_{0}\), and \(D_{ia}\). This is equivalent to consider the following Lagrangian,  
\[L_{\mathrm{V}}(t)=\langle\Psi_{\mathrm{V}}|(\hat{H}_{\mathrm{V}}-i\partial_{t })|\Psi_{\mathrm{V}}\rangle, \tag{25}\]  
and to require \(\partial L_{\mathrm{V}}/\partial D_{0}^{*}=\partial L_{\mathrm{V}}/\partial D _{ia}^{*}=0\), which derives  
\[i\partial_{t}D_{0} = \sqrt{2}\mathbf{A}\cdot\sum_{jb}\langle\phi_{j}|\hat{\mathbf{p}}|\phi_{b }\rangle D_{jb}, \tag{26a}\] \[i\partial_{t}D_{ia} = \langle\phi_{a}|\{\sum_{b}(\hat{F}_{i}+\mathbf{A}\cdot\hat{\mathbf{p}})| \phi_{b}\rangle D_{ib}+\sqrt{2}\mathbf{A}\cdot\hat{\mathbf{p}}|\phi_{i}\rangle D_{0}\}- \mathbf{A}\sum_{j}D_{ja}\cdot\langle\phi_{j}|\hat{\mathbf{p}}|\phi_{i}\rangle. \tag{26b}\]