and from t$_{n}$ to t$_{n}$$_{+1}$ ≡ t$_{n}$ + Δ t . where The same procedure in Appendices A and B can be applied to the D$_{1}$ Ansatz 36 with variational parameters ψ$_{n}$ and λ$_{qn}$ , and one can obtain For the D$_{2}$ Ansatz, α represents ψ$_{n}$ or λ$_{q}$ , and is ob- tained by Eqs. (A5) and (A6). Then, by calculating k$_{1}$ , k$_{2}$ , k$_{3}$ and k$_{4}$ , one can get the parameters α ( t$_{n}$$_{+1}$ ) for the next time step using the fourth-order Runge-Kutta method. The linear absorption spectrum of the exciton-phonon system studied in this paper is calculated by where ˆ P is the polarization operator with Here µ is the transition dipole matrix element for a single site, and | n 〉$_{ex}$ is the exciton state at the n th site | n 〉$_{ex}$ ≡ ˆ a † $_{n}$| 0 〉$_{ex}$ . Substituting Eq. (E3) into Eq. (E2) one obtains There are various numerical approaches to solve cou- pled differential equations such as Eqs. (A5) and (A6) for the D$_{2}$ Ansatz and Eqs. (B5), (B6) and (B7) for the ˜ D Ansatz. One way is to transform these equations into the Volterra integral equations 37 for minimization purposes, and then use a nonlinear optimization method such as the Newton-Raphson method to solve them. Another approach is to solve the time-dependent differential equa- tions directly. The latter method has a much higher com- putational efficiency, but it requires a higher precision in the single iterative time step since the computational er- ror may accumulate as the number of the iterative steps increases. In this paper, we use the Runge-Kutta fourth-order method 38,39 to solve Eqs. (A5) and (A6) and Eqs. (B5), (B6) and (B7). The Runge-Kutta fourth-order method is widely used to solve differential equations. Its single step error is fifth order, i.e., O (Δ t $^{5}$). The algorithm for this method can be described as follows: where the variational parameters ψ$_{n}$ ′ − $_{n}$( t ) ( n ′ = 0 , 2 , ..., N − 1) and λ$_{q}$ ( t )( q = − N/ 2+1 , − N/ 2+2 , ..., N/ 2) have the following initial values and These coupled differential equations can be numeri- cally solved by the fourth-order Runge-Kutta method, and work on this is now in progress. Appendix C: The D 1 trial state Appendix E: Linear absorption spectrum by D 2 and ˜ D Ans¨ atze wherein e − i ˆ $^{Ht}$| n 〉$_{ex}$| 0 〉$_{ph}$ can be approximated by a Davy- dov trial state, for example, by a D$_{2}$ trial state: 13 f ( t$_{n}$, α ( t$_{n}$ )) = dα ( t$_{n}$ ) dt$_{n}$ − i ˙ ψ$_{n}$ ( t ) = i 2 N − $^{1}$ψ$_{n}$ ∑ q [ ˙ λ$_{qn}$λ ∗ qn − c.c. ] + Jψ$_{n}$$_{+1}$S$_{n,n}$$_{+1}$ + Jψ$_{n}$$_{−}$$_{1}$S$_{n,n}$$_{−}$$_{1}$ − N − 1 / $^{2}$ψ$_{n}$ ∑ q g$_{q}$ω$_{q}$ [ λ$_{qn}$e iqn + c.c. ] − N − $^{1}$ψ$_{n}$ q ω$_{q}$ | λ$_{qn}$ | 2 (C1) − iN − 1 | ψ$_{n}$ | 2 ˙ λ$_{qn}$ ( t ) = + N − $^{1}$Jψ ∗ $_{n}$ψ$_{n}$$_{+1}$ ( λ$_{q,n}$$_{+1}$ − λ$_{qn}$ ) S$_{n,n}$$_{+1}$ + N − $^{1}$Jψ ∗ $_{n}$ψ$_{n}$$_{−}$$_{1}$ ( λ$_{q,n}$$_{−}$$_{1}$ − λ$_{qn}$ ) S$_{n,n}$$_{−}$$_{1}$ − N − $^{1}$| ψ$_{n}$ | $^{2}$ω$_{q}$ λ$_{qn}$ − N − 1 / $^{2}$| ψ$_{n}$ | $^{2}$g$_{q}$ ω$_{q}$e − iqn (C2) ˜ F ( ω ) = 1 π Re ∫ ∞ 0 F ( t ) e $^{iωt}$dt (E1) F ( t ) = $_{ph}$〈 0 |$_{ex}$〈 0 | ˆ P e − i ˆ Ht ˆ P $^{†}$| 0 〉$_{ex}$| 0 〉$_{ph}$ (E2) ˆ P = µ ∑ n ( | n 〉$_{ex ex}$〈 0 | + | 0 〉$_{ex ex}$〈 n | ) (E3) Appendix D: Numerical Details F ( t ) = µ 2 ∑ n ∑ m $_{ph}$〈 0 |$_{ex}$〈 m | e − i ˆ $^{Ht}$| n 〉$_{ex}$| 0 〉$_{ph}$ (E4) e − i ˆ $^{Ht}$| n 〉$_{ex}$| 0 〉$_{ph}$ ≈ ∑ n ′ ψ$_{n}$ ′ − $_{n}$( t ) | n $^{′}$〉$_{ex}$ × exp( ∑ q [ λ$_{q}$ ( t ) ˆ b † q − λ ∗ $_{q}$( t ) ˆ b$_{q}$ ]) | 0 〉$_{ph}$ (E5) ψ$_{n}$ ′ − $_{n}$( t = 0) = δ$_{n}$ ′ − $_{n}$, (E6) λ$_{q}$ ( t = 0) = 0 . (E7) k$_{1}$ = f ( t$_{n}$, α ( t$_{n}$ ))Δ t k$_{2}$ = f ( t$_{n}$ + 1 2 Δ t, α ( t$_{n}$ ) + 1 2 k$_{1}$ )Δ t k$_{3}$ = f ( t$_{n}$ + 1 2 Δ t, α ( t$_{n}$ ) + 1 2 k$_{2}$ )Δ t k$_{4}$ = f ( t$_{n}$ + Δ t, α ( t$_{n}$ ) + k$_{3}$ )Δ t α ( t$_{n}$$_{+1}$ ) = α ( t$_{n}$ ) + k$_{1}$ 6 + k$_{2}$ 3 + k$_{3}$ 3 + k$_{4}$ 6 + O (Δ t $^{5}$)