text:
# Mode Entanglement and Entangling power in Bosonic Graphs

## I introduction

Generation of quantum entanglement is an essential prerequisite for the majority of quantum information and computation tasks qip . To this aim many-body interactions have to be cleverly engineered and carefully controlled in the variety of physical systems that pretend to represent viable candidates for quantum information processors. The resulting quantum dynamics can be then characterized by means of its entangling capabilities in different ways ep ; ep\_cirac ; lea ; man . These entangling power measures provide a quantitive to assess the dynamical strength man  of a quantum evolution as physical resource for quantum information processing (QIP) protocols. In this respect the optimization of the entangling power in the set of physically available interactions clearly represents a quite relevant goal.

In this paper we will address this problem in the context of bosonic systems on lattice structures e.g., ultracold bosonic atoms in optical lattices. This kind of systems have been already rather extensively studied in the QIP related literature milb\_twowells ; bec\_lattice ; chen ; radu\_pz ; duan\_ent ; simon ; hines ; micheli ; you  in view of their experimental feasibility and inherent great flexibility that make them almost ideal candidates for studing and simulating duan\_ctrl  a plethora of important quantum many-body phenomena e.g., quantum phase transitions bec\_mott .

In the following we will adopt the view of quantum entanglement in systems of indistinguishable particles advocated in Refs virt ; ind\_ent ; gitt ; enk ; shi ; vedra  in alternative to the complementary one pursued in Refs. schli ; pask ; Li . In the former the subsystems that get entangled are provided by bosonic modes whereas in the latter are particles themselves. The mode-entanglement approach is an intrinsecally second-quantized one: the tensor product structure of the state-space necessary in order to define entanglement is provided the natural identification of the bosonic (fermionic) Fock space with a set of linear oscillators (qubits) associated with the single-particle wavefunctions.

We shall analyze systems of bosonic particles hopping over inhomogneous spatial structures bur ; zan\_graph  mathematically modelled by graphs graph . The main objective is to determine the role that the graph topology plays in the entangling power of the quantum dynamics described by a Bose-Hubbard model bh . For a given tunneling graph structure we will study how the entangling capabilities are affected by the amount of non-linearity i.e., bosons self-interactions present in the Hamiltonian. Numerical simulations as well as perturbative analytical arguments will be provided for all the graphs up to four vertices. This analysis allows to optimzes the entanglement production, from an initial product state, obtainable by switching on-and off the tunneling between the various garph vertices i.e., one-body operator control, and by varying e.g., by Fesdbach resonances, the strength of the boson-boson local coupling. We believe that the theoretical framework developed in this paper will have direct experimental relevance.

The paper is organized as follows: in Sect II is briefly recalled the Bose-Hubbard Hamiltonian and in Sect III it will be applied to the case of hopping over spatially inhomogeneus structures (graphs). In Sect IV we will describe the notion of bosonic mode-entanglement for the systems of interest and introduce the functional for measuring the entangling power of a quantum evolution. Sect V is the major one, it contains the results of numerical simulations for all the graphs up to four vertices. Finally Sect VI provides a summary and conclusions.

## II The Bose-Hubbard model

The systems we will study in this paper are all given by $`L`$ bosonic sites filled with $`N=L`$ bosonic particles. Some of the sites are linked to each other in the sense that the particles are allowed to tunnel from one site to another. The links existing between the sites will be described later (section III) by means of graphs. The particles in a site interact via a two-body interaction process (self-interaction).

The proper framework that can be used to describe the systems under study is given by the Bose-Hubbard (BH) model bh . This model has been recently use to investigate many interesting properties of different systems going from arrays of BEC where the sites form an infinite (periodic) lattice bec\_lattice , to systems in which the number of sites is finite (dimer, trimer)milb\_twowells . Moreover this model has recently been used to describe new possible schemes for quantum computing radu\_pz . According to this model the Hamiltonian of the system of $`L`$ sites and $`N`$ particles can be written as:

$$H=\underset{i}{}\epsilon _in_i^2+\underset{i}{}\gamma _in_i+\underset{<i,j>}{}\tau _{ij}(c_i^{}c_j+h.c.)$$
(1)
where: the indexes $`i,j`$ run in $`V=\{0,1,2\mathrm{},L1\}`$, the set of sites; $`c_i^{}\text{ and }c_i`$ are bosonic creation and annihilation operators that create/destroy a particle in the locale mode (site) $`i`$; $`n_i=c_i^{}c_i`$ are the corresponding occupation number operators.

The parameters $`\gamma _i`$ can be used to model the offsets in the ground state energies of the different sites. We suppose $`\gamma _i=\gamma =\text{const}i`$ and from now on we neglect this term. In fact, since the total number of particles in our systems will be held fixed, the term $`\gamma _in_i=\gamma N`$ commutes with the Hamiltonian (1) and its contribution in the dynamics of the system consists only by a global phase. The parameters $`\epsilon _i`$ control the non-linear two-body interaction between the particles in each site $`i`$, while the parameters $`\tau _{ij}`$ control the tunneling processes between the sites $`i,j`$. We consider only tunneling between next-neighbourghring sites and we suppose to be able to set the tunnel parameters $`\tau _{ij}`$ and the self-interaction parameters $`\epsilon _i`$ to some fixed desired values. These will be held fixed during the evolution, i.e. we will be intersted only on the free evolution of our systems. The latter will be described in detail in the next section.

## III GRAPHS AND HAMILTONIANS

As described in the previous section, each of the systems we consider is given by a lattice of $`L`$ vertices (sites or modes) filled with $`N=L`$ bosonic particles and we suppose that its Hamiltonian can be described in the framework of the Bose Hubbard model. We consider systems in which the self-interaction terms are all equal ($`\epsilon _i=\epsilon i`$), while for the hopping parameters we consider the case $`\tau _{ij}\{0,\tau \}i,j`$. We focus on the free evolution of the simpliest of these systems, i.e. $`N=L\{3,4\}`$.

For $`L=N=3,4`$ we consider only inequivalent systems, in the sense specified below, each characterized by a different rooted graph $`\mathrm{\Gamma }=(V,E,r)`$, where $`V=\{0,1,2,..,L1\}`$ is the set of vertices, $`E=\{(i,j):i,jV,ij\}`$is the set of edges and $`rV`$ is the root of the graph. The latter will be considered as the mode (vertex) which is relevant for the calculation of the mode entanglement. For each rooted graph $`\mathrm{\Gamma }`$ the pair $`(V,E)`$ identifies a particular Bose Hubbard hamiltonian $`H_{(V,E)}`$, since the pair $`(V,E)`$ identifies a set of sites and the way in which these sites are connected. With this Hamiltionian one can build the evolution operator $`U_{(V,E)}(t)=\mathrm{exp}(itH_{(V,E)})`$.

If the pairs $`(V,E)`$ and $`(V,E^{})`$ are such that $`EE^{}`$ they give rise to two sets of inequivalent rooted graphs since the connectivity of the two graphs they identify is different. Morover, for a given pair $`(V,E)`$ one can define a certain rooted graph by specifying the root site $`rV`$, but only some of the possible choices can be considered inequivalent. If two rooted graphs differ only for the root, they will be considered inequivalent if there is a difference in how the roots are directly connected to the other vertices of the graph or if there is a difference in the connectivity of the subgraph $`\mathrm{\Gamma }^{}=(V^{},E^{})`$ constituted by the set of vertices $`V^{}=V\backslash r`$ and the set of edges $`E^{}=\{(i,j):i,jV^{},ij\}`$.

We now give an example of the procedure used to identify the different inequivalent systems by starting with $`N=L=3`$. In this case we will consider three different inequivalent rooted graphs (see Fig. 1a).

To build these graphs we first start by labeling the root vertex with 0, the left $`(L)`$ vertex with 1 and the right $`(R)`$ vertex with 2. The choice of the labels for the vertices of the subgraph $`\mathrm{\Gamma }^{}=(V^{},E^{})`$, where $`V^{}=\{L,R\}`$ is irrelevant for our analysis (we could for example choose $`L=2`$ and $`R=1`$) since from the point of view of the entanglement of the root mode they are interchangeable and the relevant feature is how $`r`$ is directly connected to the vertices of the subgraph $`\mathrm{\Gamma }^{}`$ and the connectivity of the latter. So we can say that, for a given $`E^{}`$,the two labeling $`(0,1,2)`$ and $`(0,2,1)`$ are equivalent and consider only the case $`(r,L,R)(0,1,2)`$.

The rooted graphs $`\mathrm{\Gamma }_0`$ and $`\mathrm{\Gamma }_2`$ give rise to the same ”kind” of BH Hamiltonians, in the sense that they both have two hopping terms; the only thing that changes is the labeling of the modes:

$$\begin{array}{c}H_{\mathrm{\Gamma }_0}=\epsilon (n_r^2+n_1^2+n_2^2)+\tau (c_r^{}c_1+c_r^{}c_2+\text{(}h.c.))\\ (r1)\\ H_{\mathrm{\Gamma }_2}=\epsilon (n_r^2+n_1^2+n_2^2)+\tau (c_r^{}c_1+c_1^{}c_2+\text{(}h.c.)).\end{array}$$

But, from the point of view of the root site, in $`\mathrm{\Gamma }_2r`$ is directly linked to one site (1), whereas in $`\mathrm{\Gamma }_0r`$ is directly linked with two sites (1,2); moreover the subgraphs have different connectivity: $`E_{\mathrm{\Gamma }_2}^{}\{(1,2)\}E_{\mathrm{\Gamma }_0}^{}\mathrm{}`$.

The Hamiltonian based on $`\mathrm{\Gamma }_1`$ differs from the other two since there are three hopping terms. In this case $`r`$ is directly connected to all the other sites as in $`\mathrm{\Gamma }_0`$ but here the connectivity of the subgraph is different, $`E_{\mathrm{\Gamma }_1}^{}=\{(1,2)\}`$.

Among all the possible graphs for the case $`N=L=3`$ we do not consider disconnected graphs like the ones in Fig.1b. The first case is a particular instance of the system $`L=2`$. In the second case, since in our analysis the initial state will always be a product state belonging to the site occupation number basis, no entanglement can be genertated between the root mode and the others.

Following the same ideas one can derive the inequivalent $`\mathrm{\Gamma }`$’s also for the case $`N=L=4`$. One can start by choosing the root site and connect it to just one of the other vertices; then one has to find all the subgraphs $`\mathrm{\Gamma }^{}`$ that give rise to inequivalent rooted graphs and repeat the whole procedure for the cases in which the root is directly connected with two and three of the other vertices. The result of this process is given in Fig.2. Once selected the relevant rooted graphs $`\mathrm{\Gamma }_j`$ one can build the relative Hamiltonians $`H_{\mathrm{\Gamma }_j}(\epsilon ,\tau )`$ and the relative evolution operators $`U_{\mathrm{\Gamma }_j}(t)=\mathrm{exp}(itH_{\mathrm{\Gamma }_j}(\epsilon ,\tau ))`$. The entagling power of the latter can then be calculated as described in the next section.

## IV MODE ENTAGLEMENT

The problem of the entanglement of two-mode (component) Bose Einstein condensates (BEC) has been recently studied in different ways milb\_twowells ; micheli ; duan\_ent ; hines . In Ref. micheli  the system is composed by a two component BEC consisting of N atoms that can access two different hyperfine states $`|A`$ and $`|B`$ coupled by an electromagnetic field. In this paper it is stated that the maximally entangled state for the system is given by

$$|GHZ_N=\frac{1}{\sqrt{2}}(|A^^N+|B^^N)$$

which is the N-partite generalization of the GHZ (Greenberg- Horne-Zeilinger) state for 3 qubits

$$|GHZ_N=\frac{1}{\sqrt{2}}(|000+|111).$$

The underlying decomposition of the system is defined by considering each of the N bosons as a subsystem living in the two dimensional Hilbert space given by $`\text{span}\{|A,|B\}`$.

In this picture to each boson corresponds a qubit and the Hilbert space of the total system then corresponds to the Hilbert space of N qubits. On the other hand, as it has been pointed out in hines , the idea to consider different bosons as different two-levels subsystems and to calculate the value of the entanglement according to some suitable measure, althought mathematically possible, it is physically questionable in view of the indistinguishability of the particles. In order to have some entanglement between two systems, these have to be separately accessible from the experimental point of view.

In general one can think to speak of the entanglement between indistinguishable particles, but it is not possible to make any use of this resource since there are no physically realizable measurements on the system that can discriminate between the particles.

In milb\_twowells  one of the systems under study is a pair of tunnel-coupled BEC’s considered as a realization of a bimodal BEC. In this paper, On the contrary of what is done in Ref micheli , in order to calculate the entanglement present in the ground state of the system, the subsystems that are considered physically meaningful are the two modes. The number of bosons present in the system is always fixed. In this set up it is possible to have physical access to the subsystems. For example, one realization of the model is given by a BEC with a fixed number of bosonic particles that are trapped in a two wells potential, where each well (site) can be considered as a mode; here the modes are spatially separated. The indistinguishability of the particles does not allow to determine experimentally in which mode is a specific particle, but it is physically possible to measure the occupation number of each mode, i.e. the number of particles present at a given instant in the wells. The system can be viewed as a bipartite system for which a good measure of the entanglement exist: the Von Neumann entropy of the reduced density matrix. In this paper we follow this point of view and we make use of the notion of mode entanglement to characterize our systems.

We consider systems in which the number of modes (sites) $`L`$ is fixed and is always equal to the number of particles $`N\{3,4\}`$. The Hilbert space $`H^{(N)}`$ of the system for a fixed $`N=L`$ is a subspace of the Fock space $`h^L`$, where $`h=\text{span}\{|n\}_{n\mathrm{I}\mathrm{N}}`$. The subspace can be written as

$$H^{(N)}=\text{span}\{_{j=0}^{L1}|n_j=|\stackrel{}{n}:\stackrel{}{n}\mathrm{\Lambda }^N\}$$

where

$$\mathrm{\Lambda }^N=\{\stackrel{}{n}=(n_0,n_1,..,n_{L1}):n_j\mathrm{I}\mathrm{N}\text{ and }\underset{j=0}{\overset{L1}{}}n_j=N\}.$$

The general state of the system can then be written as $`|\psi =_{\stackrel{}{n}\mathrm{\Lambda }^N}c(\stackrel{}{n})|\stackrel{}{n}`$.

We are interested in calculating the entanglement of one of the modes, say the mode 0, with respect to the rest of the system. In practice we can always reduce to a bipartite picture and we make use of the Von Neumann entropy to calculate the entanglement of the relevant mode. To do this we need to calculate the on-site reduced density matrix. It is interesting to note that $`H^{(N)}`$ cannot be written in terms of the tensor product of the Hilbert space of the mode subsystems: in order to calculate the on-site reduced density matrix one has to remind that the $`H^{(N)}`$ is a subspace of the entire Fock space, that by definition is written as a tensor product of the mode subsystems. Thus one can borrow this tensor product structure in order to calculate the trace over all the modes but the 0 one. Let be $`V=\{0,1,..,L1\}`$ the set of sites (modes, or vertices in the graphs terminology) then the reduced density matrix of the site 0 can be written as:

$`\rho ^{(0)}`$ $`=`$ $`\text{Tr}_{V_{\backslash 0}}(|\psi \psi |)`$
$`=`$ $`{\displaystyle \underset{\stackrel{}{n}\mathrm{\Lambda }_{\backslash 0}}{}}\stackrel{}{n}|\left({\displaystyle \underset{\stackrel{}{h},\stackrel{}{k}\mathrm{\Lambda }^N}{}}c(\stackrel{}{h})c^{}(\stackrel{}{k})|\stackrel{}{h}\stackrel{}{k}|\right)|\stackrel{}{n}`$
$`=`$ $`{\displaystyle \underset{j=0}{\overset{N}{}}}|jj|\left({\displaystyle \underset{\stackrel{}{n}(j)}{}}|c(\stackrel{}{n}(j))|^2\right)`$
$`=`$ $`{\displaystyle \underset{j=0}{\overset{N}{}}}\rho _j^{(0)}|jj|`$
where

$$\mathrm{\Lambda }_{\backslash 0}=\{\stackrel{}{n}=(n_1,..,n_{L1}):n_j\mathrm{I}\mathrm{N}\},$$

and where $`\stackrel{}{n}(j)=(j,n_1,..,n_{L1})\mathrm{\Lambda }^N`$.

In view of the constrain given by the fixed number of particles $`\rho ^{(0)}`$ is diagonal and the Von Neumann entropy relative to the mode 0 can be simply written as:

$`e^{(0)}`$ $`=`$ $`S(\rho ^{(0)})=Tr(\rho ^{(0)}\mathrm{log}_2\rho ^{(0)})`$
$`=`$ $`{\displaystyle \underset{j=0}{\overset{N}{}}}\rho _j^{(0)}\mathrm{log}_2\rho _j^{(0)}.`$
That is, $`e^{(0)}`$ is the Shannon entropy of the probability distribution

$$\{\rho _0^{(0)},\rho _1^{(0)},\mathrm{}.,\rho _N^{(0)}\}.$$
(2)
In order to characterize the different systems under consideration, and that are described in detail in the next section , we use of the following functional

$$EP^{(0)}=[\underset{t[0,T]}{sup}(e^{(0)}(t))]/\mathrm{log}(N+1)$$

The expression $`sup_{t[0,T]}(e^{(0)}(t))`$ represents the maximum value of the entanglement of site 0 with respect to the rest of the system during the evolution of the system; this value has been nomalized with respect to $`\mathrm{log}(N+1)`$ that is the maximum value of the mode entanglement that can be reached by the system and that corresponds to the situation in which $`\rho ^{(0)}=(\frac{1}{N+1})I`$. We take the value of $`EP^{(0)}`$ as a measure of the entagling power which characterizes the pair $`(U_{\mathrm{\Gamma }_j},|\psi _{in})`$; what we want to compute is, for each $`\mathrm{\Gamma }_j`$, the entangling capability of the evolution operator $`U_{\mathrm{\Gamma }_j}`$ when it is applied for a given period of time to the initial state $`|\psi _{in}`$.

As far as the choice of the initial states is concerned, we decide to use two particular kind of states for both the case $`N=L=3`$ and $`N=L=4`$. The first kind of states are the uniformly occupied ones : $`|111,|1111`$. The second kind of states are the ones in which all the particles are localized in the mode $`0`$: $`|300,|4000`$. This states allow us to test the entangling power of the operators $`U_{\mathrm{\Gamma }_j}`$ in two opposite and in a sense ”extremal” situations with respect to the localization of the particles in the sites.

This choice is also motivated by the possible physical implementations of the model under study. In fact one of these possible physical realizations is given by a set of ultracold bosons placed in an optical potential. The goal of placing exactly one atom in each well of an optical lattice created by a set of counterpropagating laser beams has been recently achieved bec\_mott  thanks to the realization of a super-fluid/Mott-insulator phase transition. This has motivated us to choose, among all the possible initial states of the evolution, the ones for which the site occupation number is equal to one, $`(n_j=1,j=0,..,L1)`$. In this kind of application, in order to realize the actual systems we are interested in, one should be able to isolate the relevant sites with the required graph topology.

The choice of the states in which all the particles are localized in one single mode is motivated by the fact that the typical realization of a BEC consists in the localization of all the bosons in a trapping (parabolic) potential; with a further step it is possible to create a lattice, by means of properly sized counterpropagating laser beams. In this kind of realization one should be able to isolate the relevant sites with the required graph topology and to fill with the proper number of particles the desired mode (well, site).

In general, for the realization of our model, one should also be able to manipulate the control parameters of the model, that is to fix the self-interaction and hopping parameters to the desired values, that are the same for all the sites and all the links between them. Thus, althought the organization of the sites in graphs could be a difficult task, the realization of our model requires only a ”global” control of the parameters since it does not require to perform control operations on single sites and/or single links between sites.

## V RESULTS

Our goal is to characterize, for each $`(U_{\mathrm{\Gamma }_j},|\psi _{in})`$, the behaviour of the entangling power $`EP_{\mathrm{\Gamma }_j}^{(0)}(\epsilon ,\tau )`$ when the ratio between the self-interaction parameter and the hopping parameter $`\epsilon /\tau `$ varies from zero to a value that is much greater then one. In order to do so we fix $`\tau =1`$ (that is we measure $`\epsilon `$ in $`\tau `$ units) and let $`\epsilon `$ vary in $`[\epsilon _{max},\epsilon _{max}]`$ ($`\epsilon _{max}>0`$) with step $`\mathrm{}\epsilon `$. For each $`\epsilon [\epsilon _{max},\epsilon _{max}]`$ we simulate the full quantum evolution for $`t[0,T]`$ that the system undegoes according to the unitary operator $`U_{\mathrm{\Gamma }_j}(t)=\mathrm{exp}(tH_{\mathrm{\Gamma }_j}(\epsilon ,\tau ))`$, $`(\mathrm{}=1)`$, when the initial state is $`|\psi _{in}`$. The evolution time $`T`$ has to be taken much greater than all the eigen-periods of the syetems i.e., $`T\text{max}_iE_i^1,`$ where $`E_i`$ denote an energy eigenvalue. If this latter condition holds one has no changes in the entangling power by allowing a longer evolution time. We explicitly verified this asymptotic behaviour in our numerical simulations. For each of these evolutions we calculate the value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ .

In this way we can consider three different regimes. When $`\epsilon =0`$ there is no self-interaction between the bosonic particles but these are free to hop between the different sites according to the graph $`\mathrm{\Gamma }`$. In this regime the intial state $`|\psi _{in}`$ is not an eigenstate of the pure hopping Hamiltonian. When $`\epsilon \epsilon _{max}`$ we have $`\epsilon /\tau 1`$, the hopping between the different sites is highly reduced and the systems tend to evolve toward states which are very close to the initial ones $`|\psi _{in}=|\stackrel{}{n}`$, that in this case are eigenstates of the pure self-interaction Hamiltonians corresponding to the eigenvalues $`\epsilon _{i=0}^{L1}n_i^2`$. The last regime is the intermediate one in which there is a competition of the two processes that govern the evolution, self-interaction and hopping.

To simulate the evolution of the systems under study we use a second order Runge-Kutta routine to integrate the time dependent Shrodinger equation for the different Hamiltonians.

### V.1 N=L=3

In this subsection we illustrate the results of our simulations for the case $`N=L=3`$. In the following we refer to the figures 3, 4 and 5. For each graph $`\mathrm{\Gamma }_j,j\{0,1,2\}`$ we have computed $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ as a function of $`\epsilon `$, using $`\mathrm{}\epsilon =0.1`$. For each value of $`\epsilon `$ we evolved the system from $`t=0`$ to $`T=15`$ (arbitrary units). For the curves $`(i)`$ and the $`(iii)`$ the intial state is $`|\psi _{in}=|111`$. The curves $`(ii)`$ and the $`(iv)`$ refer to simulations in which the initial state is $`|\psi _{in}=|300`$. In the $`(i)`$ and $`(ii)`$ cases $`\epsilon `$ varies in the interval $`[0,\epsilon _{max}]`$, while in the $`(iii)`$ and $`(iv)`$ cases $`\epsilon [\epsilon _{max},0]`$; for all the different cases $`\epsilon _{max}=+20`$ and the value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ has been plotted in function of $`|\epsilon |`$.

A first interesting feature about the systems under study is that when $`\epsilon =0`$ the entagling power of the evolution operator $`U_{\mathrm{\Gamma }_j}`$ does not depend on the graph $`\mathrm{\Gamma }_j`$ but only on the initial state. In fact, $`EP_{\mathrm{\Gamma }_j}^{(0)}0.92,j`$, in the case $`|\psi _{in}=|111`$ while $`EP_{\mathrm{\Gamma }_j}^{(0)}0.91,j`$, for $`|\psi _{in}=|300`$. The entagling capability of the $`U_{\mathrm{\Gamma }_j}`$’s , when $`\epsilon =0`$, is then only slightly enhanced if the system starts in a situation in which the particles are uniformly distributed over the sites.

We now turn to analize the maximum values reached by $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ for the different sistems. A first fact that emerges from our simulations is that the maximum theoretical value of the mode entanglement $`(\mathrm{log}_24)`$ can be reached for any $`j`$ if the system starts from $`|\psi _{in}=|300`$. In fact, for an appropriate choice of the ratio $`|\epsilon |/\tau `$, regardless of the sign of the self-interaction, $`EP_{\mathrm{\Gamma }_j}^{(0)}(|300,\epsilon )>0.99`$. If $`|\psi _{in}=|111`$ this is true only in the case $`j=0`$. Another interesting feature, common to the cases just mentioned, is that $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ increases very rapidly when $`|\epsilon |`$ starts to differ from $`0`$ and passes from $`0.910.92`$ to values $`>0.97`$ in correspondance of $`|\epsilon |=0.1`$; then the value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ remains very high $`(>0.95)`$ for $`|\epsilon |<\epsilon ^{}`$, where $`\epsilon ^{}>0`$ differs for the various cases. We can the say that, in this interval of values and for the above mentioned choice of $`(\mathrm{\Gamma }_j,|\psi _{in})`$, the self-interaction between the particles, either repulsive or attractive, enhances the capability of the $`U_{\mathrm{\Gamma }_j}`$’s to entangle the mode $`0`$ with the rest of the system.

On the contrary, in the cases $`(\mathrm{\Gamma }_1,|111)`$ and $`(\mathrm{\Gamma }_2,|111)`$ the maximum theoretical value of the mode entanglement is never reached and turning on the self-interaction only slightly enhance the entangling capability of the $`U_{\mathrm{\Gamma }_j}`$’s; this happens in correspondance of values of $`|\epsilon |`$ which are smaller or of the same order of $`\tau (=1)`$.

Another interesting feature highlighted by our simulations is related to the change of sign of the self-interaction parameter $`(\epsilon \epsilon )`$ and a particular property of the graphs $`\mathrm{\Gamma }_j`$’s. If we compare the simulations $`i`$’s vs $`iii`$’s we can see that for the graphs $`\mathrm{\Gamma }_0`$ and $`\mathrm{\Gamma }_2`$

$$EP_{\mathrm{\Gamma }_j}^{(0)}(\epsilon )=EP_{\mathrm{\Gamma }_j}^{(0)}(\epsilon )\epsilon [0,20]$$

while for $`\mathrm{\Gamma }_1`$ the curves $`i`$ and $`iii`$ differ significantly; the same holds for the curves $`ii`$ and $`vi`$. The difference of the behaviour of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ in the two cases depends on a particular property of the graphs we have used to characterize our systems. These can be distinguished in bi-partite and non bi-partite graphs. A graph $`\mathrm{\Gamma }`$ is bi-partite if

$$V=AB\text{ and }(a,b)EaA\text{ and }bB.$$

Suppose for example to colour each of the vertices of a graph with one of two different colours $`A`$ and $`B`$. A graph is bi-partite if the vertices coloured with $`A`$ are linked only with vertices coloured with $`B`$, and viceversa. In our case the bi-partite graphs are $`\mathrm{\Gamma }_0`$ and $`\mathrm{\Gamma }_2`$ (note that the bi-partiteness of a graph depends only on the pair $`(V,E)`$ and is independent from the root).

If $`\mathrm{\Gamma }_j`$ is bi-partite, the Hamiltonian that results from the change of sign of $`\epsilon `$ can be written as

$$H_{\mathrm{\Gamma }_j}(\epsilon ,\tau )=PH_{\mathrm{\Gamma }_j}(\epsilon ,\tau )P^{},$$
(3)
where $`P=\mathrm{exp}(i\pi _{jA}n_j)`$ is the unitary operator that implements the canonical transformation

$$c_j(1)^{\chi _A(j)}c_j$$
(4)
($`\chi _A`$ is the characteristic function on the subgraph $`A`$). The evolution operator for a given bi-partite graph can then be expressed as

$`U_{\mathrm{\Gamma }_j}^\epsilon (t)`$ $`=`$ $`\mathrm{exp}(itH_{\mathrm{\Gamma }_j}(\epsilon ,\tau ))=P\mathrm{exp}(itH_{\mathrm{\Gamma }_j}(\epsilon ,\tau ))P^{}`$
$`=`$ $`PU_{\mathrm{\Gamma }_j}^\epsilon (t)P^{}`$
If we now apply this operator to a state belonging to the occupation number basis $`|\stackrel{}{n}`$ we obtain:

$`|\psi _\stackrel{}{n}^\epsilon (t)`$ $`=`$ $`U_{\mathrm{\Gamma }_j}^\epsilon (t)|\stackrel{}{n}=P\mathrm{exp}(itH_{\mathrm{\Gamma }_j}(\epsilon ,\tau ))\mathrm{exp}(i\varphi _\stackrel{}{n})|\stackrel{}{n}`$
$`=`$ $`{\displaystyle \underset{\stackrel{}{m}\mathrm{\Lambda }^N}{}}\mathrm{exp}(i\varphi _\stackrel{}{m})\alpha _{\stackrel{}{n}\stackrel{}{m}}(t)\mathrm{exp}(i\varphi _\stackrel{}{n})|\stackrel{}{m}`$
where $`\mathrm{exp}(i\varphi _\stackrel{}{n})`$ and $`\mathrm{exp}(i\varphi _\stackrel{}{m})`$ are, respectively, the phase factors given by the action of the operator $`P^{}`$ on the initial state $`|\stackrel{}{n}`$ and of the operator $`P`$ on the states $`|\stackrel{}{m}`$ in the sum. The terms $`\alpha _{\stackrel{}{n}\stackrel{}{m}}(t)`$ can be obtained by acting with the time reversed operator $`U_{\mathrm{\Gamma }_j}^\epsilon (t)`$ on $`P^{}|\stackrel{}{n}`$ when it is written in terms of the energy basis states, in which the Hamiltonian is diagonal, and then going back to the site occupation number basis. Since $`j`$ $`H_{\mathrm{\Gamma }_j}`$ is a real simmetric matrix there exist a real orthogonal matrix $`O`$ ($`OO^T=1`$) that diagonalize it and with which it is possible to implement the required changes of basis. We can then write

$$\alpha _{\stackrel{}{n}\stackrel{}{m}}(t)=\underset{k}{}\beta _{\stackrel{}{n}\stackrel{}{m}k}\mathrm{exp}(+itE_k),\stackrel{}{n},\stackrel{}{m}$$

where $`E_k`$ are the eigenvalues of the Hamiltonian and $`\beta _{\stackrel{}{n}\stackrel{}{m}k}`$ can be written in terms of the real entries of $`O`$. This fact is crucial when we compare the evolutions of $`|\stackrel{}{n}`$ obtained by applying $`U_{\mathrm{\Gamma }_j}^\epsilon (t)`$ and $`U_{\mathrm{\Gamma }_j}^\epsilon (t)`$ (note the sign of $`t`$). The latter gives:

$`|\psi _\stackrel{}{n}^\epsilon (t)`$ $`=`$ $`\mathrm{exp}(itH_{\mathrm{\Gamma }_j}(\epsilon ,\tau ))|\stackrel{}{n}`$
$`=`$ $`{\displaystyle \underset{\stackrel{}{m}\mathrm{\Lambda }^N}{}}\alpha _{\stackrel{}{n}\stackrel{}{m}}(t)|\stackrel{}{m}`$
where

$$\alpha _{\stackrel{}{n}\stackrel{}{m}}(t)=\underset{k}{}\beta _{\stackrel{}{n}\stackrel{}{m}k}\mathrm{exp}(itE_k),\stackrel{}{n},\stackrel{}{m}.$$

We have now the elements to explain why, when $`\mathrm{\Gamma }_j`$ is bipartite, the curves $`(i)`$ and $`(ii)`$ overlap perfectly with the curves $`iii)`$ and $`vi)`$ respectively. Since the terms $`\beta _{\stackrel{}{n}\stackrel{}{m}k}`$ are real we have that for any $`t`$ and $`\stackrel{}{n},\stackrel{}{m}`$ the amplitudes of the states $`|\stackrel{}{m}`$ in the expression of $`|\psi _\stackrel{}{n}^\epsilon (t)`$ and $`|\psi _\stackrel{}{n}^\epsilon (t)`$ satisfy:

$$|\mathrm{exp}(i\varphi _\stackrel{}{m})\alpha _{\stackrel{}{n}\stackrel{}{m}}(t)\mathrm{exp}(i\varphi _\stackrel{}{n})|^2=|\alpha _{\stackrel{}{n}\stackrel{}{m}}(t)|^2,$$
(5)
in fact $`\alpha _{\stackrel{}{n}\stackrel{}{m}}(t)=\alpha _{\stackrel{}{n}\stackrel{}{m}}^{}(t)`$. Consequently, since the value of the mode entanglement and of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ depends only on the modulus square of the probability amplitudes of the basis states $`|\stackrel{}{m}`$ in the expression of $`|\psi _\stackrel{}{n}^{\pm \epsilon }(t)`$ (see section IV), for any given $`|\psi _{in}=|\stackrel{}{n}`$ and for all the bi-partite graphs we have $`EP_{\mathrm{\Gamma }_j}^{(0)}(\epsilon )=EP_{\mathrm{\Gamma }_j}^{(0)}(\epsilon )`$. In general, this is no longer true for the non-bipartite graphs.

Another relevant and general feature of the entagling power $`\mathrm{\Gamma }`$’s and for both the selected initial state is that $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ rapidly decreses when $`|\epsilon |`$ is significantly greater than $`1(=\tau )`$ with a rate that is higer in the case $`|\psi _{in}=|300`$. This can be explained with the following argument. When $`|\epsilon |`$ becomes significantly greater then 1, the Hamiltonian of the system, $`\mathrm{\Gamma }_j`$, can be considered as composed of an unperturbed diagonal term $`H_\epsilon `$, the self-interaction one, and a perturbation term given by the hopping matrix $`H_\tau `$. Both the initial states of our computation are eigenvectors of the unperturbed term and correspond to the eigenvalues $`3\epsilon `$ ($`|111`$) and $`9\epsilon `$ ($`|300`$). As $`|\epsilon |`$ increses, the perturbation becomes smaller and the only processes relevant for the evolution are the first order ones. This means that the probability distribution $`(\rho _0^{(0)},\rho _1^{(0)},\rho _2^{(0)},\rho _3^{(0)})`$, given by the diagonal elements of $`\rho ^0`$, is higly unbalanced and this results in a value of $`e^{(0)}(t)`$ that for any $`t`$ is very small, and so it is, consequently, the value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$.

To give a more detailed explanation of this phenomenon let’s take for example $`\mathrm{\Gamma }_0`$ and suppose $`|\psi _{in}=|111`$. In general, at first order, the state $`|\psi _{in}`$ can be coupled by the perturbation to the states $`|\stackrel{}{n}^{}\chi _{|\psi _{in}}^{}`$, where

$$\chi _{|\psi _{in}}^{}=\{|\stackrel{}{n}^{}|\stackrel{}{n}^{}\mathrm{\Lambda }^N\text{and}\stackrel{}{n}^{}|H_\tau |\psi _{in}0\}.$$

In our case $`\chi _{|111}^{}=\{|021,|012,|201,|210\}`$: the modulus square of the probability amplitudes of these states contribute to the terms $`\rho _0^{(0)},\rho _2^{(0)}`$ of the distribution. The state in which $`n_0=3`$ is coupled to $`|111`$ only by second order processes and for this reason the term $`\rho _3^{(0)}\rho _0^{(0)},\rho _1^{(0)},\rho _2^{(0)}`$; furthermore, when the perturbation is very small, the terms $`\rho _0^{(0)},\rho _2^{(0)}\rho _1^{(0)}`$ and the distribution is very peaked around $`\rho _1^{(0)}`$.

The same arguments hold for $`\mathrm{\Gamma }_0`$ when $`|\psi _{in}=|300`$, but there is a significant difference between the two cases. In fact $`\chi _{|300}^{}=\{|201,|210\}`$ so that when the pertubation is small $`\rho _0^{(0)},\rho _1^{(0)}\rho _2^{(0)}\rho _3^{(0)}`$: the probability distribution is more unbalanced than in the previous case and this gives a first explanation of the fact that in the perturbative regime

$$EP_{\mathrm{\Gamma }_j}^{(0)}(|111)>EP_{\mathrm{\Gamma }_j}^{(0)}(|300),j.$$

There is another important difference between the two cases that justifies this relation: the strenght of the coupling. This, at first order, depends on the ratio $`\tau /\mathrm{}E`$, where $`\mathrm{}E`$ is the difference in energy between $`|\psi _{in}`$ and the states belonging to $`\chi _{|\psi _{in}}^{}`$. The ratio is higher for $`|\psi _{in}=|111(1/2|\epsilon |)`$ then for $`|\psi _{in}=|300(1/4|\epsilon |)`$ and this results in higer values of probability relative to the states in $`\chi _{|\psi _{in}}^{}`$ in the first case. Consequently the associated eigenvalues of $`\rho ^{(0)}`$ can be higher for $`|\psi _{in}=|111`$ thus giving rise to a probability distribution which is flatter and to higher values of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$.

We now use the perturbative approach to describe the features relative to the non-bipartite graph $`\mathrm{\Gamma }_1`$. Here the interesting thing to notice is again the role of the sign of the self-interaction, which depends on the initial state. For $`|\psi _{in}=|300`$ a repulsive (positive) self-interaction enhances the entangling capability of $`U_{\mathrm{\Gamma }_j}(t)`$. In fact, for a large part of the interval $`[0,|\epsilon _{max}|]`$,

$$EP_{\mathrm{\Gamma }_1}^{(0)}(|300,+|\epsilon |)EP_{\mathrm{\Gamma }_1}^{(0)}(|300,|\epsilon |).$$
(6)
On the other hand, when $`|\psi _{in}=|111`$ is the attractive (negative) self-interaction that favors the entangling power of the evolution operator. In fact for any $`|\epsilon |[0.5,\epsilon _{max}]`$

$$EP_{\mathrm{\Gamma }_1}^{(0)}(|111,|\epsilon |)>EP_{\mathrm{\Gamma }_1}^{(0)}(|111,+|\epsilon |).$$
(7)
This behaviour can be understood in detail at least in a perturbative picture. Let’s start from the case $`|\psi _{in}=|111`$. When $`|\epsilon |`$ is significanly greater than $`1(=\tau )`$ we can use a two-level picture to describe our system. The first level is given by $`|\psi _{in}`$, while the second one is given a uniform superposition of the $`q=|\chi _{|\psi _{in}}^{}|`$ degenerate states that are coupled to $`|\psi _{in}`$ by first order processes, that is

$$|I=(1/\sqrt{q})\underset{|\stackrel{}{n}^{}\chi _{|\psi _{in}}^{}}{}|\stackrel{}{n}^{}H_\tau |\psi _{in}.$$
(8)
The two-levels Hamiltonian can be written in general as

$$H=\left(\begin{array}{ccc}E_1& W& \\ W& E_2& \end{array}\right)$$
(9)
where

$`E_1`$ $`=`$ $`\psi _{in}|H_\epsilon |\psi _{in},`$
$`W`$ $`=`$ $`I|H_\tau |\psi _{in}=\psi _{in}|H_\tau |I^{},`$
$`E_2`$ $`=`$ $`I|H_\epsilon |I+I|H_\tau |I.`$
It is a standard result of two coupled levels systems theory that the probability of finding the system in the state $`|I`$ at time $`t`$ can be written as (Rabi’s formula)

$`P_{|\psi _{in}|I}(t)`$ $`=`$ $`|I|\mathrm{exp}[it(H_\epsilon +H_\tau )|\psi _{in}|^2`$
$`=`$ $`{\displaystyle \frac{4|W|^2}{4|W|^2+(E_1E_2)^2}}\mathrm{sin}^2(\theta t),`$
where $`\theta =\sqrt{4|W|^2+(E_1E_2)^2}/2`$. Here the important fact is the role of the term $`(E_1E_2)^2`$ in the denominator from which depends the maximum value attainable by $`P_{|\psi _{in}|I}(t)`$; this is directly linked with the maximum values of $`\rho _0^{(0)},\rho _2^{(0)}`$ and consequently with the maximum value attainable by $`EP_{\mathrm{\Gamma }_1}^{(0)}`$ in the perturbative regime. In fact it is possible to reproduce, with a good approximation, the results of our simulations by using this probability. This can be done by calculating the mode entanglement of the following state:

$$|\varphi =\sqrt{1P}|\psi _{in}+\sqrt{P}|I$$
(10)
where

$$P=P(\epsilon ,\tau )=\frac{4|W|^2}{4|W|^2+(E_1E_2)^2}$$
(11)
is the maximum attainable value of $`P_{|\psi _{in}|I}(t)`$ according to the Rabi’s formula. In the case of the system $`\mathrm{\Gamma }_0`$, the state $`|I`$ is given by a superposition of states for which $`n_0=0,2`$. The only contribution to the terms $`\rho _0^{(0)},\rho _2^{(0)}`$ of the reduced density matrix of the mode $`0`$ are given by these states. This means, as described in the previous paragraphs, that the higher are the values of these terms the ”flatter” is the probability distribution $`(\rho _0^{(0)},\rho _1^{(0)},\rho _2^{(0)},\rho _3^{(0)},)`$ and consequently the higher is the value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$. For this reason the state $`|\varphi `$ can be considered as an example of state in which the entanglement is maximum, and we can compute the entangling power for any rooted graph, and for all the values of $`\epsilon `$ in the perturbative regime as $`EP_{\mathrm{\Gamma }_j}^{(0)}=Tr(\rho ^{(0)}\mathrm{log}_2\rho ^{(0)})`$ where $`\rho ^{(0)}=\text{Tr}_{V_{\backslash 0}}(|\varphi \varphi |)`$.

By means of the two levels picture we now have the opportunity to see from another perspective the difference between bipartite and non-bipartite graphs. This result is essential to explain the difference in the behaviour of $`EP_{\mathrm{\Gamma }_1}^{(0)}`$ described by equations (6) and (7). Since both $`|\psi _{in}`$ and $`|I`$ are eigenvectors of $`H_\epsilon `$, if $`I|H_\tau |I=0`$, the difference $`(E_1E_2)`$ does not depend on $`\tau `$ and is linear in $`\epsilon `$; thus the sign of the latter does not have any effect on $`P`$. This is always true for all bipartite graphs.To demonstrate that for all bipartite graphs $`I|H_\tau |I=0`$, we first observe that, since $`|IH_\tau |\psi _{in}`$, we can always write

$$I|H_\tau |I\psi _{in}|H_\tau ^3|\psi _{in}=111|H_\tau ^3|111.$$
(12)
We can then make use of the transformation $`P=\mathrm{exp}(\pi _{jA}n_j)`$ previously described to demonstrate that for any bipartite graph $`111|H_\tau ^3|111=0`$. In order to do this we calculate the expectation value $`111|P^{}H_\tau ^3P|111`$ in two equivalent ways; in the first we let act $`P(P^{})`$ on $`|111(111|)`$ to obtain:

$$(111|P^{})H_\tau ^3(P|111)=111|H_\tau ^3|111.$$

In the second one we use the fact that for any bipartite graph $`P^{}H_\tau P=H_\tau `$; this allows to write

$$111|P^{}H_\tau ^3P|111=111|H_\tau ^3|111.$$

This implies that for any bipartite graph

$$I|H_\tau |I=111|H_\tau ^3|111=0.$$
(13)
The last result is, in general, no longer true when the graph is non-bipartite. In fact, for $`\mathrm{\Gamma }_1`$, $`I|H_\tau |I`$ is constant and proportional to $`\tau `$. This is essential to explain the difference in the behaviour of $`EP_{\mathrm{\Gamma }_1}^{(0)}`$ in correspondance of an attractive or a repulsive self-interaction. In fact, the presence in $`E_2`$ of the term proportional to $`\tau `$ renders the difference $`(E_1E_2)`$ sentive to the sign of $`\epsilon `$.

$$(E_1E_2)=(2\epsilon +3\tau ).$$

This implies that

$$P(|\epsilon |)>P(+|\epsilon |),$$

and consequently, as one can see by calculating $`EP_{\mathrm{\Gamma }_1}^{(0)}`$ from the state $`|\varphi `$,

$$EP_{\mathrm{\Gamma }_1}^{(0)}(|111,|\epsilon |)>EP_{\mathrm{\Gamma }_1}^{(0)}(|111,+|\epsilon |).$$

Since the resuts $`(12)`$ and its consequence $`(13)`$ hold also for $`|\psi _{in}=|300`$ we can use the same arguments to give a justification of $`(1)`$. Here simple calculations give:

$$E_1E_2=4\epsilon \tau .$$

The role of $`\epsilon `$ is then opposite with respect to the previous case:

$$P(+|\epsilon |)>P(|\epsilon |),$$

and consequently, in the perturbative limit

$$EP_{\mathrm{\Gamma }_1}^{(0)}(|300,+|\epsilon |)EP_{\mathrm{\Gamma }_1}^{(0)}(|300,|\epsilon |).$$

### V.2 N=L=4

In this subsection we illustrate the results of our simulations for the case $`N=L=4`$. In the following we refer to the figures 6-16. For each graph $`\mathrm{\Gamma }_j,j\{3,4,..,13\}`$ we have computed $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ as a function of $`\epsilon `$, using $`\mathrm{}\epsilon =0.1`$. For each value of $`\epsilon `$ we evolved the system from $`t=0`$ to $`T=15`$ (arbitrary units). The curves $`(i)`$ and $`(iii)`$ in each graph are the results of the simulations for $`|\psi _{in}=|1111`$. The curves $`(ii)`$ and $`(iv)`$ refer to the simulations in which the initial state is $`|\psi _{in}=|4000`$. In the $`(i)`$ and $`(ii)`$ curves $`\epsilon `$ varies in the interval $`[0,\epsilon _{max}]`$, while in the $`(iii)`$ and $`(iv)`$ curves $`\epsilon [\epsilon _{max},0]`$; for all the different cases $`\epsilon _{max}=+20`$ and the value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ has been plotted in function of $`|\epsilon |`$.

As we did in the previous section we start by observing that when $`\epsilon =0`$ the entagling power of the evolution operator $`U_{\mathrm{\Gamma }_j}`$ does not depend on the graph $`\mathrm{\Gamma }_j`$ but only on the chosen initial state. In fact $`EP_{\mathrm{\Gamma }_j}^{(0)}0.83\text{, for any }j`$, in the case $`|\psi _{in}=|1111`$, while $`EP_{\mathrm{\Gamma }_j}^{(0)}0.87\text{, for any }j`$, when $`|\psi _{in}=|4000`$. We see here a first difference with case $`N=L=3`$: the entagling capability of the $`U_{\mathrm{\Gamma }_j}`$’s , when $`\epsilon =0`$, is enhanced if the system starts with all particles localized in the mode $`0`$.

We now turn to analize the maximum values reached by $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ for the different sistems. The maximum theoretical value of the mode entanglement $`(\mathrm{log}_25)`$ is always reached $`\mathrm{\Gamma }_j`$ when the initial state is $`|4000`$. In fact, for an appropriate choice of the ratio $`|\epsilon |/\tau `$, regardless of the sign of the self-interaction, $`EP_{\mathrm{\Gamma }_j}^{(0)}(|4000)>0.99`$ . In contrast with the case $`N=L=3`$, none of the systems reaches this value when the initial state is the uniformly occupied one; the highest values of $`EP_{\mathrm{\Gamma }_j}^{(0)}(|1111,)`$ are $`0.97`$ and they are reached by the systems $`\mathrm{\Gamma }_5,\mathrm{\Gamma }_{11}`$ ( only for $`\epsilon >0`$), and $`\mathrm{\Gamma }_{10}`$.

When $`|\psi _{in}=|4000`$, $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ increases when $`|\epsilon |`$ starts to differ from $`0`$ and passes from $`0.87`$ to values $`>0.92`$ in correspondance of $`|\epsilon |=0.1`$; then the value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ remains very high for $`|\epsilon |<\epsilon ^{}`$, where $`\epsilon ^{}>0`$ differs for the various cases. We can then say that, in this interval of values the self-interaction between the particles, either repulsive or attractive, enhances the capability of the $`U_{\mathrm{\Gamma }_j}`$’s to entangle the mode $`0`$ with the rest of the system.

In the case $`|\psi _{in}=|1111`$, as mentioned above, the maximum theoretical value of the mode entanglement is never reached and the effect of turning on the self-interaction depends on the system and on the sign of the self-interaction. However, the enhancement of the entangling capability of the $`U_{\mathrm{\Gamma }_j}`$’s, when present, is found in correspondance of values of $`|\epsilon |`$ which are smaller or of the same order of $`\tau (=1)`$.

By observing the figures 6-16 we can see that also for the case $`N=L=4`$ we can distinguish between two different behaviours of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ with respect to the sign change of $`\epsilon `$. For $`j\{3,4,6,8,10\}`$ we have that for both the intial states

$$EP_{\mathrm{\Gamma }_j}^{(0)}(\epsilon )=EP_{\mathrm{\Gamma }_j}^{(0)}(\epsilon )\epsilon [0,20],$$
(14)
that is the curves $`i`$ and $`ii`$ perfectly overlap with the curves $`iii`$ and $`iv`$ respectively. On the other hand for $`j\{5,7,9,11,12,13\}`$ this is no longer true.

This fact can be explained using the same arguments we used in the previous section. In fact the $`\mathrm{\Gamma }_j,j=3,\mathrm{..13}`$, can be distinguished in bi-partite and non-bipartite graphs. For the bipartite ones we can apply to the related Hamiltonians the unitary transformation $`(4)`$ in order to get equation $`(3)`$. The latter and the subsequent reasonings hold for any $`N=L`$ and thus equation $`(5)`$ implies that $`j\{3,4,6,8,10\}`$, i.e. for all bipartite graphs equation $`(14)`$ must hold.

We now turn to examine the behaviour of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ in the pertubative regime. We first observe that for both the initial states $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ rapidly decreses when $`|\epsilon |`$ is significantly greater than $`1(=\tau )`$. This happens for the same reasons we described for the case $`N=L=3`$. The Hamiltonian, $`\mathrm{\Gamma }_j`$, is composed of an unperturbed diagonal term $`H_\epsilon `$, the self-interaction one, and a perturbation term given by the hopping matrix $`H_\tau `$. Both the initial states of our computation are eigenvectors of the unperturbed term and correspond to the eigenvalues $`4\epsilon `$ ($`|1111`$) and $`16\epsilon `$ ($`|4000`$). As $`|\epsilon |`$ increses, the perturbation becomes small and the only processes that become relevant for the evolution are the first order ones. At first order, the state $`|\psi _{in}`$ can be coupled by the perturbation to the states $`|\stackrel{}{n}^{}\chi _{|\psi _{in}}^{}`$. This implies that for $`|\psi _{in}=|1111`$ we have that $`\rho _3^{(0)},\rho _4^{(0)}\rho _0^{(0)},\rho _2^{(0)}\rho _1^{(0)}`$, while for $`|\psi _{in}=|4000`$ we have that $`\rho _0^{(0)},\rho _1^{(0)},\rho _2^{(0)}\rho _3^{(0)}\rho _4^{(0)}`$. The probability distribution $`(\rho _0^{(0)},\rho _1^{(0)},\rho _2^{(0)},\rho _3^{(0)},\rho _4^{(0)})`$ is higly unbalanced in both cases, and this results in the decreasing behaviour of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$, but it is more unbalanced in the case $`|\psi _{in}=|4000,j,`$ and this is a first explanation of the fact that in the perturbative regime

$$EP_{\mathrm{\Gamma }_j}^{(0)}(|1111)>EP_{\mathrm{\Gamma }_j}^{(0)}(|4000),j.$$
(15)
A further explanation for the last relation is given, just as for the case $`N=L=3`$, by the fact that the strenght of the first order coupling between $`|\psi _{in}`$ and the states belonging to $`\chi _{|\psi _{in}}^{}`$ is higher in the case $`|\psi _{in}=|1111,(\tau /\mathrm{}E=1/(2|\epsilon |),`$ then in the case $`|\psi _{in}=|4000,(\tau /\mathrm{}E=1/(6|\epsilon |)`$. As explained in the previous section this affects the highest reachable value for the probability related to the states belonging to in $`\chi _{|\psi _{in}}^{}`$ and consequently affects the values reached by $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ in the way described by $`(15)`$.

We finally come to the description of the role of the sign of the self-interaction when $`|\epsilon |`$ is significantly greater than $`1(=\tau )`$ in the case of the non-bipartite graphs $`\mathrm{\Gamma }_j,j\{5,7,9,11,12,13\}`$. The phenomenology is analogous to the one found for the case $`N=L=3`$. For $`|\psi _{in}=|4000`$ a repulsive (positive) self-interaction enhances the entangling capability of $`U_{\mathrm{\Gamma }_j}(t)`$, that is for a large part of the interval $`[0,|\epsilon _{max}|]`$

$$EP_{\mathrm{\Gamma }_j}^{(0)}(|4000,+|\epsilon |)EP_{\mathrm{\Gamma }_j}^{(0)}(|4000,|\epsilon |).$$
(16)
In contrast, for $`|\psi _{in}=|1111`$ it is the attractive (negative) self-interaction that for a large part of the interval $`[0,|\epsilon _{max}|]`$ enhances the entangling power of $`U_{\mathrm{\Gamma }_j}(t)`$:

$$EP_{\mathrm{\Gamma }_j}^{(0)}(|1111,|\epsilon |)>EP_{\mathrm{\Gamma }_j}^{(0)}(|1111,+|\epsilon |).$$
(17)
This behaviour can be understood by means of the same arguments exposed for the case $`N=L=3`$ that we briefly recall here. When $`|\epsilon |`$ is significantly greater than $`1(=\tau )`$, we can again use a two-levels picture to describe our systems. The first level is $`|\psi _{in}`$, while the second $`|I`$ is given by a uniform superposition of the $`q=|\chi _{|\psi _{in}}^{}|`$ degenerate states that are coupled to $`|\psi _{in}`$ by first order processes, see eq. (8). The two-levels Hamiltonian is again given by (see formula) and the probability $`P_{|\psi _{in}|I}(t)`$ is expressed by the Rabi’s formula. As we have seen for the case $`N=L=3`$ with the introduction of the state $`|\varphi `$, the maximum attainable value of this probability $`P=P(\epsilon ,\tau )`$ is directly related to the value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$, that is the higher the former the higher the latter. In this context it is again the term $`(E_1E_2)`$ that is crucial to give an explanation of the equations (16) and (17). In fact, equation (13) holds also for the case $`N=L=4`$ and for the intial states $`|1111\text{ and }|4000`$. Thus, for the bipartite graphs, we have that $`I|H_\tau |I=0`$, the term $`(E_1E_2)`$ is independent of $`\tau `$ and depends linearly on $`\epsilon `$. This implies that the probability $`P_{|\psi _{in}|I}(t)`$ does not depend on the sign of $`\epsilon `$ for the bipartite graphs $`\mathrm{\Gamma }_j,j\{3,4,6,8,10\}`$.

In the case of non-bipartite graphs this is no longer true and the dependance of $`(E_1E_2)`$ on $`\tau `$ give rise to a sensitivity to the sign of $`\epsilon `$ of $`P=P(\epsilon ,\tau )`$. Hence the different behaviour of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ described by equations (16) and (17).

### V.3 Graph Properties and Entangling Power in the Pertubative Regime

In the last part of the paper we want to deepen our analysis in the perturbative regime in order to describe the existing relations among the behaviour of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ for the different $`\mathrm{\Gamma }_j`$. As we have seen in the previous paragraphs, in the perturbative regime, this behaviour can be described by using for our systems a two levels approach. In this picture, the maximum value attainable for $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ is directly linked with $`P`$, i.e. the maximum of the probability $`P_{|\psi _{in}|I}(t)`$, see equations (8)-(11). It is possible to reproduce, with a good approximation, the results of our simulations by using this probability and by calculating the mode entanglement of the state (10), that is:

$$EP_{\mathrm{\Gamma }_j}^{(0)}=Tr(\rho ^{(0)}\mathrm{log}_2\rho ^{(0)})$$
(18)
where $`\rho ^{(0)}=\text{Tr}_{V_{\backslash 0}}(|\varphi \varphi |)`$. In order to evaluate (18) we can use analytical expressions for $`P`$ and for $`\rho ^{(0)}`$ that allow us to link $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ to some of the topological properties of the graphs. These expressions can be derived for the general case $`N=L`$.

We first start our analysis by considering $`|\psi _{in}=_{j=0}^{L1}|1_j`$; in this case the expression for P is:

$$P=\frac{16\tau ^2k}{16\tau ^2k+(2\epsilon +\tau 18l_3/2k)^2},$$
(19)
where $`16\tau ^2k=4|\tau \sqrt{2q}|=4|W|^2`$ is given by the coupling induced by the hopping term $`H_\tau `$ and is proportional to $`k=|E|=q/2=|\chi _{|\psi _{in}}^{}|/2`$, i.e. the number of links of the graph; $`2\epsilon =I|H_\epsilon |I\psi _{in}|H_\epsilon |\psi _{in}`$. The term $`\tau ^218l_3/q=\tau ^218l_3/2k=I|H_\tau |I`$ depends on $`k`$ and also on another property of the graph which is the number $`l_3`$ of closed loops of length $`3`$ present in the graph (one of such loops can be seen for example in $`\mathrm{\Gamma }_{11}`$ where the loop is given by the links between the sites $`0,1`$ and $`2`$).

We now give a brief explanation of the appearance of the number $`l_3`$ in the expression in the expression of $`P`$. As we have already pointed out in section V.A, eq. (12), $`I|H_\tau |I`$ is proportional to $`\psi _{in}|H_\tau ^3|\psi _{in}`$. The fact that this term differs from zero only if the graph is non bi-partite can be now explained in terms of the presence of loops of length $`3`$ in the graph. In fact a graph is non bi-partite if and only if it contains at least one of these loops. The presence in the graph of one of such loops results in the presence in $`H_\tau ^3`$ of terms like $`T_{i,j,h,k,m,n}=c_i^{}c_jc_h^{}c_kc_m^{}c_n`$ such that

$$\psi _{in}|T_{i,j,h,k,m,n}|\psi _{in}0.$$
(20)
For example, in the case $`\mathrm{\Gamma }_{11}`$, $`c_0^{}c_2c_2^{}c_1c_1^{}c_0`$ is one of such terms. In order to calculate $`I|H_\tau |I`$ we first have to compute the number of terms for which equation (20) holds. This can be done, for any loop of length three, as follows. We first observe that the number of possible choices for the index $`n`$ is $`3`$, that is the number of vertices involved in the loop. Since each of these is linked to the other $`2`$ vertices in the loop we have $`2`$ possible choices for $`m`$. Once $`n`$ and $`m`$ are fixed there are only two kind of terms that satisfy equation (20), that is $`a)T_{n,h,h,m,m,n}`$ and $`b)T_{h,m,n,h,m,n}`$, where $`h`$ is the third remaining vertex. Consequently, for any loop of length three we have $`12`$ terms that satisfy equation (20). By inspection we see that for the six terms of type $`a)\psi _{in}|T_{n,h,h,m,m,n}|\psi _{in}=4`$, while for the other six $`\psi _{in}|T_{h,m,n,h,m,n}|\psi _{in}=2`$. This implies that, taking in account of the fixed contribute given by each of the loops of length three present in the graph, $`I|H_\tau |I=18l_3/q`$.

We have then succeded in linking $`P`$ to the graph properties $`k`$ and $`l_3`$.The next step is to explicitly write $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ as the entanglement of the state $`|\varphi `$; to do this it is necessary to express $`\rho _0^{(0)},\rho _1^{(0)}`$ and $`\rho _2^{(0)}`$ in terms of the probability $`P`$. This can be done in a way that takes in account of the properties of the graph $`\mathrm{\Gamma }_j`$ and of the state $`|I`$:

$$\rho _1^{(0)}=(1P)+\frac{P}{q}2(kk_0)=1P\frac{k_0}{k}$$
(21)
$$\rho _0^{(0)},\rho _2^{(0)}=\frac{P}{q}k_0=\frac{Pk_0}{2k},$$
(22)
where $`q=|\chi _{|\psi _{in}}^{}|=2k`$; $`k`$ is the number of links in the graph and $`k_0`$ is the number of links starting from the root vertex. The term $`\frac{P}{q}2(kk_0)`$ in (21) takes in account of the $`2(kk_0)`$ states in $`|I`$ that are coupled by first order processes to $`|\psi _{in}`$ and for which $`n_0=1`$. We can now use the formulas just derived to describe some of the relations existing among $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ for the different rooted graphs in the case $`|\psi _{in}=_{j=0}^{L1}|1_j`$.

We first focus on the relation among the values of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$, in the perturbative regime, for rooted graphs that differ only for the choice of the root vertex. As an example we take the rooted graphs $`\mathrm{\Gamma }_5,\mathrm{\Gamma }_7`$ and $`\mathrm{\Gamma }_{11}`$; in figure (17) we have plotted the results of our simulations for the three graphs in the case $`|\psi _{in}=|1111`$ and $`\epsilon >0`$. From this figure we can see that

$$EP_{\mathrm{\Gamma }_{11}}^{(0)}>EP_{\mathrm{\Gamma }_7}^{(0)}EP_{\mathrm{\Gamma }_5}^{(0)}.$$

That is, the higher the value of $`k_0^j`$, i.e. of the coordination number of the root vertex of $`\mathrm{\Gamma }_j`$, the higher the value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$. This is an example of a feature that holds in general for any value of $`N=L`$. In fact, suppose to have

* a set of rooted graph $`A=\{\mathrm{\Gamma }_j\}`$ that differ one from the other only for the choice of the root vertex;
* $`|\psi _{in}=_{j=0}^{L1}|1_j`$;
* a fixed value of $`\epsilon `$ ($`>0\text{ or }<0`$) in the perturbative regime,

then the $`\mathrm{\Gamma }_j`$’s for which $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ is higher are those for which $`k_0^j`$ is maximum. Furthermore it is possible to say that

$$\text{ if }k_0^m>k_0^n\text{ then }EP_{\mathrm{\Gamma }_m}^{(0)}>EP_{\mathrm{\Gamma }_n}^{(0)}.$$
(23)
This can be seen from equations (19),(21),(22). In fact, for a given value of $`\epsilon `$ in the perturbative regime, $`P`$ is the same for all the rooted graphs in $`A`$, since for all them $`k`$ and $`l_3`$ are the same. On the contrary, the flatness of the probability distribution (2) and consequently the value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ depend, in the way described above, on the value of $`k_0^j`$.

We now consider a different situation in which it is meaningful to compare the values of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ for different rooted graphs. We take as example the set of graphs $`B=\{\mathrm{\Gamma }_{10},\mathrm{\Gamma }_{11},\mathrm{\Gamma }_{12},\mathrm{\Gamma }_{13}\}.`$ It is interesting to notice that the non-bipartite graphs $`\mathrm{\Gamma }_{11},\mathrm{\Gamma }_{12}`$, and $`\mathrm{\Gamma }_{13}`$ can be seen as generated from the bipartite graph $`\mathrm{\Gamma }_{10}`$ by adding respectively one, two and three links, that is by increasing the connectivity of the subgraph $`E^{}`$. It is then interesting to see what happens to the maximum attainable value of the mode entanglement when the connectivity of the subgraph changes. We refer to the figures (18) and (19), which describe the cases $`\epsilon >0`$ and $`\epsilon <0`$ respectively ($`|\psi _{in}=|1111`$). In figure (18) we can see that, for a given $`\epsilon >0`$, the increase of the connectivity of the subgraph $`E^{}`$ has the effect of reducing the value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$. On the contrary, when $`\epsilon <0`$ the effect is opposite, see figure (19). What one could conclude by looking at equations (21)-(22) is that $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ should decrease when the connectivity of the subgraph increses, regardless of the sign of the self-interaction, in view of the dependance of $`\rho _1^{(0)}`$ and $`\rho _2^{(0)}`$ on $`k_0/k`$. But this contradicts the behaviour described in figure (19). In fact, it is the role of the sign of the self-interaction in $`P`$ that allows $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ to increase with the increase of the connectivity of the sub graph $`E^{}`$. We can see here the role of the interplay between the self-interaction and the tunnelling processes.

We can now end our analysis by considering the case $`N=L`$ and $`|\psi _{in}=|N_{j=0}^{L1}|0_j.`$ As a first step we can write the analitical expression of $`P`$ in terms of the graph properties as follows:

$$P=\frac{4\tau ^2k_0N}{4\tau ^2k_0N+[2(N1)\epsilon 2\tau \frac{l_3^0}{k_0}]^2},$$
(24)
where $`k_0`$ is the number of links starting from the root vertex, while $`l_3^0`$ is the number of loops of length three that include the root vertex. Since for all the states in $`|I`$ $`n_0=N1`$, for the determination of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ in the perturbative regime the only terms of the probability distribution (2) that are relevant are $`\rho _N^{(0)}`$ and $`\rho _{N1}^{(0)}`$. These can be simply written as:

$`\rho _N^{(0)}`$ $`=`$ $`1P`$
$`\rho _{N1}^{(0)}`$ $`=`$ $`P.`$
With the help of (24) we can now see that the result (23) hold also for the state $`|\psi _{in}=|N_{j=0}^{L1}|0_j.`$ In fact, here it is dependence of $`P`$ on $`k_0`$ that determines the behaviour described by (23). Furthermor, also in this case it is possible to say that the $`\mathrm{\Gamma }_j`$’s for which $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ is higher are those for which $`k_0^j`$ is maximum.

Finally we consider again the case $`B=\{\mathrm{\Gamma }_{10},\mathrm{\Gamma }_{11},\mathrm{\Gamma }_{12},\mathrm{\Gamma }_{13}\}`$, i.e. the rooted graphs that can be can be seen as generated from the bipartite $`\mathrm{\Gamma }_{10}`$ by increasing the connectivity of the subgraph $`E^{}`$. In this case the increase of the connectivity of the subgraph plays a role only if it implies the presence in the graph of loops of length three that include the root vertex. For the set $`B`$ we have that $`l_3^0=0,1,2,3`$ for $`\mathrm{\Gamma }_{10},\mathrm{\Gamma }_{11},\mathrm{\Gamma }_{12}`$, and $`\mathrm{\Gamma }_{13}`$ respectively. The interplay between the self-interaction and the tunnelling processes plays a role which is opposite with respect to the case $`|\psi _{in}=|1111`$. In fact, in the case $`|\psi _{in}=|4000`$ we have that for a given $`\epsilon >0`$, the increase of the connectivity of the subgraph $`E^{}`$ has the effect of increasing the value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$. On the contrary, when $`\epsilon <0`$ the effect is opposite.

## VI conclusions

We have analyzed by means of numerical simulations the mode-entanglement properties of of bosonic particles hopping over graph structures. Starting from an initial unentangled state bi-partite mode-entanglement of a graph vertex with respect the rest of the graph is generated, by allowing for a finite time a tunneling along the graph edges. We characterized the entangling power of a given hopping hamiltonian by the maximum achieved during the dynamical evolution by this entanglement. We studied the dependence of this entangling power by the local self-interaction parameters for all the graphs up to four vertices and for two different natural choices of the initial state. Our results show the important role of that graph topology has in optimization of entanglement production

The main steps and results,for both cases $`N=L=3,4`$ , can be summarized as follows.

* We have classified the way in which a mode (bosonic site) can be linked to the ”rest of the world” (other bosonic sites) in terms of rooted graphs which are inequivalent from the poin of view of the mode entanglement of the reference mode (root).
* When the self-interaction term is absent ($`=0`$) the value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ does not depend on the graph choosen but only on the initial state.
* We have established for which of the configurations ($`\mathrm{\Gamma }_j,|\psi _{in}`$) the free evolution of the system can lead to the optimal value of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$.
* We have given a general (analytical) justification of the different behaviour of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ with respect to the change of sign in the selfinteraction parameter. This has been done by first observing that the graphs $`\mathrm{\Gamma }_j`$ can be classified in bi-partite and non bi-partite. We have thus given an analitical justification of the fact that

$$EP_{\mathrm{\Gamma }_j}^{(0)}(\epsilon )=EP_{\mathrm{\Gamma }_j}^{(0)}(\epsilon )\epsilon [0,20]$$

for all the bi-partite graphs. This relation, as it is highlited by our simulations, is not true for the non bi-partite ones.
* We have given a general explanation for the reduction of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ when the self-interaction term becomes sufficiently greater than the hopping term. This has been done by adopting a two-levels system perturbative approach.
* For the non bi-partite graphs, in the perturbative two-levels picture, we have given a detailed description of the following interesting features highlighted by our simulations. In the perturbative regime, the entangling power of the pair $`(U_{\mathrm{\Gamma }_j},|\psi _{in})`$ is: enhanced (decreased) by a positive (negative) self-interaction when, in the initial state, the particles are all localized in the reference mode ($`n_0=N`$); the opposite holds true when the initial state is uniformly occupied ($`n_i=1,i=1,..L`$).
* In the context of the perturbative two-levels picture we have derived analytical formulas that allow to link $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ to the topological properties of the rooted graphs; this has been done in the general case $`N=L`$ and for both the initial states $`|\psi _{in}=|N_{j=0}^{L1}|0_j`$ and $`|\psi _{in}=_{j=0}^{L1}|1_j`$. Starting from the simulations and with the help of this formulas we have compared and described the behaviour of $`EP_{\mathrm{\Gamma }_j}^{(0)}`$ for different graphs.

We believe that a natural –though not unique – experimental framework in which the theoretical analysis carried over in this paper can find application is provided by systems of cold bosonic atoms in optical lattices.

## VII acknowledgments.

We are grateful to Radu Ionicioiu for useful comments. Special thanks are due to Vittorio Penna for introducing us to the problem of bosonic lattices and useful discussions. We gratefully acknowledge financial support by Cambridge-MIT Institute Limited and by the European Union project TOPQIP