# Exciton-phonon dynamics on complex networks: Comparison between a perturbative approach and exact calculations

Saad Yalouz, Vincent J.C. Pouthier, Cyril Falvo

► To cite this version:

Saad Yalouz, Vincent J.C. Pouthier, Cyril Falvo. Exciton-phonon dynamics on complex networks: Comparison between a perturbative approach and exact calculations. Physical Review E, 2017, 96 (2), pp.022304. 10.1103/PhysRevE.96.022304. hal-01716391

HAL Id: hal-01716391

https://hal.science/hal-01716391

Submitted on 16 Mar 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Exciton-phonon dynamics on complex networks: Comparison between a perturbative approach and exact calculations

Saad Yalouz and Vincent Pouthier*

Institut UTINAM, Université Bourgogne Franche-Comté, CNRS UMR 6213, 25030 Besançon Cedex, France

Cyril Falvo

Institut des Sciences Moléculaires d'Orsay (ISMO), CNRS, Université Paris-Sud, Université Paris-Saclay, 91405 Orsay, France

(Received 29 May 2017; published 4 August 2017)

A method combining perturbation theory with a simplifying ansatz is used to describe the exciton-phonon dynamics in complex networks. This method, called PT*, is compared to exact calculations based on the numerical diagonalization of the exciton-phonon Hamiltonian for eight small-sized networks. It is shown that the accuracy of PT* depends on the nature of the network, and three different situations were identified. For most graphs, PT* yields a very accurate description of the dynamics. By contrast, for the Wheel graph and the Apollonian network, PT* reproduces the dynamics only when the exciton occupies a specific initial state. Finally, for the complete graph, PT* breaks down. These different behaviors originate in the interplay between the degenerate nature of the excitonic energy spectrum and the strength of the exciton-phonon interaction so that a criterion is established to determine whether or not PT* is relevant. When it succeeds, our study shows the undeniable advantage of PT* in that it allows us to perform very fast simulations when compared to exact calculations that are restricted to small-sized networks.

DOI: 10.1103/PhysRevE.96.022304

## I. INTRODUCTION

In molecular crystals and large molecules, the study of exciton dynamics has a long history [1] that can be traced back to the seminal works of Agranovich [2], Davydov [3], and Silbey [4,5], to name but a few. Understanding how excitons carry energy is a key step in explaining many phenomena in physics, chemistry, and biology [6]. Nowadays, this knowledge is exploited to design technological devices with specific properties. Among the many examples are artificial light harvesting systems [7], organic light-emitting diodes [8], J-aggregates with specific optical properties [9], and exciton transistors [10].

Recently, it has been suggested that exciton propagation in molecular networks may be used to perform scalable quantum computing [11–13]. For instance, the exciton is a good candidate to promote high-fidelity quantum-state transfer (QST) at the nanoscale, a fundamental task to ensure an ideal communication within and between computers [14]. Depending on how the information is encoded, different strategies have been elaborated upon that involve, for example, electronic excitons in an array of quantum dots [15,16], spin excitations in Heisenberg ferromagnets [17–19], phonons in low-dimensional crystals [20–22], and vibrons in molecular nanowires [23]. Similarly, exciton delocalization in a hyperbranched molecule defines a physical realization of a continuous-time quantum walk (CTQW) [24]. CTQWs are a very active research topic due to their potential use in developing quantum algorithms [25,26]. Their powerlessness was demonstrated to solve various problems, such as the element distinctness problem [27], the hitting problem [28], and the quantum search problem [29,30]. Consequently, CTQWs have been characterized in a great variety of networks,

including dendrimers [31], trees [32,33], Apollonian networks [34,35], fractals [36,37], small-world networks [38], and star graphs [39–43].

In that context, exciton-mediated QST or CTQWs require that the exciton keeps its coherent wavelike nature over a sufficiently long time. Unfortunately, in a realistic system, the exciton does not propagate freely. It interacts with the remaining degrees of freedom of the medium that usually form a phonon bath responsible for quantum decoherence [44]. In large molecules, the phonons act as a reservoir that is insensitive to the exciton [45], and the Born-Markov approximation is legitimate. The dynamics is well-described by a generalized master equation (GME) that accounts for the irreversible decay of the coherence. By contrast, in finite-size lattices, the phonons no longer form a reservoir [46]. The Born-Markov approximation fails to capture the exciton dynamics, and the GME approach breaks down [47]. To overcome this problem, we recently introduced a method based on the operatorial formulation of perturbation theory (PT) [48,49]. Within this method, the exciton and the phonons are treated on an equal footing. The dynamics is governed by an effective Hamiltonian that takes exciton-phonon entanglement into account: the exciton gets dressed by a virtual phonon cloud, and the phonons are dressed by virtual excitonic transitions. Quantum decoherence occurs mainly because dressed phonons and free phonons evolve differently.

Until now, PT has been successfully applied in finite-size chains [48,49]. It provides a powerful tool for describing the spectral properties of the exciton-phonon Hamiltonian, and it is particularly suitable for characterizing the exciton dynamics. However, the powerlessness of the method is intimately related to the nondegenerate nature of the exciton spectra in linear chains. In complex networks, a different situation arises since the corresponding spectra usually exhibit more or less highly degenerate energy levels. In that case, the method requires additional numerical diagonalizations that may drastically

<footnote>*vincent.pouthier@univ-fcomte.fr</footnote>
increase the numerical cost of the procedure. To overcome this difficulty, we have proposed a simplifying ansatz [50,51]. From a detailed study of the star graph, it has been observed that when the phonons are accompanied by an exciton that occupies a degenerate eigenstate, their dynamics are basically described by the unperturbed phonon Hamiltonian. In other words, the phonon dressing effect is negligible when the exciton lies in a degenerate subspace. Although such an approximation successfully captures the dynamics in a star graph, the fundamental question arises as to whether it can be generalized to other graph families.

In that context, the aim of the present paper is to check the relevance of PT when it is combined with the previous ansatz. To proceed, this approach will be compared with exact calculations based on the numerical diagonalization of the exciton-phonon Hamiltonian. We shall consider the Holstein model [52], which refers to an exciton coupled with optical phonons, i.e., a model currently used to describe excitons in molecular crystals [53]. Both approximate and exact simulations will be carried out to describe the time evolution of the excitonic coherences in various networks. These coherences are the elements of the exciton reduced density matrix that measure the ability of the exciton to develop or to maintain superpositions involving the vacuum and one-exciton states. They thus correspond to key quantities whose knowledge is required for defining QST and CTQWs. Note that in this work, which is clearly more formal than applied, we focus our attention on the standard structures usually considered in quantum graph theory (see Fig. 1). This work can thus be viewed as a necessary step to validate our approach. The next step will be to consider more realistic structures such as dendrimers or light-harvesting complexes.

The paper is organized as follows. In Sec. II, the graphs we consider are described, and the exciton-phonon Hamiltonian is defined. Then, the perturbation theory is summarized and the excitonic coherences are defined. In Sec. III, a numerical analysis is performed to describe the time evolution of the excitonic coherence using both perturbation theory and exact calculations. Finally, the results are discussed in Sec. IV.

## II. THEORETICAL BACKGROUND

### A. Exciton-phonon Hamiltonian

The various networks we consider correspond to the graphs displayed in Fig. 1. Each graph $G = \{S, \mathcal{L}\}$ consists of a set of $N$ sites $S$, and a set of links $\mathcal{L}$ that connect sites. The sites are labeled by the index $\ell = 1, \dots, N$, and the notation $(\ell, \ell')$ stands for the link between two sites $\ell$ and $\ell'$. Within these notations, the graph is completely described by the so-called adjacency matrix $A_{\ell, \ell'}$. It reduces to 1 if $\ell$ and $\ell'$ correspond to connected sites, otherwise it vanishes:

$$A_{\ell\ell'} = \begin{cases} 1 & \text{if } (\ell, \ell') \in \mathcal{L}, \\ 0 & \text{otherwise.} \end{cases} \qquad (1)$$

In that context, the exciton dynamics on a graph is modeled as follows. Each site $\ell$ is occupied by a molecular subunit whose internal dynamics is described by a two-level system with Bohr frequency $\omega_0$. Let $|\ell\rangle$ denote the state in which the $\ell$th two-level system occupies its first excited state, the other two-level systems remaining in their ground state. Similarly,

<figure><img src="image_2.png" /><figcaption>FIG. 1. Representation of the networks considered in this paper. (a) The star graph, (b) the wheel graph, (c) the H graph, (d) the hat graph, (e) the fork graph, (f) the Apollonian network, (g) the complete graph, and (h) the random complete graph.</figcaption></figure>

let $|\emptyset\rangle$ stand for the vacuum state in which all the two-level systems occupy their ground state. The exciton Hamiltonian that governs both the zero- and the one-exciton dynamics is thus defined as (within the convention $\hbar = 1$)

$$H_A = \sum_{\ell=1}^{N} \omega_0 |\ell\rangle\langle\ell| + \sum_{\ell=1}^{N} \sum_{\ell'=1}^{N} \Phi_{\ell\ell'} |\ell\rangle\langle\ell'|, \qquad (2)$$

where $\Phi_{\ell\ell'} = \Phi A_{\ell\ell'}$ is the exciton hopping matrix. It reduces to a constant $\Phi$ when $\ell$ and $\ell'$ correspond to connected sites, otherwise it vanishes. The Hamiltonian $H_A$ acts in the Hilbert space $\mathcal{E}_A$ whose dimension is equal to $N+1$. Its diagonalization yields the unperturbed one-exciton eigenstates $|\chi_k\rangle$, with $k = 1, \dots, N$, the corresponding eigenenergies being denoted $\epsilon_k$. Within this notation, the exciton Hamiltonian is finally written as

$$H_A = \sum_{k=1}^{N} \epsilon_k |\chi_k\rangle\langle\chi_k|. \qquad (3)$$

According to the Holstein model, the exciton interacts with $N$ independent oscillators localized on each site of the graph and with frequency $\Omega_0$. These oscillators form a set of optical phonons whose Hamiltonian $H_B$ acts in the Hilbert space $\mathcal{E}_B$. In terms of the standard phonon operators $a_\ell^\dagger$ and $a_\ell$, $H_B$ is
written as

$$H_B = \sum_{\ell=1}^{N} \Omega_0 a_{\ell}^{\dagger} a_{\ell}. \qquad (4)$$

Within the so-called deformation potential model [6], the exciton-phonon interaction results from the modulation of the Bohr frequency of each two-level system induced by the local oscillators. The corresponding coupling Hamiltonian is written as

$$V = \sum_{\ell=1}^{N} M^{(\ell)}(a_{\ell}^{\dagger} + a_{\ell}), \qquad (5)$$

where $M^{(\ell)} = \Delta_0 |\ell\rangle\langle\ell|$ is expressed in terms of the exciton-phonon coupling strength $\Delta_0$. Acting in $\mathcal{E}_A$ only, $M^{(\ell)}$ is diagonal in the local basis $\{|\ell\rangle\}$, the $\ell$th element representing the influence of the $\ell$th optical phonon on the $\ell$th two-level Bohr frequency. By contrast, in the exciton eigenbasis $\{|\chi_k\rangle\}$, $M^{(\ell)}$ is no longer diagonal. Its elements are defined as

$$M_{kk'}^{(\ell)} = \Delta_0 \chi_{k\ell}^* \chi_{k'\ell}, \qquad (6)$$

where $\chi_{k\ell} = \langle\ell|\chi_k\rangle$ is the exciton wave function.

The exciton-phonon system is governed by the Hamiltonian $H = H_0 + V$, where $H_0 = H_A + H_B$ is the unperturbed Hamiltonian. The corresponding Hilbert space $\mathcal{E} = \mathcal{E}_A \otimes \mathcal{E}_B$ is partitioned into two independent subspaces: $\mathcal{E} = \mathcal{E}_0 \oplus \mathcal{E}_1$. In the zero-exciton subspace $\mathcal{E}_0$, $V = 0$. The exciton-phonon eigenstates are the unperturbed states involving the tensor products between the vacuum $|\emptyset\rangle$ and the phonon number states $\{|n_\ell\rangle\} = |n_1, \dots, n_N\rangle$. They describe $n_\ell$ free phonons with energy $n_\ell \Omega_0$ localized on each site $\ell$. In the one-exciton subspace $\mathcal{E}_1$, $V$ turns on. The unperturbed states $|\chi_k\rangle \otimes \{|n_\ell\rangle\}$ are no longer eigenstates of $H$. They refer to free phonons accompanied by an exciton in state $|\chi_k\rangle$ whereas the exact eigenstates correspond to entangled exciton-phonon states.

In the following, we shall restrict our attention to the so-called nonadiabatic weak-coupling limit: $\Phi \ll \Omega_0$ (nonadiabatic limit) and $\Delta_0 \le \Omega_0$ (weak-coupling limit). Note that such a situation corresponds, for instance, to high-frequency vibrational excitons (vibrons) propagating in a molecular lattice [54]. In that case, there is no resonance between the coupled unperturbed states, and second-order PT can be used to treat the coupling $V$.

## B. Quasidegenerate perturbation theory

As discussed in great detail in Refs. [50,51], the coupling $V$ is partially removed using quasidegenerate PT [55]. To proceed, we introduce a transformation $U = \exp(S)$ that generates a new point of view in which the effective Hamiltonian $\hat{H} = U H U^\dagger$ defines a phonon-conserving Hamiltonian that is block-diagonal in the unperturbed basis. Within PT, $S$ is expanded as a Taylor series with respect to $V$ so that the diagonalization is achieved at a given order. In the original point of view, the corresponding system eigenstates involve the transformation $U^\dagger$ so that they no longer describe independent excitations but refer to entangled exciton-phonon states.

Up to second order, the transformed Hamiltonian is expressed as

$$\hat{H} = H_A + \delta H_A + \sum_{\ell=1}^{N} \sum_{\ell'=1}^{N} [\Omega_0 \delta_{\ell\ell'} + \Lambda_{\ell\ell'}] a_{\ell}^{\dagger} a_{\ell'}, \qquad (7)$$

where $\delta H_A$ and $\Lambda_{\ell\ell'}$ are operators in $\mathcal{E}_A$ whose matrix elements are defined as (in the unperturbed basis $\{|\chi_k\rangle\}$)

$$\delta H_{Ak_1k_2} = \frac{1}{2} \sum_{\ell=0}^{N} \sum_{k=0}^{N} \frac{M_{k_1k}^{(\ell)} M_{kk_2}^{(\ell)}}{\epsilon_{k_1} - \epsilon_k - \Omega_0} + \frac{M_{k_1k}^{(\ell)} M_{kk_2}^{(\ell)}}{\epsilon_{k_2} - \epsilon_k - \Omega_0}, \\ \Lambda_{\ell\ell'k_1k_2} = \frac{1}{2} \sum_{k=0}^{N} \frac{M_{k_1k}^{(\ell)} M_{kk_2}^{(\ell')}}{\epsilon_{k_1} - \epsilon_k + \Omega_0} + \frac{M_{k_1k}^{(\ell')} M_{kk_2}^{(\ell')}}{\epsilon_{k_2} - \epsilon_k - \Omega_0} \\ + \frac{1}{2} \sum_{k=0}^{N} \frac{M_{k_1k}^{(\ell')} M_{kk_2}^{(\ell)}}{\epsilon_{k_1} - \epsilon_k - \Omega_0} + \frac{M_{k_1k}^{(\ell')} M_{kk_2}^{(\ell')}}{\epsilon_{k_2} - \epsilon_k + \Omega_0}. \qquad (8)$$

In Eq. (7), $\delta H_A$ is the correction of the exciton Hamiltonian that results from the coupling with the phonons. It originates from the spontaneous emission of a phonon during which the exciton is scattered from $|\chi_{k_1}\rangle$ to $|\chi_k\rangle$. Because the energy is not conserved during this scattering process, the emitted phonon is immediately reabsorbed and the exciton is again scattered from $|\chi_k\rangle$ to $|\chi_{k_2}\rangle$. Therefore, the exciton no longer propagate freely and it is dressed by a virtual phonon cloud. The dressed exciton forms a small polaron whose dynamics is governed by the effective Hamiltonian $\hat{H}_A = H_A + \delta H_A$. The diagonalization of $\hat{H}_A$ yields the polaron eigenstates $\{|\hat{\chi}_\mu\rangle\}$, with $\mu = 1, \dots, N$, and the corresponding eigenenergies $\hat{\epsilon}_\mu$, so that

$$\hat{H}_A = \sum_{\mu=1}^{N} \hat{\epsilon}_{\mu} |\hat{\chi}_{\mu}\rangle\langle\hat{\chi}_{\mu}|. \qquad (9)$$

In Eq. (7), $\Lambda_{\ell\ell'}$ is the phonon hopping constant matrix that accounts for the correction of the phonon Hamiltonian. It has two origins. First, a phonon can be absorbed on a particular site $\ell$ giving rise to an excitonic transition. Because this transition does not conserve energy, the phonon is immediately reemitted but on another site $\ell'$. Second, a phonon localized on a site $\ell$ can favor the stimulated emission of a second phonon during which the exciton realizes a transition. But, as discussed previously, the emitted phonon is immediately reabsorbed on a second site $\ell'$. Both mechanisms are virtual processes indicating that the phonons are dressed by virtual transitions realized by the exciton.

Strictly speaking, the phonon hopping constant matrix is an operator acting in the exciton subspace $\mathcal{E}_A$. It thus provides additional modifications of the polaron eigenstates, making it even more difficult to diagonalize the matrix $\hat{H}$. To overcome these difficulties, we follow our study of the star graph and propose two simplifying approximations [50,51]. First, in the weak-coupling limit, the modifications induced by the matrices $\Lambda_{\ell\ell'}$ remain small when compared with those induced by $\delta H_A$. Therefore, up to second order in $V$, the phonon hopping constant matrices will be treated to first order in the polaronic eigenbasis $\{|\hat{\chi}_\mu\rangle\}$. Second, two situations arise depending on whether the polaron eigenstate $|\hat{\chi}_\mu\rangle$ is degenerate or not. If $|\hat{\chi}_\mu\rangle$ is associated with a nondegenerate energy level, the
influence of the phonon hopping constant matrices to lowest order is encoded in the diagonal matrix elements $\langle \hat{\chi}_{\mu} | \Lambda_{\ell\ell'} | \hat{\chi}_{\mu} \rangle$. By contrast, if $| \hat{\chi}_{\mu} \rangle$ is associated with a degenerate energy level, the first-order treatment requires an additional numerical diagonalization in the subspace involving the tensor products between the corresponding polaronic degenerate subspace and the phonon subspace. Such a procedure may present a quite important numerical cost preventing the use of PT in networks of large size. To overcome this difficulty, we disregard such effects and propose an ansatz that consists in neglecting the influence of the phonon hopping constant matrices on the degenerate polaronic eigenstates. The method combining PT with this ansatz will be called PT*. Within these approximations, the effective exciton-phonon Hamiltonian is finally rewritten as

$$ \hat{H} \approx \sum_{\mu=1}^{N} \hat{\epsilon}_{\mu} |\hat{\chi}_{\mu}\rangle\langle\hat{\chi}_{\mu}| + \hat{H}_{B}^{(\mu)} \otimes |\hat{\chi}_{\mu}\rangle\langle\hat{\chi}_{\mu}|, \quad (10) $$

where $\hat{H}_{B}^{(\mu)}$ defines the Hamiltonian that governs the phonon dynamics when the exciton lies in the polaronic state $|\hat{\chi}_{\mu}\rangle$. Therefore, when the exciton occupies a degenerate eigenstate $|\hat{\chi}_{\mu}\rangle$, the phonon dynamics is described by the unperturbed Hamiltonian so that $\hat{H}_{B}^{(\mu)} = H_B$.

This is no longer the case when the exciton lies in a non-degenerate state. In that case, by setting $\Lambda_{\ell\ell'}^{(\mu)} = \langle \hat{\chi}_{\mu} | \Lambda_{\ell\ell'} | \hat{\chi}_{\mu} \rangle$, the phonon Hamiltonian is written as

$$ \hat{H}_{B}^{(\mu)} = \sum_{\ell=1}^{N} \sum_{\ell'=1}^{N} [\Omega_0 \delta_{\ell\ell'} + \Lambda_{\ell\ell'}^{(\mu)}] a_{\ell}^{\dagger} a_{\ell'}. \quad (11) $$

According to Eq. (11), the effect of the exciton is twofold. First, it favors a shift $\Lambda_{\ell\ell'}^{(\mu)}$ of the frequency of each local oscillator $\ell = 1, \dots, N$. Then, it yields a coupling $\Lambda_{\ell\ell'}^{(\mu)}$ between distinct local oscillators so that the dressed phonons can delocalize along the network. Consequently, the vibrational normal modes that define the dressed phonons are different from the localized normal modes associated with free phonons. These dressed normal modes are obtained from the diagonalization of the phonon hopping constant matrix $\Lambda^{(\mu)}$. Such a procedure allows us to define $N$ eigenvalues $\delta\Omega_q^{(\mu)}$ and $N$ eigenvectors $\beta_q^{(\mu)}(\ell)$ labeled by the index $q = 1, \dots, N$. The index $q$ refers to a particular phonon mode with energy $\Omega_q^{(\mu)} = \Omega_0 + \delta\Omega_q^{(\mu)}$, the eigenvalues of the phonon hopping constant matrix defining the phonon frequency shifts. The dynamics of each mode is described by the canonical creation $a_q^{(\mu)\dagger}$ and annihilation $a_q^{(\mu)}$ operators, defined by

$$ \begin{aligned} a_q^{(\mu)} &= \sum_{\ell=1}^{N} \beta_q^{(\mu)*}(\ell) a_{\ell}, \\ a_q^{(\mu)\dagger} &= \sum_{\ell=1}^{N} \beta_q^{(\mu)}(\ell) a_{\ell}^{\dagger}. \end{aligned} \quad (12) $$

Within the normal mode decomposition, the phonon Hamiltonian is finally rewritten in standard form as

$$ \hat{H}_{B}^{(\mu)} = \sum_{q} \Omega_{q}^{(\mu)} a_{q}^{(\mu)\dagger} a_{q}^{(\mu)}. \quad (13) $$

### C. Excitonic coherences

To test the accuracy of PT*, we propose to study the time evolution of excitonic coherences. According to the standard theory of open quantum systems, a complete understanding of the exciton dynamics is encoded in the reduced density matrix $\sigma(t) = \text{Tr}_B[\exp(-iHt)\rho(0)\exp(iHt)]$, where $\text{Tr}_B$ is a partial trace over the phonon degrees of freedom and where $\rho(0)$ is the initial exciton-phonon density matrix [44]. In that case, the “coherences” are the off-diagonal elements $\sigma_{\ell\ell'}(t)$ that yield information about the ability of the $\ell$th two-level system to develop or to maintain a superposition between its ground state and its first excited state at a time $t$.

Without any perturbation, the network is supposed to be in thermal equilibrium at temperature $T$. We consider situations for which $\omega_0 \gg k_B T$ so that the exciton cannot be thermally excited. By contrast, the phonons form a bath described by the Boltzmann distribution $\rho_B = \exp(-\beta H_B)/Z_B$, $Z_B$ being the phonon partition function and $\beta = 1/k_B T$. To study the coherences, the network is brought in a configuration out of equilibrium in which the exciton is prepared in a state $|\psi_A\rangle = c_0|\emptyset\rangle + c_1|\ell_0\rangle$, with $|c_0|^2 + |c_1|^2 = 1$. This step is supposed to be rather fast when compared with the typical time evolution of the phonons so that the initial density matrix becomes $\rho(0) = |\psi_A\rangle\langle\psi_A| \otimes \rho_B$.

In that context, the coherences are expressed as $\sigma_{\ell\ell'}(t) = G_{\ell\ell_0}(t)\sigma_{\ell_0\ell}(0)$ with

$$ G_{\ell\ell_0}(t) = \langle \ell | \text{Tr}_B [\rho_B e^{iH_B t} e^{-iHt}] | \ell_0 \rangle. \quad (14) $$

The effective exciton propagator $G_{\ell\ell_0}(t)$ generalizes the concept of transition amplitude. It yields the probability amplitude to observe the exciton in $|\ell\rangle$ at time $t$ given that it was in $|\ell_0\rangle$ at $t=0$. Its effective nature results from the fact that the exciton interacts with the phonons during its propagation.

The effective propagator is the central object of the present study. Knowledge of this propagator provides key information on the way the phonon bath modifies the excitonic coherences either in the local basis or in the eigenbasis (by performing a change of basis). As detailed previously [48,49], it can be evaluated quite straightforwardly using PT*. To proceed, we first introduce the transformation $U$ and diagonalize $H$ in Eq. (14). Then, we define the effective phonon density matrix

$$ \rho_B^{(\mu)}(t) = \frac{1}{Z_B^{(\mu)}(t)} \exp \left[ -\beta H_B + i(H_B - \hat{H}_B^{(\mu)})t \right], \quad (15) $$

where

$$ Z_B^{(\mu)}(t) = \text{Tr}_B \exp \left[ -\beta H_B + i(H_B - \hat{H}_B^{(\mu)})t \right]. \quad (16) $$

Strictly speaking, $\rho_B^{(\mu)}(t)$ is not a density matrix since it yields complex values for the phonon population. However, it is isomorphic to the phonon Boltzmann distribution and yields averages equivalent to thermal averages with time-dependent temperatures. After simple algebraic manipulations, $G_{\ell\ell_0}(t)$ is rewritten as

$$ G_{\ell\ell_0}(t) = \sum_{\mu=1}^{N} \frac{Z_B^{(\mu)}(t)}{Z_B} e^{-i\hat{\epsilon}_\mu t} \langle\langle\ell|U_\mu^\dagger(t)|\hat{\chi}_\mu\rangle\langle\hat{\chi}_\mu|U_\mu(0)|\ell_0\rangle\rangle_{\mu,t}, \quad (17) $$
where $\langle \cdots \rangle_{\mu,t} = \text{Tr}_B[\rho_B^{(\mu)}(t) \cdots]$ and $U_\mu(t) = e^{i\hat{H}_B^{(\mu)t}} U e^{-i\hat{H}_B^{(\mu)t}}$. Expanding $U$ in a Taylor series with respect to $V$, one finally obtains the second-order expression for $G_{\ell\ell_0}(t)$. This expression is given in great detail in our previous works (see, for instance, Ref. [51]).

To check the relevance of PT\*, exact simulations are carried out for describing the time evolution of the coherences. To proceed, we have developed a code in which we explicitly build the Hamiltonian $H = H_A + H_B + V$ using a truncated (finite) basis for the phonon number states $|n_1, \dots, n_N\rangle$. The numerical process we employ can be summarized as follows. First, we set $N_P$ the maximum number of phonons on the network by fixing the value $N_P = \sum_\ell n_\ell$. In practice, depending on the temperature, we judiciously choose $N_P$ in order to correctly reproduce the phonon bath statistics (the error on the average phonon number is less than 0.1%). Thus, we ensure a good convergence of all our numerical calculations. Second, using boson symmetry properties, we numerically build the reduced set of phonon number states $|n_1, \dots, n_N\rangle$ to span the truncated Hilbert space whose dimension reduces to $\dim(\mathcal{E}_B) = (N_P + N)!/(N_P!N!)$. Finally, within this reduced basis, we explicitly calculate each contribution of the three operators $H_A, H_B$, and $V$ in the global exciton-phonon basis states formed by the tensor products between the one-exciton states $|\ell\rangle$ and the phonon number states. Without going into too much detail, the numerical process we use consists in calculating the exciton-phonon state resulting from the action of $H$ over an initial state. The real difficulty is to correctly reference the outgoing state as a function of a phonon number state and a local excitonic state. Thus, step by step, we build all the matrix elements of the full Hamiltonian $H$. In this way, we obtain a very large matrix $H$ that can be diagonalized numerically. In doing so, we get access to all the required information to efficiently simulate the time evolution of the effective propagator.

### III. NUMERICAL RESULTS

In this section, both PT\* and exact calculations are performed in order to describe the time evolution of the exciton effective propagator $|G_{\ell\ell_0}(t)|$. The following reduced parameters are used: the exciton-phonon coupling strength and the phonon frequency are fixed to $\Delta_0/\Phi = 5$ and $\Omega_0/\Phi = 100$, respectively, and the temperature is expressed in terms of the Einstein temperature $T_E = \Omega_0/k_B$ through the introduction of the parameter $\alpha = T/T_E$ with either $\alpha = 0.5$ or 1. Note that, as mentioned in Sec. II A, this set of parameters corresponds to the weak-coupling nonadiabatic limit, a situation encountered, for instance, when studying vibrational excitons in molecular networks [54]. In our numerical analysis, we consider small-sized networks with $N=6$. Exact calculations are carried out with a maximum phonon number equal to $N_P=9$ so that the size of the truncated basis reduces to 30 030. Note that for the Apollonian network, we have $N=7$. In that case, the maximum phonon number is fixed to $N_P=8$, resulting in a truncated basis whose size reaches 45 045. In the following figures, red curves (gray curves in the printed black-white versions) refer to exact calculations whereas black curves correspond to PT\* calculations. Note that the short-time

<figure><img src="image_3.png" /><figcaption>FIG. 2. Time evolution of the effective propagator on the star graph ($N=6$ and $\alpha=0.5$). (a) Short-time and (b) long-time behavior of the element $|G_{12}(t)|$. (c) Short-time and (d) long-time behavior of the element $|G_{62}(t)|$. In (a) and (c), red dots correspond to exact calculation whereas black lines refer to PT\* calculations.</figcaption></figure>

dynamics is illustrated for the star graph only (see Fig. 2), a similar behavior having been observed for the other graphs.

For the star graph, the time evolution of the propagator $|G_{12}(t)|$ is illustrated in Figs. 2(a) and 2(b). It characterizes the ability of the exciton to reach the core site $\ell=1$ at time $t$ given that it was initially created on the peripheral site $\ell=2$ [see Fig. 1(a)]. Initially equal to zero, the propagator rapidly
reaches a first maximum at time $t = 0.7\Phi^{-1}$ and its amplitude is equal to 0.44 [see Fig. 2(a)]. This peak characterizes a direct exciton transfer between the periphery and the core of the star. Then, in the short-time limit ($t < 50\Phi^{-1}$), $|G_{12}(t)|$ exhibits high-frequency oscillations whose amplitude remains smaller than 0.44. The period of these oscillations is approximately equal to $1.4\Phi^{-1}$, i.e., very close to the theoretical period $\pi\Phi^{-1}/\sqrt{N-1} \simeq 1.405\Phi^{-1}$ that describes the dynamics of a free exciton. As time increases, the maximum value of these oscillations is modulated by a slowly varying envelope [see Fig. 2(b)]. It turns out that the propagator does not vanish in the long-time limit, and different features are observed. First, the envelope exhibits incomplete quantum revivals that recur almost periodically at times $T_r, 2T_r, 3T_r, ...$ with the almost period $T_r = 880\Phi^{-1}$. Note that the amplitude of the first revival is equal to 0.32. Second, the envelope decreases until it reaches a minimum value at time $t = 2400\Phi^{-1}$ and its amplitude reaches approximately 0.12. Therefore, as time increases, the amplitude of the envelope increases again until an almost complete revival takes place at a superrevival time $T_{sr}$. In the present situation, one obtains $T_{sr} = 5600\Phi^{-1}$, the corresponding amplitude of the propagator being approximately equal to 0.44.

Figures 2(c) and 2(d) display the time evolution of the propagator $|G_{62}(t)|$ that characterizes the exciton transfer between two peripheral sites $l = 2$ and 6 [see Fig. 1(a)]. As discussed previously, the propagator initially equal to zero reaches a first maximum equal to 0.4 at time $t = 1.4\Phi^{-1}$ [see Fig. 2(c)]. Then, as time increases, the propagator exhibits high-frequency oscillations whose amplitude follows a slowly varying envelope. These oscillations vary around a nonvanishing value approximately equal to 0.2 and with a period equal to $2.8\Phi^{-1}$ [see Fig. 2(c)]. Note that the amplitude of the oscillations decreases as time increases and the propagator no longer goes through zero. Nevertheless, the slowly varying envelope supports incomplete revivals with revival time $T_r = 880\Phi^{-1}$ as well as an almost exact recurrence that arises at the superrevival time $T_{sr} = 5600\Phi^{-1}$.

As shown in Fig. 2, PT\* describes very accurately the time evolution of the effective propagator, even in the long-time limit. It perfectly reproduces the high-frequency oscillating behavior of the propagator as well as the occurrences of both the incomplete revivals and the almost exact recurrence. In addition, we have verified that this accuracy remains even at high temperature ($\alpha = 1$) and for all the matrix elements of the propagator.

For the wheel graph [see Fig. 1(b)], the time evolution of the propagator $|G_{12}(t)|$ is shown in Fig. 3(a). Initially equal to zero, $|G_{12}(t)|$ reaches a first maximum equal to 0.4 at time $t = 0.6\Phi^{-1}$. Then, it shows high-frequency oscillations that vary around a nonvanishing value approximately equal to 0.15. The period of these oscillations is equal to $1.2\Phi^{-1}$. As for the star graph, these oscillations are modulated by a slowly varying envelope that exhibits incomplete revivals. These revivals arise almost periodically, with the corresponding almost period being approximately equal to $T_r = 1000\Phi^{-1}$. Nevertheless, the amplitude of the envelope varies from one revival to another, and no exact recurrence is observed. In that case, the figure clearly shows that PT\* provides a very good estimate of the propagator. It yields a perfect description of the oscillating

<figure><img src="image_4.png" /><figcaption>FIG. 3. Time evolution of the effective propagator elements (a) $|G_{12}(t)|$ and (b) $|G_{62}(t)|$ on the wheel graph ($N = 6$ and $\alpha = 0.5$).</figcaption></figure>

behavior of the propagator as well as of the shape of the envelope. Note that a very slight discrepancy arises between the numerical calculations and PT\* calculations in the very long-time limit.

A different situation occurs when one considers the propagator $|G_{62}(t)|$ whose time evolution is displayed in Fig. 3(b). $|G_{62}(t)|$ increases from its vanishing initial value to reach a maximum equal to 0.63 at time $t = 15.5\Phi^{-1}$. Then, it exhibits high-frequency oscillations whose amplitude decreases until the time becomes equal to $800\Phi^{-1}$. As time increases, the oscillations vary around a nonvanishing value approximately equal to 0.2 for the exact simulation and equal to 0.35 for the PT\* approach. In the two cases, the period of high oscillations is the same, namely $\sim 1.2\Phi^{-1}$. The amplitude of these oscillations is modulated by a slowly varying envelope that exhibits small-amplitude revivals. These revivals recur almost periodically with an almost period approximately equal to $T_r = 1000\Phi^{-1}$. For this coherence, PT\* yields a quite good estimate of the propagator in the short-time limit, but it fails to reproduce the time evolution of the propagator at long times. Although the shape of the envelope is more or less well captured, PT\* clearly overestimates the value of the propagator. We have verified that such a discrepancy affects all the propagator elements that involve two peripheral sites of the wheel graph. Moreover, it is enhanced as the temperature increases.

As illustrated in the following three figures, PT\* provides a very good estimate of the time evolution of the effective propagator for both the H graph (Fig. 4), the hat graph (Fig. 5), and the fork graph (Fig. 6). It describes very accurately the high-frequency oscillations of the propagator, and it perfectly accounts for the occurrence of revivals. This accuracy remains
<figure><img src="image_5.png" /><figcaption>FIG. 4. Time evolution of the effective propagator elements (a) $|G_{12}(t)|$ and (b) $|G_{65}(t)|$ on the H graph ($N = 6$ and $\alpha = 0.5$).</figcaption></figure>

even at high temperature ($\alpha = 1$) and for all the matrix elements of the propagator.

For the H graph [see Fig. 1(c)], we consider two different elements of the propagator. The element $|G_{12}(t)|$ describes the exciton transfer between the peripheral sites $\ell = 2$ with degree 1 and the site $\ell = 1$ with degree 3 (the degree of a

<figure><img src="image_6.png" /><figcaption>FIG. 5. Time evolution of the effective propagator elements (a) $|G_{12}(t)|$ and (b) $|G_{62}(t)|$ on the hat graph ($N = 6$ and $\alpha = 0.5$).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Time evolution of the effective propagator elements (a) $|G_{16}(t)|$ and (b) $|G_{63}(t)|$ on the fork graph ($N = 6$ and $\alpha = 0.5$).</figcaption></figure>

site is the number of links that emanate out from that site). As illustrated in Fig. 4(a), $|G_{12}(t)|$ first increases from zero to reach a maximum amplitude of 0.58 at time $t = 0.95\Phi^{-1}$. Then, it exhibits high-frequency oscillations whose mean period is $1.5\Phi^{-1}$. These oscillations are modulated by a slowly varying envelope that first decreases. Then, the envelope increases again until an almost complete revival takes place at a superrevival time $T_{sr} = 7300\Phi^{-1}$. The time evolution $|G_{65}(t)|$ is displayed in Fig. 4(b). It characterizes the exciton transfer between two peripheral sites $\ell = 2$ and 6. After reaching its first maximum equal to 0.75 at time $t = 2.1\Phi^{-1}$, $|G_{65}(t)|$ shows high-frequency oscillations whose amplitude follows a slowly varying envelope. These oscillations vary around a nonvanishing value approximately equal to 0.5. Note that at time $T_{sr} = 7600\Phi^{-1}$, a very pronounced recurrence arises.

For the hat graph [see Fig. 1(d)], the element $|G_{12}(t)|$, which characterizes the transfer between the top site $\ell = 1$ and the base site $\ell = 2$, basically behaves as the element $|G_{12}(t)|$ of the H graph, except that no pronounced revivals are observed [see Fig. 5(a)]. For the element $|G_{62}(t)|$, a remarkable feature arises since the high-frequency oscillations of the propagator are modulated by a slowly varying envelope that scales as a sinelike function [Fig. 5(b)]. Therefore, revivals recur periodically according to an almost period $T_r = 2550\Phi^{-1}$.

As shown in Fig. 6, similar features are observed for the fork graph [see Fig. 1(e)]. First, the time evolution of $|G_{16}(t)|$, which refers to a transition between two peripheral nodes with degree unity, is quite similar to that of the element $|G_{12}(t)|$ of the hat graph [see Fig. 6(a)]. Then, for the element $|G_{63}(t)|$ [Fig. 6(b)], which describes a transition between the central node $\ell = 3$ with degree 4 and the peripheral node $\ell = 6$, one recovers the occurrence of a slowly varying envelope that scales as a sinelike function.
<figure><img src="image_8.png" /><figcaption>FIG. 7. Time evolution of the effective propagator elements (a) $|G_{72}(t)|$ and (b) $|G_{12}(t)|$ on the Apollonian network ($N = 7$ and $\alpha = 0.5$).</figcaption></figure>

For the Apollonian network [see Fig. 1(f)], a different situation arises since the accuracy of PT\* depends on the propagator we consider. First, Fig. 7(a) displays the evolution of the element $|G_{72}(t)|$, which describes the exciton transfer between the core site $\ell = 7$, whose degree is equal to 6, and the peripheral site $\ell = 2$ with degree 5. As time increases, the propagator, which is initially equal to zero, reaches a first maximum equal to 0.47 at time $t = 3\Phi^{-1}$. Then, it shows high-frequency oscillations that vary around a nonvanishing value approximately equal to 0.15. The mean period of these oscillations is $1.2\Phi^{-1}$. These oscillations are modulated by a slowly varying envelope that exhibits incomplete revivals. These revivals arise almost periodically, the corresponding almost period being approximately equal to $T_r = 1250\Phi^{-1}$. Nevertheless, the amplitude of the envelope varies from one revival to another, and no exact recurrence is observed. In that case, PT\* yields a very good estimate of the propagator, a situation encountered each time one considers a matrix element of the propagator that involves the central node $\ell = 7$.

A different situation occurs when one considers the propagator $|G_{12}(t)|$ whose time evolution is displayed in Fig. 7(b). This propagator refers to the excitonic transfer between the two peripheral sites $\ell = 2$ and 1 [see Fig. 1(f)]. From its initial value equal to zero, $|G_{12}(t)|$ increases to reach its maximum value equal to 0.65 at time $t = 25\Phi^{-1}$. As time increases, it exhibits high-frequency oscillations whose amplitude decreases until the time becomes equal to $980\Phi^{-1}$. Then, the oscillations vary around 0.2 for exact calculations and 0.3 for PT\* calculations. The corresponding period $\sim 1\Phi^{-1}$ is the same for both approaches. The amplitude of these oscillations is modulated by a slowly varying envelope that exhibits small-amplitude revivals. These revivals recur almost

<figure><img src="image_9.png" /><figcaption>FIG. 8. Time evolution of the effective propagator elements (a) $|G_{62}(t)|$ and (b) $|G_{66}(t)|$ on the complete graph ($N = 6$ and $\alpha = 0.5$).</figcaption></figure>

periodically with an almost period approximately equal to $T_r = 1250\Phi^{-1}$.

As for the wheel graph, PT\* allows us to capture the shape of the envelope, but it clearly overestimates the value of the propagator. Such a discrepancy affects all the propagator elements that involve the sites of the Apollonian network that differ from the central site. In addition, it is enhanced as the temperature increases, so that PT\* leads to completely erroneous results.

As shown in Fig. 1(g), the complete graph involves $N = 6$ equivalent sites whose degree reduces to the number of nodes. Therefore, the propagator is characterized by two different kinds of elements, i.e., those involving distinct sites (transfer amplitude) and those involving the same site (survival amplitude). In that context, the time evolution of the transfer amplitude $G_{12}(t)$ is shown in Fig. 8(a). Initially equal to zero, the propagator rapidly increases to reach a first peak value equal to 0.33 at time $t = 0.55\Phi^{-1}$. Then, it develops high-frequency oscillations whose mean period is approximately equal to $1.1\Phi^{-1}$ and whose amplitude is modulated by a slowly varying envelope. This envelope behaves as an almost periodic function that exhibits two kinds of revivals. Indeed, small-amplitude revivals take place at times $T_r, 2T_r, 3T_r, ...$ with $T_r = 1250\Phi^{-1}$, their amplitude being approximately equal to 0.2. By contrast, almost exact recurrences arise regularly according to the characteristic time $T_{sr} = 2500\Phi^{-1}$, their amplitude being approximately equal to 0.3. As shown in Fig. 8(a), PT\* partially captures the behavior of the propagator. More precisely, it perfectly accounts for its short-time evolution. It predicts the occurrence of high-frequency oscillations with a mean period $\sim 1.05\Phi^{-1}$ that basically corresponds to the period obtained from exact
<figure><img src="image_10.png" /><figcaption>FIG. 9. Time evolution of the effective propagator elements (a) $|G_{62}(t)|$ and (b) $|G_{66}(t)|$ on the inhomogeneous complete graph ($N = 6$ and $\alpha = 0.5$).</figcaption></figure>

calculations. However, PT\* fails to reproduce the shape of the envelope. It favors the occurrence of a periodic sinelike envelope whose period is equal to $1250\Phi^{-1}$. Consequently, PT\* accounts for the occurrence of the complete revivals, at least over the time scale shown in the figure. Conversely, it clearly overestimates the amplitude of the incomplete revivals.

The behavior of the survival amplitude is displayed in Fig. 8(b). Initially equal to unity, it exhibits high-frequency oscillations whose mean period is equal to $1.05\Phi^{-1}$. These oscillations are modulated by a sinelike envelope that shows an almost periodic behavior, the almost period being approximately equal to $T_{sr} = 2600\Phi^{-1}$. As a result, almost complete recurrence arises according to the characteristic time $T_{sr}$ for which the survival amplitude reaches a value close to unity. In that case, the figure reveals that PT\* completely fails to reproduce the time evolution of the survival amplitude. Although it predicts high-frequency oscillations that match the shape of the propagator in the short-time limit, it clearly overestimates the propagator. It yields a nonvanishing slowly varying envelope whose amplitude varies periodically around an average value equal to 0.83, the period being approximately equal to $1250\Phi^{-1}$. Therefore, PT\* predicts the occurrences of exact revivals that recur periodically.

Finally, let us consider the inhomogeneous complete graph in which the site energies have been chosen randomly [see Fig. 1(h)]. Each site energy corresponds to a random variable uniformly distributed over the energy range $[\omega_0 - 0.5\Phi, \omega_0 + 0.5\Phi]$. The behavior of the transfer amplitude $G_{12}(t)$ is displayed in Fig. 9(a) for a given set of site energies. In the short-time limit, the propagator rapidly increases from its vanishing initial value. It thus reaches a maximum value equal to $0.35$ at time $t = 0.55\Phi^{-1}$. Then it exhibits high-frequency

oscillations whose mean period is approximately equal to $1.1\Phi^{-1}$. These oscillations are modulated by a slowly varying envelope whose shape evolves as time increases. This envelope behaves as a damped sinelike function whose period is equal to $1250\Phi^{-1}$, indicating that revivals recur periodically. However, as time increases, the oscillating nature of the envelope tends to disappear. The time evolution of the survival amplitude is shown in Fig. 9(b). Initially equal to unity, it exhibits high-frequency oscillations whose period is $\sim 1\Phi^{-1}$. These oscillations vary around an average value approximately equal to 0.5. They are modulated by a slowly varying envelope that exhibits a sinelike behavior. As a result, almost complete recurrence arises according to the superrevival time $T_{sr} = 1250\Phi^{-1}$ and for which the survival amplitude reaches a value close to unity. As shown in Fig. 9, PT\* describes very accurately the time evolution of the effective propagator. This accuracy remains at high temperature ($\alpha = 1$) and for all the matrix elements of the propagator.

## IV. DISCUSSION

The numerical results show that the accuracy of the method PT\* depends on the nature of the graph. Three different situations have been identified. First, for most graphs, i.e., for the star graph, the H graph, the fork graph, the hat graph, and the complete random graph, PT\* provides a very accurate description of the excitonic dynamics. This method perfectly reproduces the oscillating behavior of the effective exciton propagator over short, intermediate, and long time scales. In addition, it accounts for the occurrence of both incomplete revivals and almost exact recurrences. For these graphs, the fundamental point is that PT\* works very well for all temperatures and regardless of the element of the propagator under study. Then, for the Wheel graph and the Apollonian network, the relevance of PT\* depends on the nature of the elements of the excitonic propagator, and two situations arise. These graphs are organized around a central site, and they exhibit a specific symmetry. We have shown that PT\* is particularly suitable for describing the time evolution of the matrix elements of the propagator that involve the central site. Conversely, for the other elements, PT\* fails to reproduce the dynamics. It overestimates the excitonic coherence, a feature that is enhanced as the temperature increases. Finally, for the complete graph, PT\* breaks down. The method provides satisfactory results over short time scales, only. For intermediate and long time scales, it overestimates the excitonic coherence, regardless of the element of the propagator that we consider, this effect being strongly exacerbated as the temperature increases.

In a general way, the relevance of PT\* is intimately related to the validity of our ansatz presented in Sec. II. According to this ansatz, we assume that the influence of the phonon hopping constant matrices $\Lambda_{\ell\ell'}$ on the degenerate polaronic eigenstates can be neglected. To check this hypothesis, we have realized a detailed study of the phonon hopping constant matrices in order to determine whether or not the residual exciton-phonon coupling induced by those matrices can be disregarded. To proceed, for degenerate polaronic eigenstates, we define a measure of the coupling strength as the spectral radius $\rho_D$ of the restriction of the matrix $\Lambda_{\ell\ell'}^{\mu\mu'} = \langle \hat{x}_\mu | \Lambda_{\ell\ell'} | \hat{x}_{\mu'} \rangle$ to the
<table><caption>TABLE I. Spectral radius in $\Phi$ units.</caption><thead><tr><th></th><th>ρ<sub>D</sub></th><th>ρ<sub>ND</sub></th></tr></thead><tbody><tr><td>Star graph</td><td>7.5 × 10<sup>-5</sup></td><td>6.7 × 10<sup>-3</sup></td></tr><tr><td>Wheel graph</td><td>4.5 × 10<sup>-3</sup></td><td>2.5 × 10<sup>-2</sup></td></tr><tr><td>H graph</td><td>2.5 × 10<sup>-7</sup></td><td>5.3 × 10<sup>-2</sup></td></tr><tr><td>Fork graph</td><td>2.5 × 10<sup>-5</sup></td><td>6.0 × 10<sup>-3</sup></td></tr><tr><td>Apollonian network</td><td>5.0 × 10<sup>-3</sup></td><td>5.3 × 10<sup>-3</sup></td></tr><tr><td>Complete graph</td><td>1.6 × 10<sup>-2</sup></td><td>5.0 × 10<sup>-3</sup></td></tr></tbody></table>

degenerate subspace under study. Note that the spectral radius of a square matrix is the largest element among the absolute values of the eigenvalues of that matrix [56]. In the same way, in the case of a nondegenerate state $|\hat{\chi}_{\mu}\rangle$, the coupling strength corresponds to the largest eigenvalues $|\delta\Omega_q^{(\mu)}|$ of the matrix $\Lambda_{\ell\ell'}^{\mu}$ (see Sec. II B). This quantity defines the matrix spectral radius $\rho_{ND}$. Within these notations, one thus expect that $PT^*$ is relevant provided that $\rho_D \ll \rho_{ND}$.

The pertinence of $PT^*$ can thus be interpreted as follows. First, for both the hat graph and the random complete graph, the excitonic spectrum does not show degeneracy. In that case, it is not necessary to invoke our ansatz since $PT^*$ corresponds to standard perturbation theory. It thus works very well provided that the exciton-phonon coupling remains sufficiently weak, as demonstrated previously [48,49]. Then, for the star graph, the H graph, and the fork graph, the excitonic spectrum exhibits degeneracy. However, two situations arise. If one considers the elements of the propagator that involve a site whose degree is greater than 1, as well as the site 1 for the fork graph, it behaves as if the excitonic dynamics was confined in the subspace entirely generated by the nondegenerate eigenstates. Therefore, once again, $PT^*$ is equivalent to PT so that it provides a perfect description of the excitonic dynamics. Conversely, if one considers the elements of the propagator that involve the other sites, the degenerate excitonic states participate in the dynamics. However, for these three graphs, the relation $\rho_D \ll \rho_{ND}$ is satisfied, as shown in Table I. Consequently, $PT^*$ perfectly reproduces the time evolution of the propagator.

A different situation appears for the Wheel graph and for the Apollonian network whose excitonic spectra exhibit degenerate energy levels. In both cases, the spectral radius $\rho_{ND}$ and $\rho_D$ are of the same order of magnitude so that our ansatz breaks down (see Table I). Nevertheless, two situations arise depending on the nature of the propagator under study. If one considers the elements of the propagator that involve the central site, the excitonic dynamics is confined in the subspace generated by the nondegenerate eigenstates. Therefore, as discussed previously, $PT^*$ is equivalent to PT so that it perfectly describes the dynamics. By contrast, if one considers the elements of the propagator that involve the other sites, the degenerate excitonic states participate in the dynamics. Since the ansatz is no longer valid, $PT^*$ fails to reproduce the temporal evolution of the propagator.

Finally, for the complete graph with $N$ sites, the excitonic spectrum shows an $(N-1)$-fold degenerate eigenenergy that plays a key role in the dynamical processes. However, the analysis of the phonon hopping constant matrices reveals that

$\rho_D$ is one order of magnitude larger than $\rho_{ND}$ (see Table I). In this particular case, regardless of the elements of the propagator under study, our ansatz is no longer valid so that $PT^*$ cannot be used to describe the excitonic properties. At this point, let us mention that $PT^*$ gives very accurate results when a weak disorder arises in the complete graph. The disorder breaks the symmetry, and the degeneracy is lifted so that $PT^*$ becomes equivalent to PT.

These different results clearly show that $PT^*$ is a powerful tool for studying the excitonic dynamics on a graph, as long as the exciton-phonon coupling remains weak. Although it cannot be generalized, it seems to be adapted to a large number of graphs. Moreover, an analysis of the spectral radius of the phonon hopping constant matrices in the degenerate excitonic subspaces allows us to determine whether or not the method is relevant.

To conclude, we must highlight the undeniable advantage $PT^*$ has over exact calculations. Indeed, the perturbative approach allows us to perform very fast calculations, even for graphs whose size is important. Conversely, because of the very important size of the Hamiltonian matrix $H$, exact calculations can only be carried out on small graphs. And even in this case, the computing times become very long when compared to those provided by $PT^*$ simulations. Let us take an example: we consider here a graph with $N=6$ sites, and we want to simulate $10^5$ times steps. On the one hand, the use of exact calculations (with a maximum phonon number $N_p=9$) requires a CPU time approximately equal to 35 h using an OpenMP INTEL MKL subroutine to diagonalize numerically the Hamiltonian matrix $H$ and distributing the computational load over 16 cores. On the other hand, without any optimization processes and with only one core of our laptop (3.6 GHz Intel Core i7 processor), the $PT^*$ simulation only required 12 s of computational time. Thus, the perturbative approach appears to be $10^4$ times faster than exact simulation in that context. Therefore, this method is a very efficient way to obtain good results in a very short time.

## V. CONCLUSION

In this paper, a method involving standard perturbation theory combined with a simplifying ansatz has been introduced to describe the dynamics of the exciton-phonon system in a complex molecular network. A detailed numerical analysis has been carried out to check the relevance of this method, which is called $PT^*$. $PT^*$ has been compared with exact calculations based on the numerical diagonalization of the exciton-phonon Hamiltonian for eight small-sized networks. Special attention has been paid to characterize the time evolution of the effective propagator. This propagator is a measure of the excitonic coherence, and it characterizes the propagation of excitonic superpositions involving the vacuum and one-exciton states.

It has been shown that the accuracy of $PT^*$ depends strongly on the nature of the network, so that three different situations have been identified. First, for most graphs (star graph, H graph, fork graph, hat graph, and complete random graph), $PT^*$ yields a very accurate description of the exciton dynamics. It perfectly reproduces the oscillating behavior of the effective propagator, and it accounts for the occurrence of
both incomplete and almost exact quantum recurrences. PT* works very well at all temperatures and for all the excitonic coherences we have considered. Then, for the Wheel graph and the Apollonian network, we have shown that the relevance of PT* depends on the nature of the element of the propagator. These graphs are organized around a central site so that they exhibit a specific symmetry. Therefore, PT* perfectly accounts for the time evolution of the matrix elements of the propagator that involve the central site. By contrast, for the other elements, PT* fails to reproduce the dynamics, a feature that is enhanced as the temperature increases. Finally, for the complete graph, the PT* method breaks down.

We have shown that these different behaviors originate in the interplay between the degenerate nature of the excitonic energy spectrum and the strength of the exciton-phonon interaction. Therefore, a criterion has been established to

determine whether or not PT* is relevant. Consequently, although it cannot be generalized, it turns out that PT* is adapted to treat a large number of graphs. This is a fundamental point due to its undeniable advantage over exact calculations. The perturbative approach allows us to perform very fast calculations, even for graphs whose size is important, whereas exact calculations are restricted to small graphs because of the large size of the exciton-phonon Hamiltonian matrix. The superiority of PT* will be exploited in forthcoming works to describe exciton-mediated QTS or CTQWs in large networks.

## ACKNOWLEDGMENTS

Calculations have been performed thanks to the computing resources of the Mésocentre de Calcul, a Regional Computing Center at Université de Franche-Comté.

[1] G. D. Mahan, *Many-Particle Physics* (Kluwer Academics, New York, 2000).

[2] V. Agranovich, Zh. Eksp. Teor. Fiz. **37**, 430 (1959) [Sov. Phys. JETP **10**, 307 (1960)].

[3] A. S. Davydov, *Theory of Molecular Excitons* (Plenum, New York, 1971).

[4] R. Silbey, J. Chem. Phys. **46**, 4029 (1967).

[5] M. Grover and R. Silbey, J. Chem. Phys. **54**, 4843 (1971).

[6] V. May and O. Kuhn, *Charge and Energy Transfer Dynamics in Molecular Systems* (Wiley-VCH Verlag, Berlin, 2000).

[7] V. S. Lin, S. G. DiMagno, and M. J. Therien, Science **264**, 1105 (1994).

[8] G. L. Ingram, C. Nguyen, and Z.-H. Lu, Phys. Rev. Appl. **5**, 064002 (2016).

[9] J. Knoester and V. M. Agranovich, *Thin Films and Nanostructures: Electronic Excitations in Organic Based Nanostructures* (Elsevier, Amsterdam, 2003).

[10] P. Andreakou, S. V. Poltavtsev, J. R. Leonard, E. V. Calman, M. Remeika, Y. Y. Kuznetsova, L. V. Butov, J. Wilkes, M. Hanson, and A. C. Gossard, Appl. Phys. Lett. **104**, 091101 (2014).

[11] O. Mulken and A. Blumen, Phys. Rep. **502**, 37 (2011).

[12] V. Pouthier, J. Phys.: Condens. Matter **24**, 445401 (2012).

[13] H. Kamada and H. Gotoh, Semicond. Sci. Technol. **19**, S392 (2004).

[14] C. H. Bennet and D. P. DiVincenzo, Nature (London) **404**, 247 (2000).

[15] D. Loss and D. P. DiVincenzo, Phys. Rev. A **57**, 120 (1998).

[16] F. L. Semião, K. Furuya, and G. J. Milburn, New J. Phys. **12**, 083033 (2010).

[17] S. Bose, Phys. Rev. Lett. **91**, 207901 (2003).

[18] S. Bose, Contemp. Phys. **48**, 13 (2007).

[19] M. Christandl, N. Datta, A. Ekert, and A. J. Landahl, Phys. Rev. Lett. **92**, 187902 (2004).

[20] M. Plenio, B. Hartley, and J. Eisert, New J. Phys. **6**, 36 (2004).

[21] M. B. Plenio and F. L. Semião, New J. Phys. **7**, 73 (2005).

[22] C. Gollub, Ph.D. thesis, Ludwig Maximilian University of Munich, 2009.

[23] V. Pouthier, Phys. Rev. B **85**, 214303 (2012).

[24] O. Mülken, V. Bierbaum, and A. Blumen, J. Chem. Phys. **124**, 124905 (2006).

[25] A. M. Childs, Phys. Rev. Lett. **102**, 180501 (2009).

[26] S. E. Venegas-Andraca, Quant. Inf. Proc. **11**, 1015 (2012).

[27] A. Ambainis, SIAM J. Comput. **37**, 210 (2007).

[28] A. M. Childs, E. Farhi, and S. Gutmann, Quant. Inf. Proc. **1**, 35 (2002).

[29] E. Farhi, J. Goldstone, and S. Gutmann, Theory Comput. **4**, 169 (2008).

[30] A. M. Childs and J. Goldstone, Phys. Rev. A **70**, 022314 (2004).

[31] V. Pouthier, Phys. Rev. E **90**, 022818 (2014).

[32] P. Rebentrost, M. Mohseni, I. Kassal, S. Lloyd, and A. Aspuru-Guzik, New J. Phys. **11**, 033003 (2009).

[33] S. R. Jackson, T. J. Khoo, and F. W. Strauch, Phys. Rev. A **86**, 022335 (2012).

[34] A. L. Cardoso, R. F. S. Andrade, and A. M. C. Souza, Phys. Rev. B **78**, 214202 (2008).

[35] X.-P. Xu, W. Li, and F. Liu, Phys. Rev. E **78**, 052103 (2008).

[36] Z. Darázs, A. Anishchenko, T. Kiss, A. Blumen, and O. Mülken, Phys. Rev. E **90**, 032113 (2014).

[37] E. Agliari, A. Blumen, and O. Mülken, Phys. Rev. A **82**, 012305 (2010).

[38] O. Mülken, V. Pernice, and A. Blumen, Phys. Rev. E **76**, 051125 (2007).

[39] S. Salimi, Ann. Phys. (NY) **324**, 1185 (2009).

[40] X.-P. Xu, J. Phys. A **42**, 115205 (2009).

[41] A. Ziletti, F. Borgonovi, G. L. Celardo, F. M. Izrailev, L. Kaplan, and V. G. Zelevinsky, Phys. Rev. B **85**, 052201 (2012).

[42] A. Anishchenko, A. Blumen, and O. Mülken, Quant. Inf. Proc. **11**, 1273 (2012).

[43] V. Pouthier, Quant. Inf. Proc. **14**, 491 (2015).

[44] H. P. Breuer and F. Petruccione, *The Theory of Open Quantum Systems* (Oxford University Press, New York, 2007).

[45] S. M. Barnett and S. Stenholm, Phys. Rev. A **64**, 033808 (2001).

[46] M. Esposito and P. Gaspard, Phys. Rev. E **68**, 066112 (2003).

[47] V. Pouthier, J. Phys.: Condens. Matter **22**, 385401 (2010).

[48] V. Pouthier, Phys. Rev. B **83**, 085418 (2011).

[49] V. Pouthier, J. Chem. Phys. **134**, 114516 (2011).

[50] S. Yalouz and V. Pouthier, Phys. Rev. E **93**, 052306 (2016).

[51] S. Yalouz, C. Falvo, and V. Pouthier, Quant. Inf. Proc. **16**, 143 (2017).

[52] T. Holstein, Ann. Phys. (NY) **8**, 325 (1959); **8**, 343 (1959).
[53] Y. Tanimura and S. Mukamel, in *Ultrafast Dynamics of Chemical Systems*, edited by J. D. Simon (Kluwer, Dordrecht, 1994), p. 327.

[54] V. Pouthier, J. C. Light, and C. Girardet, *J. Chem. Phys.* **114**, 4955 (2001).

[55] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, *Atoms-Photon Interactions: Basic Processes and Applications* (Wiley, New York, 1992).

[56] P. D. Lax, *Functional Analysis* (Wiley-Interscience, New York, 2002).