# Ripple-like instability in the simulated gel phase of finite size
phosphocholine bilayers

Vivien Walter  Céline Ruscher  Adrien Gola  Carlos M. Marques  Olivier
Benzerara  Fabrice Thalmann  Department of Chemistry, King’s College London,
Britannia House, 7 Trinity Street, SE1 1DB, London, United Kingdom  Institut
Charles Sadron, CNRS and University of Strasbourg, 23 rue du Loess, F-67034
Strasbourg, France

###### Abstract

Atomistic molecular dynamics simulations have reached a degree of maturity
that makes it possible to investigate the lipid polymorphism of model bilayers
over a wide range of temperatures. However if both the fluid $L_{\alpha}$ and
tilted gel $L_{\beta^{\prime}}$ states are routinely obtained, the
$P_{\beta^{\prime}}$ ripple phase of phosphatidylcholine lipid bilayers is
still unsatifactorily described. Performing simulations of lipid bilayers made
of different numbers of DPPC (1,2-dipalmitoylphosphatidylcholine) molecules
ranging from 32 to 512, we demonstrate that the tilted gel phase
$L_{\beta^{\prime}}$ expected below the pre-transition cannot be obtained for
large systems ($>$ 94 DPPC molecules) through common simulations settings or
temperature treatments. Large systems are instead found in a disordered gel
phase which display configurations, topography and energies reminiscent from
the ripple phase $P_{\beta^{\prime}}$ observed between the pretransition and
the main melting transition. We show how the state of the bilayers below the
pretransition can be controlled and depends on thermal history and conditions
of preparations. A mechanism for the observed topographic instability is
suggested.

###### keywords:

Lipid bilayers , Molecular dynamics simulations , Phase transition

##  1 Introduction

Lipid membranes are fundamental components of living organisms, for their
pivotal role in the structure and the biochemistry of the cells. Their
properties are for a large part the consequence of the organisation of the
phospholipid molecules that compose them to a large extent. This organisation
is best revealed by experimenting with artificial lipid bilayers [1, 2].
Through a perfect control of the molecular composition, an extended range of
geometries and the possibility to insert membrane proteins and all sort of
molecules, artificial bilayers have become a standard tool in modern
biophysics [3, 4].

A striking property of pure phospholipids bilayers is to exhibit a number of
thermodynamic transitions upon temperature changes, meaning that they can be
found in several phases including a dense structured $L_{\beta^{\prime}}$
tilted gel phase at low temperature or a disordered fluid phase at higher
temperature $L_{\alpha}$ [3, 4, 5, 6]. Decades ago, a new phase have been
identified in the most common sort of phospholipid, the phosphocholines (PC
lipids) [7, 8, 3]. This new phase, specific to the PC lipids and called the
ripple phase $P_{\beta^{\prime}}$, is characterised by important corrugation
along the bilayer and the alternation between interdigitated leaflets and non-
interdigitated ones. The ripple phase is observed experimentally in pure
bilayers of 1,2-dipalmitoylphosphatidylcholine (DMPC) between 16 and 24∘C,
1,2-dipalmitoylphosphatidylcholine (DPPC) between 34 and 41∘C and
1,2-distearoylphosphatidylcholine (DSPC) between 51 and 55∘C [9, 3].

This structure is sometimes presented as resulting from an alternation of gel
and fluid lipid configurations [10]. There is no consensus regarding the
cause[11, 12, 13, 14, 15], and the experimental structure is still subject to
detailed investigations [16]. The interest in the ripple phase has grown over
the years, and intense research in its nature are made jointly in experimental
and numerical sciences [17, 18, 19].

In a recent paper, Khakbaz and Klauda investigated the phase transition of
DMPC and DPPC [20]. They reported formation of a structure resembling the
$P_{\beta^{\prime}}$ phase for DMPC bilayers in a range of temperature, while
for DPPC only a transition from the tilted to the fluid phase was observed in
bilayers composed of 70 lipids. More recently, we investigated the phase
transition of DPPC bilayers using Machine Learning algorithms [21]. In
membranes composed of 212 lipids, we observed a transition at 315 K from a
fluid phase to a condensed disordered gel phase similar to the ripple phase
that persists well below the pre-transition temperature of 307K and the
$L_{\beta^{\prime}}$ tilted phase was never obtained directly from a quench of
the $L_{\alpha}$ state.

In this article, we provide a detailed analysis of the nature of the phase of
DPPC simulated with the CHARMM36 force field below its melting temperature. We
investigate the state of DPPC bilayers at 288 K, for different system sizes
and different thermal routes. While the $L_{\beta^{\prime}}$ phase was
observed for small systems, we found that a disordered gel phase, reminiscent
from the $P_{\beta^{\prime}}$ phase, occurs in larger systems. In particular,
we noticed the formation of corrugations whose amplitude increases with the
system sizes investigated. By applying different thermal treatments, we
characterize the metastability of DPPC and show that disordered gel is the
preferred phase of simulated DPPC at 288K, in the thermodynamic limit.

##  2 Material & Methods

###  2.1 System Description

DPPC bilayers were obtained using the CHARMM-GUI online Membrane Builder [22,
23, 24, 25]. The size of the membrane was controlled via the number of lipids,
ranging from 32 to 512 and with equal amounts of lipids in both leaflets. The
bilayers were hydrated with water blocks of 10 nm on each side to prevent any
interaction of the leaflets through the PBCs along the Z-axis. The exact
composition of each system used is given in the Supplementary Informations.

The constructed system were minimized in energy using a steepest descent
algorithm and equilibrated by running two NVT simulations at 288 K for 10 ps,
with respectively a 0.001 ps and a 0.002 ps step; and two NPT simulations at
288 K and 1 bar, with a 0.002 ps step for respectively 100 ps and 1 ns.

###  2.2 Simulation Runs

Unless specified, the following conditions and parameters have been used for
all simulations. All simulations were performed using GROMACS 2016.4 [26, 27]
along with the CHARMM-36 all-atom force-field [28] (June 2015 version). The
force field parameters for the lipid molecules were provided directly by
CHARMM-GUI [29, 30].

All the molecular dynamics simulations used the leap-frog integration
algorithm [31] with a time step set to 2 fs. The temperature was controlled
during the simulation using a Nosé-Hoover thermostat [32, 33] with a
correlation time of $\tau_{T}=0.4$ ps, and a Parrinello-Rahman semi-isotropic
barostat [34, 35] set to 1 bar in all directions was applied to the system
(correlation time $\tau_{P}=2.0$ ps, compressibility $4.5\times 10^{-5}$
bar-1).

Lipid and water molecules were separately coupled to the thermostat. Following
GROMACS recommendations for the CHARMM-36 all-atom force field, a Verlet cut-
off scheme on grid cells was used with a distance of 1.2 nm, and non-bonded
interactions cut-offs (Van der Waals and Coulombic) were also set to 1.2 nm.
Fast smooth Particle-Mesh Ewald electrostatics was selected for handling the
Coulombic interactions, with a grid spacing of 4 nm. A standard cut-off scheme
with a force-switch smooth modifier at 1.0 nm was applied to the Van der Waals
interactions. We did not account for long range energy and pressure
corrections, and constrained all the hydrogen bonds of the system using the
LINCS algorithm.

Molecular dynamics production were run for a minimum of 50 ns. When
temperature treatments occurred, the systems were simulated for another 50 ns
to allow for the bilayer to reach equilibrium. Regardless of the length of the
simulation, the analysis were always performed on the last 25 ns of the
simulations.

###  2.3 Analysis

The areas per lipid of the systems were measured using two different methods:
(i) measure by projection of the bilayer in the XY plane of the box, noted
$A_{H}^{p}$, and (ii) measure by meshing of the water-lipid interfaces of each
leaflet of the bilayer, written $A_{H}^{m}$. The projected area per lipid was
measured directly in Gromacs by measuring the area of the simulation box in
the XY plane and by dividing it by the number of lipids per leaflet. The
meshed area per lipid was computed using Ovito 2.9. To do so, the surface of
the water blocks at the water-lipid interface was meshed using a probe sphere
radius of 6 and a smoothing level of 30 after removing the lipid molecules.
The area of the meshing was then divided by the total number of lipids in the
bilayer.

The amplitude and the period of the corrugations were measured using the
meshing of the water-lipid interface collected for the measurement of the area
per lipid. The position of the vertices of the meshing were interpolated using
SciPy [36] on a (200,200) uniform grid of the size of the simulation box, and
the height of these 200$\times$200 points were analysed. This grid was
directly used to generate the height maps. The amplitude of the corrugations,
represented by the root mean square (RMS) height of the points on the grid,
were calculated by subtracting the mean and then measuring the standard
deviation of the height $\sigma(h)$. The period of the corrugations were
calculated by computing the azimuthal power spectrum density (PSD) of the
points of the grid and determining the frequency of the highest peak of the
PSD.

To measure the local chain tilt direction, the positions of the atoms of the
tails of the lipids were extracted from the simulation files using MDAnalysis
[37, 38]. The vectors from the first carbon atom of each tail to all the
carbon atom of their respective tails were computed, before calculating the
mean vector. The vector field was then created by binning the membrane in
$n\times m$ squares and averaging over time the tilt direction vectors found
in each square $(X_{i},Y_{i})$. As a consequence, the vector field account for
the average direction of the tails found at this location on the membrane over
time and not to the average direction of a lipid over time.

The enthalpy of the systems were extracted from the simulation using directly
the tool provided with Gromacs, gmx energy. The systems used strictly had the
exact same number of atoms, for both lipids and water molecules, to prevent
changes due to the system composition. The (tilted) gel - fluid phase
transition enthalpy was measured using the method from [39]. Briefly, the
total enthalpy of a system simulated at different temperatures over a wide
range, here 283 to 358 K, was collected. The effects of the temperature
increase on the enthalpy, besides any phase transition, were removed from the
measurement by subtracting the affine baseline measured in each phase. The
values were then divided by the total number of lipids in the system, and the
transition enthalpy was calculated by integrating the difference between the
fitted gel and fluid baseline over the range of temperature at which the
transition occurs (here 308 to 318 K, accounting for Gromacs accuracy in
setting a temperature). For the tilted to disordered gel phase transition, the
enthalpy was measured by bringing the two systems at the exact same
temperature. The difference in enthalpy measured between these two systems was
then divided by the total number of lipids in the system.

##  3 Results

###  3.1 Corrugation formation and characterization

We considered DPPC bilayers prepared with CHARMM-GUI at temperature $T=288$K
whose size is ranging from 32 to 512 lipids. For small bilayers made of 64
lipids or less, we observed the formation of a smooth $L_{\beta^{\prime}}$
tilted gel phase, in agreement with the numerous experiments and simulations
reported in the literature [40]. The geometry of these bilayers remained
stable over time, even after a 50 ns simulation. When larger systems $(>94)$
lipids were concerned, we noticed the formation of a corrugation deforming
significantly the leaftlets (cf Figure 1).

Figure 1: DPPC bilayer made of 212 lipids simulated at 288 K right after
construction and equilibration. 2.5 nm thick slices of the system along each
axis are shown on the bottom, highlighting the corrugations of the membrane.
Top and bottom lipid leaflets are coloured respectively in cyan and orange
while the phosphorus atoms are shown as silver beads. Hydrogen atoms were
removed from the screenshot to improve readability. The PBC box is displayed
using the blue plain lines.

The corrugations takes place along both simulations box axis ($x$ and $y$). As
already observed and investigated in details for DMPC by Khakbaz and Klauda
[20], the molecular configurations of the lipids do not appear uniform along
the corrugations, and the typical stretched tails of the gel phase seem to
turn into the typical disordered tail configurations of the fluid phase in the
thin portions of the corrugations where interdigitation happens. These
variations in configuration can be highlighted by computing and mapping the
local segment order parameter $S_{\text{mol}}$ of each atom of each lipid (cf
Figure 2(a)).

|
---|---

(a)

(b)

Figure 2: (Left) Topography of the upper leaflet of (Top) 256 and (Bottom) 64
lipid bilayers, highlighting the corrugations observed in the membrane. Both
contour plot shares the same color code and contour line scale, in Å. The
white dashed line show the system before replication. (Right) Representation
of slices of the lipid membranes on the left, with a color code on the tail
atoms corresponding to their order parameters.

We systematically characterized the corrugation and probed its evolution with
system size. To quantify it in a more specific way, we investigated the
topographic elevation function for each leaflet (see Figure 2(b)) from which
we defined the amplitude of the corrugation as the root mean square (RMS)
height $h_{\text{RMS}}$ of the two water-lipid interfaces of the bilayers. The
heights of the interfaces were obtained by meshing the water surface and
removing the mean height of each leaflet. In these circumstances,
$h_{\text{RMS}}$ is equal to the standard deviation of the heights
$\sigma(h)$. The leaflet corrugation amplitudes are identical for both
leaflets even though cross section pictures might suggest otherwise (Figure
3(a)). The corrugation amplitude increases with system size but saturates for
large systems (Figure 3(b)). The longitudinal period of the corrugations
always coincide with the periodic boundary conditions (PBC, see Figure 3(c)).
The power spectrum densities (PSD) used to measure the periods are given in
the Supplementary Materials.

|  |
---

(a)

(b)

(c)

Figure 3: Effect of the size of the simulation box on the geometry of
bilayers. (a) 2.5 nm thick slices of the systems made with different amount of
lipid molecules and simulated at 288 K. Measurement of the (b) amplitude, or
$h_{\text{RMS}}$ and (c) period of the corrugations observed in different
systems versus the number of lipids. The evolution of the size of the box is
shown on top of the period. Dashed lines are respectively the cubic root and
square root fits of the data. Figure 4: (Top) Differences between the area per
lipid measured by surface meshing and by projection in the XY plane $\alpha$
(in %) and (Bottom) evolution of each type of area per lipid for different
system sizes. Blue dashed line on the bottom graph is the average projected
area per lipid.

To better quantify the nature of the corrugation, we investigated the area per
lipid using two methods: i) _projection_ in the $(XY)$ plane of the box, i.e
the area $A_{H}^{p}$ of the box divided by the number of lipid; ii) _meshing_
of the water lipid interface that enables the determination of the interfacial
area $A_{H}^{m}$. We found that the difference between interfacial and
projected area per lipid

$\displaystyle\alpha=\frac{A_{H}^{m}-A_{H}^{p}}{A_{H}^{p}}$ (1)

provides an excellent characterization of the corrugated phase (Figure 4).
Indeed, $\alpha$ is three times larger for corrugated systems ($\alpha=14\pm
2\%$) than for tilted ones ($\alpha=5\pm 2\%$). The differentiation using
usual methods such as the area per lipid or the tail order parameter achieved
less significant results (cf Supplementary materials).

The local chain tilt direction in the XY plane was also determined and mapped
(Figure 5(a)). As expected, tilts are uniform in the $L_{\beta^{\prime}}$
state and present a random short-range correlation in the fluid $L_{\alpha}$
state. This was confirmed by the distribution of the angles with respect to
the X-axis, as a function of the height on the leaflet. In the corrugated
systems, the local chain tilts displays long-range variations. The angle
distribution of the corrugated systems shows that small heights have random
orientation similar to the fluid, while small heights are subject to less
variations.

|
---|---

(a)

(b)

Figure 5: (a) Orientation in the XY plane of the tails of each of the 256
lipids forming a DPPC bilayer (Top) in the disordered gel phase, (Middle) in
the tilted gel phase and (Bottom) in the fluid phase. Only one leaflet of the
membrane is shown for each system. (b) Respective scatter distribution of the
angles between the tails of the lipids and the X-axis of the membrane, as a
function of the height on the leaflet. The measurement were performed on both
top and bottom leaflets, respectively colored in cyan and in orange.

We also probed the sensitivity of the corrugation to modifications of space
group. A simple change from cubic to rectangular simulation boxes with
$L_{x}=2L_{y}$ shows commensurate corrugations periods. Moreover we carry out
simulations in hexagonal and monoclinic simulation boxes counting 128 lipids.
Both of these systems were found to be corrugated, with modulation vectors
directed along the PBC/crystallographic directions (Figure 6).

---

(a)

(b)

Figure 6: Topography of bilayers (a) made of 188 DPPC molecules in a
simulation box with $L_{x}=2L_{y}$ and (b) made of 128 DPPC molecules in an
hexagonal simulation box, both simulated at 288 K.

To release residual stress on the bilayer configuration that could have been
brought by the semi-isotropic barostat, an anisotropic barostat was applied to
the systems, with 1 bar and a $4.5\times 10^{-5}$ bar-1 compressibility set
along each axes, and the pressure crossed terms set to 0 for both pressure and
compressibility. We found that releasing residual anisotropic stress did not
modify the structure of the bilayers. After this simulation run, the system
was still found in a corrugated state, as shown in Figure 7, with a relative
area increase $\alpha$ of 11.9 $\pm$ 0.8% and a RMS corrugation amplitude of
4.89 Å comparable to the amplitude before the run (5.18 Å). We therefore
conclude from this result that both tilted gel and corrugated states behave as
a cohesive, solid state on the simulated time scales. They also display
significant residual static stresses of 0.8 and 1.0 bars respectively in the x
and the y directions. Based on these characterizations we refer to the
corrugated state as the ”disordered gel state”, noted here $L^{d}_{\beta}$, to
make a clear distinction between these configurations and the ripple phase
$P_{\beta^{\prime}}$ that has only be experimentally reported above the pre-
transition temperature.

Finally one can wonder whether the appearance of the corrugation is restricted
to DPPC or DMPC. To answer this question, we probed the effect of the tail and
the head groups by considering the longer-tailed DSPC and the ethanolamine-
based DPPE. As shown in Figure 8, the large DSPC systems were found in the
disordered gel state while the DPPE systems remained in the expected
homogeneous tilted gel state. The respective area differences of these systems
are 12 $\pm$ 1 and 5 $\pm$ 1 %.

Figure 7: Slice of the 212 DPPC molecule system obtained after being simulated
at 288 K for 50 ns with an anisotropic barostat. Figure 8: Profiles of
bilayers made of 64 or 212 lipids, either with DPPC, DSPC or DPPE, after a 50
ns simulation at 288 K. Only the DPPE membrane remained in the tilted gel
phase when constructed with a large number of lipids.

###  3.2 Influence of thermal history

The previous observations suggest that disordered $L^{d}_{\beta}$ and tilted
$L_{\beta^{\prime}}$ gel states are two competing states whose appearance seem
to be correlated to the system size for freshly thermalised systems generated
with CHARMM-GUI. In what follows we probe the sensitivity of the stability of
these phases to different routes of thermal treatments.

We focus our analysis on two system sizes: 64 lipids, which has been found in
the $L_{\beta^{\prime}}$ phase, and 256 lipids which has been found in the
$L^{d}_{\beta}$. Both systems were subjected to the following thermal
treatment: starting from $T=288K$, systems were annealed at $T=358$K in the
fluid phase. These systems are respectively named pc64-A and pc256-A, for
annealing. They were then cooled down to $T=288$K in two different ways:
either with a brutal fast cooling, named here quenched, or with a slow gradual
cooling of $1$K/ns that we denote gentle cooling. The quenched systems are
noted pc64-AQ and pc256-AQ (annealing-quenching), and the cooled systems
called pc64-AC and pc256-AC (annealing-cooling).

We first notice that whatever the thermal history, large systems were always
found in the disordered gel state as shown in Figure 9. Gentle cooling of
small systems allowed them to recover the $L_{\beta^{\prime}}$ phase while
quenching lead to $L^{d}_{\beta}$ phase. As in this latter case, either phases
could be obtained depending on the thermal history, we conclude that the
disordered gel is metastable with respect to the tilted gel for small systems.
Furthermore, we can note that both the gently cooled and the quenched systems
have a final projected area per lipid $A_{H}^{p}$ (respectively 50.1 $\pm$0.7
and 50.8 $\pm$ 0.8 Å2) close to the average projected area of 50 Å2 found for
systems of all sizes, hence correcting the odd value of 52.3 Å2 found before
temperature treatment (cf Figure 4).

As large systems of 256 lipids were never spontaneously found in the
$L_{\beta^{\prime}}$ state, we decided to force them into this state by
duplicating along both $x$ and $y$ directions the 64 $L_{\beta^{\prime}}$
system. The flat tilted duplicated system was found to be stable on the
simulation time scale, suggesting that $L_{\beta^{\prime}}$ could be either
stable or metastable in large systems too. The nature of the relevant stable
thermodynamic phase for large systems remains an open question, while our
simulations clearly favor the disordered gel state $L^{d}_{\beta}$.
Measurements of the difference in area $\alpha$ can be found in the
Supplementary Materials.

Figure 9: Synopsis of the different thermal trajectories, showing the systems
obtained for different cooling rates as well as the systems they originated
from. Only the results from the 256 DPPC system obtained by replication of the
64 lipid system are shown for the quenching and cooling experiments. Results
for the 256 DPPC system constructed by CHARMM-GUI, as well as the intermediate
systems, are shown in the Supplementary Materials.

###  3.3 Thermodynamics of tilted and disordered states

Since the small systems made of 64 lipids can be prepared and controlled to
reach all the observed phases, we used them to compare the energetic
properties of the respective gel phase. Having in mind the idea of a complex
underlying potential energy landscape composed of several minima located at
different energy levels, we probe how far energetically the disordered gel
phase stand from the tilted gel phase minimum. Therefore we performed an
energy minimization using conjugate gradient to remove thermal fluctuations on
the tilted and disordered gel configurations. We found $\Delta
E_{P}=E_{P}(\text{disordered})-E_{P}(\text{tilted})=64$kJ/mol, meaning that if
both states are metastable, the tilted system is the most stable one. In
addition, we also computed the difference in enthalpy between the disordered
and the tilted thermalised gel phases. To this aim, we used the same initial
states, after quenching or gentle cooling but without energy minimisation.
Since the two systems share the same atomic compositions (same number of
lipid, water molecules and constraints), the difference in enthalpy (kinetic
and potential energy according to the force field, plus the $PV$ contribution)
at a given temperature should be characteristic of the enthalpy difference
between the two states. We found an enthalpy difference of $12\pm 4$kJ/mol at
$288$K. This value is higher than the one reported experimentally for the
calorimetric gel to ripple phase pretransition, 4.6 kJ/mol [4] but obtained at
a higher temperature. We therefore repeated the measure after rising both
systems at a temperature close to the pretransition temperature, namely 305 K,
and then found a difference in enthalpy of $3.9\pm 0.4$ kJ/mol now compatible
with experimental measurements. It is essential to note here that, in order to
obtain a disordered system configuration at 305 K, the system had first to be
annealed at 375 K instead of 358 K, meaning that the difference in temperature
required for quenching should be at least of 70 K in order to obtain the
$L^{d}_{\beta}$ phase (see Supplementary materials).

The enthalpy of the $L^{d}_{\beta}$ $\longrightarrow L_{\alpha}$ transition
was also measured in our small systems. This was performed by simulating the
64 DPPC system at temperatures ranging from 283 to 358 K and by removing the
change in enthalpy due to the change in temperature (data shown in
Supplementary Materials). The transition enthalpy was found equal to 27.3
kJ/mol, which is also comparable with the experimental reported values circa
32.2 kJ/mol [4]. We can conclude that despite the approximations, and the
absence of quantum corrections to the bond vibrations contributions, the
CHARMM36 force field thermodynamic predictions seem in good agreement with
experimental observations.

##  4 Discussion

A careful inspection of the corrugations shows that the topographic modulation
is imposed by the periodicity of the simulation box. The corrugation shows
finite system size dependence, as its amplitude increases up to sizes of the
order of $L_{x}=8$ nm where the modulation saturates to a value of 7 $\pm$ 1
Å. Assuming a sine-like corrugation, this maximum RMS amplitude measured can
be converted into a peak-to-peak amplitude of 2.0 $\pm$ 0.2 nm similar to the
previously published values of 2.4 nm for the DPPC or even 1.8 nm in DMPC
bilayers [41, 16]. The saturation of the corrugation amplitude is expected,
given that it can only reach a fraction of the total membrane thickness. The
associated lateral length 8 nm can therefore be considered as a lower bound of
the instability characteristic longitudinal length scale.

We believe that the relative difference between the interfacial area and the
projected one, $\alpha$, introduced in equation (1) can be taken as a relevant
order parameter for the transition between the $L_{\beta^{\prime}}$ and the
disordered gel phases $L^{d}_{\beta}$. The latter being reminiscent from
$P_{\beta^{\prime}}$ ripple phase, $\alpha$ could be seen as a critical
parameter to investigate the existence of the ripple phase and discriminate it
from the tilted gel or from the fluid phase.

However, unlike experiments, the numerical instability occurs along two
orthogonal directions, or along the hexagonal axes. Non-square boxes fail to
select only one modulation direction. We nevertheless think that the numerical
corrugation instability is related to the experimental ripple instability, as
also suggested by the dependence of the presence of corrugation to the
chemical nature of the heads and tails of the lipids. Indeed, experiments have
shown that the ripple phase is specific from the phosphocholine (PC) lipids
[40, 42].

Another striking observation is the insensitivity of large systems (256
lipids) to thermal treatment which have been systematically ended in the
disordered gel state. On the opposite small systems can alternate between both
phases. However, when thermalised at low temperatures from CHARMM-GUI or
slowly quenched, they end up preferentially into the tilted gel state,
suggesting that in the range of temperatures investigated the
$L_{\beta^{\prime}}$ state is thermodynamically favored. By contrast the
$L_{\beta^{\prime}}$ state is never selected spontaneously by larger systems.
For the latter, MD suggests that the disordered gel state $L^{d}_{\beta}$ is
the most preferred state for all temperatures below melting. Metastability and
kinetic effects are certainly significant and may hide to true nature of the
stable phase.

Assuming that the difference between tilted and disordered gels has something
to do with the pretransition, we found an enthalpy difference 3 times larger
than the premelting latent heat at 288K, and of the same order of magnitude at
305K, close to the observed experimental transition. This points towards the
relevance of disordered gel state $L_{\beta}^{d}$ as a ripple state analogue.
Moreover the higher enthalpy of the disordered gel with respect to the tilted
gel is consistent with a $L_{\beta^{\prime}}\longrightarrow$ $L^{d}_{\beta}$
$\longrightarrow L_{\alpha}$ sequence of transitions as temperature is
increased.

The reason of the outcome of a ripple instability below the melting
temperature is nothing but obvious. Our simulations point out to a competition
between an homogeneous tilted state and an inhomogeneous corrugated state. The
transition between these states is discontinuous. The corrugated state is not
very tilted, and partially melted, or disordered, and interdigitated. We
suggest now a possible mechanism explaining the observed situation. The tilted
phase can be understood as the result of a frustration between lipid
headgroups which try to increase their exposure to water in the interface
region, lipid chains which try to reach an optimal packing density as a result
of cohesive forces, and chain stiffness for which the introduction of gauche
dihedral angles is unfavorable at low temperatures. Tilt allows lipid to
optimize simultaneously those three constraints. On the other hand, melting a
lipid chain enables the release of the constraint acting on the chain
stiffness, and makes it possible to increase the hydration free-energy by
reducing the membrane thickness. Below melting temperature, thermodynamics
makes it unfavorable to melt all the lipid molecules. However, some local
disordering of the lipids may still be favorable, increasing the hydration of
the headgroups without need of spending too much energy in melting the chains.

Based upon those considerations, we designed a simple one dimensional lipid
chain model that supports the idea that in an temperature range just below
melting, the homogeneous tilted state energetically unstable with respect to a
local corrugation of the bilayer, see Figure 10. Details of the
parameterisation of the model and further results are presented in the
Supplementary Information section. We conclude that the thickness modulation
might indeed be caused by a subtle interplay between headgroup hydration,
hydrophobic chain packing, trans-gauche isomerisation and tilt elasticity
energy terms.

Figure 10: Schematic representation of the mechanical instability occurring in
a simple one dimensional lipid chain model, where the tilted state has a
higher energy than a locally disordered state. The gain in energy originates
mostly from the hydration term. Blue solid rods stand for gel state, red
dashed rods for fluid state.

##  5 Conclusion

We have successfully demonstrated in this work how the size of the simulation
box influences the ripple-like instability in a PC membrane simulated with the
Charmm36 force field at low temperature, where it is usually expected to be in
the tilted gel $L_{\beta^{\prime}}$ phase. This unexpected organisation, which
we called the disordered gel phase $L^{d}_{\beta}$, does not appear in small
systems, which is consistent with the results from Khakbaz and Klauda [20] as
well as with our results in a previously published paper [21]. The energy and
geometry analyses demonstrated that this disordered gel phase has a lot in
common with the $P_{\beta^{\prime}}$ phase. Furthermore, this instability was
not observed with PE lipids, in agreement with experimental findings. For
small systems, we found ways of preparing the system in either tilted gel or
disordered gel states by acting on the thermal treatment. More work is needed
to determine whether the Charmm36 force-field can describe a one-dimensional
$P_{\beta^{\prime}}$ spatial thickness modulation, with the right periodicity,
and at which temperatures. Our work suggests that simulations will have to be
guided to the desired structure. Finally, we conclude that the ripple
instability looks like a generic mechanism adopted by the phosphatidylcholines
lipids to increase the headgroup hydration while still satisfying the packing
constraints, at the expense of a mild cost in disordering/melting a small
fraction of the chains.

## Acknowledgements

V. W. warmly thanks T.E. de Oliveira for helping with simulations set-up. A.G.
acknowledges partial support from the Investissements d’Avenir program
“Développement de l’Économie Numérique” through the SMICE project. The authors
gratefully acknowledge support from the high performance cluster (HPC)
Equip@Meso from the University of Strasbourg, through grants no. G2019A131C
and G2020A140C.

## References

  * Mouritsen [2005] O. G. Mouritsen, Life-as a matter of fat: the emerging science of lipidomics, Springer, 2005\.
  * Dimova and Marques [2019] R. Dimova, C. M. Marques (Eds.), The Giant Vesicle Book, 1st ed., CRC Press, Taylor and Francis, 2019.
  * Marsh [2013] D. Marsh, Handbook of Lipid Bilayers, 2nd ed., CRC Press, Boca Raton, 2013\.
  * Cevc and Marsh [1987] G. Cevc, D. Marsh, Phospholipid Bilayers. Physical Principles and Models, John Wiley & Sons, New-York, 1987.
  * Mabrey and Sturtevant [1976] S. Mabrey, J. M. Sturtevant,  Investigation of phase transitions in lipids and lipid mixtures by high sensitivity differential scanning calorimetry,  Proceedings of the Natural Academy of Sciences USA 73 (1976) 3862–3866.
  * Heimburg [2007] T. Heimburg, Thermal Biophysics of Membranes, Wiley-VCH, 2007.
  * Tardieu et al. [1973] A. Tardieu, V. Luzzati, F. C. Reman,  Structure and polymorphism of the hydrocarbon chains of lipids: A study of lecithin-water phases,  Journal of Molecular Biology 75 (1973) 719–733.
  * Sun et al. [1996] W. J. Sun, S. Tristram-Nagle, R. M. Suter, J. F. Nagle,  Structure of the ripple phase in lecithin bilayers.,  Proceedings of the National Academy of Sciences 93 (1996) 7008–7012. doi:10.1073/pnas.93.14.7008.
  * Lewis et al. [1987] R. N. A. H. Lewis, N. Mak, R. N. McElhaney,  A differential scanning calorimetric study of the thermotropic phase behavior of model membranes composed of phosphatidylcholines containing linear saturated fatty acyl chains,  Biochemistry 26 (1987) 6118–6126. doi:10.1021/bi00393a026.
  * Heimburg [2000] T. Heimburg,  A model for the lipid pretransition: coupling of ripple formation with the chain-melting transition,  Biophysical Journal 78 (2000) 1154–1165.
  * Doniach [1979] S. Doniach,  A thermodynamic model for the monoclinic (ripple) phase of hydrated phospholipid bilayers,  The Journal of Chemical Physics 70 (1979) 4587–4596. doi:10.1063/1.437292.
  * Lubensky and MacKintosh [1993] T. C. Lubensky, F. C. MacKintosh,  Theory of “ripple” phases of lipid bilayers,  Physical Review Letters 71 (1993) 1565–1568. doi:10.1103/physrevlett.71.1565.
  * Fournier [1998] J.-B. Fournier,  Coupling between membrane tilt-difference and dilation: A new “ripple” instability and multiple crystalline inclusions phases,  Europhysics Letters (EPL) 43 (1998) 725–730. doi:10.1209/epl/i1998-00424-4.
  * Misbah et al. [1998] C. Misbah, J. Duplat, B. Houchmandzadeh,  Transition to ripple phases in hydrated amphiphiles,  Physical Review Letters 80 (1998) 4598–4601. doi:10.1103/physrevlett.80.4598.
  * Sengupta et al. [2001] K. Sengupta, V. A. Raghunathan, Y. Hatwalne,  Role of tilt order in the asymmetric ripple phase of phospholipid bilayers,  Physical Review Letters 87 (2001). doi:10.1103/physrevlett.87.055705.
  * Akabori and Nagle [2015] K. Akabori, J. F. Nagle,  Structure of the dmpc lipid bilayer ripple phase,  Soft Matter 11 (2015) 918–926.
  * de Vries et al. [2005] A. H. de Vries, S. Yefimov, A. E. Mark, S. J. Marrink,  Molecular structure of the lecithin ripple phase,  Proceedings of the National Academy of Sciences 102 (2005) 5392–5396. doi:10.1073/pnas.0408249102.
  * Lenz and Schmid [2007] O. Lenz, F. Schmid,  Structure of symmetric and asymmetric ”ripple” phases in lipid bilayers,  Physical Review Letters 98 (2007).
  * Debnath et al. [2014] A. Debnath, F. M. Thakkar, V. Kumaran, P. K. Maiti, K. G. Ayappa,  Laterally structured ripple and square phases with one and two dimensional thickness modulations in a model bilayer system,  Soft Matter 10 (2014) 7630–7637.
  * Khakbaz and Klauda [2018] P. Khakbaz, J. B. Klauda,  Investigation of phase transitions of saturated phosphocholine lipid bilayers via molecular dynamics simulations,  Biochimica et Biophysica Acta (BBA) - Biomembranes 1860 (2018) 1489 – 1501.
  * Walter et al. [2020] V. Walter, C. Ruscher, C. M. Marques, O. Benzerara, F. Thalmann,  A machine learning study of the two states model for lipid bilayer phase transitions,  Phys. Chem. Chem. Phys. (2020) –.
  * Jo et al. [2007] S. Jo, T. Kim, W. Im,  Automated builder and database of protein/membrane complexes for molecular dynamics simulations,  PLoS ONE 2 (2007) 880.
  * Jo et al. [2008] S. Jo, T. Kim, V. G. Iyer, W. Im,  Charmm-gui: A web-based graphical user interface for charmm,  Journal of Computational Chemistry 29 (2008) 1859–1865.
  * Jo et al. [2009] S. Jo, J. B. Lim, J. B. Klauda, W. Im,  Charmm-gui membrane builder for mixed bilayers and its application to yeast membranes,  Biophysical Journal 97 (2009) 50–58.
  * Wu et al. [2014] E. L. Wu, X. Cheng, S. Jo, H. Rui, H. K. Song, E. M. Davila-Contreras, Y. Qi, J. Lee, V. Monje-Galvan, R. M. Venable, J. B. Klauda, W. Im,  Charmm-gui membrane builder toward realistic biological membrane simulations,  Journal of Chemical Theory and Computation 35 (2014) 1997–2004.
  * Berendsen et al. [1995] H. Berendsen, D. van der Spoel, R. van Drunen,  Gromacs: A message-passing parallel molecular dynamics implementation,  Computer Physics Communications (1995).
  * Abraham et al. [2015] M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, E. Lindahl,  Gromacs: High performance molecular simulations through multi-level parallelism from laptops to supercomputers,  SoftwareX (2015).
  * Best et al. [2012] R. B. Best, X. Zhu, J. Shim, P. E. M. Lopes, J. Mittal, M. Feig, A. D. MacKerell Jr.,  Optimization of the additive charmm all-atom protein force field targeting improved sampling of the backbone phi, psi and side-chain khi1 and khi2 dihedral angles,  Journal of Chemical Theory and Computation 8 (2012) 3257–3273.
  * Brooks et al. [2009] B. R. Brooks, C. L. Brooks III, A. D. MacKerell Jr, L. Nilsson, R. J. Petrella, B. Roux, Y. Won, G. Archontis, C. Bartels, S. Boresch, A. Caflisch, L. Caves, Q. Cui, A. R. Dinner, M. Feig, S. Fischer, J. Gao, M. Hodoscek, W. Im, K. Kuczera, T. Lazaridis, J. Ma, V. Ovchinnikov, E. Paci, R. W. Pastor, C. B. Post, J. Z. Pu, M. Schaefer, B. Tidor, R. M. Venable, H. L. Woodcock, X. Wu, W. Yang, D. M. York, M. Karplus,  Charmm: The biomolecular simulation program,  Journal of Computational Chemistry 30 (2009) 1545–1614.
  * Lee et al. [2016] J. Lee, X. Cheng, J. M. Swails, M. S. Yeom, P. K. Eastman, J. A. Lemkul, S. Wei, J. Buckner, J. C. Jeong, Y. Qi, S. Jo, V. S. Pande, D. A. Case, C. L. Brooks III, A. D. MacKerell Jr, J. B. Klauda, W. Im,  Charmm-gui input generator for namd, gromacs, amber, openmm, and charmm/openmm simulations using the charmm36 additive force field,  Journal of Chemical Theory and Computation 12 (2016) 405–413.
  * Hockney et al. [1974] R. W. Hockney, S. P. Goel, J. W. Eastwood,  Quiet high-resolution computer models of a plasma,  Journal of Computational Physics 14 (1974) 148–158.
  * Nosé [1984] S. Nosé,  A molecular dynamics method for simulations in the canonical ensemble,  Molecular Physics 52 (1984) 255–268.
  * Hoover [1985] W. G. Hoover,  Canonical dynamics: Equilibrium phase-space distributions,  Physical Review A 31 (1985) 1695.
  * Nose and Klein [1983] S. Nose, M. L. Klein,  Constant pressure molecular dynamics for molecular systems,  Molecular Physics 50 (1983) 1055–1076.
  * Parinello and Rahman [1998] M. Parinello, A. Rahman,  Polymorphic transitions in single crystals: A new molecular dynamics method,  Journal of Applied Physics 52 (1998) 7182.
  * Virtanen et al. [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. Jarrod Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, S. . . Contributors,  Scipy 1.0: Fundamental algorithms for scientific computing in python,  Nature Methods 17 (2020) 261–272.
  * Michaud‐Agrawal et al. [2011] N. Michaud‐Agrawal, E. J. Denning, T. B. Woolf, O. Beckstein,  Mdanalysis: A toolkit for the analysis of molecular dynamics simulations,  Journal of Computational Chemistry 32 (2011) 2319–2327.
  * Gowers et al. [2016] R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, O. Beckstein,  Mdanalysis: A python package for the rapid analysis of molecular dynamics simulations,  Proceedings of the 15th Python in Science Conference (2016) 98–105.
  * de Oliveira et al. [2020] T. E. de Oliveira, F. Leonforte, L. Nicolas-Morgantini, A.-L. Fameau, B. Querleux, F. Thalmann, C. M. Marques,  Fluid bilayer phase in aqueous mixtures of fatty alcohol and cationic surfactant,  Phys. Rev. Research 2 (2020) 013075. doi:10.1103/PhysRevResearch.2.013075.
  * Cunningham et al. [1998] B. A. Cunningham, A.-D. Brown, D. H. Wolfe, W. P. Williams, A. Brain,  Ripple phase formation in phosphatidylcholine: Effect of acyl chain relative length, position, and unsaturation,  Physical Review E 58 (1998).
  * de Vries et al. [2005] A. H. de Vries, S. Yefimov, A. E. Mark, S. J. Marrink,  Molecular structure of the lecithin ripple phase,  PNAS 102 (2005) 5392–5396.
  * Katsaras et al. [2000] J. Katsaras, S. Tristram-Nagle, Y. Liu, R. L. Headrick, E. Fontes, P. C. Mason, J. F. Nagle,  Clarification of the ripple phase of lecithin bilayers using fully hydrated, aligned samples,  Physical Review E 61 (2000).


2021

[2]Jean Christophe Tremblay

1]Department of Physics, Indian Institute of Technology Bombay, Powai, Mumbai,
400076, India

2]CNRS-Université de Lorraine, UMR 7019, , 1 Bd Arago, Metz, 57070, France

# Charge Migration in Heterocyclic Five-Membered Rings

Sucharita Giri  Gopal Dixit<EMAIL_ADDRESS><EMAIL_ADDRESS>[  [

###### Abstract

This contribution presents numerical simulations of $N$-electron dynamics in
heterocyclic five-membered ring molecules to shed light on the effect of
molecular symmetry on charge migration. Laser-driven dynamics is studied using
the hybrid time-dependent density functional theory/configuration methodology,
and the ensuing field-free charge migration is investigated by means of
transient electronic flux density maps. Our results demonstrate that the
charge migration in aromatic rings is sensitive to the presence of heteroatoms
such as oxygen and nitrogen. Their presence within the ring induces
significant modifications of the character in the ground and low-lying
electronic states, which is imprinted in the charge migration mechanism.

###### keywords:

attosecond, charge migration, symmetry reduction, electronic flux density

##  1 Introduction

Recent developments in laser technologies facilitate the generation of
ultrashort and intense laser pulses. The availability of such laser pulses
allows studying ultrafast processes in molecules on their natural timescale
krausz2009attosecond . During laser-molecule interaction, electrons are the
first to respond to the action of an ultrashort laser pulse. As a result of
the interaction, a superposition of several electronic states in the form of
an electronic wave packet is created, which triggers an electronic charge
distribution into motion. The migration of electronic charge around the atoms
of the molecular scaffold affects the outcomes of chemical reactions and
biological functions cederbaum1999ultrafast ; kuleff2005multielectron ;
breidbach2003migration ; remacle2006electronic ; breidbach2005universal ;
nisoli2017attosecond ; bredtmann2014x . Thus, inducing and controlling the
charge migration in molecules have emerged as one of the cornerstone topics of
attosecond science in recent years folorunso2021molecular ;
kraus2015measurement ; calegari2014ultrafast ; he2022filming .

Charge migration is connected to the transiently-evolving electronic flux
density (EFD) as a result of the quantum version of the continuity equation
sakurai2006advanced . The EFD is a vector field, which provides mechanistic
information about the spatial distribution of the charge migration along the
direction of flow. Over the years, EFD has emerged as an important tool to
gain detailed insights into the underlying mechanism of ultrafast charge
migration in systems ranging from molecules to nano-materials
nagashima2009electron ; giri2020time ; okuyama2009electron ; giri2021imaging ;
giri2022probing ; shao2020gal ; hermann2016ultrafast ; pohl2019imaging ;
shao2021bio ; sobottka2020tuning .

Charge migration of valence electrons and corresponding EFD in ring-shaped
molecules have received ample attention in the last few years. Charge
migration in magnesium porphyrin has been studied by means of EFD
nam2020monitoring ; tremblay2021time ; barth2006periodic ;
barth2006unidirectional ; koksal2017effect . Control over aromaticity in
benzene is understood by analysing charge migration ulusoy2011correlated .
Multidirectional angular EFD with a pincer-type flow during charge migration
in benzene has been observed, which further highlights that a significant
fraction of electrons is inactive during the flow. Furthermore, it has been
found that the nature and direction of EFD can be controlled by selective
excitation using optimised laser pulses hermann2016multidirectional ;
jia2017quantum ; hermann2017attosecond ; hermann2020probing . At this
juncture, it is natural to ponder how charge migration and the corresponding
EFD will differ if one reduces the symmetry of the molecule by transiting from
six-membered benzene to five-membered ring-shaped molecules. This is the main
aim of the present work.

Both cyclopentene and cyclopentadiene are homocyclic and nonaromatic
pentagonal molecules of carbon atoms. Moreover, both molecules are ionic in
nature, which makes them less attractive for the present study. On the other
hand, pyrrole, furan, and oxazole are five-membered ring-shaped neutral
molecules. Thus, these three molecules appear suitable to serve our purpose to
explore the role of symmetry reduction on charge migration and EFD. These
three heterocyclic aromatic molecules are present in many important biological
systems as a fundamental unit and have enormous applications demingos2021first
. As such, the static electronic properties of pyrrole, furan, and oxazole
have received significant attention from theoreticians
christiansen1999electronic ; roos2002theoretical ; lee1996molecular ;
serrano1993theoretical ; wan2000electronic ; burcl2002study ;
palmer1995electronic ; christiansen1998electronic . Further, a few non-
equilibrium dynamics in three molecules have been documented in Refs.
cao2016excited ; geng2020time ; carrascosa2021mapping ; barbatti2010non ;
barbatti2006nonadiabatic ; yong2021ultrafast .

The present work provides a detailed and comparative study of charge migration
of valence electrons in pyrrole, furan, and oxazole. In the following, we
demonstrate how the charge migration is sensitive to the symmetry reduction,
i.e., progressing from benzene with $\mathcal{D}_{\textrm{6h}}$ symmetry to
these five-membered ring-shaped molecules with $\mathcal{C}_{\textrm{2v}}$ and
$\mathcal{C}_{\textrm{s}}$ symmetries. Furthermore, we will discuss how the
nature of the electronic charge and flux densities during charge migration is
altered by the presence of heteroatoms, i.e., nitrogen in pyrrole, oxygen in
furan; and nitrogen and oxygen in oxazole.

##  2 Computational Methods

Charge migration is simulated by solving numerically the time-dependent
Schrödinger equation

$i\hbar\frac{\partial}{\partial
t}|\Psi(t)\rangle=\mathcal{H}(t)|\Psi(t)\rangle,$ (1)

where $\mathcal{H}(t)$ is the time-dependent Hamiltonian describing an
$N$-electron system interacting with an ultrashort laser pulse, and treated in
the semi-classical dipole approximation as

$\mathcal{H}(t)=\mathcal{H}_{0}-\hat{\mu}\cdot\mathbf{F}(t).$ (2)

Here, $\mathcal{H}_{0}$ is the field-free molecular Hamiltonian within the
clamped nuclei approximation, $\hat{\mu}$ is the dipole moment operator, and
$\mathbf{F}(t)$ is the electric field of the incident laser. An electronic
wave-packet is created in a molecule upon application of an ultrashort laser.
Within the hybrid time-dependent density functional theory-configuration
interaction (TDDFT/CI) method, the $N$-electron wave-packet is expressed as a
linear combination of many-body excited electronic states

$\Psi(\mathbf{r}^{N},t)=\sum_{k=0}^{N_{\textrm{states}}}C_{k}(t)\Phi_{k}(\mathbf{r}^{N}),$
(3)

with expansion coefficients $C_{k}(t)$. $N_{\textrm{states}}$ is the number of
many-electron states, $\Phi_{k}(\mathbf{r}^{N})$, necessary to describe the
wave-packet and its response in the presence of an ultrashort, intense laser
pulse. In the CI Singles picture, each $N$-electron excited state is written
as a linear combination of singly-excited configuration state functions (CSFs)
as

$\Phi_{k}(\mathbf{r}^{N})=D_{0,k}\Phi_{0}(\mathbf{r}^{N})+\sum_{ar}D_{a,k}^{r}\Phi_{a}^{r}(\mathbf{r}^{N}).$
(4)

The first term on the right-hand side represents the ground state Slater
determinant $\Phi_{0}(\mathbf{r}^{N})$, obtained from Kohn-Sham density
functional theory. The second term contains CSFs describing single excitations
from an occupied molecular orbital “a” to a virtual molecular orbital “r”. In
TDDFT/CI, the expansion coefficients $D_{0,k}$ and $D_{a,k}^{r}$ are obtained
using linear-response TDDFT in the Tamm-Dancoff approximation. This
combination was shown to provide good energetics while keeping the
computational cost low klinkusch2016 .

To analyse the laser-induced charge migration, the time-evolution of the one-
electron density $\rho(\mathbf{r},t)$ and the associated EFD $j(\mathbf{r},t)$
are reconstructed from the $N$-electron wavepacket dynamics. These two
observable quantities are connected via the quantum continuity equation as

$\frac{\partial}{\partial
t}\rho(\mathbf{r},t)=-\vec{\nabla}\cdot\mathbf{j}(\mathbf{r},t).$ (5)

Here, $\rho(\mathbf{r},t)$ is the expectation value of the one-electron
density operator as

$\displaystyle\hat{\rho}\left(\mathbf{r}\right)=\sum_{k=1}^{N}\delta\left(\mathbf{r}-\mathbf{r}_{k}\right).$
(6)

The latter also gives access to the electron density in a given state,
$\rho_{m}(\mathbf{r})=\langle\Phi_{m}|\hat{\rho}\left(\mathbf{r}\right)|\Phi_{m}\rangle$,
and the transition densities
$\rho_{mn}=\langle\Phi_{m}|\hat{\rho}\left(\mathbf{r}\right)|\Phi_{n}\rangle$.
Similarly, the EFD can be written in the following operator form

$\displaystyle\hat{j}\left(\mathbf{r}\right)=\frac{1}{2}\sum_{k=1}^{N}\left(\delta(\mathbf{r}-\mathbf{r}_{k})\hat{p}_{k}+\hat{p}^{\dagger}_{k}\delta(\mathbf{r}-\mathbf{r}_{k})\right),$
(7)

where $\hat{p}_{k}=-i\hbar\vec{\nabla}_{k}$ is the momentum operator of
electron $k$. The transient flux density maps are then the expectation values
of this operator.

To achieve the convergence of the laser-induced charge migration simulations,
the $N_{\textrm{states}}=$ 24, 22, and 32 lowest-lying bound excited states
are used for pyrrole, furan and oxazole, respectively. The CAM-B3LYP
functional and aug-cc-pVQZ basis sets on all atoms are used to compute these
electronic states 04:YTH:cam ; dunning1989ccpvxz . All electronic structure
calculations are performed with the Gaussian16 program frisch2016gaussian  and
the open-source toolbox detCI@ORBKIT is used post-process the information
obtained from the quantum chemistry package and to calculate all the required
observable hermann2016orbkit ; pohl2017open ; hermann2017open . The dynamical
simulations are performed using an in-house implementation of the TDCI method
08:TKS:rhoTDCI ; 10:TKKS:resp ; 11:TKKS:irhoTDCI .

The time evolution of the many-electron population dynamics in the three
molecules is shown for a 20 fs sine-squared pulse linearly polarised along the
$y$-axis. The carrier frequency of 6.26 eV, 6.26 eV, and 6.37 eV and peak
intensity of 3.7$\times$1013 W/cm2, 2.8$\times$1013 W/cm2, and 7$\times$1013
W/cm2 pulses are used for inducing charge migration in pyrrole, furan, and
oxazole, respectively.

##  3 Results and Discussion

Figure 1: Depiction of the electronic structure in the five-membered ring
molecules. 5${}^{\text{th}}$, 2${}^{\text{nd}}$, and 1${}^{\text{st}}$ excited
states of pyrrole, furan, and oxazole, respectively are represented from left
to right. (a) Ground state electron density (isocontours equally separated by
0.001 a.u.). (b) Natural transition orbitals. Orange (cyan) surfaces represent
the particle (hole) densities, respectively for isosurface value of 0.005. (c)
Difference electron densities between target excited and ground electronic
states, $\rho_{m}(\mathbf{r})-\rho_{0}(\mathbf{r})$. Violet(yellow) colour
represents regions of density increase (depletion) for isosurface value of
+0.005 (-0.005).

To compare the charge migration in pyrrole, furan, and oxazole on an equal
footing, a few conditions are considered. The first one is a static
consideration, where the target excited states are chosen to have a similar
character. This can be achieved by exciting states that, e.g., create nodal
structures around the heteroatoms. Further, in order to populate a target
excited state, it must be optically accessible. That is, the intensity of the
transition from the ground to target excited state (proportional to the
transition dipole moment $\mu^{2}$) must be high. From a dynamical
perspective, charge migration should be faster than the timescale of nuclear
rearrangement, ideally in the sub-femtosecond regime, yet accessible with a
few-femtosecond laser pulse. The target excited states satisfying both
constraints will thus have similar energies and transition dipoles. Further,
understanding the dynamics of the electronic density and EFD after such a pump
pulse is greatly simplified if the wavepacket is composed of only a few
components in the excited state manifold. This is achieved by exciting only a
small fraction of electronic populations and optimizing the laser parameters
(field intensity and carrier frequency), thereby leading to a controlled
degree of charge migration.

Figure 2: Laser-induced population dynamics in heterocyclic five-membered ring
molecules triggered by 20 fs sine-squared pulse linearly polarised along
$y$-axis. The label “resp” indicates the residual population of all other
electronic states, which acts as a non-linear response during the excitation.
The three molecules are shown in ball-and-stick representation, which are
projected in the $yz$-plane throughout this work. White, grey, blue, and red
colours represent hydrogen, carbon, nitrogen and oxygen, respectively. (a)
Pyrrole: carrier frequency resonant with the energy of the 5${}^{\text{th}}$
electronic state (6.26 eV); peak intensity of 3.7$\times$1013 W/cm2. (b)
Furan: carrier frequency resonant with the energy of the 2${}^{\text{nd}}$
electronic state (6.26 eV); peak intensity of 2.8$\times$1013 W/cm2. (c)
Oxazole: carrier frequency resonant with the energy of the 1${}^{\text{st}}$
electronic state (6.37 eV); peak intensity of 7$\times$1013 W/cm2.

The charge distribution corresponding to the ground electronic state of the
three molecules is presented in Fig. 1(a). All three molecules are projected
in the $yz$-plane throughout this work. Pyrrole and furan are symmetric with
respect to reflection about the $z$-axis that bisects molecule through the
heteroatom and the opposing CC bond, and the ground state densities both
transform according to the totally symmetric point IRREP of the
$\mathcal{C}_{\textrm{2v}}$ point group. The presence of the nitrogen atom
(blue) appears to favour the delocalisation of the electrons over the whole
ring in the $\pi$-orbitals. The more strongly electronegative oxygen atom
(red) leads to a more compact ground state density in furan (central panel),
which is more concentrated at the heteroatom. Due to the presence of both
oxygen and nitrogen atoms on the two sides of the pentagon, oxazole is not
symmetric with respect to a reflection in the plane and retains only the
operations of the $\mathcal{C}_{\textrm{s}}$ point group.

To identify interesting target state in each molecule, the natural transition
orbitals (NTOs) of  all states up to the ionisation threshold have been
inspected. The NTOs of the lowest excited state accessible via a linearly
polarised field along the $y$ axis are shown in Fig. 1(b). These correspond to
the fifth, second, and first excited states in pyrrole, furan, and oxazole,
respectively. The transition intensities of these target excited states are
much higher than the transition intensities from the ground to all other
excited states of the respective molecule. NTOs offer a compact depiction of
electronic transitions from a reference state in terms of particle and hole
densities martin2003natural . In all three molecules, the excitation indeed
creates nodal structure around the heteroatoms but the excitation characters
are radically different. The hole density (cyan) in pyrrole is found to be
delocalised on both sides of the ring on the carbon atoms, with a clear $\pi$
character and a node along the $y$-axis. The particle density (orange) is more
strongly localised close to the atoms of the ring, with a $\pi$ character
around the nitrogen atom. The structure is very different in furan, where the
hole density is localised around the carbon atoms and exhibits a stronger in-
plane bonding character (reminiscent of a $\sigma$ character, in opposition to
an out-of-plane $\pi$-character). The same character is found in the particle
density, which is delocalised along the C-O-C fragment, as well as through
space between the two C=C bonds. Both excitations in pyrrole and furan retain
the structure of the $\mathcal{C}_{\textrm{2v}}$ point group. The NTO
densities in oxazole are more similar to the furan case, although it is the
hole density that has the through-space structure.

Figure 3: Electronic flux density and time-dependent difference electronic
density associated with charge migration in pyrrole at $\tau$, 2$\tau$, and
6$\tau$. Here, $\tau=\textrm{T}/8$ with $\textrm{T}=660$ attoseconds as the
characteristic timescale of the charge migration in pyrrole. Ground-state
electron density is subtracted from the time-dependent difference electronic
density at each subsequent time-step. The colour of the streamline arrows
varies from white to black according to their increasing intensity. The violet
(yellow) contours present the enhancement (depletion) of the electronic
density with respect to the ground state density. For better visualisation,
the projection of pyrrole in the $yz$-plane is superposed in ball-stick
representation.

The difference densities between the ground and target states are also shown
in Fig. 1(c). The strong $\pi$-type characters, i.e., the out-of-plane bonding
characters of the excitation in pyrrole is noticed. In opposition to $\pi$
character, a clear nodal structure in the plane of the molecule analogous to
the $\sigma$ character of the excitations in furan and oxazole can be readily
recognised. More importantly, the difference densities reveal the symmetry of
the excited states, which are found to transform according to the B2 IRREP of
the $\mathcal{C}_{\textrm{2v}}$ point group for pyrrole and furan. These
important differences should lead to markedly different charge migration
dynamics after laser excitation.

We employ a linearly polarised pulse along the $y$ direction to induce the
charge migration. The resultant population dynamics are presented in Fig. 2.
The pulses for the three molecules have comparable peak intensities, with
similar carrier frequencies tuned at the respective target transitions. The
peak intensities of the pulses, in excess of 1013W/cm2, are tuned to obtain a
similar population in the target states, hence leading to a similar contrast
in the charge migration patterns. As evident from the figure, the target
excited states reach around $4\%$ population in all three cases. The values of
the transition dipole ($\mu^{2}$) to the target excited states are 0.9898
$ea_{0}$, 1.1071 $ea_{0}$, and 1.0048 $ea_{0}$ for pyrrole (fifth), furan
(second), and oxygen (first), respectively. This compares favourably to all
other transitions with values lower than 0.25 $ea_{0}$ in all three molecules.

The population in all off-resonant electronic states are termed as “resp” in
the figure, as they mediate the non-linear response of the electronic density
to the strong external electric field. The electronic response is found to be
more important in oxazole than in either pyrrole or furan. During the pulse,
the population of these non-resonant states is exactly zero for pyrrole and
furan, but non-zero in oxazole. This is because of the strong transition
dipole moments to other excited states in oxazole, which stems from its
different chemical nature compared to the other two molecules. The field-free
charge migration sets into motion at time zero, after the pulse. Since the
target excited states in the three systems have similar excitation energies
with respect to the ground electronic state, the charge migration will have a
similar oscillation period around 650 as.

To shed light on the effect of symmetry on field-free charge migration on its
characteristic timescale, we analyse snapshots of the EFD and of the
associated time-dependent difference electronic density. For the latter, the
contribution from the electronic ground state is subtracted from the wave
packet at each time-step.

Figure 4: Same as Fig. 3 for furan with $\textrm{T}=660$ attoseconds as the
characteristic timescale.

Figure 3 presents the EFD and the difference density at three time steps for
pyrrole. The overall behaviour of the difference densities at all times
reflects that they are mostly located around the atoms, which indicates charge
localization. This corresponds well to the picture offered by the density
difference in the left panel of Fig. 1(c). The strong charge localisation
around the atoms should not affect strongly the bonding properties in pyrrole.
Consequently, it is unlikely that the excitation chosen would lead to any bond
breaking, ring-opening or other photochemical reaction. Moreover, the
difference densities are positive and negative in alternate carbon atoms,
whereas nitrogen is always surrounded by positive electron density. At $\tau$
and 2$\tau$, the positive contour around nitrogen is tilted towards the
positive $y$-axis as evident from the figure. This situation is reversed at
6$\tau$. The difference densities around carbon atoms are anti-symmetric about
$y$-axis at all times. The overall intensity of the difference density reduces
from $\tau$ to $2\tau$. There is a phase shift of ${\pi}/{2}$ from $2\tau$ to
$6\tau$, which is in accordance with the half time-period of the charge
migration as evident from the figure. The characteristic timescale form charge
migration in pyrrole is $\textrm{T}=660$ attoseconds for the specific
superposition state created here.

The electron flow can be understood better in terms of the EFD. As evident
from Fig. 3, the direction of the electron flow is mostly around the atoms.
The streamlines start from the atomic positions and terminate around the bonds
at all times. This reflects the atom-to-bond charge migration in pyrrole,
which corresponds well to the picture offered by the NTO in Fig. 1(b). Indeed,
the particle density is observed to be strongly localised around the atoms.
Consequently, the streamlines appear as piecewise disjoint lines of flow
around each atom. From a symmetry point of view, there is no general electron
flow from left to right or from upward to downward. At $\tau$, the flow around
the bottom two carbon atoms is upward inside and downward outside of the
pyrrole ring. Similar to the difference densities, the streamlines follow the
same phase reversal from $2\tau$ to $6\tau$. The EFD transforms according to
the totally symmetric IRREP of the $\mathcal{C}_{\textrm{2v}}$ point group at
all times.

At this juncture, it is natural to ask how the above findings will alter if
one replaces nitrogen by a more electronegative atom, oxygen. To explore the
question further, let us analyse EFD and the difference density for furan at
three time steps as presented in Fig. 4. The electronic ground state density
and the NTO of the target transition in Fig. 1 already hint at a significant
change in the charge migration dynamics. Unlike for pyrrole, the positive and
negative contours of the difference electronic densities are distributed not
only around the atoms but also around the bonds, and the latter is higher in
intensity. They are anti-symmetric about the $y$-axis except around the oxygen
atom. These findings are in agreement with the density difference profile
observed in Fig. 1(c).

The bond strengthening (weakening) shown by this increase (decrease) in the
electronic densities around the bonds is primarily visible on the two sides of
the oxygen atom. In contrast to pyrrole, these will affect the bond strength
and potentially cause bonds to break, resulting in ring opening on either side
of the heteroatom. Exciting furan with short, intense laser pulses of the form
used in this work is therefore more likely to trigger photochemical processes
on longer timescales. Interestingly, the difference density around the oxygen
is negative at all times, reflecting the larger electronegativity of oxygen
compared to the carbon atoms of the ring. Since the oxygen substitution
preserves all operations of the $\mathcal{C}_{\textrm{2v}}$ point group, the
density difference upon charge migration in furan retains the same B2 IRREP as
for pyrrole.

From the EFD for furan displayed in Fig. 4, the EFD reveals a swirling type of
motion between the atoms of the ring. The streamlines show a flow downward
through the bonds between the oxygen and its nearest carbon atoms, and going
upward to the hydrogen atoms at time steps $\tau$ and $2\tau$. The direction
of the flow is completely reversed at 6$\tau$. These patterns are consistent
with the NTO picture in Fig. 1(b), which show a stronger in-plane bonding
character and more nodal structure in the plane of the molecule for the target
excitation. Despite the apparent complexity of the flow patterns, they
transform according to the totally symmetric IRREP of the
$\mathcal{C}_{\textrm{2v}}$ point group at all times. It can be concluded that
the stronger in-plane bonding character of the target excitation in a ring
bearing a more electronegative atom does not affect the symmetry of the charge
migration.

Figure 5: Same as Fig. 3 for oxazole with $\textrm{T}=650$ attoseconds as the
characteristic timescale.

For the charge migration in pyrrole and furan, both molecules belong to
$\mathcal{C}_{\textrm{2v}}$ point group. The symmetry of the five-membered
ring can be reduced by substitution of two carbon atoms by both nitrogen and
oxygen. To understand this further symmetry reduction, let us explore the
charge migration in oxazole, which has $\mathcal{C}_{\textrm{s}}$ symmetry.

As a result of the symmetry reduction, the charge migration in oxazole changes
drastically as two carbon atoms are replaced by one nitrogen and one oxygen
atoms on the opposite sides of the pentagon ring (see Fig. 5). Due to the
presence of the different atoms in the ring, the difference electron densities
are entirely different in comparison to pyrrole and furan, and do not exhibit
any symmetric behaviour upon rotation with respect to the $z$-axis or
reflection in the $xz$-plane. The regions of density increase/depletion are
not specific to the bonds or atoms. The dumbbell structure observed at all
time step in the electronic density different around nitrogen atom hints at a
strong contribution of a $p$-type orbital in the molecular plane. At $\tau$,
the negative contours are present around the bonds and positive contours are
situated in the space between nitrogen, oxygen and carbon atoms. In general,
there is no strong correlation between the difference density plots during
charge migration and the ones reported in Fig. 1(c).

To have a better understanding of the charge migration mechanism in oxazole,
let us analyse the time-resolved EFD after laser excitation, as depicted in
Fig. 5. The multiple directions of the flux densities provide a picture of the
charge migration that is consistent with the NTO densities depicted in Fig.
1(b). A swirling motion is observed, separated in bottom and top contributions
that barely exchange electrons at all times. This is intriguing, since the
excitation using a $y$-polarised pulse would likely create nodal structure
along this direction. As such, a left-right separation of the charge migration
should be expected, as was the case for pyrrole and furan. From the NTOs in
Fig. 1(b), it appears that the top part of the particle density only migrates
towards the hole density located between the two heteroatoms. We attribute
this feature to their electronegativity, that prefer to retain the excited
electron localised in this fragment. On the other hand, the bottom part of the
charge migration sees a transfer of the particle localised on the CO bond to
the CC bond across the ring. This part of the charge migration operates
through space rather than following the bonds, in stark contrast with the more
symmetrical pyrrole and furan. The bonds connecting to the oxygen atom
experience the strongest bond strength oscillations, whereas bonds connected
to the nitrogen atom see less pronounced fluctuations. Much like the case of
furan, it can be inferred from the electronic density and EFD snapshots that
oxazole could probably undergo further photochemical reactions.

##  4 Conclusion

In summary, we have investigated the role of symmetry reduction on charge
migration in selected heterocyclic five-membered ring molecules. Time-
dependent transition electron densities and electronic flux densities were
used to unravel the mechanism of charge migration induced by selective laser
pulses. To compare the charge migration on three molecules on equal footing,
transitions with similar energetics were used to reach target excited states
with high transition intensities and nodal structure close to the heteroatoms.
It has been found that laser-induced charge migration and the corresponding
EFD are significantly different in the three molecules in spite of these
constraints. Moreover, variations in electronic densities and EFDs reveal bond
strengthening or weakening upon excitation by light, which is useful to infer
photochemical potential in these molecules. For the laser pulses used in this
work, it was inferred that pyrrole is less likely to undergo ring-opening
reactions than furan and oxazole.

The presence of nitrogen in pyrrole and oxygen in furan make the charge
migration in both molecules significantly different. In pyrrole, a signature
of charge localisation is observed during charge migration as the difference
electronic densities are positioned around the atoms. In contrast, electron
density depletion/increase regions are delocalised around the atoms and around
the bonds during light-induced charge migration in furan, which indicates the
bond-to-atom charge transfer mechanism. The charge migration patterns observed
in both cases are totally symmetric. For oxazole, the difference electronic
densities documenting the charge migration dynamics display no symmetry
outside of the $yz$-plane as a result of the lower symmetry, from
$\mathcal{C}_{\textrm{2v}}$ in furan and pyrrole to
$\mathcal{C}_{\textrm{s}}$. In this case, charge migration takes the form of a
swirling motion through space.

We believe that our results on attosecond charge migration in heterocyclic
five-membered rings will motivate further theoretical investigations in other
types of molecular systems, where studies of structure-dynamics relationships
are still few and far apart. Further, state-of-the-art emerging experimental
techniques, such as attosecond transient absorption spectroscopy
kraus2015measurement ; calegari2014ultrafast , high-harmonic generation
spectroscopy folorunso2021molecular ; he2022filming ; dixit2018control ;
chandra2019experimental , or time-resolved x-ray diffraction dixit2012imaging
; dixit2014theory , to name but a few, will provide potent probes for our
findings in heterocyclic five-membered ring molecules in the coming future.

Acknowledgments

S. G. acknowledges Council of Scientific and Industrial Research (CSIR) for
senior research fellowship (SRF). G. D. acknowledges support from Science and
Engineering Research Board (SERB) India (Project No. MTR/2021/000138).

## Declarations

The datasets generated during and/or analysed during the current study are
available from the corresponding author on reasonable request.

## References

  * (1) F. Krausz and M. Ivanov, “Attosecond physics,” Reviews of Modern Physics, vol. 81, no. 1, p. 163, 2009.
  * (2) L. S. Cederbaum and J. Zobeley, “Ultrafast charge migration by electron correlation,” Chemical Physics Letters, vol. 307, no. 3-4, pp. 205–210, 1999.
  * (3) A. I. Kuleff, J. Breidbach, and L. S. Cederbaum, “Multielectron wave-packet propagation: General theory and application,” J. Chem. Phys., vol. 123, no. 4, p. 044111, 2005.
  * (4) J. Breidbach and L. Cederbaum, “Migration of holes: Formalism, mechanisms, and illustrative applications,” J. Chem. Phys., vol. 118, no. 9, pp. 3983–3996, 2003.
  * (5) F. Remacle and R. D. Levine, “An electronic time scale in chemistry,” Proceedings of the National Academy of Sciences, vol. 103, no. 18, pp. 6793–6798, 2006.
  * (6) J. Breidbach and L. Cederbaum, “Universal attosecond response to the removal of an electron,” Phys. Rev. Lett., vol. 94, no. 3, p. 033901, 2005.
  * (7) M. Nisoli, P. Decleva, F. Calegari, A. Palacios, and F. Martín, “Attosecond electron dynamics in molecules,” Chemical Reviews, vol. 117, no. 16, pp. 10760–10825, 2017.
  * (8) T. Bredtmann, M. Ivanov, and G. Dixit, “X-ray imaging of chemically active valence electrons during a pericyclic reaction,” Nature Communications, vol. 5, no. 1, pp. 1–7, 2014.
  * (9) 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,” Physical Review Letters, vol. 126, no. 13, p. 133002, 2021.
  * (10) P. M. Kraus, B. Mignolet, D. Baykusheva, A. Rupenyan, L. Hornỳ, E. F. Penka, G. Grassi, O. I. Tolstikhin, J. Schneider, F. Jensen, et al., “Measurement and laser control of attosecond charge migration in ionized iodoacetylene,” Science, vol. 350, no. 6262, pp. 790–795, 2015.
  * (11) F. Calegari, D. Ayuso, A. Trabattoni, L. Belshaw, S. De Camillis, S. Anumula, F. Frassetto, L. Poletto, A. Palacios, P. Decleva, et al., “Ultrafast electron dynamics in phenylalanine initiated by attosecond pulses,” Science, vol. 346, no. 6207, pp. 336–339, 2014.
  * (12) L. He, S. Sun, P. Lan, Y. He, B. Wang, P. Wang, X. Zhu, L. Li, W. Cao, P. Lu, et al., “Filming movies of attosecond charge migration in single molecules with high harmonic spectroscopy,” Nature Communications, vol. 13, no. 1, pp. 1–9, 2022.
  * (13) J. J. Sakurai, Advanced quantum mechanics.  Pearson Education India, 2006.
  * (14) K. Nagashima and K. Takatsuka, “Electron-wavepacket reaction dynamics in proton transfer of formamide,” J. Phys. Chem. A, vol. 113, no. 52, pp. 15240–15249, 2009.
  * (15) S. Giri, A. M. Dudzinski, J. C. Tremblay, and G. Dixit, “Time-dependent electronic current densities in chiral molecules,” Phys. Rev. A, vol. 102, no. 6, p. 063103, 2020.
  * (16) M. Okuyama and K. Takatsuka, “Electron flux in molecules induced by nuclear motion,” Chem. Phys. Lett., vol. 476, no. 1-3, pp. 109–115, 2009.
  * (17) S. Giri, J. C. Tremblay, and G. Dixit, “Imaging charge migration in chiral molecules using time-resolved x-ray diffraction,” Phys. Rev. A, vol. 104, no. 5, p. 053115, 2021.
  * (18) S. Giri, J. C. Tremblay, and G. Dixit, “Probing the effect of molecular structure saddling on ultrafast charge migration via time-resolved x-ray diffraction,” Phys. Rev. A, vol. 106, p. 033120, 2022.
  * (19) J. Shao, V. Pohl, L. E. Marsoner Steinkasserer, B. Paulus, and J. C. Tremblay, “Electronic current mapping of transport through defective zigzag graphene nanoribbons,” J. Phys. Chem. C, vol. 124, no. 43, pp. 23479–23489, 2020\.
  * (20) G. Hermann and J. Tremblay, “Ultrafast photoelectron migration in dye-sensitized solar cells: Influence of the binding mode and many-body interactions,” J. Chem. Phys., vol. 145, no. 17, p. 174704, 2016.
  * (21) V. Pohl, L. E. Marsoner Steinkasserer, and J. C. Tremblay, “Imaging time-dependent electronic currents through a graphene-based nanojunction,” J. Phys. Chem. Lett., vol. 10, no. 18, pp. 5387–5394, 2019.
  * (22) J. Shao, B. Paulus, and J. C. Tremblay, “Local current analysis on defective zigzag graphene nanoribbons devices for biosensor material applications,” J. Comput. Chem., vol. 42, pp. 1475–1485, 2021.
  * (23) S. Sobottka, M. Nößler, A. L. Ostericher, G. Hermann, N. Z. Subat, J. Beerhues, M. Behr-van der Meer, L. Suntrup, U. Albold, S. Hohloch, et al., “Tuning ptii-based donor–acceptor systems through ligand design: Effects on frontier orbitals, redox potentials, uv/vis/nir absorptions, electrochromism, and photocatalysis,” Chemistry–A European Journal, vol. 26, no. 6, pp. 1314–1327, 2020.
  * (24) Y. Nam, J. R. Rouxel, J. Y. Lee, and S. Mukamel, “Monitoring aromatic ring-currents in mg-porphyrin by time-resolved circular dichroism,” Phys. Chem. Chem. Phys., vol. 22, no. 45, pp. 26605–26613, 2020.
  * (25) J. C. Tremblay, V. Pohl, G. Hermann, and G. Dixit, “Time-resolved imaging of correlation-driven charge migration in light-induced molecular magnets by x-ray scattering,” Faraday Discussions, vol. 228, pp. 82–103, 2021.
  * (26) I. Barth and J. Manz, “Periodic electron circulation induced by circularly polarized laser pulses: Quantum model simulations for mg porphyrin,” Angewandte Chemie International Edition, vol. 45, no. 18, pp. 2962–2965, 2006\.
  * (27) I. Barth, J. Manz, Y. Shigeta, and K. Yagi, “Unidirectional electronic ring current driven by a few cycle circularly polarized laser pulse: quantum model simulations for mg- porphyrin,” Journal of the American Chemical Society, vol. 128, no. 21, pp. 7043–7049, 2006.
  * (28) K. Köksal and F. Koç, “The effect of twisted light on the ring-shaped molecules: the manipulation of the photoinduced current and the magnetic moment by transferring spin and orbital angular momentum of high frequency light,” Computational and Theoretical Chemistry, vol. 1099, pp. 203–208, 2017.
  * (29) I. S. Ulusoy and M. Nest, “Correlated electron dynamics: How aromaticity can be controlled,” Journal of the American Chemical Society, vol. 133, no. 50, pp. 20230–20236, 2011.
  * (30) G. Hermann, C. Liu, J. Manz, B. Paulus, J. F. Perez-Torres, V. Pohl, and J. C. Tremblay, “Multidirectional angular electronic flux during adiabatic attosecond charge migration in excited benzene,” J. Phys. Chem. A, vol. 120, no. 27, pp. 5360–5369, 2016.
  * (31) 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., vol. 482, pp. 146–159, 2017.
  * (32) G. Hermann, C. Liu, J. Manz, B. Paulus, V. Pohl, and J. C. Tremblay, “Attosecond angular flux of partial charges on the carbon atoms of benzene in non-aromatic excited state,” Chem. Phys. Lett., vol. 683, pp. 553–558, 2017.
  * (33) G. Hermann, V. Pohl, G. Dixit, and J. C. Tremblay, “Probing electronic fluxes via time-resolved x-ray scattering,” Phys. Rev. Lett., vol. 124, no. 1, p. 013002, 2020.
  * (34) P. G. Demingos, N. M. Balzaretti, and A. R. Muniz, “First-principles study of carbon nanothreads derived from five-membered heterocyclic rings: thiophene, furan and pyrrole,” Phys. Chem. Chem. Phys., vol. 23, no. 3, pp. 2055–2062, 2021.
  * (35) O. Christiansen, J. Gauss, J. F. Stanton, and P. Jorgensen, “The electronic spectrum of pyrrole,” J. Chem. Phys., vol. 111, no. 2, pp. 525–537, 1999\.
  * (36) B. O. Roos, P.-Å. Malmqvist, V. Molina, L. Serrano-Andrés, and M. Merchán, “Theoretical characterization of the lowest-energy absorption band of pyrrole,” J. Chem. Phys., vol. 116, no. 17, pp. 7526–7536, 2002.
  * (37) S. Y. Lee and B. H. Boo, “Molecular structures and vibrational spectra of pyrrole and carbazole by density functional theory and conventional ab initio calculations,” J. Phys. Chem., vol. 100, no. 37, pp. 15073–15078, 1996\.
  * (38) L. Serrano-Andres, M. Merchan, I. Nebot-Gil, B. O. Roos, and M. Fulscher, “Theoretical study of the electronic spectra of cyclopentadiene, pyrrole, and furan,” Journal of the American Chemical Society, vol. 115, no. 14, pp. 6184–6197, 1993.
  * (39) J. Wan, J. Meller, M. Hada, M. Ehara, and H. Nakatsuji, “Electronic excitation spectra of furan and pyrrole: Revisited by the symmetry adapted cluster–configuration interaction method,” J. Chem. Phys., vol. 113, no. 18, pp. 7853–7866, 2000.
  * (40) R. Burcl, R. D. Amos, and N. C. Handy, “Study of excited states of furan and pyrrole by time-dependent density functional theory,” Chem. Phys. Lett., vol. 355, no. 1-2, pp. 8–18, 2002.
  * (41) M. H. Palmer, I. C. Walker, C. C. Ballard, and M. F. Guest, “The electronic states of furan studied by vuv absorption, near-threshold electron energy-loss spectroscopy and ab initio multi-reference configuration interaction calculations,” Chem. Phys., vol. 192, no. 2, pp. 111–125, 1995\.
  * (42) O. Christiansen and P. Jørgensen, “The electronic spectrum of furan,” Journal of the American Chemical Society, vol. 120, no. 14, pp. 3423–3430, 1998\.
  * (43) J. Cao, Z.-Z. Xie, and X. Yu, “Excited-state dynamics of oxazole: A combined electronic structure calculations and dynamic simulations study,” Chem. Phys., vol. 474, pp. 25–35, 2016.
  * (44) T. Geng, J. Ehrmaier, O. Schalk, G. W. Richings, T. Hansson, G. Worth, and R. D. Thomas, “Time-resolved photoelectron spectroscopy studies of isoxazole and oxazole,” J. Phys. Chem. A, vol. 124, no. 20, pp. 3984–3992, 2020\.
  * (45) A. M. Carrascosa, M. Yang, H. Yong, L. Ma, A. Kirrander, P. M. Weber, and K. Lopata, “Mapping static core-holes and ring-currents with x-ray scattering,” Faraday Discussions, vol. 228, pp. 60–81, 2021.
  * (46) M. Barbatti, J. Pittner, M. Pederzoli, U. Werner, R. Mitrić, V. Bonačić-Kouteckỳ, and H. Lischka, “Non-adiabatic dynamics of pyrrole: Dependence of deactivation mechanisms on the excitation energy,” Chem. Phys., vol. 375, no. 1, pp. 26–34, 2010.
  * (47) M. Barbatti, M. Vazdar, A. J. Aquino, M. Eckert-Maksić, and H. Lischka, “The nonadiabatic deactivation paths of pyrrole,” J. Chem. Phys., vol. 125, no. 16, p. 164323, 2006.
  * (48) H. Yong, S. M. Cavaletto, and S. Mukamel, “Ultrafast valence-electron dynamics in oxazole monitored by x-ray diffraction following a stimulated x-ray raman excitation,” J. Phys. Chem. Lett., vol. 12, no. 40, pp. 9800–9806, 2021\.
  * (49) S. Klinkusch and J. C. Tremblay, “Resolution-of-identity stochastic time-dependent configuration interaction for dissipative electron dynamics in strong fields,” J. Chem. Phys., vol. 144, no. 18, p. 184108, 2016.
  * (50) T. Yanai, D. Tew, and N. Handy, “A new hybrid exchange-correlation functional using the coulomb-attenuating method (cam-b3lyp),” Chem. Phys. Lett., vol. 393, pp. 51–57, 2004.
  * (51) J. Dunning and H. Thom, “Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen,” J. Chem. Phys., vol. 90, p. 1007, 1989.
  * (52) M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, et al., “Gaussian 16, revision a. 03, gaussian,” Inc., Wallingford CT, 2016.
  * (53) G. Hermann, V. Pohl, J. C. Tremblay, B. Paulus, H. C. Hege, and A. Schild, “Orbkit: A modular python toolbox for cross-platform postprocessing of quantum chemical wavefunction data,” J. Comp. Chem., vol. 37, no. 16, pp. 1511–1520, 2016.
  * (54) V. Pohl, G. Hermann, and J. C. Tremblay, “An open-source framework for analyzing n-electron dynamics. i. multideterminantal wave functions,” J. Comput. Chem., vol. 38, no. 17, pp. 1515–1527, 2017.
  * (55) G. Hermann, V. Pohl, and J. C. Tremblay, “An open-source framework for analyzing n-electron dynamics. ii. hybrid density functional theory/configuration interaction methodology,” J. Comput. Chem., vol. 38, no. 28, pp. 2378–2387, 2017.
  * (56) J. C. Tremblay, T. Klamroth, and P. Saalfrank, “Time-dependent configuration-interaction calculations of laser-driven dynamics in presence of dissipation,” J. Chem. Phys., vol. 129, pp. 084302:1–8, 2008.
  * (57) J. C. Tremblay, P. Krause, T. Klamroth, and P. Saalfrank, “The effect of energy and phase relaxation on dynamic polarizability calculations,” Phys. Rev. A, vol. 81, pp. 063420:1–10, 2010.
  * (58) J. C. Tremblay, S. Klinkusch, T. Klamroth, and P. Saalfrank, “Dissipative many-electron dynamics of ionizing systems,” J. Chem. Phys., vol. 134, p. 044311, 2011.
  * (59) R. L. Martin, “Natural transition orbitals,” J. Chem. Phys., vol. 118, no. 11, pp. 4775–4777, 2003.
  * (60) G. Dixit, Á. Jiménez-Galán, L. Medišauskas, and M. Ivanov, “Control of the helicity of high-order harmonic radiation using bichromatic circularly polarized laser fields,” Phys. Rev. A, vol. 98, no. 5, p. 053402, 2018.
  * (61) S. Chandra, I. N. Ansari, G. Dixit, F. Lepine, and A. Bhattacharya, “Experimental evidence of sensitivity of the high harmonic generation to hydrogen bonding,” J. Phys. Chem. A, vol. 123, no. 24, pp. 5144–5149, 2019\.
  * (62) G. Dixit, O. Vendrell, and R. Santra, “Imaging electronic quantum motion with light,” Proceedings of the National Academy of Sciences, vol. 109, no. 29, pp. 11636–11640, 2012.
  * (63) G. Dixit, J. M. Slowik, and R. Santra, “Theory of time-resolved nonresonant x-ray scattering for imaging ultrafast coherent electron motion,” Phys. Rev. A, vol. 89, no. 4, p. 043409, 2014.