Annihilation of positrons from 22Na in novae

N. N. Chugai1 and A. D Kudryashov2
1Institute of Astronomy, Russian Academy of Sciences, Pyatnitskaya ul. 48,
Moscow, 119017 Russia e-mail<EMAIL_ADDRESS>
2Russian Institute of scientific and technical Information, Russian Academy of
Sciences, Moscow

Abstract

We explore for the first time effects of the magnetic field on the escape of
22Na positrons and on the flux evolution of annihilation 511 keV line in
novae. It is shown that for the white dwarf magnetic field of $\sim 10^{6}$ G
the field of the expanding nova shell is able to significantly impede
positrons escape and increase the time of the nova emission in 511 keV up to
hundreds days. Keywords: stars — novae — nucleosynthesis

##  1 Introduction

Nova phenomenon is caused by a thermonuclear runaway (TNR) in a hydrogen shell
accumulated on a white dwarf (WD) by an accretion in a binary system with a
normal companion (usually red dwarf). Despite a general understanding of the
nova outburst, widely accepted models for major processes involved in the nova
phenomenon are lacking. Particularly, a mixing mechanism of the WD matter with
the hydrogen shell that is crucial for the TNR is still debatable. Moreover,
it is not clear, whether the nova ejects less or more material than
accumulated between consecutive outbursts (Epelstain et al. 2007; Starrfield
et al. 2020). Uncertainties remain also in the description of the shell
ejection and shell morphology.

Given the theoretical status, a significant role in the theoretical progress
belongs to observations. This is supported by two fundamental facts laid in
the basis of the nova theory: (i) novae are binary system containing the WD
(Kraft 1964); (ii) nova shell material is enriched by a factor of 100 with CNO
elements (Mustel & Baranova 1966), which indicates a crucial role of mixing
between accreting and WD matter for the nova energetics (Sparks et al. 1976).
It should be emphasised that apart from novae related to the TNR on the CO WD
there is an abundant class of novae related to the TNR on ONe dwarfs
(Starrfield et al. 1986). This conclusion is prompted by the detection of Ne
overabundance in spectra of some novae (Ferland & Shields 1978). Noteworthy,
according to the interesting conjecture (Shara & Prialnik 1994) neon novae
could arise from CO WD with outer layers significantly enriched by Ne and Mg
due to the preceded hydrogen accretion with the rate of $\sim 10^{-6}$
$M_{\odot}$ yr-1.

The nova TNR results in the synthesis of radioactive isotopes which decay with
half life from hundreds of seconds (13N) till several years (22Na). As far
back as half century Clayton & Hoyle (1974) emphasised the importance of
observations of gamma-lines from novae. The isotope composition of nova shell
depends on WD type (CO or ONe) and on some poorly fixed parameters,
particularly, the WD mass, the accretion rate and the admixing of the WD
material into the hydrogen shell. A detection of gamma-lines from novae would
permit us to determine the mass of a certain isotope and thus to obtain useful
constraints on outburst models.

Highly prospective for gamma-line observations is the isotope 22Na with the
half life of 2.6 years that provides a possibility to apply a long integration
time. The outburst theory predicts a high 22Na abundance in neon novae, $\sim
10^{-3}$ by mass (Denissenkov et al. 2014; Kudryashov 2019), which is by a
three order higher compared to that of CO novae. The 22Na decay is accompanied
by the emission of the 1275 keV gamma-line, the 511 keV annihilation line, and
three-photon annihilation continuum. Five neon novae have been observed with
the gamma-ray telescope COMPTEL resulting in the upper limit of the 22Na mass
between $3\cdot 10^{-8}$ $M_{\odot}$ and $2\cdot 10^{-6}$ $M_{\odot}$ based on
the $2\sigma$ upper limit of the flux in the 1275 keV line (Iyudin et al.
1995). These upper limits do not contradict to theoretical 22Na abundance in
current models of neon novae provided a total ejecta mass is $<3\cdot 10^{-5}$
$M_{\odot}$. Subsequently, Iyudin (2010) reported on the COMPTEL detection of
1275 keV line at the confidence level of 4$\sigma$ from the slow nova V723 Cas
of 1975, which is rather unexpected result for these category of novae.

The model luminosity of novae in the annihilation 511 keV line of 22Na has
been addressed by Gomez et al. (1998) who demonstrated that already a week
after the TNR the nova shell becomes transparent for positrons that escape
without annihilation. If this is the case, the possibility to detect 511 keV
line from 22Na in novae becomes doubtful. On the other hand, the conclusion
about the high transparency of the nova shell for 22Na positrons is made
neglecting a possible magnetic field in the shell. It may well be that a
moderate value of the magnetic field would be enough to essentially impede the
positrons escape and thus to prolong the time of a significant nova luminosity
in 511 keV line. In this case the detection possibility for 511 keV line from
neon novae would become feasible.

This paper for the first time explores the issue, to which extent the
inclusion of the magnetic field could affect the escape of 22Na positrons from
the shell and the evolution of 511 keV line flux from the neon nova. The shell
model with the frozen-in magnetic field and the prescription for the
annihilation are considered in Section 2. The computations of the 511 keV flux
evolution for different values of relevant parameters are presented in Section
3.

##  2 Model overview

###  2.1 Nova shell with frozen-in magnetic field

The issue of the shell structure and kinematics for different subclasses of
novae currently is not completely clear. The most natural is the assumption
that the ejecta forms in the regime of the optically thick wind (Ruggles &
Bath 1979; Kato & Hachisu 2007). The idea is based on the fact that the major
stage of the TNR is longer compared to the hydrodynamic time at the WD
surface. Following this arguments we assume that the ejected shell is produced
by the wind with the constant mass loss rate and constant velocity of $v\sim
1500$ km s-1 that forms at the outer boundary of the expanded envelope with
the radius $r_{1}\sim 2GM/v^{2}\sim 10^{10}$ cm, where $M\sim 1$ $M_{\odot}$
is the WD mass. The convenient wind parameter is the wind kinetic luminosity
$L_{w}=0.5wv^{3}$, where $w=\dot{M}/v=4\pi r^{2}\rho$ is the wind density
parameter. The $L_{w}$ value should be of the order of the Eddington
luminosity, $1.26\cdot 10^{38}(M/M_{\odot})$ erg s-1. For the standard model
we adopt $L_{w}=10^{38}$ erg s-1. It takes thus $t_{w}=M_{s}/\dot{M}\sim 26$
days for the wind to form ejecta with the mass of $M_{s}\sim 10^{-5}$
$M_{\odot}$.

For the adopted wind velocity the ejecta perturbation from the binary
components are insignificant. Indeed, for the typical orbital period of 3.5 h,
the WD mass of $M_{1}=1$ $M_{\odot}$, and the red dwarf mass of $M_{2}=0.5$
$M_{\odot}$ one gets the large semi-axis $a={9\cdot 10^{10}}$ cm, and orbital
velocities of components $v_{1}=150$ km s-1 and $v_{2}=300$ km s-1. The
maximal variation of the wind velocity caused by the WD motion is about 10%
close to the orbital plane. The variation of the density and velocity caused
by the gravitation from the secondary is also small due to high wind velocity.
Below we adopt that the shell expands spherically with the constant velocity.
Apropos, the neon nova V382 Vel (Takeda & Diaz 2019) shows the spherical shell
expanding with the velocity of 1200 km s-1. A possible deviation from the
spherical symmetry (cf. Chomiuk et al 2020) is discussued in the final
section.

In our picture the hydrogen shell before the outburst resides in the magnetic
field of the WD with the average strength of $B_{0}=10^{6}$ G. This value is
rather moderate compared to the field of polars and intermediate polars of
$3\cdot 10^{6}-2\cdot 10^{8}$ G. (Ferrario et al. 2020). The assumption on the
magnetic field value is also consistent with the magnetic momentum of the WD
of classical nova 1934 DQ Her, viz., $(1-3)\cdot 10^{32}$ G cm3 (Patterson
1994). For the WD radius of $\sim 4\cdot 10^{8}$ cm this value corresponds to
the magnetic field on the WD surface of $(1-3)\times 10^{6}$ G.

The TNR is accompanied by a vigorous convection all over the hydrogen shell
(Casanova et al. 2018; Starrfield et al. 2020). The adopted magnetic field
($\sim 10^{6}$ G) is weak in the energetic and dynamic respect, so the
convection unavoidably results in the creation of the entangled field in the
convective shell with the thickness of $h\approx 250$ km in the case of a ONe
WD of 1.25 $M_{\odot}$ (Casanova et al. 2018). We neglect both the
amplification of the magnetic field by the convection and the field
dissipation at the small turbulent scales. The average magnetic field in the
convective shell is assumed to be equal to the average field on the WD surface
$B_{0}=10^{6}$ G.

Our scenario suggests that the nova shell forms in three stages: (i) the TNR
accompanied by the convection, (ii) the hydrogen shell expansion, and (iii)
the wind outflow. The envelope expansion at the stage (ii) from the WD radius
of $r_{0}=4\cdot 10^{8}$ cm up to the radius of $r_{1}=10^{10}$ cm results in
the weakening of the frozen-in magnetic field down to
$B_{1}=B_{0}(r_{0}h/r_{1}^{2})$ that is adopted as the initial field of the
wind at the level of its formation. The subsequent evolution of the field in
the wind proceeds differently for the radial ($B_{r}$) and tangential
($B_{t}$) components. The radial field decreases as
$B_{r}=B_{1}(r_{1}/r)^{2}$, whereas the tangential as
$B_{t}=B_{1}r_{1}^{2}/(r\Delta r)$, where $\Delta r=vt_{w}$ is the shell
thickness that remains constant. At the late time
$t\gg\mbox{max}(r_{1}/v,\Delta r/v)$ the entangled field $B$ in the shell is
dominated by the tangential components, so we adopt $B=B_{t}$. We neglect the
shell radial spread-out due to thermal velocities and ignore a possible shell
compression by the fast rarefied wind at the nebular stage of the nova
evolution ($t\gtrsim 100$ d).

###  2.2 Positron annihilation and escape

The 22Na isotope decays via the positron emission with the probability of 0.9
and the $e^{-}$ capture with the probability of 0.1. In both cases the
daughter nucleus 22Ne ends up in the excited state with a subsquent emission
of 1275 keV gamma-quantum (Firestone et al. 1999). Positron energies are
distributed in the range of 0 - 546 keV. Our description of positron
annihilation takes into account the numerical result that the probability of a
free annihilation is $p_{f}=0.1$ (Krannell et al. 1976; Leising and Clayton
1987), whereas 90% of positrons annihilate via the positronium formation at
thermal energies.

A positron slows down via ionization losses that depend on the matter
composition and the ionization degree. We take into account only hydrogen (X =
0.7) and helium and adopt that the hydrogen ionization fraction is 0.5, while
helium is neutral. In this case the fraction of free electron per barion is
$y_{f}=0.35$ and the fraction of bound electrons is $y_{b}=0.5$. The positron
slowdown is calculated according to the Bethe formula for ionization losses
that is splitted into terms related to neutrals and free electrons (Ahlen
1980)

$\frac{dE}{d\mu}=-\frac{4\pi
e^{4}N_{A}}{mc^{2}\beta^{2}}\left(y_{b}\ln{\frac{2mc^{2}\beta^{2}\gamma^{2}}{I}}+y_{f}\ln{\frac{2mc^{2}\beta^{2}\gamma^{2}}{\hbar\omega_{0}}}\right)\,,$
(1)

where $\mu$ is the surface density along the positron path (g cm-2), $e$ is
the elementary charge, $N_{A}$ is the Avogadro number, $m$ is the electron
mass, $c$ is the speed of light, $\beta=u/c$ ($u$ is the positron velocity),
$\gamma$ is the positron Lorentz factor, $I=29$ eV is the average ionization
potential for the adopted mixture of hydrogen and helium and the ionization
fraction, $\omega_{0}$ is the plasma frequency.

The positron range for a given energy is calculated taking into account the
energy distribution (Xiao et al. 2018). As an illustration, for the initial
energy of 200 keV and 400 keV the positron range at $t=100$ d ($\omega_{0}$
depends on the shell age) is $R=0.012$ g cm-2 and $R=0.035$ g cm-2,
respectively. These values differ from the average range of 0.1 g cm-2
estimated earlier (Leising & Clayton 1987). The difference presumably is
related to the adopted ionization degree. Indeed, we find for the neutral gas
the comparable value $R\sim 0.1$ g cm-2. The slowdown due to free electrons is
efficient because of the significant energy loss in Coulomb collisions.

The escape probability is determined by the ratio of the total surface density
along the random walk $\mu$ (g cm-2) in the shell and the range:
$p_{esc}=(1-p_{f})\exp{(-\mu/R)}$. The $\mu$ value is assumed to be a product
of the total geometrical length of the random walk $l$ and the density at the
radius of the positron creation. An alternative version, based on the average
density of the shell, produces the similar result. With the probability of
($1-p_{esc}$) the positron radiatively recombines with the branching ratio of
1/4 and 3/4 into the singlet and triplet states, respectively, in line with
statistical weights. The subsequent kinetics depends on the balance between
annihilation rates $t(^{1}S)=1.2\cdot 10^{-10}$ s, $t(^{3}S)=1.38\cdot
10^{-7}$ s (Karshenboim 2004), on the one hand, collisional transition rates
${}^{1}S\leavevmode\nobreak\ \rightleftarrows\leavevmode\nobreak\ ^{3}S$ and
the photoionization rate, on the other hand; these processes are taken into
account. The photoionization might affect the ratio between two- and three-
photon annihilation at the early epoch. However, already the day after the
TNR, the photoionization rate $P$ in the wind turns out already small enough,
$Pt(^{3}S)<0.1$, so at the subsequent stage the ratio between the two-photon
and three-photon annihilation does not deviate from 1/3.

A frozen-in entangled magnetic field impedes the positrons escape and this
should change the evolution of the annihilation line from the 22Na decay
compared to the case without magnetic field. The positron diffusion in the
entangled field is unlike the case of a uniform field. In the latter case the
diffusion across the magnetic field is strongly inhibited in the absence of
collisions and field perturbations, whereas along the field the positrons can
propagate freely. In the entagled field the diffusion accros the field can
occur because at some coherent length $l_{\mbox{\tiny B}}$ the field
significantly changes that results in the positron drift accross the field.
The diffusion in the entagled field looks therefore like an excursion along
the field on the scale $l_{\mbox{\tiny B}}$ and a drift accross the field.
Although there are approaches to the description of the diffusion in the
entangled field based on the scale $l_{\mbox{\tiny B}}$ and Larmor radius
(Narayan & Medvedev 2001; Chandran & Maron 2004) we do not see a possibility
to apply these approaches without three-dimensional modelling of the magnetic
field and the positron random walk.

Here we use the diffusion treatment in terms of the average Bohm diffusion
coefficient $D_{\mbox{\tiny B}}=(1/3)R_{\mbox{\tiny L}}u$, where
$R_{\mbox{\tiny L}}$ is the Larmor radius and $u$ is the positron velocity.
For the diffusion coefficient we adopt parametric representation $D=\xi
D_{\mbox{\tiny B}}$ with the coefficient $\xi\geq 1$. This approach with $\xi$
in the range of 1… 10 is adopted, e.g., for the diffusion in the chaotic field
when addressing the particle acceleration in supernova remnants (Marcowith et
al. 2006).

The characteristic diffusion time for the positron in the expanding shell with
the width $d$ is $t_{dif}\sim d^{2}/D$ while the total surface density
ignoring a possible annihilation is $\mu=\rho ut_{dif}$. This value combined
with the positron range $R$ due to ionization losses determines the avearge
escape probability of positrons $p_{esc}$ introduced above. In the adopted
model of the nova shell and the magnetic field evolution for the fixed initial
value at the WD surface the probability $p_{esc}$ is determined by the
parameter $\xi$ only. The adopted form of the diffusion coefficient is
equivalent to a random walk with the step $l$ being proportional to the shell
radius $R_{s}$. Indeed, in our model $l\propto R_{\mbox{\tiny L}}\propto
B^{-1}\propto R_{s}$. It is interesting that in this case the diffusion
coefficient formally corresponds to the approximation $\lambda=l/R_{s}=const$
used earlier in the description of the positrons diffusion in the entangled
field of SN Ia ejecta (Churazov & Khabibulin 2018).

##  3 Results

A standard model suggests the shell with the mass $M_{s}=10^{-5}$ $M_{\odot}$,
created by the wind with the velocity $v=1500$ km s-1 and the kinetic
luminosity $L_{w}=10^{38}$ erg s-1. The mass fraction of 22Na in the shell is
$10^{-3}$ and this isotope is assumed to be mixed in the shell homogeneously
with the constant fraction in any point. The wind presumably forms at the
radius $r_{1}=10^{10}$ cm; the result is weakly sensitive to the $r_{1}$
value. The average magnetic field at the WD surface is $B_{0}=10^{6}$ G.

The flux of directly escaped 511 keV quanta from the 22Na decay in nova at the
distance of 1 kpc is shown for $\xi=1$, 10, and 100, as well as for the case
without magnetic field (Fig. 1). The plot also shows the flux in the 1275 keV
line of 22Na. As expected, the magnetic field essentially changes the
evolution of the annihilation line, so its detection on the time scale of one
year becomes feasible. Remarkably, in the absence of the magnetic field the
flux in the 511 keV line does not reach the maximal possible value because of
the early positrons escape. For models with the magnetic field and $\xi=1$,
10, and 100 the fraction of escaped positrons is 0.41, 0.64, and 0.8,
respectively.

To remind, we do not take into account the magnetic field amplification by the
convection that can be significant. Indeed, the MHD 3D-modelling of the solar
convection demonstrates that a weak magnetic field can be amplified to the
value corresponding the Alfven velocity $V_{\tiny A}\sim 0.45u$ (Cataneo et
al. 2003). This relation suggests that the convection during the nova outburst
with the convective velocity $u\sim 100$ km s-1 (Casanova et al. 2018) in the
layer with the density of $\sim 3\cdot 10^{2}$ g cm-3 (Denissenkov et al.
2014) is able to amplify the field upto $\sim 10^{8}$ G. For $B_{0}>10^{6}$ G
positrons trapping will be even more effective.

The effect of other parameters — the shell mass, the shell velocity, and the
width of the convective shell — in the evolution of the annihilation line from
22Na decay is shown in Figure 2. As a fiducial model we take the standard
model with $\xi=10$. The twice as large shell mass increases the flux and the
time of the flux maximum. The lower expansion velocity by a factor of 1.5
increases the positron trapping and also decreases the time of the flux
maximum. The increase of the width of the convective shell results in the
larger field at the late time and therefore more efficient trapping. We do not
demonstrate the case with the kinetic wind luminosity of $2\cdot 10^{38}$ erg
s-1 because the only effect of this variation is the decrease of the maximum
age. Versions with the different magnetic field are not shown as well because
effect of smaller/larger value of the initial field is similar to the effect
of larger/smaller value of the $\xi$ parameter by the same factor. This
suggests that the effect of the field and the diffusion coefficient is
actually determined by the combined parameter $B_{0}/\xi$.

##  4 Discussion and conclusions

The paper has been aimed at the study of effects of the magnetic field in the
nova shell on the escape of positrons from the 22Na decay and on the flux
evolution in the annihilation line of 511 keV. We present the model for the
positron annihilation in a nova shell that takes into account a magnetic
field. The modelling shows that for the sensible value of the average magnetic
field on the WD surface of $10^{6}$ G, the frozen-in magnetic field is able to
essentially impede the positrons escape and thus to prolong the luminosity of
the annihilation line from 22Na up to hundreds days.

Results, however, depend not only on the strength of magnetic field, but on
the diffusion coefficient as well. The dimensionless diffusion coefficient
$\xi$ used above is in the range $1\leq\xi\leq 100$, that includes the range
of $1\leq\xi\leq 10$ in models of the diffusive shock particle acceleration
constrained by the non-thermal X-ray brightness at the limb of young supernova
remnants (Marcowith et al. 2006).

Yet there remains an open question: to which extent an assumption of the
homogeneous chaotic magnetic field is applicable to nova shells. It may well
be that the turbulent convection results in the intermittent chaotic field
with regions of strong and weak field. If this the case and the filling factor
of regions with the weak field is significant, the diffusion and escape of
positrons could be more efficient compared to the homogeneous case. Partially
this situation is taken into account by the large value of $\xi=100$. However,
actually this does not resolve the issue of a possible intermittency of the
chaotic field and the problem requires separate study. It is noteworthy that a
simultaneous detection of 1275 keV and 511 keV lines will become a valuable
tool for the diagnostics not only of the synthesised 22Na mass but also of the
parameter $B_{0}/\xi$ that determines the positron escape.

Shells of historical novae have clumpy structure, e.g., nova 1970 FH Ser (Gill
& O’Brien 2000), and this, at first glance, is in conflict with our model of
the spherical shell. However, the effect of the clumpy structure on the 22Na
positron escape depends on the stage at which this clumpiness does emerge.
There is no reason for the clumpiness to appear at the stage of the optically
thick wind. More naturally to admit that the shell experinces a fragmentation
significantly later under the dynamical effect of a high velocity rarefied
wind (Cassatella et al. 2004). The high velocity wind is formed by the remnant
of the hot envelope of the WD due to the radiative acceleration. A signature
of the wind with the velocity of $\sim 3000$ km s-1, versus the main shell
velocity of $\sim 1800$ km s-1, has been observed in nova V1974 Cyg
(Cassatella et al. 2004). Based on the WR-star wind one expects that the
kinetic luminosity of the high-velocity wind ($L_{w,f}$) should be $\sim
10^{-2}$ of the radiative luminosity of the WD (Sander et al. 2020), i.e.,
$L_{w,f}\sim 10^{36}$ erg s-1 in our case. One can estimate the fragmentation
time assuming the fast wind velocity of 3000 km s-1. The major fragmentation
mechanism is the Rayleigh-Taylor instability initiated by the shell
acceleration driven by the fast wind. The fragmentation time scale exceeds the
initial phase of the acceleration determined by the shell compression due to
the shock propagation. For the shell mass of $M_{s}=10^{-5}$ $M_{\odot}$ and
the thickness $d=M_{s}/w$ formed by the wind with the kinetic luminosity
$L_{w}=10^{38}$ erg s-1 the compression by the fast wind takes the time
$t_{c}\sim d/D$, where $D$ is the shock wave speed determined from the
momentum conservation $D=(\rho_{f}/\rho)^{1/2}(v_{f}-v)=53$ km s-1, $\rho_{f}$
is the density of the fast wind, $v_{f}$ is its velocity, while the rest of
values are related to the nova shell. In the considered case $t_{c}\sim 2$
years. This time refers to the initial stage of the nova shell fragmentation
that requires somewhat longer time, of the order of $t_{f}\sim 4t_{c}$ (Klein
et al. 1994), i.e., about 8 years. The found fragmentation time indicates that
the approximation of the spherical shell on the time scale of one year
generally does not contradict to the observed clumpy shells of historical
novae.

Although novae have been considered as a possible source of Galactic
positrons, in effect their contribution to the Galactic 511 keV gamma-line is
unlikely significant even for 100% positron escape. With the Galactic nova
rate of 50 yr-1, neon novae fraction of 1/4, nova shell mass of $10^{-5}$
$M_{\odot}$, and 22Na mass fraction of $10^{-3}$ the Galactic production rate
of positrons by novae is only of $2\cdot 10^{41}$ s-1 that is by a two order
lower compared to the Galactic annihilation rate of $5\cdot 10^{43}$ s-1
(Siegert et al. 2016).

It remains unclear, whether WD of novae are always magnetic with the field
$\gtrsim 10^{6}$ G? The recent study of the cataclismic variables from the
Gaia DR2 within 150 pc (Pala et al. 2020) shows that 36% of them posess
megagauss fields. Given the fact that the field of $\sim 10^{6}$ G is a border
value for the cataclismic variables to be classified as an intermediate polar
with the accretion on the poles, a significant fraction of WD in cataclismic
variables could have magnetic fields $<10^{6}$ G and thus not show signatures
of the intermediate polar. Yet in this case the magnetic field can be
nevertheless strong enough to impede 22Na positrons escape. For example, for
the field of $10^{5}$ G and $\xi=10$ the situation with the trapping of 22Na
positrons in effect is equivalent to the case of $B_{0}=10^{6}$ G É $\xi=100$
(Fig. 2) that demonstrates the efficient positron annihilation on the time
scale of $\sim 100$ days.

One of us (N.C.) is grateful to Eugene Churazov and Alexander Getling for
stimulating discussions.

References

Ahlen S.P., Reviews of Modern Physics 52, 121 (1980)
Casanova J., José J., Shore S., Astron. Astrophys. 619, 121 (2018)
Cassatella A., Lamers H. J. G. L. M., Rossi C., et al., Astron. Astrophys.
420, 571 (2004)
Cattaneo F., Emonet T., Weiss N., Astrophys. J. 588, 1183 (2003)
Chandran B.D.G., Maron J.L., Astrophys. J. 602, 170 (2004)
Chomiuk L., Metzger B.D., Shen K.J.), arXiv201108751 (2020)
Churazov E., Khabibulin I., Mon. Not. R. Astron. Soc. 480, 1394 (2018)
Clayton D.D., Hoyle F., Astrophys. J. 187, L101 (1974)
Crannell C.J., Joyce G., Ramaty R., Astrophys. J. 210, 582 (1976)
Denissenkov P.A., Truran J. W., Pignatari M., et al., Mon. Not. R. Astron.
Soc. 442, 2058 (2014)
Epelstain N., Yaron O., Kovetz A., Prialnik D., Mon. Not. R. Astron. Soc. 374,
1449 (2007)
Ferland G.J., Shields G.A., Astrophys. J. 226, 172 (1978)
Ferrario L., Wickramasinghe D., Kawka A., Advances Space Res. 66, 1025 (2020)
Firestone R.B., Shirley V.S., Baglin C.M. et al., Table of Isotopes (John
Wiley and Sons, New York, 1999)
Gill C.D., O’Brien T.J., Mon. Not. R. Astron. Soc. 314, 175 (2000)
Gomez-Gomar J., Hernanz M., Jose J., Isern J., Mon. Not. R. Astron. Soc. 296,
913 (1998)
Iyudin A.F., Astron. Rep. 54, 611 (2010)
Iyudin A.F., Bennett K., Bloemen H., et al., Astron. Astrophy. 300, 422 (1995)
Karshenboim S.G., International Journal of Modern Physics A 19, 3879 (2004)
Kato M., Hachisu I., Astrophys. J. 657, 1004 (2007)
Klein R.I., McKee C.F., Colella P., Astrophys. J. 420, 213 (1994)
Kraft R.P., Astrophys. J. 139, 457 (1964)
Kudryashov A.D., INASAN Science Reports 3, 205 (2019)
Leising M.D., Clayton D.D., Astrophys. J. 323, 159 (1987)
Marcowith A., Lemoine M., Pelletier G.), Astron. Astrophys. 453, 193 (2006)
Mustel E.R, Baranova L.I., Soviet Astron. 10, 388 (1966)
Narayan R., Medvedev M.V., Astrophys. J. 562, L129 (2001)
Pala A.F., Gänsicke B.T., Breedt E., et al., Mon. Not. R. Astron. Soc. 494,
3799 (2020)
Patterson J., Publ. Astron. Soc. Pacific 106, 209 (1994)
Ruggles C.L.N., Bath G.T., Astron. Astrophys. 80, 97 (1979)
Sander A.A.C., Vink J.S., Hamann W.-R., Mon. Not. R. Astron. Soc. 491, 2206
(2020)
Shara M.M., Prialnik D., Astron. J. 107, 1542 (1994)
Siegert T., Diehl R., Vincent A.C., et al., Astron. Astrophys. 595, 25 (2016)
Sparks W.M., Starrfield S., Truran J.W., Astrophys. J. 208, 819 (1976)
Starrfield S., Bose M., Iliadis C. et al., Astrophys. J. 895, 70 (2020)
Starrfield S., Sparks W.M., Truran J.W., Astrophys. J. 303, 186 (1986)
Takeda L., Diaz M., Publ. Astron. Soc. Pacific 131, 54205 (2019)
Xiao H., Hajdas W., Wu B., et al., Astroparticle Physics, 103, 74 (2018)

Figure 1:  Flux of gamma-quanta in 1275 keV line (thin line) and in the
annihilation line 511 keV from the 22Na decay in ONe nova shell at the
distance of 1 kpc. Numbers next to lines indicate the value of $\xi$. Dotted
line shows the model without magnetic field.

Figure 2:  Flux in the annihilation line 511 keV from the 22Na decay in the
standard model with $\xi=10$ (thick line) and in models with the variation of
a certain parameter (shell mass, velocity, and thickness of the convective
shell on the WD surface).


# Exploring the Mechanisms Behind Non Aromatic Fluorescence with Density
Functional Tight Binding Method

Gonzalo Díaz Mirón  [  Carlos R. Lien-Medrano  [  Debarshi Banerjee  [  Uriel
N. Morzan  [  Michael A. Sentef  [  Ralph Gebauer  [  Ali Hassanali  [
<EMAIL_ADDRESS>
###### Abstract

Recent experimental findings reveal non-conventional fluorescence emission in
systems devoid of conjugated bonds or aromatic compounds, termed Non-Aromatic
Fluorescence (NAF). This phenomenon is exclusive to aggregated or solid
states, remaining absent in monomeric solutions. Previous studies on small
system models in vacuum emphasize that the carbonyl stretching mode, along
with the strong interaction of Short Hydrogen Bonds (SHBs), remains the
primary vibrational mode explaining NAF in these systems. This work
investigates NAF in the crystal structure of L-pyroglutamine-ammonium,
comparing it with the related non-fluorescent amino acid L-glutamine.
Employing Density Functional Tight Binding (DFTB) theory in combination with
Non-Adiabatic Molecular Dynamics (NAMD) and the Mixed Quantum/Molecular
Mechanics (QM/MM) approach, we simulate larger model systems and consider the
effects of the surrounding environment. Our results on more realistic systems
align with previous reports using higher-level theories. Moreover, our study
demonstrates the efficacy of DFTB in elucidating the intricacies inherent to
NAF, underscoring the qualitative insight into this novel optical property.
Beyond contributing to the ongoing exploration of NAF, this work forms the
basis for delving further into the complexities characterizing non-aromatic
systems within the fluorescence paradigm.

###### keywords:

Non-aromatic fluorescence, Density Functional Tight-Binding, Non adiabatic
molecular dynamics

ICTP]Condensed Matter and Statistical Physics, The Abdus Salam International
Centre for Theoretical Physics, Trieste, Italy Bremen]Institute for
Theoretical Physics and Bremen Center for Computational Materials Science,
University of Bremen, 28359 Bremen, Germany ICTP]Condensed Matter and
Statistical Physics, The Abdus Salam International Centre for Theoretical
Physics, Trieste, Italy DFUBA]Instituto de Fisica de Buenos Aires, Facultad de
Ciencias Exactas y Naturales, Universidad de Buenos Aires, Argentina
Bremen]Institute for Theoretical Physics and Bremen Center for Computational
Materials Science, University of Bremen, 28359 Bremen, Germany [MPSD] Max
Planck Institute for the Structure and Dynamics of Matter, Center for Free-
Electron Laser Science (CFEL), Luruper Chaussee 149, 22761 Hamburg, Germany
ICTP]Condensed Matter and Statistical Physics, The Abdus Salam International
Centre for Theoretical Physics, Trieste, Italy ICTP]Condensed Matter and
Statistical Physics, The Abdus Salam International Centre for Theoretical
Physics, Trieste, Italy NAF,DFTB,NAMD

##  1 Introduction

Over the last decade, a growing body of experimental evidence has revealed the
emergence of non-conventional fluorescence signals in biomolecular systems. In
contrast to our standard spectroscopic intuition attributing fluorescence in
biological systems to aromatic or conjugated groups, recent intriguing
experiments have shown that non-aromatic systems exhibit inherent absorption
and fluorescence in the UV-Visible range1, 2, 3, 4. These fluorescent signals
appear to be enhanced in aggregated forms, analogous to the aggregation-
induced emission observed in organic materials5, 6 Noteworthy examples include
amyloid structures7, 8, 9, individual amino acids10, 11, compounds derived
from amino acids12, 10, poly-amidoamines13 and small organic molecules14, 15.

While the detailed molecular mechanism underlying Non-Aromatic Fluorescence
(NAF) remains an open question, several possible explanations have been
proposed1, 5. These include electron delocalization along peptide bonds in
ordered secondary structures3, the constriction of the carbonyl stretching
mode caused by intense local interactions such as strong hydrogen bonds16, 12,
proton transfer along short hydrogen bonds17 and inter-residue charge-transfer
excitations followed by charge recombination4. These molecular phenomena are
believed to introduce low lying fluorescent electronic states in the visible
regime arising from the supramolecular assemblies that do not rely on the
presence of $\pi$-delocalized electrons. Since fluorescence emission involves
the de-activation from electronic excited states, it requires a complex
interplay involving the dynamics of both nuclear and electronic degrees of
freedom which are notoriously difficult to model.

Over the last few years, our group has been leading theoretical and
computational efforts to investigate the spectroscopic origins of this
phenomenon9, 16, 12, 18, 1. We have recently investigated NAF in model amyloid
proteins9, 18 and L-glutamine/L-glutamine-derived16, 12 amino-acids combining
both experiments and theory. Employing Trajectory Surface Hopping (TSH)19, 20,
simulations at a Time Dependent Density Functional Theory (TD-DFT) level21, we
have shown that carbonyl groups (CO) appear to play a critical role in
modulating the observation of NAF in a wide variety of systems.

Upon an initial photo-excitation to the first excited state $S_{1}$, which is
mainly localized around the CO group, strong inter-molecular interactions
involving hydrogen bonds within certain biological aggregates hinder the
stretching of the CO bond. This leads to a restriction of the non-radiative
relaxation pathways, subsequently boosting the fluorescence quantum yield. The
CO stretching mode is, at the same time, strongly coupled to other vibrational
modes involving the amide bond which appear to be modulated by local
environmental effects. Our previous research has been based on dimer model
systems carved out of crystal structures; in these studies we have accounted
for the steric effect of the surrounding molecules by incorporating specific
constraints in our simulations12, 9, 16. The computational demands associated
with simulating the photophysics of these systems make TSH calculations
extremely intensive. Furthermore, the inherent variability in the electronic
relaxation pathways can introduce slow convergence problems. Consequently, a
substantial number of trajectories becomes essential to derive meaningful and
reliable results. We have recently employed hybrid Quantum-Classical
Mechanical (QM/MM) approaches to enable a realistic description of the
environment effects on the fluorescence.16 However, the computationally
tractable sizes of the quantum mechanical region are severely limited. These
limitations have sparked our interest in exploring alternative methods to
effectively tackle these challenges. In this regard, Density Functional Tight-
Binding (DFTB)22, 23, an approximate method based on Density Functional Theory
(DFT), has emerged as a promising approach. DFTB achieves computational
efficiency through a combination of the two-center approximation, pre-
tabulated integrals, and the use of a minimal basis set24, 25, 26. This
accelerates DFTB calculations by a factor of 2-3 compared to traditional
DFT25. Moreover, the pre-tabulation of integrals at DFT level ensures the
precision of DFTB in predicting electronic properties for molecular, nano, and
periodic systems, effectively emulating DFT27, 28, 29, 30, 31, 32. This
framework has also been extended to include the modeling of excited state
dynamics leading to the Time-Dependent Density Functional Tight-Binding method
(TD-DFTB). These extensions include both the frequency33, 34(TD-DFTB) and the
time domain35, 36 (real-time TD-DFTB). In addition, TD-DFTB has been coupled
to the Ehrenfest approach35, 36 and TSH37, 38, 39 to perform non-adiabatic
dynamics.

The goal of this work is to comprehensively assess the potential of the TD-
DFTB method in investigating the photophysics underlying NAF phenomenon
employing the DFTB+ code40. The primary focus of our investigation is on
amino-acid crystals which have recently been studied in our group using both
ab-initio electronic methods (TD-DFT) and experiments. It was found that
L-glutamine, characterized by regular hydrogen bond interactions, is a non-
fluorescent system. However, following a chemical transformation into L-pyro-
ammonium (pyroglutamic acid complexed with an ammonium ion), the system
undergoes a transition to a short hydrogen bond configuration and displays a
remarkable fluorescent signal. (Panels a and b of Figure 1 show the molecular
structure of these systems.). By comparing and contrasting the mechanisms
obtained with TD-DFTB to our previous work, we demonstrate that the DFTB
approach correctly captures the non-radiative decay pathways occurring in
L-glutamine as well as the origins of fluorescence in L-pyro-ammonium. With
this validation in hand, we use TD-DFTB to understand the role of
environmental effects using a QM/MM approach allowing for an accurate
characterization of the experimentally observed Stoke-Shift between the
excitation and emission spectra in this system. We also highlight aspects of
non-radiative decay mechanisms which are different from our previous studies
paving the way for future directions and developments needed to use Tight-
Binding approaches to study non-aromatic fluorescence.

##  2 Theoretical Approaches

The objective of this section is to provide a background on the various
theories and techniques employed in our methodologies to study non-aromatic
fluorescence. In this section, we will summarize some of the most important
aspects of DFTB, as well as its Time-Dependent extension to the frequency
domain (TD-DFTB), with a special focus on the implementations made in the
DFTB+ code40. Additionally, we will describe the key aspects of non-adiabatic
dynamics in the context of Trajectory Surface Hopping (TSH), with a specific
emphasis on the different analytical expressions used to calculate the non-
radiative probability decay.

###  2.1 Density Functional Tight-Binding (DFTB and TD-DFTB)

DFTB equations are derived from DFT by expansion of the total energy in a
Taylor series of the electron density fluctuations $\delta\rho$ around a
reference density $\rho^{0}$40. The chosen reference density is usually a
summation of overlapping, spherical and noninteracting atomic charge
densities. Up to the third order, the total energy can be approximated as:41

$E^{DFT}[\rho^{0}+\delta\rho]\approx
E^{0}[\rho^{0}]+E^{1\text{st}}[\rho^{0},\delta\rho]+E^{2\text{nd}}[\rho^{0},(\delta\rho)^{2}]+E^{3\text{rd}}[\rho^{0},(\delta\rho)^{3}].$
(1)

The initial term $E^{0}$ can be represented by summing pairwise repulsive
potential energies $E^{rep}_{AB}$. The truncation up to first order in the
energy can be expressed as

$E^{DFTB1}=E^{0}+E^{1\text{st}}=\sum_{A>B}E_{AB}^{rep}+\sum^{occ.}_{i}n_{i}\langle\psi_{i}|\mathcal{H}[\rho^{0}]|\psi_{i}\rangle,$
(2)

where $\langle\psi_{i}|$ is the $i^{\text{th}}$ molecular orbital in a linear
combination of atomic orbitals (LCAO), $n_{i}$ is the occupation and
$\mathcal{H}[\rho^{0}]$ is the Hamiltonian operator in a two-center
approximation.

In the second-order DFTB (DFTB2) method, the $E^{2nd}$ term is incorporated
through a monopole approximation23, where the electron density fluctuations
are approximated as a sum of atomic Slater-type spherical charge densities:

$E^{DFTB2}=E^{DFTB1}+\frac{1}{2}\sum_{AB}\gamma_{AB}\Delta q_{A}\Delta q_{B}.$
(3)

Here, $\Delta q_{A}$ is the Mulliken charge on atom A and $\gamma_{AB}$
represents the electron interaction of two Slater-type spherical charge
densities on atoms A and B. The inclusion of the second term requires a self-
consistent solution since the Mulliken charges depend on the molecular
orbitals. This provides a more accurate description of charge-transfer and
charge-charge interactions.

DFTB3 allows the chemical hardness of an atom to change with its charge states
by including third order terms:

$\displaystyle E^{DFTB3}=E^{DFTB2}+\frac{1}{3}\sum_{AB}\Gamma_{AB}\Delta
q_{A}^{2}\Delta q_{B},$ (4)

where $\Gamma_{AB}$ is the derivative of $\gamma_{AB}$ with respect to the
$\Delta q_{A}$.

Corrections for both van-der-Waals42, 43 and hydrogen-bonding interactions44,
45 were also implemented in the method owing to the fact that these types of
weak interactions are not accurately captured by DFT based methods including
DFTB.

Within the DFTB+40 code there are several available approaches to determine
excited states. We employed the Casida formalism21 which within the framework
of DFTB, follows the same ansatz as in the original version of TD-DFT34, 33.
The electronic excitation energy $\omega_{I}$ in TD-DFTB can be obtained by
solving the following eigenvalue problem:

$\displaystyle\Omega\mathbf{F}_{I}=\omega^{2}_{I}\mathbf{F}_{I}.$ (5)

Here, $\Omega$ is the response matrix and its elements are

$\displaystyle\Omega_{ij\sigma,kl\tau}=\delta_{ik}\delta_{jl}\delta_{\sigma\tau}\omega^{2}_{ij}+2\sqrt{\omega_{ij}}K_{ij\sigma,kl\tau}\sqrt{\omega_{kl}},$
(6)

where $\omega_{ij}=\epsilon_{j}-\epsilon_{i}$, $i,k$ and $j,l$ are occupied
and unoccupied KS orbitals respectively, $\sigma$ and $\tau$ are spin indices,
whereas $K$ is the so-called coupling matrix. Within the monopole
approximation, the coupling matrix elements $K_{ij\sigma,kl\tau}$ adopt simple
expressions, reducing the computational cost.33

Moreover, the incorporation of long-range corrections into the two-electron
integrals have also been implemented in the DFTB+ code46 with the aim to
provide a reduction in the self-interaction error in DFT (and also present in
DFTB) yielding a correct behaviour of the asymptotic potential. This
implementation along with a new reparametrization of the Slater-Koster
parameters enables a more accurate description of the excited states, mainly
for those which involves charge transfer excitations.47

###  2.2 Trajectory Surface Hopping (TSH)

Trajectory Surface Hopping (TSH) is one of the most widely employed techniques
to model non-adiabatic dynamics20. The non-radiative relaxation is a non-
equilibrium process and therefore the dynamics strongly depends on the initial
conditions. For this reason, the simulation of a large number of thermally
distributed independent trajectories, with varying initial conditions becomes
essential48, 49. In each TSH trajectory, the electronic states are propagated
using quantum mechanics, while the nuclear motion is handled classically
employing the forces coming from a single potential energy surface from the
quantum method. For more information on the theoretical and practical details,
the reader is referred to previous studies20, 50. Bellow, we highlight
specific details relevant to our present work.

The TSH workflow involves propagating the system on distinct potential energy
surfaces. At each timestep, the non-radiative probability is computed and a
stochastic algorithm is employed to decide which potential energy surface the
system will proceed along. The most common expression for the transition
probability is the fewest switches surface hopping (FSSH) proposed by Tully19:

$P_{ij}(t)=-2\int_{t}^{t+\delta
t}dt^{\prime}\dfrac{C_{i}(t^{\prime})C_{j}(t^{\prime})^{*}g_{ij}(t^{\prime})}{C_{i}(t^{\prime})C_{i}(t^{\prime})^{*}}$
(7)

The $g_{ij}$ terms are the non adiabatic coupling (NAC) elements, one of the
key components in non-adiabatic dynamics. Several alternatives approaches have
been developed in order to avoid the explicit determination of the non
adiabatic coupling elements51, 52, 53, since sometimes there are no analytical
expressions available or because they are computationally expensive to
determine for large systems. One of the most popular in this context is the
Landau-Zener (LZSH) approximation54:

$P_{ij}(t)=exp\left(-\dfrac{\pi}{2\hbar}\sqrt{\dfrac{\Delta
E_{ij}^{3}(t)}{\dfrac{d^{2}}{dt^{2}}\Delta E_{ij}(t)}}\right)$ (8)

where $\Delta E_{ij}$ is the adiabatic energy gap55, 50. Within the context of
our TD-DFTB simulations, we employ the LZSH approximation instead of Tully’s
method due to its easy implementation into our computational framework.

##  3 Computational Details

In this section, we summarize the computational protocols applied to
L-glutamine (L-gln) and L-pyroglutamine ammonium (L-pyro-amm) to simulate the
ground-state molecular dynamics in periodic boundary conditions, the
determination of the absorption spectra, the excited states surface hopping
dynamics and finally the QM/MM simulations with DFTB/TD-DFTB. In the present
study, we employed the DFTB+ package40 for the electronic structure
calculations. The Genesis56, 57 software was used to perform classical and
QM/MM simulations.

###  3.1 Ground State DFTB Molecular Dynamics

Ab-initio molecular dynamics simulations of both L-gln and L-pyro-amm were
conducted using DFTB in the ground-state in periodic boundaries conditions
(PBC) using the gamma point. The initial coordinates for both systems were
obtained from previously resolved crystallographic structures (CSD 116931658
and 198155112). We employed third order corrections in the Hamiltonian with
the corresponding Slater Koster parameters 3ob59 as well as D3 Grimme
corrections42. Molecular dynamics simulations were conducted for a total of
600 ps for each system using a timestep of 0.5 fs and within the NVT ensemble
using the Berendsen thermostat with a coupling of 0.1 $ps^{-1}$. For all the
analysis, the first 20 ps were discarded as equilibration. As will be shown
later, we compare some of the ground-state structural properties obtained from
DFTB to those obtained with DFT. Details of these simulations can be found in
the Computational Methods Section of a previous work in our group12.

###  3.2 Absorption Spectra

After conducting the DFTB simulation in the ground state for L-gln and L-pyro-
amm, we proceeded to extract frames at intervals of 500 fs. A total of 1000
frames were used to compute the absorption spectra for each system with the
aim to allow for adequate structural fluctuations in the crystal and their
subsequent coupling to electronic transitions.

These calculations were carried out within the TD-DFTB method. To account for
the correct description of the excited states, we also incorporated Long Range
Corrections (LC) with the appropriate Slater-Koster files47. A total of 100
states were calculated for each snapshot employing the Casida algorithm. The
final absorption spectrum for each system was generated by averaging the
frames, and each individual band was broadened using a Gaussian function with
a Full Width at Half Maximum (FWHM) of 5 nm.

###  3.3 Excited State Molecular Dynamics in Gas Phase

TSH simulations were performed on dimer model systems carved out from the
crystal structures of L-gln and L-pyro-amm. The dimer models were carefully
selected to preserve the Short Hydrogen Bonds (SHBs) in L-pyro-amm and the
normal Hydrogen Bonds (HBs) in L-gln. These simulations were conducted with
the dimer in the gas phase. To account for the crystallographic environment,
we applied atom constraints on specific selected atoms (see Figure S1 in the
Supplementary Information).

Our protocol to generate the initial conditions consists of performing 6 ns
long classical molecular dynamics simulations using the Generalized Amber
Force Field (GAFF)60 for both model systems with a timestep of 1 fs and a
thermostat constant of 1 $ps^{-1}$. From these long simulations, a total of
600 configurations were selected with a spacing of 10 ps. For each of these
600 configurations, ground-state simulations were continued using DFTB as
described earlier in section 3.1, for 2 ps each at 300K.

TSH simulations were started with the system in the electronic state $S_{1}$
corresponding to the initial photo-absorption. The dynamics were then
performed in the NVE ensemble for 250 fs using a timestep of 0.5 fs with TD-
DFTB as described in section 3.2. At each timestep, the probability of non-
radiative decay was calculated using the LZSH algorithm (Equation 8).

As noted earlier in the Introduction, one of the goals of this work is to
compare and contrast the underlying non-radiative decay mechanisms associated
with non-aromatic fluorescence obtained with TD-DFT in our earlier work16 to
TD-DFTB. It is therefore important to recall two important differences.
Firstly, the use of TD-DFTB expands the statistics of TSH trajectories by a
factor of 3 (600 vs 200) in the case of L-gln and a factor of 6 for L-pyro-amm
(600 vs 100). Secondly, as already indicated earlier, the TD-TBDFT hopping
events between different states are performed using LZSH (Equation 8) while in
our previous TD-DFT simulations, were conducted with FSSH (Equation 7).
Despite this methodological difference, the comparison between the two
approaches remains valid, as both theories have been extensively tested and
compared in several works61, 62, 63, demonstrating consistent results.

###  3.4 QM/MM Simulations

In addition to simulations of the dimer in the gas phase, we also conducted
the QM/MM simulations in order to understand the role of the environment64,
65. This method involves partitioning the system into two distinct regions,
namely a QM and MM, and describing the coupling between these regions using
the electrostatic embedding approach66, 65, 67. To investigate the radiative
and non-radiative process in both systems, we selected two different QM sub-
systems for L-pyro-amm, the dimer and unit cell, while for L-gln we studied
the unit cell.

To construct the system, we began by replicating the unit cell of the crystal
structure in all three directions resulting in a supercell comprising 375 unit
cells for both L-gln and L-pyro-amm. The subsequent step involved a classical
equilibration at room temperature for 10 ns with a timestep of 1 fs. Classical
parameters for bonded and Lennard-Jones were adopted from GAFF60 and for
treating the electrostatics we used Mulliken charges which were extracted from
optimized crystal structures using DFTB in PBC (see section 3.1). A total of
10 snapshots were extracted from the simulations.

After selecting the QM region, QM/MM ground state dynamics simulation of 25 ps
each were conducted with the same simulation protocol for each snapshot. From
the last 20 ps of these simulations, we extracted 20 frames for each
trajectory to calculate the absorption spectra using TD-DFTB in the crystal
environment (see section 3.2). A total of 200 snapshots were utilized for the
absorption spectra.

In addition, we performed TSH simulations for each of the QM sub-system with
the aim to study the emission as well as the photo-physical relaxation of the
$S_{1}$ state of these systems. TSH simulations started from the $S_{1}$
electronic state employing the initial conditions generated in the previous
step for a total of 1 ps, with a timestep of 0.5 fs within the NVE ensemble. A
total of 10 TSH dynamics were conducted for each system using the same
methodology described section 3.3.

##  4 Results and Discussion

###  4.1 Ground State Properties in Molecular Crystals

Before analyzing the optical properties of these non-aromatic systems with TD-
DFTB, it is necessary first to verify the quality of the hydrogen bonds (HBs)
description, in particular the difference in the distance and its strength for
the L-Pyro-Ammonium (L-pyro-amm) and L-Glutamaine (L-gln) in the crystal
structures, since this property plays an important role in the fluorescence12,
9. Figure 1 shows the probability distribution function (PDF) obtained from
ab-initio ground state molecular dynamics with PBC (see section 3.1) for the
typical descriptors employed to describe the HBs including the proton transfer
coordinate (PTC) and bond distance (BD) between the two heavy atoms among
which proton transfer potentially occurs. BD is defined as the distance
between O and N atoms in the case of L-gln (panel a and c) and the distance
between the two O atoms in the case of L-pyro-amm (panel b and d). The PTC is
defined as the distance difference between the O-H and N-H in L-gln (panel e)
and between the two O-H involved in L-pyro-amm (panel f).

Figure 1: Unit cell crystallographic geometries for L-gln (panel a) and
L-pyro-amm (panel b), where the blue circles highlight the HBs and SHBs,
respectively. Histograms of BD (panels c and d) and PTC (panels e and f)
obtained from an ab-initio ground state dynamics with PBC using DFT (blue
line) and DFTB (red lines)

Overall, we observe that DFTB works very well at reproducing the properties of
the HBs in both L-gln and L-pyro-amm. Specifically, the bond lengths
associated with the SHBs are shorter than those of normal HBs comparing the
two systems ($2.4-2.6$ Åvs $2.7-3.0$ Å). Another important feature of the SHBs
in L-pyro-amm was the proton transfer, which has been observed in our previous
DFT simulations12. This is reflected in a double-well potential along the PTC.
Panels e and f show that DFTB accurately reproduces this feature for the SHBs
and the single asymmetric well of the normal HBs. Given the nature of the DFTB
method along with its various approximations, we find the agreement very
encouraging.

###  4.2 Absorption Spectra in Molecular Crystals

With the ground-state configurations obtained from the ab-initio MD, we next
turn to extracting the absorption spectra. Figure 2 shows the absorption
spectra for L-gln and L-pyro-amm obtained experimentally12 (left panel) as
well as the calculated with TD-DFTB (right panel). As one can observe, an
important difference between the two systems is the extent of the red-
tail/edge excitation that is significantly more pronounced in the case of
L-pyro-amm. In the case of L-gln, the absorption peaks are below the 250 nm
and outside the window of experimentally detectable wavelengths. On other hand
for L-pyro-amm, there is a very pronounced red-edge excitation tail ranging
from 250-400 nm.

Figure 2: Absorption UV-Vis Spectra for L-Glutamine and L-Pyro-Ammonium. Left
panel shows the experimental data, extracted from reference 12. Right panel
shows the calculated average spectra using TD-DFTB. Red lines correspond to
L-pyro-amm system while blue lines for L-gln.

The right panel of Figure 2 depicts the average calculated absorption spectra
obtained for both systems using 1000 snapshots each, at the level of TD-DFTB.
Interestingly, we observe that TD-DFTB reproduces the correct experimental
trends comparing L-gln and L-pyro-amm. Analysing the transitions involved in
the UV-Vis energy range (330-370 nm) we found that the main molecular orbitals
that explain these transitions involve significant contributions coming from
the carbonyl bonds and the SHB (see Figure S2 in the Supplementary
Information). It is also worth noting that this prediction is also fully
consistent with the results that were previously obtained at TD-DFT level in
our previous work12.

###  4.3 Excited State TSH Simulations in Gas phase

The two molecular systems (L-gln and L-pyro-amm) display substantially
different photophysical properties. Specifically in L-pyro-amm, non-radiative
decay pathways are effectively shut off16 thereby enhancing the excited state
lifetime and subsequently increasing the fluorescence yield. From the TSH
simulations on both systems, we extracted the average populations of the
$S_{1}$ and $S_{0}$ electronic states over time which are illustrated in
Figure 3 for the dimer systems in vacuum. The solid lines come from the TD-
DFTB simulations in this work, while the dotted lines are obtained from our
previous TD-DFT calculations16. From the evolution of the electronic states,
we can observe that TD-DFTB captures the essential photophysics of both
systems. In L-gln dimers, an ultra-fast non-radiative decay is evident, with a
characteristic lifetime of $\sim$90 fs for TD-DFT and $\sim$50 fs with TD-
DFTB. Conversely, L-pyro-amm exhibits a much longer $S_{1}$ lifetime with both
theory levels. These findings are again, fully consistent with the
experimental findings12 where fluorescence is only observed in the L-pyro-amm
molecular crystal.

Figure 3: Average populations of the electronic states obtained from TSH
trajectories for L-gln (left panel) and L-pyro-amm (right panel) dimer systems
in vacuum. Results using TD-DFT are depicted as dotted lines while those
obtained with TD-DFTB in solid lines. The evolution of $S_{1}$ and $S_{0}$ are
denoted in red and blue lines, respectively.

In our previous studies employing TD-DFT, we observed that the transition from
$S_{1}\to S_{0}$ occurs via a conical intersection (CoIn) where the essential
vibrational modes coupled to the process involves the extension of the
carbonyl bond, proton transfer along the SHB, and finally, the deplanarization
of the amide group. Examining each of these modes in the ground and excited
state as well as near the CoIn, a very specific structural signature becomes
apparent. Figure 4 depicts the histograms for each of these modes obtained
from TSH simulations employing both TD-DFT (top) and TD-DFTB (bottom) theory
levels for L-gln dimer.

Figure 4: Structural features characterizing ground and excited states as well
as CoIn crossing for L-gln system. Molecular structures showing the parameters
to describe the carbonyl stretching (panel a), the proton transfer (panel d)
and the deplanarization (panel g). Probability Distribution Function (PDF) for
the CO distance (panels b and c), PTC (panels e and f) and Deplanarization
modes (panels h and i) for L-gln dimers in vacuum obtained from TSH
trajectories using TD-DFT (middle panels) and TD-DFTB (lower panels). All PDFs
were normalized independently.

The comparison between TD-DFTB and TD-DFT reveals excellent agreement both in
terms of the overall trends observed but also in the specific details of how
the modes change during non-radiative decay. This is a clear evidence of the
high degree of distortions experienced by the system upon excitation, in which
the minimum on both potential energy surfaces $S_{1}$ and $S_{0}$ are shifted
between each other, a feature that is captured by the Tight Binding (TB)
approach. In particular, we observe that an additional stretching for CO
respect to $S_{1}$ distribution is necessary for the system to reach the CoIn
(green lines in panel b and c). Similar behavior is also observed for the PTC
(panels e and f). The primary difference arises in relation to the
deplanarization mode. In TD-DFT (panel h), we observe a higher degree of
deplanarization at the CoIn compared to the distribution of this mode in the
$S_{0}$ state. In contrast, although the distributions obtained with TD-DFTB
reflect the same trends (panel i), the deplanarization mode appears to be much
more constrained compared to our previous findings. Despite the differences,
it should be underscored that, based in our previous work16, the
deplanarization mode is less critical for the photophysics of both L-pyro-amm
and L-gln.

We next turn to comparing the behavior of the non-radiative decay modes that
are important for L-pyro-amm. Since this system is fluorescent, extracting the
decay mechanisms that emerge from the excited state simulations is
particularly challenging due to the fact that there are fewer statistics of
reactive events. In our previous study16 for example, out of a total of 100
trajectories, 4 displayed non-radiative decay from the excited to the ground-
state. Employing excited states simulations for L-pyro-amm as well as model
peptides inspired from the amyloid proteins, we demonstrated that a key
ingredient for enhancing excited state lifetime is facilitated when the
position of the minima for the $S_{0}$ and $S_{1}$ states overlap, which
effectively makes the conical intersection much less accessible16.

Figure 5: Structural features characterizing ground and excited states as well
as CoIn crossing for L-pyro-amm system. Molecular structures showing the
parameters to describe the carbonyl stretching and the proton transfer are in
panel a and d, respectively. Probability Distribution Function (PDF) for the
CO distance (panels b and c), PTC (panels e and f) for L-pyro-amm dimers in
vacuum obtained from TSH trajectories using TD-DFT (middle panels) and TD-DFTB
(lower panels). All PDFs were normalized independently.

In Figure 5 we show the most relevant features, namely the carbonyl stretching
and the proton transfer modes obtained from the TSH simulations at TD-DFT and
TD-DFTB levels. Specifically, we are able to reproduce the significant overlap
in the distribution of the two modes between the $S_{1}$ and $S_{0}$
electronic states (red and blue lines). These findings indicate that the
energy minimum in the potential energy surfaces for each of these modes do not
exhibit significant displacement, suggesting that these modes are restricted
in L-pyro-amm. This feature is also fully consistent with typical aromatic
fluorescent systems68, where the restriction arises from the high electron
delocalization in the planar arrangement of these compounds. In the case of
non-aromatic systems, like the ones we are tackling, the strong interactions
due to SHBs induce this restriction on different modes.

It is important to highlight when comparing the top and bottom panels of
Figure 5, that due to the reduced computational complexity, we are able to run
significantly more TSH trajectories compared to our previous TDDFT simulations
(600 vs 100). DFTB thus permits the observation of many more non-radiative
decays yielding statistics which allow the construction of histograms of
various modes near the CoIn, for instance the green distribution in the case
of TD-DFTB represents 122 decays (panels c and f) and the green vertical lines
for TD-DFT represents only 4 decays (panels b and e). The examination of the
distribution of these modes at the CoIn reveals that the mechanism observed in
L-pyro-amm is similar to that observed in the L-gln system. Specifically, we
find with the DFTB approach that the CO stretching mode undergoes two distinct
signatures at the CoIn. One exhibits a more subtle stretching in comparison to
the distribution of $S_{1}$, while the other shows a more significant
stretching of approximately $1.6\AA$ (green line in panel c in Figure 5). Our
previous TDDFT simulations sample a few events between these two peaks.

The effect of having more statistics with DFTB is observed more clearly in the
PTC distributions (panels e and f) which display more structure compared to
those obtained with TD-DFT. For example, we can capture the double-well
character of the PTC (blue line in panel f) with TD-DFTB, while with TD-DFT we
observe one main peak along with a shoulder (blue line panel e). A distinction
emerges when comparing the behavior of proton transfer between the excited and
ground states in TD-DFTB. Notably, in the excited state, there is a
significant reduction in proton transfer events, evidenced by the presence two
peaks located at $\pm 0.5$ Åwith difference intensity in the PDF (depicted by
the red line in panel f). A similar trend is observed in TDDFT, albeit with
less intensity. In this case, a narrow distribution is noted in $S_{1}$
compared to the ground state, with the hydrogen positioned equidistant from
the two oxygen atoms involved in the SHB (illustrated by the red line in panel
e). Additionally, we observe that the SHB experiences a deformation at the
CoIn (green line in panel f in Figure 5), transforming into a transient normal
HB. This is evident from the peaks at higher values of approximately $\pm
2.5\AA$, which is also present in TD-DFT (green bars in panel e). Although we
cannot rule out the possibility that some of these differences arise from the
underlying potential energy surface along various vibrational modes, the
limited statistics with TD-DFT are more likely to be the source of the
differences we observe. However, we can see that the mechanism predicted by
TD-DFTB is in excellent agreement with our previous works12, 16.

###  4.4 Optical Properties in Realistic Biological Environments

The preceding results provide evidence that DFTB is a very viable approach for
studying the non-radiative decay mechanisms in organic molecular crystals
implicated in non-aromatic fluorescence. In the following section, we turn to
using TD-DFTB to study the optical properties of L-gln and L-pyro-amm where
the environment is included in a more realistic manner within a QM/MM
framework.

####  4.4.1 L-Pyro-Ammonium

We constructed two model QM regions embedded within an MM crystallographic
environment including a dimer and the unit cell as shown in the Figure 6.

Figure 6: Dimer and Unit cell systems in the QM/MM simulations. The QM region
is depicted in ball and stick while the MM in line representations. The real
systems were replicated in the three directions; for the sake of clarity, here
we are showing a 2D representation.

The experimental data used to compare the optical properties in this system
was the excitation and emission spectra, as reported in ref 12. Therefore, to
allow for a one-to-one comparison, we determined the excitation spectra
instead of the absorption. To compute the absorption spectra, one has to
include all the electronic states in the UV-Vis range (250-400 nm), while to
determine the excitation spectra in general one includes only the first few
transitions. We chose to compute the excitation spectra using only the first
transition ($S_{0}\to S_{1}$). This decision is grounded in the well-
established mirror image rule69, which states that if the emission spectrum of
a fluorophore mirrors its absorption or excitation spectra, implying that the
probability of the $S_{1}\to S_{0}$ transition matches that of the $S_{0}\to
S_{1}$ transition. Indeed we observe that this feature is observed in the
experimental spectra shown in the right panel of Figure 7.

Figure 7: Excitation (solid lines) and Fluorescence (dashed lines) spectra for
L-pyro-amm. Left panel shows the average spectra obtained from QM/MM
simulations using the dimer (red line) and the unit cell (blue line) in the QM
region. Right panel shows the experimental spectra (extracted from reference
12, red line) and the average calculated from the QM/MM simulations for the
unit cell in the QM region, where the calculated spectra were shifted to match
the maximum in the experimental fluorescence spectra.

Left panel in figure 7 shows the results for both QM systems in L-pyro-amm in
the QM/MM simulation. We can observe that both systems show similar excitation
maximum (approx. 215 nm) and emission (approx. 290 nm). When we compared these
values with the experimental data (red lines in the right panel) we observe
notable differences. This finding is a well know fact that TD-DFT and
therefore in its approximated variant TD-DFTB, tend to exhibit a notable shift
in the absolute magnitude of the transition energies when comparing with
experimental measurements70, 71. An additional factor contributing to the
discrepancy between our computational outcomes and experimental findings lies
in the inherent limitations of the QM/MM approach in representing the
environmental effects. As demonstrated in Section 4.2, our methodology
achieves outstanding agreement when employing ab-initio ground state
simulations with periodic boundary conditions. Regrettably, this level of
concordance is not replicated in our QM/MM simulations. However, as we will
elucidate in the following, these disparities persist not only in the ground
state but also in excited state simulations. Consequently, focusing on energy
differences becomes a more pertinent and meaningful avenue for comparison. We
have thus determined the Stokes Shift, defined as the difference between the
maximum of the emission band and the absorption or excitation spectra the
results of which are summarized in Table 1. Notably, we find remarkable
agreement between the experimental results and the two model systems for the
Stokes shift. It is also interesting to observe that the dimer system appears
to fully capture the energetics associated with the fluorescence Stokes Shift
which reinforces the use of the dimer model systems in disentangling the non-
radiative decay mechanisms.

System | Absorption/Emission Max. [nm] | FWHMa [nm] | Stokes Shift [nm]
---|---|---|---
Dimer | 211/298 | 120 | 87
Unit Cell | 218/295 | 121 | 77
Experimentb | 354/441 | 92 | 87

Table 1: Stokes Shifts for the different models in the QM/MM simulations

  * •

$a$ FWHM denotes the Full Width at Half Maximun for the fluorescence spectra.

  * •

$b$ Experimental results were extracted from the reference 12.

The right panel of Figure 7 presents the excitation and fluorescence spectra
for both experimental and calculated using the unit cell of L-pyro-amm in the
QM/MM simulations. The calculated spectra were adjusted to align with the peak
of the emission spectra. Impressively, our calculations are in rather modest
agreement with the experimental data. The deviation of 10 nm in the Stokes
Shift is attributed to the fact that the mirror image rule isn’t entirely
applicable in this context. In particular, slight shoulders are noticeable in
the experimental excitation spectra (solid red line), which are absent in the
emission (dotted red line). This is an indication that the excitation spectra
involves more that only one transition from the ground state. Consequently, we
computed the excitation spectra by incorporating an increasing number of
transitions to accommodate this behavior. The detailed results are illustrated
in Figure S3 and table S1 of the Supplementary Information.

####  4.4.2 L-Glutamine

In the case of L-gln, the non-fluorescent system, our objective was to
investigate the non-radiative relaxation mechanism from the $S_{1}$ electronic
state. It’s important to note that the purpose of this section is to gain
insights into the distinctions compared to the dimer in vacuum, rather than
providing a quantitative analysis. Consequently, we conducted simulations for
10 trajectories. Figure 8 depicts the evolution of the potential energy
surfaces ($S_{1}$ and $S_{0}$) along with the pertinent modes coupled to the
CoIn, as discussed earlier, for two TSH trajectories in which we observed the
non-radiative transition employing QM/MM simulations where the QM region is
the unit cell.

Figure 8: Temporal evolution of the most relevant features for two TSH
trajectories in which we observed non-radiative decay in QM/MM simulation for
the unit cell of L-gln system. Top panels show the evolution of the potential
energies surfaces $S_{1}$ (red line) and $S_{0}$ (blue line), where the black
dot indicates the current surface along the time. The other panels show the
evolution of CO mode, Proton transfer Coordinate (PTC) and the Deplanarization
mode as a function of time. Vertical green lines indicate the time of the
transition $S_{1}\to S_{0}$ (CoIn). Black lines show the average of each modes
after the CoIn.

From both trajectories, it is evident that the qualitative findings align with
those discovered in the dimer case of L-gln in vacuum. Notable features
include the stretching of the CO bond, the proton transfer and the
deplanarization of the amide group in the proximity of the Conical
Intersections (CoIn). It’s worth noting that the characteristic lifetimes for
this non-radiative decay are now longer in the QM/MM simulations, at 306 and
421 fs, as compared to the vacuum dimer (approximately 50 fs). This difference
is likely due to the increased complexity in the unit cell system, offering
more degrees of freedom and thus a broader range of relaxation pathways.
However, as discussed previously, the photophysical mechanism governing the
$S_{1}\to S_{0}$ transition remains consistent across both the unit cell QM/MM
and vacuum dimer scenarios.

An interesting distinction arises after the system encounters the CoIn. As
illustrated in Figure 8, all relevant modes persist in an excited vibrational
state on the $S_{0}$ surface. In our simulations, we chose to employ the NVE
ensemble, avoiding the use of a thermostat. This decision was motivated by the
fact that introducing a thermostat could potentially alter the velocities of
the atoms, resulting in a different pattern of vibrational relaxation in the
excited states. Such alterations might be attributed to the system’s
interaction with the thermal bath rather than its inherent photophysical
behavior. Therefore, achieving a comprehensive understanding of vibrational
relaxation would necessitate longer simulation times than the scope of this
methodology allows.

Up to this point, our primary emphasis regarding the fluorescence property has
been centered on exploring the non-radiative decay within L-gln and L-pyro-
amm, utilizing a wide variety of model systems and methodologies. The
rationale behind this approach was to gain insights into the emission process,
given the competitive nature of both radiative and non-radiative pathways.
However, as previously discussed, not all the trajectories in L-gln resulted
in a non-radiative decay to the ground state. Thus, an interesting point of
comparison lies in the emission probabilities between the L-gln and L-pyro-amm
systems. This involves calculating the distribution of oscillator strengths
derived from TSH simulations, where non-radiative decay was not observed.
Consequently, the system maintains its evolution on the $S_{1}$ potential
energy surface. The resulting distribution is illustrated in Figure 9.

Figure 9: Oscillator strength distribution for the transition $S_{1}\to S_{0}$
for L-gln and L-pyro-amm systems, employing the unit cell in the QM/MM
simulations. The inset shows the same distribution at higher oscillator
strength.

As we can discern from Figure 9, both compounds display a peak close to zero
oscillator strength. This aligns with the fact that the L-gln system lacks
fluorescence emission. In the case of L-pyro-amm, its fluorescence is
relatively weak in intensity, several orders below the typical aromatic
compounds1, yet it remains detectable through conventional spectroscopic
methods. Furthermore, the plot reveals that as the oscillator strength
increases, the probability of emission is notably higher for L-pyro-amm
compared to L-gln, which strongly concurs with experimental findings12.
Moreover, a minor population of L-pyro-amm conformations demonstrates
significantly elevated oscillator strength values (refer to the inset in
Figure 9). Although this population is minor, it is evidently greater than
what’s observed for the non-fluorescent L-gln system.

##  5 Conclusions

In this study, we focused on exploring the phenomenon of Non-Aromatic
Fluorescence (NAF) in L-pyro-amm and its comparison with the non fluorescent
L-gln system, leveraging the computational efficiency offered by Density
Functional Tight-Binding Theory (DFTB). Our application of DFTB theory yielded
a robust structural and chemical description of the hydrogen bonds (HB)
present in these structures. The correct behavior of the HBs in the L-gln
system, as well as the characteristic double-well structure in Short Hydrogen
Bonds (SHB) present in L-pyro-amm, was successfully reproduced. These
properties were recognized for their pivotal role in NAF across various
systems.

The integration of Trajectory Surface Hopping (TSH) with TD-DFTB in our study
provided a mutual consensus of results with higher-level theories,
exemplifying the versatility and efficacy of this combined approach. The
computational efficiency inherent in DFTB enabled us to attain additional
statistical significance, particularly advantageous in comprehending the less-
frequent non-radiative decay in the L-pyro-amm fluorescence system. Despite
these strengths, subtle differences surfaced. Agreement in modes like carbonyl
and proton transfer with our previous work highlighted the consistency of key
features. On the other hand, modes such as the de-planarization display more
differences compared to the previous TD-DFT findings. The precise origin of
this observed difference could be an inherent limitation of the DFTB method or
arise as a consequence of our non adiabatic dynamics approach, specifically
the utilization of the LZ method. In this method, the focus in determining non
radiative probabilities is on the energy gap between the states. As soon as
the system reaches a minimum energy gap, the algorithm initiates a transition
between potential energy surfaces. In contrast, FSSH method incorporates non
adiabatic coupling vectors, which significantly influence the hopping
criteria, taking into account both the energy gap and the vectors.
Consequently, our results exhibit a slightly accelerated decay to the ground
state with LZ compared to those in FSSH, reducing the likelihood of observing
the minor contribution from the deplanarization mode. Clearly, this
discrepancy necessitates further investigation in future studies.

In addition, our investigation extended to larger systems within the QM
regions and incorporating environmental effects through the QM/MM approach,
reveals the presence of a consistent non-radiative mechanism in the L-gln
system. The scalability of our simulations to encompass more degrees of
freedom in a larger system elucidates that the observed decay mechanisms
persist but manifest over extended time. The computational efficiency afforded
by DFTB method proves indispensable in capturing these prolonged processes,
which would be impractical with traditional ab-initio methods. In the case of
L-pyro-amm, the introduction of environmental factors induces a discernible
blue-shift in the absorption/excitation spectra. Despite this, a detailed
analysis between ground and excited states simulations allows for the faithful
reproduction of the experimentally observed stoke-shift in this system.

GDM, DB and AH thank to the European Commission for funding on the ERC Grant
HyBOP 101043272. GDM, UNM and AH also acknowledges to CINECA supercomputing
for the resource allocation (project NAFAA-HP10B4ZBB2 and V-CoIns-HP10BY0AET).

## References

  * Morzan et al. 2022 Morzan, U. N.; Díaz Mirón, G.; Grisanti, L.; González Lebrero, M. C.; Kaminski Schierle, G. S.; Hassanali, A. Non-Aromatic Fluorescence in Biological Matter: The Exception or the Rule? _The Journal of Physical Chemistry B_ 2022, _126_ , 7203–7211
  * Ma et al. 2021 Ma, Y.; Zhang, H.; Wang, K.; Cao, D.; Wang, K.; Guan, R.; Zhou, C. The bright fluorescence of non-aromatic molecules in aqueous solution originates from pH-induced CTE behavior. _Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy_ 2021, _254_ , 119604
  * Shukla et al. 2004 Shukla, A.; Mukherjee, S.; Sharma, S.; Agrawal, V.; Kishan, K. R.; Guptasarma, P. A novel UV laser-induced visible blue radiation from protein crystals and aggregates: scattering artifacts or fluorescence transitions of peptide electrons delocalized through hydrogen bonding? _Archives of biochemistry and biophysics_ 2004, _428_ , 144–153
  * Del Mercato et al. 2007 Del Mercato, L. L.; Pompa, P. P.; Maruccio, G.; Torre, A. D.; Sabella, S.; Tamburro, A. M.; Cingolani, R.; Rinaldi, R. Charge transport and intrinsic fluorescence in amyloid-like fibrils. _Proceedings of the National Academy of Sciences_ 2007, _104_ , 18019–18024
  * Tang et al. 2021 Tang, S.; Yang, T.; Zhao, Z.; Zhu, T.; Zhang, Q.; Hou, W.; Yuan, W. Z. Nonconventional luminophores: characteristics, advancements and perspectives. _Chemical Society Reviews_ 2021, _50_ , 12616–12655
  * Zhou et al. 2016 Zhou, Q.; Cao, B.; Zhu, C.; Xu, S.; Gong, Y.; Yuan, W. Z.; Zhang, Y. Clustering-triggered emission of nonconjugated polyacrylonitrile. _Small_ 2016, _12_ , 6586–6592
  * Pinotsi et al. 2013 Pinotsi, D.; Buell, A. K.; Dobson, C. M.; Kaminski Schierle, G. S.; Kaminski, C. F. A label-free, quantitative assay of amyloid fibril growth based on intrinsic fluorescence. _ChemBioChem_ 2013, _14_ , 846–850
  * Pansieri et al. 2019 Pansieri, J.; Josserand, V.; Lee, S.-J.; Rongier, A.; Imbert, D.; Sallanon, M. M.; Kövari, E.; Dane, T. G.; Vendrely, C.; Chaix-Pluchery, O.; others Ultraviolet–visible–near-infrared optical properties of amyloid fibrils shed light on amyloidogenesis. _Nature photonics_ 2019, _13_ , 473–479
  * Grisanti et al. 2020 Grisanti, L.; Sapunar, M.; Hassanali, A.; Došlić, N. Toward understanding optical properties of amyloids: a reaction path and nonadiabatic dynamics study. _Journal of the American Chemical Society_ 2020, _142_ , 18042–18049
  * Arnon et al. 2021 Arnon, Z. A.; Kreiser, T.; Yakimov, B.; Brown, N.; Aizen, R.; Shaham-Niv, S.; Makam, P.; Qaisrani, M. N.; Poli, E.; Ruggiero, A.; others On-off transition and ultrafast decay of amino acid luminescence driven by modulation of supramolecular packing. _Iscience_ 2021, _24_
  * Homchaudhuri and Swaminathan 2001 Homchaudhuri, L.; Swaminathan, R. Novel absorption and fluorescence characteristics of l-lysine. _Chemistry Letters_ 2001, _30_ , 844–845
  * Stephens et al. 2021 Stephens, A. D.; Qaisrani, M. N.; Ruggiero, M. T.; Díaz Mirón, G.; Morzan, U. N.; González Lebrero, M. C.; Jones, S. T.; Poli, E.; Bond, A. D.; Woodhams, P. J.; others Short hydrogen bonds enhance nonaromatic protein-related fluorescence. _Proceedings of the National Academy of Sciences_ 2021, _118_ , e2020389118
  * Tsai et al. 2011 Tsai, Y.-J.; Hu, C.-C.; Chu, C.-C.; Imae, T. Intrinsically fluorescent PAMAM dendrimer as gene carrier and nanoprobe for nucleic acids delivery: bioimaging and transfection study. _Biomacromolecules_ 2011, _12_ , 4283–4290
  * Guan et al. 2020 Guan, R.; Dong, B.; Xu, C.; Zhang, H.; Cao, D.; Lin, W. A strategy to construct fluorescent non-aromatic small-molecules: hydrogen bonds contributing to the unexpected fluorescence. _Chemical Communications_ 2020, _56_ , 4424–4427
  * Fang et al. 2018 Fang, M.; Yang, J.; Xiang, X.; Xie, Y.; Dong, Y.; Peng, Q.; Li, Q.; Li, Z. Unexpected room-temperature phosphorescence from a non-aromatic, low molecular weight, pure organic molecule through the intermolecular hydrogen bond. _Materials Chemistry Frontiers_ 2018, _2_ , 2124–2129
  * Mirón et al. 2023 Mirón, G. D. et al. The carbonyl-lock mechanism underlying non-aromatic fluorescence in biological matter. _Nature Communications_ 2023, _14_
  * Pinotsi et al. 2016 Pinotsi, D.; Grisanti, L.; Mahou, P.; Gebauer, R.; Kaminski, C. F.; Hassanali, A.; Kaminski Schierle, G. S. Proton Transfer and Structure-Specific Fluorescence in Hydrogen Bond-Rich Protein Structures. _Journal of the American Chemical Society_ 2016, _138_ , 3046–3057, PMID: 26824778
  * Grisanti et al. 2017 Grisanti, L.; Pinotsi, D.; Gebauer, R.; Schierle, G. S. K.; Hassanali, A. A. A computational study on how structure influences the optical properties in model crystal structures of amyloid fibrils. _Physical Chemistry Chemical Physics_ 2017, _19_ , 4030–4040
  * Tully 1990 Tully, J. C. Molecular dynamics with electronic transitions. _The Journal of Chemical Physics_ 1990, _93_ , 1061–1071
  * Barbatti 2011 Barbatti, M. Nonadiabatic dynamics with trajectory surface hopping method. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2011, _1_ , 620–633
  * Casida 2009 Casida, M. E. Time-dependent density-functional theory for molecules and molecular solids. _Journal of Molecular Structure: THEOCHEM_ 2009, _914_ , 3–18
  * Cui and Elstner 2014 Cui, Q.; Elstner, M. Density functional tight binding: Values of semi-empirical methods in an ab initio era. _Physical Chemistry Chemical Physics_ 2014, _16_ , 14368–14377
  * Elstner et al. 1998 Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties. _Physical Review B_ 1998, _58_ , 7260–7268
  * Koskinen and Mäkinen 2009 Koskinen, P.; Mäkinen, V. Density-functional tight-binding for beginners. _Computational Materials Science_ 2009, _47_ , 237–253
  * Elstner and Seifert 2014 Elstner, M.; Seifert, G. Density functional tight binding. _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 2014, _372_ , 20120483
  * Gaus et al. 2014 Gaus, M.; Cui, Q.; Elstner, M. Density functional tight binding: application to organic and biological molecules. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2014, _4_ , 49–61
  * Gaus et al. 2013 Gaus, M.; Goez, A.; Elstner, M. Parametrization and benchmark of DFTB3 for organic molecules. _Journal of Chemical Theory and Computation_ 2013, _9_ , 338–354
  * Gruden et al. 2017 Gruden, M.; Andjelkovic, L.; Jissy, A. K.; Stepanović, S.; Zlatar, M.; Cui, Q.; Elstner, M. Benchmarking density functional tight binding models for barrier heights and reaction energetics of organic molecules. _Journal of Computational Chemistry_ 2017, _38_ , 2171–2185
  * Poidevin et al. 2023 Poidevin, C.; Duplaix-Rata, G.; Costuas, K.; Fihey, A. Evaluation of tight-binding DFT performance for the description of organic photochromes properties. _The Journal of chemical physics_ 2023, _158_ , 074303
  * Fihey et al. 2015 Fihey, A.; Hettich, C.; Touzeau, J.; Maurel, F.; Perrier, A.; Köhler, C.; Aradi, B.; Frauenheim, T. SCC-DFTB parameters for simulating hybrid gold-thiolates compounds. _Journal of Computational Chemistry_ 2015, _36_ , 2075–2087
  * Maghrebi et al. 2023 Maghrebi, K.; Chantrenne, I.; Messaoudi, S.; Frauenheim, T.; Fihey, A.; Lien-Medrano, C. R. Rapid Access to Ground- and Excited-State Properties of Gold Nanoclusters Coated with Organic Ligands: Evaluation of the DFTB Method Performance. _Journal of Physical Chemistry C_ 2023, _144_
  * Dolgonos et al. 2010 Dolgonos, G.; Aradi, B.; Moreira, N. H.; Frauenheim, T. An improved self-consistent-charge density-functional tight-binding (SCC-DFTB) set of parameters for simulation of bulk and molecular systems involving titanium. _Journal of Chemical Theory and Computation_ 2010, _6_ , 266–278
  * Niehaus 2009 Niehaus, T. A. Approximate time-dependent density functional theory. _Journal of Molecular Structure: THEOCHEM_ 2009, _914_ , 38–49
  * Niehaus et al. 2001 Niehaus, T. A.; Suhai, S.; Della Sala, F.; Lugli, P.; Elstner, M.; Seifert, G.; Frauenheim, T. Tight-binding approach to time-dependent density-functional response theory. _Physical Review B_ 2001, _63_ , 085108
  * Bonafé et al. 2020 Bonafé, F. P.; Aradi, B.; Hourahine, B.; Medrano, C. R.; Hernández, F. J.; Frauenheim, T.; Sánchez, C. G. A real-time time-dependent density functional tight-binding implementation for semiclassical excited state electron–nuclear dynamics and pump–probe spectroscopy simulations. _Journal of Chemical Theory and Computation_ 2020, _16_ , 4454–4469
  * Uratani and Nakai 2021 Uratani, H.; Nakai, H. Scalable Ehrenfest Molecular Dynamics Exploiting the Locality of Density-Functional Tight-Binding Hamiltonian. _Journal of chemical theory and computation_ 2021, _17_ , 7384–7396
  * Stojanovic et al. 2017 Stojanovic, L.; Aziz, S. G.; Hilal, R. H.; Plasser, F.; Niehaus, T. A.; Barbatti, M. Nonadiabatic dynamics of cycloparaphenylenes with TD-DFTB surface hopping. _Journal of Chemical Theory and Computation_ 2017, _13_ , 5846–5860
  * Mitric et al. 2009 Mitric, R.; Werner, U.; Wohlgemuth, M.; Seifert, G.; Bonacic-Koutecky, V. Nonadiabatic dynamics within time-dependent density functional tight binding method. _The Journal of Physical Chemistry A_ 2009, _113_ , 12700–12705
  * Posenitskiy et al. 2019 Posenitskiy, E.; Rapacioli, M.; Lepetit, B.; Lemoine, D.; Spiegelman, F. Non-adiabatic molecular dynamics investigation of the size dependence of the electronic relaxation in polyacenes. _Physical Chemistry Chemical Physics_ 2019, _21_ , 12139–12149
  * Hourahine et al. 2020 Hourahine, B.; Aradi, B.; Blum, V.; Bonafé, F.; Buccheri, A.; Camacho, C.; Cevallos, C.; Deshaye, M.; Dumitrică, T.; Dominguez, A.; others DFTB+, a software package for efficient approximate density functional theory based atomistic simulations. _The Journal of chemical physics_ 2020, _152_
  * Gaus et al. 2011 Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method (SCC-DFTB). _Journal of Chemical Theory and Computation_ 2011, _7_ , 931–948
  * Brandenburg and Grimme 2014 Brandenburg, J. G.; Grimme, S. Accurate modeling of organic molecular crystals by dispersion-corrected density functional tight binding (DFTB). _The journal of physical chemistry letters_ 2014, _5_ , 1785–1789
  * Miriyala and Řezáč 2017 Miriyala, V. M.; Řezáč, J. Description of non-covalent interactions in SCC-DFTB methods. _Journal of Computational Chemistry_ 2017, _38_ , 688–697
  * Rezac 2017 Rezac, J. Empirical self-consistent correction for the description of hydrogen bonds in DFTB3. _Journal of Chemical Theory and Computation_ 2017, _13_ , 4804–4817
  * Rezac et al. 2009 Rezac, J.; Fanfrlik, J.; Salahub, D.; Hobza, P. Semiempirical quantum chemical PM6 method augmented by dispersion and H-bonding correction terms reliably describes various types of noncovalent complexes. _Journal of Chemical Theory and Computation_ 2009, _5_ , 1749–1760
  * Niehaus and Della Sala 2012 Niehaus, T. A.; Della Sala, F. Range separated functionals in the density functional based tight-binding method: Formalism. _physica status solidi (b)_ 2012, _249_ , 237–244
  * Lutsker et al. 2015 Lutsker, V.; Aradi, B.; Niehaus, T. A. Implementation and benchmark of a long-range corrected functional in the density functional based tight-binding method. _The Journal of chemical physics_ 2015, _143_
  * Subotnik et al. 2016 Subotnik, J. E.; Jain, A.; Landry, B.; Petit, A.; Ouyang, W.; Bellonzi, N. Understanding the surface hopping view of electronic transitions and decoherence. _Annual review of physical chemistry_ 2016, _67_ , 387–417
  * Subotnik and Shenvi 2011 Subotnik, J. E.; Shenvi, N. Decoherence and surface hopping: When can averaging over initial conditions help capture the effects of wave packet separation? _The Journal of chemical physics_ 2011, _134_
  * Crespo-Otero and Barbatti 2018 Crespo-Otero, R.; Barbatti, M. Recent advances and perspectives on nonadiabatic mixed quantum–classical dynamics. _Chemical reviews_ 2018, _118_ , 7026–7068
  * Nikitin 1968 Nikitin, E. Theory of non-adiabatic transitions. recent developement of the Landau-Zener (Linear) model. _Chemische elementarprozesse_ 1968, 43–77
  * Hayashi et al. 2009 Hayashi, S.; Tajkhorshid, E.; Schulten, K. Photochemical reaction dynamics of the primary event of vision studied by means of a hybrid molecular simulation. _Biophysical journal_ 2009, _96_ , 403–416
  * Zhu and Nakamura 1993 Zhu, C.; Nakamura, H. The two-state linear curve crossing problems revisited. III. Analytical approximations for Stokes constant and scattering matrix: Nonadiabatic tunneling case. _The Journal of chemical physics_ 1993, _98_ , 6208–6222
  * Zener 1932 Zener, C. Non-adiabatic crossing of energy levels. _Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character_ 1932, _137_ , 696–702
  * Xie et al. 2019 Xie, W.; Sapunar, M.; Došlić, N.; Sala, M.; Domcke, W. Assessing the performance of trajectory surface hopping methods: Ultrafast internal conversion in pyrazine. _The Journal of chemical physics_ 2019, _150_
  * Kobayashi et al. 2017 Kobayashi, C.; Jung, J.; Matsunaga, Y.; Mori, T.; Ando, T.; Tamura, K.; Kamiya, M.; Sugita, Y. GENESIS 1.1: A hybrid-parallel molecular dynamics simulator with enhanced sampling algorithms on multiple computational platforms. 2017
  * Jung et al. 2015 Jung, J.; Mori, T.; Kobayashi, C.; Matsunaga, Y.; Yoda, T.; Feig, M.; Sugita, Y. GENESIS: a hybrid-parallel and multi-scale molecular dynamics simulator with enhanced sampling algorithms for biomolecular and cellular simulations. _Wiley Interdisciplinary Reviews: Computational Molecular Science_ 2015, _5_ , 310–323
  * Cochran and Penfold 1952 Cochran, W.; Penfold, B. R. The crystal structure of L-glutamine. _Acta Crystallographica_ 1952, _5_ , 644–653
  * Gaus et al. 2013 Gaus, M.; Goez, A.; Elstner, M. Parametrization and benchmark of DFTB3 for organic molecules. _Journal of Chemical Theory and Computation_ 2013, _9_ , 338–354
  * Wang et al. 2004 Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. _Journal of computational chemistry_ 2004, _25_ , 1157–1174
  * Xie et al. 2020 Xie, W.; Holub, D.; Kubar, T.; Elstner, M. Performance of Mixed Quantum-Classical Approaches on Modeling the Crossover from Hopping to Bandlike Charge Transport in Organic Semiconductors. _Journal of Chemical Theory and Computation_ 2020, _16_ , 2071–2084
  * Xie and Domcke 2017 Xie, W.; Domcke, W. Accuracy of trajectory surface-hopping methods: Test for a two-dimensional model of the photodissociation of phenol. _The Journal of Chemical Physics_ 2017, _147_
  * Suchan et al. 2020 Suchan, J.; Janoš, J.; Slavicek, P. Pragmatic approach to photodynamics: Mixed Landau–Zener surface hopping with intersystem crossing. _Journal of Chemical Theory and Computation_ 2020, _16_ , 5809–5820
  * Crespo et al. 2003 Crespo, A.; Scherlis, D. A.; Martí, M. A.; Ordejón, P.; Roitberg, A. E.; Estrin, D. A. A DFT-based QM-MM approach designed for the treatment of large molecular systems: Application to chorismate mutase. _The Journal of Physical Chemistry B_ 2003, _107_ , 13728–13736
  * Senn and Thiel 2009 Senn, H. M.; Thiel, W. QM/MM methods for biomolecular systems. _Angewandte Chemie International Edition_ 2009, _48_ , 1198–1229
  * Morzan et al. 2018 Morzan, U. N.; Alonso de Armino, D. J.; Foglia, N. O.; Ramirez, F.; Gonzalez Lebrero, M. C.; Scherlis, D. A.; Estrin, D. A. Spectroscopy in complex environments from QM–MM simulations. _Chemical reviews_ 2018, _118_ , 4071–4113
  * Boulanger and Harvey 2018 Boulanger, E.; Harvey, J. N. QM/MM methods for free energies and photochemistry. _Current Opinion in Structural Biology_ 2018, _49_ , 72–76
  * Strickler and Berg 1962 Strickler, S.; Berg, R. A. Relationship between absorption intensity and fluorescence lifetime of molecules. _The Journal of chemical physics_ 1962, _37_ , 814–822
  * Albani 2011 Albani, J. R. _Structure and dynamics of macromolecules: absorption and fluorescence studies_ ; Elsevier, 2011
  * Caricato et al. 2010 Caricato, M.; Trucks, G. W.; Frisch, M. J.; Wiberg, K. B. Electronic transition energies: A study of the performance of a large range of single reference density functional and wave function methods on valence and Rydberg states compared to experiment. _Journal of Chemical Theory and Computation_ 2010, _6_ , 370–383
  * Leang et al. 2012 Leang, S. S.; Zahariev, F.; Gordon, M. S. Benchmarking the performance of time-dependent density functional methods. _The Journal of chemical physics_ 2012, _136_