3.4.1. Population Dynamics

Figure 2 presents the transient time evolution of *ρ*00(*t*) during the first period of modulation by changing the time interval *δt* while keeping the number of divisions *N* constant at *N* = 40. We set parameters as *λ* = 0.01, *ωc* = 3*<sup>ω</sup>*0, and *h*¯ *ω*0 = 25 meV which shows the relaxation time *τ* ¯ *r* ≡ *ω*0*τr* ≈ 5. (The value of *h*¯ *ω*0 is chosen to be the same as the typical value for a molecular junction in Reference [28].) Setting the initial condition of the two-level system with the effective inverse temperature as *β* ¯ *s* = *h*¯ *ω*0 *kBTs* = *β* ¯(0)(≈1.07) corresponding to the stationary state for the initial temperature setting, we plot the time dependence of the population in the lower state, *ρ*00(*t*). The population *ρ*00(*t*) under the adiabatic approximation (Figure 2, red line) shows that the relevant system quickly approaches the stationary state corresponding to the temperature setting in each time interval. Setting the interval *δt* to be much larger than the relaxation time as in Figure 2A corresponding to the lower modulation frequency Ω = 0.3 THz, we find that the relevant system mostly follows the temperature modulation as the stationary state is approached, which shows the feature close to the adiabatic approximation. With decreasing interval *δt* (Figure 2B,C), we find that the relevant system does not follow the temperature modulation thus exhibiting nonadiabaticity.

**Figure 2.** Time dependence of the population in the lower state of the two-level system with changing modulation frequency : (**A**) Ω = 0.3 THz, (**B**) Ω = 1 THz, and (**C**) Ω = 5 THz with *s* = 0.01, *ωc* = 3*<sup>ω</sup>*0, *h*¯ *ω*0 = 25 meV , and *N* = 40. The time variable is scaled with *ω*0 as ¯*t* = *<sup>ω</sup>*0*t*.
