# Quantum self-trapping on a star graph

Vincent Pouthier

► To cite this version:

Vincent Pouthier. Quantum self-trapping on a star graph. Physical Review E, 2022, 105 (4), pp.044304. 10.1103/PhysRevE.105.044304. hal-03466644

HAL Id: hal-03466644

https://hal.science/hal-03466644

Submitted on 6 Dec 2021

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
Quantum self-trapping on a star graph

Vincent Pouthier*

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

(Dated: December 6, 2021)

The attractive Bose-Hubbard model is applied for describing the two-exciton dynamics in a non-linear quantum star graph. When the excitons are created on the core of the star, it is shown that the interplay between the complex architecture of the network and the nonlinearity favors the occurrence of a real quantum self-trapping. Quite weak in the small nonlinearity limit, this self-localization is enhanced as the nonlinearity increases. This feature originates in the restructuring of the two-exciton eigenstates whose localized nature intensifies with the nonlinearity. Nevertheless, the quantum self-trapping is never complete since it is impossible to localize the entire exciton density, even in the strong nonlinearity limit.

PACS numbers:

# I. INTRODUCTION

Nonlinearity-induced energy localization in classical lattices has been intensively studied during the last four decades. This concept can be traced back to the seminal works of Davydov devoted to the vibrational energy flow in proteins¹. Using a quasi-classical approximation, Davydov suggested that a vibrational exciton propagates according as a soliton mechanism, a solution of the Nonlinear Schrodinger (NLS) equation in the continuum limit². In the mid-1980’s, lattice effects were introduced through the analysis of the discrete version of NLS. This equation revealed the occurrence of a remarkable phenomenon known as self-trapping³: the local accumulation of energy remains trapped where it has been created. Later, Sievers and Takeno⁴ showed that self-trapping is a special example of more general solutions called discrete breathers⁵. In classical anharmonic lattices, discrete breathers correspond to time-periodic and spatially localized solutions which result from the interplay between discreteness and nonlinearity.

In the quantum regime, a different situation occurs. Indeed, the quantum equivalent of the discrete NLS equation is the Bose version of the Hubbard model. This model has been used to study a great variety of situations such as Bose-Einstein condensates⁶, photonic quantum computers⁷ or vibrations in molecular lattices⁸. The Bose-Hubbard model describes interacting bosonic excitations (called excitons in the following of the text) evolving on a lattice with translational invariance. Therefore, the Bloch theorem applies so that the corresponding eigenstates cannot localize the energy. Nevertheless, the nonlinearity is responsible for the occurrence of specific states called multi-exciton bound states⁸⁻¹⁷. A bound state corresponds to the trapping of several excitons over only a few neighboring sites, with a resulting energy which is less than the energy of excitons lying far apart. The distance separating the excitons is small, so that they behave as a single particle delocalized along the lattice with a well-defined momentum. Although bound states cannot localize the energy because they must share

the symmetry of the translation operator, they take a very long time to tunnel from one lattice site to another. In other words, the initial creation of several excitons on a single site produces a localization of the energy over a timescale that increases with both the nonlinearity and the exciton number. This localized behavior, known as the quantum signature of the classical self-trapping, disappears in the long time limit due to the non-vanishing dispersion of the bound state energy band.

At present, because the occurrence of bound states in lattices with translational invariance is relatively well understood, our aim is to investigate what happens in complex networks one encounters in graph theory.

Indeed, it has been suggested recently that exploiting the motion of a single exciton in a complex network is a promising way for performing scalable quantum computing¹⁸. For instance, in a dendrimer, the propagation of an exciton corresponds to a physical realization of a continuous time quantum walk (CTQW)¹⁹. Extensively studied during the past few years, CTQW has become a very popular research subject due to its potential use in quantum information processing²⁰⁻²². For example, a CTQW on a complex network provides a natural way for performing efficient quantum searches in the spirit of the well-known Grover’s algorithm²³⁻²⁵. Consequently, CTQW and single exciton dynamics have been investigated in a great variety of networks such as extended dendrimers¹⁹,²⁶, binary and glued trees²⁷,²⁸, Apollonian networks²⁹, fractal networks³⁰,³¹, sequentially growing networks³² and star graphs³³⁻³⁹, to cite but a few examples.

Although some nonlinear effects in complex networks have been investigated using the discrete NLS equation⁴⁰,⁴¹, little is known about the quantum case. Therefore, in this paper, the concept of CTQW is extended to the case of several quantum walkers, i.e. to the case of several excitons moving on a complex network according to a Bose-Hubbard model. More precisely, to show that nonlinearity-induced quantum self-trapping may occur in complex networks, we will consider the situation in which two excitons are initially created
on the core of a star graph. The star graph is one of the most regular structures in graph theory. Organized around a central core, it exhibits the local tree structure of irregular and complex networks. However, its topology remains sufficiently simple so that analytical calculations can be carried out. The influence of both the exciton numbers and the graph architecture will be investigated in forthcoming works.

The paper is organized as follows. In Sec. II the star graph is introduced and the exciton Bose-Hubbard Hamiltonian is defined. Then the time dependent Schrodinger equation is established and the relevant ingredients required for characterizing the dynamics are described. The problem is solved numerically in Sec. III and the results are finally discussed and interpreted in Sec. IV.

## II. THEORETICAL BACKGROUND

### A. Model Hamiltonian

The star graph $S_N$ we consider is shown in Fig. 1. It corresponds to a tree that involves $N$ branches that emanate out from a central core. The central core, labeled by the index $\ell = 0$, is connected to $N$ branch sites $\ell = 1, ..., N$. Each site $\ell$ is occupied by a molecular subunit whose internal dynamics is described by an anharmonic oscillator. Let $b_{\ell}^{\dagger}$ and $b_{\ell}$ denote the corresponding standard boson operators. Within these notations, the exciton Hamiltonian is the Bose-Hubbard model defined as (with the convention $\hbar = 1$)

$$H = \sum_{\ell=0}^{N} \omega_0 b_{\ell}^{\dagger} b_{\ell} - Ab_{\ell}^{\dagger} b_{\ell}^{\dagger} b_{\ell} b_{\ell} + \sum_{\ell=1}^{N} \Phi(b_{0}^{\dagger} b_{\ell} + b_{\ell}^{\dagger} b_{0}), \quad (1)$$

where $\omega_0$ is the internal frequency of each oscillator, $\Phi$ represents the hopping constant between the core site and each branch site, and $A$ is the nonlinearity responsible for an attractive interaction between the excitons.

<figure><img src="image_2.png" /><figcaption>FIG. 1: The star graph.</figcaption></figure>

To characterize the exciton dynamics, the corresponding time dependent Schrodinger equation has to be solved. Since the Hamiltonian $H$ conserves the number of excitons, this can be achieved by using the number states

method¹². To proceed, the Hilbert space $E$ is partitioned into independent subspaces as $E = E_0 \oplus E_1 \oplus E_2 \oplus ..., E_v$ refers to the $v$-exciton subspace. The dimension of $E_v$ is equal to the number of ways for distributing $v$ indistinguishable quanta onto $N+1$ sites, i.e. $(v+N)!/(v!N!)$. Within this representation, the Hamiltonian is block-diagonal, each block corresponding to a particular exciton number.

In this paper, we focus our attention on the two-exciton dynamics. To proceed, a useful basis set to generate the entire $E_2$ subspace is given by the normalized and symmetrized states $|\ell; \ell'|$ with $\ell = 0, ..., N$ and $\ell' = \ell, ..., N$. A particular vector $|\ell; \ell'|$ characterizes two excitons located onto the sites $\ell$ and $\ell'$, respectively, as

$$|\ell; \ell'| = \begin{cases} b_{\ell}^{\dagger} b_{\ell'}^{\dagger} |\emptyset\rangle & \text{if } \ell' > \ell \\ \frac{1}{\sqrt{2}} b_{\ell}^{\dagger 2} |\emptyset\rangle & \text{if } \ell = \ell', \end{cases} \quad (2)$$

where $|\emptyset\rangle$ stands for the vacuum state. Note that the dimension of $E_2$ reduces to $(N+1)(N+2)/2$.

### B. Schrodinger Equation

In the local basis $|\ell; \ell'|$, the two-exciton quantum state is expanded as

$$|\Psi(t)\rangle = \sum_{\ell=1}^{N} \sum_{\ell'=\ell}^{N} \Psi_{\ell\ell'}(t) |\ell; \ell'\rangle. \quad (3)$$

Therefore, the time dependent Schrodinger equation depends on the nature of the basis vectors involved in so that different situations occur. When the two excitons are located on the core of the graph ($\ell = \ell' = 0$), the Schrodinger equation is expressed as

$$i\dot{\Psi}_{00} = (2\omega_0 - 2A)\Psi_{00} + \sum_{\ell=1}^{N} \sqrt{2}\Phi\Psi_{0\ell}. \quad (4)$$

When the first exciton is located on the core whereas the second exciton belongs to the periphery of the graph ($\ell = 0$ and $\ell' = 1, ..., N$), the Schrodinger equation is written as

$$i\dot{\Psi}_{0\ell'} = 2\omega_0\Psi_{0\ell'} + \sqrt{2}\Phi(\Psi_{00} + \Psi_{\ell'\ell'}) \\ + \Phi(\Psi_{1\ell'} + \dots + \Psi_{\ell'-1\ell'} + \Psi_{\ell'\ell'+1} + \dots \Psi_{\ell'N}). \quad (5)$$

When the two excitons occupy the same site of the periphery of the star ($\ell = \ell' = 1, ..., N$), the Schrodinger equation is written as

$$i\dot{\Psi}_{\ell\ell'} = (2\omega_0 - 2A)\Psi_{\ell\ell'} + \sqrt{2}\Phi\Psi_{0\ell'}. \quad (6)$$

Finally, when the two excitons are far apart and far from the core ($\ell = 1, ..., N$ and $\ell' = \ell + 1, ..., N$), the Schrodinger equation is expressed as

$$i\dot{\Psi}_{\ell\ell'} = 2\omega_0\Psi_{\ell\ell'} + \Phi(\Psi_{0\ell} + \Psi_{0\ell'}). \quad (7)$$
Eqs. (4)–(7) reveal the equivalence between the two-exciton dynamics and the dynamics of a single fictitious particle moving quantum mechanically on the complex graph displayed in Fig. 2. Within this equivalence, the two-exciton wave function $\Psi_{\ell\ell'}(t)$ can be viewed as the wave function of the fictitious particle. According to Eqs. (4)–(7), its dynamics is described by a tight-binding Hamiltonian characterized by self-energies located on each site and hopping matrices which couple different sites. The nonlinearity is responsible for the occurrence of defects in the graph leading to a shift of the corresponding self-energies. This equivalence allows us to simplify the resolution of the Schrodinger equation. Indeed, when the two excitons are initially created on the core of the star, Fig. 2 clearly shows the symmetry of the problem. It turns out that the wave function $\Psi_{\ell\ell'}(t)$ exhibits only 4 different values since both $\Psi_{0\ell}(t)$ and $\Psi_{\ell\ell}(t)$ are $\ell$ independent $\forall \ell = 1, N$ and $\Psi_{\ell\ell'}(t)$ is constant $\forall \ell = 1, \dots, N$ and $\forall \ell' = \ell + 1, \dots, N$.

<figure><img src="image_3.png" /><figcaption>FIG. 2: Equivalence between the two-exciton Schrodinger equation and the dynamics of a single fictitious particle moving quantum mechanically on a complex network (see the text).</figcaption></figure>

Consequently, Eqs. (4)–(7) can be solved by performing the following change of variables

$$ 
\begin{aligned}
\Psi_{00}(t) &= \Psi_{00}(t) \\
\chi_{0p}(t) &= \frac{1}{\sqrt{N}} \sum_{\ell=1}^{N} \Psi_{0\ell}(t) \\
\chi_{pp}(t) &= \frac{2}{\sqrt{N(N-1)}} \sum_{\ell=1}^{N} \sum_{\ell'=\ell+1}^{N} \Psi_{\ell\ell'}(t) \\
\Psi_{pp}(t) &= \frac{1}{\sqrt{N}} \sum_{\ell=1}^{N} \Psi_{\ell\ell}(t).
\end{aligned} 
\quad (8) 
$$

With these new variables, the Schrodinger equation reduces to a system of four equations expressed as

$$ 
\begin{aligned}
i\dot{\Psi}_{00} &= (2\omega_0 - 2A)\Psi_{00} + \sqrt{2N}\Phi\chi_{0p} \\
i\dot{\chi}_{0p} &= 2\omega_0\chi_{0p} + \sqrt{2N}\Phi\Psi_{00} + \sqrt{2}\Phi\Psi_{pp} \\
&\quad + \sqrt{2(N-1)}\Phi\chi_{pp} \\
i\dot{\chi}_{pp} &= 2\omega_0\chi_{pp} + \sqrt{2(N-1)}\Phi\chi_{0p} \\
i\dot{\Psi}_{pp} &= (2\omega_0 - 2A)\Psi_{pp} + \sqrt{2}\Phi\chi_{0p}.
\end{aligned} 
\quad (9) 
$$

Finally, to remove the coupling between $\chi_{0p}$ and $\chi_{pp}$, one introduces the new variables $\chi_{\pm} = (\chi_{0p} \pm \chi_{pp})/\sqrt{2}$ so that the Schrodinger equation is written as

$$ 
\begin{aligned}
i\dot{\Psi}_{00} &= (2\omega_0 - 2A)\Psi_{00} + \sqrt{N}\Phi(\chi_+ + \chi_-) \\
i\dot{\chi}_+ &= \epsilon_+\chi_+ + \sqrt{N}\Phi\Psi_{00} + \Phi\Psi_{pp} \\
i\dot{\chi}_- &= \epsilon_-\chi_- + \sqrt{N}\Phi\Psi_{00} + \Phi\Psi_{pp} \\
i\dot{\Psi}_{pp} &= (2\omega_0 - 2A)\Psi_{pp} + \Phi(\chi_+ + \chi_-),
\end{aligned} 
\quad (10) 
$$

with $\epsilon_{\pm} = 2\omega_0 \pm \sqrt{2(N-1)}\Phi$.

<figure><img src="image_4.png" /><figcaption>FIG. 3: Energy diagram of the 4-level system in which the two-exciton dynamics is confined.</figcaption></figure>

To summarize, when the two excitons are initially created on the core of the star, the quantum dynamics is isomorphic to that of a 4-level system whose energy diagram is shown in Fig. 3. This system involves the four orthogonal and normalized states $\{|\hat{\Psi}_{00}\rangle, |\hat{\chi}_+\rangle, |\hat{\chi}_-\rangle, |\hat{\Psi}_{pp}\rangle\}$ defined as

$$ 
\begin{aligned}
|\hat{\Psi}_{00}\rangle &= |0; 0\rangle \\
|\hat{\chi}_{\pm}\rangle &= \frac{1}{\sqrt{2}}(|\hat{\chi}_{0p}\rangle \pm |\hat{\chi}_{pp}\rangle) \\
|\hat{\Psi}_{pp}\rangle &= \frac{1}{\sqrt{N}} \sum_{\ell=1}^{N} |\ell; \ell\rangle,
\end{aligned} 
\quad (11) 
$$

with

$$ 
\begin{aligned}
|\hat{\chi}_{0p}\rangle &= \frac{1}{\sqrt{N}} \sum_{\ell=1}^{N} |0; \ell\rangle \\
|\hat{\chi}_{pp}\rangle &= \frac{2}{\sqrt{N(N-1)}} \sum_{\ell=1}^{N} \sum_{\ell'=\ell+1}^{N} |\ell; \ell'\rangle.
\end{aligned} 
\quad (12) 
$$
The state $|\hat{\Psi}_{00}\rangle$, whose energy is equal to $2\omega_0 - 2A$, defines an excitonic pair located on the core of the graph whereas the state $|\hat{\Psi}_{pp}\rangle$, with the same energy, describes an excitonic pair delocalized over the periphery of the star. These two pair states are coupled with the states $|\hat{\chi}_{\pm}\rangle$ that define superpositions of states involving excitons far from each other. Within this representation, the Hamiltonian of the 4-level system will be diagonalized numerically. The knowledge of the corresponding eigenvalues $\epsilon_{\mu}$ and eigenvectors $|\phi_{\mu}\rangle$, with $\mu = 1, ..., 4$, will allow us to solve the Schrodinger equation Eq. (10) and to compute the two-exciton quantum state as

$$|\Psi(t)\rangle = \Psi_{00}(t)|\hat{\Psi}_{00}\rangle + \chi_{+}(t)|\hat{\chi}_{+}\rangle + \chi_{-}(t)|\hat{\chi}_{-}\rangle + \Psi_{pp}(t)|\hat{\Psi}_{pp}\rangle, \quad (13)$$

with the initial condition $|\Psi(t = 0)\rangle = |\hat{\Psi}_{00}\rangle$.

### C. Observables

From the knowledge of both the two-exciton eigenstates and the two-exciton time dependent wave function, different observables can be computed.

First, we introduce the survival probability of the initial state $P_S(t) = |(0; 0|\Psi(t)\rangle|^2$. It defines the probability to observe the two excitons on the core of the star at time $t$, and characterizes the network memory of the initial localized state as

$$P_S(t) = |\Psi_{00}(t)|^2. \quad (14)$$

Then, information about the way the energy is distributed along the star is given by the expectation value of the population operator $\Pi_{\ell}(t) = \langle\Psi(t)|b_{\ell}^{\dagger}b_{\ell}|\Psi(t)\rangle$. The exciton population represents a key observable to describe the energy flow between the core and the periphery. It allows us to discriminate between both energy localization and delocalization. Therefore, in terms of the two-exciton wave function, the population at the core site and at time $t$, is expressed as

$$\Pi_0(t) = 2|\Psi_{00}(t)|^2 + \frac{1}{2}|\chi_+(t) + \chi_-(t)|^2. \quad (15)$$

Note that by symmetry, the exciton population at a peripheral site $\ell \neq 0$ is defined as $\Pi_{\ell}(t) = (2 - \Pi_0(t))/N$.

Finally, the physics of the two excitons is encoded in their eigenstates $|\phi_{\mu}\rangle$, with $\mu = 1, ..., 4$, whose nature strongly depends on the model parameters. Therefore, to characterize this feature, we define $P_{\mu,\alpha}$, with $\alpha = 00, +, -, pp$, as the weight of the eigenstate $|\phi_{\mu}\rangle$ in the basis $\{|\hat{\Psi}_{00}\rangle, |\hat{\chi}_{+}\rangle, |\hat{\chi}_{-}\rangle, |\hat{\Psi}_{pp}\rangle\}$.

### III. NUMERICAL RESULTS

In this section, the previous formalism is applied for describing the quantum dynamics when the two excitons

are initially created on the core of the star graph. Note that the convention $\Phi = 1$ will be used.

When $A = 0$ and $N = 9$, Fig. 4 reveals that both the exciton density and the survival probability correspond to periodic functions whose period is equal to $T_0 = 1.05\Phi^{-1}$. The exciton density varies between 0 and 2 (Fig. 4a) whereas the survival probability oscillates between 0 and 1 (Fig. 4b). Note that both functions reach their maximum value simultaneously. By contrast, we have verified that the probability to occupy the state $\chi_{0p}$ exhibits oscillations that ranges between 0 and 0.5 and whose period is equal to $T_0/2$ (note drawn). These results indicate that a perfect delocalization arises so that a coherent energy transfer takes place between the core and the periphery of the graph.

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

<figure><img src="image_6.png" /><figcaption>FIG. 4: Time evolution of (a) the exciton density and (b) the survival probability for $A = 0$ and $N = 9$.</figcaption></figure>

When $A = 1.0$ and $N = 9$, a different behavior occurs, as illustrated in Fig. 5. Initially equal to 2, the exciton density decreases as time increases, until it reaches a value equal to 0.085 at time $t = 0.525\Phi^{-1}$ (Fig. 5a). Then it increases until it reaches a value approximately equal to 1.68 at time $t = 1.035\Phi^{-1}$. Such a behavior
continues so that the density exhibits fast oscillations whose corresponding period is approximately equal to $T_0$. Nevertheless, these oscillations are modulated by a slowly varying envelope which prevents the density to vary between 0 and 2. Indeed, over the timescale considered in Fig. 5b, the minimum value of the density is equal to 0.084. However, recurrences take place almost periodically over a timescale defined by the period $T_1 = 10.8\Phi^{-1}$. Therefore, at time $t = T_1$, the density reaches a local maximum equal to 1.985. Similarly, at time $t = 2T_1$, the density reaches another local maximum equal to 1.940.

<figure><img src="image_7.png" /><figcaption>FIG. 5: Time evolution of (a) the exciton density and (b) the survival probability for $A = 1.0$ and $N = 9$.</figcaption></figure>

As illustrated in Fig. 5b, the survival probability behaves similarly, and it shows fast oscillations modulated by a slowly varying envelope. But a fundamental difference occurs. Although the exciton density no longer vanishes, the survival probability still reaches extremely low values. Over the present timescale, the minimum value of the survival probability is equal to $5 \times 10^{-5}$. These results indicate that although a coherent energy transfer still arises between the core and the periphery of the graph, approximately 4.25% of the initial energy remains

always trapped on the core site where it has been created. The difference between the density and the survival amplitude indicates that this trapping does not exclusively correspond to a confinement in the pair state localized on the core.

<figure><img src="image_8.png" /><figcaption>FIG. 6: Time evolution of (a) the exciton density and (b) the survival probability for $A = 3.0$ and $N = 9$.</figcaption></figure>

This is no longer the case when $A = 3.0$ and $N = 9$, as shown in Fig. 6. Indeed, the exciton density still exhibits fast oscillations, whose period is approximately equal to $T_0$, modulated by a slowly varying envelope (Fig. 6a). Initially equal to 2, recurrences take place almost periodically so that the density reaches local maximum equal to 1.98 and 1.95 at $t = 11.32\Phi^{-1}$ and $24.48\Phi^{-1}$, respectively. However, over the timescale shown in Fig. 6a, the minimum value of the density is equal to 0.47 indicating that more than 23.5% of the initial energy is trapped on the core site. This trapping effect now involves the pair state localized on the core site since the survival probability no longer vanishes as shown in Fig. 6b. Instead, it exhibits a minimum value equal to 0.053 over the timescale considered here.

As shown in Fig. 7 for $A = 10$ and $N = 9$, this localization process is clearly enhanced as $A$ increases. The
<figure><img src="image_9.png" /><figcaption>FIG. 7: Time evolution of (a) the exciton density and (b) the survival probability for $A = 10.0$ and $N = 9$.</figcaption></figure>

exciton density now behaves as a slowly varying function that supports a high-frequency small-amplitude modulation (Fig. 7a). The low-frequency component, whose period is approximately equal to $6.4\Phi^{-1}$, scales as a sine-like function that vary typically between 2 and 1. More precisely, over the timescale considered here, the maximum value of the exciton density is equal to 1.99 and it arises at time $t = 6.43\Phi^{-1}$. Similarly, the minimum value of the exciton density is equal to 1.16 and it arises at time $t = 28.48\Phi^{-1}$. The key point is that the survival probability behaves as the exciton density, and we have verified that the relation $\Pi_0(t) \approx 2P_S(t)$ is almost satisfied for all time (Fig. 7b). In other words, the energy transfer is now mediated by a pair state. Moreover, an important self-trapping arises since more than 50% of the initial energy stay localized on the core of the star.

In Fig. 8, a special attention is paid for describing the minimum value of both the exciton density on the core of the star $\Pi_0^*$ (Fig. 8a) and the survival probability $P_S^*$ (Fig. 8b). These observables have been extracted over a timescale equal to $100\Phi^{-1}$ for $N$ equal to 5, 10, 15 and

<figure><img src="image_10.png" /><figcaption>FIG. 8: A dependence of the minimum value of (a) the exciton density and (b) the survival probability over the timescale $100\Phi^{-1}$ and for different $N$ values.</figcaption></figure>

Equal to zero when $A=0$, the minimum value of the exciton density $\Pi_0^*$ increases as $A$ increases, indicating that a localization occurs on the core as soon as the nonlinearity turns on (Fig. 8b). However, the strength of this localization strongly depends on the nonlinearity and different regimes take place. For small $A$ values, $\Pi_0^*$ slowly increases with $A$ indicating that a quite weak localization arises. It depends quadratically on the nonlinearity provided that $A$ remains smaller than 2. Note that, the larger the size of the star $N$ is, the smaller is the minimum value of the exciton density. For larger $A$ values, $\Pi_0^*$ increases faster with the nonlinearity indicating that the localization is enhanced. It approximately scales linearly with $A$ according to a slope that is almost $N$ independent. Such a behavior persists until $\Pi_0^*$ reaches a maximum value that strongly depends on the size of the star. The larger the size of the star is, the larger is the maximum value reached by $\Pi_0^*$. This maximum is equal to 0.89, 1.33, 1.52, and 1.62 for $N = 5, 10, 15$ and $20$, respectively. In other words, the convergence of $\Pi_0^*$ towards a $A$ independent maximum value reveals that a strong, but incomplete, localization arises.
As shown in Fig. 8b, the minimum value of the survival probability $P_S^*$ behaves as $\Pi_0^*$ provided that $A$ is sufficiently important. Indeed, for small $A$ values, a fully different behavior occurs since $P_S^*$ remains at zero until $A$ reaches a critical value $A^*$. This critical value increases with the size of the star. It is approximately equal to 1.8, 2.2, 2.6 and 3.0 for $N = 5, 10, 15$ and $20$, respectively. For $A > A^*$, $P_S^*$ increases with the nonlinearity until it reaches a maximum and becomes $A$ independent. This maximum strongly depends on the size of the star and it is equal to 0.44, 0.66, 0.75, and 0.80 for $N = 5, 10, 15$ and $20$, respectively.

<figure><img src="image_11.png" /><figcaption>FIG. 9: A dependence of the two-exciton eigenenergies in the four-dimensional active subspace for $N = 9$.</figcaption></figure>

To understand these features, let us focus our attention on the nature of the two-exciton eigenstates. The $A$ dependence of the two-exciton eigenenergies in the four-dimensional active subspace are displayed in Fig. 9 for $N = 9$. When $A = 0$, the energy spectrum exhibits two non-degenerate states whose energies are defined as $2\omega_0 \pm 6\Phi$. In addition, the spectrum shows the energy $2\omega_0$ that is two-fold degenerate. Note that this result perfectly matches with the one-exciton properties of the star graph that supports two non-degenerate states whose energy is $\omega_{\pm} = \omega_0 \pm \sqrt{N}\Phi^{37}$. When two excitons are present, the corresponding eigenenergies are thus $2\omega_+, 2\omega_-, \omega_+ + \omega_- \text{ and } \omega_- + \omega_+$, as observed in Fig. 9.

When $A$ turns on, two distinct behaviors take place. First, the energy of the two lowest-energy states $\mu = 1$ and $\mu = 2$ decreases when $A$ increases. It is straightforward to show that $\epsilon_2 = 2\omega_0 - 2A$ and $\epsilon_1$ scales similarly for large $A$ values. By contrast, the energies $\epsilon_3$ and $\epsilon_4$ are slowly varying functions of the nonlinearity which slightly decreases as $A$ increases.

Fig. 10 shows the $A$ dependence of the weight $P_{\mu,\alpha}$, with $\alpha = 00, +, -$ and $pp$, of the eigenstates $|\phi_\mu\rangle$ in the basis $\{|\tilde{\Psi}_{00}\rangle, |\hat{\chi}_+\rangle, |\hat{\chi}_-\rangle, |\tilde{\Psi}_{pp}\rangle\}$.

When $A = 0$, Fig. 10a reveals that the lowest-energy

<figure><img src="image_12.png" /><figcaption>FIG. 10: A dependence of the two-exciton eigenstates in the four-dimensional active subspace for $N = 9$.</figcaption></figure>

state $\mu = 1$ corresponds basically to a superposition involving $|\hat{\chi}_-\rangle$ (69.44%) and $|\tilde{\Psi}_{00}\rangle$ (25%). As $A$ increases, an important restructuring arises. The weight of the state $|\hat{\chi}_-\rangle$ decreases whereas the weight of the states $|\tilde{\Psi}_{00}\rangle$ and $|\tilde{\Psi}_{pp}\rangle$ increases. Finally, for strong $A$ values, the state $\mu = 1$ almost localizes in the state $|\tilde{\Psi}_{00}\rangle$ (89.5%), exhibiting a rather small contribution of the state $|\tilde{\Psi}_{pp}\rangle$ (0.09%).

As illustrated in Fig. 10b, the structure of the state $\mu = 2$ is quite surprising. Indeed, when $A = 0$, it corresponds to a degenerate state fully delocalized over the 4 basis states. However, when $A$ turns on, the degeneracy is removed. Therefore, the state $\mu = 2$ becomes $A$ independent. It almost localizes in the state $|\tilde{\Psi}_{pp}\rangle$ (90.0%), exhibiting a rather small contribution of the state $|\tilde{\Psi}_{00}\rangle$ (0.10%).
When $A=0$, the state $\mu=3$ defines the second degenerate state that is delocalized over the 4 basis states (Fig. 10c). Then, for small $A$ values, it basically involves $|\hat{\Psi}_{00}\rangle$ and $|\hat{\chi}_{\pm}\rangle$. However, an important restructuring arises as $A$ increases. The weight of the states $|\hat{\chi}_{-}\rangle$ increases whereas the weight of the remaining states decreases. Consequently, for strong $A$ values, the state $\mu=3$ almost reduces to $|\hat{\chi}_{-}\rangle$ whose weight reaches 99.65% for $A=30$.

As shown in Fig. 10d, when $A=0$ the state $\mu=4$ corresponds basically to a superposition involving $|\hat{\chi}_{+}\rangle$ (69.44%) and $|\hat{\Psi}_{00}\rangle$ (25%). As $A$ increases, the weight of the states $|\hat{\chi}_{+}\rangle$ increases whereas the weight of the remaining states decreases. As a result, the state $\mu=3$ almost reduces to $|\hat{\chi}_{+}\rangle$ for strong $A$ values. For instance, the weight of the state $|\hat{\chi}_{+}\rangle$ represents 99.71% for $A=30$.

## IV. DISCUSSION

Our numerical results reveal a fundamental feature. Indeed, in a star graph, the Bose-Hubbard model allows the occurrence of a real quantum self-trapping as soon as the nonlinear parameter turns on. Such a phenomenon is quite remarkable because it does not appear in lattices with translational invariance. In other words, it is the interplay between the complex architecture of the network and the nonlinearity that gives rise to the self-localization of the energy. Nevertheless, the quantum self-trapping is not perfect in the sense that even for a very strong nonlinearity, it is impossible to localize the entire exciton density on the core site.

The key point is that the self-trapping phenomenon remains quite subtle because it strongly depends on the strength of the nonlinearity. As shown in Fig. 3, such a sensitivity results from the position of the energy level of the excited localized pair state $|\hat{\Psi}_{00}\rangle$ with respect to the other energy levels. This position controls the nature of the eigenstates that, in turn, govern the dynamics. This feature allows us to introduce a critical nonlinearity $A^* = \sqrt{2(N-1)\Phi}/2$ for which there is a resonance between the pair states $|\hat{\Psi}_{00}\rangle$ and $|\hat{\Psi}_{pp}\rangle$ and the state $|\hat{\chi}_{\pm}\rangle$ (see Fig. 3). Note that this value basically corresponds to the critical value identified in Fig. 8b. Indeed, in a rather good agreement with the numerical results, one obtains $A^*/\Phi = 1.41, 2.12, 2.64$ and 3.08 for $N=5, 10, 15$ and 20, respectively.

When the nonlinearity is rather weak, i.e. when $A < A^*$, the localized pair state $|\hat{\Psi}_{00}\rangle$ interacts with all the other basis states. As a consequence, the initial state that follows the creation of the two excitons on the core site decomposes almost over all the eigenstates. Due to the nature of these eigenstates, such a situation favors the transfer of the localized pair state $|\hat{\Psi}_{00}\rangle$ towards the other basis states $|\hat{\chi}_{\pm}\rangle$ and $|\hat{\Psi}_{pp}\rangle$ giving rise to the delocalization of the exciton density from the core to the periphery of the star. Nevertheless, this density exhibits

two contributions, i.e. $\Pi_0(t) = 2|\Psi_{00}(t)|^2 + |\chi_{0p}(t)|^2$. The first contribution involves the population $|\Psi_{00}(t)|^2$ of the pair state localized on the core whereas the second contribution refers to the population $|\chi_{0p}(t)|^2$ of the state involving an exciton localized on the core and an exciton delocalized over the periphery. These two populations vary as time elapses and they vanish for specific times indicating that a delocalization arises. But the key point is that these two populations never vanish simultaneously. As a consequence, the exciton density on the core of the star is always greater than zero resulting in a self-localization of the energy.

As the nonlinearity increases, the self-trapping is enhanced. Nevertheless, its origin begins to change as one reaches the resonance for $A = A^*$.

When $A > A^*$, the quantum self-trapping originates in the restructuring of the eigenstates whose localized nature is intensified. More precisely, a key role is played by the two lowest-energy eigenstates $\mu=1$ and $\mu=2$ that basically correspond to superpositions involving the pair states. The eigenstate $\mu=1$ refers to the pair state localized on the core site and it exhibits a small contribution of the pair state delocalized over the periphery. In turn, the state $\mu=2$ refers to the pair state delocalized over the periphery and it exhibits a small contribution of the pair state localized on the core site. Consequently, the initial creation of the two excitons on the core site mainly excites the eigenstate $\mu=1$, as well as the state $\mu=2$ to a lesser extent. Therefore, the population $|\Psi_{00}(t)|^2$ oscillates around an important mean value whereas the population $|\Psi_{pp}(t)|^2$ oscillates around a quite small mean value. A rather strong self-trapping arises that mainly corresponds to the localization of the energy in the pair state localized on the core of the star. But this self-trapping is never complete because the population of the pair state delocalized over the periphery never vanishes.

To interpret that the eigenstate restructuring favors the self-trapping, let us apply to the star graph the model developed in Refs$^{16,17}$. This model is based on the numerical observation that the dynamics is confined in a relevant subspace provided that $A$ is sufficiently strong. The dimension of the problem is thus reduced so that the model provides a simple view of the dynamics in which pair states play the central role. According to Fig. 2, the relevant subspace is generated by the number states $|0;0\rangle$, $|0;\ell\rangle$ and $|\ell;\ell\rangle$, with $\ell=1, \dots, N$, only. Disregarding the influence of the remaining states, the representation of $H$ reduces to a tight-binding model on an extended star graph. Both states $|0;0\rangle$ and $|\ell;\ell\rangle$ exhibit the same energy $\epsilon_{00} = \epsilon_{pp} = 2\omega_0 - 2A$ whereas the energy of the states $|0;\ell\rangle$ is $2\omega_0$.

Therefore, when the two excitons are created on the core site, the dynamics reduces to that of pair states whose properties are modified due to their interactions with the intermediate states $|0;\ell\rangle$. According to standard perturbation theory, these interactions are responsible for the following features. First, due to the coupling with the $N$ states $|0;\ell\rangle$, the self-energy of the
<figure><img src="image_13.png" /><figcaption>FIG. 11: Effective Hamiltonian that govern the pair states dynamics in the strong A limit. Full circles stand for the renormalized self-energy $\tilde{\epsilon}_{pp}$ of the peripheral pair states $|\ell; \ell\rangle$ whereas the open circle defines the renormalized self-energy $\tilde{\epsilon}_{00}$ of the localized pair states $|0; 0\rangle$. Dashed lines represent the effective coupling $J = -\Phi^2/A$ between the pair state localized on the core and the peripheral pair states (see the text).</figcaption></figure>

localized pair state $|0; 0\rangle$ is renormalized and becomes $\tilde{\epsilon}_{00} = 2\omega_0 - 2A - N\Phi^2/A$. Second, because the peripheral pair state $|\ell; \ell\rangle$ interacts with one intermediate state $|0; \ell\rangle$, its self-energy is also modified according to $\tilde{\epsilon}_{pp} = 2\omega_0 - 2A - \Phi^2/A$. Then, an effective coupling arises between $|0; 0\rangle$ and each pair state $|\ell; \ell\rangle$ whose intensity is $J = -\Phi^2/A$.

Consequently, it is as if the excitonic pair behaves as a single particle moving on a the star graph shown in Fig. 11, this graph exhibiting a defect on its core site. In that case, different strategies have been developed for describing the corresponding eigenstates$^{33-35}$. Here, we take advantage of the fact that the graph possesses discrete rotational symmetry. It remains invariant under the discrete rotation of angle $\theta_0 = 2\pi/N$ and centered on the core site. The diagonalization of the Hamiltonian is thus greatly simplified when one works with the Bloch basis$^{35}$ that involves the local state $|0; 0\rangle = |\hat{\Psi}_{00}\rangle$ and $N$ orthogonal Bloch states $|\hat{\Psi}_k\rangle$ ($k = 1, ..., N$) defined as

$$|\hat{\Psi}_k\rangle = \frac{1}{\sqrt{N}} \sum_{\ell=1}^{N} e^{ikl\theta_0} |\ell; \ell\rangle \quad (16)$$

Within the Bloch basis, the graph exhibits two kinds of eigenstates. First, the spectrum shows the $(N-1)$-fold degenerate eigenenergy $\tilde{\epsilon}_{pp}$, the corresponding eigenstates being the $N-1$ Bloch states $|\hat{\Psi}_k\rangle$, with $k = 1, ..., N-1$. These states do not play any role in the present situation. Second, the graph supports two eigenstates that correspond to superpositions involving the state $|0; 0\rangle = |\hat{\Psi}_{00}\rangle$ and the Bloch state $|\hat{\Psi}_{k=N}\rangle = |\hat{\Psi}_{pp}\rangle$ that is uniformly distributed over the periphery of the star. These eigenstates, which govern the dynamics when two excitons are created on the core site, are defined as

$$|\phi_1\rangle = (\sqrt{N}|\hat{\Psi}_{00}\rangle + |\hat{\Psi}_{pp}\rangle)/\sqrt{N+1} \quad (17) \\ \epsilon_1 = 2\omega_0 - 2A - (N+1)\Phi^2/A,$$

and

$$|\phi_2\rangle = (|\hat{\Psi}_{00}\rangle - \sqrt{N}|\hat{\Psi}_{pp}\rangle)/\sqrt{N+1} \quad (18) \\ \epsilon_2 = 2\omega_0 - 2A.$$

In a perfect agreement with our numerical results, Eqs. (18) and (19) define the two lowest-energy states that govern the dynamics in the strong A limit. The eigenstate $\mu = 1$ defines a pair state mainly localized on the core site whereas the state $\mu = 2$ is mainly delocalized over the periphery. The nature of these two states is due to the presence of a defect on the star graph that describes the pair state dynamics. This defect originates in the energy correction of the pair states that differs depending on whether the pair occupies the core or the periphery of the star. In other words, because the degree of the core site is $N$ times larger than the degree of a peripheral site, the correction of the self-energy of the localized pair state $|0; 0\rangle$ is $N$ times larger than that of the pair state $|\ell; \ell\rangle$. This feature favors the localization of the eigenstates that gives rise to the quantum self-trapping phenomenon observed in the previous section.

Form the knowledge of Eqs. (18) and (19), standard quantum mechanism calculations allow to get an analytical expression of the survival probability to observe the excitons in the localized pair state $|0; 0\rangle$ at time $t$. This probability is written as

$$P_S(t) = 1 - \frac{4N}{(N+1)^2} \sin^2 \left( \frac{(N+1)\Phi^2 t}{2A} \right) \quad (19)$$

As observed in the numerical section, the survival probability oscillates around an average value $\bar{P}_S = 1 - 2N/(N+1)^2$. This value increases as $N$ increases indicating that the degree of the central core enhances the quantum self-trapping effect. Note that in a perfect agreement with the observations in Fig. 8b, Eq. (19) yields a minimum value of the survival probability equal to 0.44, 0.67, 0.76 and 0.82 for $N = 5, 10, 15$ and 20, respectively.

## V. CONCLUSION

In the present paper, the Bose-Hubbard model has been used to analyze the energy transfer in a nonlinear quantum star graph. Within this model, the dynamics is controlled by two relevant parameters, i.e. the nonlinearity $A$ and the degree $N$ of the core site. When two excitons are initially created on the core site, it has been shown that a real quantum self-trapping occurs. Such a phenomenon is quite surprising because it does not appear in lattices with translational invariance. In fact, the self-localization of the energy results from the interplay between the complex architecture of the network and the nonlinearity. Rather weak in the small nonlinearity limit, the self-trapping is enhanced as the nonlinearity increases due to the restructuring of the two-exciton eigenstates whose localized nature intensifies. Nevertheless, the quantum self-trapping is never complete since it is
impossible to localize the entire exciton density, even in
the strong nonlinearity limit.

This work, which can be viewed as a first step, falls
within a more general framework devoted to the study
of nonlinear quantum complex networks. The next steps
will concern the characterization of different features ex-
pected to play a crucial role. For instance, it would be
interesting to analyze the influence of the initial condi-
tions, i.e. to create initially either an excitonic pair on the

periphery of the star or two excitons lying far apart. The
self-trapping could compete with degeneracy-induced lo-
calization that arises in star graphs³⁷,³⁸. Then, we know
that the exciton number is also a key ingredient that
may enhances self-trapping effects. Finally, the natural
pursuit of these researches would be to investigate what
happens in more complex networks such as extended star
graph, dendrimers or glued trees, to cite but a few exam-
ple[s].

* Electronic address: vincent.pouthier@univ-fcomte.fr

1 A.S. Davydov and N.I. Kisluka, Phys. Status Solidi **59**, 465 (1973) ; Zh. Eksp. Teor. Fiz **71**, 1090 (1976) ; Sov. Phys. JETP **44**, 571 (1976).

2 A.C. Scott, Phys. Rep. **217**, 1 (1992).

³ J.C. Eilbeck, P.S. Lomdahl, and A.C. Scott, Physica D **16**, 318 (1985).

⁴ A.J. Sievers and S. Takeno, Phys. Rev. Lett. **61**, 970 (1988).

⁵ S. Flach and C.R. Willis, Phys. Rep. **295**, 181 (1998).

⁶ G. Kalosakas, A.R. Bishop, and V.M. Kenkre, Phys. Rev. A **68**, 023602 (2003).

⁷ S. Yalouz, B. Senjean, F. Miatto, and V. Dunjko, Quantum **5**, 572 (2021).

⁸ J.C. Kimball, C.Y. Fong, and Y.R. Shen, Phys. Rev. B **23**, 4946 (1981).

⁹ F. Bogani, G. Cardini, V. Schettino, and P.L. Tasselli, Phys. Rev. B **42**, 2307 (1990).

10 E. Wright, J.C. Eilbeck, M.H. Hays, P.D. Miller, and A.C. Scott, Physica D **69**, 18 (1993).

11 L. Bernstein, J.C. Eilbeck, and A.C. Scott, Nonlinearity **3**, 293 (1990).

12 A.C. Scott, J.C. Eilbeck, and H. Gilhoj, Physica D **78**, 194 (1994).

13 V. Pouthier, Phys. Rev. E **68**, 021909 (2003).

14 J. Dignac, J.C. Eilbeck, M. Salerno, A.C. Scott, Phys. Rev. Lett. **93**, 025504 (2004).

15 L. Proville, Phys. Rev. B **71**, 104306 (2005).

16 C. Falvo, V. Pouthier, and J.C. Eilbeck, Physica D **221**, 58 (2006).

17 V. Pouthier, Phys. Rev. B **76**, 224302 (2007).

18 O. Mulken and A. Blumen, Phys. Rep. **502**, 37 (2011).

19 O. Mulken, V. Bierbaun and A. Blumen, J. Chem. Phys. **124**, 124905 (2006).

20 A.M. Childs, Phys. Rev. Lett. **102**, 180501 (2009).

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

22 A.M. Childs, E. Farhi, and S. Gutmann, Quantum Inf. Proc. **1**, 35 (2002).

23 E. Farhi and S. Gutmann, Phys. Rev. A **57**, 2403 (1998).

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

25 L. K. Grover, Phys. Rev. Lett. **79**, 325 (1997).

26 V. Pouhtier, J. Chem. Phys. **139**, 234111 (2013).

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

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

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

30 Z. Darazs, A. Anishchenko, T. Kiss, A. Blumen, and O. Mulken, Phys. Rev. E **90**, 032113 (2014).

31 E. Agliari, A. Blumen, and O. Mulken, Phys. Rev. A **82**, 012305 (2010).

32 O. Mulken, M. Dolgushev, and M. Galiceanu, Phys. Rev. E **93**, 022304 (2016).

33 S. Salimi, Ann. Phys. **324**, 1185 (2009).

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

35 X.P. Xu, Phys. Rev. E **79**, 011117 (2009).

36 A. Anishchenko, A. Blumen, and O. Mulken, Quantum Inf. Process **11**, 1273 (2012).

37 V. Pouthier, Quantum Inf. Process **14**, 491 (2015).

38 V. Pouthier, Quantum Inf. Process **14**, 3139 (2015).

39 S. Yalouz, V. Pouthier, and C. Falvo, Quantum Inf. Process **16**, 143 (2017).

40 G.P. Tsironis, Phys. Lett. A **375**, 676 (2011).

41 G.P. Tsironis, Phys. Lett. A **375**, 1304 (2011).