# Charge migration and attosecond solitons in conjugated organic molecules

François Mauger, Aderonke S Folorunso, Kyle A Hamer, Cristel Chandre,

Mette B Gaarde, Kenneth Lopata, Kenneth J Schafer

► To cite this version:

Francois Mauger, Aderonke S Folorunso, Kyle A Hamer, Cristel Chandre, Mette B Gaarde, et al. Charge migration and attosecond solitons in conjugated organic molecules. Physical Review Research, 2022, 4 (1), pp.013073. 10.1103/PhysRevResearch.4.013073. hal-02866922v3

HAL Id: hal-02866922

https://hal.science/hal-02866922v3

Submitted on 13 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Charge migration manifests as attosecond solitons in conjugated organic molecules

François Mauger¹, Aderonke S. Folorunso², Kyle A. Hamer¹, Cristel

Chandre³, Mette B. Gaarde¹, Kenneth Lopata²,⁴, and Kenneth J. Schafer¹

¹Department of Physics and Astronomy, Louisiana State University, Baton Rouge, Louisiana 70803, USA

²Department of Chemistry, Louisiana State University, Baton Rouge, Louisiana 70803, USA

³CNRS, Aix Marseille Univ, I2M, Marseille, France

⁴Center for Computation and Technology, Louisiana State University, Baton Rouge, Louisiana 70803, USA

(Dated: January 13, 2022)

Charge migration is the electronic response that immediately follows localized ionization or excitation in a molecule, before the nuclei have time to move. It typically unfolds on sub-femtosecond time scales and most often corresponds to dynamics far from equilibrium, involving multi-electron interactions in a complex chemical environment. While charge migration has been documented experimentally and theoretically in multiple organic and inorganic compounds, the general mechanism that regulates it remains unsettled. In this work we use tools from nonlinear dynamics to analyze charge migration that takes place along the backbone of conjugated hydrocarbons, which we simulate using time-dependent density functional theory. In this electron-density framework we show that charge migration modes emerge as attosecond solitons and demonstrate the same type of solitary-wave dynamics in both simplified model systems and full three-dimensional molecular simulations. We show that these attosecond-soliton modes result from a balance between dispersion and nonlinear effects tied to time-dependent multi-electron interactions. Our soliton-mode mechanism, and the nonlinear tools we use to analyze it, pave the way for understanding migration dynamics in a broad range of organic molecules. For instance, we demonstrate the opportunities for chemically steering charge migration via molecular functionalization, which can alter both the initially localized electron perturbation and its subsequent time evolution.

# I. INTRODUCTION

The movement of electrons and holes in matter regulates many physical and chemical processes such as chemical reactions, photosynthesis and photovoltaics, and charge transfer [1, 2]. For the electrons, these dynamics can reach down to the Angstrom and attosecond spatio-temporal scales, and are commonly referred to as charge migration (CM) [3–8]. At the fastest time scales, CM is understood as the purely electronic-driven dynamics that takes place before nuclei have time to move. It can be the precursor for many of the down-stream processes mentioned above [2, 9–11], and therefore a means of understanding and ultimately steering them with the goal of charge-directed reactivity [12, 13].

The study of molecular CM is a formidable endeavor. Experimental studies require coherent probes with attosecond resolution [2, 6]. Such ultrafast probes have been enabled by the continuous progress in x-ray and table-top infrared sources over the past few decades [1]. For example, recent experiments have investigated how fast CM couples to the nuclear degrees of freedom, leading to structural changes [14, 15]. Theoretical investigations of CM necessitate models with multiple interacting electrons [16, 17], even when nuclear motion is ignored [18–22]. Furthermore, systematic studies of the underlying mechanisms responsible for regulating CM involve the analysis of systems with a large number of coupled degrees of freedom.

In this work, we leverage tools from nonlinear dynamics to study (sub-)femtosecond field-free CM unfolding in conjugated organic molecules following the

<figure><img src="image_2.png" /><figcaption>FIG. 1. Charge migration (CM) manifests as attosecond solitons in conjugated organic molecules: (top) Schematic of our investigations of CM in the $\pi$ system of a conjugated carbon chain. We initially put a localized hole at one end of the conjugated system and then (colormap) track the field-free CM electron dynamics it induces in the molecule.</figcaption></figure>

sudden creation of a localized hole in the system, a paradigm for site specific ionization. This is motivated by the fact that nonlinear dynamics [23, 24] has developed general-purpose methods for tackling and understanding the structure of high-dimensional phase spaces. Indeed, nonlinear dynamical analyses have been instrumental in many areas of atomic, molecular and optical science, including in transition-state theory of chemical reactions [25–27] and strong-field physics [28–31].

We show that periodic modes of CM that feature a localized hole that travels back and forth along the molecule's conjugated chain, as sketched in figure 1, emerge as attosecond solitons. These solitary waves rep-
resent a balance between dispersion and nonlinear effects that are driven by time-dependent multi-electron interactions. Previous studies have invoked the beating between a few molecular orbitals as the underlying mechanism that regulates CM [6, 9, 18]. Instead, here we propose a different view of the CM dynamics, directly in the time domain. Our attosecond-soliton picture provides a generic mechanism for sustained CM motions in organic molecules that does not rely on nonphysical molecular orbitals to describe CM dynamics. Notably, we find that the same molecule can support several of these CM modes with periods varying by several hundred attoseconds. We discuss the implication of our results for future theoretical and experimental CM studies, including the opportunities for chemically steering CM by molecular functionalization, both in creating the initially localized electron hole and for its subsequent time evolution.

The paper is organized as follows: Section II introduces the theoretical and computational framework we use to study CM dynamics. Section III proceeds with a detailed nonlinear dynamical analysis of CM modes in model conjugated hydrocarbons for which we can perform extensive simulations. Specifically, we investigate the phase-space structures that enable periodic CM modes and establish the solitary-wave mechanism that supports such dynamics. Then, in section IV, we build on those results to provide evidence for the same type of attosecond-soliton CM mode in full quantum-chemistry computations of the bromohexatriyne molecule. We also discuss experimental implications of our findings. Section V concludes the paper and discusses implication of our results for future CM studies.

## II. THEORETICAL AND COMPUTATIONAL FRAMEWORK

We consider field-free CM in a conjugated hydrocarbon molecular cation following the sudden creation of a localized one-electron hole in the system, e.g., as would result from site-specific ionization. Our general approach is sketched in figure 1: we initially put the localized hole at one end of the conjugated $\pi$ system and study how the ensuing CM dynamics moves the electron/hole density across the entire molecule. The delocalized $\pi$ structure of conjugated hydrocarbons has previously been shown to serve as the backbone for CM in these compounds [22]. To demonstrate the generality of our finding, we consider two levels of theory: (1) a reduced model of the $\pi$ system for which we can perform extensive computations and detailed analyses, and (2) high-performance quantum-chemistry simulations in a real molecule, where we leverage what we learned in the reduced case. The reduced $\pi$ system corresponds to a one-dimensional (1D) model of an alkene. It emulates the carbon chain conjugated structure by matching the bond lengths between the various C centers of the full 3D molecules. We detail our model $\pi$ system in Appendix 1.

In all simulations we use time-dependent density-functional theory (TDDFT) with fixed nuclei. TDDFT has been shown to successfully explain CM experiments and to reproduce correlated wave-function CM simulations [7, 19]. More broadly, TDDFT computations are commonly used in atto/femtosecond science for their ability to systematically handle large molecules with many active and correlated electrons [32, 33]. For a molecular cation with $N-1$ active electrons the TDDFT dynamics, in the Kohn-Sham (KS) formalism [34], is given by the system of one-particle equations

$$i\partial_t\phi_k(\vec{r};t) = \hat{\mathcal{H}}_{\text{eff}}[\rho](\vec{r})\phi_k(\vec{r};t), \quad (1)$$

together with the one-body density

$$\rho(\vec{r};t) = \sum_{k=1}^{N-1} |\phi_k(\vec{r};t)|^2. \quad (2)$$

Following Pauli's exclusion principle, the KS orbitals $\{\phi_k\}_k$ are orthonormal wave function in their spin and space coordinates. The one-body density provides a real-space representation of the electronic-charge distribution in the molecule. From it, the hole density $\rho_h(\vec{r};t) = \rho_N(\vec{r}) - \rho(\vec{r};t)$ is computed by taking the difference with the neutral's ground-state density $\rho_N(\vec{r})$. We stress the distinction between three related but different elements of TDDFT: (i) The KS orbitals of equation (1) are the dynamical variables of the TDDFT problem; in themselves, KS orbitals do not correspond to physically observable quantities. (ii) Molecular orbitals (MO) are time-independent wave functions; we use them as a spatial basis set to build initial conditions in the model $\pi$ system. (iii) The one-body density of equation (2) is the physical observable we are ultimately interested in for CM analyses. We provide additional details about our TDDFT simulations in Appendix 2.

In what follows, we investigate the influence of the initial hole localization on its subsequent CM dynamics. Specifically, we systematically and continuously vary the initial degree of localization of the one-electron hole around one end of the $\pi$-system structure in both the 1D and full 3D simulations. By convention we set 0% localization to correspond to the cation ground state and 100% to correspond to a fully localized hole, as sketched at the top of figure 1. Because of the different implementation of the 1D and 3D TDDFT simulations of Eqs. (1,2) (see Appendix 2), we use different avenues to introduce a variably localized initial the hole in model $\pi$ system and full TDDFT computations. We provide details about our initial hole configurations in Appendix 3. In short, for the model $\pi$ system, we build the variably-localized hole configuration with a linear combination of a few occupied and unoccupied MOs of the corresponding cation. To build the initial hole for the full TDDFT simulations, we use constrained DFT (cDFT) [35], which combines standard energy minimization techniques with the constraint of having a certain amount of the hole localized around a specific center in the molecule. This approach
allows us to consistently impose the initial-hole configu-
ration without involving *ad hoc* MO mixing. The cDFT
approach was shown to be successful in a previous study
of CM in halocarbons [22]. Following the initialization
of the hole, we compute the subsequent TDDFT dy-
namics in the full, unrestricted TDDFT framework of
equations (1,2), for both the model system and the real
molecule. Taken together, these equations correspond to
a nonlinearly coupled system, which we study using tools
from nonlinear dynamics.

III. NONLINEAR DYNAMICAL ANALYSIS

Figure 2 summarizes the result of our nonlinear analysis of CM dynamics in a reduced model of the conjugated π system of alkene/alkyne hydrocarbons. Specifically, we investigate the influence of the initial hole localization on the CM dynamics it induces in the molecule (dark orange markers). The more the initial hole is localized in the molecule the further away its electronic structure is from the ground-state distribution. This translates into a higher level of molecular excitation, which we indicate along the lower x axis. For each initial hole configuration we compute the subsequent field-free CM dynamics as given by equations (1,2). We then extract the main frequency component of any motion that moves the hole density between the two ends of the molecule – shown as the right y axis in the figure. We aim to obtain a global picture of phase space through the dependence of this frequency on the initial conditions. In nonlinear dynamics, this approach corresponds to a frequency-map analysis (FMA) [36, 37] which, for instance, has been very successful in celestial mechanics to understand the dynamics

<figure><img src="image_3.png" /><figcaption>FIG. 2. A localized one-electron hole initially introduced in the π system of a conjugated hydrocarbon leads to CM modes whose periods vary with the initial degree of localization. While varying the degree of localization of the initial hole, or equivalently the excitation energy of the molecule, we track the main frequency/period of the CM motion that moves the particle hole throughout the entire chain (dark orange dots). We compare these with direct computations of periodic CM modes in the same system (light grey dots).</figcaption></figure>

inside the solar system [38, 39]. In particular, FMA can
be used to discriminate between chaotic motions, with
a strong sensitivity to the initial conditions, and regular
ones, among them periodic motions.

A. Frequency-map analysis

For this work, we configure our FMA to focus on CM
dynamics unfolding within the first few tens of femtosec-
onds following the initial localized-hole creation. We pro-
vide further details regarding our implementation of the
FMA in Appendix 4. The results of the FMA, as shown
in figure 2, reveal several striking features. First it shows
that, as a trend, the period of the CM increases (fre-
quency decreases) with increasing molecular excitation
energy. In other words, initially more localized holes
lead to slower CM. Second, this increase is irregular and
proceeds in an almost step-wise fashion via a series of
plateaus. Within each of these plateaus, the period of
the CM is essentially independent of the excitation en-
ergy or, equivalently, of the details of the way in which
the hole was initialized.

Next, we look at the type of dynamics associated with
the plateaus in the FMA of figure 2. We have found
that these dynamics are best understood by looking at
the density contribution from the KS orbital in which we
introduce the initial localized-hole perturbation. In sec-
tion III B below we explain how this KS-orbital density
matches the CM motion of the hole we are ultimately
interested in. For illustration, in figure 3 (a,b) we show
the temporal evolution of two sample densities – see the
matching labels 1 and 2 in figure 2. Both panels reveal
qualitatively similar motions where the initial electron
density in the KS orbital, instead of spreading, remains
localized in space and periodically propagates through
the entire π system like a particle.

<figure><img src="image_4.png" /><figcaption>FIG. 3. Examples of charge-migration modes in a conjugated hydrocarbon model with a particle electron/hole that periodically moves through the π system. (a,b) Electron density contribution from the KS-orbital channel in which we introduce the initial localized-hole perturbation – see text. (c) Hole density associated with panel (b). The overhead labels ① and ② indicate the plateaus of figure 2 from which these are taken.</figcaption></figure>
Aside from the difference in the periods, a second quantitative difference between panels (a) and (b) of figure 3 is that (b) exhibits a “cleaner” propagation of the electron density along the molecule. This difference in the quality of the CM motion can be explained with the extended FMA of figure 4, by looking at the next-to-leading frequency components in the CM signal. In figure 4, in addition to the main frequency component $\Omega_1$ we identify several overtones corresponding to an anharmonic component $\Omega_2$ and two harmonic ones $2\Omega_1$ and $3\Omega_1$ (seen above $\approx 3.5$ eV). Notably, the plateau labeled 1 has both $\Omega_1$ and $\Omega_2$ components, which manifests as a faint periodicity in the evolution of the electron density. On the other hand, in plateau 2, only multiples of the lowest frequency are visible, which translates into a very clean quasi-periodic migration. Generally speaking, figure 4 shows that we have a single fundamental CM frequency plus its harmonics, and hence clean CM motion, for initial hole localization above about 60% in the model $\pi$ system. We have checked that we observe qualitatively similar results in the other plateaus of the FMA.

<figure><img src="image_5.png" /><figcaption>FIG. 4. Extended frequency-map analysis (FMA) for the model $\pi$ system showing the leading frequencies of the charge-migration dynamics. The component labeled $\Omega_1$ corresponds to the main frequency shown in figure 2. The curves labeled $2\Omega_1$ and $3\Omega_1$ mark FMA components with twice and three times the frequency value of the main component $\Omega_1$, respectively – note that below $\approx 3.5$ eV, the $3\Omega_1$ curve plots outside of the range of display. The labels ① and ② indicate the plateaus from which the sample CM-mode of figure 3 are taken. We also show the result of the FMA with the linearized-dynamics approximation of $\Delta\hat{H}_{\text{eff}} = 0$ in equation (3), which reduces the dynamics to the beating between molecular orbitals.</figcaption></figure>

## B. Kohn-Sham orbitals vs hole dynamics

The single KS-orbital density contributions shown in figure 3 (a,b) are important for understanding how the CM dynamics is organized, however, recall that they do

not correspond to physically observable quantities. For that, in figure 3 (c) we show the hole density associated with (b). Clearly, the two panels exhibit a similar pattern as the hole density mostly stays localized in space while it propagates through the $\pi$ system. Overall, the contributions from the other KS-orbital channels to the hole density show up mostly as higher-frequency patterns localized around the individual atomic centers. We understand the correspondence between the selected KS-orbital and the hole densities as the result of the orthonormality of the KS orbitals in equation (1): The localized density of panel (c) corresponds to an unpaired KS-orbital channel, which therefore contributes a single electron to the one-body density of equation (2). The other KS orbitals, which are forced to “stay away” from it via the orthogonality condition, all correspond to fully filled KS channels and therefore contribute two electrons to the density. In the end, the orthogonality between the unpaired and paired KS orbitals leads to the matching hole density in panel (c).

## C. Solitary-wave charge migration mechanism

To better understand the mechanism that regulates the particle-like hole dynamics shown in figure 3, we formally decompose the TDDFT Hamiltonian operator of equation (1) into its linear and nonlinear parts

$$ \hat{H}_{\text{eff}}[\rho(t)] = \hat{H}_{\text{eff}}[\rho_{\text{GS}}] + \Delta\hat{H}_{\text{eff}}[\rho(t)], \quad (3) $$

where $\rho_{\text{GS}}$ is the ground-state one-body density of the molecular cation. Modeling the CM as a beating of MOs amounts to neglecting the nonlinear part $\Delta\hat{H}_{\text{eff}}[\rho(t)]$ given that, in the basis of MOs, $\hat{H}_{\text{eff}}[\rho_{\text{GS}}]$ is a diagonal matrix with the MO energies on its diagonal. We find that the linearized MO-beating approximation yields qualitatively and quantitatively different results from the full TDDFT dynamics. To illustrate this, in figure 4 we compare the extended FMA for the full TDDFT dynamics, including time-dependent electron-electron interactions, with its linearized approximation $\Delta\hat{H}_{\text{eff}} = 0$ in equation (3), which reduces the dynamics to the beating between MOs. The figure shows a clear disagreement between the two results, both in terms of the shapes of the frequency components and in the number of them.

Intuitively, we understand the inadequacy of the MO-beating picture as follows: MOs are delocalized over the entire $\pi$ system and thus, in order to obtain a tightly localized electron density over a portion of the chain, one needs the coherent superposition of multiple MO wave functions. Then, the mismatch in the energy spacing between these MOs would lead to a decoherence, and thus spread, of the electronic density, which we do not observe in our simulations. This shows that the dispersion associated with the linear part of $\hat{H}_{\text{eff}}[\rho(t)]$ in equation (3) is cancelled by nonlinear effects associated with the nonlinear component $\Delta\hat{H}_{\text{eff}}[\rho(t)]$, ultimately leading to the non-dispersing solitary-wave dynamics shown in figure 3.
## D. Periodic charge-migration soliton modes

To conclude our nonlinear analysis of the model $\pi$ system, we return to figure 2. While the FMA of figure 2 sheds light on the CM dynamics unfolding from the initial hole we impose in the conjugated system, it only explores a narrow portion of the phase space that would otherwise be accessible to CM motions in general. Note also that one could design different ways to generate the initial hole than what we have chosen here, each potentially exploring a different portion of the phase space. To more fully explore the phase space in our model system, we complement our FMA with a direct systematic search for periodic CM modes that exhibit a traveling solitary wave similar to figure 3 (b,c), in the full parameter space of the model $\pi$ system. We discuss our strategy for finding these periodic CM modes in Appendix 5. We show the period and excitation-energy of these periodic CM modes with light grey markers in figure 2. As a whole, these results give a very clear picture of how the dynamics of particle-like CM is organized in the phase space of the $\pi$ system: The “ladder” of extended plateaus again shows that essentially the same CM period can be observed over a wide range of excitation energy, often spanning more than one eV. We also recover the general trend that slower CM modes are only available for more excited electronic configurations of the target – note the lack of periodic CM modes in the upper-left corner of the figure.

In figure 2, the comparison between the periodic CM modes and the FMA – light grey and dark orange markers, respectively – is stunning: All the plateaus in the latter match a set of periodic modes in the former. In other words, these periodic CM modes form the dynamical skeleton that regulates the CM we observed when creating an *ad hoc* localized hole at one end of the $\pi$ system. The “cleanliness” of that motion – e.g., figure 3 (a) vs (b) – depends on how close the initial condition puts the initial electronic configuration to a suitable periodic CM mode it can mimic. Altogether this suggests a two-pronged approach to CM studies: (i) asking whether the molecule of interest supports periodic soliton CM modes. If so, then (ii) identify how to tailor the ionization process to access the CM mode(s) of interest. We revisit this idea in section IV below.

## E. Functionalizing the carbon chain

The bare conjugated-carbon chain we have considered so far is, of course, highly symmetric and therefore an impractical system for generating a localized hole at one end in a realistic experimental scenario. Instead, one can consider functionalizing the chain by attaching a functional group at one of its ends, as sketched in figure 1. In simulations, we do so both for the model $\pi$ system and full TDDFT simulations (see section IV). Specifically, in this paper we consider the example of using a halogen

functional group.

Experiments and simulations have shown that strong field ionization can create localized ionization on the halogen center [32] in halocarbons and that the resulting valence hole can then migrate through the rest of the carbon-chain conjugated system [7, 22]. In the model $\pi$ system, we can emulate a halogen function by putting a single atomic center at one end of the conjugated-carbon chain. When varying the properties of this atomic center we have found that, in order for the electron/hole perturbation to migrate between the function and the carbon chain as it does in our original CM example of figure 1, their respective orbitals must hybridize when forming the overall compound’s electronic structure. This finding is consistent with observations in full TDDFT simulations of functionalized alkynes using different halogens [22]. Intuitively, such delocalized hybridized orbitals provide a bridge for the electron density to move between the function and carbon-chain parts of the molecule.

To conclude this section, we illustrate a second potential use of molecular functionalization, as a means to steer CM dynamics. In the same type of halo-functionalized model $\pi$ system as previsouly, we now start the hole on the chain side of the compound. Figure 5 compares the

<figure><img src="image_6.png" /></figure>

FIG. 5. Molecular functionalization can be used to steer charge migration (CM): Adding different functional groups to the same carbon chain can (top) inhibit or (bottom) enhance CM modes in the rest of the $\pi$ system. Here we vary the functionalization by having different distances between the halogen and the chain: 4 a.u. in panels (a,c) and 3 a.u. in (b,d). The figure compares CM dynamics after initiating the variably localized hole on the chain side of the molecule. Panels (a,b) correspond to extended FMAs similar to figure 4. For the initial condition marked by the vertical blue lines, panels (c,d) show two sample CM dynamics as in figure 3. In panel (b), like in figure 4, the labels $2\Omega_1$ and $3\Omega_1$ indicate FMA frequency components with twice and three times the fundamental frequency $\Omega_1$, respectively.
CM dynamics for two different functionalization configurations, here controlled with the distance between the function-atom and chain. Panel (a) shows a fuzzy FMA that lacks the clean solitary-wave CM-mode dynamics of the chain alone, as illustrated in panel (c). On the other hand, panel (b) exhibits clean frequency components with plateaus similar to that of the chain on its own. This is confirmed in the sample CM dynamics of panel (c) where the particle density travels almost periodically along the conjugated backbone and only transiently “leaks” to the function when the migration reaches the lower end of the chain. In both cases, the halogen function acts as a bias for the rest of the $\pi$ system by (bottom) enhancing or (top) inhibiting its CM modes.

#### IV. CHARGE MIGRATION IN “REAL” CONJUGATED MOLECULES

We now leverage the results of our nonlinear analysis of CM dynamics in the model $\pi$ system and apply them to “real” molecules and to experimental considerations. Though we cannot systematically explore the phase space in the full molecule in the way we did in the reduced system, the FMA analysis is quite revealing in terms of what we have already observed. We provide details of the TDDFT simulations, initial-hole configuration, and FMA computations in bromohexatriyne in Appendices 2, 3, and 4, respectively.

In figure 6 we show the result of an FMA analysis of CM in the singly-ionized bromohexatriyne molecular cation computed with the high-performance full-dimensional *ab initio* quantum-chemistry TDDFT package NWChem [40, 41]. Here we vary the initial portion of the electron hole localized on the Br atom using cDFT as discussed above and, after releasing the hole-localization constraint and following the time evolution, we compute the main frequency of CM dynamics

<figure><img src="image_7.png" /><figcaption>FIG. 6. A similar charge-migration frequency/period analysis as in figure 2 can be applied to full-dimension, all-electron, quantum-chemistry TDDFT simulations in bromohexatriyne (inset molecule). Orange dots mark the main frequency of the CM motion that moves the hole throughout the entire molecule – between the Br and final C≡C-H groups. We added the horizontal dashed grey lines by hand to highlight the two plateaus we observe in the frequency map.</figcaption></figure>

that move the density between the Br and final C≡C-H groups in the target – orange markers in figure 6. The results are strikingly similar to the simplified $\pi$-system model: As a trend, the period of the CM increases with the initial-hole localization (excitation energy). This increase happens in a step-wise fashion over ranges of initial conditions that lead to essentially the same CM period. Note the two distinctive plateaus between about 1.7 eV and 3.3 eV, and 4.3 eV and 6 eV, highlighted with the dashed grey lines in the figure. We have checked that the hole dynamics in these plateaus resembles that of the model $\pi$ system – see reference [22], and its supplemental movie, for sample CM motions in full simulations of halocarbons.

A notable quantitative difference between the FMA in the model $\pi$ system of figure 2 and in the full bromohexatriyne of figure 6 is that the latter has apparently fewer plateaus, and that the plateaus extend over a wider range of excitation energies. We speculate on two possible mechanisms for this. (i) *ionization condition*: By using an energy-minimization condition, the cDFT might be more “gentle” than the *ad hoc* initial-hole configuration we use in the model $\pi$ system. This would lead the FMA in bromohexatriyne to effectively explore a relatively narrower portion of phase space and thus “see” fewer of the CM-mode plateaus that the molecule can support. In an experimental context this speaks to the idea that, by properly tuning the initial parameters, one can have an ionization process that yields a relatively large bandwidth of excitation energy but still induces a consistent CM response. (ii) *structural condition*: Compared to the rest of the alkyne group, the Br function is effectively a heterogeneous center. So, while the alkyne and Br components properly hybridize to form the molecular cation, this mixing might not be compatible with all the CM-mode plateaus of the bare C chain alone. In turn this would result in a net reduction of the number of plateaus supported by the molecule and thus accessible to the FMA. This second explanation speaks to the idea of chemical control of CM motions, where a functional group added to a conjugated organic system can act as a bias and alter the ability of the rest of the molecule to support periodic CM modes – see the discussion around figure 5. More broadly, the possibility for altering a molecule’s ability to support CM dynamics through changes in its electronic configuration opens the door to doing so through external knobs like laser fields: The laser would selectively switch sustained CM “on” or “off” in the target sample by enabling or preventing its CM modes.

#### V. CONCLUSION AND OUTLOOK

In conclusion, we have performed detailed analyses of CM in conjugated organic molecules using tools from nonlinear dynamics (FMA and periodic-motion analysis). In the density picture we showed that periodic
CM modes, with a hole traveling back-and-forth through the $\pi$ system in a particle-like manner, emerge as solitary waves. This mechanism is fundamentally different from the few-orbital beating pictures that have previously been employed [6, 9, 18] and is driven by time-dependent multi-electron interactions. Obviously, for a detailed quantitative prediction of CM modes, like their precise period or other metrics associated with the CM dynamics [22], one needs a detailed modeling of the molecular system. However, we showed that a lower level of theory, using a simplified conjugated model, qualitatively reproduces key features of the full system. In an extension of the results shown here we have found similar features when further simplifying our model $\pi$ system and neglecting exchange and correlation interactions in the simulations (see Appendix 6). This suggests that the particle-like CM motions we identify emerge as a result of the dynamical mean-field interaction alone. Long-range many-electron interactions are a hallmark of conjugated organic molecules, which makes the possibility for sustained CM motions widely available in these systems – at least until the onset of nuclear motion.

Surprisingly, our analysis reveals that the same molecule can support several soliton CM modes with very different periods. For the model $\pi$ system shown in figure 2, the identified CM periods vary over a range of about 1 fs. In the full bromohexatriyne simulations of figure 6, the two identified CM modes are 500 as apart. Those modes are characterized by regions of the parameter space over which different initial conditions lead to the same almost-periodic CM motion, all essentially with the same period. In other words, the overall CM dynamics is not so much determined by the details of the individual electronic degrees of freedom – nor the competition between them – than it is driven by their collective response. This type of collective behavior is usually referred to as synchronization and has been identified throughout physics, engineering and biology [42–44]. We also note that these findings are compatible with the fully correlated Schrödinger-equation formalism, where field-free dynamics can only take discrete frequency values associated with the energy difference between the all-electron excited-state wave functions, and that more excited levels tend to be more closely spaced leading to longer periods. In comparison, our attosecond-soliton mechanism provides a novel way to understand CM dynamics directly in the time domain and real space. In all cases, the type of parametric robustness of the CM dynamics we have identified is essential for experimental applications as it provides a robustness of the migration dynamics against uncertainties in the way the hole might be created. The results presented in this work, along with the analysis tools that we employed, can help provide important perspectives for the design of future CM studies.

# ACKNOWLEDGMENTS

We thank L.F. DiMauro and R.R. Jones for helpful discussions on this work. This work was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Award No. DE-SC0012462.

# APPENDICES

## 1. Reduced $\pi$-system model

We systematically build our reduced $\pi$ system from one-dimensional (1D) carbon chains, denoted $(C_2)_n$, with $n$ the number of pairs of “C” centers. These chains are the 1D analog to 3D alkenes, without the hydrogen centers. We use soft-Coulomb [45] effective potentials to describe electron interactions, parameterized as $V_{sc}[Z, a](x) = -\frac{Z}{\sqrt{x^2+a^2}}$, with $Z$ the effective charge and $a$ the softening parameter. For electron-electron interactions we take a softening parameter $a_{ee} = \sqrt{2}$. In the model $\pi$ system, each C center contributes one electron with $Z_C = 1$. We set the respective positions for the atomic centers to emulate conjugation-like bonding by using slightly different inter-atomic distances within and between $C_2$ pairs, respectively set to 2.5 a.u. and 2.7 a.u. We chose these distances to be comparable to the ones from full-dimensional (functionalized) alkenes. Finally, for electron-ion interactions we select the softening parameter $a_C = 1$ for which we consistently find reasonable electronic properties throughout the $(C_2)_n$ family, both in terms of molecular orbital energies and shapes.

For simplicity, we only show results for the $(C_2)_4$ system. We systematically find similar results – multiple plateaus in the FMA associated with periodic CM modes with a traveling particle-like hole, longer CM periods accessible only to more excited configurations, etc. – for the other members of the $(C_2)_n$ family.

To functionalize the model $\pi$ system we add a single atomic “X” center with $Z_X = 2$ at one end of the molecule, meant to emulate a halogen group. We then use the softening parameter $a_X$ and distance between the function and chain to tune the coupling/hybridization between the two.

## 2. (TD)DFT simulations

For the 1D model $\pi$ system, we use spin-restricted TDDFT with local-density approximation (LDA) Slater-exchange and correlation potentials [46] and an average-density self-interaction correction (ADSIC) [47]. Numerically, we discretize the model $\pi$ system and simulate its TDDFT dynamics of Eq. (1) on a Cartesian grid.

For full-dimensional ab-initio TDDFT simulations in bromohexatriyne [22], we use an all-electron, spin polarized level of theory with the hybrid PBE0 functional,
cc-pVDZ basis set for the H and C atoms and Stuttgart RLC ECP for Br, as implemented in the NWChem package [40, 41].

### 3. Variable initial hole configuration

We initialize all simulations with a one-electron hole in the conjugated system of the molecular cation.

For the model $\pi$ system, we take an intuitive approach to building the initial condition using a single parameter to control the initial hole localization. In this system, the molecular cation has one unpaired KS orbital, while the other ones are fully filled. In the unpaired KS-orbital channel, we linearly mix a one-electron wave function localized around an end C dimer with the delocalized cation's ground state highest-occupied molecular orbital; the initial localization of the hole is controlled by this mixing coefficient. We then reorthonormalize the remaining paired KS orbitals. In the end, the difference in the number of electrons contributed by the unpaired and paired KS orbital yields the relative deficit of electronic density over one of the final C dimers.

Figure 7 illustrates how the initial electron and hole configuration depends on the localization parameter in the model $\pi$ system. The colormaps stop at 100% initial hole localization because our initial-condition mixing scheme has an upper bound for the excitation energy it can induce in the cation. On the other hand, we do not have the same upper limit when directly searching for periodic CM modes in the full parameter space, and we can therefore achieve the higher excitation energies shown in figure 2.

For full (3D) TDDFT simulations, the use of the constrained DFT (cDFT) method [35] was described in reference [22]. In reference [22] all CM simulations were initialized with the constraint of having exactly one electron hole on the halogen center. Here, instead, we keep the overall one-electron hole but vary the amount of that hole that is constrained to be localized around the Br center. Figure 8 illustrates how the initial hole configuration varies with the DFT constraint in bromohexatriyne computations. Note that the cDFT algorithm enables us to impose over 100% initial localization of the hole on the Br end.

### 4. Frequency-map analysis (FMA)

Generally speaking, an FMA tracks the variation of the main frequency component(s) associated with a dynamical process as a function of a continuously-varied initial condition [36, 37]. While we use slightly different implementation of the FMA for the analysis of the model $\pi$ system and full (3D) TDDFT simulations, the underlying idea is the same in both cases: To analyze CM motions across the entire molecule, we construct a complex-valued scalar signal $s_{\text{FMA}}(t) = s_l(t) + is_r(t)$ by

<figure><img src="image_8.png" /><figcaption>FIG. 7. Variation of the initial electron/hole configuration we use for the frequency-map analysis (FMA) of the model $\pi$ system. (a-c) Samples of the initial electron and hole density distributions along the molecular backbone for 20%, 55%, and 85% localization, respectively. (d) Initial electron density in the unpaired Kohn-Sham (KS) orbital we use as initial condition in the FMA of figure 2. (e) Initial hole density associated with (d). Note that panels (d,e) use the same colormaps as in figure 3.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. Like in figure 7 (e), variation of the initial hole configuration for the frequency-map analysis (FMA) of bromohexatriyne of figure 6.</figcaption></figure>

computing the amount of electron/hole density around the left and right ends of the molecule, respectively $s_{l,r}$. For the FMA itself, we use the time interval $10 \text{ fs} \le t \le 30 \text{ fs}$: We discard the first $10 \text{ fs}$ to avoid transient effects associated with the sudden introduction of the localized perturbation in the system. We have checked that including them in the analysis has only cosmetic effects
and does not change our results and conclusions.

For the model $\pi$ system (figure 2), we compute $s_{l,r}$ by projecting the electron density contribution for the unpaired KS orbital over the end $C_2$ dimers. We choose this projection for consistency with the results of figure 3 (a,b) where we showed that the underlying organization of the CM dynamics is best apparent in this KS-orbital channel. We have checked that we obtain essentially the same results when projecting the entire hole density on the same end dimers.

For the full TDDFT simulations (figure 6), we compute $s_{l,r}$ by simply integrating the hole density around the Br (left side of the median plane to the Br-C bond) and around the $-C\equiv C-H$ (right side of the median plane to the final $\equiv C-C\equiv$ bond) groups, respectively.

## 5. Periodic CM mode search

The TDDFT system of equations (1,2) has an infinite number of degrees of freedom, which makes it impracti-

cal for a direct computation of periodic CM modes. Instead, we employ a two-step approach: First, we restrict the dynamics to few occupied and unoccupied MOs of the corresponding cation and use a nonlinear solver – here the Levenberg-Marquardt Method as implemented in MATLAB® – to find periodic motions in that restricted space. Then, we use those restricted periodic modes as initial conditions in unrestricted TDDFT simulations and check that they still correspond to periodic motions. We stress that, with this approach, the restricted computations only serve as an intermediary to determining periodic CM modes in the same TDDFT framework we use for our other CM simulations.

## 6. $\pi$-system model without exchange-correlation interactions

Figure 9 illustrates how one obtains qualitatively very similar results when neglecting exchange-correlation interactions in the model $\pi$ system – compare to figure 4.

[1] F. Krausz and M. Ivanov, Attosecond physics, Rev. Mod. Phys. **81**, 163 (2009).

[2] F. Calegari, A. Trabattoni, A. Palacios, D. Ayuso, M. C. Castrovilli, J. B. Greenwood, P. Decleva, F. Martín, and M. Nisoli, Charge migration induced by attosecond pulses in bio-relevant molecules, J. Phys. B: At. Mol. Opt. Phys. **49**, 142001 (2016).

[3] J. Breidbach and L. S. Cederbaum, Universal attosecond response to the removal of an electron, Phys. Rev. Lett. **94**, 033901 (2005).

[4] S. Lünnemann, A. I. Kuleff, and L. S. Cederbaum, Ultrafast charge migration in 2-phenylethyl-n,n-dimethylamine, Chem. Phys. Lett. **450**, 232 (2008).

[5] O. Smirnova, Y. Mairesse, S. Patchkovskii, N. Dudovich, D. Villeneuve, P. Corkum, and M. Y. Ivanov, High harmonic interferometry of multi-electron dynamics in molecules, Nature **460**, 1476 (2009).

[6] F. Calegari, D. Ayuso, A. Trabattoni, L. Belshaw, S. De Camillis, S. Anumula, F. Frassetto, L. Poletto, A. Palacios, P. Decleva, J. Greenwood, F. Martín, and M. Nisoli, Ultrafast electron dynamics in phenylalanine initiated by attosecond pulses, Science **346**, 336 (2014).

[7] P. Kraus, B. Mignolet, D. Baykusheva, A. Rupenyan, L. Horný, E. Penka, G. Grassi, O. Tolstikhin, J. Schneider, F. Jensen, L. Madsen, A. Bandrauk, F. Remacle, and H. Wörner, Measurement and laser control of attosecond charge migration in ionized iodoacetylene, Science **350**, 790 (2015).

[8] A. I. Kuleff, N. V. Kryzhevoi, M. Pernpointner, and L. S. Cederbaum, Core ionization initiates subfemtosecond charge migration in the valence shell of molecules, Phys. Rev. Lett. **117**, 093002 (2016).

[9] F. Remacle and R. D. Levine, An electronic time scale in chemistry, Proceedings of the National Academy of Sciences **103**, 6793 (2006).

<figure><img src="image_10.png" /></figure>

FIG. 9. Same extended frequency-map analysis (FMA) as in figure 4, but here we have neglected the exchange-correlation term in our TDDFT simulation – i.e., electron-electron interactions are reduced to the mean-field Hartree potential only. The grey curves are the FMA components for the TDDFT simulations, labeled by the frequency values $\Omega_1$, its multiples $2\Omega_1$ and $3\Omega_1$, and $\Omega_2$. The green horizontal lines show the result of an FMA starting from the same initial conditions but using the linearized-dynamics approximation of $\Delta\hat{H}_{\text{eff}} = 0$ in equation (3), which reduces the dynamics to the beating between molecular orbitals.

[10] F. Lépine, M. Y. Ivanov, and M. J. J. Vrakking, Attosecond molecular dynamics: fact or fiction?, Nat. Photon. **8**, 1749 (2014).

[11] V. Despré, A. Marciniak, V. Loriot, M. C. E. Galbraith, A. Rouzée, M. J. J. Vrakking, F. Lépine, and A. I. Kuleff, Attosecond hole migration in benzene molecules surviving nuclear motion, J. Chem. Phys. Lett. **6**, 426 (2015).
[12] R. Weinkauf, E. W. Schlag, T. J. Martinez, and R. D. Levine, Nonstationary electronic states and site-selective reactivity, J. Phys. Chem. A **101** (1997).

[13] F. Remacle, R. Levine, and M. Ratner, Charge directed reactivity: a simple electronic model, exhibiting site selectivity, for the dissociation of ions, Chem. Phys. Lett. **285**, 25 (1998).

[14] M. Lara-Astiaso, M. Galli, A. Trabattoni, A. Palacios, D. Ayuso, F. Frassetto, L. Poletto, S. De Camillis, J. Greenwood, P. Decleva, I. Tavernelli, F. Calegari, M. Nisoli, and F. Martán, Attosecond pump-probe spectroscopy of charge dynamics in tryptophan, J. Phys. Chem. Lett. **9**, 4570 (2018).

[15] E. P. Månsson, S. Latini, F. Covito, V. Wanie, M. Galli, E. Perfetto, G. Stefanucci, H. Hübener, U. De Giovannini, M. C. Castrovilli, A. Trabattoni, F. Frassetto, L. Poletto, J. B. Greenwood, F. Légaré, M. Nisoli, A. Rubio, and F. Calegari, Real-time observation of a correlation-driven sub 3 fs charge migration in ionised adenine, Comm. Chem. **4**, 73 (2021).

[16] M. Vacher, M. J. Bearpark, M. A. Robb, and J. P. Malhado, Electron dynamics upon ionization of polyatomic molecules: Coupling to quantum nuclear motion and decoherence, Phys. Rev. Lett. **118**, 083001 (2017).

[17] V. Despré, N. V. Golubev, and A. I. Kuleff, Charge migration in propiolic acid: A full quantum dynamical study, Phys. Rev. Lett. **121**, 203002 (2018).

[18] D. Ayuso, A. Palacios, P. Decleva, and F. Martín, Ultrafast charge dynamics in glycine induced by attosecond pulses, Phys. Chem. Chem. Phys. **19**, 19767 (2017).

[19] A. Bruner, S. Hernandez, F. Mauger, P. M. Abanador, D. J. LaMaster, M. B. Gaarde, K. J. Schafer, and K. Lopata, Attosecond charge migration with TDDFT: Accurate dynamics from a well-defined initial state, J. Phys. Chem. Lett. **8**, 3991 (2017).

[20] D. Jia, J. Manz, B. Paulus, V. Pohl, J. C. Tremblay, and Y. Yang, Quantum control of electronic fluxes during adiabatic attosecond charge migration in degenerate superposition states of benzene, Chem. Phys. **482**, 146 (2017).

[21] E. Perfetto, D. Sangalli, A. Marini, and G. Stefanucci, Ultrafast charge migration in XUV photoexcited phenylalanine: A first-principles study based on real-time nonequilibrium Green's functions, J. Phys. Chem. Lett. **9**, 1353 (2018).

[22] A. S. Folorunso, A. Bruner, F. Mauger, K. A. Hamer, S. Hernandez, R. R. Jones, L. F. DiMauro, M. B. Gaarde, K. J. Schafer, and K. Lopata, Molecular modes of attosecond charge migration, Phys. Rev. Lett. **126**, 133002 (2021).

[23] J. Masoliver and A. Ros, Integrability and chaos: the classical uncertainty, Eur. J. Phys. **32**, 431 (2011).

[24] A. Lichtenberg and M. Lieberman, *Regular and Chaotic Dynamics*, Vol. 38 (Springer-Verlag, New York, 1992).

[25] W. H. Miller, Spiers memorial lecture quantum and semiclassical theory of chemical reaction rates, Faraday Discuss. **110**, 1 (1998).

[26] T. Bartsch, R. Hernandez, and T. Uzer, Transition state in a noisy environment, Phys. Rev. Lett. **95**, 058301 (2005).

[27] S. Kawai, A. D. Bandrauk, C. Jaffé, T. Bartsch, J. Palacián, and T. Uzer, Transition state theory for laser-driven reactions, J. Chem. Phys. **126**, 164306 (2007).

[28] G. van de Sand and J. M. Rost, Irregular orbits generate higher harmonics, Phys. Rev. Lett. **83**, 524 (1999).

[29] F. Mauger, C. Chandre, and T. Uzer, Strong field double ionization: The phase space perspective, Phys. Rev. Lett. **102**, 173002 (2009).

[30] F. Mauger, C. Chandre, and T. Uzer, Recollisions and correlated double ionization with circularly polarized light, Phys. Rev. Lett. **105**, 083002 (2010).

[31] W. Becker, X. Liu, P. J. Ho, and J. H. Eberly, Theories of photoelectron correlation in laser-driven multiple atomic ionization, Rev. Mod. Phys. **84**, 1011 (2012).

[32] P. Sándor, A. Sissay, F. Mauger, M. W. Gordon, T. T. Gorman, T. D. Scarborough, M. B. Gaarde, K. Lopata, K. J. Schafer, and R. R. Jones, Angle-dependent strong-field ionization of halomethanes, J. Chem. Phys. **151**, 194308 (2019).

[33] D. R. Tuthill, F. Mauger, T. D. Scarborough, R. R. Jones, M. B. Gaarde, K. Lopata, K. J. Schafer, and L. F. DiMauro, Multidimensional molecular high-harmonic spectroscopy: A road map for charge migration studies, J. Mol. Spectro. **372**, 111353 (2020).

[34] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. **140**, A1133 (1965).

[35] H. Eshuis and T. van Voorhis, The influence of initial conditions on charge transfer dynamics, Phys. Chem. Chem. Phys. **11**, 10293 (2009).

[36] J. Laskar, Frequency analysis for multi-dimensional systems. Global dynamics and diffusion, Physica D **67**, 257 (1993).

[37] J. Laskar, Introduction to frequency map analysis, In: Hamiltonian systems with three or more degrees of freedom, in *Hamiltonian Systems with Three or More Degrees of Freedom*, edited by C. Simó (Springer Netherlands, 1999) p. 134.

[38] J. Laskar and P. Robutel, The chaotic obliquity of the planets, Nature **361**, 608 (1993).

[39] J. Laskar, F. Joutel, and P. Robutel, Stabilization of the Earth's obliquity by the Moon, Nature **361**, 615 (1993).

[40] K. Lopata and N. Govind, Modeling fast electron dynamics with real-time time-dependent density functional theory: Application to small molecules and chromophores, J. Chem. Theory Comput. **7**, 1344 (2011).

[41] E. Aprà, E. J. Bylaska, W. A. de Jong, N. Govind, K. Kowalski, T. P. Straatsma, M. Valiev, H. J. J. van Dam, Y. Alexeev, J. Anchell, V. Anisimov, F. W. Aquino, R. Atta-Fynn, J. Autschbach, N. P. Bauman, J. C. Becca, D. E. Bernholdt, K. Bhaskaran-Nair, S. Bogatko, P. Borowski, J. Boschen, J. Brabec, A. Bruner, E. Cauët, Y. Chen, G. N. Chuev, C. J. Cramer, J. Daily, M. J. O. Deegan, T. H. Dunning, M. Dupuis, K. G. Dyall, G. I. Fann, S. A. Fischer, A. Fonari, H. Früchtl, L. Gagliardi, J. Garza, N. Gawande, S. Ghosh, K. Glaesemann, A. W. Götz, J. Hammond, V. Helms, E. D. Hermes, K. Hirao, S. Hirata, M. Jacquelin, L. Jensen, B. G. Johnson, H. Jónsson, R. A. Kendall, M. Klemm, R. Kobayashi, V. Konkov, S. Krishnamoorthy, M. Krishnan, Z. Lin, R. D. Lins, R. J. Littlefield, A. J. Logsdail, K. Lopata, W. Ma, A. V. Marenich, J. Martin del Campo, D. Mejia-Rodriguez, J. E. Moore, J. M. Mullin, T. Nakajima, D. R. Nascimento, J. A. Nichols, P. J. Nichols, J. Nieplocha, A. Otero-de-la Roza, B. Palmer, A. Panyala, T. Pirojsirikul, B. Peng, R. Peverati, J. Pitner, L. Pollack, R. M. Richard, P. Sadayappan, G. C.
Schatz, W. A. Shelton, D. W. Silverstein, D. M. A. Smith, T. A. Soares, D. Song, M. Swart, H. L. Taylor, G. S. Thomas, V. Tipparaju, D. G. Truhlar, K. Tsemekhman, T. Van Voorhis, A. Vázquez-Mayagoitia, P. Verma, O. Villa, A. Vishnu, K. D. Vogiatzis, D. Wang, J. H. Weare, M. J. Williamson, T. L. Windus, K. Woliński, A. T. Wong, Q. Wu, C. Yang, Q. Yu, M. Zacharias, Z. Zhang, Y. Zhao, and R. J. Harrison, NWChem: Past, present, and future, J. Chem. Phys. 152, 184102 (2020).

[42] J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Rítort, and R. Spigler, The Kuramoto model: A simple paradigm for synchronization phenomena, Rev. Mod. Phys. 77, 137 (2005).

[43] J. Feldhaus, J. Arthur, and J. Hastings, X-ray free-electron lasers, J. Phys. B: At. Mol. Opt. Phys. 38, S799 (2005).

[44] A. Pikovsky and M. Rosenblum, Synchronization, Schol-arpedia 2, 1459 (2007).

[45] J. Javanainen, J. H. Eberly, and Q. Su, Numerical sim-ulations of multiphoton ionization and above-threshold electron spectra, Phys. Rev. A 38, 3430 (1988).

[46] N. Helbig, J. I. Fuks, M. Casula, M. J. Verstraete, M. A. L. Marques, I. V. Tokatly, and A. Rubio, Density functional theory beyond the linear regime: Validating an adiabatic local density approximation, Phys. Rev. A 83, 032503 (2011).

[47] C. Legrand, E. Suraud, and P.-G. Reinhard, Comparison of self-interaction-corrections for metal clusters, J. Phys. B: At. Mol. Opt. Phys. 35, 1115 (2002).