Quantum Entanglement and State-Transference in Fenna–Matthews–Olson Complexes: A Post-Experimental Simulation Analysis in the Computational Biology Domain

Abstract

Fenna-Mathews-Olson complexes participate in the photosynthetic process of Sulfur Green Bacteria. These biological subsystems exhibit quantum features which possibly are responsible for their high efficiency; the latter may comprise multipartite entanglement and the apparent tunnelling of the initial quantum state. At first, to study these aspects, a multidisciplinary approach including experimental biology, spectroscopy, physics, and math modelling is required. Then, a global computer modelling analysis is achieved in the computational biology domain. The current work implements the Hierarchical Equations of Motion to numerically solve the open quantum system problem regarding this complex. The time-evolved states obtained with this method are then analysed under several measures of entanglement, some of them already proposed in the literature. However, for the first time, the maximum overlap with respect to the closest separable state is employed. This authentic multipartite entanglement measure provides information on the correlations, not only based on the system bipartitions as in the usual analysis. Our study has led us to note a different view of FMO multipartite entanglement as tiny contributions to the global entanglement suggested by other more basic measurements. Additionally, in another related trend, the initial state, considered as a Förster Resonance Energy Transfer, is tracked using a novel approach, considering how it could be followed under the fidelity measure on all possible permutations of the FMO subsystems through its dynamical evolution by observing the tunnelling in the most probable locations. Both analyses demanded significant computational work, making for a clear example of the complexity required in computational biology.

1. Introduction

Photosynthesis is a key process intended to obtain energy from sunlight, and is vital for the survival of life on Earth. A series of physical and metabolic reactions occur in specialized structures to absorb, transfer, and store light energy efficiently within photosynthetic organisms. The importance of micro-organisms in anthropic activities has been highlighted nowadays as a third generation of biofuels, particularly those coming from photosynthetic ones [1], especially cyanobacteria [2], which are the main target of the current study due to the important role played by their dependence on temperature [3]. Among those structures, Light-Harvesting Complexes (LHC) are responsible for light capture and energy transference into the Reaction Centres (RC) [4], where electron transference takes place following the photosynthetic process. Due to their almost $100\%$ efficiency, LHCs have been the subject of study, particularly in the quantum realm [5,6,7,8].

The current studies on LHCs cannot proceed without meeting of several necessary technological developments and achieving collaboration among different disciplines. An array of proteins in which the chlorophyll molecules are embedded to transfer light energy in the further processes within photosynthesis, LHCs appear to perform tasks concerning the most basic rules of matter, where our measurement instruments are not yet completely able to reach in understanding such mechanisms more deeply. Thus, a multidisciplinary approach is necessary to combine experimental, theoretical, and computational approaches. Such combined efforts have allowed us to unravel the mechanism behind the aforementioned energy transfer. Photonics has provided the basis for advanced femtosecond laser techniques to trace the excitation energy transfer (EET) dynamics [9,10], while biochemistry and genetics have contributed to the methods used to purify, characterize, and modify the protein structures necessary in spectroscopy essays [11,12,13,14,15,16]. While theoretical methods and computational models allow us to simulate the light absorption behaviour in different structural features [17], this is only affordable through computational biology, which provides observable data to compare with experimental mesoscopic observations. In fact, the computer atomistic modelling of LHC comprises multiscale methods involving dynamical and spectroscopic properties regarding excitation and site energies, excitonic couplings, and spectral densities determining the exciton dynamics and spectroscopic properties in the key aspects of photosynthesis [18].

Accordingly, LHCs are the focus of a multidisciplinary effort to understand their global behaviour. The growing interest in these structures is centred on their biological study, potential photovoltaic applications, and the natural structures that control certain quantum computing processes at room temperature [19,20]. All those processes lie in the common frontier of biology, chemistry, and physics. Biophysics, a discipline requesting all necessary physical theories and methods to understand how a biological system works, remains open to the most elusive phenomena behind life. There, the role of computational biology completes the analysis coming from the experimental data and further theoretical analysis.

Particularly, the Fenna–Matthews–Olson (FMO) complex is a photosynthetic trimer antenna present in anoxygenic photobacteria as the Green Sulfur Bacteria (GSB) [8]. In each monomer, it holds eight Bacteriochlorophylls-a photopigments (BChls-a) inside a protein scaffolding structure. These antennas capture photons, transferring light energy from FMO to RC through their eight BChls. The absorption and transfer mechanisms are not yet completely understood [21]. While thermal excitation in the complete structure alone could possibly be sufficient to understand this basic phenomenon and its efficiency, other analyses suggest the crucial intervention of quantum features to explain it. Thus, quantum biology is a biophysics subdiscipline accounting for the importance of quantum mechanical effects in living systems. FMOs are sufficiently simple life systems to understand the quantum behaviour involved within protein structures. Because certain topics involving the current research, particularly those related to quantum mechanics, could become difficult for inexperienced readers, we provide more friendly descriptions of the meaning of several concepts, phenomena, and equations throughout the presentation.

There exists an interesting series of experimental outcomes leading to a deeper understanding. A water-soluble pigment–protein complex was the first to be analysed by X-Ray spectroscopy [22,23]. In addition to its chemical structure, the identification of several photopigments [24] (the BChls) was attained by electron microscopy, which provided detailed and exact knowledge of their structure [5], distances, and orientations as well as their protein supporting structure. Such knowledge led to the theoretical settlement of the dipole–dipole interactions among BChls and precise quantification of the involved Hamiltonian [9,25] (without considering other secondary time-dependent interactions). This allowed the possibility of inquiring about other more subtle features of this structure using computer simulations [26]. Thus, compared with other observational entities, certain key constants rule their precise behaviour [27]. Through the corresponding computational analysis, at least for quantum purposes, several models analysing their dynamics have been implemented to reach more precise predictions fitting the facts revealed by spectroscopy [28,29]. Detailed research efforts concerning the FMO complex knowledge and their understanding can be found in [30].

Figure 1 shows several views of the FMO complex [31]: (a) a whole view exhibiting its trimer structure, (b) one of the monomers with its protein scaffolding, (c) the same monomer free of its scaffolding, and (d) a rotation of the last monomer, showing the eight BChls in their commonest numbering. The interactions among the BChls are mainly dipolar electric. Note that each BChls has a specific orientation in terms of its dipolar momentum and that the eighth BChl is located inside the monomers (closer to the chlorosome baseplate, the region where the light energy comes in); for this reason, it was the last to be discovered and initial analyses did not consider it. In fact, BChls 3 and 4 are the nearest to the cytoplasmic membrane and RC. On the other hand, BChl 1, BChl 6, and BChl 8 are the closest to the chlorosome baseplate, where the light energy enters into FMO [5,8,32,33]. In fact, the BChl positions and orientations are interestingly associated with their own excitation energy range and their order in the energy transfer path, as observed in spectroscopic studies. Thus, these BChls are mainly responsible for the absorption energy face (the outer antenna towards blue), while those delivering the energy gathering and linked to the RC are shifted towards red [5]. The FMO complex is a photosynthetic component in all species of GSB. Of these species, Chlorobaculum tepidum [34] and Prosthecochloris aestuarii are the most deeply studied, exhibiting slight structural differences among certain strains belonging to each species.

The FMO became the first LHC structure to be mapped under spectroscopy [22,23]. The spectroscopic analysis demonstrated the existence of long-lived quantum coherence inside the FMO complex [35,36,37] together with its role and importance for its high energy transfer efficiency. Nevertheless, it is a matter of debate [38,39] between coherent quantum (e.g., Lindblad equation) and incoherent non-quantum (e.g., Förster theory) approaches to explain the efficient light-harvesting process useful for photosynthesis as well as their different timescales, questioning the relevance of quantum coherence in the process [38,40,41]. On the other hand, several studies have stated that FMO requires both coherent and incoherent energy transfers [36,42,43,44]. Together, theoretical physics and chemistry are developing a better understanding of the mechanisms involved under a quantum perspective by enabling computer simulations in the light-harvesting tracking [35,45,46].

In a simplified interpretation of the phenomenon, BChls interact among them through electrical forces. However, their protein scaffolding structure can be understood as springs between each pair of BChls exerting elastic forces. Thus, light excitation is transmitted through the entire elastic structure. Nevertheless, in the beginning, such transmission has a quantum nature, leaving room for slower classical and mechanical transmission [30]. Thus, the faster quantum mechanical excitation is ruled by features such as state superposition, coherence, and entanglement, which are explained below.

Quantum entanglement constitutes an important feature behind protein structures through the quantum coherence exhibited in the FMO. This phenomenon appears to be responsible for the set of correlated BChl behaviours that jointly transfer the light energy excitation under a process emulating the well-known quantum tunnelling effect. A clearer and more intuitive depiction of such phenomena is properly discussed in the development. This tunnelling effect and its understanding could become useful in the comprehension of the FMO as a quantum channel inside a biological system, as it has been tracked by spectroscopic methods [47]. However, entanglement is an elusive quantum property under the absence of a final absolute measure to comprise its richness and complexity. To inquire about such phenomena, condensed-matter models, including quantum mechanics, should be implemented by proper computer simulations considering their complex structure as a part of their protein scaffolding. Thus, mathematical modelling boosted by experimental observations in spectroscopy on biological samples has allowed post-experimental analyses to simulate those quantum features. In particular, and beyond the biological interest in these simple systems, both quantum entanglement and state transference through each monomer structure are interesting features in quantum information and processing.

The aim of the current research has two main directions inside the computational biology domain, sharing the frontier with quantum biology, first through the analysis of more complex entanglement models that are commonly recurred in the contemporary literature, and second to provide a certain insight into the state transfer through the excitation of light energy crossing the complex structure of the surrounding protein scaffolding. Section 2.1 synthetically presents the basic details of the quantum mechanical excitation modelling through each monomer inside the FMO. Such modelling allows us to simulate the time evolution of the excitation through a non-Markovian method. The obtained outcomes are the building blocks of the following sections. Section 2.3 contains a short review of entanglement measures commonly used or suggested in the study of FMO entanglement. We present a precise timely description of them along with their theoretical relations, followed by an analysis of specific examples with their similarities and discrepancies. Despite most of these measures referring to bipartition-based models, which involve the overall set of BChls, we introduce the maximum overlap with respect to the closest separable state, a foreseeable and promising quantum measure more appropriate for this phenomenon. Nevertheless, it requires considerable computing effort. We obtain the entanglement analysis of several important groups of BChls as a function of time and temperature. The fourth section proposes a model to track the state transference over time as departing from the initial condition stated in the second section to reach the tunnelling effect in the excitation. The final section concludes the main findings.

2. Results and Discussion

2.1. Fundamentals of FMO Dynamics Simulation Using the Hierarchical Equations of Motion

In this section, a brief review concerning the theoretical aspects relevant to modelling the excitation evolution inside the FMO is discussed, as well as other complementary references regarding each fundamental topic. Thus, we first discuss the main considerations about excitation among BChls together with their scaffold structure, then the quantum equations in the open quantum systems realm, where this problem should be analysed, in addition to their computational addressing. With these elements, it is possible to model the excitation evolution under specific initial conditions of excitation, which are be discussed as well. Finally, traditional measures of entanglement inside the group of BChls of each monomer are briefly studied; we improve upon such models in the next section.

When energy transmission through the FMO is settled, the excitation as vibrational modes through the BChls has a mainly quantum nature. This implies that several modes between each pair coexist at the same time in superposition, a feature present exclusively in quantum systems. This generates the energy levels and characteristic transitions depicted in Figure 2a regarding the excitation of each BChl, which means an increase and decrease from/to the ground energy level of each one. The phenomenon can be observed in the experimental spectroscopy analysis depicted below.

2.1.1. FMO Structural Model

The FMO complex is a nanoscale water-soluble protein. It has a molecular weight of 150,000 Da, a maximum diameter of $8.3\mathrm{nm}$ [22,23], and an average distance between nearest-neighbour BChls of 12 Å [13]. The proven approach to model the excitation evolution of BChls inside FMO first considers the group of eight BChls as the core quantum system [30]. In experimental spectroscopy analysis, it has been observed that each BChl exhibits a two-level system behaviour where in the main only one BChl is excited at a time, with bi-excitonic or higher states being rare [48]; this effect is due to and reinforced by the dipole blockade behaviour [49]. Thus, a preferred basis on which to depict the evolution is the occupation basis $|k\rangle \equiv |0_1...0_{k-1}{1}_k{0}_{k-2}...0_8\rangle$; alternatively, the excitonic basis $|{ϵ}_i\rangle ,i=1,...,8$ corresponds to the free Hamiltonian eigenvalues of the system formed by the eight BChls [40].

As the knowledge of precise energies inside and together with the protein structure becomes a theoretical challenge, it is commonly obtained by fitting the optical spectra measured with spectroscopic techniques considering a femtoseconds range [9]. The Hamiltonians are reported in the following subsections and in Appendix A. In addition, the model should consider the interaction with the protein scaffolding. Despite its structural complexity, it can be considered as a continuous media stating a thermal bath, thereby moving the problem into an open system one. Such systems include the bath environment through quantum master equations (Lindblad [28,45] or Redfield [26,28]). More accurate approximations are reached using non-Markovian approaches such as the Hierarchical Equations of Motion (HEOM) [27,50,51]. Finally, to make the simulation implementation easier for the computer systems evolution, the superoperator method converts such problems into a finite set of first-order differential linear equations [30]. We briefly deal with such approaches in the next subsections.

2.1.2. Interactions inside of FMO among Their Structural Components

The Hamiltonians for each monomer in the FMO are obtained by combining spectroscopy measurements and theoretical models. Figure 2a comprises several of the observed series of excitation dynamics. Main energy entries occur in the nearest BChls to the baseplate (1, 6, and 8) [5]. A cascade energy transition begins with the privileged paths shown there. In the more recent studies, BChl 8 appears as the main entry, despite it not being considered in the original analysis (where BChl 1 appeared as the main) because it was discovered later [23,31]. The final energy transfer moves on the RC through BChls 3 and 4 [51].

The full modelling for the interactions in each monomer considers the dipole–dipole interaction among BChls [28]. In spite of this, the structure becomes modified by the protein scaffolding stating new equilibrium positions [30]. Such a development introduces a relocalization energy $H_{\mathrm{reloc}}$ and a set of oscillatory energies around those new equilibrium points $H_S$ for the system of the set of eight BChls. This Hamiltonian has been reported for C. tepidum [24] and P. aestuarii [5]. We include such $H_S$ Hamiltonians in Appendix A. In addition, protein scaffolding (considered as a thermal Bath) has its own energy $H_B$ depending on several external parameters. Finally, other main terms consider the interaction between the System (the set of eight BCHls) and the Bath. Thus, the main Hamiltonian for the entire problem becomes $H=H_S+H_{\mathrm{reloc}}+H_B+H_{S-B}$. Figure 2b depicts graphically the last set of interactions on the monomer $H_S$ for the group of BChls $H_{\mathrm{reloc}}$ for the relocalization energy due to De protein folding, $H_B$ for the state of the phononic bath, and $H_{S-B}$ for the interaction between the main system and the bath.

2.1.3. Lindblad and Redfield Equations for an FMO Monomer

We are interested in the joint behaviour of BChls, which is depicted as a statistical mixture of states through the density matrix written in terms of whether the occupation or the excitonic bases. Clearly, a first approximation is the Liouville–von Neumann equation:

$$\begin{array}{c}\hfill i\hslash \frac{\partial \rho }{\partial t}=[H_S,\rho ]\end{array}$$

(1)

The time dependence of the set of BChls imposes a natural statistical mixture of quantum states. Nevertheless, $\rho ={\rho }_S$, and hence $H=H_S$, is unable to capture the interaction with the protein scaffolding, which is a complex and dense surrounding structure, which suggests treatment as a thermal bath. Otherwise, if $\rho$ represents the entire system, it is not possible to isolate ${\rho }_S$, as in the further time-evolution both the BChls system and the bath become entangled regardless of the initial condition $\rho ={\rho }_S\otimes {\rho }_B$. Concerning our interest in ${\rho }_S$ and the complexity of such a scaffolding, a better approximation is provided by the open quantum system approach, where the protein scaffolding is modelled as a thermal bath. Thus, the corresponding master equations provide a reasonable approximation. For instance, the Lindblad master equation [52]:

$${\dot{\rho }}_S=-\frac{i}{\hslash }[H_S,{\rho }_S]+\frac{1}{{\hslash }^2}\sum _{\alpha }\left(L_{\alpha }{\rho }_S{L}_{\alpha }^{†}-\frac{1}{2}\{{\rho }_S,L_{\alpha }^{†}{L}_{\alpha }\}\right)$$

(2)

where $L_{\alpha }$, the so-called Lindblad operators, should integrate the interaction between the system and the bath. Despite this, such operators do not have a direct expression with the physical arrangement, instead forming a generic basis for the bath. Nevertheless, it has been used as an approximation of the FMO dynamics [43,45,53]. Alternatively, the Redfield master equation [54]:

$${\dot{\tilde{\rho }}}_S=-\frac{i}{\hslash }[H_S,{\rho }_S]-\frac{1}{{\hslash }^2}\sum _{\alpha }[{\Sigma }_{\alpha },{\mathsf{\Lambda }}_{\alpha }{\rho }_S-{\rho }_S{\mathsf{\Lambda }}_{\alpha }^{†}]$$

(3)

with ${\tilde{\Sigma }}_i,{\tilde{\mathsf{\Lambda }}}_i$ being operators which physically depict the bath coupling. This approach has been extensively used to model FMO behaviour [26,28,29].

In both cases, the later master equations are constructed to depict the main system (the set of BChls in each monomer) considering external interactions with the phononic bath (the protein scaffolding). While the Lindblad master equation is mainly constructed on a mathematical basis depicting the external interaction, the Redfield one involves a clearer and closer relationship with the involved physical systems. Nevertheless, there is criticism of such approaches concerning their Markovian nature. In fact, it is believed that the FMO dynamics depend not only on its present state but also on the past states, exhibiting a non-Markovian nature [55]. We deal with such a requirement in the following subsection.

2.1.4. Hierarchical Equation of Motion Model and the Superoperator Approach of Master Equations

The Hierarchical Equation of Motion (HEOM) [50] is a non-Markovian model based on a non-perturbative approach. It was developed to reproduce the evolution of $\rho$ for quantum dissipative systems such as FMO. After the following technical presentation, we discuss a more friendly interpretation of the HEOM method applied to our current problem. HEOM is constructed as a recursive set of equations considering D previous temporal stages as a bath memory. Each $\rho$ stage is labelled with a vector $\mathbf{n}$ of dimension N (in this case $N=8$, the total number of BCHls). Then, ${\rho }_{\mathbf{n}}$ with $\mathbf{n}=(0,..,0)$ corresponds to the real density matrix of the BChls system ${\rho }_S$. Other ${\rho }_{\mathbf{n}}$ are auxiliary matrices regarding the entire global system memory. This set of matrices ${\rho }_{\mathbf{n}}$ has a degree $s=0,1,...,D$ related to all vectors $\mathbf{n}=({\mathrm{n}}_1,{\mathrm{n}}_2,...,{\mathrm{n}}_N)$, where ${0\le \mathrm{n}}_k\in {\mathrm{Z}}^{+}\cup \left\{0\right\}$ with the property ${\sum }_{k=1}^N{\mathrm{n}}_k=s$.

HEOM comprises the main features of the Markovian models as Redfield or Lindblad quantum master equations by including the Hamiltonian $H_S$ and the relocalization energy correction $H_{\mathrm{reloc}}$. In addition, the trapping from BChls 3 and 4 observed in the spectroscopy analysis, together with terms comprising the non-Markovian features, read as follows [30,50]:

$$\begin{array}{ccc}\hfill \dot{{\rho }_{\mathbf{n}}}& =& -\frac{i}{\hslash }[H_S+\sum _{k=1}^N{\lambda }_k|k\rangle \langle k|,{\rho }_{\mathbf{n}}]-\sum _{k=1}^N{\mathrm{n}}_k{\gamma }_k{\rho }_{\mathbf{n}}-r_{\mathrm{trap}}(\left\{\right|3\rangle \langle 3|,{\rho }_{\mathbf{n}}\}+\left\{\right|4\rangle \langle 4|,{\rho }_{\mathbf{n}}\})\hfill \\ & & +i\sum _{k=1}^N{\Delta }_k([V_k,{\rho }_{{\mathbf{n}}_{k+}}]+{\mathrm{n}}_k([V_k,{\rho }_{{\mathbf{n}}_{k-}}]-i\frac{\beta \hslash {\gamma }_k}{2}\{V_k,{\rho }_{{\mathbf{n}}_{k-}}\}))\hfill \end{array}$$

(4)

in the above expression, k is the Boltzman constant and $\beta =kT$; moreover, $V_k=|k\rangle \langle k|$, ${\Delta }_k={\lambda }_k/\beta {\hslash }^2{\gamma }_k$ (${\lambda }_k$ is the relocalization coefficient for $k-$th BChl [30]), and ${\gamma }_k$ can be considered as the bath cutoff frequency characterizing its non-Markovian nature [56]. Thus, ${\gamma }_k$ is the interaction strength between the bath and the $k-$th BChl; because BChls 3 and 4 drive the energy on the RC with a trapping rate $r_{\mathrm{trap}}$, it is represented in (4). There, if $\mathbf{n}$ is a vector of order s, then ${\mathbf{n}}_{k\pm }$ is a vector of order $s\pm 1$ identical to $\mathbf{n}$ but with the component ${\mathrm{n}}_k$ increased (decreased) by one. Together, ${\rho }_{\mathbf{n}}\left(0\right)$ for $\mathbf{n}\ne \mathbf{0}$ are initially settled as the zero matrix. Finally, ${\mathbf{n}}_{k+}=\mathbf{0}$ is selected as the cut-off for $\mathbf{n}$ of order D.

Initial conditions ($t=0$) for the FMO system correspond to ${\rho }_S$ unentangled with the bath. Because it has been observed in spectroscopic studies that BChls 1, 6, and 8 are the antennas of FMO, ${\rho }_S\left(0\right)={\rho }_{\mathrm{LOC}}^k=|k\rangle \langle k|$ with $k=1,6,8$ is widely considered as an initial condition in the literature [5,51,57]. Despite this, another kind of initial condition corresponds to the Förster Resonance Energy Transfer (FRET) [41], a feasible initial condition related to the fluorescence, a mechanism of energy transfer between light-sensitive molecules as the BChls. This suggests energy transference from a main BChl into the entire system as a function of their proper energies [42]:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{FRET}}^i=\sum _{k=1}^N|\langle {ϵ}_k|i\rangle {|}^2|{ϵ}_k\rangle \langle {ϵ}_k|,i=1,...,8\end{array}$$

(5)

where $\left\{\right|{ϵ}_k\rangle |k=1,...,N\}$ are, as explained before, the eigenvalues of $H_S$. Commonly, $i=8$ is used as the main entry of light energy into each monomer of FMO.

Clearly, the HEOM approach has a more complex interpretation due to its non-Markovian basis; it is pictorially summarized in Figure 3. Through several $D+1$ time steps, the complete system (with the BChls represented with the ball and stick models and the bath as protein alpha-helices and beta-sheets in green) evolves ruled by the electric interaction $H_S$ between BChls, in agreement with a theoretical analysis bas on the experimental spectroscopy observation and their mutual interaction. Interactions are shown with rays in different colours representing the quantum excitation modes. This representation can be useful in understanding the entanglement concept in the model. In the upper formula in Figure 3, in a first approximation [30], this interaction shifts the positions of the BChls, providing a relocalization energy $H_{\mathrm{reloc}}={\sum }_{k=1}^N{\lambda }_k|k\rangle \langle k|$ (the shifting of the electrostatic point of equilibrium from red dots to yellow ones). Such an energy reference can be compared with that generated by gravity on a vertical spring by shifting the equilibrium position. In our model, negative trapping energy is introduced in the form of $-r_{\mathrm{trap}}\left\{\right|3\rangle \langle 3|+|4\rangle \langle 4|,{\rho }_{\mathbf{n}}\}$ (note that below this term already includes the non-Markovian HEOM approach), reflecting the final delivery to the RC (the last formula in the figure). In other approaches [57], the RC is modelled as an additional system in the form of a sink, reflecting the excitement from the BChls system. Additionally, the effect of bath energy on the main system is introduced by $H_B=-{\sum }_{k=1}^N{\mathrm{n}}_k{\gamma }_k{\rho }_n$ (the second formula in the figure). Finally, the HEOM approach relates the quantum state in each time step through the Hamiltonian in the third formula to obtain a future quantum state through the derivative of ${\rho }_{\mathbf{n}}$ on the left side in (4).

Because HEOM requires a numerical approach, as is generally the case in all master equations, a computational procedure is convenient. In fact, the HEOM and open quantum system master Equations (2)–(4) can be transformed into a set of differential linear equations. This is known as the superoperator-supervector procedure [30]. In this procedure, the density matrix is considered as a supervector with two indexes, while each pair of matrices multiplying it on each side is considered a tensor product of supervectors, thereby generating a superoperator:

$${\left(A\rho B\right)}_{ij}=\sum _{k,l}{A}_{ik}{\rho }_{kl}{B}_{lj}=\sum _{k,l}{(A\otimes B^T)}_{ij,kl}{\rho }_{kl}$$

(6)

This rule can be applied in each term of (2)–(4), thereby jumping $\rho$ as a supervector of size $N^2=64$ on the right side of each term through other operators and becoming multiplied by a superoperator $\kappa$ of size $N^2\times N^2$ obtained by summing the tensor products appearing in each master equation:

$${\dot{\overrightarrow{\rho }}}_S={\kappa }_{\mathrm{HEOM}}{\overrightarrow{\rho }}_S$$

(7)

Thus, such a single transformation allows us to change any quantum master equation into a linear single matrix differential equation. This procedure allows an easy computer algorithm to b implemented to obtain the evolution of ${\rho }_S$. This strategy is simply a method of translating the HEOM equations into a single linear differential problem which can be easily numerically solved. A more detailed implementation of this procedure and of the HEOM method is reported in [30].

2.1.5. Complementary and Operative Information of the BChls Dynamics Simulation

Certain operative aspects should be stated before the performance of the primary simulation to reach the oscillation dynamics for the eight BChls. Because of the important role of spectroscopy in the FMO dynamics, the constant parameters involved, such as $r_{\mathrm{trap}},{\lambda }_k$, and ${\gamma }_k$, which are typically expressed in energy units, are expressed in spectroscopy units instead (see Section 3).

Following the last statement, the constant $r_{\mathrm{trap}}=1{\mathrm{ps}}^{-1}$ is responsible from the trap efficiency of the BChls 3 and 4 energy towards the end of the process. This value is approximated by spectroscopy observations. Note that in our model we omit the final entrapment by the RC (as an example, see the modelling by [57]); thus, it is expected that energy finally becomes gathered by this pair of BChls. Reorganization energies are ruled by ${\lambda }_k$. This constant reflects the architecture of BChls inside FMO, and is responsible for the decoherence strength. Despite differentiated values for each BChl being expected, our actual knowledge only allows us to estimate them in the range $[35,65]{\mathrm{cm}}^{-1}$ [58,59]. Both extremes are commonly analysed in the literature by taking each value for the overall BChls. Similarly, the bath cutoff frequency ${\gamma }_k$ (the inverse of the bath coherence time-scale) is experimentally approximated at around $50{\mathrm{cm}}^{-1}$ [56,59].

An external physical parameter is the operation temperature T. The entire behaviour of FMO dynamics is influenced by it, including equilibrium populations, quantum coherence time, and energy capture efficiency. Extreme values used to observe the FMO behaviour range from experimental spectroscopy analysis temperatures ($77\mathrm{K}$) to hot water temperatures around 333–343 K (60–70 °C) near the thermal vents where bacteria survives in absence of sunlight [12], living from the dim glow at a depth of $2500\mathrm{m}$. For this reason, a range of $[77,347]\mathrm{K}$ is commonly considered of interest in the simulation analysis. The following HEOM simulations were performed for a depth $D=4$. Despite higher values for D having been performed and reported [27], our results remain in the percentual variation range of ${10}^{-4}$.

2.2. Quantum Dynamics Modelling of BChls Excitation in the FMO Complex

In this subsection, we graphically report the data which are the basis for our entanglement and state transfer analysis. Operative details are summarized in Section 3.

Using the HEOM method, the parameters $r_{\mathrm{trap}}=1{\mathrm{ps}}^{-1},{\gamma }_k=50{\mathrm{cm}}^{-1}$, and the pair of values ${\lambda }_k=35,65{\mathrm{cm}}^{-1}$ (the same values for each BChl, $k=1,...,8$), a set of computer simulations using $D=4$ between $t\in [0,15]\mathrm{ps}$ were performed for P. aestuarii and C. tepidum. A wide group of representative temperatures were used: $T=77,131,185,239,293,347\mathrm{K}$. The use of $\rho \left(0\right)={\rho }_{\mathrm{FRET}}^8$ as an initial condition was implemented in all simulations. Several of these outcomes have been presented (P. aestuarii) in a limited previous analysis of single bipartite entanglement and coherence [60,61]. This work widely extends their analysis by considering another bacteria specie, and particularly by doing so under a more formal and diverse analysis of entanglement.

The meaning of outcomes comprised in $\rho \left(t\right)$ should be interpreted remembering that the evolution process produces a statistical mixture of quantum states, each one being a superposition of excitation modes between pairs of BCHls. Then, the diagonal entries of $\rho \left(t\right)$, ${\rho }_{ii}\left(t\right),i=1,...,8$ represent the probability of each BChl becoming the excited one, considering that only one becomes excited at time t. Outcomes are plotted in Figure 4 (P. aestuarii) and Figure 5 (C. tepidum). In both cases, subfigures (a–f) correspond to $\lambda =35{\mathrm{cm}}^{-1}$ and subfigures (g–l) to $\lambda =65{\mathrm{cm}}^{-1}$. Only the populations ${\rho }_{ii},i=1,...,8$ are shown as functions of $t\in [0,10]\mathrm{ps}$ (the initial part of the simulated dynamics). In each group, the temperatures are ordered from lower to higher. The population for each BChl is presented in colour, in agreement with the colour legend on the right side. An important aspect not shown here (previously analysed [61]) is the increase in decoherence, which means the absence of quantum features through time. In fact, for higher times than $15\mathrm{ps}$ the coherence (the diagonal-off entries in $\rho \left(t\right)$) drops to zero, reaching a classical statistical mixture of physical states (populations without quantum properties).

The effect of the balance on the final populations ${\rho }_{33}$ and ${\rho }_{44}$ is notable when the temperature rises. In addition, the reduction of coherence timelapse is notable with its increase (asymptotic behaviour of populations). The effect of ${\lambda }_k$ occurs in combination with T, noting that a higher ${\lambda }_k$ value accelerates the dynamics for the lower temperatures and slows down with the higher ones. We consider the entire data of those numerical simulations in order to inquire more deeply into the entanglement phenomenon among BCHls under several models of entanglement measures, as well as both semi-locally and globally. The development and outcomes are presented in the next section.

2.3. Entanglement and Quantum Coherence Measures

The eight BChls in each monomer of FMO evolved according to the HEOM stated in (4) and the initial condition ${\rho }_{\mathrm{FRET}}$ are considered. The states of such sets of BChls become mixed and entangled in general. Entanglement is a quantum feature that is not present in classical systems. While classical systems exhibit certain independence among physical dynamic variables or observables, such variables can only be weakly related through functional relations. Instead, certain groups of quantum observables can be strongly related in a superposition of such correlations. This means that upon measurement of just one of those observables, the remainder becomes strictly determined. Thia property belongs to the quantum state, not to the system itself. When no correlations are present among the related observables, we can say that the quantum state is separable.

Nevertheless, this apparent single feature becomes highly complex to quantify when the number of subsystems increases. For instance, for a pair of systems exhibiting two levels (with observables as non-excited or excited for each one, as BChls are), there is a single type of entanglement when considering both systems together. For the same case, three systems could exhibit correlations by pairs or for all three together. For eight BChls, each one excited or not, the type of correlations increases dramatically, by pairs, by thirds, etc., and then again through subgroups having correlations. Thus, the quantification of entanglement becomes highly complex, as there exists no clear measure comprising this richness.

Only scarce measures of bipartite pure systems are well known (those not involving statistical mixtures of quantum states), despite there being certain useful measures for a limited group of subsystems. For mixed states, the situation is more intricate around the quantification of entanglement, as there are both statistical and quantum correlations in the same system. The success of bipartite systems lies in their simplicity as well as in that of the Peres—Horodecki criterion [62,63] for a pair of subsystems. Clearly, such a criterion could be applied to larger systems by considering bipartitions, which means splitting the system into two parts, each one possibly containing more than a single system. Such a measure provides the correlation only between the subsystems, not inside each one.

For this reason, the most common measures considered to quantify entanglement in each monomer of the FMO are based on bipartitions. In this section, we pursue an analysis of the entanglement inside the FMO to characterize it as a function of parameters or their biological properties.

2.3.1. Analysis under Several Measures of Global and Bipartite Entanglement

In [30], an introductory analysis was presented for models with $N=7$; however, a deeper analysis can be performed using improved entanglement measures. In the literature, the more immediate and recurrent measures correspond with the concurrences [26,30,60,64]:

$$\begin{array}{ccc}\hfill {\mathcal{C}}_{\left\{kl\right\}}& \equiv & \mathcal{C}\left({\rho }_{\left\{kl\right\}}\right)=2\left|{\rho }_{kl}\right|\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill {\mathcal{C}}_{\left\{k\right\}}& \equiv & \mathcal{C}\left({\rho }_{\left\{k\right\}}\right)=2\sqrt{{\rho }_{kk}(1-{\rho }_{kk})}\hfill \end{array}$$

(9)

both obtained by partial tracing and using the criterion developed in [27,63]. These measures, the concurrences, are simply mathematical functions quantifying the bipartite entanglement in a specific order, i.e., zero for separable states and one for maximal entanglement, with tighter correlations. They are a kind of measure belonging to the monotones of entanglement, a wider group of ordered measures ranging from separability (minimum) to maximal entanglement (maximum).

In the previous formulas, ${\rho }_{\left\{A\right\}}$ is the tracing of $\rho$ on the entire system, with the exception of those subsystems in the set A. Both measures report entanglement in the range $[0,1]$, 0 for separability and 1 for maximal entanglement. For these reasons, sometimes the squares of (8) and (9) are considered as the concurrence definitions. In general, they provide certain local information about entanglement (between two BChls or between one and the rest, respectively); however, they remain limited. In particular, the local paired concurrence ${\mathcal{C}}_{\left\{kl\right\}}$ becomes naturally related with the coherence [65] of the system (using the measure $l_1$-norm): $C_{l_1}\left(\rho \right)={min}_{\sigma \in \Gamma }{\sum }_{ij}|{\rho }_{ij}-{\sigma }_{ij}|=2{\sum }_{i<j}\left|{\rho }_{ij}\right|$ [66]. Quantum coherence is understood as the stage at which a quantum system exhibits a superposition of classical physical states; otherwise, the system exhibits decoherence. For a mixed state, it may be noted that when only a statistical mixture of states remains, as in our case, it means that $\rho \left(t\right)$ becomes practically a diagonal matrix in the occupation basis (their diagonal-off entries are practically zero). In that sense, both ${\mathcal{C}}_{\left\{kl\right\}}$ and $C_{l_1}\left(\rho \right)$ respectively function as local and global entanglement measurements [27].

Other partitioned or global entanglement measures have been implemented in the entanglement analysis of FMO. For instance, in [67], logarithmic negativity was employed by defining a bipartition $k,N-k$ of BChls $A|B$: $\{{\mathrm{BChl}}_{i_1},...,{\mathrm{BChl}}_{i_k}\}\to A=\{i_1,...,i_k\}$ and $\{{\mathrm{BChl}}_{i_{k+1}},...,{\mathrm{BChl}}_{i_N}\}\to B=\{i_{k+1},...,i_N\}$. The expression for the logarithmic negativity for this bipartition considering the single excitation states model becomes

$$\begin{array}{ccc}\hfill E_{A|B}^{\mathrm{LN}}\left[\rho \right]& =& {log}_2\left(1+2\sqrt{\sum _{i\in A}\sum _{j\in B}{\left|{\rho }_{ij}\right|}^2}\right)={log}_2\left(1+\sqrt{E_{A|B}^{\mathrm{N}}\left[\rho \right]}\right)\hfill \end{array}$$

(10)

$$\begin{array}{ccc}& & \mathrm{with}:E_{A|B}^{\mathrm{N}}\left[\rho \right]\equiv \sum _{i\in A}\sum _{j\in B}{\mathcal{C}}_{\left\{ij\right\}}^2\hfill \end{array}$$

(11)

which is a bipartitioned entanglement measurement. Note that the inner double sum on the partition elements inside the square root, $E_{A|B}^{\mathrm{N}}\left[\rho \right]$, provides a cumulative entanglement measure on each partition, the negativity.

This measure clearly remains in the interval $[0,\left(\genfrac{}{}{0pt}{}{N}{k}\right)]$. The sum over all partitions can provide a global measurement achieved with little computing effort (note that in the partitions for $k=[N/2]=N/2$ only half of them are different). In any case, this measure accounts for an additive amount of overall paired entanglement among all BChls. In other approaches, the relative entropy of entanglement has been used [68]:

$$\begin{array}{ccc}\hfill E^{\mathrm{RE}}\left[\rho \right]& =& -\sum _{i=1}^N{log}_2{\rho }_{ii}-S\left(\rho \right)\hfill \end{array}$$

(12)

$$\begin{array}{ccc}& & \mathrm{with}:S\left(\rho \right)=-\mathrm{Tr}\left(\rho {log}_2\rho \right)\hfill \end{array}$$

(13)

providing a global entanglement measure, which regards the information contained in the state. Another global measure is the Meyer–Wallach one for the system of BChls. It is a monotone measure [69] considered as an entanglement measure for several two-level systems together. For N BChls, the Meyer–Wallach entanglement measure is based on the entanglement of each BChls with the remainder, as follows:

$$\begin{array}{ccc}\hfill E^{\mathrm{MW}}\left[\rho \right]& =& \frac{1}{N}\sum _{k=1}^{N}2(1-\mathrm{Tr}\left[{\rho }_{\left\{k\right\}}^2\right])=\frac{1}{N}\sum _{k=1}^{N}4{\rho }_{kk}(1-{\rho }_{kk})=\frac{1}{N}\sum _{k=1}^N{\mathcal{C}}_{\left\{k\right\}}^2\hfill \end{array}$$

(14)

becoming the average quadratic concurrence ${\mathcal{C}}_{\left\{k\right\}}^2$ of each BChl. Here, as depicted before, ${\rho }_{\left\{k\right\}}$ is the reduced density matrix for the k-th BChl obtained after tracing out all the remaining qubits [30,60]:

$$\begin{array}{ccc}\hfill {\rho }_{\left\{k\right\}}& =& {\mathrm{Tr}}_{\left\{k\right\}}\left(\rho \right)=(1-{\rho }_{kk})|0\rangle \langle 0|+{\rho }_{kk}|1\rangle \langle 1|\hfill \end{array}$$

(15)

In this interpretation, it can be visualized as an average measure of how each BChl remains entangled with its BChl partners. Nevertheless, this measure has inherent deficiencies, as it is unable to distinguish fully non-separable states from those with separable subsystems. Despite this, it is useful for obtaining a global vision of entanglement in the group of N BChls [70].

An alternative approach to provide a global measure of entanglement is a weighted average of entanglement entropies upon all bipartitioned groups of BCHls [70,71]:

$$\begin{array}{ccc}\hfill E_{\tau }^{\mathrm{WAE}}\left[\rho \right]& =& \frac{1}{[N/2]}\sum _{m=1}^{[N/2]}{\tau }^{\left(m\right)}\left[\rho \right]\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill {\tau }^{\left(m\right)}\left[\rho \right]& =& \frac{1}{N^{\left(m\right)}}\sum _{n=1}^{N^{\left(m\right)}}\tau \left(n\right)\left[\rho \right]\hfill \end{array}$$

(17)

Each m-bipartition splits the set of N BChls in two groups, for instance, $\{{\mathrm{BChl}}_{i_1},...,{\mathrm{BChl}}_{i_m}\}$ and $\{{\mathrm{BChl}}_{i_{m+1}},...,{\mathrm{BChl}}_{i_N}\}$. An $m-$bipartition has $N^{\left(m\right)}$ different elements, provided by $\left(\genfrac{}{}{0pt}{}{N}{m}\right)$ if $m\ne N/2$ and $\frac{1}{2}\left(\genfrac{}{}{0pt}{}{N}{m}\right)$ if $m=N/2$. In these expressions, $[·]$ is the floor function. Summing the total number of bipartitions, we obtain ${\sum }_{m=1}^{[N/2]}{N}^{\left(m\right)}=2^{N-1}-1$. Despite $\tau \left(n\right)\left[\rho \right]$ being any bipartite entanglement measure, the normalised negativity$\tau \left(n\right)\left[\rho \right]$ [63] has commonly been used [70]:

$$\begin{array}{ccc}\hfill \tau \left(n\right)\left[\rho \right]& =& \frac{2}{2^m-1}\sum _i\left|{\lambda }_i\right|\hfill \end{array}$$

(18)

which is the normalised sum of the negative eigenvalues ${\lambda }_i$ of the partial transpose matrix associated with each bipartition. Despite being global, each weighted contribution ${E_{\tau }^{\mathrm{WAE}}}^{\beta }\left[\rho \right]={[N/2]}^{-1}{N^{\left(m\right)}}^{-1}\tau \left(n\right)$ can be considered as a specific entangled measure for the corresponding m-bipartition $\beta =i_1,...,i_m|i_{m+1},...,i_N$. It is a more complex partitioned entanglement measure than (10), requiring a more powerful computer effort. Despite this, the main contribution of this measure is the weighted interpretation to reach an entanglement average; nevertheless, the exchangeable normalised negativity becomes more difficult to manage.

The last formulation can be alternatively applied using $E_{A|B}^{\mathrm{N}}\left[\rho \right]$ instead of $\tau \left(n\right)\left[\rho \right]$, thereby defining $E_{E^{\mathrm{N}}}^{\mathrm{WAE}}\left[\rho \right]$. While this measure uses the same set of bipartitions, it is a more economic measure in terms of computer calculations. In a further development, the complete set of bipartitions $\mathcal{P}$ is assumed in the following natural order:

$$\begin{array}{ccc}\hfill \mathcal{P}=\{{\mathcal{P}}_1,{\mathcal{P}}_2,{\mathcal{P}}_3,{\mathcal{P}}_4\}& \equiv & \left\{\right\{\left\{1\right\},\left\{2\right\},...,\left\{8\right\}\},\hfill \\ & & \left\{\right\{1,2\},...,\{7,8\left\}\right\},\hfill \\ & & \left\{\right\{1,2,3\},...,\{6,7,8\left\}\right\},\hfill \\ & & \left\{\right\{1,2,3,4\},...,\{1,6,7,8\left\}\right\}\}\hfill \end{array}$$

(19)

In turn, each element of ${\mathcal{P}}_m$ is labelled as ${\mathcal{P}}_m^n,n=1,...,N^{\left(m\right)}$. Then, the last construction, $E_{E^{\mathrm{N}}}^{\mathrm{WAE}}\left[\rho \right]$, can be made regarding only the main part prescription $A={\mathcal{P}}_m^n$ of order $m=1,2,...,[N/2]=4$, as the remaining part B is easily obtained from each A as its complement $B=A^C$. Note that ${\mathcal{P}}_4$ only contains the different bipartitions, one-half of the total. Thus, the number of elements in each order of bipartitions is $8,28,56,35$ for $k=1,...,4$, respectively, for a total of 127 bipartitions. Regarding the entire presentation at this point, it is noticeable that the entanglement measures presented before are related to the simplest measures in (8) and (9) as an additive or averaged quantity in each case, based on bipartitions, mainly because of the Peres–Horodecki criterion.

Thus, in the following analysis, we have initially compared $E_{A|B}^{\mathrm{N}}\left[\rho \right]$ (averaged as in $E_{E^{\mathrm{N}}}^{\mathrm{WAE}}\left[\rho \right]$) and $E^{\mathrm{MW}}\left[\rho \right]$ in Figure 6 (P. aestuarii) and Figure 7 (C. tepidum) in $[0,5]$ ps, where the entanglement values become more important. Note that $E_{A|B}^{\mathrm{N}}\left[\rho \right]$ is a direct summative entanglement property based on the entire possible bipartitions settled by $\mathcal{P}$ and measured by the paired squared concurrences ${\mathcal{C}}_{\left\{kl\right\}}$ between single BChls belonging to each group of the bipartition. Thus, this cumulative quantity can be represented as piled to obtain a global entanglement measurement considering the entire set of bipartitions, as in these figures. This allows us to identify the contribution of each bipartition and each bipartition order at each time in the evolution. Moreover, the last piling is not only the simple sum of each $E_{A|B}^{\mathrm{N}}\left[\rho \right]$; instead, each contribution is re-scaled as in $E_{\tau }^{\mathrm{WAE}}$, except more precisely, using $E_{A|B}^{\mathrm{N}}\left[\rho \right]$ as $\tau \left(n\right)$. While, $E^{\mathrm{MW}}\left[\rho \right]$ is a global entanglement measure based on bipartitions, in this case it uses the single average of the entanglement squared concurrence of each BChl with the remaining system. Clearly, in the complexity of the entanglement measurement for multipartite systems, both measures provide different information about this feature; however, they both try to provide a global entanglement measurement.

Panels (a–c) show the case for ${\lambda }_k=35{\mathrm{cm}}^{-1}$ and (d–f) show the case for ${\lambda }_k=65{\mathrm{cm}}^{-1}$. Each set corresponds to $T=77,181,293\mathrm{K}$. In all cases, ${\gamma }_k=50{\mathrm{cm}}^{-1}$. $E^{\mathrm{MW}}\left[\rho \right]$ is a global entropy measure obtained from the sum of entropies between each BChl and its remainder. It is shown by the dashed orange line, and should be read on the right scale. On the other hand, $E_{A|B}^{\mathrm{N}}\left[\rho \right]$ is a bipartite measure between each bipartition in (19). It is shown (scaled as in the average of $E_{E^{\mathrm{N}}}^{\mathrm{WAE}}\left[\rho \right]$) for each set of bipartitions with orders $k=1$ (red), $k=2$ (green), $k=3$ (blue), and $k=4$ (grey), and should be read on the left scale.

Each bipartition corresponds to a line from below (darker) to above (lighter), in agreement with (19). They are shown piled on each other in an additive way, with each one scaled by the factor $\frac{1}{4N^{\left(k\right)}}$, as for $E_{\tau }^{\mathrm{WAE}}\left[\rho \right]$. Thus, the gap between lines provides each individual measure integrated into one global measure. The top grey line provides the global measure $E_{E^{\mathrm{N}}}^{\mathrm{WAE}}\left[\rho \right]$. Despite the different scales, note that for larger values of T this global measure reaches its maximum at the beginning of the dynamics. Despite these measures being monotones and having been re-scaled, we can note the different information provided by each one; for example, $E^{\mathrm{MW}}\left[\rho \right]$ is in a sense an average measure about the even filling of populations (${\rho }_{ii}=\frac{1}{N}=\frac{1}{8},i=1,...,8$): $E^{\mathrm{MW}}\left[\rho \right]=\frac{4}{8}(1-{\sum }_{i=1}^8{\rho }_{ii}^2)\le \frac{7}{16}=0.4375$. The maximum value is almost reached by P. aestuarii around 1ps, under $T=77\mathrm{K}$, and C. tepidum almost reaches a value between $2-3$ps. In another issue, $E_{E^{\mathrm{N}}}^{\mathrm{WAE}}\left[\rho \right]$ exhibits larger values for certain bipartitions if the gap between lines are noticeable (see the upper one for such gap). Observing Figure 6 and Figure 7 panels (a) and (d), such gaps are present for BChl 8 in the beginning. Afterwards, they are apparently more spaced as a consequence of the averaging re-scaling. Note the dramatic differences imposed by the increase of T. In general, there is no apparently privileged order k for the bipartitions in terms of a larger entanglement provided by this measure.

In general, despite both being based in bipartitions, it could be said that $E_{E^{\mathrm{N}}}^{\mathrm{WAE}}\left[\rho \right]$ comprises more entanglement cases than $E^{\mathrm{MW}}\left[\rho \right]$ due to being mainly oriented to the entanglement of single BChls with the remaining system. Regarding the maximum observed in the last measure in comparison with the minimum observed in $E_{E^{\mathrm{N}}}^{\mathrm{WAE}}\left[\rho \right]$, this shows that each BChl remains entangled with the remainder at the beginning of the dynamics. Instead, a tiny entanglement between all possible bipartition appears followed by a remarkable increase, which then fades with the increase in temperature.

This brief review provides an illustrative vision of entanglement measures, as well as those composed of single measures such as ${\mathcal{C}}_{\left\{kl\right\}}$ and ${\mathcal{C}}_{\left\{k\right\}}$. Note that $E^{\mathrm{MW}}\left[\rho \right]$ and $E_{E^{\mathrm{N}}}^{\mathrm{WAE}}\left[\rho \right]$ were selected as two bipartite emblematic measures regarding the entanglement between each BChl and its environment as well as between each pair of subgroups of BChls.

2.3.2. The Maximum Fidelity with Respect to the Closest Separable State: An Authentic Multipartite Entanglement Measure

Until now the entanglement of the BChls has been characterized in terms of measures based on the bipartitions of the system. However, in this section we discuss the entanglement quantification of subsystems composed of three components in terms of fidelity with respect to the closest separable state. This is a subtle distinction inside the multipartite entanglement complexity, because while bipartition-based entanglement measures refer to entanglement regarding the overall set of BChls, subgroups of them can instead possibly involve quantum correlations independently of the remaining BChls. For the measure, this quantity should be able to detect fully separable states between the BChls in a subgroup of all of them.

We begin with several definitions. A density matrix $\sigma$ is called fully separable if it can be written as the following convex combination [72]:

$$\sigma =\sum q_k{\sigma }_k^{\left(1\right)}\otimes \cdots \otimes {\sigma }_k^{\left(N\right)},\sum _k{q}_k=1,q_i\ge 0,$$

(20)

where the density operator ${\sigma }_k^{\left(i\right)}$ acts on the space of the i-th BChl. For an arbitrary multipartite mixed state $\rho$, the maximum fidelity with respect to the closest separable state is defined as

$$\mathsf{\Lambda }\left(\rho \right)=\underset{\sigma \in S}{max}{\mathrm{Tr}}^2\left(\sqrt{\sqrt{\rho }\sigma \sqrt{\rho }}\right),$$

(21)

where S stands for the set of fully-separable density matrices [73]. The quantity (21) is bounded between 0 and 1, and the largest value is attained for states of the form (20). Note that the closer $\mathsf{\Lambda }$ is to zero, the more entangled the state is. For this reason, sometimes it is useful to use instead the geometric measure of entanglement defined as $E_G\left(\rho \right)=1-\mathsf{\Lambda }\left(\rho \right)$ [74]. Thus, this measure is free of the partitioned principle to understand entanglement. Instead, it recurs to the separability concept and a metric approach measuring the minimum geometrical distance of each state to the closest separable state through the fidelity measure. This natural measure is in the scale $[0,1]$, from the closest state to the farthest one to the separable state subspace in the Hilbert space regarded.

In the case of two-qubit mixed systems, $\mathsf{\Lambda }$ can be evaluated in terms of concurrence [74]. However, a closed expression in the multipartite scenario is not known even for pure states, and one has to rely on numerical procedures. Our purpose is to analyse the maximum fidelity with respect to the closest separable state $\mathsf{\Lambda }$ for subsystems compounded of three BChls, whose density matrix is in general mixed. Accordingly, we implemented a numerical algorithm which first computed a generalized Schmidt decomposition for multi-qubit pure states, then optimized over all the possible decompositions of a mixed state using Uhlmann’s theorem [73].

Representative results of the $\mathsf{\Lambda }$ time-evolution are shown in Figure 8 and Figure 9 for several values of the temperature. In each case, panels (a,c) correspond to the BChls 3, 6, and 7, commonly the three with the largest ending populations, for ${\lambda }_k=35{\mathrm{cm}}^{-1},65{\mathrm{cm}}^{-1}$, respectively, while panels (b,d) show the subsystem compounded with the energy-exit BChls, 3 and 4, together with the main entry energy BChl, 8. This correlates the main entry and the supposed main exits in the spectroscopy observations. In this case, ${\lambda }_k=35{\mathrm{cm}}^{-1},65{\mathrm{cm}}^{-1}$ with respect to the panels. In both cases, the three-qubit reduced density matrix is computed as ${\rho }_{\{j,k,\ell \}}={\mathrm{Tr}}_{\{j,k,\ell \}}\left(\rho \right)$. As before, the partial trace is taken over the complement of the set $\{j,k,\ell \}$. In all the cases, ${\gamma }_k=50{\mathrm{cm}}^{-1}$ and the time-interval shown corresponds with when the more relevant entanglement behaviour happens. We first note that the initial state ${\rho }_{\mathrm{FRET}}^8$ (5) is not fully-separable, as it is the convex combination of the eigenvectors $|{ϵ}_k\rangle$, which in general are entangled. Thus, in general, any partition of the initial state attains a non-vanishing entanglement degree. As a reference, the $\mathsf{\Lambda }$-values for the eigenvectors are shown in

Table 1. The maximum fidelity with respect to the closest separable state of the eigenvectors $|{ϵ}_k\rangle$ for the BChls P. aestuarii (centre) and C. tepidum (right).

$|{\mathit{ϵ}}_{\mathit{k}}\rangle$	${\mathbf{\Lambda }}_{\mathit{P}. \mathit{aestuarii}}\left(\right|{\mathit{ϵ}}_{\mathit{k}}\rangle \langle {\mathit{ϵ}}_{\mathit{k}}\left|\right)$	${\mathbf{\Lambda }}_{\mathit{C}. \mathit{tepidum}}\left(\right|{\mathit{ϵ}}_{\mathit{k}}\rangle \langle {\mathit{ϵ}}_{\mathit{k}}\left|\right)$
$|{ϵ}_1\rangle$	0.9888	0.9909
$|{ϵ}_2\rangle$	0.7050	0.7683
$|{ϵ}_3\rangle$	0.9287	0.9327
$|{ϵ}_4\rangle$	0.7159	0.8604
$|{ϵ}_5\rangle$	0.5719	0.9397
$|{ϵ}_6\rangle$	0.8315	0.8275
$|{ϵ}_7\rangle$	0.8110	0.9454
$|{ϵ}_8\rangle$	0.8292	0.9235

.

For later times, the state is in general close to a fully-separable state. This is particularly evident in Figure 9a,c, where the values in the vertical axis fluctuate around the worked numerical precision. Nevertheless, in the remaining cases the model allows us to produce a more entangled state depending on the temperature. For instance, in Figure 8a,c, the higher the temperature, the more entangled the BChl state of 3, 6, and 7 becomes, reaching the lowest $\mathsf{\Lambda }$ for $T=77\mathrm{K}$ at around $9.05\mathrm{ps}$ and $3.2\mathrm{ps}$ for ${\lambda }_k=35{\mathrm{cm}}^{-1}$ and ${\lambda }_k=65{\mathrm{cm}}^{-1}$, respectively. In the case of BChls 3, 4, and 8, panels b,d of Figure 8 and Figure 9 show dependence of the entanglement on both the temperature and ${\lambda }_k$.

2.3.3. Entanglement Measures Using Bipartitions against Authentic Multipartite Measures in FMO

It is well known that in the multipartite case there is no universal measure to quantify entanglement. In this sense, each one of the proposed measures provides different information on the non-local correlations. Accordingly, no unique maximally entangled state can be identified in general [75]. The present analysis has been focused on three-qubit states, which stand for the simplest multipartite entanglement scenario and represent, to our knowledge, one of the first steps in the characterization of entanglement in the FMO systems beyond measures based on bipartitions.

The current measure, the maximum fidelity with respect to the closest separable state (or alternatively the geometric measure of entanglement), has not been used before in the analysis of FMO entanglement, and genuinely shows the clear existence of entanglement among three BCHls (used here as an example because of the algorithmic complexity of higher subsystems). This measure is not directly shown by other measures previously used in the literature, which in fact provide other entanglement features, typically only involving bipartitions such as those constructed in Section 2.3.1.

Despite being tiny (near to one), again, a cumulative effect could be expected for the bipartition-based versions. For instance, $E_{A|B}^N\left[\rho \right]$ reflects the entanglement of several k-sized bipartions in an additive way (Figure 6 and Figure 7). For P. aestuarii, they clearly grow at the end before decreasing, particularly for the lowest temperatures. Those for $k=1$ provide higher contributions, mainly due to the lower number of cases. For C. tepidum, this behaviour grows again for a second time at the very end. Nevertheless, for both species the main contribution grows as k grows. Instead, the global measure $E^{\mathrm{MW}}\left[\rho \right]$ exhibits a single peak behaviour in its growth, commonly delayed with respect to $E_{A|B}^N\left[\rho \right]$. This measure, unlike the previous ones, only considers the entanglement of each BChl, with the rest in agreement with (14). In both measures and both species, a remainder entanglement appears due to the tiny correlations between bipartitions or the entire sets of BChls. When the maximum fidelity with respect to the closest separable state is observed between BChls of a subset of the BChls in the monomer, the tiny entanglement there is verified (although we show only the authentic entanglement among the three emblematic groups of BChls).

Thus, there is a global climax when all single entanglements of each BChl with the remaining ones is reached (provided by $E^{\mathrm{MW}}\left[\rho \right]$, the red dashed line in Figure 6 and Figure 7). Such a climax does not equate to that reached when most bipartitions become entangled on average (indicated by $E_{E^{\mathrm{N}}}^{\mathrm{WAE}}\left[\rho \right]$ in Figure 6 and Figure 7). The last condition appears prevalent at the end of the dynamics in different strengths as a function of T. In another trend, authentic tripartite entanglement for BChls 3, 6, and 7 and 3, 4, and 8 initially increases, then decreases again during the dynamics ($\mathsf{\Lambda }\left[\rho \right]$ decreases) as a function of T and at different timescales, as Figure 8 and Figure 9 suggest. All of these are observed behaviours resulting from the different entanglement measures used.

In any case, it is noticeable that a fundamental feature of GBS, as entanglement is, can be finally tracked with the convergence of experimental techniques, theoretical physics, mathematical modelling, and computational biology [76]. This remarkable effort to understand the more basic quantum behaviour behind the photosynthetic processes for such bacteria exhibits the importance of multidisciplinary research merging biology, physics, mathematics, and computer science [30].

2.4. State Transference Analysis

This section analyses the state’s transference through the dynamic process. In fact, the initial state ${\rho }_{\mathrm{FRET}}^8$ is distributed through the localization states with the three central populations: ${\rho }_{8,8}\left(0\right)=0.707,{\rho }_{1,1}\left(0\right)=0.216,{\rho }_{2,2}\left(0\right)=0.065$ for P. aestuarii and ${\rho }_{8,8}\left(0\right)=0.703,{\rho }_{5,5}\left(0\right)=0.080,{\rho }_{1,1}\left(0\right)=0.065$ for C. tepidum. The analysis tracks such states distributed on the set of BChls in a different arrangement. Thus, by considering all BChls permutations, we have

$$\begin{array}{c}\hfill \mathsf{\Pi }=\left\{{\pi }_i\right|i=1,...,8!\}=(\{1,2,...,7,8\},\{1,2,...,8,7\},...,\{2,1,3,...8\},...,\{8,7,...,1\})\end{array}$$

(22)

ordered place by place, which means that each subgroup of permutations is ordered first by its first position inside, again by the second, and so on. There are $8!=40,320$ different permutations. Then, we obtain the fidelity ${\mathcal{F}}_{{\pi }_i}$ for each permutation:

$$\begin{array}{c}\hfill {\mathcal{F}}_{{\pi }_i}\left(t\right)={\mathrm{Tr}}^2\left(\sqrt{\sqrt{{\rho }_{\mathrm{FRET}}^8}\rho \left(t\right)\sqrt{{\rho }_{\mathrm{FRET}}^8}}\right)\end{array}$$

(23)

and as function of t through the evolution. In each t, we obtain the maximum fidelity on the entire set of BChls permutations: ${\mathcal{F}}_{\mathsf{\Pi }}\left(t\right)={\max }_{\mathsf{\Pi }}{\mathcal{F}}_{{\pi }_i}\left(t\right)$. Clearly, ${\mathcal{F}}_{\mathsf{\Pi }}\left(t\right)$ occurs in each t for a different permutation with the time interval. As is well-known, fidelity is a positive valued measure providing a comparison between two states, pure or mixed, obtained as the inner product between them. When both states fit perfectly, this measure becomes one; otherwise, if states are orthogonal (which is understood as the most different possible state), it drops to zero. By tracking the time that represents the best fitting of the initial state as the energy excitation moves into the whole system, we can try to identify which permutation of BChls best fits that initial state at each time. In a sense, this is a discrete comparison of the continuous forms of a wavefunction when tunnelling through a potential obstacle. The outcomes of such tracking are analysed in the next subsections.

2.4.1. Comparative: Input versus Output States

The analysis outcomes are presented by remarking on the permutation ${\pi }_i$ nearest to the initial state as measured by ${\mathcal{F}}_{\mathsf{\Pi }}\left(t\right)$. Thus, Figure 10 (P. aestuarii) and Figure 11 (C. tepidum) show the state transference dynamics in $t\in [0,12]\mathrm{ps}$ for the representative temperatures $T=77,185,293\mathrm{K}$. In each case, panels (a–c) correspond to ${\lambda }_k=35{\mathrm{cm}}^{-1}$ and panels (d–f) correspond to ${\lambda }_k=65{\mathrm{cm}}^{-1}$. In all cases, ${\gamma }_k=50{\mathrm{cm}}^{-1}$. In each plot, each different coloured segment corresponds to the closest permutation carrying out a similar state to the initial ${\rho }_{\mathrm{FRET}}^8$. Colours are reported in the colour bar at the top. The scale reports the initial BChls in the permutation in agreement with the order settled in (22), and each unitary interval in the bar comprises $7!$ permutations.

In each panel, the dominant asymptotic permutation is reported in red on the top right. Thus, considering that the initial permutation is $\{1,2,...,8\}$, we observe that for P. aestuarii the initial dominant excitations in BChls 8, 1, and 2 finally move on $(4,7,5)$, $(7,3,6)$, $(3,6,7)$, or $(4,5,7)$, respectively, depending on the set of parameters $T,{\lambda }_k$ (those BChls in the original positions of the initial excitations). For C. tepidum, the initial dominant excitations in BChls 8, 5, and 1 finally move on $(6,4,5)$, $(4,6,5)$, or $(3,2,5)$, respectively.

Interesting facts arise here. First, note the similarities between the corresponding panels in Figure 10 and Figure 11. In addition, the speedup in the process is imposed by temperature, together with a slower recovery into an similar state to the initial one (${\mathcal{F}}_{\mathsf{\Pi }}$ lowering). A similar effect is produced by the increase of the ${\lambda }_k$ value. Note that for the lower temperature the similarity transitions are more acute. All cases in the dynamics exhibit a kind of excitation dispersion followed by a regrouping on other BChls, which is more uniform for C. tepidum than for P. aestuarii.

2.4.2. Analysis of Main Excitation Paths

In an alternative view, we provide a global vision of the analysis of the whole fidelities for all permutations. Figure 12 (P. aestuarii) and Figure 13 (C. tepidum) show the fidelities calculated for the entire $8!$ permutations in the horizontal axis, evolving on the vertical one for the first 8ps. Each fidelity is shown as a point in agreement with the colour bar scale at the bottom, with 0 (red) and 1 (blue) representing the fidelity. In each figure, panels (a–c) correspond to ${\lambda }_k=35{\mathrm{cm}}^{-1}$ and $T=77,181,293\mathrm{K}$ and (d-f) to ${\lambda }_k=65{\mathrm{cm}}^{-1}$ for the same corresponding T. In all cases, ${\gamma }_k=50{\mathrm{cm}}^{-1}$. Points are horizontally arranged as a function of the permutation using a representative scale showing the integers 1–8, the first number in the permutation; thus, each unitary space contains $7!$ permutations, meaning that the bluest vertical traces show the most likely permutations to which the initial excitation shape is moving.

Plots should be interpreted as follows. Note the permutation 1 ($1,\dots ,8$ in fact) at $t=0$ ps on the left in blue, indicating the initial state. This initial excitation evolves by dispersion to reach a similar structure on other BChls or permutations. Red zones indicate fewer likelihood distributions where that initial excitation shape could be located, while the bluest or greenest ones show higher likelihood. Again, a speedup can be observed for the largest ${\lambda }_k$ with the lowest T. The similarity between the plots for larger temperatures is again observed. Upper green or blue bands denote the main state transference to certain BChls (3, 4, 7), as indicated in the previous subsection, despite one permutation being optimal, as indicated in Figure 10 and Figure 11. Despite not being identical, this tracing analysis of the initial state distribution provides insight into the more feasible paths being followed by excitation through the set of eight BChls when energy arrives crossing the structure. Note the fast decay of the eighth BChl in the lead of this transition. In addition, it should be noted that for P. aestuarii (reddest) the final state transference is higher, which is preserved through fewer BChl combinations (green or blue) than for C. tepidum (greenest/yellowest), where the exit similarity is partial (in agreement with Figure 10 and Figure 11). Nevertheless, this better-defined final transference is dependent on the temperature (lower) and more noticeable for ${\lambda }_k=65{\mathrm{cm}}^{-1}$. This analysis suggests that the tunnelling is clearly limited by the temperature increase. Nonetheless, this analysis shows that for lower temperatures the tunnelling phenomenon appears strongly if ${\lambda }_k$ increases, otherwise a dispersion effect appears.

3. Materials and Methods

3.1. Data Collection

Experimental data introduced in the analysis are achieved through the Hamiltonians reported in Appendix A for both species, P. aestuarii and C. tepidum. Those Hamiltonians were obtained by combining theory and experiment in spectroscopy including the eighth BChls. Other key structural parameters involved as $r_{\mathrm{trap}}=1{\mathrm{ps}}^{-1},{\lambda }_k=35,65{\mathrm{cm}}^{-1}$ and ${\gamma }_k=50{\mathrm{cm}}^{-1}$ are considered departing from the best spectroscopy observations available. Such experimental information is introduced in the HEOM modelling to reach a time-dependent depiction of the density matrix $\rho \left(t\right)$ for the BChls system. Using spectroscopy units implied that energy is reported in ${\mathrm{cm}}^{-1}$ using the conversion factor $200\pi \hslash c$. In addition, the conversion between ${\mathrm{cm}}^{-1}$ and ${\mathrm{s}}^{-1}$ is reached through the factor $200\pi c$.

3.2. Time-Dependent Density Matrix

HEOM method was implemented under the supervector-superoperator procedure depicted in Section 2.1.4. Implementation of HEOM method was programmed using Mathematica. The basic simulations to reach $\rho \left(t\right)$ with HEOM, a depth $D=4$ was used with processing power in the range of 3.4 GHz and implementing parallel computing with twelve processing cores. Each computer process with a fixed set of parameters took around 10 h considering time simulations of $15\mathrm{ps}$ to reach the equilibrium of populations in ${\rho }_S$.

Using the last HEOM implementation with $D=4$, together with the parameters $r_{\mathrm{trap}}=1{\mathrm{ps}}^{-1},{\gamma }_k=50{\mathrm{cm}}^{-1}$, and the pair of values ${\lambda }_k=35,65{\mathrm{cm}}^{-1}$ ($k=1,...,8$), a set of computer simulations between $t\in [0,15]\mathrm{ps}$ were performed for both species, P. aestuarii and C. tepidum. Such a set was obtained by considering a wide group of representative temperatures: $T=77,131,185,239,293,347\mathrm{K}$. In addition, the overall analysis considered the use of $\rho \left(0\right)={\rho }_{\mathrm{FRET}}^8$ as the initial condition. Each simulation was stored as an individual file of complex-valued data regarding each entry of $\rho \left(t\right)$ for the times considered in $t\in [0,15]\mathrm{ps}$. The numerical time-step for calculations being used was $\delta t={10}^{-4}\mathrm{ps}$, despite, only data with a time-step of $\Delta t={10}^{-3}\mathrm{ps}$ were stored. Those data synthetically and partially represented in Figure 4 and Figure 5 (only the populations ${\rho }_{ii}\left(t\right)$ are plotted there) became the input for further analyses.

3.3. Entanglement Evolution Analysis

Regarding the $\rho \left(t\right)$ values stored, several entanglement measures were obtained for each corresponding time. For $E_{A|B}^{\mathrm{N}}\left[\rho \right]$ (and then for $E_{E^{\mathrm{N}}}^{\mathrm{WAE}}\left[\rho \right]$) and $E^{\mathrm{MW}}\left[\rho \right]$ own Mathematica programming procedures were constructed regarding the overall bipartitions needed. Those procedures let to reach the last measures for the corresponding discrete times. Similarly, to reach each $\mathsf{\Lambda }\left({\rho }_{jk\ell }\left(t\right)\right)$ considered in the analysis, the corresponding Uhlmann’s procedure [73] was programmed using Matlab.

3.4. State Transference Evolution Analysis

State transference analysis has required to transform each $\rho \left(t\right)$ stored by exchanging the order of BChls in agreement with each element of the set $\mathsf{\Pi }$. The procedure was directly constructed in Mathematica to then get ${\mathcal{F}}_{{\pi }_i}\left(t\right)$ for each element in the set, and in each discrete time for our original data in $\rho \left(t\right)$. Finally, the maximum value, ${\mathcal{F}}_{\mathsf{\Pi }}\left(t\right)$ in each time is easily obtained to become represented in Figure 10 and Figure 11. Instead, Figure 12 and Figure 13 are complete representations of each ${\mathcal{F}}_{{\pi }_i}\left(t\right)$ through the time.

4. Conclusions

Quantum features are undoubtedly present in living systems. We are still learning how these features set important aspects of life survival. With even higher complexity than chemical systems, the task is very intricate because it needs to pass through several size orders. In such circumstances, computational biology is advisable as an important element gathering efforts from several disciplines, including physics, chemistry, mathematics, and computer science, leading to specialized computer simulations that can describe the meaningfulness of these features for life.

For FMO, energy transition through LHC exhibits clear quantum features such as entanglement among BChls, together with state tunnelling crossing the set of BChls to reach the RC. The analysis of the behaviour of these quantum features requires experimental work developed in the spectroscopy domain as well as computer modelling to perform simulations of the complex processes involved in their dynamics. Such combined approaches are settled in the terrain of computational biology to obtain insight into relevant features. The current work extends several traditional approaches to analyse the quantum features of entanglement and state transference.

By introducing alternative or more elucidated entanglement measures for the set of eight subsystems, a clearer understanding of the phenomenon is pursued. Note that due to the current limitations of entanglement measures for larger mixed states, an approach based on bipartitions becomes useful; however, each type of measure provides different information. In addition, authentic multipartite entanglement beyond bipartite becomes present, as our analysis based on the fidelity with respect to the closest separable state has shown.

In this sense, multipartite entanglement must be understood from different viewpoints beyond the bipartite approach. The complexity of this quantum phenomenon is reflected in the fact that several entanglement measures have been proposed, each one considering particular aspects of non-locality. In the present work, we have considered the maximum overlap with respect to the closest separable state as a natural entanglement measure, as it provides information concerning the likeness of a given mixed state with one that is fully separable. We focused on two particular cases: the three BChls with the largest populations, and BChls 3, 4, and 8, which are related to the main BChl exit and entry. Whereas most of the measures used in the literature are computed as the averages of bipartite measures, this function is able to quantify the entanglement of subsystems compounded with three or more qubits. In this context, we considered partitions of three chromophores and computed $\mathsf{\Lambda }\left[\rho \right]$ using a numerical algorithm [73].

Nevertheless, the same study could be accomplished for any number of BChls traced out the general state. Our numerical simulations showed that the state of the system is close to being in a fully-separable state at most times, even though there are cases where a non-vanishing degree of entanglement is achieved. Our study represents a step forward in the analysis of the multipartite entanglement of three or more BChl subsystems.

For state transference, models based on the tracking of the initial state distribution through BChls have allowed us to identify a process with fidelities above $60\%$ for the worst cases, sometimes reaching nearly $80\%$. In particular, the approach taken in this work to identifying the initial state in all possible permutations of the site basis elements has allowed us to identify its transition through the evolved state in this complex feature, again strongly supported by a computational approach regarding a post-experimental simulation analysis on the biological FMO system.

Superposition, entanglement, and tunnelling are key aspects of the quantum world. Recently, it has been noticed that entangled systems enlarge the coherence time needed for quantum process [71,77]. Thus, by extending the entanglement, it appears that coherence becomes preserved over larger times, which provides the possibility of computing the energy transference more accurately. In such research, each time larger system sizes are used, implementing an atom in a superposition of hundreds of locations, the quantum coherence becomes extended on the scale of seconds. If a similar mechanism happens in living systems the related quantum effects could play a key role in daily survival.