# Non-classicality of molecular vibrations activating electronic dynamics at
room temperature

Edward J. O’Reilly  Alexandra Olaya-Castro<EMAIL_ADDRESS>Department of
Physics and Astronomy, University College London, Gower Street, London, WC1E
6BT United Kingdom.

###### Abstract

Advancing the debate on the possible advantages of quantum effects in
biomolecular systems requires illustrations of biomolecular prototypes that
can exhibit and use non-classical phenomena to aid performance. Here we
investigate non-classical fluctuations in molecular dimers present in
photosynthetic light-harvesting antennae. We show that exciton-vibrational
resonances enable non-classical occupation fluctuations of the driven
collective high-energy vibration at room temperature. This phenomenon is a
manifestation of strong modulation of low-phonon occupation levels through the
transient formation of hybrid exciton-vibration states. For the prototype
considered, direct quantitative relations seem to hold between non-
classicality and optimal exciton population transfer.

The experimental investigation of quantum effects in photosynthetic energy
transfer Engel07 ; Panitchayangkoon07 ; Calhoun09 ; Collini10 ; Harel12  has
inspired exciting research towards asserting the relevance of non-classical
behaviour in the biological function of the molecular components of living
organisms. Theoretical works have made emphasis on the quantumness exhibited
by electronic degrees of freedom Whaley11  and on the key role of
environmental fluctuations Plenio08 ; Rebentrost09 . In order to understand
the full range of possible non-classical phenomena present during room-
temperature dynamics of electronic excitations in light-harvesting antennae it
is important to expand the analysis to the molecular motions that play a
fundamental role in driving exciton transport Renger01 . Research in this
direction is timely given that experiments Richards12 ; Turner12  and
theoretical works Womick11 ; Chin12a ; Chin12b ; Kreisbeck12 ; Kolli12 ;
Christensson12 ; Tiwari12  are highlighting the importance of concerted
electronic and vibrational dynamics in energy transfer. Moreover, recent
experimental techniques able to manipulate vibrational states Biggs09 ;
Biggs12  and probe their quantum properties may indeed provide experimental
access to the non-classicality of vibrational modes in biomolecular systems.
Here we consider a prototype molecular dimer with an structured exciton-phonon
interaction and investigate non-classical behaviour of molecular vibrations
that steer energy dynamics and exciton coherences in this molecular unit.
While previous studies have considered coherent exciton-vibration dynamics in
molecular dimersKuhn96 ; Chin12b ; Christensson12 , it is unclear how this
deviates from classical (wave-like) coherent behaviour and if so whether there
are quantitative relationships between non-classicality and effective energy
transfer.

The prototype of interest consists of two chromophores coupled by weak
electronic interactions and subject to the influence of both a low-energy
thermal background and a discrete high-frequency intramolecular vibration. The
key aspect of this molecular vibration is that its energy matches closely the
exciton gap but it is much larger than the thermal energy scale. Dimers of
this type naturally occur in several light-harvesting proteins Womick11 ;
Doust05  and are of particular interest as they represent systems where,
although counter-intuitive, exciton coherence beating has been probed Wong12 .
Insights into the importance of resonances between exciton energy mismatch and
vibrational frequencies for efficient energy transfer can be gained from
Förster theory Forster59 . However, the wider implications of such resonances
for optimal spatio-temporal distribution of energy Kolli12 ; Perdomo10 , for
manipulation of exciton coherence Chin12a ; Kolli12 ; Chin12b  and for non-
adiabatic excitation dynamics Tiwari12  have only started to be clarified. An
issue that remains to be tackled, is the extent to which such resonances can
expose the non-classical features of the molecular motions activating
transport. In the dimer investigated here, commensurate energies of exciton
gap and discrete vibrations imply that after input of exciton population,
effective energy distribution proceeds by driving the collective relative
displacement mode out of thermal equilibrium towards a non-classical state.
Non-classicality of the state of this relative motion is unambiguously
indicated by sub-Poissonian occupation statistics Mandel79  or alternatively
by strong modulation of low-phonon occupations Klyshko96  and is a
manifestation of transient collective coherent behaviour between excitons and
vibrations. Low-energy thermal fluctuations stabilize energy distribution to
the lowest exciton state before destroying non-classicality of the collective
mode assisting the process. Meanwhile, coherent interaction between electronic
degrees of freedom and such a mode manifests itself in the beating of
excitonic coherences with frequency components that vary depending on the
coupling to the thermal bath. When time-averaging over a relevant time scale,
we find direct quantitative relations between population transfer and the
degree of non-classicality as indicated by the average of negative values of
the Mandel $Q-$parameter Mandel79 . This relationship holds even when the
electronic system is initialized in a statistical mixture of exciton states.
Hence, efficient photosynthetic energy transport could take advantage of non-
classicality.

The field of quantum optics has developed a solid framework to quantify the
non-classicality of bosonic fields Davidovich96 . It therefore provides useful
conceptual and quantitative tools to investigate non-classicality of the
harmonic vibrational degrees of freedom of interest in this work. From the
perspective of quantum optics, non-classicality arises if the state of the
system of interest cannot be expressed as a statistical mixture of coherent
states defining a valid probability measure, this then leads to non-positive
or highly singular values of the Glauber-Sudarshan $P$-function or the Wigner
functionScully97 . However, experimental characterization of these
distributions is not always a straightforward task. Alternatively, it is of
interest to investigate signatures of non-classicality. For a single-mode
field, a sub-Poissonian occupation number distribution quantified by the
Mandel’s $Q-$parameter Mandel79  is a signature of deviations from classical
behaviour. The Mandel $Q-$parameter characterizes the departure of the
occupation number distribution $P(n)$ from Poissonian statistics through the
inequality

$\displaystyle
Q=\frac{\langle\hat{n}^{2}\rangle-\langle\hat{n}\rangle^{2}}{\langle\hat{n}\rangle}-1<0$
(1)

where $\langle\hat{n}\rangle$ and $\langle\hat{n}^{2}\rangle$ denote the first
and second moments of the bosonic number operator $\hat{n}$ respectively.
Vanishing $Q$ indicates Poissonian number statistics where the mean of
$\hat{n}$ equals its variance (such is the case of a coherent state of light).
For a chaotic thermal state one finds that $Q=\langle\hat{n}\rangle>0$
indicating that particles are ‘bunched’. The departure from classical (wave-
like) behaviour is characterised by $Q<0$ or particle number distributions
where the variance is less than the mean, as is the case for Fock states. The
occupation fluctuations associated with a sub-Poissonian distribution
originate from the quantized nature of the field and cannot be described in
terms of an underlying _classical_ stochastic process. A criteria that
directly connects non-classical occupation fluctuations and the negativity of
a quasiprobability distribution is an inequality indicating modulation of
nearest number occupationsKlyshko96 . In this work we use the above framework
to quantify non-classicality of the undamped vibrational modes supporting
exciton transport in the molecular dimer of interest. We consider each
chromophore with an excited electronic state of energy $\varepsilon_{i}$
strongly coupled to a quantized vibrational mode of high-frequency
($\omega_{\textrm{vib}}>K_{B}T$) described by the bare Hamiltonians
$H_{\textrm{el}}=~{}\sum_{i=1,2}\varepsilon_{i}\sigma^{+}_{i}\sigma^{-}_{i}$
and
$H_{\textrm{vib}}=~{}\omega_{\textrm{vib}}(b_{1}^{\dagger}b_{1}+~{}b_{2}^{\dagger}b_{2})$
respectively. Inter-chromophore coupling is generated by dipole-dipole
interactions of the form
$H_{\textrm{d-d}}=V(\sigma^{+}_{1}\sigma^{-}_{2}+\sigma^{+}_{2}\sigma^{-}_{1})$.
The electronic excited states interact with their local vibrational
environments with strength $g$, linearly displacing the corresponding mode
coordinate, $H_{\textrm{el-
vib}}=g\sum_{i=1,2}\sigma^{+}_{i}\sigma^{-}_{i}(b^{\dagger}_{i}+b_{i})$. In
the above $b^{\dagger}_{i}$($b_{i}$) creates (annihilates) a phonon of the
vibrational mode of chromophore $i$ while $\sigma^{\pm}_{i}$ creates or
destroys an electronic excitation at site $i$. The eigenstates $|X\rangle$ and
$|Y\rangle$ of $H_{\textrm{el}}+H_{\textrm{d-d}}$ denote exciton states with
energy splitting given by $\Delta
E=\sqrt{(\varepsilon_{1}-\varepsilon_{2})^{2}+4V^{2}}$. Transformation into
collective mode coordinates shows that only the mode corresponding to the
relative displacement of the modes with the bosonic operators
$b^{({\dagger})}=(b^{({\dagger})}_{1}-b^{({\dagger})}_{2})/\sqrt{2}$, couples
to the excitonic system. The centre of mass mode decouples from the electronic
degrees of freedom. The occupation of the relative displacement mode of
frequency $\omega_{\textrm{vib}}$ is denoted by $n$ and the non-classical
properties of such mode is the subject of interest in this article.

Our prototype system corresponds to the central dimer of the phycoerythrin 545
(PE545) antennae found in cryptophyte marine algae Doust04 . Exciton coherence
beating involving this dimer has been probed Collini10  and a structured
spectral density has been phenomenologically determined via linear
spectroscopy Novoderezhkin10 . The electronic parameters for the dimer taken
from reference Novoderezhkin10  are $\varepsilon_{1}=19574\textrm{cm}^{-1}$,
$\varepsilon_{2}=18532\textrm{cm}^{-1}$ and $V=92.2\textrm{cm}^{-1}$. The
intramolecular mode that is separated from the set of vibrational degrees of
freedom has frequency $\omega_{\textrm{vib}}=1111\textrm{cm}^{-1}$ being
quasi-resonant with the exciton energy splitting $\Delta
E=1058.2\textrm{cm}^{-1}$. The coupling strength of each chromophore to the
this mode is $g=\omega_{\textrm{vib}}\sqrt{0.0578}=267.1\textrm{cm}^{-1}$
obtained from Doust04 ; Novoderezhkin10 , and the thermal energy scale at room
temperature is $k_{\rm B}T\approx 200\textrm{cm}^{-1}$. We consider the
influence of the low energy thermal fluctuations via a thermal bath $H_{B}$
described by a continuous distribution of harmonic modes, each linearly and
independently coupled to the excited electronic states $H_{\textrm{el}}$ with
coupling strength characterised by a Drude spectral density
$J(\omega)=(\lambda\Omega_{c}\omega/\pi)/{(\omega^{2}+{\Omega_{c}}^{2}})^{2}$
with reorganization energy $\lambda=110\textrm{cm}^{-1}$ and cut-off frequency
$\Omega_{c}=100\textrm{cm}^{-1}$ Hoda12 , hence the strength of the system-
bath interaction is parametrized by $\sqrt{\lambda\Omega_{c}}$ Prior10 . Since
$J(\omega_{\textrm{vib}})\ll g^{2}$, this does not over-account for the
influence of the resonant mode on the electronic system. Importantly, although
the electronic coupling $V$ is the smallest energy scale in the prototype
dimer, the exciton dynamics and phonon fluctuations here investigated will not
be properly described by Förster theory. In order to accurately account for
the effects of the thermal bath we have extended the hierarchical expansion of
the dynamics Tanimura89 ; Ishizaki09 ; Shi09  to explicitly include the
quantum interactions between electronic excitations and the selected mode.
This approach allows us to account accurately for the effects of the low-
energy thermal bath. In calculating dynamics we consider an initial state, in
which the upper energy exciton eigenstate is rapidly excited while both high-
energy quantized vibrations and low-energy thermal modes are initially in
thermal equilibrium at room temperature. Truncating the quantized modes at an
occupation Fock level $M=6$ adequately describes both reduced dynamics of the
collective quantized mode and the electronic dynamics. A description of the
numerical methodology can be found in the appendix.

Figure 1:  a, Representation of the pigments and protein environment of the
PE545 complex (Protein Data Bank ID code 1XG0Doust04 ). The central dimer is
highlighted. Also shown are the energy levels of the exciton-vibration states
used to describe energy transfer $\left|{X,n}\right\rangle$. b, Coherent
dynamics of the exciton-vibration dimer. Oscillatory population of the lowest
exciton population $\rho_{YY}(t)$ (thick blue curve) is accompanied by
oscillations of the inter-exciton coherence in the ground-state of the
vibrational mode $|\rho_{X0-Y0}(t)|$ (thin red curve). c, Mandel $Q$-parameter
of the vibrational mode for biological electronic coupling and d for zero
electronic coupling.

Hamiltonian dynamics and non-classicality.- We first illustrate the quantum
coherent dynamics of excitation when coupled to the undamped quasiresonant
high-energy vibrational modes in the absence of the low-energy thermal bath.
This allows us to identify relevant time scales in the evolution of a state
with no initial excitonic superpositions:
$\rho(t_{0})=\left|{X}\right\rangle\langle
X|\otimes\varrho^{\textrm{th}}_{\textrm{vib}}$ which written in the basis of
exciton-vibration states of the form $|X,n\rangle$ becomes
$\rho(t_{0})=\sum_{n}P_{th}(n)\left|{X,n}\right\rangle\langle X,n|$. Here $n$
denotes the phonon occupation number of the collective mode while $P_{th}(n)$
denotes the thermal occupation of such level. We investigate the population of
the lowest excitonic state $\rho_{YY}(t)=\sum_{n}\langle
Y,n|\rho(t)|Y,n\rangle$, the absolute value of the coherence
$\rho_{X0-Y0}(t)=\langle X,0|\rho(t)|Y,0\rangle$ which denotes the inter-
exciton coherence in the ground-state of the collective vibrational mode, and
the non-classicality given by negative values of $Q(t)$. Hamiltonian evolution
generates coherent transitions from states $|X,n\rangle$ to $|Y,n+1\rangle$
(see Figure 1(a)) with a rate $f$ that depends on the exciton delocalization
($|V|/\Delta\epsilon$), the coupling to the mode $g$, and the phonon
occupation $n$ i.e. $f\simeq g(2|V|/\Delta\varepsilon)\sqrt{(n+1)/2}$. Since
$\omega_{\textrm{vib}}\gg K_{\rm B}T$ the ground state of the effective mode
is largely populated, such that the Hamiltonian evolution of the initial state
is largely dominated by the evolution of the state $|X,0\rangle$. This implies
that the energetically close exciton-vibration state $|Y,1\rangle$ becomes
coherently populated at a rate $f\simeq g(2|V|/\Delta\epsilon)$, leading to
the oscillatory pattern observed in the probability of occupation
$\rho_{YY}(t)$ as illustrated in Figure 1 (b). The low-frequency oscillations
of the dynamics of $\rho_{YY}(t)$ cannot be assigned to the exciton or the
vibrational degrees of freedom alone as expected from quantum-coherent
evolution of the exciton-plus-effective mode system. For instance, if the mode
occupation is restricted to at most $n=1$, the period of the amplitude of
$\rho_{YY}(t)$ is given approximately by the inverse of $(1/2)(\sqrt{(\Delta
E-\alpha)^{2}+(2g^{2}4V^{2}/\alpha\Delta E)}-(\Delta E+\alpha))$ with
$\alpha^{2}=2g^{2}+\omega_{\textrm{vib}}^{2}$ and
$(2g^{2}4V^{2}/\alpha^{2}\Delta E^{2})\ll 1$. Coherent exciton population
transfer is then accompanied by coherent dynamics of the inter-exciton
coherence $|\rho_{X0-Y0}(t)|$ with the main amplitude modulated by the same
low-frequency oscillations of $\rho_{YY}(t)$ and a superimposed fast
oscillatory component of frequency close to $\omega_{\textrm{vib}}$ as can be
seen in Figure 1 (b). This fast driving component arises from local
oscillatory displacements: when $V\simeq 0$ the time evolution of each local
mode is determined by the displacement operator with amplitude
$\alpha(t)={2g}(1-\exp(-i\omega_{\textrm{vib}}t))/\omega_{\textrm{vib}}$
Palma96 . As the state $|Y,1\rangle\langle Y,1|$ is coherently populated, the
collective quantized mode is driven out of equilibrium towards a non-classical
state in which selective occupation of the first vibrational level takes
place. This manifests itself in sub-Poissonian phonon statistics as indicated
by negative values of $Q(t)$ shown in Figure 1(c). Similar phenomena have been
described in the context of electron transport in a nanoelectromechanical
system Merlo08 ; Cavaliere08 . Importantly, such sub-Poissonian statistics
arises only when the electronic interaction is finite. For comparison, Figure
1(d) shows that if $V=0$, a local excitation drives the mode towards a
classical thermal displaced state with super-Poissonian statistics ($Q(t)>0$).
In short, non-classicality of the relative displacement mode quasiresonant
with the excitonic transition arises through the formation of exciton-
vibration states.

Energy and coherence dynamics - Interplay between vibration-activated dynamics
and thermal fluctuations leads to two distinct regimes of energy transport.
For our consideration of weak electronic coupling, the coherent transport
regime is determined approximately by $\sqrt{\lambda\Omega_{c}}\leq
2g|V|/\Delta\epsilon$. Here population of the low-lying exciton state is
dominated by coherent transitions between exciton-vibration states and the
rapid short-time growth of $\rho_{YY}(t)$ in this regime can be traced back to
coherent transfer of population from $|X,0\rangle\langle X,0|$ to
$|Y,1\rangle\langle Y,1|$. At longer time scales thermal fluctuations induce
incoherent transitions from $|X,0\rangle\langle X,0|$ to $|Y,0\rangle\langle
Y,0|$ with a rate proportional to $\sqrt{\lambda\Omega_{c}}$, thereby
stabilizing population of $|Y,1\rangle\langle Y,1|$ to a particular value (see
Figure 2(a)). In contrast, for $\sqrt{\lambda\Omega_{c}}>2g|V|/\Delta\epsilon$
population transfer to $\rho_{YY}(t)$ is incoherent although the dynamics is
still accompanied by generation and coherent evolution of excitonic coherences
(see Figures 2(g) and 2(h)). In this incoherent transport regime
$\rho_{YY}(t)$ has a slow but continuous rise reflecting the fact that the
incoherent transitions from $|X,0\rangle\langle X,0|$ to $|Y,0\rangle\langle
Y,0|$ now contribute significantly to the transport. However, transfer to
$\rho_{YY}(t)$ is always more efficient with the quasiresonant mode than in
the situations where only thermal-bath induced transitions are considered (see
dashed lines in Figures 2(d) and 2(g)). The underlying reason is that the
system is transiently evolving towards a thermal configuration of exciton-
vibration states. Hence, in both regimes excitation transfer to the lowest
exciton state is achieved by transiently sharing energy with the relative
displacement mode. The transition from coherent to incoherent exciton
population dynamics is also marked by the onset of energy dissipation of the
exciton-vibration system as shown in figures 2(c),2(f) and 2(i). While
population transfer is coherent, energy dissipation into the thermal bath is
prevented (Fig. 2(c) and 2(f)).

The beating patterns of the coherence $\rho_{X0-Y0}(t)$ reveals the structured
nature of the exciton-phonon interaction and witnesses whether there is
coherent exciton-vibration evolution. Because of this the frequency components
of such oscillatory dynamics vary depending of the coupling to the thermal
bath as can be seen in Figures 2(b), (e) and (h) . In the coherent regime
$\rho_{X0-Y0}(t)$ follow exciton populations with the main amplitude modulated
by the same relevant energy difference between exciton-vibration states (see
Figures 2(b) and (e)). In contrast, in the regime of incoherent transport the
short-time oscillations of $\rho_{X0-Y0}(t)$ (between $t=0$ and $t=0.1$ps)
arise from electronic correlations due to bath-induced renormalization of the
electronic Hamiltonian Silbey84  and as such, the frequency of these short-
time oscillations correspond to transitions of renormalized excitons. As the
main amplitude is damped, this exciton coherence retains the superimposed
driving at a frequency $\omega_{\textrm{vib}}$. The dynamical features
presented in Figures 2(g) and 2(h) confirm previous findings based on
perturbative calculations Kolli12  and agrees with the time-scales of the
exciton coherence beating reported in the two-dimensional spectroscopy of
cryptophyte algae Wong12 .

Figure 2: The dynamics of $\rho_{YY}(t)$ (top row), $|\rho_{X0,Y0}(t)|$
(middle row) and Energy of the exciton vibration system (bottom row) for three
interaction strengths to the low-energy thermal bath: a-c:
$\lambda=6\textrm{cm}^{-1}$, d-f $\lambda=20\textrm{cm}^{-1}$ and g-i:
$\lambda=110\textrm{cm}^{-1}$

Non-classicality - Interaction with a large thermal environment leads to the
emergence of classicality. Indeed, as the coupling of the electronic degrees
of freedom to the thermal background is increased, the maximum non-
classicality of the collective mode as indicated by the most negative value of
$Q(t)$ is reduced. However, the time interval for which the mode remains in a
non-classical state is not a monotonic function of $\lambda$. For moderate
values of $\lambda$, the relative displacement mode spends longer periods in
states with non-classical fluctuations i.e. periods for which $Q(t)<0$ as seen
in Figures 3(a)-(b). Hence relatively medium coupling to a thermal background
assists selective population of state $|Y,1\rangle$ thereby stabilizing non-
classicality at a particular level. An important feature of the results
presented in Figures 3(a)-3(c) is therefore that non-classical features
surviving the picosecond time scale are present across the whole range of
thermal bath couplings. This is consistent with a slow decay of the exciton-
vibration coherence $\rho_{X0,Y1}$ (not shown). Interestingly, the non-
classical properties of the collective mode resemble non-classicality of
thermal states (completely incoherent states) that are excited by a single
quanta Zavatta07 . The non-classical properties predicted by $Q(t)$ also agree
with the non-classicality predicted by a parameter proposed by Klyshko
Klyshko96  as is described in appendix B.

Figure 3: a-c Dynamics of the Mandel $Q$-parameter for
$\lambda=6,20,110\textrm{cm}^{-1}$. Figure 4:  Time integrated averages of
exciton population $\rho_{YY}(t)$, inter-exciton coherent in the lowest
vibrational level $|\rho_{X0-Y0}(t)|$ and non-classicality as quantified by
$Q(t)\Theta[-Q(t)]$ (blue, red and green respectively) as functions of
coupling to the thermal background by fixing environment cut-off frequency
$\Omega_{c}=100\textrm{cm}^{-1}$ and varying reorganization energy $\lambda$.

Non-classicality is witnessed by exciton population dynamics. In order to
illustrate this, we now investigate quantitative relations between non-
classicality and energy transport by considering relevant integrated averages
in the time scale of coherent exciton-vibration dynamics denoted by $\tau$.
For the parameters chosen, this time scale is about half a picosecond and
corresponds to the time in which excitation energy would be distributed away
to other chromophores Kolli12 . The time integrated averages over $\tau$ are
defined as: $\langle
F[\rho(t)]\rangle_{\tau}=\frac{1}{\tau}\int_{0}^{\tau}\textrm{d}t~{}F[\rho(t)]~{},$
where $F[\rho(t)]$ corresponds to the exciton population $\rho_{YY}(t)$ and
the non-classicality of the collective quantized mode through periods of sub-
Poissonian statistics $Q(t)\Theta[-Q(t)]$ as functions of the coupling to the
bath $\lambda$. As shown in Figure 4 the average exciton population and non-
classicality follow a similar non-monotonic trend as a function of the
coupling to the thermal bath, indicating a direct quantitative relation
between efficient energy transfer in the time scale $\tau$ and the degree of
non-classicality. The appearance of maximal point in the average non-
classicality as a function of the system-bath coupling indicates that the
average quantum response of the relative displacement mode to the impulsive
electronic excitation, is optimal for a small amount of thermal noise.

The above results indicate that the degree of delocalization of the initial
exciton state is key to enabling non-classical fluctuations of the collective
quantized vibration. This suggests that statistical mixtures of excitons with
finite purity can still trigger such non-classical response. To illustrate
this we now consider mixed initial states of the form
$\rho(t_{0})=\varrho_{\textrm{ex}}\otimes\varrho^{\textrm{th}}_{\textrm{vib}}$
where $\varrho_{\textrm{ex}}=r\left|{X}\right\rangle\langle
X|-(1-r)\left|{Y}\right\rangle\langle Y|$ with $1/2\leq r\leq 1$. The
associated linear entropy quantifying the mixedness of the initial exciton
states is given by $S_{L}=2r(r-1)$. The time-averaged non-classicality (Figure
5 (a)) and average population transfer (Figure 5 (b)) follow similar
decreasing monotonic trends for mixed states. These findings then suggest
direct quantitative relationships between non-classicality and exciton
population transport in the relevant time scale.

Figure 5:  The time integrated averages of exciton population
$\rho_{YY}(t)-\rho_{YY}(0$) and non-classicality as quantified by
$Q(t)\Theta[-Q(t)]$ as functions of the linear entropy $S_{L}=2r(r-1)$ for
$0<r<1/2$. Curves are for different system-bath couplings.

Outlook.- Our work shows that in molecular dimers operating in the weak
electronic coupling regime, exciton-vibrational resonances can trigger a non-
classical response of a quasiresonant collective high-energy vibration thereby
taking advantage of quantum fluctuations to efficiently populate a target
exciton state. These non-classical fluctuations are observed for a variety of
initial exciton conditions including statistical mixtures suggesting that such
non-classicality may be triggered even under incoherent input of electronic
population. These quantum phenomena could be experimentally accessible by
transient coherent ultrafast phonon spectroscopy Waldermann08  which is
sensitive to low-phonon populations of high-energy vibrations Lee11 . The
framework here discussed could also be used to give more insight into the
possible non-classical response of vibrational modes in other chemical sensors
Turin96  which are conjectured to operate through weak electronic
interactionsBrookes07  to sense vibrational frequencies of the order of a
thousand wavenumbers Franco11 .

From a more general perspective, our work illustrates how exciton-vibrational
resonances allow for rich dynamics of energy and coherences thereby providing
a mechanism for adjusting the rates of entropy production and output energy.
In this context, the associated vibrational dynamics could be understood as an
internal quantum mechanism for thermodynamic control Erez08 . Understanding
how quantum phenomena affect the thermodynamics of energy conversion in
biomolecules could provide an unifying framework to quantify the advantages of
quantum phenomena in a variety of biomolecular transport processes
Novoderezhkin04 ; Hay12 .

Acknowledgements.- We thank Greg Scholes and Rienk van Grondelle for
discussions, and Avinash Kolli for support on numerical calculations in the
initial stages of this project. Financial support from the Engineering and
Physical Sciences Research Council (EPSRC UK) is gratefully acknowledged.

## References

  * [1] Gregory S. Engel, Tessa R. Calhoun, Elizabeth L. Read, Tae-Kyu Ahn, Tomàš Mancǎl, Yuan-Chung Cheng, Robert E. Blankenship, and Graham R. Fleming.  Evidence for wavelike energy transfer through quantum coherence in photosynthetic systems.  Nature, 446(7137):782–786, Apr 12 2007.
  * [2] Gitt Panitchayangkoon, Dugan Hayes, Kelly A. Fransted, Justin R. Caram, Elad Harel, Jianzhong Wen, Robert E. Blankenship, and Gregory S. Engel.  Long-lived quantum coherence in photosynthetic complexes at physiological temperature.  Proceedings of the National Academy of Sciences, 107(29):12766–12770, 2010.
  * [3] Tessa R. Calhoun, Naomi S. Ginsberg, Gabriela S. Schlau-Cohen, Yuan-Chung Cheng, Matteo Ballottari, Roberto Bassi, and Graham R. Fleming.  Quantum coherence enabled determination of the energy landscape in light-harvesting complex II.  The Journal of Physical Chemistry B, 113(51):16291–16295, 2009\.
  * [4] Elisabetta Collini, Cathy Y. Wong, Krystyna E. Wilk, Paul M. G. Curmi, Paul Brumer, and Gregory D. Scholes.  Coherently wired light-harvesting in photosynthetic marine algae at ambient temperature.  Nature, 463(7281):644–647, Feb 4 2010.
  * [5] Elad Harel and Gregory S. Engel.  Quantum coherence spectroscopy reveals complex dynamics in bacterial light-harvesting complex 2 (LH2).  Proceedings of the National Academy of Sciences, 2012.
  * [6] K. Birgitta Whaley, Mohan Sarovar, and Akihito Ishizaki.  Quantum entanglement phenomena in photosynthetic light harvesting complexes.  Procedia Chemistry, 3(1):152 – 164, 2011.  22nd Solvay Conference on Chemistry.
  * [7] M B Plenio and S F Huelga.  Dephasing-assisted transport: quantum networks and biomolecules.  New Journal of Physics, 10(11):113019, 2008.
  * [8] Patrick Rebentrost, Masoud Mohseni, Ivan Kassal, Seth Lloyd, and Alán Aspuru-Guzik.  Environment-assisted quantum transport.  New Journal of Physics, 11(3):033003, 2009.
  * [9] Thomas Renger, Volkhard May, and Oliver Kühn.  Ultrafast excitation energy transfer dynamics in photosynthetic pigment-protein complexes.  Physics Reports, 343(3):137–254, 2001.
  * [10] G. H. Richards, K. E. Wilk, P. M. G. Curmi, H. M. Quiney, and J. A. Davis.  Coherent vibronic coupling in light-harvesting complexes from photosynthetic marine algae.  The Journal of Physical Chemistry Letters, 3(2):272–277, 2012.
  * [11] Daniel B. Turner, Rayomond Dinshaw, Kyung-Koo Lee, Michael S. Belsley, Krystyna E. Wilk, Paul M. G. Curmi, and Gregory D. Scholes.  Quantitative investigations of quantum coherence for a light-harvesting protein at conditions simulating photosynthesis.  Phys. Chem. Chem. Phys., 14:4857–4874, 2012.
  * [12] Jordan M. Womick and Andrew M. Moran.  Vibronic enhancement of exciton sizes and energy transport in photosynthetic complexes.  The Journal of Physical Chemistry B, 115(6):1347–1356, 2011.
  * [13] A. W. Chin, S. F. Huelga, and M. B. Plenio.  Coherence and decoherence in biological systems: principles of noise-assisted transport and the origin of long-lived coherences.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 370(1972):3638–3657, 2012.
  * [14] A. W. Chin, J. Prior, R. Rosenbach, F. Caycedo-Soler, S. F. Huelga, and M. B. Plenio.  The role of non-equilibrium vibrational structures in electronic coherence and recoherence in pigment-protein complexes.  Nature Physics, 2013.
  * [15] Christoph Kreisbeck and Tobias Kramer.  Long-lived electronic coherence in dissipative exciton dynamics of light-harvesting complexes.  The Journal of Physical Chemistry Letters, 3(19):2828–2833, 2012\.
  * [16] Avinash Kolli, Edward J. O’Reilly, Gregory D. Scholes, and Alexandra Olaya-Castro.  The fundamental role of quantized vibrations in coherent light harvesting by cryptophyte algae.  Journal of Chemical Physics, 137(17), Nov 7 2012.
  * [17] Niklas Christensson, Harald F. Kauffmann, Tõnu Pullerits, and Tomàš Mancǎl.  Origin of long-lived coherences in light-harvesting complexes.  The Journal of Physical Chemistry B, 116(25):7449–7454, 2012.
  * [18] Vivek Tiwari, William K. Peters, and David M. Jonas.  Electronic resonance with anticorrelated pigment vibrations drives photosynthetic energy transfer outside the adiabatic framework.  Proceedings of the National Academy of Sciences, 2012.
  * [19] Jason D. Biggs and Jeffrey A. Cina.  Using wave-packet interferometry to monitor the external vibrational control of electronic excitation transfer.  The Journal of Chemical Physics, 131(22):224101, 2009.
  * [20] Jason D. Biggs and Jeffrey A. Cina.  Studies of impulsive vibrational influence on ultrafast electronic excitation transfer.  The Journal of Physical Chemistry A, 116(7):1683–1693, 2012.
  * [21] Oliver Kühn, Thomas Renger, and Volkhard May.  Theory of exciton-vibrational dynamics in molecular dimers.  Chemical Physics, 204(1):99 – 114, 1996.
  * [22] Alexander B. Doust, Ivo H. M. van Stokkum, Delmar S. Larsen, Krystyna E. Wilk, Paul M. G. Curmi, Rienk van Grondelle, and Gregory D. Scholes.  Mediation of ultrafast light-harvesting by a central dimer in phycoerythrin 545 studied by transient absorption and global analysis.  The Journal of Physical Chemistry B, 109(29):14219–14226, 2005\.
  * [23] Cathy Y. Wong, Richard M. Alvey, Daniel B. Turner, Krystyna E. Wilk, Donald A. Bryant, Paul M. G. Curmi, Robert J. Silbey, and Gregory D. Scholes.  Electronic coherence lineshapes reveal hidden excitonic correlations in photosynthetic light harvesting.  Nature Chemistry, 5(396-404), 2012.
  * [24] Th. Förster.  10th spiers memorial lecture. transfer mechanisms of electronic excitation.  Discuss. Faraday Soc., 27:7–17, 1959.
  * [25] Alejandro Perdomo, Leslie Vogt, Ali Najmaie, and Alàn Aspuru-Guzik.  Engineering directed excitonic energy transfer.  Applied Physics Letters, 96(9):093114, 2010.
  * [26] L. Mandel.  Sub-poissonian photon statistics in resonance fluorescence.  Opt. Lett., 4(7):205–207, Jul 1979.
  * [27] D.N. Klyshko.  Observable signs of nonclassical light.  Physics Letters A, 213(1-2):7 – 15, 1996.
  * [28] Luiz Davidovich.  Sub-poissonian processes in quantum optics.  Rev. Mod. Phys., 68:127–173, Jan 1996.
  * [29] M.O. Scully and S. Zubairy.  Quantum Optics.  Cambridge University Press, 1997.
  * [30] Alexander B. Doust, Christopher N.J. Marai, Stephen J. Harrop, Krystyna E. Wilk, Paul M.G. Curmi, and Gregory D. Scholes.  Developing a structure-function model for the cryptophyte phycoerythrin 545 using ultrahigh resolution crystallography and ultrafast laser spectroscopy.  Journal of Molecular Biology, 344(1):135 – 153, 2004.
  * [31] Vladimir I. Novoderezhkin, Alexander B. Doust, Carles Curutchet, Gregory D. Scholes, and Rienk vanÂ Grondelle.  Excitation dynamics in phycoerythrin 545: Modeling of steady-state spectra and transient absorption with modified redfield theory.  Biophysical Journal, 99(2):344 – 352, 2010.
  * [32] Hoda Hossein-Nejad, Alexandra Olaya-Castro, and Gregory D. Scholes.  Phonon-mediated path-interference in electronic energy transfer.  The Journal of Chemical Physics, 136(2):024112, 2012.
  * [33] Javier Prior, Alex W. Chin, Susana F. Huelga, and Martin B. Plenio.  Efficient simulation of strong system-environment interactions.  Phys. Rev. Lett., 105:050404, Jul 2010.
  * [34] Yoshitaka Tanimura and Ryogo Kubo.  Time evolution of a quantum system in contact with a nearly gaussian-markoffian noise bath.  Journal of the Physical Society of Japan, 58(1):101–114, 1989.
  * [35] Akihito Ishizaki and Graham R. Fleming.  Unified treatment of quantum coherent and incoherent hopping dynamics in electronic energy transfer: Reduced hierarchy equation approach.  The Journal of Chemical Physics, 130(23):234111, 2009.
  * [36] Qiang Shi, Liping Chen, Guangjun Nan, Rui-Xue Xu, and YiJing Yan.  Efficient hierarchical liouville space propagator to quantum dissipative dynamics.  The Journal of Chemical Physics, 130(8):084105, 2009.
  * [37] G. Massimo Palma, Kalle-Antti Suominen, and Artur K. Ekert.  Quantum computers and dissipation.  Proceedings: Mathematical, Physical and Engineering Sciences, 452(1946):567–584, 1996.
  * [38] M Merlo, F Haupt, F Cavaliere, and M Sassetti.  Sub-poissonian phononic population in a nanoelectromechanical system.  New Journal of Physics, 10(2):023008, 2008.
  * [39] F Cavaliere, G Piovano, E Paladino, and M Sassetti.  Phonon distributions of a single-bath mode coupled to a quantum dot.  New Journal of Physics, 10(11):115004, 2008.
  * [40] Robert Silbey and Robert A. Harris.  Variational calculation of the dynamics of a two level system interacting with a bath.  The Journal of Chemical Physics, 80(6):2615–2617, 1984.
  * [41] Alessandro Zavatta, Valentina Parigi, and Marco Bellini.  Experimental nonclassicality of single-photon-added thermal light states.  Phys. Rev. A, 75:052106, May 2007.
  * [42] F. C. Waldermann, Benjamin J. Sussman, J. Nunn, V. O. Lorenz, K. C. Lee, K. Surmacz, K. H. Lee, D. Jaksch, I. A. Walmsley, P. Spizziri, P. Olivero, and S. Prawer.  Measuring phonon dephasing with ultrafast pulses using raman spectral interference.  Phys. Rev. B, 78:155201, Oct 2008.
  * [43] K. C. Lee, M. R. Sprague, B. J. Sussman, J. Nunn, N. K. Langford, X.-M. Jin, T. Champion, P. Michelberger, K. F. Reim, D. England, D. Jaksch, and I. A. Walmsley.  Entangling macroscopic diamonds at room temperature.  Science, 334(6060):1253–1256, 2011.
  * [44] Luca Turin.  A spectroscopic mechanism for primary olfactory reception.  Chemical Senses, 21(6):773–791, 1996.
  * [45] Jennifer C. Brookes, Filio Hartoutsiou, A. P. Horsfield, and A. M. Stoneham.  Could humans recognize odor by phonon assisted tunneling?  Phys. Rev. Lett., 98:038101, Jan 2007.
  * [46] Maria Isabel Franco, Luca Turin, Andreas Mershin, and Efthimios M. C. Skoulakis.  Molecular vibration-sensing component in drosophila melanogaster olfaction.  Proceedings of the National Academy of Sciences, 108(9):3797–3802, 2011.
  * [47] Noam Erez, Goren Gordon, Mathias Nest, and Gershon Kurizki.  Thermodynamic control by frequent quantum measurements.  Nature, 452(7188):724 – 727, 2008.
  * [48] Vladimir I. Novoderezhkin, Andrey G. Yakovlev, Rienk van Grondelle, and Vladimir A. Shuvalov.  Coherent nuclear and electronic dynamics in primary charge separation in photosynthetic reaction centers: a redfield theory approach.  The Journal of Physical Chemistry B, 108(22):7445–7457, 2004.
  * [49] Sam Hay and Nigel S. Scrutton.  Good vibrations in enzyme-catalysed reactions.  Nature Chemistry, 4(3):161–168, Mar 2012.

##  Appendix A Hierarchical expansion of exciton-vibration dynamics

Here we provide details of the model for a system consisting of $N$ electronic
degrees of freedom each interacting with a quantized vibrational mode and the
method used to accurately calculate the dynamics of these exciton-vibration
systems. The system Hamiltonian
$H_{S}=H_{\textrm{el}}\otimes\mathbbm{1}_{\textrm{vib}}+\mathbbm{1}_{\textrm{el}}\otimes
H_{\textrm{vib}}+H_{\textrm{el-vib}}$, is defined by

$H_{\textrm{el}}=\sum_{i=1}^{N}\varepsilon_{i}\sigma^{+}_{i}\sigma^{-}_{i}+\sum_{i,j\leq
i}^{N}V_{ij}(\sigma^{+}_{i}\sigma^{-}_{j}+\sigma^{+}_{j}\sigma^{-}_{i})$ (2)
$H_{\textrm{vib}}=\sum_{i=1}^{N}\omega_{\textrm{vib}}^{i}b^{\dagger}_{i}b_{i}$
(3)

and

$H_{\textrm{el-
vib}}=g_{i}\sigma^{+}_{i}\sigma^{-}_{i}\otimes(b^{\dagger}_{i}+b_{i})~{}.$ (4)

Here, the operators $\sigma^{\pm}_{i}$ create/annihilate electronic
excitations at site $i$, while $b^{\dagger}_{i}$ ($b_{i}$) create (annihilate)
phonons associated with the vibrational mode of site $i$. $\varepsilon_{i}$
denotes the excitation energies and $V_{ij}$ the strength of interactions
between sites. $\omega_{\textrm{vib}}^{i}$ is the frequency of the vibrational
mode at site $i$ and $g_{i}=\sqrt{S_{i}}\omega_{\textrm{vib}}^{i}$ the
coupling of the mode to its site. $S_{i}$ is the Huang-Rhys factor or mean
number of phonons in the ‘polaron cloud’ formed when the mode is fully
displaced. The the high energy modes considered here, we expand the quantized
vibrational modes of the system in the basis of single phonon (Fock) states
(e.g. $b_{i}=\sum_{n}^{M}(n+1)^{1/2}\left|{n}\right\rangle\langle n+1|$)
truncating at $M$ phonons. This system is then coupled to a continuum of
bosonic modes (the thermal background)
$H_{B}=\sum_{\textbf{k}}\omega_{\textbf{k}}b_{\textbf{k}}^{\dagger}b_{\textbf{k}}$,
with an interaction

$H_{I}=\sum_{\textbf{k}}\sum_{i}g_{\textbf{k},i}(\sigma^{+}_{i}\sigma^{-}_{i}\otimes\mathbbm{1}_{\textrm{vib}})(b_{\textbf{k}}^{\dagger}+b_{\textbf{k}})~{},$
(5)

such that a mode with frequency $\omega_{\textbf{k}}$ is linearly coupled to
the electronic excitation at site $i$ with a strength $g_{\textbf{k,i}}$ as
specified by a spectral density
$J_{i}(\omega)=\sum_{\textbf{k}}|g_{\textbf{k,i}}|^{2}\delta(\omega-\omega_{\textbf{k}})$.
We assume the bath is spatially uncorrelated and identical for each site.

An infinite hierarchy of coupled differential equations[34, 35, 36] is used to
express the dynamics of the exciton-vibration system density matrix
$\rho_{S}(t)$. These read as,

$\displaystyle\partial_{t}\tilde{\rho}_{\boldsymbol{n}}(t)$
$\displaystyle=-i/{\hbar}[H_{S},\tilde{\rho}_{\boldsymbol{n}}(t)]-\sum_{i=1}^{N}\sum_{k=0}^{K}n_{ik}v_{ik}\tilde{\rho}_{\boldsymbol{n}}(t)$
$\displaystyle{-i\sum_{i=1}^{N}\left(\sum_{k=1}^{\infty}\frac{c_{ik}}{\nu_{ik}}-\sum_{k=1}^{K}\frac{c_{ik}}{\nu_{ik}}\right)\left[Q_{i},\left[Q_{i},\tilde{\rho}_{\boldsymbol{n}}(t)\right]\right]}$
$\displaystyle-i\sum_{i=1}^{N}\sum_{k=0}^{K}\sqrt{(n_{ik}+1)|c_{ik}|}\left[Q_{i},\;\tilde{\rho}_{\boldsymbol{n_{ik}^{+}}}\right]$
$\displaystyle-i\sum_{i=1}^{N}\sum_{k=0}^{K}\sqrt{n_{ik}/|c_{ik}|}\left(c_{ik}\,Q_{i}\tilde{\rho}_{\boldsymbol{n_{ik}^{-}}}-c_{ik}^{*}\,\tilde{\rho}_{\boldsymbol{n_{ik}^{-}}}Q_{i}\right)~{}.$
(6)

for the case of a Drude spectral density and in a rescaled formulation[36]
which admits their efficient adaptive time-step numerical integration.
Operators $\tilde{\rho}_{\boldsymbol{n}}$ are indexed by the multi-index
$\boldsymbol{n}$ which has entries $n_{ik}$ and $n=\sum_{ik}n_{ik}$.
$\boldsymbol{n_{ik}^{\pm}}$ is $\boldsymbol{n}$ with entry $n_{jk}\rightarrow
n_{jk}\pm 1$. The reduced density matrix of the system $\rho_{S}$ is given by
the operator $\tilde{\rho}_{\boldsymbol{0}}$. The coefficients $\nu_{ik}$ and
$c_{ik}$ are the Matsubara frequencies and coefficients appearing in the
(truncated) exponential decomposition $\sum_{k=0}^{K}c_{ik}e^{-\nu_{ik}t}$ of
the bath correlation function applied in Ref [35]. In the above terms $k>K$
are truncated with a Markovian approximation scheme [35]. This results in
accurate dynamics for $K=1$ at high temperatures ($\beta\Omega_{c}<1$). In the
present exciton-vibration case the system-bath coupling operator takes the
form $Q_{i}=\sigma_{i}^{+}\sigma_{i}^{-}\otimes\mathbbm{1}_{\textrm{vib}}$.

##  Appendix B $B_{0}$ parameter

Modulations of adjacent phonon number occupation as quantified by negative
$B_{n}$, guarantees negative regions of a quasi-probability distribution[27]
as is a potentially an experimentally accessible quantity. The non-
classicality indicated by sub-Poissonian statistics is confirmed by negative
$B_{0}=2P(0)P(2)-P(1)$. This arises from selective population of the $P(1)$
level of the relative displacement vibrational mode.

Figure 6: $B_{0}$ for a) $\lambda=0\textrm{cm}^{-1}$, b)
$\lambda=6\textrm{cm}^{-1}$, c) $\lambda=20\textrm{cm}^{-1}$ and d)
$\lambda=110\textrm{cm}^{-1}$


# Habits of Mind: Reusing Action Sequences for Efficient Planning

Noémi Éltető1<EMAIL_ADDRESS>& Peter Dayan1,2
<EMAIL_ADDRESS>
1Computational Neuroscience, MPI for Biological Cybernetics 2University of
Tübingen

###### Abstract

When we exercise sequences of actions, their execution becomes more fluent and
precise. Here, we consider the possibility that exercised action sequences can
also be used to make planning faster and more accurate by focusing expansion
of the search tree on paths that have been frequently used in the past, and by
reducing deep planning problems to shallow ones via multi-step jumps in the
tree. To capture such sequences, we use a flexible Bayesian action chunking
mechanism which finds and exploits statistically reliable structure at
different scales. This gives rise to shorter or longer routines that can be
embedded into a Monte-Carlo tree search planner. We show the benefits of this
scheme using a physical construction task patterned after tangrams.

Keywords: Planning; Sequence models; Bayesian nonparametrics; Action chunking;

## Introduction

The law of exercise holds that animals have a tendency to repeat or
perseverate on frequent choices, even when this reduces reward (?, ?, ?). This
leveraging of structure in behaviour complements the leveraging of structure
in the environment. ? (?) suggested that perseveration arises from a trade-off
between reward maximization and policy complexity, with repetition making
policies more compact. This theory was focused on isolated actions; however,
animals also pick up on sequences of action. The action chunking process that
underlies skill or habit learning (?, ?) has also been observed in multi-step
decision-making tasks, for instance with rats repeating whole trajectories of
choices (?, ?). A hallmark of such action chunks is that they are performed
quickly and uninterrupted by state assessment, that is, open-loop.

The law of exercise may apply not only to the execution of actions in multi-
step problems, but also to the planning thereof. ? (?) tested how human
participants simplified a deep planning problem. One of the three heuristics
they found was stochastic memoization, in which participants became
increasingly prone to repeat previous action paths in their entirety or chunks
of them, even though these might not be optimal. In their interpretation, this
marks a transfer from flexible but costly computation to the less flexible but
cheaper reliance on past experience. A similar observation was made in a task
where participants planned and constructed complex structures that were
required to be stable in a virtual environment simulating rigid-body physics
(?, ?). During training, participants’ construction procedures became
increasingly stereotyped.

The trouble with sequences of actions is that their number tends to grow
exponentially with the length of the sequence. This makes methods that try to
be comprehensive (?, ?) computationally infeasible. Indeed, we know of no
scalable or cognitively plausible model of action chunk reuse in planning.
Here, in keeping with a general trend towards internalized decision-making (?,
?), we examine the benefits of using a non-parametric Bayesian sequence model
that can infer sparse chunks of various lengths, depending on the statistics
of the sequences executed in past episodes. We use this model to augment
Monte-Carlo Tree Search (MCTS) in two ways: using single actions that it
predicts as being likely, to bias one-step tree expansions, and proposing
popular multi-step expansions. We test our model on a physical construction
task called the Sticky Tangram.

## MCTS-with-HABITS

To examine the benefits of biasing planning towards previously learned action
sequences, we used a Bayesian infinite sequence model (a hierarchical
Dirichlet process) to augment Monte Carlo Tree Search with open-loop, multi-
action, node expansions. We call the new planning algorithm MCTS-with-HABITS
(Figure 1).

Figure 1: Schematic of the MCTS-with-HABITS model. The planner builds a search
tree with nodes representing states and edges representing different possible
actions, marked by colors (left box). The planner traverses the tree by
choosing actions that ultimately lead to more simulated wins and that are more
predictable by the HABITS sequence model (upper right box). At the leaf state
(marked by the bold node) several potential actions (marked by dashed edges)
are considered that lead to states (marked by dashed nodes) that have not been
added to the tree yet. The red and the green actions are primitive actions
whose winning values (marked as pies) are the same. Yet, the green action has
a higher likelihood under the sequence model, when conditioned on the action
trace from the root to the leaf node. Therefore, the green action will be
preferred over red. The green action is also predictably followed, given the
context of the action trace, by the blue action. Therefore, the chunk
generator (lower right box) appends the blue action to the green one.
Conditioned on the previous action trace and the blue action, no further
action is strongly predicted and the chunk generation is terminated. Thus we
extend the available action set by the green-blue chunk. Note that the
primitive green action, from which we unrolled a chunk, is also kept as an
alternative. If the action chunk is selected by the planner then it jumps over
the state marked by a grey node, not considering the available actions from
that state, also marked by grey. This effective stunting results in a
different value estimate for the chunk green-blue (marked by the larger pie in
the node) compared to the primitive action green; in this example, the chunk
will be preferred by the tree policy.

Consider the search tree in Figure 1, representing a deterministic planning
problem in which nodes correspond to states and the edges correspond to
actions. In Monte Carlo Tree Search (MCTS; ? (?)), the tree is traversed from
the root node to a leaf by repeated action choices based on their simulated
expected value. This decision mechanism is called the tree policy. A leaf is
any node that has a potential child whose value has not been estimated. Once a
leaf is reached, it is expanded by considering the actions that are available
in that state and choosing one of them. Then, the value of the node is
estimated by simulating a random rollout using a uniform random policy. The
reward of the rollout (in our case of a puzzle, $1$ for a correct solution;
$0$ otherwise) is backpropagated to the selected node and its ancestors. The
four ingredients – selection, expansion, simulation, and backpropagation –
constitute one step of the MCTS algorithm. Since it estimates the value of
nodes via simulation, MCTS focuses the search on promising paths without the
need for handcrafted distance heuristics as in A*. In order to balance this
exploitative focus (from the fraction $\frac{w_{t}}{n_{t}}$ of wins $w_{t}$
from $n_{t}$ tries from a node on the $t^{\text{th}}$ MCTS step) with
exploration of potentially even better states, a bonus is added to actions
leading to lesser tried nodes ($\sqrt{\frac{lnN_{t}}{n_{t}}}$, where $N_{t}$
is the total number of simulations run from the parent node; as in the UCT
algorithm).

To this conventional tree policy, we add a habit value term, that is, the
scaled likelihood of the action $a$ planned to get to a node, making its total
propensity be:

$Q(a)=\dfrac{w_{t}}{n_{t}}+c*\sqrt{\dfrac{lnN_{t}}{n_{t}}}+h*p_{t}(a_{i}=a|a_{0}:a_{i-1})$
(1)

where $c$ is the exploration coefficient, here fixed to 1; $h$ is the habit
coefficient; $a_{i}$ is the action directly leading to the node; and
$a_{0}:a_{i-1}$ is the path of all preceding planned actions from the root,
where $i$ indicates the current depth in the tree. Note that the win
proportion and exploration terms are node-dependent, while the habit term is
node-independent. Probability $p_{t}(a_{i}=a|a_{0}:a_{i-1})$ comes from a
Bayesian non-parametric sequence model that flexibly combines the predictive
power of dependencies at variable depths (?, ?). The model was recently used
to explain motor skill learning in humans (?, ?)). Here, we adapt it to the
task of hierarchical state augmentation (HABITS) in which it provides a
flexible window onto the action path leading to a state, for determining the
likelihood of the next action in the search (upper right box of Figure 1). Its
predictions are based on previous action sequence outputs of the planner.
Formally (temporarily dropping the iteration index $t$), the model is an
hierarchical Dirichlet process (HDP):

$a_{i}\sim\mathbf{G}_{\mathbf{u}_{i}}\sim
HDP(\alpha,\mathbf{G}_{\pi(\mathbf{u}_{i})})$ (2)

where $\mathbf{G}_{\mathbf{u}_{i}}$ is the vector of action probabilities
given the context $\mathbf{u}_{i}$ of some previous actions before the $i$-th
action; $\alpha$ is a strength parameter, controlling the resemblance to the
base distribution and, therefore, the speed of sequence learning; and the base
distribution $\mathbf{G}_{\pi(\mathbf{u}_{i})}$, which is the probability
distribution over actions given the truncated context $\pi(\mathbf{u}_{i})$
containing all but the earliest action. Applying Equation 2 recursively
performs a weighted smoothing across the action probabilities conditioned on
preceding action chunks of shrinking sizes. For more details and algorithms,
see ? (?). In sum, the sequence model flexibly augments the states with a
weighted window onto the previous actions, efficiently utilizing the ones that
have predictive power.

Assume that the HABITS module has been trained on action sequence outputs of
the planner from previous episodes (irrespectively of the quality of those
plans). Then, the HABITS module interacts with the planner in two ways. It
allows for biasing the one-step search towards predictable action sequence
paths through the habit term in Equation 1. Moreover, by unrolling the
sequence model, we also extend the node’s action repertoire by action chunks
(lower right box in Figure 1). Action chunks were created by sampling and
attaching successor actions while the conditional entropy of the distribution
over the next action was lower than the open-loop threshold parameter
$\omega$. When the action chunks were considered by the tree policy, their
habit value was determined by the predictability of the first action in the
chunk. By choosing an action chunk, the agent jumped deeper in the tree
without adding the intermediate node(s) to the tree and evaluating them,
effectively stunting all the other branches that would have grown from the
intervening nodes. As such, the MCTS-with-HABITS can be viewed as a hybrid
between open- and closed-loop. The complexity of the chunks proposed by the
HABITS scaled with the complexity of the sequences that the agent executed in
past episodes, opportunistically saving as much planning cost as possible.

Our experiments tested the unique benefits of the two separate mechanisms: the
sequence bias in the one-step search, realized by the habit term in the tree
policy, and the open-loop action chunks, realized by unrolling the sequence
model until the open-loop entropy threshold $\omega$ is reached. Therefore, we
compared three variants of our model: MCTS-with-HABITSfull, utilizing both
one-step biases and open-loop action chunks; MCTS-with-HABITSopen-loop,
utilizing no one-step bias for the actions but having the alternative of
unrolling primitive actions into chunks; and MCTS-with-HABITSone-step having
no action chunk alternatives but utilizing sequence predictions to bias its
action choices; along with the control model, MCTSvanilla. For brevity, we
will refer to MCTS-with-HABITSfull as MCTS-with-HABITS. Within the scope of
this paper, we did not optimize the model parameters, but selected values for
demonstration (Table 1).

Table 1: Parameter settings for four alternative model variants. The node budget $b$, exploration coefficient $c$ and sequence learning parameter $\alpha$ were fixed across the model variants, while we varied the habit coefficient $h$ and the open-loop threshold $\omega$. Model name | $b$ | $c$ | $\alpha$ | $h$ | $\omega$
---|---|---|---|---|---
MCTS-with-HABITSfull | 50 | 1 | 1 | 5 | 1.5
MCTS-with-HABITSopen-loop | 50 | 1 | 1 | 0 | 1.5
MCTS-with-HABITSone-step | 50 | 1 | 1 | 5 | 0
MCTSvanilla | 50 | 1 | - | 0 | 0

Apart from the labelled additions to the tree policy, the MCTS algorithm was
used conventionally, looping through the node selection, expansion, rollout
with uniform random policy, and backpropagation of a binary value (success or
failure) to the selected node and its ancestors.

## Results

### Task

We tested our model on a physical construction task inspired by the classical
tangram puzzle, versions of which were recently employed by ? (?) and ? (?).
Our version is called the Sticky Tangram (Figure 2A). The agent had to plan a
construction that matched a target silhouette, using the seven building blocks
in the inventory. The blocks could be placed into a grid space which
discretized and constrained the state space. The coarse discretization is
essential in order to keep the puzzle in a moderate complexity regime. The
blocks had to be placed in such a sequence that they touched each other, i.e.
floating blocks were not allowed.

Figure 2: (A) Schematic of an example trial of the Sticky Tangram task. The
silhouette (grey area) is the target that should be built up using the
building blocks (shown in different colors). The first chosen block has to be
placed on the floor (yellow line); the following blocks have to be placed such
that they touch at least a unit segment of a previous block’s edge (hence the
’stickiness’). In the duplet condition, two blocks always occurred together,
in the same relative position; in the triplet condition, three blocks formed
such a chunk. Throughout the paper, these two particular chunks will be used
as examples. (B) The full search tree of the example problem in (A). The leaf
node in red brackets is the only solution, the other leaf nodes represent dead
ends.

A state included the target silhouette and the action history of the relative
positions of the previously placed building blocks (Figure 2B). Adding a block
to the construction constituted a primitive action; adding several blocks at
once constituted an action chunk. In any state, only valid actions were
considered – that is, block placements that were inside the borders of the
target silhouette and that resulted in no disjoint parts. Note that this
applied to action chunks as well, meaning that actions that were unrolled into
chunks that would have violated the rules of the game were discarded.

### Experiments

### Reusing variable-size chunks

We randomly generated silhouettes that were composed of four building blocks
(Figure 2A). In the duplet condition, two building blocks formed a chunk such
that they would always occur together, in the same relative position and had
to be placed in a fixed order. In the triplet condition, three formed such a
chunk. The unchunked blocks occurred in unpredictable relative positions.
’Chunky silhouettes’ comprised the duplet and two random unchunked building
blocks in the duplet condition or the triplet and one random unchunked
building block in the triplet condition.

In the chunky silhouettes, the marginal probabilites of the chunked blocks
were $1$, but in the unchunked blocks they were $<1$. That is because the
elements of the chunks were always present in the chunky silhouettes, per
definition, but the other blocks were not. Since we wanted to study the
advantage of sequence structure in planning problems, we ensured that the
overall marginal probabilities of the blocks were equal, and a chunked block
was more predictable only conditioned on the other elements of the chunk.
Therefore, we introduced ’random silhouettes’ comprising only unchunked
blocks. A 4:3 ratio of the chunky and random silhouettes in the training set
ensured that the marginal probability distribution of the blocks was uniform –
that is, the agent used all the blocks equally frequently for the correct
solutions.

Since the chunky silhouettes were generated using more constraints than the
random ones (i.e. the chunk constraint), on average, the full search trees of
chunky silhouettes are less complex. We quantified the problem complexity as
the median number of nodes that the MCTSvanilla evaluated until finding the
solution. Then we sampled sets of the chunky and random silhouettes such that
they were matched by their complexity. We restricted the problem complexity to
be $\leq 50$. The example silhouette shown in Figure 2A and B has a full
search tree of 18 nodes and a low complexity of 12.

We trained four MCTS variants on 19 trials (a randomized sequence of chunky
and random silhouettes). During training, all model variants used a budget of
50 nodes, enabling nearly perfect performance due to the silhouettes having
complexity values $\leq 50$ (Figure 3A and E). In the case of the three MCTS-
with-HABITS model versions, the sequence module was trained on the action
sequence output of the planner. We ran 32 simulations for both the duplet and
triplet conditions, from different random seeds, and averaged the results.

All model variants showed nearly perfect performance given a node budget
sufficient for the problems ($>96\%$ in the duplet condition and $>94\%$ in
the triplet condition; Figure 3A and E), albeit with different patterns of
chunk use. The sequence model of the MCTS-with-HABITS model versions gradually
extracted the structure in the action sequences generated by the planner, that
is, the action duplets and triplets in the respective conditions. The MCTS-
with-HABITS and MCTS-with-HABITSopen-loop models began to propose action
chunks and these chunks were used more often in the correct solutions
throughout the training, both in the duplet (effect of trial:
$F(1,696)=293.70$, $p<.001$; Figure 3B) and in the triplet conditions (effect
of trial: $F(1,712)=29.51$, $p<.001$; Figure 3F). In both conditions, the
MCTS-with-HABITS was more disposed to include chunks in its plans (model
effect: $F(1,696)=207.58$, $p<.001$ for the duplet condition;
$F(1,712)=150.39$, $p<.001$ for the triplet condition). Problems that were
solved in four primitive steps by the MCTSvanilla (Figure 3C and G) were
solved by the MCTS-with-HABITS in only three steps in the duplet condition
(Figure 3D) and in two steps in the triplet condition (Figure 3 H). Note that
action chunks were reused in spite of the same states never reoccurring across
episodes since the sequence model proposes chunks in a state-independent
manner. In sum, the MCTS-with-HABITS model found and utilized the largest
exploitable chunks that reoccurred in planned action sequences.

Figure 3: Emergence of action sequence reuse in the duplet condition (A-B) and
the triplet condition (C-D). Error bars and bands indicate the $95\%$ CI.
(A)(E) All model variants had nearly perfect performance both in the duplet
and triplet conditions. (B)(F) The two model versions that were permitted to
use open-loop action chunks gradually picked up on the structure in the action
sequences and became more likely to use an action chunk for the solution. In
the triplet condition (F), chunk use emerged even faster than in the duplet
condition (B). (C)(G) Examples of a correct solutions by the MCTSvanilla.
(D)(H) Examples of solutions by the MCTS-with-HABITS model, using an action
duplet (D) and an action triplet in (H). The black polygons mark chunks.

### Solving high complexity problems with less compute

After the training trials, we gradually reduced the node budget from the
initial value of 50 to 12, 8, 5, and 1. In order to compare performance across
these node budgets, we froze the sequence model at the test trials. At each
reduced node budget value, the agents solved four silhouettes. We ran 32
simulations for both the duplet and triplet conditions separately, from
different random seeds, and averaged the results.

When solving random silhouettes, all model variants’ performance dropped to
similar degrees as a function of node budget restriction (effect of budget:
$F(1,1872)=124.58$, $p<.001$; effect of model: $F(1,1872)=0.42$, $p<.73$;
Figure 4A and C). In the case of the chunky silhouettes, when previously
learned action sequences were appropriately reusable, the resource-constrained
model performance depended on whether the model was allowed to use one-step
biases, action chunks, or both; this held both in the duplet (effect of model:
$F(1,1244)=17.03$, $p<.001$; Figure 4B) and triplet conditions (effect of
model: $F(1,1188)=57.16$, $p<.001$; Figure 4D).

Figure 4: Model performance in the face of node budget restriction in the test
phase of Experiment 1, in the duplet condition (A-B) and in the triplet
condition (C-D). Error bars indicate the $95\%$ CI. (A)(C) Performance on
random silhouettes decreased for all model variants to similar degrees due to
the node budget restriction. (B)(D) For chunky silhouettes, the performance of
the MCTS-with-HABITSopen-loop and MCTS-with-HABITSone-step was more resilient
to the node budget restriction than that of the MCTSvanilla. The full MCTS-
with-HABITS model showed the best performance in the face of node budget
restriction.

We then analyzed the planning performance in the resource-restricted regime of
12, 8, and 5 nodes, as a function of problem complexity. In the case of the
random silhouettes, the different model variants’ performance scaled with the
problem complexity similarly, both in the duplet (model*complexity
interaction: $F(6,384)=1.72$, $p=.11$; Figure 5A) and triplet conditions
(model*complexity interaction: $F(6,364)=.52$, $p=.78$; Figure 5C). However,
performance on the chunky silhouettes scaled differently with complexity
across the model variants, in both conditions (model*complexity interaction in
duplet condition: $F(6,744)=1.89$, $p=.07$; Figure 5B; model*complexity
interaction in triplet condition: $F(6,764)=4.07$, $p<.001$; Figure 5D). This
suggests that the effective complexity of problems was reduced by relying on
past action sequence statistics.

Figure 5: Model performance as a function of problem complexity in the test
phase of Experiment 1, in the duplet condition (A-B) and in the triplet
condition (C-D). Error bars indicate the $95\%$ CI. Color coding is identical
to that of Figure 4.

### Solving ambiguous problems

We trained the model variants on chunky silhouettes that enforced an action
chunk – that is, the chunked blocks had to be placed in fixed relative
positions and order for the correct solution. Here, we test the model variants
on so-called ’sequence-ambiguous silhouettes’ that have multiple solutions
allowing for placing the chunked blocks in different orders (Figure 6A and E).
Crucially, such sequence-breaking was unlikely under the sequence model,
making the sequence-breaking alternative a lower-value plan.

Indeed, solutions containing chunks were preferred by both MCTS-with-HABITS
and MCTS-with-HABITSopen-loop models (Figure 6C and G). We allowed for a
flexible node budget and measured the number of node evaluations performed
until the solution was reached (Figure 6D and H). In both the duplet and
triplet conditions, the number of evaluated nodes was significantly different
between the model variants, and the lowest in the case of the MCTS-with-HABITS
(model effect: $F(3,156)=25.86$, $p<.001$ and model effect: $F(3,156)=42.99$,
$p<.001$. To conclude, our model prefers to reuse past action sequences even
when other alternatives are viable and this confers it the computational
benefit of planning with smaller trees.

Figure 6: Action chunk reuse in sequence-ambiguous problems in the duplet
(A-C) and triplet conditions (D-F). Error bars indicate the $95\%$ CI. (A)(D)
Example silhouettes whose correct solutions allowed for flexible ordering of
the chunk elements. The digits indicate the possible serial positions of
blocks in different solutions. The blue digits indicate the serial positions
of the blocks in solutions that preserved the chunked order (same example
chunks as in Figure 3). (B)(E) Solutions containing chunks were preferred by
both model variants that were enabled to generate chunks. (C)(F) The MCTS-
with-HABITSfull model evaluated the least number of nodes until finding a
solution.

## Discussion

Discovering the building blocks of a problem is key to its efficient solution.
The behavior of animals suggest that they form action chunks that they
flexibly employ for solving novel tasks. We proposed that such actions chunks
can be usefully integrated into planning as well, effectively reducing the
depth of problems. We introduced the MCTS-with-HABITS, a model that finds and
leverages action sequence patterns at variable scales and integrates them into
the planning process in order to save computational costs.

We tested the model on a physical construction task that was constrained such
that certain blocks were placed in predictable positions relative to each
other. Our method used a flexible Bayesian sequence model, which can find
predictable action sequences of different sizes in past episodes, as induced
by the task constraints, which it can then reuse for the planning problem at
hand. We allowed the sequence model to influence the one-step search inherent
to MCTS, and also to propose multi-step node expansions to the planner. Both
mechanisms gave the model an edge in performance over the baseline model when
tested under resource-constraints. We showed that this was related to the
model being more economical with its node budget and searching deeper on
predictable paths via the multi-step expansions.

Various extensions to the method are possible. For instance, in our current
algorithms, the learned action sequences only influenced MCTS’ tree policy,
that is, the search among the nodes whose values were evaluated. One could
also use the sequence model to inform the rollout policy, potentially
improving the simulation step. Equally, it would be possible to store
information about the complexity of subsequent planning with nodes, and use
this to adjust the use of the sequence model depending on overall demands of
planning space or time.

In the current work, the construction task that our model was tested on was
more complex than typical psychology experiments but less complex than the
ones many real-life applications pose. Testing the model on more complex
problems and problems with partially transferable structure will be revealing
about the scalability of our approach. It would be particularly interesting to
look at stochastic problems, and those with intermediate rewards, which
exploit more of the power of MCTS.

Our model extracts and reuses action patterns in past behavior. Another class
of approaches focuses on past state visitation patterns. ? (?) proposed that
so-called bottleneck states can be identified at rare state transitions and
that they are useful subgoals that link subproblems. Several methods based on
the successor representation (?, ?, ?) use the eigendecomposition of state
succession patterns to carve the environment into subproblems. However, these
approaches require the agent to have explored the entire state space
extensively, a requirement that often cannot be met in machine learning and
robotics.

Applying the law of exercise to the realm of planning was inspired by the
cognitive science literature. How model-free (?, ?) or value-free (?, ?)
statistics interact with model-based reasoning has been a long-standing
interest in the field. Here we proposed a hybrid model that flexibly adjusts
the weighting of value-based and value-free assessments by relying on
predictably repeated chunks of behavior. Such sequence reuse yields a solution
for optimizing the direction and depth of search (?, ?, ?), reminiscent of
humans’ resource-rational strategies (?, ?). Future experiments should
investigate whether humans indeed use habits for planning, whether they
utilize such habits of mind adaptively in order to save planning costs, and,
if so, whether their planning is well described by the model presented here.

## Acknowledgments

NE and PD were funded by the Max Planck Society. PD was also funded by the
Alexander von Humboldt Foundation. The authors thank Tankred Saanum for
helpful conversations.

## References

  * Bapst et al.Bapst et al. Bapst, V., Sanchez-Gonzalez, A., Doersch, C., Stachenfeld, K., Kohli, P., Battaglia, P., Hamrick, J.  (2019).  Structured agents for physical construction.  In _International conference on machine learning_ (pp. 464–474).
  * Callaway et al.Callaway et al. Callaway, F., van Opheusden, B., Gul, S., Das, P., Krueger, P. M., Griffiths, T. L., Lieder, F.  (2022).  Rational use of cognitive resources in human planning.  _Nature Human Behaviour_ , _6_(8), 1112–1125.
  * Daw, Niv, DayanDaw et al. Daw, N. D., Niv, Y., Dayan, P.  (2005).  Uncertainty-based competition between prefrontal and dorsolateral striatal systems for behavioral control.  _Nature neuroscience_ , _8_(12), 1704–1711.
  * DayanDayan Dayan, P.  (2012).  How to set the switches on this thing.  _Current opinion in neurobiology_ , _22_(6), 1068–1074.
  * Dezfouli BalleineDezfouli Balleine Dezfouli, A., Balleine, B. W.  (2012).  Habits, action sequences and reinforcement learning.  _European Journal of Neuroscience_ , _35_(7), 1036–1051.
  * Dhawale, Wolff, Ko, ÖlveczkyDhawale et al. Dhawale, A. K., Wolff, S. B., Ko, R., Ölveczky, B. P.  (2021).  The basal ganglia control the detailed kinematics of learned motor skills.  _Nature Neuroscience_ , _24_(9), 1256–1269.
  * Éltető, Nemeth, Janacsek, DayanÉltető et al. Éltető, N., Nemeth, D., Janacsek, K., Dayan, P.  (2022).  Tracking human skill learning with a hierarchical bayesian sequence model.  _bioRxiv_.
  * GershmanGershman Gershman, S. J.  (2020).  Origin of perseveration in the trade-off between reward and complexity.  _Cognition_ , _204_ , 104394.
  * Huys et al.Huys et al. Huys, Q. J., Lally, N., Faulkner, P., Eshel, N., Seifritz, E., Gershman, S. J., … Roiser, J. P.  (2015).  Interplay of approximate planning strategies.  _Proceedings of the National Academy of Sciences_ , _112_(10), 3098–3103.
  * Kocsis SzepesváriKocsis Szepesvári Kocsis, L., Szepesvári, C.  (2006).  Bandit based monte-carlo planning.  In _European conference on machine learning_ (pp. 282–293).
  * Kuperwajs MaKuperwajs Ma Kuperwajs, I., Ma, W. J.  (2021).  Planning to plan: a bayesian model for optimizing the depth of decision tree search.  In _Proceedings of the annual meeting of the cognitive science society_ (Vol. 43).
  * Machado, Barreto, Precup, BowlingMachado et al. Machado, M. C., Barreto, A., Precup, D., Bowling, M.  (2021).  Temporal abstraction in reinforcement learning with the successor representation.  _arXiv preprint arXiv:2110.05740_.
  * Machado et al.Machado et al. Machado, M. C., Rosenbaum, C., Guo, X., Liu, M., Tesauro, G., Campbell, M.  (2017).  Eigenoption discovery through the deep successor representation.  _arXiv preprint arXiv:1710.11089_.
  * W. McCarthy, Kirsh, FanW. McCarthy et al. McCarthy, W., Kirsh, D., Fan, J.  (2020).  Learning to build physical structures better over time.  In _Cogsci._
  * W. P. McCarthy, Mattar, Kirsh, FanW. P. McCarthy et al. McCarthy, W. P., Mattar, M. G., Kirsh, D., Fan, J. E.  (2021).  Connecting perceptual and procedural abstractions in physical construction.  In _Proceedings of the annual meeting of the cognitive science society_ (Vol. 43).
  * McGovern BartoMcGovern Barto McGovern, A., Barto, A. G.  (2001).  Automatic discovery of subgoals in reinforcement learning using diverse density.
  * Miller, Shenhav, LudvigMiller et al. Miller, K. J., Shenhav, A., Ludvig, E. A.  (2019).  Habits without values.  _Psychological review_ , _126_(2), 292.
  * Riegler, Polani, SteuberRiegler et al. Riegler, B., Polani, D., Steuber, V.  (2021).  Embodiment and its influence on informational costs of decision density—atomic actions vs. scripted sequences.  _Frontiers in Robotics and AI_ , _8_ , 535158.
  * Sezener, Dezfouli, KeramatiSezener et al. Sezener, C. E., Dezfouli, A., Keramati, M.  (2019).  Optimizing the depth and the direction of prospective planning using information values.  _PLoS computational biology_ , _15_(3), e1006827.
  * TehTeh Teh, Y. W.  (2006).  A hierarchical bayesian language model based on pitman-yor processes.  In _Proceedings of the 21st international conference on computational linguistics and 44th annual meeting of the association for computational linguistics_ (pp. 985–992).
  * ThorndikeThorndike Thorndike, E. L.  (1898).  Animal intelligence: An experimental study of the associative processes in animals.  _The Psychological Review: Monograph Supplements_ , _2_(4), i.