When the exact factorization meets conical intersections...

Capturing nuclear dynamics through conical intersections is pivotal to understand the fate of photoexcited molecules. The concept of a conical intersection, however, belongs to a specific definition of the electronic states, within a Born–Huang representation of the molecular wavefunction. How would these ultrafast funneling processes be translated if an exact factorization of the molecular wavefunction were to be used? In this article, we build upon our recent analysis [B.F.E. Curchod, F. Agostini, J. Phys. Chem. Lett. 8, 831 (2017)] and address this question in a broader perspective by studying the dynamics of a nuclear wavepacket through two types of conical intersections, differing by the strength of their underlying diabatic coupling. Our results generalize our previous findings by (i) showing that the time-dependent potential energy surface smoothly varies, both in time and in position, between the corresponding diabatic and adiabatic potentials, with sometimes more complex features if interferences are observed, (ii) highlighting the non-trivial behavior of the time-dependent vector potential and the fact that it cannot be gauged away in general, and (iii) justifying some approximations employed in the derivation of a mixed quantum/classical scheme based on the exact factorization.


Introduction
Conical intersections (CIs) [1][2][3][4][5][6][7][8] are often invoked to interpret relaxation processes undergone by photoexcited molecules. CIs are prototypical examples of the breakdown of the Born-Oppenheimer (BO) approximation, as they represent efficient funnels [9][10][11] for population transfer between electronic states, mediated by nuclear motion. They are regions of configuration space where the adiabatic potential energy surfaces (PESs) are degenerate and exhibit, within the so-called branching place, a typical double-cone shape. These curious features of adiabatic PESs have been widely studied in the chemical physics literature, not only for their role in nonadiabatic processes but also for the effect that the related Berry phase (a topological phase) has on adiabatic phenomena, for instance occurring purely in the electronic ground state [2][3][4][5][6][7][8][11][12][13][14][15][16][17][18][19][20]. The concept of CI appears as consequence of the particular choice in our description of the electronic system at a given nuclear configuration. In the diabatic basis, CIs and Berry-phase effects disappear. Unfortunately, the diabatic basis is not rigorously defined as the set of eigenstates of an electronic hermitian operator [21][22][23], meaning that sometimes more chemical arguments are needed to define (quasi-)diabatic states. Despite the obvious challenges, both fundamental and numerical, encountered when dealing with CIs, the adiabatic representation is commonly employed to perform trajectorybased on-the-fly nonadiabatic molecular dynamics simulations of photoinduced phenomena (sometimes coupled with a local diabatization [24]). Perhaps, this choice is related to the fact that an alternative to the adiabatic framework, more rigorous than switching to the diabatic representation, has not yet been exhaustively investigated. It is maybe important to realize at this stage that the description of nonadiabatic processes proposed in the previous paragraphs -and its corresponding vocabularydirectly emanates from a Born-Huang representation of the molecular wavefunction. Such representation offers the well-known picture of nuclear amplitudes evolving on coupled potential energy surfaces and transferring from one electronic state to the other in coupling regions, possibly at a conical intersection in the adiabatic representation. Hence, instead of questioning the electronic representation itself within this picture, one could prefer to change the overall representation for the molecular wavefunction and see how the deactivation of an electronically-excited molecule takes place. In the following, we will make use of the exact factorization (EF) of the molecular wavefunction [25,26], and this paper aims at discussing such alternative theoretical framework in the context of nonadiabatic phenomena through conical intersections.
Recently [27], we have employed the EF to analyze nuclear dynamics at a CI in a model potential for the photoisomerization of retinal. We have concluded that no peculiar, i.e., singular, behavior that can be traced back to the CI arises either in the evolution of the nuclear wavepacket or in the time-dependent potentials that drive its evolution. This study was the first detailed analysis in the time domain of the properties of the EF in connection to CIs. (The EF can also be performed in the time-independent picture [28][29][30][31][32][33][34][35], and analyses of CI in this framework were also proposed [36][37][38][39].) In the framework of the EF, nuclear dynamics is generated by a Hamiltonian containing scalar and vector potentials that are time-dependent, since they represent the excited-state (dynamical) effect of the electrons. Formally, this structure is equivalent to an adiabatic picture, where the static adiabatic PES replaces the time-dependent scalar potential, for instance the ground-state PES, and the time-dependent vector potential reduces to the, once again static, nonadiabatic coupling vectors. The relation is not only formal, as in the limit of the electron-nuclear mass ratio µ = m/M tending to zero [40][41][42], the adiabatic quantities are recovered from the exact-factorization properties. However, it seems that degeneracies, conical shapes, or singularities are not so easily encountered for finite values of the electron-nuclear mass ratio. Hence, is the EF a particularly suitable framework to address "dynamics at conical intersections"? To address this question, we are going to analyze the dynamics around a CI using a model system, and for different topologies of the PESs. We will identify general features and similarities in the properties of the time-dependent potentials, that we believe will represent essential points of reflection for future, and possibly more practical, method developments.
In the following, we briefly introduce the EF in Section 2.1, present the model object of our analysis in Section 2.2, and provide computational details in Section 2.3. The following Sections are devoted to the comparison, based on the simulated dynamics, between the adiabatic and diabatic framework, arising from the Born-Huang expansion of the molecular wavefunction, which is presented in Section 3.1, and the EF framework, discussed in Section 3.2. Conclusions are reported in Section 4.

Theoretical background 2.1 Exact factorization
We introduce here the basics of the exact factorization of the electron-nuclear wavefunction [25,26]. In this framework, we use a specific Ansatz for the molecular wavefunction that departs from the most commonly employed Born-Huang representation, namely The molecular wavefunction, Ψ (r, R, t), is the solution of the time-dependent Schrödinger equation with Hamil-tonianĤ =T n +Ĥ BO . It contains the nuclear kinetic energy,T n , and the BO Hamiltonian,Ĥ BO , which is the sum of the electronic kinetic energy and all interaction potentials. The symbols r, R stand for electronic and nuclear positions, respectively. In the EF Ansatz, χ(R, t) is the nuclear wavefunction, whereas Φ R (r, t) is an electronic factor, that depends parametrically on the nuclear positions and satisfies the partial normalization condition (PNC) The symbol · r indicates an integration over electronic coordinates only. The PNC is pivotal to interpret |χ(R, t)| 2 as the marginal probability of finding the nuclear configuration R at time t, and |Φ R (r, t)| 2 as the conditional probability of finding the electronic configuration r at time t for a given nuclear configuration R.
The existence of the wavefunctions in equation (1) has been proved in references [25,26], as well as their uniqueness up to a (R, t)-dependent gauge transformation, where θ(R, t) is some real function of the nuclear coordinates and time.
Equations of motion for Φ R (r, t) and χ(R, t) result from the stationary variations [43] of the quantum mechanical action and read Nn ν=1 Importantly, the PNC was enforced by means of Lagrange multipliers [44,45]. As discussed in the Introduction, the nuclear equation (6) is clearly a time-dependent Schrödinger equation that contains time-dependent vector A ν (R, t) and scalar (R, t) potentials. The termÛ coup en [Φ R , χ] in equation (5) is coined electron-nuclear coupling operator [46], (R, t) is the time-dependent potential energy surface (TDPES) [47][48][49][50][51][52], and A ν (R, t) the time-dependent vector potential (TDVP) [27], All these terms act as an exact coupling between the electrons and the nuclei. The electron-nuclear coupling operatorÛ coup en [Φ R , χ] depends on the nuclear wavefunction and acts on the parametric dependence of Φ R (r, t) as a differential operator. Hence, this "pseudo-operator" includes the coupling to the nuclear subsystem beyond the parametric dependence found in the BO Hamiltonian. χ(R, t) is a genuine nuclear wavefunction since it leads to an N -body nuclear density, and an N -body nuclear current-density, which reproduce the true nuclear N -body density and current-density [26] obtained from the full wavefunction Ψ (r, R, t).

Models
The CI model employed in this work is based on the work by Villani et al [53]. The 2-by-2 Hamiltonian is given in the diabatic (D) basis, with "potential energy" terms This model describes the crossing of two similar parabola, one being slightly displaced both in the x direction and in energy. Adiabatic potential energies are produced by diagonalizing the diabatic potential matrix, and will lead to the appearance of a CI whose shape can be altered by modifying the off-diagonal elements H 01 (x, y) and H 10 (x, y). The parameter γ, introduced in the definition of the off-diagonal elements, can be used for this purpose by tuning the strength of the diabatic coupling between the two parabola. Two different regimes will be used in this work: weak coupling with γ = 0.01 and strong coupling with γ = 0.08. In the weak regime, the diabatic states are only weakly coupled and, upon diagonalization, give rise to a localized conical intersection (Fig. 1, upper panel). On the other hand, the intense interaction between the diabatic states in the strong-coupling regime results in adiabatic PESs differing importantly from the diabatic ones ( Fig. 1, lower panel), i.e., from the diagonal elements of the potential energy matrix. A cut through the x coordinate at fixed y = 0.0 bohr (upper panel of Fig. 2) reveals that, along this particular axis, the strong and weak coupling lead to similar adiabatic PESs. The picture dramatically changes by slightly deviating from the axis of the conical intersection and performing the cut at y = 0.333 bohr. In this second cut (middle panel of Fig. 2), the adiabatic PESs around the CI are only moderately altered as compared to the diabatic ones for the weak-coupling regime (dotted lines). Conversely, the strong diabatic coupling leads to a significant change in the topology of the adiabatic PESs with the formation of a barrier in S 1 and a deep minimum in S 0 (dashed lines). Such a difference in the adiabatic PESs resulting from the two different coupling regimes is further visible by performing a cut through y, this time, for a fixed value of x = 3.0 bohr (lower panel of Fig. 2). Once more, the adiabatic surfaces in the weak-coupling case only slightly deviate from the diabatic behavior (continuous lines), whereas a pronounced double minimum can be observed for the adiabatic ground-state surface (dashed lines) for the strong-coupling regime.
In the following, we will investigate the dynamics of a nuclear wavepacket in both situations. Even if the nuclear wavepacket dynamics is initiated in the same way for the two regimes (see Sects. 2.3 and 3), the overall dynamics will strongly differ as a result of the coupling strength between the surfaces. These differences were extensively discussed in the context of the adiabatic and diabatic representations [53,54], and we will here present the viewpoint of the EF.

Computational details
The results presented in this article were obtained in a two-step process. First, we solve the full timedependent Schrödinger equation numerically using a splitoperator formalism [55] in a diabatic basis and the 2-by-2 Hamiltonian described in Section 2.2. The initial (diabatic) nuclear wavepacket is taken as Gaussian, with widths σ x = 0.15 and σ y = 0.197, and is centered at R ini = (x ini , y ini ) = (2.0, 0.0). The dynamics is initiated in the second adiabatic state (S 1 ). We then reconstruct the TDPES and the TDVP from the time-dependent nuclear wavefunctions χ l (R, t), with l = 1, 2, using the relation In other words, if the full wavefunction is expanded in the diabatic basis, and considering the electronic states being orthogonal and normalized, the nuclear densities is simply reconstructed as the sum of "projected densities" |χ l (R, t)| 2 . In the following, we fix the gauge by imposing that the phase of the nuclear wavefunction S(R, t) = 0, i.e., the nuclear wavefunction given by is real and non-negative in this gauge. The coefficients of the electronic wavefunction in the diabatic basis are then given as χ l (R, t)/χ(R, t), by virtue of the EF Ansatz equation (1). The TDPES and the TDVP can be decomposed as well in the diabatic basis. To this end, let us first give the explicit expression of the TDPES, that follows from its definition in equation (8), namely The gauge transformations (Eq. (3)) only affect the last term on the right-hand side. Therefore, the TDPES can be decomposed in three gauge-invariant (GI) contributions, , and one gauge-dependent part (GD), GD (R, t). We express them in the diabatic basis for actual calculations, namely as if the electronic time-dependent wavefunction is expanded in the diabatic basis, as the property that nonadiabatic coupling vectors are identically zero in the diabatic basis avoids incurring in singularities at the CI when computing the TDPES and the TDVP. We can write: where the symbol H kl (R) indicates the elements of the electronic Hamiltonian (the potential matrix in Eq. (10)) in the diabatic basis and l, k label the states (see Eq. (11)). The vector potential is given by While the TDPES and TDVP are expanded here in a diabatic basis, it is critical to note that they do not depend on any particular choice of electronic representations.

Results and discussion
In the first part of this Section, we discuss the nonadiabatic dynamics of the two models presented above in a purely Born-Huang perspective, and we owe to translate them in an EF picture in the second part.

Born-Huang picture
As described in Section 2.2, the initial total wavefunction for both models is given by with χ H00 (R, t 0 ) a Gaussian wavepacket on the diabatic surface H 00 at R = (2.0, 0.0) with no initial momentum. In an adiabatic picture, this initial wavefunction represents a nuclear wavepacket initiated on the excited adiabatic state S 1 (see Fig. 2, upper panel for the relation between adiabatic/diabatic states). The nuclear wavepacket immediately starts moving towards the CI as a result of the steep gradient in the x direction of the S 1 PES in this region (Fig. 2), both in the weak-and the strong-coupling case. In the weak-coupling regime, where the CI is spatially very localized, the wavepacket enters the nonadiabatic region after 500 au, leading to a rapid transfer of nearly 85% of the S 1 population to S 0 (Fig. 3, continuous lines). In the strong-coupling case, the picture is slightly different as the coupling region extends more around the R CI = (3.0, 0.0) region. This results in an earlier transfer of the population towards S 0 that stabilizes at around 93% (Fig. 3, dashed lines). Perhaps more interesting is the effect of the coupling strength on the spatial extent of the nuclear wavepacket over time, reflecting the shape of the adiabatic potential surfaces. In the weak-coupling regime (upper panel of Fig. 4), the wavepacket extends on the x coordinate and hits the CI (t = 1000 au), transferring population to the lower state (grey contour lines in Fig. 4), while a small contribution stays on the S 1 state (red contour lines in Fig. 4). No important distortions of the wavepacket are observed, besides a splitting between the S 0 and S 1 contributions at later times (t = 2000 au). After the passage through the CI, the S 1 contribution immediately enters a repulsive part of the adiabatic PES, while the contribution evolving on S 0 further continues its relaxation towards the S 0 minimum (see the upper panel of Fig. 2). The picture is vastly different for the strong-coupling regime. As observed in the lower panel of Figure 4, the wavepacket spreads more towards the y direction as a result of the distortions of the adiabatic PESs (Fig. 2), even before reaching the exact location of the CI due to the early transfer of population towards S 0 . Upon passage through the CI (t = 1000 au), the amplitudes on the different states strongly interfere, leading to a complex oscillatory structure of the nuclear wavepacket at later times (t = 2000 au). We notice here that, as it is evident in Figure 4 at t = 1000 au the nuclear density itself does not have a node at the CI. Only its "projected" contribution onto S 1 presents a node at the CI, consequence of the fact that at the exact location of the CI, population is funnelled from S 1 to S 0 . In this sense, nuclear dynamics does not reveal any signature of the CI.

Exact-factorization picture
In the following, we will study how the critical quantities of the EF, namely the TDPES and the TDVP, portray the nonadiabatic dynamics in the two coupling regimes. Our analysis focuses on the two time snapshots, t = 1000 au and t = 2000 au, highlighted in Section 3.1.

Time-dependent vector potential -A(R, t)
Let us first observe how the TDVP, A(R, t) = A(x, y, t) in the present case, behaves at the time of the passage through the CI (t = 1000 au). In the chosen gauge, this 2-dimensional vector field is where α = x, y, and Ψ (R, t)|∂ α |Ψ (R, t) r (divided by the nuclear mass) is the Cartesian α-component of the nuclear velocity field. The additional term [25] found in the definition of the TDVP and depending on the gradient of the nuclear phase S(R, t) is identically zero in the chosen gauge.
In the weak-coupling regime (Fig. 5, upper panel), the TDVP exhibits a simple behavior when the wavepacket reaches the nonadiabatic region: it mainly points towards larger x, with a strength increasing along x, since the nuclear wavepacket itself, after being initiated on S 1 , moves rapidly towards the CI region without spreading significantly in the y-direction. We recall that in the chosen gauge the vector potential is the nuclear velocity field, thus it should not come as a surprise that it illustrates the evolution of the nuclear wavepacket. The TDVP resulting from the nonadiabatic dynamics in the strong-coupling regime shows a very different behavior at t = 1000 au, and reflects the complex nuclear dynamics around the nonadiabatic region. In contrast with the weak-coupling case, the components of A(R, t = 1000) reflect the nuclear extension in the y-direction (as described in Sect. 3.1) and its strength anticipates the structure resulting from wavepacket interferences.
After the passage through the CI (t = 2000 au), the TDVP in the weak-coupling regime still mostly points towards larger x (upper panel of Fig. 6), but now with some variations in the region 3.5 < x < 4.5 bohr. Additionally, in the region 4.0 < x < 4.5 bohr, the magnitude of A(R, t) is lower than at larger values of x, reflecting the change in behavior of the nuclear wavepacket: in the Born-Huang picture, we observed that the nuclear wavepacket component on S 0 overtakes the one on S 1 at t = 2000 au, which suffers the repulsive part of S 1 . Therefore, the nuclear component evolving on S 1 , localized in the region 4.0 < x < 4.5 bohr, moves slower than the one on S 0 , that spreads along larger values of x. In the region 4.5 < x < 6.0 bohr, the magnitude of A(R, t) starts decreasing, reflecting how the wavepacket slows down when it enters a more repulsive part of the potential energy surface. An even more complex behavior is observed for the TDVP in the strong-coupling regime (lower panel of Fig. 6). The magnitude of A(R, t) rapidly varies and its direction reflects the generation of the complex interference structure in the nuclear wavepacket following the transition through the nonadiabatic region. Owing to its definition (Eq. (21)), the TDVP is expected to diverge at nodal points of the nuclear density. In our simulation we observe, in fact, that the modulus of the vector potential increases in very localized regions, that are related to the minima of the nuclear density contribution (almost fully) propagating on S 0 . Such regions of low density arise from the interfering portions of the S 0 wavepacket, thus eventual singularities of the vector potential are, in this case, clearly unrelated to the CI, which we recall is localized at (3.0, 0.0), well outside the range of (x, y) values shown in Figure 6.
It is very important to stress that the EF does not lead to an absence of geometric phases [38]. We can compute the circulation of the vector potential as γ(t) = 1 C A(R, t) · dR, along a closed path C. According to this expression, we observe that the value of γ(t), the geometric phase, is in general non-zero, independently of the fact that the closed path C encircles or not a CI. This observation is in agreement with the findings of reference [38] still in the context of the EF but in the time-independent case. More importantly, the value of γ(t) depends on the path, and is not quantized, that is, γ(t) can have any value from 0 to 2π. Such an observation is critical for the EF as it implies that A(R, t) cannot be gauged away. The TDVP has, in general, a non-zero curl, meaning that it cannot be written as the gradient of a scalar function.

Time-dependent potential energy surface -(R, t)
The TDPES is maybe one of the most exciting features of the EF, as it allows to clearly interpret the behavior of a nuclear wavepacket in nonadiabatic conditions. For instance, different shapes of the TDPES connected by "abrupt" changes (or steps) are symptoms of a spatial splitting of the nuclear density; strong oscillations, on the other hand, while difficult to clearly decompose, are a signature of interferences occurring between two or more portions of the nuclear wavepacket interacting in a nonadiabatic region. The TDPES encodes signatures of the dynamics, and helps us in identifying the type of nonadiabatic effects most strongly affecting the behavior of the nuclear wavefunction.
In what follows, we first focus our attention on GI1 (R, t), whose expression in the diabatic basis in given in equation (15). We stress once again that the TDPES does not depend on a particular representation for the electronic system, as further highlighted later, and its expression in the diabatic basis is only used here for numerical convenience.
At t = 1000 au, in the weak-coupling regime, GI1 (R, t) simply bridges the S 1 adiabatic PES (for x < x CI ) with the S 0 one (for x > x CI ), as shown in Figure 7 (upper panel). In this sense, the first GI part of the TDPES exhibits a diabatic behavior, connecting smoothly the two BO PESs in the direction of the evolution of the nuclear wavepacket. Furthermore, GI1 (R, t) does not have singularities or discontinuities in the region of the CI. The GI1 (R, t) component of the TDPES obtained during the strong-coupling dynamics highlights an important feature of the TDPES: its behavior is not tied to a diabatic or an adiabatic picture, but can accommodate both worlds naturally. The part of GI1 (R, t) close to the CI exhibits a similar diabatic character as in the weak-coupling case, but this character evolves to a more adiabatic one when we look at y values away from the CI. In these regions, the strong diabatic coupling results in important distortions of the adiabatic PESs that now differ largely from their diabatic counterparts (see Fig. 2, lower panel). Hence, GI1 (R, t) exhibits both a diabatic and a pronounced adiabatic character at the same time. This is shown in Figure 7 (lower panel), but it appears more clearly in the cut along the y coordinate in the axis of the CI, x = 3.0 bohr (Fig. 8, lower panel), where GI1 (R, t) exhibits a mixed adiabatic/diabatic etiquette.
The relation between GI1 (R, t) and the location and character associated to the nuclear wavepacket is explicit in the cuts shown in Figure 8. In some regions, the nuclear wavefunction has a mixed character, namely the two projected contributions of the nuclear density on the adiabatic states are non-zero in the same region. This situation occurs at t = 1000 au in the weak-coupling regime, as shown in Figure 8 (upper panel). In such cases, GI1 corresponds to a mean field potential, between (S0)

BO (R) and
(S1) BO (R). In the more complex case of strong coupling, the nuclear wavepacket near the CI has a mix electronic character, reflected in an average potential. Further away from the CI, the nuclear density is dominated by the S 0 contribution, and the TDPES resembles now the ground-state adiabatic PES.
In Figure 8, both adiabatic contributions to the nuclear density present a peculiar behavior at the CI, in the weakcoupling and in the strong-coupling regime. The S 1 part has (what seems to be) a node, while the S 0 part has a sharp peak, but the two "singular" contributions cancel each other perfectly, yielding a smooth total density [56]. At the CI, the S 1 population is identically zero, due to the infinitely large nonadiabatic coupling to S 0 .
Let us now focus on the remaining two GI terms of the TDPES. We first recall the analytic expression, showing that GI2 (R, t) is non-negative, while GI3 (R, t) is non-positive (it is defined as minus the squared value of the vector potential). In a previous work [27] as well as in the weak-coupling regime studied here, we observe that the sum of these contributions is almost zero over the whole range of x and y. This is shown in the cut along the y coordinates in Figure 8 (upper panel). In the strong-coupling case, GI2 (R, t) and GI3 (R, t) do not seem to exactly cancel each other in the region around the CI (Fig. 8, lower panel). Further studies that prove and justify analytically this behavior are indeed necessary to predict a general trend. However, the numerical validation presented here provides sufficient evidence to support some of the approximations introduced in the derivation of the coupled-trajectory mixed quantum-classical (CT-MQC) algorithm [57][58][59][60]. CT-MQC is a numerical scheme that solves within a quantum-classical approximation the electronic and nuclear equations of the EF, equations (5) and (6). In CT-MQC, the TDPES has simply been approximated as GI1 (R, t) + GD (R, t), with the aim to avoid calculations of second-order derivatives of the electronic wavefunction (present in GI2 ) while keeping the gauge invariance of the neglected contributions to the TDPES (neglect of GI3 in relation to the neglect of GI2 ). According to our numerical observations, the combined effect of GI2 (R, t) and GI3 (R, t) does not have a strong R-dependence (in some cases it even appears to be zero). Therefore, their effect on the classical force computed from the gradient of the TDPES can be considered, to a good approximation, negligible. As stated above, further studies are required to investigate the properties of GI2 (R, t) and GI3 (R, t), and to identify the consequences on nuclear observables, e.g., on the nuclear kinetic energy [26]. This analysis further justifies the focus on the behavior of GI1 (R, t), as it is clearly the GI contribution to the TDPES that mostly affects the dynamics. Additionally, in our choice of gauge, the GD part of the TDPES GD (R, t) (see Eq. (18)) is mostly constant for both coupling regimes and would lead to a rigid shift of the GI parts of the TDPES, as it is shown in both panels of Figure 8 for a cut along the y axis for a fixed value of the x coordinate.
Let us now add up all the GI contributions to the TDPES at t = 2000 au (Fig. 9). 1 As observed previously in our analysis of the TDVP, the full TDPES exhibits features that drive the complex dynamics of the nuclear wavepacket. Hence, the pattern leading to a splitting of the nuclear wavepacket in the weak-coupling regime -explained in a Born-Huang picture by the S 1 nuclear component separating from the S 0 nuclear component -is clearly visible from the TDPES (Fig. 9, upper panel). The TDPES is composed of two main regions (3.7 < x < 4.6 bohr and 4.6 < x < 5.5 bohr), with a pronounced change of behavior at their interface. A sharp repulsive potential in the x direction composes the first region, supplemented by a central repulsive component at y = 0 bohr. Conversely, the second region only shows a slowly increasing potential towards larger x value. This strong variation of the TDPES is responsible for the splitting of the nuclear wavepacket into two components. Even though the dynamics is more involved in the strongcoupling case, especially because of interferences observed in the S 0 wavepacket after the passage through the CI, the TDPES still clearly modulates the shape of the nuclear density, as shown in Figure 9 (lower panel). Regions where the TDPES is large are associated to regions where the nuclear density is small (tending towards a node), whereas the series of minima observed in the TDPES (at around x = 3.9 and for y varying between -0.5 and 0.3) creates a multi-peaked nuclear density. The oscillatory features of the TDPES in the strong-coupling case can be interpreted as a two-dimensional generalization of our previous analysis [51] on the effect of interferences on the TDPES. From Figure 2 (middle panel) it is evident that when the nuclear wavepacket moving on S 1 reaches the CI, it transfers population to S 0 . However, due to the shape of the adiabatic PES, the incoming density is partially trapped in a potential well. The result is that the wavepacket in S 0 interferes with the wavepacket still incoming from S 1 and transfers some population back to S 1 . This statement is also validated by the results presented in Figure 3 (dashed lines), where we observe that after 1500 au the population of S 0 slightly decreases and the population of S 1 slightly increases. Interferences observed at t = 2000 au in the strong-coupling regime are therefore the effect of nonadiabatic interferences, propagated over time, that occur at the coupling region between the S 0 and the S 1 wavepackets.

Conclusions
The theoretical framework of the exact factorization of the electron-nuclear wavefunction has been employed to investigate the nuclear dynamics at conical intersections. The time-dependent potential energy surface and the time-dependent vector potential have been analyzed as indicators of the nonadiabatic effects influencing nuclear relaxation through a region of strong coupling between electronic states. We have pointed out key and general features of the time-dependent potentials in different coupling regimes, proving the general validity of our previous observations [27] based on the photoisomerization process of a retinal model.
Despite the fact that the potentials of the exact factorization are generalizations of the adiabatic potential energy surfaces and of the nonadiabatic coupling vectors (all static properties of the electronic system), they do not show any singular behavior usually viewed as signature of the presence of conical intersections.
The time-dependent potential energy surface presents features that, somehow, adapt to the dynamics. Upon comparison with diabatic and adiabatic potential energy surfaces, we observed that the time-dependent potential energy surafce can show properties that connect in space these two representations at a given time along the dynamics. But the TDPES can also develop properties that switch from one representation to the other over time. The nature of the time-dependent potential energy surface can, however, be more involved, and in fact we observed mean-field shapes, as well as strongly oscillatory features, symptoms of interferences between wavepackets associated to different electronic states.
The time-dependent vector potential is smooth along the whole studied dynamics, in particular in the region of coupling when the population transfer takes place between electronic states. In the chosen gauge, the vector potential equals the nuclear velocity field; thus even though it is a gauge-dependent quantity, the time-dependent vector potential has in our study an important physical interpretation and indicates how the nuclear density evolves. The vector potential is far from being a trivial property, in the sense that, in general, it is not a curl-free vector field, and it cannot be gauged away. This observation should not come as a surprise, since there is no particular reason for the nuclear velocity field to be irrotational over all configuration space. Being not irrotational, the circulation of the vector potential along a closed path in nuclear space is, once again in general, non-zero: it yields a non-trivial, path-dependent, and not quantized, geometric phase.
The exact factorization offers very rich information about excited-state molecular processes, which departs from the picture that most standard approaches provide. Our studies have always as main target the development of simpler and more efficient ways to perform simulations of nonadiabatic phenomena. Several advances in this direction have already been made. We mentioned for example the development of a new quantum-classical approach based on the exact-factorization equations that properly captures quantum decoherence effects. Another interesting development around the exact factorization is the recent proposition [61] to treat the electronic problem within a density-functional framework. Therefore, we strongly believe that future developments around the exact factorization might lead to approaches where the difficulties related to electronic-structure representations can be relieved, but by still preserving a clear physical and chemical pictures of the processes of interest.