# Magnetic and singlet phases in the three-dimensional periodic Anderson Model

Wiliam S. Oliveira Instituto de Física, Universidade Federal do Rio de Janeiro
Cx.P. 68.528, 21941-972 Rio de Janeiro RJ, Brazil Thereza Paiva Instituto de
Física, Universidade Federal do Rio de Janeiro Cx.P. 68.528, 21941-972 Rio de
Janeiro RJ, Brazil Richard T. Scalettar Department of Physics and Astronomy,
University of California, Davis, CA 95616 USA Natanael C. Costa Instituto de
Física, Universidade Federal do Rio de Janeiro Cx.P. 68.528, 21941-972 Rio de
Janeiro RJ, Brazil

###### Abstract

Heavy fermion materials are compounds in which localized $f$-orbitals
hybridize with delocalized $d$ ones, leading to quasiparticles with large
renormalized masses. The presence of strongly correlated $f$-electrons at the
Fermi level may also lead to long-range order, such as magnetism, or
unconventional superconductivity. From a theoretical point of view, the
“standard model” for heavy fermion compounds is the Periodic Anderson Model
(PAM). Despite being extensively scrutinized, its thermodynamic properties in
three-dimensional lattices have not been carefully addressed by unbiased
methodologies. Here we investigate the 3D PAM employing state-of-the-art
finite temperature auxiliary field quantum Monte Carlo simulations. We present
the behavior of the kinetic energy, the entropy, the specific heat, and the
double occupancy as functions of the temperature and the hybridization
strength. From these quantities, and by the analysis of the spin-spin
correlation functions, we investigate the occurrence of magnetic phase
transitions at finite temperatures, and determine the phase diagram of the
model, including the behavior of the Néel temperature as a function of the
external parameters.

###### pacs:

71.10.Fd, 02.70.Uu

##  I Introduction

Real-space Quantum Monte Carlo (QMC) studies of itinerant electron and
electron-phonon models, such as the Hubbard, Holstein, Su-Schrieffer-Heeger
(SSH), and Periodic Anderson Hamiltonians, have often focused on one and two
dimensions Raczkowski  _et al._ (2020); Paleari  _et al._ (2021); Xing  _et
al._ (2021); Hu  _et al._ (2019). In the case of the Hubbard model, this
choice is partly driven by direct implications in quasi-2D compounds, e.g.
layered geometries relevant to cuprate materials and $d$-wave
superconductivity Scalapino (2012). Similarly, other 2D materials/lattices,
such as graphene and transition metal dichalcogenides, may exhibit interesting
topological properties or long-range ordered phases Kotov  _et al._ (2012);
Manzeli  _et al._ (2017). However, the choice of lower dimensional lattices is
also a practical one: with the system sizes limited to ${\cal
O}(10^{2}$-$10^{3})$ sites, reliable finite-size scaling analysis becomes
challenging in 3D geometries. Unfortunately, due to the Mermin-Wagner theorem,
lower dimension may also lead to restrictions for the occurrence of long-range
ordered phases. This is the case of Hubbard-like models, which at half-filling
any possible long-range antiferromagnetic phase should occur only at zero
temperature. 111The reduced (Ising) order parameter symmetry of the Holstein
case does support a finite $T$ transition in 2D. Studies within dynamical
mean-field theory (DMFT), dynamical cluster approximation, or diagrammatic
Monte Carlo avoid this issue, while introducing other approximations or
limitations.

A significant exception is the great numerical effort devoted to investigating
the 3D (cubic) Hubbard Hirsch (1987); Scalettar  _et al._ (1989); Staudt  _et
al._ (2000); Fuchs  _et al._ (2011) model, the 3D (cubic) Holstein Cohen-Stead
_et al._ (2020), and (cubic perovskite) SSH models Cohen-Stead  _et al._
within QMC simulations. These studies have established quantitatively the
dependence of their critical temperatures at half-filling – $T_{N}$ (Néel) to
the Hubbard model, and $T_{\rm cdw}$ (CDW) to the Holstein model. Such a
difficult task is alleviated by the increasing computational resources, with
parallelization protocols, as well as the development of new methodologies
(e.g., machine learning or Langevin accelerations). Although this leap to 3D
geometries is a significant accomplishment, so far it has been carried out
mostly for single-band models.

As a paradigm for two-band systems, the periodic Anderson model (PAM)
considers the interplay of localized and itinerant fermions. Unlike the half-
filled 3D Hubbard model, where one has ground state antiferromagnetic even for
$U_{f}/t\to 0$, in the PAM there is a finite critical ratio of interband to
conduction electron hybridization $V/t$, separating antiferromagnetic and
singlet phases at $T=0$. Thus, one has a substantially more rich phase
diagram. However its quantitative critical boundaries have not been determined
by unbiased methodologies at the present moment. Early attempts to describe
the properties of the 3D PAM by QMC was performed in Refs.McMahan  _et al._ ,
1998; Huscroft  _et al._ , 1999; Paiva  _et al._ , 2003, in which the behavior
of the spin-spin correlation functions was examined for different values of
$V/t$ and $U_{f}/t$, but for small lattices and without a finite-size scaling
analysis. The main goal of this paper is to bridge this gap by extending work
on real space QMC in 3D to the PAM, and comparing it with the well-known 2D
case.

As the standard model for heavy-fermion compounds, the study of the 3D PAM is
a significant step towards understanding the thermodynamic properties of these
materials Bernhard and Lacroix (1999). In particular, the application of
pressure changes the conduction electron bandwidth $W=12\,t$ and hybridization
$V$ between localized and conduction bands (with a lesser effect on on-site
repulsion). These changes affect the nature of both high and low-temperature
responses, from the transition temperatures to transport properties.
Identifying these energy scales is crucial since, given a density-functional
theory parameterization of $W$ and $V$, our phase transitions can be used for
quantitative benchmarking of heavy fermion properties – such as the Cerium
volume collapse transition Held  _et al._ (2001). The present paper is
organized as follows: In Sec. II we present the Periodic Anderson model, and
highlight the main features of the QMC method together with the quantities of
interest. Our results are presented in Sec. III, while the main conclusions
are summarized in Sec. IV.

##  II Model and Method

The PAM Hamiltonian reads

$\displaystyle\hat{\mathcal{H}}=$
$\displaystyle-t\sum_{\langle\mathbf{i},\mathbf{j}\rangle,\sigma}\left(d_{\mathbf{i}\sigma}^{\dagger}d_{\mathbf{j}\sigma}^{\phantom{\dagger}}+\mathrm{h.c.}\right)-V\sum_{\mathbf{i},\sigma}\left(d_{\mathbf{i}\sigma}^{\dagger}f_{\mathbf{i}\sigma}^{\phantom{\dagger}}+\mathrm{H.c.}\right)$
$\displaystyle-\mu\sum_{\mathbf{i},\sigma,\alpha}n^{\alpha}_{\mathbf{i}\sigma}+U_{f}\sum_{\mathbf{i}}\left(n^{f}_{\mathbf{i}\uparrow}-\frac{1}{2}\right)\left(n^{f}_{\mathbf{i}\downarrow}-\frac{1}{2}\right)\,\,,$
(1)

where the sums run over a three-dimensional simple cubic lattice, with
$\langle{\bf i,j}\rangle$ denoting nearest-neighbor sites, and $\alpha=d$ or
$f$. The first term on the r.h.s. of Eq.(II) corresponds to the hopping of
conduction d-electrons, while the last one describes the Coulomb repulsion on
localized $f$-orbitals. The hybridization between these two orbitals is given
by an on-site inter-orbital hopping $V$ and the filling of the lattice is set
by adjusting the chemical potential $\mu$. Hereafter, we set the energy scale
as being units of the hopping integral, $t$, and define $U_{f}/t=6$. We have
assumed there is no difference in the single particle $d$ and $f$ levels, that
is, $\epsilon_{f}-\epsilon_{d}=0$.

We investigate the thermodynamic properties of the 3D PAM by performing
finite-temperature auxiliary-field QMC simulations; namely, determinant
QMCBlankenbecler  _et al._ (1981); Hirsch (1985); White  _et al._ (1989). The
DQMC method is an unbiased technique commonly used to investigate tight-
binding Hamiltonians: it maps a $d$-dimensional interacting system to a
noninteracting ($d$+1)-dimensional one, with the additional imaginary-time
coordinate $0\leq\tau\leq\beta$, where $\beta$ is the inverse temperature.
Within this approach, one separates the one-body ($\hat{\mathcal{K}}$) and
two-body ($\hat{\mathcal{P}}$) pieces in the partition function by using the
Trotter-Suzuki decomposition, i.e. by defining $\beta=L_{\tau}\Delta\tau$,
with $L_{\tau}$ being the number of imaginary-time slices, and $\Delta\tau$
the discretization grid. Then

$\displaystyle{\cal Z}$ $\displaystyle={\rm
Tr}\,e^{-\beta\hat{\mathcal{H}}}={\rm
Tr}\,\big{[}\big{(}e^{-\Delta\tau(\hat{\mathcal{K}}+\hat{\mathcal{P}})}\big{)}^{L_{\tau}}\big{]}$
$\displaystyle\approx{\rm
Tr}\,\big{[}e^{-\Delta\tau\hat{\mathcal{K}}}e^{-\Delta\tau\hat{\mathcal{P}}}e^{-\Delta\tau\hat{\mathcal{K}}}e^{-\Delta\tau\hat{\mathcal{P}}}\cdots\big{]},$
(2)

with an error proportional to $(\Delta\tau)^{2}$, but being exact in the limit
$\Delta\tau\to 0$. The resulting partition function is rewritten in quadratic
(single-body) form through a discrete Hubbard-Stratonovich transformation
(HST) on the two-body terms, $e^{-\Delta\tau\hat{\mathcal{P}}}$. This HST
introduces discrete auxiliary fields with components on each of the space and
imaginary-time lattice coordinates, which are sampled by Monte Carlo
techniques. In this work we choose $\Delta\tau=0.1$, so that the error from
the Trotter-Suzuki decomposition is less than, or comparable to, statistical
errors from the Monte Carlo sampling. More details about the method are
discussed in Refs. dos Santos, 2003; Gubernatis  _et al._ , 2016; Becca and
Sorella, 2017, and references therein.

Although DQMC is unbiased, its low-temperature application is restricted to
systems with particle-hole or other symmetries, owing to the minus-sign
problem Loh  _et al._ (1990); Troyer and Wiese (2005); Mondaini  _et al._
(2022). For this reason, our focus is on half-filling, $\mu=0$, where the sign
problem is absent. Fortunately, this density is of considerable interest, both
because of the strong magnetic order favored by commensurate filling, and by
the materials for which half filling is appropriate (e.g., the undoped parent
compounds of the cuprate superconductors).

In order to examine signatures for the formation of local momentum and
antiferromagnetic long-range order, we first investigate the thermodynamic
properties, namely, the internal energy and specific heat per site,

$e(T,V)=\frac{1}{N}\langle\hat{\mathcal{H}}\rangle\leavevmode\nobreak\ ,$ (3)
$c(T)=\frac{1}{N}\frac{d\langle\hat{\mathcal{H}}\rangle}{dT}\leavevmode\nobreak\
,\,\,$ (4)

respectively. We also investigate the magnetic properties by performing
calculations of the antiferromagnetic spin structure factor, defined as the
Fourier transform of the spin-spin correlation functions for $f$ orbitals,

$S_{ff}(\boldsymbol{q})=\frac{1}{3N}\sum_{\mathbf{i},\mathbf{j}}e^{i\mathbf{q}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})}\,\langle\vec{S}_{f\mathbf{i}}\cdot\vec{S}_{f\mathbf{j}}\rangle,$
(5)

with $\boldsymbol{q}=(\pi,\pi,\pi)$. For singlet formation, we examine a
correlator function

$C_{\textit{fd}}=\frac{1}{3N}\sum_{\mathbf{i}}\langle\vec{S}_{f\mathbf{i}}\cdot\vec{S}_{d\mathbf{i}}\rangle.$
(6)

Here, we define the fermionic spin operators as
$\vec{S}_{f\mathbf{i}}=(f^{\dagger}_{\uparrow\mathbf{i}}\,\,f^{\dagger}_{\downarrow\mathbf{i}})\,\vec{\sigma}\,\begin{pmatrix}f_{\uparrow\mathbf{i}}\\\
f_{\downarrow\mathbf{i}}\end{pmatrix}$, with $\vec{\sigma}$ being the Pauli
spin matrices, while the 1/3 factor in Eqs. (5) and (6) is for the average
value along one spin component.

We extract critical behavior by means of the antiferromagnetic correlation
ratio

$R_{c}(L)=1-\frac{S_{ff}(\textbf{q}+\delta\textbf{q})}{S_{ff}({\textbf{q}})}$
(7)

with $|\delta\textbf{q}|=\frac{2\pi}{L}$. This is a renormalization-group
invariant observable, which crossing points for different lattice sizes
provides the critical points Kaul (2015); Darmawan  _et al._ (2018).

In addition to these equal-time correlation functions, we also calculate
appropriate unequal-time quantities, including the magnetic susceptibility

$\chi(\mathbf{q})=\chi_{dd}(\mathbf{q})+2\chi_{df}(\mathbf{q})+\chi_{ff}(\mathbf{q}),$
(8)

with

$\chi_{\alpha\gamma}(\mathbb{q})=\frac{1}{3N}\sum_{\mathbf{i},\mathbf{j}}\int_{0}^{\beta}d\tau\bigg{\langle}{\vec{S}}_{\alpha,\mathbf{j}}(\tau)\cdot{\vec{S}}_{\gamma,\mathbf{i}}(0)\bigg{\rangle}e^{i\mathbf{q}\cdot(\mathbf{r}_{i}-\mathbf{r}_{j})}\,,$
(9)

($\alpha,\gamma=d$, or $f$) where we examine the uniform and staggered cases,
$\boldsymbol{q}=(0,0,0)$ and $(\pi,\pi,\pi)$, respectively.

Figure 1: (Color online) (a) Energy as a function of temperature for different
$V/t$, and the (b) specific heat as a function of the temperature for
$V/t=1.4$; closed symbols represent data from numerical differentiating the
energy values, whereas dashed line represents the differentiation of an
exponential fit to the energy Paiva  _et al._ (2001). (c) Specific heat from
the fit procedure for the same values of $V/t$ as in (a). All data are for
$6\times 6\times 6$ lattices.

Finally, we also evaluate Nuclear Magnetic Resonance (NMR) quantities, such as
the time relaxation rate for the $f$-electrons. For the low-frequency limit of
the dynamic susceptibility, it is defined as (see, e.g. Ref. Curro, 2009)

$T_{1,ff}^{-1}=\gamma^{2}k_{B}T\lim_{\omega\to
0}\sum_{\mathbf{q}}A^{2}(\mathbf{q})\,\frac{\chi_{ff}^{{}^{\prime\prime}}(\mathbf{q},\gamma)}{\hbar\omega},$
(10)

where $A^{2}(\mathbf{q})$ is the square of the Fourier transform of the
hyperfine interaction, and $\gamma$ is the gyromagnetic ratio. The latter is
related to the nuclear magnetic moment by $\gamma\hbar=g\mu_{N}\sqrt{I(I+1)}$,
with $\mu_{N}$ being the nuclear magneton, $g$ the nuclear $g$-factor, and $I$
the nuclear spin. $T_{1,ff}^{-1}$ quantifies a characteristic time in which a
component of the nuclear spin (of a given $f$-orbital) reaches equilibrium
after an external perturbation (magnetic field pulse). It is a dynamical (real
frequency) quantity whose numerical evaluation usually requires an analytic
continuation of the imaginary-time DQMC data. Instead, we follow Ref. Randeria
_et al._ , 1992, which performs an approximation to this procedure, leading to

$\displaystyle\frac{1}{T_{1}T}=\frac{1}{\pi^{2}T^{2}}\frac{1}{N}\sum_{\mathbf{i}}\bigg{\langle}S_{f,\mathbf{i}}(\tau=\beta/2)S_{f,\mathbf{i}}(0)\bigg{\rangle}\leavevmode\nobreak\
.$ (11)

##  III Results

###  III.1 Thermodynamic properties

We start our analysis by discussing the thermodynamic properties of the system
– unless otherwise indicated, the following results are obtained for $6\times
6\times 6$ lattice (i.e., 432 sites). First, we investigate the internal
energy, Eq. (3), displayed as a function of $T/t$ in Fig. 1 (a). The
derivative of the internal energy with respect to the temperature yields the
specific heat, Eq. (4), whose behavior should indicate the occurrence of phase
transitions. Figure 1 (b) displays $c(T)$ for $V/t=1.4$, exhibiting a two-peak
structure. Since numerical differentiation for a few data points is usually
noisy, one may perform a nonlinear fit of the QMC energy points, and
differentiate the fitted curve, as shown by the dashed line in Fig. 1 (b).
Here we define such a fit on an exponential basisPaiva  _et al._ (2001), by
the function $e_{\rm fit}(T)=a_{0}+\sum_{n=1}^{M}a_{n}{\rm exp}\,(-\beta
n\Delta)$, with a cut-off in $M$ (fixing $M\leq 10$ is enough). The agreement
indicates the rubustness of the fitting approach, which we repeat for other
values of hybridization $V/t$, whose results are shown in Fig. 1 (c).

Figure 2: (Color online) (a) Singlet correlator and (b) antiferromagnetic
structure factor as functions of the hybridization for different temperatures
for $6\times 6\times 6$ lattices. The inset in (b) shows the antiferromagnetic
structure factor as a function of inverse temperature for hybridizations
$V/t=0.9$ and $V/t=1.7$.

There are several important features concerning the behavior of $c(T)$ to be
emphasized. First, a high-temperature hybridization-independent peak at
$T/t\sim 1.45$ is evident, which is due to the formation of local moments
Paiva  _et al._ (2001). We also find a low-temperature peak, that is pushed to
higher temperatures and decreases in intensity as $V/t$ increases. This low-
temperature peak is also seen for the two and three-dimensional Hubbard model
Duffy and Moreo (1997); Paiva  _et al._ (2001); Staudt  _et al._ (2000), and
usually is associated to collective spin-wave excitations. In particular, for
the three-dimensional Hubbard model this low-temperature peak position signals
the Néel temperature since the ground state is AFM. For the 3D PAM, on the
other hand, previous calculations for the specific heat were unable to resolve
the low-$T$ peak for values of $V/t$ below the AF-singlet quantum critical
point McMahan  _et al._ (1998); Huscroft  _et al._ (1999), as we have done in
this paper. As one expects a phase transition from the antiferromagnetic to
the singlet phase, the position of the low-$T$ peak may not uniquely identify
the Néel temperature for the 3D PAM. Therefore, to probe whether the low-$T$
peaks indicate the presence of long-range spin correlations, we next
investigate the magnetic properties of the system.

###  III.2 Magnetic properties

We proceed to examine the singlet correlator, $C_{fd}$: a measure of the on-
site spin-spin correlation between an f-moment and the spin of a d-electron.
We find this quantity is always negative, pointing to antiferromagnetic
coupling and singlet formation that increases in strength as $V/t$ increases.
As displayed in Fig. 2(a), $C_{fd}$ has a rapid increase in modulus around
$V/t\approx 1.3$, providing an initial indication of a change into a singlet
phase.

In order to go beyond this rough indication from the singlet correlator, we
also analyze the $f$-orbital AFM structure factor, Eq. (5), shown in Fig.
2(b). $S_{ff}$ grows rapidly as $V/t$ decreases from $V/t\sim 2$, especially
for larger $\beta$, indicating the formation of long range AF order on the $f$
orbitals as the singlet correlations diminish (Fig. 2(a)). However, after
reaching a maximum at $V/t\sim 1.3$, $S_{ff}$ begins falling. Naively, the
$f$-orbital magnetic behavior is symmetric about the maximum. However,
examining the temperature dependence of $S_{ff}$ reveals a profoundly
different story. Specifically, for $V/t=0.9$, $S_{ff}$ while seemingly small,
grows very rapidly as $\beta$ increases – see the inset of Fig. 2(b). This is
in sharp contrast to the behavior at $V/t=1/7$, where $S_{ff}$ saturates at a
small value as the temperature is lowered (increasing $\beta$), as also shown
in the inset of Fig. 2(b). This, together with the fact that the singlet
correlator is small for $V/t\sim 0.9$ suggests that the low $V/t$ region still
has AF order, but that the transition temperature $T_{N}$ is below the
smallest temperature simulated $T_{N}<t/12$.

Figure 3: (Color online) Total uniform susceptibility (a) and staggered
susceptibility in the f sites (b) for different values of $V/t$ on $6\times
6\times 6$ lattices.

To further understand the magnetic behavior we turn to the uniform and
staggered magnetic susceptibilities, shown in Figs. 3(a) and (b),
respectively. For high temperatures, both uniform and staggered
susceptibilities show the expected $1/T$ Curie behavior of free spins, while
for low temperatures and small hybridizations (e.g., $V/t=1.2$ and 1.3), the
uniform susceptibility $\chi(0)$ displays a clear peak, followed by a strong
suppression as $T$ decreases. We recall that, in antiferromagnetic systems,
such a peak of $\chi(0)$ occurs at the Néel temperature. Here, as our system
may exhibit both AFM and singlet phases, the peak turns into a cusp pushed to
higher temperatures for larger $V/t$, blurring the identification of $T_{N}$.

As will be further discussed later, such a cusp in $\chi(0)$ occurs due to the
emergence of the spin-singlet phase, and gives the energy scale for the spin
gap formation, $\Delta_{s}$. Such a gap in the spin excitations is clearly
noticed when analyzing the staggered susceptibility, $\chi_{ff}(\pi,\pi,\pi)$,
displayed in Fig. 3(b). In the thermodynamic limit, the staggered
susceptibility should diverge at $T_{N}$, while $\chi_{ff}(\pi,\pi,\pi)\to 0$
in presence of a spin gap. Indeed, as shown in Fig. 3(b) the susceptibility
displays a fast increase with decreasing $T/t$ for $V/t\lesssim 1.4$, while is
strongly suppressed for larger hybridizations.

Figure 4: (Color online) Spin-lattice relaxation rate of the 3D PAM as a
function of temperature, for different hybridizations and $6\times 6\times 6$
lattices.

Early DQMC studies of two-dimensional PAM have shown evidence of a QCP at
$V_{c}/t\approx 1.0$, separating an antiferromagnetically (AFM) ordered ground
state from a spin-singlet Vekić  _et al._ (1995). The problem was recently
revisited by DQMC and Dynamical Vertex Approximation analyses, which show a
similar value for the QCP Hu  _et al._ (2017); Schäfer  _et al._ (2019). One
would expect that some finite temperature quantities should reflect the
characteristics of the different low $T$ phases, even for the 2D case (where
the occurrence of long-range ordered AFM phase is forbidden at finite
temperatures). Indeed, this is the case of NMR measurements, in particular the
spin-lattice relaxation rate, Eq. (10), which may signal the existence of a
spin gap Mendes-Santos  _et al._ (2017). Previous DQMC studies of the 2D PAM
have shown that, at the critical point $V_{c}$, the $1/T_{1}T$ is almost
constant at low temperatures, while increasing (reducing) for $V<V_{c}$
($V>V_{c}$)Costa  _et al._ (2019). Figure 4 displays the behavior of the
relaxation rate as a function of temperature for different values of $V/t$,
for the 3D PAM (linear size $L=6$). Notice that, within the AFM phase
($V/t=1.2$ or $1.3$) $1/T_{1}T$ approaches a finite nonzero value as $T/t\to
0$, consistent with the absence of a spin gap, i.e. the presence of spin-wave
excitations. On the other hand, for larger $V$, $1/T_{1}T$ decreases
monotonically when $T$ is lowered, reflecting a spin-gapped ground state.
Interestingly, the change in behavior of $1/T_{1}T$ occurs around $V/t\sim
1.4$, where both results from specific heat and susceptibilities support the
disappearance of AFM order.

In order to estimate and track the spin-gap $\Delta_{s}$ opening, it is useful
to analyze the unequal-time spin-spin correlation functions,
$C(\mathbf{i},\mathbf{j},\tau)=\langle
S^{z}_{\mathbf{j}}(\tau)S^{z}_{\mathbf{i}}(0)\rangle$. In particular, in the
presence of a spin-gap, the asymptotic behavior of the unequal-time staggered
structure factor becomes

$\displaystyle\lim_{\tau\to\infty}\frac{1}{N}\sum_{j}(-1)^{|\mathbf{i}-\mathbf{j}|}C(\mathbf{i}=0,\mathbf{j},\tau)\propto\exp(-\Delta_{s}\tau)\leavevmode\nobreak\
,$ (12)

from which $\Delta_{s}$ may be extracted. Following this procedure, Fig. 5
exhibits our estimation for $\Delta_{s}$ as a function of the hybridization,
and at different temperatures. For large hybridization values, we notice that
the gap increases with decreasing temperature, remaining finite at low-$T$.
For small hybridizations, on the other hand, the gap is small and decreases
with decreasing temperature; that is, one has a spin gapless phase, in line
with our previous expectation of AFM order. The change in such behavior occurs
for $V/t\approx 1.5$, where the spin-gap is temperature independent; in fact,
as shown in the inset of Fig. 5, the extrapolation to $T/t\to 0$ for
$V/t=1.45$ leads to a vanishing gap. In summary, the results presented in
Figs. 3-5 point to a QCP at $V_{c}/t=1.45\pm 0.05$, separating an AFM phase at
low hybridization from a singlet phase at large $V/t$.

Figure 5: (Color online) The spin gap $\Delta_{s}$ as a function of
hybridization for different inverse temperatures. As is apparent, there are
two distinct behaviors as we move from high to low hybridization, crossing a
temperature-independent critical region around $V/t\sim$ 1.5. Inset: The spin-
gap as a function of temperature for $V/t=1.45$.

Although we have estimated the location of the QCP, the Néel temperatures (for
$V<V_{c}$) are ill-defined by the specific heat and homogeneous susceptibility
behavior. A thorough probe of $T_{N}$ is provided by a finite-size scaling
analysis of the correlation ratio, Eq. (7), a renormalization-group invariant,
and shown in Fig. 6 as a function of temperature, for different hybridizations
and system sizes. The crossing of the $Rc(L)$ for different lattice sizes
provides the critical temperature. Since the correlation ratio shows no abrupt
features, the paramagnetic metallic to antiferromagnetic insulator finite
temperature phase transition is more likely to be continuous. For the range of
$V/t$ analyzed we have $T_{N}/t\simeq 0.08-0.10$. However, for smaller
hybridizations $T_{N}$ will decrease, as previously discussed in the context
of Fig. 2(b). Going to such low $T/t$ is challenging to achieve in three-
dimensional geometries. It is worth noting that this computational limitation
also affects the determination of $T_{N}$ in the 3D single band Hubbard model,
where $T_{N}\sim t\,e^{-a\sqrt{t/U}}$ is very small at weak couplingSlater
(1951); Hirsch (1987); Scalettar  _et al._ (1989); Nagaosa (1999); Fuchs  _et
al._ (2011); Kozik  _et al._ (2013).

Figure 6: (Color online) Determinant quantum Monte Carlo results for the AF
correlation ratio, as function of inverse temperature, for different lattice
sizes and hybridizations, at fixed $U_{f}/t=6$. The crossing points separate
the paramagnetic phase and Nèel phases.

We conclude by comparing the energy scales of our results with those from the
literature. Early QMC attempts to investigate the three-dimensional PAM were
performed in Refs. McMahan  _et al._ , 1998; Huscroft  _et al._ , 1999; Paiva
_et al._ , 2003, providing hints that such an AFM long-range order requires
very low temperatures to be achieved. By ‘low-temperature’ we mean by
comparison with the 3D Hubbard model Staudt  _et al._ (2000); Hirschmeier  _et
al._ (2015); Khatami (2016), whose maximum of the Néel temperature is
$T_{N}/t\approx 0.33$, and occurs at $U_{f}/t\approx 8$. Here, we provide
$T_{N}$ for the 3D PAM, showing that these energy scales are indeed much lower
than those from the Hubbard model222Interestingly, the Néel temperatures
obtained by DMFT Held  _et al._ (2000); Schäfer  _et al._ (2019) are
compatible with our QMC ones.. This difference is expected due to the weak
indirect coupling between $f$-electrons mediated by conduction ones – first
addressed by DMFT studies Held  _et al._ (2000). In view of this, it is
expected that the value of $V_{c}/t$ should change for different values of
$U_{f}$, but the maxima of $T_{N}$ should keep similar values as those
presented for $U_{f}/t=6$.

###  III.3 Transport properties

In the ground state, the AFM and singlet phases are both insulators. However,
at finite temperatures, a crossover from a metallic to insulating behavior
should occur at an appropriate energy scale. Therefore, as a complementary
study, we also probe such a crossover as the temperature is lowered, showing
that its features can be identified from the behavior of thermodynamic
quantitiesKim  _et al._ (2020).

Figure 7: (Color online) Kinetic energy $K$ as a function of temperature for
different values of $V/t$ for $6\times 6\times 6$ lattices. The inset shows
$K$ as a function of $T/t$ for $V/t=1.4$ and two different values of imaginary
time discretization: $\Delta\tau=0.05$ (squares) and $\Delta\tau=0.10$
(triangles).

We first examine the $c$-electrons kinetic energy per particle,

$\big{\langle}K\big{\rangle}=-\frac{t}{N}\bigg{\langle}\sum_{\langle\mathbf{i},\mathbf{j}\rangle,\sigma}(d_{\mathbf{i},\sigma}^{\dagger}d_{\mathbf{j},\sigma}+H.c)\bigg{\rangle},$
(13)

whose behavior is presented in Fig. 7 as a function of temperature, for
different values of hybridization. For all hybridization values, the following
behavior is noticed: at high temperatures, $\big{\langle}K\big{\rangle}$ is
reduced as temperature is lowered, reaching a minimum for a given energy scale
$T_{\rm min(\langle K\rangle)}$ and, therefore, increasing as $T\to 0$ when
$T<T_{\rm min(\langle K\rangle)}$. The sign change in
$\frac{\partial\big{\langle}K\big{\rangle}}{\partial T}$, from positive to
negative, is consistent with an evolution from metallic to insulator behavior,
respectively. In order to emphasize the absence of imaginary time
discretization errors, the inset in Fig. 7 shows the kinetic energy per
particle for $V/t=1.4$ as a function of temperature for two different values
of $\Delta\tau$, from which is clear that our choice of $\Delta\tau$ does not
affect the results of $\big{\langle}K\big{\rangle}$. This result of the
kinetic energy (suggesting an insulator for $T<T_{\rm min(\langle
K\rangle)}$), combined with the fact that $\frac{\partial}{\partial
T}\chi(0)\leq 0$ at lower temperatures (consistent with a Pauli metal for
$d$-electrons + local moment Curie behavior for $f$-electrons, instead of a
spin-gapped insulator; see Fig. 3), point out to a bad metal phase.

To further emphasize this bad metallic phase, we also investigate the dc
conductivity,

$\sigma_{dc}=\frac{\beta^{2}}{\pi}\Lambda_{xx}(\mathbf{q=0},\tau=\beta/2),$
(14)

in which

$\Lambda_{xx}(\mathbf{q},\tau)=\langle
j_{x}(\mathbf{q},\tau)j_{x}(-\mathbf{q},0)\rangle,$ (15)

with $j_{x}(\mathbf{q},\tau)$ being the Fourier transform of the unequal-time
current-current correlation functions

$j_{x}(\mathbf{i},\tau)=\mathrm{e}^{\tau\mathcal{H}}\left[it\sum_{\sigma}\left(d_{\mathbf{i}+\mathbf{x}\sigma}^{\dagger}d_{\mathbf{i}\sigma}^{\phantom{\dagger}}-d_{\mathbf{i}\sigma}^{\dagger}d_{\mathbf{i}+\mathbf{x}\sigma}^{\phantom{\dagger}}\right)\right]\mathrm{e}^{-\tau\mathcal{H}};$
(16)

see, e.g., Refs.Trivedi  _et al._ , 1996; Denteneer  _et al._ , 1999; Mondaini
_et al._ , 2012. The behavior of $\sigma_{dc}$ as a function of temperature is
displayed in Fig. 8, which shows $\sigma_{dc}\to 0$ as $T\to 0$ for all $V/t$,
consistent with an insulating phase ground state for the entire range of
hybridizations. However, this insulating feature – i.e.
$\frac{\partial\sigma_{dc}}{\partial T}<0$ – changes depending on the
temperature, with $\sigma_{dc}$ exhibiting a maximum at $T_{\rm
max(\sigma_{dc})}$. Interestingly, the values of $T_{\rm max(\sigma_{dc})}$
are consistent with $T_{\rm min(\langle K\rangle)}$, both signaling a
crossover into a bad metallic phase.

Figure 8: (Color online) $dc$ conductivity as a function of the temperature,
for different values of $V/t$ and $6\times 6\times 6$ lattices.

##  IV Conclusions

In this work, we have investigated the thermodynamic, magnetic, and transport
properties of the three-dimensional Periodic Anderson Model through unbiased
DQMC simulations. For the thermodynamic properties, we examined the specific
heat for different hybridizations. We find a broad, high-temperature peak that
is independent of $V/t$, while the low-$T$ peak, on the other hand, is
strongly $V$-dependent. This low-$T$ peak is associated with the emergence of
the AFM phase for $V<V_{c}$. However, it does not give precise information
about the singlet phase for $V>V_{c}$ 333We expect that, for larger lattice
sizes, the low-$T$ peak would be enhanced for $V<V_{c}$, and suppressed for
$V>V_{c}$, emphasizing the transition..

Figure 9: (Color online) Energy scales for the three dimensional Periodic
Anderson Model for different values of $V/t$. For comparison, the Néel
temperature $T_{N}$ of the 3D Hubbard model, when plotted versus on-site
interaction $U_{f}/t$, exhibits a broad maximum with $T_{N}/t\sim 0.3$ in the
range $6\lesssim U_{f}/t\lesssim 10$.

To get more quantitative insight into the singlet formation, we investigated
the AF structure factor, and correlator function for different $V/t$ and
temperatures. In particular, the AF structure factor pointed out the energy
scale for the quantum critical point. This quantum phase transition between an
AFM state and a singlet state was obtained by analyzing the spin gap, the NMR
time relaxation rate, and the staggered magnetic susceptibility, which agree
with $V_{c}/t=1.45\pm 0.05$. The Néel temperature within the AFM phase was
probed by the antiferromagnetic correlation ratio and lies in the range
$T/t\simeq 0.08-0.10$ for all $V/t$ investigated. Finally, by investigating
the transport properties, we find out that there is a region in which the
system behaves as a bad metal, i.e. it presents insulator features for the
kinetic energy (namely $\frac{\partial\langle K\rangle}{\partial T}<0$), while
keeping metallic ones for the homogeneous susceptibility (namely
$\frac{\partial\chi}{\partial T}\leq 0$ – a Pauli metal for $d$-electrons +
local moment Curie behavior for $f$ ones).

In summary, our work presents detailed finite temperature results for the
three-dimensional Periodic Anderson Model in a simple cubic lattice, allowing
the quantitative determination of different energy scales. To summarize our
findings, Fig. 9 presents these energy scales, and the resulting finite
temperature phase diagram. In particular, we show [i] the quantum critical
point, QCP, [ii] the temperature of the maxima for the homogeneous
susceptibility, $\chi(0)_{\rm max}$, [iii] the minima for the $d$-electrons
hopping energy, $K_{\rm min}$, [iv] the temperature of the maxima of the
conductivity $\sigma_{\rm max}$, [v] the low-temperature specific heat peak
position $c(T)_{\rm max}$, and finally [vi] the Néel temperatures $T_{N}$
obtained from a Finite Size Scaling analysis. Knowing precise values for these
energy scales at half-filling is an essential step towards examining other
more complex facets of this model which emerge with doping, including the
occurrence of $d$-wave superconductivity due to frustration effects Wu and
Tremblay (2015), or the enhancement of magnetism due to depletion Assaad
(2002); Costa  _et al._ (2018); Jiang (2020). Indeed, Fig. 9 provides detailed
data for the PAM across an interesting part of its phase diagram encompassing
the AF-singlet QCP, which may help guide future theoretical and experimental
work on heavy fermion materials.

## ACKNOWLEDGMENTS

The authors are grateful to the Brazilian Agencies Conselho Nacional de
Desenvolvimento Científico e Tecnológico (CNPq), Coordenação de
Aperfeiçoamento de Pessoal de Ensino Superior (CAPES), Fundação Carlos Chagas
de Apoio à Pesquisa do Estado do Rio de Janeiro (FAPERJ), and Instituto
Nacional de Ciência e Tecnologia de Informação Quântica (INCT-IQ) for funding
this project. N.C.C. acknowledges financial support from CNPq, Grant No.
313065/2021-7, and from FAPERJ – Fundação Carlos Chagas Filho de Amparo à
Pesquisa do Estado do Rio de Janeiro –, Grant No. 200.258/2023
(SEI-260003/000623/2023). T. P. acknowledges financial support from CNPq Grant
Nos. 403130/2021-2 and 308335/2019-8, and FAPERJ Grant Nos. E-26/200.959/2022
and E-26/210.100/2023. The work of R.S. was supported by the grant DE-
SC-0014671, funded by the U.S. Department of Energy, Office of Science.

## References

  * Raczkowski  _et al._ (2020) Marcin Raczkowski, Robert Peters, Thị Thu Phùng, Nayuta Takemori, Fakher F. Assaad, Andreas Honecker, and Javad Vahedi, “Hubbard model on the honeycomb lattice: From static and dynamical mean-field theories to lattice quantum Monte Carlo simulations,” Phys. Rev. B 101, 125103 (2020).
  * Paleari  _et al._ (2021) G. Paleari, F. Hébert, B. Cohen-Stead, K. Barros, RT. Scalettar, and G. G. Batrouni, “Quantum Monte Carlo study of an anharmonic Holstein model,” Phys. Rev. B 103, 195117 (2021).
  * Xing  _et al._ (2021) Bo Xing, Wei-Ting Chiu, Dario Poletti, R. T. Scalettar, and George Batrouni, “Quantum Monte Carlo simulations of the 2D Su-Schrieffer-Heeger model,” Phys. Rev. Lett. 126, 017601 (2021).
  * Hu  _et al._ (2019) Danqing Hu, Jian-Jun Dong, and Yi-feng Yang, “Hybridization fluctuations in the half-filled periodic Anderson model,” Phys. Rev. B 100, 195133 (2019).
  * Scalapino (2012) D. J. Scalapino, “A common thread: The pairing interaction for unconventional superconductors,” Rev. Mod. Phys. 84, 1383–1417 (2012).
  * Kotov  _et al._ (2012) Valeri N. Kotov, Bruno Uchoa, Vitor M. Pereira, F. Guinea, and A. H. Castro Neto, “Electron-electron interactions in graphene: Current status and perspectives,” Rev. Mod. Phys. 84, 1067–1125 (2012).
  * Manzeli  _et al._ (2017) S. Manzeli, D. Ovchinnikov, D. Pasquier, O.V. Yazyev, and A. Kis, “2D transition metal dichalcogenides,” Nature Reviews Materials 2, 17033 (2017).
  * Note (1) The reduced (Ising) order parameter symmetry of the Holstein case does support a finite $T$ transition in 2D.
  * Hirsch (1987) J. E. Hirsch, “Simulations of the three-dimensional Hubbard model: Half-filled band sector,” Phys. Rev. B 35, 1851–1859 (1987).
  * Scalettar  _et al._ (1989) R. T. Scalettar, D. J. Scalapino, R. L. Sugar, and D. Toussaint, “Phase diagram of the half-filled 3D Hubbard model,” Phys. Rev. B 39, 4711–4714 (1989).
  * Staudt  _et al._ (2000) R Staudt, M. Dzierzawa, and A. Muramatsu, “Phase diagram of the three-dimensional Hubbard model at half filling,” Eur. Phys. J. B 17, 411–415 (2000).
  * Fuchs  _et al._ (2011) Sebastian Fuchs, Emanuel Gull, Lode Pollet, Evgeni Burovski, Evgeny Kozik, Thomas Pruschke, and Matthias Troyer, “Thermodynamics of the 3D Hubbard model on approaching the Néel transition,” Phys. Rev. Lett. 106, 030401 (2011).
  * Cohen-Stead  _et al._ (2020) B. Cohen-Stead, Kipton Barros, ZY Meng, Chuang Chen, R. T. Scalettar, and G. G. Batrouni, “Langevin simulations of the half-filled cubic Holstein model,” Phys. Rev. B 102, 161108 (2020).
  * (14) Benjamin Cohen-Stead, Kipton Barros, Richard Scalettar, and Steven Johnston, “A hybrid Monte Carlo study of bond-stretching electron-phonon interactions and charge order in the bismuthate family of superconductors,” cond-mat/2208.02339 .
  * McMahan  _et al._ (1998) A.K. McMahan, C. Huscroft, R.T. Scalettar, and E.L. Pollock, “Volume-collapse transitions in the rare earth metals,” Journal of Computer-Aided Materials Design 5, 131–162 (1998).
  * Huscroft  _et al._ (1999) Carey Huscroft, A. K. McMahan, and R. T. Scalettar, “Magnetic and thermodynamic properties of the three-dimensional periodic Anderson hamiltonian,” Phys. Rev. Lett. 82, 2342–2345 (1999).
  * Paiva  _et al._ (2003) Thereza Paiva, Gökhan Esirgen, Richard T. Scalettar, Carey Huscroft, and A. K. McMahan, “Doping-dependent study of the periodic Anderson model in three dimensions,” Phys. Rev. B 68, 195111 (2003).
  * Bernhard and Lacroix (1999) Ben Hur Bernhard and Claudine Lacroix, “Thermodynamics of the Anderson lattice,” Phys. Rev. B 60, 12149–12154 (1999).
  * Held  _et al._ (2001) K. Held, A. K. McMahan, and R. T. Scalettar, “Cerium volume collapse: Results from the merger of dynamical mean-field theory and local density approximation,” Phys. Rev. Lett. 87, 276404 (2001).
  * Blankenbecler  _et al._ (1981) R. Blankenbecler, D. J. Scalapino, and R. L. Sugar, “Monte Carlo calculations of coupled boson-fermion systems. I,” Phys. Rev. D 24, 2278–2286 (1981).
  * Hirsch (1985) J. E. Hirsch, “Two-dimensional Hubbard model: Numerical simulation study,” Phys. Rev. B 31, 4403–4419 (1985).
  * White  _et al._ (1989) S. R. White, D. J. Scalapino, R. L. Sugar, E. Y. Loh, J. E. Gubernatis, and R. T. Scalettar, “Numerical study of the two-dimensional Hubbard model,” Phys. Rev. B 40, 506–516 (1989).
  * dos Santos (2003) Raimundo R. dos Santos, “Introduction to quantum Monte Carlo simulations for fermionic systems,” Brazilian Journal of Physics 33, 36 – 54 (2003).
  * Gubernatis  _et al._ (2016) J Gubernatis, N Kawashima, and P Werner,  _Quantum Monte Carlo Methods: Algorithms for Lattice Models_ (Cambridge University Press, Cambridge, England, 2016).
  * Becca and Sorella (2017) Federico Becca and Sandro Sorella,  _Quantum Monte Carlo approaches for correlated systems_ (Cambridge University Press, Cambridge, England, 2017).
  * Loh  _et al._ (1990) E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White, D. J. Scalapino, and R. L. Sugar, “Sign problem in the numerical simulation of many-electron systems,” Phys. Rev. B 41, 9301–9307 (1990).
  * Troyer and Wiese (2005) Matthias Troyer and Uwe-Jens Wiese, “Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations,” Phys. Rev. Lett. 94, 170201 (2005).
  * Mondaini  _et al._ (2022) Rubem Mondaini, Sabyasachi Tarat, and Richard T Scalettar, “Quantum critical points and the sign problem,” Science 375, 418–424 (2022).
  * Kaul (2015) Ribhu K. Kaul, “Spin nematics, valence-bond solids, and spin liquids in $\mathrm{SO}(n)$ quantum spin models on the triangular lattice,” Phys. Rev. Lett. 115, 157202 (2015).
  * Darmawan  _et al._ (2018) Andrew S. Darmawan, Yusuke Nomura, Youhei Yamaji, and Masatoshi Imada, “Stripe and superconducting order competing in the Hubbard model on a square lattice studied by a combined variational Monte Carlo and tensor network method,” Phys. Rev. B 98, 205132 (2018).
  * Paiva  _et al._ (2001) Thereza Paiva, R. T. Scalettar, Carey Huscroft, and A. K. McMahan, “Signatures of spin and charge energy scales in the local moment and specific heat of the half-filled two-dimensional Hubbard model,” Phys. Rev. B 63, 125116 (2001).
  * Curro (2009) N.J. Curro, “Nuclear magnetic resonance in the heavy fermion superconductors,” Reports on Progress in Physics 72, 026502 (2009).
  * Randeria  _et al._ (1992) Mohit Randeria, Nandini Trivedi, Adriana Moreo, and Richard T. Scalettar, “Pairing and spin gap in the normal state of short coherence length superconductors,” Phys. Rev. Lett. 69, 2001–2004 (1992).
  * Duffy and Moreo (1997) Daniel Duffy and Adriana Moreo, “Specific heat of the two-dimensional Hubbard model,” Phys. Rev. B 55, 12918–12924 (1997).
  * Vekić  _et al._ (1995) M. Vekić, J. W. Cannon, D. J. Scalapino, R. T. Scalettar, and R. L. Sugar, “Competition between Antiferromagnetic Order and Spin-Liquid Behavior in the Two-Dimensional Periodic Anderson Model at Half Filling,” Phys. Rev. Lett. 74, 2367–2370 (1995).
  * Hu  _et al._ (2017) Wenjian Hu, Richard T. Scalettar, Edwin W. Huang, and Brian Moritz, “Effects of an additional conduction band on the singlet-antiferromagnet competition in the periodic Anderson model,” Phys. Rev. B 95, 235122 (2017).
  * Schäfer  _et al._ (2019) T. Schäfer, A. A. Katanin, M. Kitatani, A. Toschi, and K. Held, “Quantum criticality in the two-dimensional periodic Anderson model,” Phys. Rev. Lett. 122, 227201 (2019).
  * Mendes-Santos  _et al._ (2017) T. Mendes-Santos, N. C. Costa, G. Batrouni, N. Curro, R. R. dos Santos, T. Paiva, and R. T. Scalettar, “Impurities near an antiferromagnetic-singlet quantum critical point,” Phys. Rev. B 95, 054419 (2017).
  * Costa  _et al._ (2019) N. C. Costa, T. Mendes-Santos, T. Paiva, N. J. Curro, R. R. dos Santos, and R. T. Scalettar, “Coherence temperature in the diluted periodic Anderson model,” Phys. Rev. B 99, 195116 (2019).
  * Slater (1951) J.C. Slater, “Magnetic effects and the Hartree-Fock equation,” Physical Review 82, 538 (1951).
  * Nagaosa (1999) Naoto Nagaosa,  _Quantum field theory in strongly correlated electronic systems_ (Springer Science & Business Media, 1999).
  * Kozik  _et al._ (2013) E. Kozik, E. Burovski, V. W. Scarola, and M. Troyer, “Néel temperature and thermodynamics of the half-filled three-dimensional Hubbard model by diagrammatic determinant Monte Carlo,” Phys. Rev. B 87, 205102 (2013).
  * Hirschmeier  _et al._ (2015) Daniel Hirschmeier, Hartmut Hafermann, Emanuel Gull, Alexander I. Lichtenstein, and Andrey E. Antipov, “Mechanisms of finite-temperature magnetism in the three-dimensional Hubbard model,” Phys. Rev. B 92, 144409 (2015).
  * Khatami (2016) Ehsan Khatami, “Three-dimensional Hubbard model in the thermodynamic limit,” Phys. Rev. B 94, 125114 (2016).
  * Note (2) Interestingly, the Néel temperatures obtained by DMFT Held  _et al._ (2000); Schäfer  _et al._ (2019) are compatible with our QMC ones.
  * Held  _et al._ (2000) K. Held, C. Huscroft, R. T. Scalettar, and A. K. McMahan, “Similarities between the Hubbard and periodic Anderson models at finite temperatures,” Phys. Rev. Lett. 85, 373–376 (2000).
  * Kim  _et al._ (2020) Aaram J. Kim, Fedor Simkovic, and Evgeny Kozik, “Spin and charge correlations across the metal-to-insulator crossover in the half-filled 2D Hubbard model,” Phys. Rev. Lett. 124, 117602 (2020).
  * Trivedi  _et al._ (1996) Nandini Trivedi, Richard T. Scalettar, and Mohit Randeria, “Superconductor-insulator transition in a disordered electronic system,” Phys. Rev. B 54, R3756–R3759 (1996).
  * Denteneer  _et al._ (1999) P. J. H. Denteneer, R. T. Scalettar, and N. Trivedi, “Conducting phase in the two-dimensional disordered Hubbard model,” Phys. Rev. Lett. 83, 4610–4613 (1999).
  * Mondaini  _et al._ (2012) Rubem Mondaini, K. Bouadim, Thereza Paiva, and Raimundo R. dos Santos, “Finite-size effects in transport data from quantum Monte Carlo simulations,” Phys. Rev. B 85, 125127 (2012).
  * Note (3) We expect that, for larger lattice sizes, the low-$T$ peak would be enhanced for $V<V_{c}$, and suppressed for $V>V_{c}$, emphasizing the transition.
  * Wu and Tremblay (2015) Wei Wu and A.-M.-S. Tremblay, “$d$-wave superconductivity in the frustrated two-dimensional periodic Anderson model,” Phys. Rev. X 5, 011019 (2015).
  * Assaad (2002) F. F. Assaad, “Depleted Kondo lattices: Quantum Monte Carlo and mean-field calculations,” Phys. Rev. B 65, 115104 (2002).
  * Costa  _et al._ (2018) N. C. Costa, M. V. Araújo, J. P. Lima, T. Paiva, R. R. dos Santos, and R. T. Scalettar, “Compressible ferrimagnetism in the depleted periodic Anderson model,” Phys. Rev. B 97, 085123 (2018).
  * Jiang (2020) Mi Jiang, “Enhanced tendency towards $d$-wave pairing and antiferromagnetism in a doped staggered periodic Anderson model,” Phys. Rev. B 102, 085119 (2020).


# Computational study of optical absorption spectra of helicenes as applied to
strain sensing

Veera Sundararaghavan  Department of Aerospace Engineering, University of
Michigan, Ann Arbor MI 48109, U.S.A<EMAIL_ADDRESS>Vikas Varshney
Materials and Manufacturing Directorate, Air Force Research Laboratory,
Wright-Patterson Air Force Base, OH 45433, U.S.A.  Davide Simone  Materials
and Manufacturing Directorate, Air Force Research Laboratory, Wright-Patterson
Air Force Base, OH 45433, U.S.A.

###### Abstract

Helicenes, a class of organic molecules consisting of ortho-fused benzene
rings in a spring-like configuration have found several interesting
applications in nonlinear optical materials and opto-electronic devices. Under
the action of strain, i.e., via mechanical stretching or compression, the
optical absorption spectra of helicenes change which can be employed for
strain sensing. The present study presents a detailed investigation of the
optical absorption spectra of helicenes using density functional theory along
with calculations of the changes in the spectra during mechanical axial
stretching or compression of helicenes. The electronic band gap followed a
non-symmetric parabolic form with the amount of applied strain. A lowering of
the gap in stretched or compressed helicenes compared to the pristine helicene
was observed. The compressed state shows a smaller energy gap compared to
tension for the same strain magnitude. A detailed inspection of the optical
absorption spectra shows that compressive states show significantly lower
absorption at higher optical energies (shorter wavelengths) which can provide
greater sensitivity to the strain measurement.

###### keywords:

Helicenes, Optical absorption spectra, circular dichroism, strain

## Introduction

Helicenes comprise a class of organic molecules in which planar cyclic rings,
such as benzene, thiophene, and their combination thereof, etc., are fused at
their ortho–positions, creating a helical structure [1]. Even though these
molecules are devoid of ‘classical’ chiral centers, they exhibit chiral
characteristics that arise from their intrinsic molecular architecture (i.e.,
handedness of the helix), leading to their unique optical properties, such as
exceptionally large optically rotatory dispersion (ORD) [2]. Due to their
unique properties, helicenes and their functional derivatives have been
employed in a wide variety of chemical disciplines such as catalysis [3, 4],
polymer chemistry [5, 6], biological sciences [7], non-linear optics [8, 9],
and optoelectronics [10, 11, 12].

In addition to their ability to greatly influence polarized light, helicenes
present themselves as potential sensing entities due to their intrinsic
spring-like characteristics, where a change in either electronic, optical, or
piezoelectric response, could be examined while deforming the helicene either
under compression or elongation. To date, there has been very limited non-
theoretical research towards investigating the different types of responses
due to external mechanical perturbations or stimuli in helicenes.
Investigation of charge transport in helicenes was carried out [13, 14],
indicating that current along the helicene molecule depends on its length (or
number of fused rings) and degree of deformation. Specifically, under
stretching or compression along the helical axis, a non-monotonic (U-shaped
curve) current–strain relationship was observed by Guo et al. [13]. In a
different study, Vacek and co-workers computed several orders of modulation in
on-off ratios in conductance, when helicenes were utilized as a molecular
switch [14].

In the context of modulating optical features, while experimental optical
absorption and circular dichroism spectra of [5]–[9] helicenes have been
reported, the effect of tensile and compressive deformation of the helicenes
on their optical absorption is yet to be studied. In principle, while it is
possible to isolate and strain individual helicene molecules via optical
tweezers [15, 16], it is challenging to systematically investigate the effect
of strain on opto-spectral response for individual helicene molecules via
experiment. Theoretical investigation of such responses and identifying
helicene characteristics that maximize the response sensitivity is key for
their successful incorporation as in-situ, non-invasive sensing components in
various applications.

Computational chemistry provides a powerful and complementary tool for
assessing the physical and chemical properties of organic molecules, often
allowing direct comparison of predictions with experimental observations. In
this work, we focus on modeling the electronic characteristics via density
functional theory (DFT) and computation of excited states in the form of
UV–Vis and electronic circular dichroism (ECD) spectra of pristine and
strained helicenes using time-dependent density functional theory (TDDFT) to
explore their feasibility as strain sensing entities.

Figure 1: Comparison of the HOMO–LUMO gap of helicenes compared to its planar
phenacenes and acene counterparts (shown in the left) using DFT-LDA method.
The x–axis indicates the number of aromatic rings (the ‘n’ in [n]-helicene)

(a) (b)

Figure 2:  Decrease in HOMO–LUMO gap with an increase in helicene length. (a)
The trend follows Kuhn’s equation (shown as inset), plotted here for even
numbered helicenes. (b) Trend zoomed for shorter helicenes. Inset in (a) shows
the zigzag length of even helicenes used in the formula as highlighted in
yellow.

## Results

### Electronic gap Prediction

Fig. 1 shows the predicted gap (from local density approximation (LDA))
between the highest occupied molecular orbital (HOMO) and lowest unoccupied
molecular orbital (LUMO) for helicenes and its planar analogues (phenacenes
and acenes) as a function of the number of benzene units within the modeled
structures. Although LDA under-predicted the true HOMO-LUMO gap (i.e., HF
predicted gap), the observed trends are in good agreement with what has been
observed in experimental data as aggregated in Ref. [17]. It should be noted
that the similar HOMO-LUMO gaps of helicenes and phenacenes are due to their
cross-conjugated nature of the sp2 bonds in the aromatic rings relative to the
conjugated structure of acenes. We find that the gap of helicenes is
marginally lower than that of phenacenes and is primarily attributed to the
non-planar nature of helicenes and the deviation of bond lengths associated
with external aromatic bonds of the helix versus internal aromatic bonds [18].
On the contrary, no cross-conjugation occurs in acenes, leading to sharper
decay in the gap with increasing acene length.

Interestingly, the HOMO-LUMO gap for all the structures was predicted to
decrease with increasing length. Due to the conjugated nature of studied
structures, it is imperative that the $\pi$-electrons are involved in
electronic transitions between the HOMO and LUMO. As the length of the
helicene and its planar analogue increases, the conjugated network on the
molecule becomes increasingly delocalized, and hence, the transition from HOMO
to LUMO is expected to require progressively lower energy, as can be seen from
the computed gap.

The decrease of energy gap as a function of length for conjugated organic
chains is embodied in Kuhn’s equation [19] as follows:

$\Delta E=\frac{h^{2}}{8mL^{2}}(N+1)+V_{o}(1-\frac{1}{N})$ (1)

where $N=4n+2$ is the number of $\pi$ electrons, $h$ is the Planck’s constant,
$m$ is the electron mass, $L$ is the zig–zag length (as described next), and
$V_{0}$ is the HOMO–LUMO gap of an infinitely long helicene, which is a
fitting parameter. The zig–zag length, $L$, of a representative even helicene
is shown in yellow in the inset in Fig. 2(a) and is equal to 3$n$+3 times the
C-C bond length of 1.4 Å and follows the bond conjugation. Kuhn’s equation is
a single parameter fit based on an appropriate choice of $V_{0}$, the
HOMO–LUMO gap of an infinitely long helicene. The value of $V_{0}$ is taken to
be $0.06$ Ha ($1.6$ eV) in the fit shown in Fig. 2(a).

Band gap of planar analogs have been studied in the past in the form of
zero–dimensional (finite length) and one-dimensional (infinite length)
graphene nano ribbons [20]. Among these 1D ribbons, while armchair ribbons are
either calculated as semiconducting or metallic depending on their width,
zigzag ribbons are found to be metallic for all widths [21]. 1D armchair
ribbons follow a rule that the nanoribbon is metallic if the number of rows of
carbon (width-wise) is 3$M$-1, where $M$ is an integer. Phenacenes and acenes
analogs studied in our study can respectively be treated as finite length arm-
chair (4 rows of carbon atoms) and zig-zag nanoribbons. It can be seen that
acene analogs decay toward a zero HOMO-LUMO gap with increasing length,
corroborating with the metallic nature of infinitely long metallic
nanoribbons. Similarly, phenacenes, which have four rows of carbon atoms along
ribbon-width are expected to have a non–zero band gap at infinite lengths. The
reported DFT band gap for phenacenes of infinite length is 2.6 eV[20].
Helicenes are expected to have a lower HOMO-LUMO gap and the estimated value
of $1.6$ eV appears reasonable. Note that Kuhn’s equation does not account for
features such as non–planarity of molecules other than through the chosen
parameter $V_{0}$. Fig. 2(b) shows that odd numbered $n$–helicenes have
similar energy gaps as ($n$+1) even helicenes although one would expect the
gap to be slightly higher owing to the shorter zig–zag length (3$n$+2) times
the C-C bond length for the same $V_{0}$.

(a) (b) (c)

Figure 3: Comparison of predicted optical absorption with experimentally
observed values. The bars indicate the absorption computed at the first 25
excitation frequencies. (a) [4]-helicene (b) [5]-helicene (c) [6]-helicene.
Experimental data from Ref. [22, 23].

### Prediction of Optical Absorption Spectra

Next, we discuss and compare the predicted optical absorption with the
experimental literature in Figs. 3–4. The vertical bars in these plots
indicate the absorption intensity computed at the first 25 excitation
frequencies and the data is directly compared against published experimental
spectra for helicene molecules of various lengths reported in Refs. [22, 23,
24, 25].

Fig(s). 3a & 3b compare our predictions against the experimental data of Clar
and Stewart (1952) [22] for [4]-helicene and [5]-helicene, respectively. The
figures also contain the experimental spectra of the planar phenacene
counterpart (solid line) which, as previously shown in Fig. 1, shows an
optical gap at a marginally lower wavelength (higher energy) compared to the
helicene. The [6]-helicene optical absorption data is compared against another
experiment by Newman et al. [23] in Fig. 3(c). The optical absorption spectra
of helicenes contain two prominent peaks near the HOMO–LUMO gap called the
para ($p$) band and the $\beta$ band. The $p$-bands are related to fused
benzene rings, while $\beta$ peaks are associated with higher energy, high
intensity B type excitations [26, 27, 28]. Both peak wavelengths are well
captured by our simulations. The experimental data also shows that helicene
$\beta$ peak is shifting to lower energy with respect to its phenacene
analogue. This shift has been previously explained to be the result of greatly
increased overlapping of H atoms due to the twisting and the subsequent loss
in resonance energy [22]. With [4]–, [5]– and [6]–helicenes, coplanar overlap
of end benzene rings is not possible.

Optical spectra has been analyzed in terms of electronic transitions in Table
1 below. Only the three states with the lowest energy (longest wavelengths)
are analyzed. The configurations with highest contributions to the excited
state are given in a notation having commonality across molecules. The HOMO-
LUMO transition is always considered 1-1* regardless of system size. HOMO-1 to
LUMO+1 is 2-2* etc. The percentage contribution of these orbitals to the
excited state are indicated. The oscillator and rotatory strengths (in
$10^{-40}$ cgs) are provided. The electronic excitations of $C_{2}$-symmetric
helicenes can be classified according to the two types of polarized
transitions: 1-1* and 2-2* are polarized along the short $C_{2}$-axis and are
of the ‘a’ type, while 1-2* and 2-1* are ‘b’-type excitations are polarized
along the long axis (orthogonal to $C_{2}$). In even helicenes, the lowest
energy excitations tend to be of the ‘b’ type and for odd helicenes, these are
of the ‘a’ type, which is consistent with past calculations[26, 29]. In the
context of helicenes, the three excited states indicated correspond to
$\alpha$, p and $\beta$ bands respectively. The $\alpha$ peak always consists
of the low energy low intensity (‘L’ type) CI states. The $\beta$ bands are of
the high energy high intensity (‘B’ type) and is prominent compared to the
other two bands in helicenes of lower lengths. The prominence of the $\beta$
peak reduces with higher helicenes as described in Ref. [27], which we
postulate arises from end benzene rings axially overlapping due to a complete
pitch of rotation. The overlap gives rise to an enhanced probability of intra-
molecular pi-stacking interactions leading to lower energy electronic
transitions.

(a) (b) (c)

Figure 4: Comparison of predicted optical absorption with experimentally
observed values for (a) [7]-helicene (b) [8]-helicene (c) [9]-helicene.
Experimental data from [24, 25]

Fig. 4(a) compares the optical gap and absorption spectra of [7]–, [8]–, and
[9]– helicenes, which involve a direct and increasing overlap of benzene rings
(due to complete pitch), with experiments reported in Refs [24, 25]. The
figure shows excellent agreement with experimental spectra and highlights the
decreasing trend of lower energy optical gap with increasing helicene length.
The trend is also consistent with the Kuhn’s equation which predicts lower
energy gaps due to extended delocalization as a function of helicene length.

(a) (b)

Figure 5: Comparison of predicted spectra against data from [30, 24, 25] (a)
UV-Vis spectra (b) ECD spectra.

(a) (b) (c)

Figure 6: (a) Axial straining of [9]–helicene. (b) The change in HOMO-LUMO gap
as a function of strain for [9]-helicene, the data is plotted against data
from Guo et al (2015) for [12]-helicene. (c) The HOMO and LUMO orbitals
(plotted together) for 40% compression, 25% and 100% tension.

The optical absorption spectra for [7]–, [8]–, and [9]– helicenes are
aggregated together along with experimental data from Refs [24, 25] in Fig.
5(a). In addition to optical absorption spectra, we also calculated the ECD
spectra for [6]– to [8]–helicene and compare the findings against experimental
data in Ref. [30]. The predicted as well as experimental ECD spectra are shown
in Fig. 5(b). We find that the predicted ECD spectra compare favorably with
experiments. It should be noted that due to the calculation of only 25 excited
states, the higher wavelength (low energy) results can be more accurately
compared against the experiments. Helicenes with even-numbered rings typically
show a higher peak differential absorption of left and right circularly
polarized light as compared to odd–helicenes as seen in both experiments and
DFT simulations. The good agreement between measured and predicted optical
spectra (as discussed in Figs 3 and 4) validated the level of
accuracy/fidelity of first-principle calculations necessary to carry out
helicene deformation (compression and elongation) calculations, which we
discuss next.

[n] | %contribution | CI parentage | Wavelength (nm) | Oscillator strength | Rotatory strength
---|---|---|---|---|---
4 | 52,47 | 2-1*,1-2* | 339.5 | 0 | -0.1
| 88,10 | 2-2*,1-1* | 330.59 | 0.02 | -19
| 48,45 | 1-2*,2-1* | 279.82 | 0.78 | 112
5 | 40,22 | 1-1*,2-2* | 342.13 | 0 | -3
| 53,32 | 2-1*,1-1* | 326.34 | 0.03 | 88
| 49,20 | 1-2*,2-2* | 305.19 | 0.17 | 197
6 | 57,40 | 2-1*,1-2* | 358.35 | 0 | 2
| 22,77 | 2-2*,1-1* | 343.64 | 0 | -0.6
| 55,39 | 1-2*,2-1* | 315.59 | 0.42 | 686
7 | 62,36 | 1-1*,2-2* | 377.41 | 0 | 3
| 95,2.7 | 2-1*,1-2* | 361.79 | 0.03 | 274
| 78,17 | 1-2*,3-1* | 349.44 | 0.05 | 291
8 | 60,28 | 1-2*,2-1* | 392.47 | 0 | 25
| 62,21 | 1-1*,2-2* | 373.21 | 0.01 | 112
| 66,21 | 2-1*,1-2* | 368.05 | 0.08 | 702
Table 1: Calculated spectral data (%contribution corresponding to CI
parentage, corresponding wavelength, oscillator and rotatory strengths (in
$10^{-40}$ cgs) are indicated. Only top three long wavelength transitions are
considered.

### Effect of Applied Strain on HOMO–LUMO gap & Spectra Modulation

In this section, we discuss how the tensile and compressive axial deformation
affects the optical absorption spectra of helicene molecules towards exploring
their potential as embedded sensing modalities using the representative
example of [9]–helicene. In this context, we start our discussion with the
visualization of different strained states of [9]–helicene as well as the
estimation of change in respective HOMO–LUMO gap with respect to an undeformed
helicene as shown in Fig(s). 6(a) & (b). It can be seen from the figure that
the change in gap follows a non-monotonic, asymmetric, inverted parabola-like
form with respect to the amount of applied strain. For both large compressive
and tensile deformation, a lower band gap is observed as compared to the
pristine, un-deformed helicene. The parabolic form is unique to helical
graphene structure and its analogs. For example, under bending as shown in Ref
[31], bandgap of semiconducting nanotubes monotonically decreases with
bending, while that of metallic nanotubes would increase with bending. In
helicenes, the LUMO energies are sensitive to straining and follow a similar
parabolic form to strain as the HOMO–LUMO gap itself. This unique trend in
LUMO energies can be associated with the twisted structure of helicenes. On
compression, the gap decreases due to the higher overlap of orbitals and
subsequent loss of resonance energy. Fig. 6(c)(top) depicts the HOMO and LUMO
orbitals for helicene axially strained under compression (40% strain)
depicting overlap of orbitals from the center towards the benzene units at the
free ends.

The behavior on stretching is more interesting, as one would expect the gap to
increase towards those what is observed in planar phenacenes. This is indeed
the case until a critical strain of 25%. Beyond this strain, although HOMO
energies continue to increase upon tensile straining as expected, the LUMO
energies show an opposing trend, i.e., decrease upon further tensile
straining. In a previous study by Guo et al. [13] on [12]-helicene, a critical
strain of 23.53% was reported where the inflection occurs towards a decreasing
gap under increased tension (curve depicted in Fig. 6(b)) .

In Fig. 6(c), we show the behavior of orbital wave functions at the critical
strain of 25% where the highest gap is realized. At this strain, overlap of
orbitals out-of-plane of the benzene rings as seen at 40% compression is
completely removed, leading to an increase in the HOMO-LUMO gap. Beyond this
critical strain, the HOMO-LUMO gap begins to decrease. The HOMO-LUMO orbital
structure is illustrated for the case of 100% tensile strain in Fig.
6(c)(bottom). In Ref. [13], it was found that stretching beyond 23.53% tensile
strain tends to localize the LUMO wave function toward the inner helix.
However, that study enforced proportional elongation of all helicene atoms
before optimization which could have biased the structure towards a geometry
with a homogeneous strain distribution. In our study, we see that the geometry
optimization after straining leads to inhomogeneous straining of helicenes
(see Fig 6(a)). At 100% strain, the free ends take up significant deformation
via bending of the chain, while the central helix is subject to a relatively
smaller amount of stretching. Alternating inner helix bonds are highly
deformed during tension, with a few bonds stretching towards the $sp^{3}$ bond
failure length of 1.54 Å(at 100% strain) making conjugation along the inner
helix unlikely. Fig 7 shows the inner helix bond lengths at 100% strain at
which the bond failure length close to 1.54 Å is reached. Bond lengths at 75%
and 125% strains are also shown as comparison (also see Fig S1 in
supplementary for additional bond lengths along inner helix).

A different explanation of lower gap upon elongation compared to Ref. [13] can
be provided in terms of the outer helix bond lengths in strained helicenes. An
$sp^{2}$ carbon in an alkene is 1.34 Å in length, an $sp^{3}$ carbon is 1.46 Å
in length; an $sp^{2}$ benzene bond is 1.40 Å in length. The inner helix of
the helicene has bond lengths of 1.45 Å similar to $sp^{3}$ carbon and is thus
not very aromatic or conjugated. The outer helix has alternating bond lengths
of 1.36 Å and 1.43 Å and are more conjugated and aromatic than the inner helix
[18]. A decrease in the HOMO-LUMO gap implies better conjugation of bonds.
However, one would expect bond lengths to increase during tensile straining
greater than 25% lowering the strength of bonds and affecting conjugation in
the structure. On inspection of the structures, it is seen that length of
longer bonds in the outer helix decreases towards 1.40 Å with increasing
strain. Figure 8 shows the decrease of certain outer helix bond lengths,
indicated with arrows, towards 1.4 Å with straining and the accompanying
change in the distribution of the highest occupied orbital wavefunctions. Fig.
6(c)(bottom) depicts the HOMO and LUMO orbitals for helicene under tension (at
100% strain) showing that orbitals continue to localize with improved
conjugation along the outer helix compared to the 25% strain case.
Equalization of these outer helix bond lengths towards the $sp^{2}$ benzene
bond length contributes to lower band gap through improved conjugation in the
outer helix.

Additional analysis of the change in electronic transitions upon straining has
been included in Table 2. Although the change in the two lowest energy
excitations is small across the large change in strains, the $\beta$ peak
(third entry) shows a larger sensitivity, with a red shift seen under
compression and a blue shift under tension. This can be attributed to the
improved out-of-plane conjugation observed under compression, lowering the
excitation energy. The strength of the $\beta$ peak is also seen to
significantly increase under tension allowing differentiation of the strain
state.

Figure 7: HOMO and LUMO orbitals (together) are indicated at three strain
levels around failure strain.

Figure 8: Change of few outer helix bond lengths towards 1.4 Å, the $sp^{2}$ benzene bond length, at higher strains tends to improve conjugation along these bonds. The bonds indicated with arrows. Strain levels indicated below each case. Strain | %contribution | CI parentage | Wavelength (nm) | Oscillator strength | Rotatory strength
---|---|---|---|---|---
0 | 57,30 | 1-1*,2-2* | 404.4 | 0 | -0.4
| 88,8 | 2-1*,1-2* | 395.55 | 0.04 | 398
| 60,33 | 1-2*,3-1* | 381.2 | 0.02 | 271
-40 | 74,20 | 1-1*,3-2* | 404.82 | 0 | -20
| 83,13 | 2-1*,1-2* | 403.64 | 0 | 22
| 92,2 | 3-1*,1-2* | 398.4 | 0.03 | 358
50 | 54,37 | 1-1*,2-2* | 407.97 | 0 | -0.01
| 66,27 | 2-1*,1-2* | 387.7 | 0.03 | 116
| 64,27 | 1-2*,2-1* | 369.72 | 0.3 | 1197
Table 2: Calculated spectral data for [9]-helicene as a function of strain
(other entries similar to table 1). Only the top two main contributing
configurations are indicated.

## Discussion

The spring-like configuration of Helicene molecules leads to a smaller HOMO-
LUMO gap than that of its planar phenacene analogs. As the length of the
helicene chain increases, the conjugated network on the molecule becomes more
spread out, and hence, the transition from HOMO to LUMO requires progressively
lower energy, following Kuhn’s equation. The computed optical absorption
spectra of selected helicene molecules are found to compare well against
published experimental data. The spectra contain two prominent peaks near the
HOMO–LUMO gap called the para ($p$) band and the $\beta$ band. The helicene
peaks shift to lower energy with respect to its phenacene analogue as a result
of increased overlapping of atoms due to the twisting.

When helicenes are strained, the changes in the HOMO-LUMO gap followed a non-
symmetric parabolic trend with the amount of straining. Compressed helicenes
have a lower gap due to increased conjugation along the helical axis due to
chain overlap. The inflection in the HOMO–LUMO gap with tensile strain can be
explained as a transition from increasing HOMO-LUMO gap due to decreased
overlap to decreasing HOMO-LUMO gap due to improved conjugation as few outer
helix bonds reach the $sp^{2}$ benzene bond length during stretching. The
observed parabolic trend poses an issue for strain sensing as the fact that
for the same measured optical gap, there are two possible strain states of
opposite manner (either compression or elongation). This can be seen in the
optical absorption spectra plotted for different strain states in Fig. 9(a)
where selected compressive and tensile strain states show similar spectra near
the optical gap. However, over a larger range of wavelengths, compressive
states show lower absorption at higher optical energies (lower wavelengths),
resolving the issue of differentiating the sign of strains. Fig. 9(b) and Fig
10(a) compare the absorption spectra (linear scale) for tension versus
compression states. It is seen that the para ($p$-band) and $\beta$ peaks are
decreased due to the overlap of orbitals in compression states. Similar
behavior is also seen in the ECD spectra (Fig. 10b). Such a feature would
allow for a higher sensitivity measurement of strain states, coupled with
optical gap measurements.

(a) (b)

Figure 9: (a) Optical absorption spectra for selected strain states (a)
comparison at wavelengths near the optical gap (b) Comparison over a larger
set of wavelengths for 50 % tension and 40% compression.

(a) (b)

Figure 10: (a) Optical absorption compared for 10% compression or tension
states (b) Comparison of the ECD spectra for these strain states.

These insights provide promise for helicenes as a potential strain sensing
molecule when used in composite materials. Additional investigations on how
the optical gap as well as optical spectra gets modulated with respect to
different helicene parameters such as helicene length, heteroatom
functionalization (S, N, P atoms), ring substituents, and electronic doping,
would be valuable in highlighting their capability for future strain sensing
applications.

## Methods

All structures were generated and optimized using COMPASS force-field in
Forcite code [32] to be used as starting structures for first-principle
calculations. Preliminary band gap calculations for various helicenes (as a
right handed helix) as well as their planar phenacene and acene analogs (see
Fig. 1) were carried out using DFT calculations with local density
approximation with Vosko-Wilk-Nusair[33](LDA-VWN) exchange correlation. In
order to generate optimized geometries of various helicene structures for
TDDFT calculations, Hartree-Fock (HF) theory was employed as it provides exact
exchange correlation. Furche et al (2000) [29] presented the first TDDFT
simulations of circular dichroism spectra of [4]-[7] carbohelicenes using
geometry relaxed using Hatree-Fock theory. They found that DFT calculations
with GGA and hybrid functionals lead to inferior results through comparison
with experimental molecular structure. Double numeric basis sets with
polarization functions were used with an energy convergence criteria of 1.0e-4
Ha to obtain the optimized structures of different helicenes studied here.

Optical absorption simulations of helicenes were carried out using TDDFT
method in Gaussian-09 code (frequency space implementation is described in
detail in Ref. [34]). The first 25 excitation frequencies were computed to
model the optical spectra. The B3LYP exchange correlation functional has been
used based on its past success in modeling polyaromatic hydrocarbons[35, 36,
37]. In the context of carbohelienes, Nakai (2012) [30] performed a comparison
of different exchange functionals and found that the amount of exact exchange
in the DFT functionals plays an important role. B3LYP was found to reasonably
reproduce the experimental CD spectrum of [6] helicene, although slightly
underestimating the experimental CD intensities. The correlation consistent
polarized valence double zeta (ccPVDZ) basis set was chosen, which has been
successfully employed in the past for polyaromatic hydrocarbons [38] including
carbohelicenes[39] in combination with B3LYP exchange correlation functional.
The smoothed response was computed using a peak broadening parameter (half
width of the peak at half height) of 0.233 eV. To model the effect of helicene
deformation on optical properties, the helicenes were strained along the
helical axis direction by displacing the end aromatic carbon atoms of the
helical structure by a displacement corresponding to the needed strain level
(in both elongation and compression). Then, the structure was re-optimized
using HF theory (as detailed above) but with the two end-atoms fixed. The
optical absorption spectra of the strained structure were then computed using
TDDFT.

## References

  * [1] Urbano, A.  Recent developments in the synthesis of helicene-like molecules.  _Angewandte Chemie International Edition_ 42, 3986–3989 (2003).
  * [2] Gingras, M.  One hundred years of helicene chemistry. part 3: applications and properties of carbohelicenes.  _Chemical Society Reviews_ 42, 1051–1095 (2013).
  * [3] Narcis, M. J. & Takenaka, N.  Helical-chiral small molecules in asymmetric catalysis.  _European Journal of Organic Chemistry_ 2014, 21–34 (2014).
  * [4] Yavari, K. _et al._ Helicenes with embedded phosphole units in enantioselective gold catalysis.  _Angewandte Chemie_ 126, 880–884 (2014).
  * [5] Shen, Y. & Chen, C.-F.  Helicenes: synthesis and applications.  _Chemical reviews_ 112, 1463–1535 (2012).
  * [6] Dhbaibi, K. _et al._ Chiral diketopyrrolopyrrole-helicene polymer with efficient red circularly polarized luminescence.  _Frontiers in chemistry_ 8, 237 (2020).
  * [7] Xu, Y. _et al._ (p)-helicene displays chiral selection in binding to z-dna.  _Journal of the American Chemical Society_ 126, 6566–6567 (2004).
  * [8] Nuckolls, C. _et al._ Circular dichroism and uv- visible absorption spectra of the langmuir- blodgett films of an aggregating helicene.  _Journal of the American Chemical Society_ 120, 8656–8660 (1998).
  * [9] Wigglesworth, T. J., Sud, D., Norsten, T. B., Lekhi, V. S. & Branda, N. R.  Chiral discrimination in photochromic helicenes.  _Journal of the American Chemical Society_ 127, 7272–7273 (2005).
  * [10] Yang, Y., Da Costa, R. C., Fuchter, M. J. & Campbell, A. J.  Circularly polarized light detection by a chiral organic semiconductor transistor.  _Nature Photonics_ 7, 634–638 (2013).
  * [11] Fuchter, M. J. _et al._ [7]-helicene: a chiral molecular tweezer for silver (i) salts.  _Dalton Transactions_ 41, 8238–8241 (2012).
  * [12] Shi, L. _et al._ Synthesis, structure, properties, and application of a carbazole-based diaza [7] helicene in a deep-blue-emitting oled.  _Chemistry–A European Journal_ 18, 8092–8099 (2012).
  * [13] Guo, Y.-D., Yan, X.-H., Xiao, Y. & Liu, C.-S.  U-shaped relationship between current and pitch in helicene molecules.  _Scientific reports_ 5, 16731 (2015).
  * [14] Vacek, J., Chocholoušová, J. V., Stará, I. G., Starỳ, I. & Dubi, Y.  Mechanical tuning of conductance and thermopower in helicene molecular junctions.  _Nanoscale_ 7, 8793–8802 (2015).
  * [15] Neuman, K. C. & Nagy, A.  Single-molecule force spectroscopy: optical tweezers, magnetic tweezers and atomic force microscopy.  _Nature methods_ 5, 491–505 (2008).
  * [16] Capitanio, M. & Pavone, F. S.  Interrogating biology with force: single molecule high-resolution measurements with optical tweezers.  _Biophysical journal_ 105, 1293–1303 (2013).
  * [17] Malloci, G., Cappellini, G., Mulas, G. & Mattoni, A.  Electronic and optical properties of families of polycyclic aromatic hydrocarbons: A systematic (time-dependent) density functional theory study.  _Chemical Physics_ 384, 19–27, DOI: 10.1016/j.chemphys.2011.04.013 (2011).
  * [18] Schulman, J. M. & Disch, R. L.  Aromatic character of [n] helicenes and [n] phenacenes.  _The Journal of Physical Chemistry A_ 103, 6669–6672 (1999).
  * [19] Kuhn, H.  A quantum-mechanical theory of light absorption of organic dyes and similar compounds.  _The Journal of chemical physics_ 17, 1198–1212 (1949).
  * [20] Shemella, P., Zhang, Y., Mailman, M., Ajayan, P. M. & Nayak, S. K.  Energy gaps in zero-dimensional graphene nanoribbons.  _Applied Physics Letters_ 91, 042101 (2007).
  * [21] Son, Y.-W., Cohen, M. L. & Louie, S. G.  Energy gaps in graphene nanoribbons.  _Physical review letters_ 97, 216803 (2006).
  * [22] Clar, E. & Stewart, D.  Aromatic hydrocarbons. lxiii. resonance restriction and the absorption spectra of aromatic hydrocarbons1.  _Journal of the American Chemical Society_ 74, 6235–6238 (1952).
  * [23] Newman, M. S., Darlak, R. S. & Tsai, L. L.  Optical properties of hexahelicene.  _Journal of the American Chemical Society_ 89, 6191–6193 (1967).
  * [24] Martin, R., Flammang-Barbieux, M., Cosyn, J. & Gelbcke, M.  1-synthesis of octa-and nonahelicenes. 2-new syntheses of hexa-and heptahelicenes. 3-optical rotation and ord of heptahelicene.  _Tetrahedron Letters_ 9, 3507–3510 (1968).
  * [25] Flammang-Barbieux, M., Nasielski, J. & Martin, R.  Synthesis of heptahelicene (1) benzo [c] phenanthro [4, 3-g ]phenanthrene.  _Tetrahedron Letters_ 8, 743 – 744 (1967).
  * [26] Buss, V. & Kolster, K.  Electronic structure calculations on helicenes. concerning the chirality of helically twisted aromatic systems.  _Chemical physics_ 203, 309–316 (1996).
  * [27] Vander Donckt, E., Nasielski, J., Greenleaf, J. & Birks, J.  Fluorescence of the helicenes.  _Chemical Physics Letters_ 2, 409–410 (1968).
  * [28] Weigang Jr, O. & Dodson, P. T.  Low-temperature circular dichroism of hexahelicene.  _The Journal of Chemical Physics_ 49, 4248–4250 (1968).
  * [29] Furche, F. _et al._ Circular dichroism of helicenes investigated by time-dependent density functional theory.  _Journal of the American Chemical Society_ 122, 1717–1724 (2000).
  * [30] Nakai, Y., Mori, T. & Inoue, Y.  Theoretical and experimental studies on circular dichroism of carbo [n] helicenes.  _The Journal of Physical Chemistry A_ 116, 7372–7385 (2012).
  * [31] Farajian, A. A., Yakobson, B. I., Mizuseki, H. & Kawazoe, Y.  Electronic transport through bent carbon nanotubes: nanoelectromechanical sensors and switches.  _Physical Review B_ 67, 205423 (2003).
  * [32] Accelrys.  Inc., materials studio software (2010).
  * [33] Vosko, S. H., Wilk, L. & Nusair, M.  Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis.  _Canadian Journal of physics_ 58, 1200–1211 (1980).
  * [34] Stratmann, R. E., Scuseria, G. E. & Frisch, M. J.  An efficient implementation of time-dependent density-functional theory for the calculation of excitation energies of large molecules.  _The Journal of chemical physics_ 109, 8218–8224 (1998).
  * [35] Cardia, R., Malloci, G., Mattoni, A. & Cappellini, G.  Effects of tips-functionalization and perhalogenation on the electronic, optical, and transport properties of angular and compact dibenzochrysene.  _The Journal of Physical Chemistry A_ 118, 5170–5177 (2014).
  * [36] Cardia, R. _et al._ Electronic and optical properties of hexathiapentacene in the gas and crystal phases.  _Physical Review B_ 93, 235132 (2016).
  * [37] Dardenne, N. _et al._ Tuning optical properties of dibenzochrysenes by functionalization: A many-body perturbation theory study.  _The Journal of Physical Chemistry C_ 121, 24480–24488 (2017).
  * [38] Martin, J. M.  The vibrational spectra of corannulene and coronene. a density functional study.  _Chemical physics letters_ 262, 97–104 (1996).
  * [39] Schaack, C. _et al._ Helicene monomers and dimers: chiral chromophores featuring strong circularly polarized luminescence.  _Chemistry–A European Journal_ 25, 8003–8007 (2019).

## Acknowledgements

This research was supported in part by the Air Force Research Laboratory
Materials and Manufacturing Directorate, through the Air Force Office of
Scientific Research Summer Faculty Fellowship Program, Contract Numbers
FA8750-15-3-6003 and FA9550-15-0001. This research was supported in part
through computational resources and services provided by Advanced Research
Computing at the University of Michigan, Ann Arbor.

## Author contributions statement

V.V., D.S. and R.M. conceived the research idea(s), V.S. conducted the
simulations(s), V.V., V.S. and D.S. analysed the results. All authors reviewed
the manuscript.

## Competing interests

The author(s) declare no competing interests.

## Supplementary Data

A decrease in the HOMO-LUMO gap implies better conjugation of bonds. However,
one would expect bond lengths to increase during tensile straining greater
than 25% lowering the strength of bonds and affecting conjugation in the
structure. In Figure S1, alternating inner helix bonds are indeed highly
deformed during tension, with a few bonds stretching towards the $sp^{3}$ bond
failure length of 1.54 Å(at 100% strain) making conjugation along the inner
helix unlikely.

(a) (b)

Figure S1: Failure strain in [9]-helicene: Bond lengths in the inner helix
upon stretching. Close to 100% strain, a bond length of 1.54 Å(failure point
of sp3 bonds) is reached.