Parametric Mapping of Quantum Regime in Fenna–Matthews–Olson Light-Harvesting Complexes: A Synthetic Review of Models, Methods and Approaches

: Developments in ultrafast-spectroscopy techniques have revealed notably long-lived quantum coherence between electronic states in Fenna–Matthews–Olson complex bacteriochlorophylls, a group of molecules setting a nanoscale structure responsible of the coherent energy transfer in the photosynthetic process of green sulfur bacteria. Despite the experimental advances, such a task should normally be complemented with physical computer simulations to understand its complexity. Several methods have been explored to model this quantum phenomenon, mainly using the quantum open systems theory as a first approach. The traditional methods used in this approach do not take into account the memory effects of the surroundings, which is commonly approximated as a phonon bath on thermal equilibrium. To surpass such an approximation, this article applies the Hierarchical Equations of Motion method, a non-markovian approach also used to analyze the dynamic of such a complex, for the modeling of the system evolution. We perform a parametric analysis about some physical features in the quantum regime involved during the quantum excitation process in order to get a comprehension about its non-trivial dependence on operation parameters. Thus, the analysis is conducted in terms of some relevant physical parameters in the system to track the complex global behavior in aspects as coherence, entanglement, decoherence times, transference times, and efficiency of the main process of energy capturing. As a complementary analysis from the derived outcomes, we compare those features for two different species as a suggestive possible roadmap to track genetic differences in the photosynthetic performance of the complex through its biological nature.


Introduction
Photosynthesis is a fundamental process for life on Earth through which light energy coming from the sun is absorbed and used by photosynthetic organisms to perform necessary metabolic reactions. Its success is the result of millions of years of evolution during which those organisms have developed highly specialized structures able to capture, transfer and store light energy by a series of physical and chemical transformations. An important group of these structures are the light-harvesting complexes (LHC), responsible for capturing photons and transferring their energy to the reaction centers (RC), where charge separation takes place. Since their early discovery, LHCs have captured the attention of researchers due to their remarkable nearly 100% photochemical quantum yield [1][2][3][4].
The study of photosynthesis continued as theories on electron transfer in biological systems matured [5,6], and the Fenna-Matthews-Olson (FMO) complex took the spotlight as it became the first mapped LHC structure [7,8]. Multiple disciplines are actively collaborating in the race for unraveling the mechanism behind this energy transfer efficiency and its large range of possible applications that range from photovoltaics to quantum computing [9,10]. New and improved information about the structure and its behaviour continues to emerge [1,11,12]. While theoretical physics has provided the basis for advanced femtosecond laser techniques necessary for tracing excitation energy transfer (EET) dynamics [13,14], biochemical studies on the species [15,16] and genetic methods are being used for purifying, characterizing and even modifying the protein structure [17][18][19][20]. Additionally, improving X-ray crystallographic procedures have been key for obtaining high resolution electron density maps of protein complexes which allow insight into existing interactions within the macro molecule [7,8,21].
Though the presence of long-lived quantum coherence in the FMO complex has been confirmed [22][23][24], its role and relevance for the high energy transfer efficiency properly in the central photosynthetic processes within the protein is still a matter of debate [25,26]. The conflict originates from the comparison between coherent (e.g., Lindblad equation) and incoherent (e.g., Förster theory) theories which appear to yield similar approximations since coherence decay occurs on different timescale as the energy transfer, thus questioning the actual importance of quantum coherence [21,25]. Many authors have stated that the observed efficiency in energy transfer within light-harvesting complexes such as the FMO require both coherent and incoherent energy transfer [23,[27][28][29].
Although the presence of environmental noise and the acceleration of decoherence may seem counterproductive, this transition might be necessary to assure fast unidirectional transfer towards the reaction center [23,30,31]. Going further, theoretical physics and chemistry are developing a better understanding of the mechanisms involved from a quantum perspective that enables their simulation and allows tracking of possible control parameters still for further technological applications exploding the quantum features on a biological substrate. This may lead to the hopeful application of this system as a biological resource of quantum processing, particularly due to its room temperature performance [22,30,32].
This merge of disciplines has resulted in the introduction of physical concepts into a biological context for the wide study of photosynthesis. For instance, in solid-state physics the exciton refers to a collective electronic excitation with well-defined quantum properties in a molecular crystal. Though lacking the same ordered and static structure, photosynthetic proteins also share excitation among their pigment molecules in a more ordered manner than pure pigments in solution [21]. Thus, these systems have been settled to be somewhere between crystals and fluids [21]. The complexity of the system has required the use of well-known developments of other disciplines as tentative approaches to the problem. A diagram of the main scientific branches involved in the study of the FMO complex is depicted in Figure 1 remarking just some emblematic works contributing to the understanding of the FMO complex. This map shows the complexity and richness of the problem which should be addressed in the understanding and collaborative research around it.
In fact, the first widely studied light-harvesting complex was the FMO pigment-protein. This photosynthetic antenna holds bacteriochlophyll a photopigments within its protein scaffold. Bacteriochlorophyll a (commonly abbreviated as BChl a but not widely used in this work) is only one type of the different photopigments (also known as chromophores) being present in the majority of anoxygenic photobacteria [4]. These bacteriochlorophylls act as antennas capturing photons and transferring their energy through the complex to the bacteria reaction center. For the reader's ease, in the remaining development of this article, we will be using only the abbreviation BChl as reference to any bacteriochlorophyll a chromophores in the complex. The FMO complex is found in all species of green sulfur bacteria (see Section 2.1). The study of the FMO complex is a clear example of a multidisciplinary effort to understand its entire behavior departing from its biological role, chemical structure and physical processes involved, among them, the quantum properties of such a biological system. Excitation energy transfer in GSB within the FMO and including other photosynthetic complexes can be pictured as a multi-layer cascade process for transfer of energy first coming from the environmental photons ( Figure 2). There, the FMO could be conceived as a multipartite funnel with a spring scaffolding where energy is collected at a top large surface to excite some of its components (the pigments) and then being transferred and concentrated towards a sink, while the process is damped by its scaffolding structure to reboot the FMO. The energy transfer happening in this funnel involves the migration of excited states being sequentially shared by the pigments. In order to bring this analogy closer to reality we must think of this light-harvesting funnel as an arrangement of pigment particles (the BChls) with dipole like charge distributions and springs holding them to the surrounding structure and connecting them to one another. The photons (energy packets) from which the energy is being absorbed and transferred can be visualized as a rain of pebbles from different sizes (their energy or frequency). To know how many photons are absorbed by a pigment particle one must only know the photon flux and the effective cross-sectional area of the pigment [4]. This is not a physical size but rather an effective size considering factors such as the wavelength of excitation. As energy transfer occurs in this system, we can imagine the springs movement causing slight changes to the original arrangement at the time they partially damp the process. It can be understood as the changes in the protein structure as the transfer occurs, which we label as the "reorganization energy", λ. The driving force for energy transfer is more or less equal but in opposite direction to the reorganization energy [21]. Key quantities are the electronic coupling between the BChls, the electron-phonon coupling (reorganization energy λ k and time scale γ −1 k ), the temperature T, and the disorder [21]. Finally, such an energy transfer boosts the reduction/oxidation chemical reactions taking place in the RC, involving other biological components and continuing the photosynthetic process. This loony picture ( Figure 2) is useful to capture the general process formally depicted below in a much more technical description. Experimental work has addressed those problems related with the understanding of the chemical structure of this nano-structure. In the quantum domain, two-dimensional spectroscopy has been able to capture the quantum beating of its inner components in terms of strengths and times of absorption-emission in order to quantify its quantum energy levels for the entire structure whose main interactions are electrical due to the dipole momentum of the BChls and mechanical between BChls and the protein scaffolding structure. Despite the fact that analysis is not in vivo and then developed to very low temperatures, it means under very different operation conditions that still this analysis together with theoretical modeling allows to advise the quantum mechanical properties of its structure. Then, such a structure can be computer-modeled to extrapolate its behavior under realistic operation conditions. That is the reason to improve the quality and the extent of such computer simulations for the FMO based on more faithful models.
The aim of this article is three-fold. Firstly, to synthetically present a road map of well-stated methods to understand and to simulate the FMO light-harvesting complex of green sulfur bacteria since a quantum point of view. Such comprehension could set a fingerprint to characterize biochemical structures being present in nature while they still exhibit genetic differentiation. In addition, the immersion into the FMO simulation research could be complex and confuse in spite of the large amounts and diversity of the literature written in the last two decades. Thus, the article sets a certain roadmap transiting the main stages in their study and simulation. In other trend in parallel, our main goal is to develop a parametric analysis of the quantum nature of the complex which has been partially attempted in some previous works as a secondary task in them. We analyze here a more extensive numeric and mathematical analysis of the non-monotonic and multivariable quantum behavior for the FMO complex as function of certain characteristic structural parameters of operation. Such methods have a feedback from several theoretical and experimental techniques ranging from biochemistry to spectroscopy to genetic mapping, performing a multidisciplinary approach in their analysis. Through this multidisciplinary approach, we set the most common setup for the simulation analysis moving forward into discussion, focusing and deepening on some novel approaches in terms of the parametric dependence of certain generic quantum markers associated with the efficiency, long-lived coherence and quantum entanglement sustenance in the complex. Thus, in parallel, the article develops and presents an extensive analysis of the parametric dependence on some of those structural properties which has been barely determined through experimental analysis, conforming the entire quantum regime during the excitation. A third exploratory aspect, in spite the parametric analysis, is the relation among the predictable behavior depending on such structural parameters with the genetic study or inclusively the manipulation through similar species or related strains. We devote certain room to discuss this interesting issue as an epilogue. Thus, the second section presents a structural description of FMO complexes and their main elements necessary for its modeling as well as an introductory description of their main mechanisms to boost quantum excitation and light-harvesting. The third section presents the generalities of several approaches in the modeling of quantum interactions inside of each monomer of the FMO complex from data reported for the Chlorobaculum tepidum bacteria, treating the protein structure as an open quantum system, particularly through Redfield and Lindblad markovian equations. This section is closed with the presentation of the Hierarchical Equation of Motion (HEOM) method which extends the most traditional equations to model open quantum systems into a non-markovian approach. This method will be used for the further development of analysis in the article. The fourth section discusses the simulation details in terms of physical and computational elements to perform the base simulations together with the outcomes through a range of selected physical parameters, setting a multivariable dependence of remarkable properties in the quantum functioning of FMO complex. The section discusses the parametric dependence of coherence, efficiency and transference, together with the entanglement created in the process for the Chlorobaculum tepidum bacteria as instance. Finally, the fifth section discusses and analyzes openly, in the terms of the previous presentation, the extension for other bacteria and strains about possible genetic methods of discrimination based on such parametric characterization on their quantum fingerprint, as well as its implication for their possible genetic manipulation. The final section exposes the conclusions of this work and sets the questions leading towards further research.

Structure of the Fenna-Matthews-Olson Light-Harvesting Complex
This section is intended to familiarize the reader with the biological roots of the problem and its relation to quantum mechanics. It begins by describing the microorganisms where the FMO LHC is found and the conditions under which they grow. Then a technical description of the structural and chemical characteristics of the complex is presented. Finally, the quantum dynamics describing the interactions governing the behaviour of the complex are presented for further sections.

Green Sulfur Bacteria and the FMO Complex
The FMO protein is found in green sulfur bacteria (GSB) and chloroacidobacteria (excluding filamentous anoxygenic phototrophs) that contains a Type I reaction center [4]. Sulfur bacteria have evolved to use sulfur in higher amounts as an electron donor or electron acceptor for energy production, taking advantage of its properties as redox reaction mediator [15]. These prokaryotic organisms belong to the Chlorobiceae family. See Figure 3 including the excellent microscopic approaches reported in [16]. Photosynthetic sulfur bacteria are divided in green sulfur bacteria and purple sulfur bacteria due to the photosynthetic pigments being present in each one. Their distribution and growth are determined by the light availability and sulfide concentrations in their media [15]. GSB have been found in a variety of environments ranging from sulfur springs and deep-sea hydrothermal vents to the anoxic hypolimnia of lakes, where they grow under either high or low but stable sulfide concentrations [15,16]. The depths at which these bacteria are found, where there is little if any access to light, and the vast range of temperatures at which they grow and adapt, contributes to the possibility that the efficiency of their photosynthetic mechanisms take advantage of quantum effects at relatively high temperatures. Green sulfur bacteria are anoxygenic photosynthetic sulfur bacteria that use light as energy source for carbon fixation, mainly using hydrogen sulfide (H 2 S) as electron donor [15]. These bacteria are obligate anaerobic photolithotrophs that only grow under strict anoxic conditions using (depending on species) sulfide, elemental sulfur, thiosulfate, molecular hydrogen or even reduced iron and other organic substances as electron donors for anoxygenic photosynthesis [15]. The electrons from the reduced form of sulfur are used for CO 2 fixation via the reverse tricarboxylic acid cycle, while the oxidation of sulfide results in the formation of sulfur globules deposited outside the cell [4,15]. GSB have adapted to survive the low light intensity conditions of their habitats [15].
Initially known as the "bacteriochlorophyll a protein", the Fenna-Matthews-Olson (FMO) protein gets its name from Roger Fenna and Brian Matthews, who first determined its structure [7,8], and John Olson, who discovered the protein [33]. This LHC is located between the chlorosome (another LHC) and the reaction center as shown in the representation included in Figure 4. It participates in the excitation energy transfer process during photosynthesis by funneling that excitation from the baseplate of the chlorosome to the reaction center [4]. This is possible due to its pigment-protein couplings which tune the optical properties of the complex [1]. . Schematic representation of the spatial arrangement of the Chlorosome, the FMO protein and the reaction center in the cytoplasmic membrane. Light is absorbed by photopigments in the chlorosome and its excitation energy is transferred down to the baseplate which then transfers it to the FMO complex and sequentially to the reaction center.

Structure of the FMO Complex
The FMO is a water-soluble protein with a molecular weight of 150,000 Da [7,8] with a maximum diameter around of 8.3 nm. This rare property has facilitated its crystallization for high resolution spectroscopic studies [8,12]. Trimer structure of FMO includes monomer sequences working together still with certain independence to gather the light excitation. Each monomer is a ∼360 amino acid sequence of the complex which folds into a "bag" containing eight bacteriochlorophyll a molecules with an average distance between nearest-neighbor pigments of 12 Å [17] (Figure 4). The FMO structure, shown in Figure 5, consists of a homotrimeric protein complex, which subunits are related by a 3-fold symmetry axis. The side chains of the amino acids conforming the protein structure interact with one another through Van der Waals forces, dipole-dipole interactions, ions and hydrogen bonding creating structures called α-helices, β-sheets and random coils resulting in the monomer shown in Figure 5B. The β-sheet ribbons compose the large surface wall of the complex, and its amphipathic nature (possessing both hydrophilic and hydrophobic properties) provides shielding for the non-polar chlorophyll core from the surroundings, while α-helices and random coils make up the contact region between subunits, Figure 5A. These nano-structures have an important effect in BChl interactions [14]. It was initially reported that each of its subunits contains seven bacteriochlorophylls [7], but an eighth BChl molecule was later discovered between subunits, bringing the total number of BChls in the protein to 24 [1,11,12,34]. The established nomenclature for numbering the pigments is still the originally proposed by Fenna and Matthews [7], to which the eighth BChl was added. Following this nomenclature, BChl 3 and BChl 4 are the nearest to the cytoplasmic membrane while BChl 1, BChl 6 and BChl 8 are closer to the chlorosome baseplate (as seen in Figure 5) [1,4,34,35].
The reader must remember that the FMO complex is present in all species of GSB, being C. tepidum and P. aestuarii the two most widely studied cases. Genetic diversity between species and strains (genetic variants of the species) becomes in slight structural differences of their FMO complex. To fix the analysis, the model calculations performed in this article are based on some experimental reported values for C. tepidum [13] despite we tackle a compared view for P. aestuarii in the closing section. Section 6 will discuss a deeper analysis on the effects of these structural differences and provide a comparison between FMO complexes among certain species and strains.

Interactions within the FMO Complexes and Its Environment
The folding of the protein will result in the formation of local environments in which the amino acids may be subject to different conditions thus modifying their behaviour. These localized environments influence the optical transition energies of the pigments [1]. These interactions are known as pigment-protein interactions, and their resulting influence on the whole quantum site energies in the FMO complex are schematized in Figure 6 without a detailed quantum description which will be addressed in the development below. The position and orientation of the chromophores is highly related to their excitation energy and their role in the energy transfer path, shifting the absorption energy of the pigments facing the outer antenna towards blue compared to those linked to the reaction center shifted towards red [1].
Determining the site energies of each BChl molecule in the complex protein structure can prove to be very challenging. For this reason, site energies are often treated as parameters that are determined from the fit of optical spectra [13]. Experimental spectroscopy is able to discriminate emissions of excited BChls in short time window intervals, but they are still limited to some few femto-seconds. In addition, such two-dimensional spectroscopy experimental techniques require working at very low temperatures (typically around 77 • K) thus detecting quantum beats in the BChls. In this sense, spectroscopy has been successful to get, together with theoretical modeling, a concrete quantification of the Hamiltonian corresponding to dipole-dipole interactions and some limited deviations due to the protein scaffolding of the set of BChls inside FMO complex [36]. From a remarkable interest is the novel inclusion of the eighth BChl in the experimental and FMO computer modeling literature [1,37]. Those Hamiltonians will be discussed in the following subsections and reported in Appendix A.
Those main experimentally obtained energies show that the excited BChl 8 and its strongly excitonic coupled BChl 1 are the two main entrance points for excitation energy coming from the baseplate [1]. BChl 6 provides an alternative route for the excitation energy to the reaction center [38]. Thus, for the studies based only on the first seven BChls, authors commonly consider two pathways for the initial excitation energy: BChl 1 and 6. The process of relaxation of the excitation energy is funnelled to BChl 3 along the branches formed by BChl 2 and by BChl 4 to BChl 7 as seen in Figure 6 [1]. The latter being completed within 500 fs, while BChl 2 appears to limit the overall relaxation/equilibration time to about 1.5 ps [1]. In any case, the short relaxation times and the large energy gap between the entrance BChl 8 and the sink BChl 3 prevents recombination and leaking [1].
Due to their structural complexity and behaviour, proteins are systems with thousands of degrees of freedom [39]. On the other hand, the functional subsystem (the active site where the BChl cofactors are bound) involves only a few quantum states [39]. The change in the quantum state as the transition from the ground state to the excited electronic state is associated with a change in the electric dipole momentum of the subsystem [39]. The polar residues contained in the protein and its highly polar solvent (water) surroundings result in a strong interaction between the functional subsystem and its environment. Protein itself may undergo structural and electrostatic changes depending on its environmental conditions, meaning that the environment must be included in the model for a correct approximation to the system [1,14,39].

Quantum Modelling of FMO Quantum Excitation
This section states the basic quantum modeling principles and the most relevant approaches that have been applied to model the FMO behavior. The interactions within the complex itself and its environment are described together with the main assumptions considered in the beginning of the section. The Hamiltonian describing the system dynamics is introduced for further treatment.
Concept of quantum open systems is established considering the environment contribution to the system through quantum master equations (as Lindblad and Redfield ones) as the first attempts to solve the problem. This section then deepens into the hierarchical equations of motion (HEOM) method as a more accurate approximation together with the superoperator approach to the computer solving of master equations. Finally, the effects of varying parameters of the model on coherence and entanglement are discussed through the section.
The quantum mechanics approach to the modeling of FMO complexes is first given through a sufficiently accurate Hamiltonian H S reproducing the main interactions in the complex. It means, those dipole-dipole electrical components among BChls considered as system S. This main Hamiltonian is built by the individual excitation energies along the diagonal and the dipole-dipole coupling terms in the off-diagonal positions. Diagonalization of this Hamiltonian yields to a set of eigenenergies corresponding to the eigenstates in the exciton basis [21]. The energies of these excitons correspond to the observable transitions in a linear absorption spectrum, nevertheless, due to homogeneous and inhomogeneous broadening the energies of the entire seven or eight excitons have never been experimentally measured, although an enormous effort in theoretical work has been done to model the system [35], departing from the main interactions obtained by spectroscopy in order to be predicted. Thus, this interacting set of BChls states a quantum system with well-defined energy levels.

Main Dipole Interactions among BChls
Changes in the absorption spectrum of chromophores are usually different for excitonic coupled chromophores compared to the isolated ones [21]. Both in the ground and excited state, the molecule has a permanent dipole momentum, µ g and µ e , respectively, which can differ significantly for asymmetric pigments like bacteriochlorophylls [21]. In the presence of a linear electric field F E , the energy levels E g and E e for the ground and excited states change (together with the difference ∆E between them). Those changes are: where ∆E = E e − E g , ∆E 0 = E 0 e − E 0 g , ∆µ = µ e − µ g , and θ is the angle between F E and ∆µ. The superscript 0 refers to the case F E = 0. We avoid the report of the spatial form of fields and forces F E for the dipole-dipole interactions because it is common in the electromagnetism literature [40]. The main aspect is highlight that the knowledge about the FMO structure is essential to understand the energy transfer dynamics and their performance between chromophores. Nevertheless, the protein itself may undergo structural and electrostatic changes depending on its environment conditions. Then, to understand the EET dynamics within the molecule, one must first to understand the coherent distribution of excitation among the different BChl sites. This is usually achieved taking two assumptions into consideration: (1) BChl sites are modeled as two-level systems, and (2) only one site can be excited at a time for EET to occur, making negligible the probabilities of bi-excitonic and other higher states compared to the single excitation [41]. Those assumptions can be attributed to the dipole blockade, an effect in which excitation of one site shifts other sites out of resonance due to the interaction energy being added to or subtracted from, the excitation energy for attractive and repulsive interactions, respectively [41,42]. The strongest forces are exerted by the electric dipole interaction, arising from the closely spaced arrangement of the BChl sites [41,42]. Since only exciton temporal dynamics have been experimentally observed through spectroscopy techniques, site excitation dynamics still must be indirectly modeled [35,41,[43][44][45] in order to get an entire description of the Hamiltonian.

Construction of the Entire Hamiltonian to Model FMO Complex
Because of the experimental analysis, each BChl i becomes in two most common excitation states denoted by |0 i (the ground state) and |1 i (the first excited state). Such most common states for each single BChl set a two level restricted system for each BChl. Together, the tensor product of them could set a basis for the entire set of N BChls (N = 7, 8). This notable property has raised the attention of quantum information science on the complex. Nevertheless, as it has been proved [41], the most common global state corresponds to the superposition of single excited BChls, |k ≡ |0 1 0 2 ...1 k ...0 N , called the occupation or site basis. Then, H S is written as: where | i are the eigenstates of energy in the BChls system, the excitonic basis. Those Hamiltonians have been reported in the literature for N = 7 [13] and N = 8 [1], both are included in Appendix A. While, the structure surrendering such BChls, the protein structure scaffolding such chemical structures ( Figure 5B) has a lower level of interactions but still meaningful on the BChls. For instance, as we will discuss below, characteristic wavelengths exciting the set of BChls through dipole-dipole interactions represent energies of order of 1.2 × 10 4 cm −1 which are compared with the reorganization corrections due to the interactions with the protein scaffolding in the order of just 35∼65 cm −1 (using the common spectroscopy units of energy in cm −1 discussed in the following subsections). Thus, it is commonly approached as a phononic medium. It means quantum vibratory modes on the mechanical structure exchange energy with the set of BChls. Due to the complexity of such a scaffolding structure, it is normally modeled as a continuum medium of oscillators working additionally as a thermal bath, B, which barely reproduces the whole phenomena present there. The Hamiltonian for the bath H B is then considered as a set of oscillators of mass m α for each BChl. In terms of their positions q k,α and momenta p k,α [46], their global Hamiltonian is: which is written in terms of the modes ω α of such oscillators. Each oscillator has the possibility to switch the BChl state |k in H S . Due to the interaction with each BChl, such oscillators become displaced into a new equilibrium position q k,α = q k,α − d k,α . Thus, under such an interaction: where the following relations were considered: m α c 2 =hω α , c = ω α λ α (for the mass, the frequency and the wavelength of each (k, α) phononic mode). Together, d k,α = d k,α /λ α is the dimensionless displacement from each equilibrium configuration. Then, the first term represents the bath energy. The second term is a reorganization Hamiltonian modifying the diagonal elements of H S . Third term is the interaction between the system and the bath (with a linear response). By defining S k = |k k| and B k = − ∑ αh ω α d k,α q k,α , we arrive at: Be careful with the difference between λ α , the phonon wavelenght, and λ k , the reorganization energy for each BChl site excited state which is mainly related with the resilience of each BChl under the bath interaction. Note the bilinear form for H S−B as in (9), where S i , B i are operators stating a basis for each subsystem S and B, respectively (particularly S i is an orthonormal basis of N 2 matrices with S N 2 the only with trace different from zero in {S i |i = 1, 2, ..., N 2 }). Thus, the complete Hamiltonian for each monomer in the FMO complex becomes (considering H S → H S ⊗ 1 B ): Then, this Hamiltonian is not separable, thus able to generate entanglement among the BChls and the bath together. Hamiltonian is intended to obey for the entire density matrix ρ T (including the BChls system S plus the thermal bath S + B) the Schrödinger-Liouville equation: by switching into the interaction picture between the system and the thermal bath (H S−B ): the Schrödinger-Liouville equation becomes in this picture:ρ

LHCs As Open Quantum Systems
Master equations of open quantum systems are commonly obtained by tracing the dynamics of the thermal bath system approximating thus the behaviour of the main system in the Schrödinger-Liouville equation. Integrating (13):ρ (14) and then substituting last equation in (13), we obtain:ρ Tracing the bath and assuming as approximationρ S = Tr B (ρ T ) (barely valid only if both systems become weakly entangled), which it is additionally based on certain approximate stability for the bath: we get by using the commutator properties:ρ (16) which under the Born-Markov assumption stating thatρ S (t) only depends on the current state of ρ S (t), it becomes:˙ρ it is a simplified departing point to derive the Lindblad and Redfield master equations (brief derivations of them are included in Appendix B and Appendix C, respectively). In the next sections we report some interesting and successful efforts to use those quantum master equations in the understanding of the FMO complex simulation.

Lindblad Approach to FMO Complexes Modelling
Lindblad master equation in its diagonalized operators form: represents the interaction of one quantum system with an external bath through the generic operators S j or L α as algebraic basis of the systems (see Appendix B). For the FMO complexes, despite L α still comprising the physical information, it becomes limited due to the impossibility to reproduce faithfully more details of the system and bath interaction. A representative proposal performed by [47] in order to adjust this master equation to the FMO problem has been described by a reduced density matrix ρ for the system of BChls. By including terms in the Schrödinger-Liouville equation for the decoherence due to environmental interaction, recombination of excitation on a common ground state, and transfer of excitation to a sink (explicitly representing the reaction center): captures the basic environmental effects used in several previous studies [28,30,31]. By introducing directly some operators in the site basis, it has been shown that the description of the Lindblad equation accounts for the main features of the dynamics [30,36]: where a site-independent dephasing rate is given by γ and {, } is the anticommutator. Similarly, an irreversible transfer of excitation from the m−th chromophore to the sink |s (an additional subsystem in the modeling, while m commonly will come m = 3, 4 BChls) is represented by: where κ denotes the sink rate. The irreversible loss of excitation due to recombination is given by an analogous term: with a site-independent recombination rate Γ with the ground state (represented by |0 ). Note that the need for an additional sink |s and the ground state |0 should include two additionally dimensions giving a total of N + 2 sites in that model [47]. Note that each additional term in the modeling fits to the form of the generic term 1 The relevant environmental parameters are almost universal in the related literature: the dephasing strength γ, recombination rate Γ and the sink rate κ. Standard values are inferred from experiments, the sink rate κ = 1 ps −1 , the recombination rate Γ = 1 ns −1 and the dephasing rate at room temperature γ = 300 cm −1 given in their traditional spectroscopy unities transforming between cm −1 to s −1 via the factor 200πc or between cm −1 to J through the factor 200πhc. The dephasing rate, being a product of the temperature and the derivative of the spectral density, can be estimated by using the experimentally determined parameters of the spectral density [13,47]. This value approximately agrees with the optimal dephasing rate at which transfer is most efficient [28,31].

Redfield Approach to FMO Complexes Modelling
The Redfield master equation:ρ with Σ i and Λ i comprising the physical information (see Appendix B) of coupling between the system and the environment. This master equation has been used successfully to model physically the quantum FMO dynamics [36,43] and adaptations of it has been used to analyze and to propose notable quantum processing applications based on FMO complexes [10]. A concrete implementation to such a system is given explicitly by [48] in the form of the Bloch-Redfield equation: where V † H S V = diag({ i }) and |a n can be any basis of states of the quantum system. If particularly the eigenstates of the system Hamiltonian H S are chosen as basis |a n = | n , then the matrices V become a unitary operator. There, C jk (ω) is the spectral function of the environment defining the noise spectrum as well as the strength of spatial correlations q jk ,q jk between sites j and k. The system operators S j define which part of the system couples to the noise environment.

HEOM Method Approach in FMO Complexes Models
The Hierarchical Equations of Motion (HEOM) method [49] is a non-markovian approach including the features of previous models, thus considered better than the previous ones which are markovian. It is considered and constructed as D recursive equations thus considering D previous temporal stages for the memory of the bath. Then, such stages are labeled each one with a vector n of dimension N considering each chromophore and its interaction with the bath. Thus, only ρ n with n = (0, ..., 0) is the real density matrix of the system ρ S while other ρ n are just auxiliary matrices regarding the memory. This set of matrices have an associated degree s = 0, 1, ..., D corresponding to all vectors n = (n 1 , n 2 , ..., n N ) with 0 ≤ n k ∈ Z + ∪ {0} such as ∑ N k=1 n k = s. Thus, the following method gathers the main features than its markovian relatives, Redfield or Lindblad quantum master equations. In fact, it includes the main Hamiltonian H S , the relocalization correction, the trapping from BChls 3 and 4, together with other recursive non-markovian terms obtained from its derivation [49] which includes a treatment for the equilibrium: where k is the Boltzman constant and β = kT. Together, V k = |k k| and ∆ k = λ k βh 2 γ k (note particularly the reorganization term). In addition, if n is a vector of order s, then n k± is the vector of order s ± 1 similar to n with only their component n k increased (decreased) by one. γ k is the bath cutoff frequency reflecting the non-markovian nature of the bath, the rate that information flows from the system to the bath [50]. Thus, γ k comprises the interaction strengths between the bath and each BChl. The initial condition for such a system is the corresponding to ρ S at t = 0 being unentangled with the bath (typically ρ S (0) = |k k| with k = 1, 6 if N = 7 or k = 8 if N = 8 in agreement with the behavior of the set of BChls in each monomer). Last model has been considered by several authors [38,51] for N = 7. Because it has been observed that BChls 1, 6 and 8, works as FMO antennas while BChls 3 and 4 drive the energy oscillations to the RC at the trapping rate, r trap . Thus, the first ones are considered as initial values in the equations and the second ones within a term in (28). While, other ρ n (0) for n = 0 are typically settled as the zero matrix. In the same way, if n = 0 is considered in (28), note that n k− does not appear in the equation because the coefficient n k = 0. Otherwise, if n is of order D then we will take n k+ = 0 as cut-off. At this point, the main models to afford the computer simulations for the FMO have been considered to arrive in the central master equation used in our further development: HEOM method (still valid in the non-markovian regime). Those methods are used to supply complementary aspects not yet affordable with experimental techniques in order to extend their operative behavior under natural conditions not directly observed in the experimental realm. A remarkable aspect is that methods as Redfield and still Lindblad (or inclusively the Secular approximation not considered here [36]) have reported simulation outcomes completely in agreement with experimental spectroscopy studies in terms of the excitation times, the unchained BChls excitation ordering, together with the prediction of population rates for the site occupation for the main BChls [36,43,47,48]. Those computer simulation techniques show oscillations in the dynamics of the exciton populations persisting even at temperatures over T = 300 • K.

Superoperator Version of Master Equations
In order to take computational advantage, master equations as (18), (23) and still (28) can be written simpler by switching them into the so called superoperator approach. In fact, by applying the rule: in each term of (18), (23) and (28), we can jump ρ on the right through other operators by switch it into a supervector of size N 2 (if ρ is a N × N matrix). The remaining operators in the expressions become comprised in a superoperator κ (M) (M = L, R or H for Lindblad, Redfield and HEOM, respectively) of size N 2 × N 2 obtained through the tensor product ⊗ as: Note that this single transformation allows to transform any of the three quantum master equations into a linear single matrix differential equation, easy to implement computationally despite its size. Appendix D develops the operative process applied to the HEOM method, thus providing a tool to translate it into a matrix differential equation.

Numerical Solving of HEOM Models and Analysis of Physical Measures to Feature Parametrically the Quantum Regime in the FMO Complex
As it was previously mentioned in Section 2.2 and they are illustrated in Figure 4, there are eight BChl molecules embedded in the protein scaffold of each FMO monomer. Even though this has been known and confirmed since over a decade [1,11,12], most of the publications related to the theoretical calculations of energy dynamics within the FMO have only considered the seven originally identified BChl sites [14,41,43,47]. This is due to the fact that some models were published before the eight BChl was located, but also more recent studies continue considering only seven BChls in their models under the premises that the eighth BChl is weakly bound to each monomer and is relatively far away from the other seven pigments [47]. We will come back to analyze briefly this fact at the end of the article. It is believed that this BChl increases the overall energy-transfer efficiency of the complex acting as a bridge between the baseplate pigments and the core BChls [12]. In this study, we deliberately decided to consider first only the core seven BChls as we aim to compare our analysis with those previously reported. Such tracks consider the BChls 1 or 6 as detonators of dynamics. As we will see, those initial conditions rise two different process, each one corresponding to each side (left-right) of interactions depicted in Figure 6. On other report, we have already considered the eighth BChl in the dynamics following the entire model [52]. Thus, after to present the different approaches and lines of research in the study of FMO complex dynamics in the quantum regime, we arrive at the main topic of our analysis, the parametric behavior of such dynamics. The extensive analysis has been developed departing from the initial condition ρ(0) = |1 1| instead ρ(0) = |6 6| also, because the first is considered as the main entry of excitation energy for many authors. Nevertheless, in the final section, we will come back briefly to this other pathway considering in addition the case of eight BChls.

Methods, Extent of Analysis, and Source Data
The HEOM method reproduces the behavior of one monomer on the FMO complex under the descriptions before depicted. We are considering three different parameters which are important in the description of such behavior to perform a parametric analysis of the quantum features involved in the excitation process. Thus, in our analysis, a series of HEOM simulations were performed by sweeping the parameters of temperature (T), relocation energy (λ) and dephasing rate (γ) to simulate a wide range of conditions for the site excitation energy transfer in the FMO pigment-protein.

Computer Resources in the Modelling
The high complexity of the HEOM method is translated to the computational resources required to solve the system. We used three commercial personal computers with an average of processing power in the range of 2.7∼3.4 GHz using parallel computing with at least four processing cores. They were tasked to process 150 different combinations of the previous parameters. For such a reason, the hierarchy depth was fixed to D = 3. Each process took between 6 and 8 h to be completed until the stage of equilibrium where the BChls populations became asymptotically unchanged. Out of this 150 combinations, only 121 yielded to adequate converging results based on an analysis of convergence using D = 4 as comparison with deviations constrained by 10 −3 . Depth higher values would have allowed higher precision as described in [51] but it would have required additional computational power competing with the extended analysis of parameters.

Ranges and Extent of the Parametric Analysis
While the constant r trap = 1 ps −1 is maintained fix in the analysis despite it is directly responsible from the trap efficiency for the BChls 3 and 4. Because we do not introduce the RC explicitly as a quantum system (as an instance in [47]), certain BChls (as the before mentioned) gather the final excitation upon the equilibrium.
The parameter λ k for the reorganization energies (5) reflecting the architecture of BChls inside FMO complex is responsible of the decoherence strength. Note there are differentiated λ k values for each BChl, despite the analysis in the literature commonly uses a unique value for all of them in order to describe the dependence on their main working value for such elements. In addition, the experimental knowledge about their concrete values is still poor and not completely conclusive. The most accepted value experimentally for λ k is approximately 35 cm −1 , then we will analyze the behavior in the interval [35,75] cm −1 to include other values discussed in the literature (particularly λ k = 65 cm −1 ) [53,54].
The parameter γ k , as it was previously discussed, is the bath cutoff frequency (the inverse of the bath coherence time-scale, the memory time of the environment). γ k has been experimentally determined around of 50 cm −1 . Here, we analyze it in the interval [30,50] cm −1 [50,54]. For that reason, in the analysis sections we refer to them as λ and γ simply. As it was previously stated, note all of those constants are normally expressed in energy units, but in spectroscopy they are traditionally reported in cm −1 through the converting factor 200πhc, while similarly, the switch between cm −1 and s −1 is reached through the factor 200πc.
Finally, the temperature T is an important factor in the equilibrium populations, the efficiency and in the coherence time of quantum coherence (a phenomenon expected to occur at very low temperatures avoiding that excitation not will be rapidly extinguished for the environment noise at high temperatures). Despite in the FMO the actual quantum phenomena occur at room temperature, experimental spectroscopy analysis are performed to relatively low temperatures around 77 • K.
Notably, green sulfur bacteria have been discovered living near from a black smoker at the coast of Mexico in the Pacific Ocean at its limit operation temperature of 60 • C at a depth of 2500 m. At such a depth, no sunlight penetrates, but it still lives from the dim glow of the thermal vent [16]. In the current analysis, we will consider the range [77, 347] • K.

Data Source Methodology
Data source for the analysis were obtained from the computer simulations of different combinations of parameters as it was previously depicted. HEOM method delivered time-series of matrices ρ from the initial condition until getting convenient stabilization values for the populations ρ ii in each combination. Those data are the basis or our analysis (and other complementary at the end of the article extending the main one) by tracking the behavior of several quantum features as entanglement, stability populations, coherence, together with decoherence and transfer times, as well as their relations. Such features are presented and developed in the following sections.

Parametric Dynamics of Populations during the BChls Excitation
The interest on the parametric analysis of some remarkable properties of the quantum regime lays on the precise knowledge for those parameters, whose combination sets notable differences in such quantum properties. Moreover, a set of well differentiated values for those parameters are reported in the literature for several species and strains, some of them exhibiting tiny differences in their compared behavior. Thus, parametric differences are becoming important for issues as: bio-inspired resources in quantum information, genetic manipulation of bacteria modifying such behavior, development of control on the FMO complex as technological resource, etc. Thus, a more deep understanding of the parametric dependence on the FMO complex quantum behavior is beginning a new era to manipulate biological resources for quantum applications. To the best of our knowledge, no previous analysis has been performed with such extent on those parameters despite interesting analysis has been performed as in [55] with the purpose of getting control on this nano-structure for quantum processing purposes.
The phononic interactions with the bath are crucial, they are resonance modes whose wavelengths are not easily observed directly from the experiment instead to its associated reorganization energy λ k in (5). An interesting issue about this parameter is the fact that while the absorption wavelength peaks λ absorption for the BChls in the FMO complex of green sulfur bacteria ranges approximately from 830 nm to 840 nm as they are observed in the spectroscopy absorption spectrum [56,57], instead, the characteristic excitonic wavelength λ α for the largest phononic frequency (α = 0) can be estimated as λ 0 ≈ 2πc ω 0 = 2πhc 2λ k = 142.9 µm (for λ k = 35 cm −1 changing the spectroscopy unities discussed previously), a very larger value compared with the characteristic absorption wavelengths but still meaningful to produce the decoherence, thus limiting the time of the process then re-initializing the process. Figure 7 shows some emblematic results obtained for the population evolution of each of the 7 BChls being present in the FMO monomer simply plotted against time, as it is customary in the most works [38,43,51]. With an excitation initially located only in the BChl 1, ρ(0) = |1 1|, Figure 7 shows only the first 4 ps despite simulations considered more extended time periods until the equilibrium. It corresponds to the period immediately after excitation energy is captured by coupled BChls 1 and 2. It is important to note that there is another possible path through BChl 6 initial excitation (see Figure 6), but here, we are specifically modeling an excitation coming through BChl 1. One can appreciate the excitation declines in pigments 1 and 2 while it is transferred towards BChls 3 and 4, in accordance with the energy flows presented in Figure 6. Lower temperatures yield longer sustained coherent oscillation periods although these are not directly related with the population oscillations seen close to t = 0. All plots shown in Figure 7 belong to the fixed parameter γ = 50 cm −1 contrasting the stabilization times between λ = 35 cm −1 and λ = 65 cm −1 at three distant temperatures (77 • K, 185 • K and 293 • K) which presented general variations around of 1/100 ps advisable in the length of the initial excitation time. In the following section we gather all simulations to get information about the parametric behaviour of those processes inside of each FMO monomer.
Some comparative outcomes should be commented at this point in spite the most approaches go through the dynamics simulation presenting similar plots as those of Figure 7, with variants in the modeling, initial conditions, complementary quantum systems (as a generic ground state or a sink representing the RC as in [47]) or still the master equation being used. Thus, outcomes included in the Figure 7 are consistent with those presented in [51,52] using HEOM method but exploring different initial conditions or considering insight parametric dependence with N = 8 BChls. Initial oscillations in the main interactive BChls 1 and 2 are notably present in those models. Equilibrium populations are comparable due they depend mainly from T. Outcomes generated using Redfield master equation could undergo to the failure of Tr(ρ) = 1 as it has been remarked in [36], while the HEOM method, under the correct convergence, exhibits a correct behavior inherited from is formulation [49]. As it was previously remarked, the use of Lindblad master equations comprising the main observable features in the FMO dynamics, successfully reproduces it in [47], not showing the extended initial oscillations as product the markovianity and closer to the Secular approximation compared in [36,58] with the Lindblad formulation combined with Redfield in the analysis of a single dimer interaction inside of the FMO complex.

Quantum Populations Reached in the Equilibrium Depending of T, λ and γ
The population evolution data obtained from the parametric analysis were simulated until the stabilization time was reached (the time when system populations reach asymptotic values transferring all energy to the final BChls). Then, asymptotic populations ρ ii ∞ ≡ lim t→∞ ρ ii (t) were fitted through the method of maximum likelihood to be represented by the surface contours shown in Figure 8, with complementary information reported in the Table 1. The average square distance between the simulation data and the outcoming approximation function falls in the range from the extreme values ∆ 2 ρ 22∞ = 2.2 × 10 −5 (best) to ∆ 2 ρ 33∞ = 1.6 × 10 −3 (worst). Unlike Figure 7, Figure 8 represents the final state of the system when this stabilization time was reached. Panels A, B and C, respectively correspond to populations ρ 11 , ρ 33 and ρ 44 , showing the populations with the greatest range of variability (other populations with less variability are included together in the fourth panel, D-G). Their corresponding insets show alternatively the continuous parametric dependence on the parameters for each final state in color in the common scale 0 (red) to 1 (blue). BChl 3 is clearly the pigment more populated by the time of system stabilization in all cases. Each plot in Figure 8 exhibits a series of surface contours for each approximated function with the ρ ii ∞ ≡lim t→∞ ρ ii (t) values given in color in the side scale on the left. The green, yellow and red mesh lines drawn on each contour surface in the Figure 8 corresponds to cuts of constant values of T (green), γ (yellow) and λ (red), respectively. These lines allow us to trace the effect of the each other two parameters on the equilibrium populations. The parametric analysis allows the detection of behavioural trends on the populations while parameters change. As instance, for BChl 1, the highest population ρ 11 ∞ is achieved at a parameter combination of the smallest values of T and γ, but highest values of λ (see the lower inset). The lowest values of ρ 11 ∞ are located for lowest temperatures in general. Notably, for ρ 33 ∞ , in contraposition with BChl 1, the highest populations occur for low T values, boosted slightly for the lowest λ and γ values, showing evidence about the FMO efficiency at lower temperatures (meaning that excitation is moving from BChls 1 and 2 to BChl 3 and 4 more efficiently). Other populations show a similar 2-fold behavior ( Figure 8D, E, F and G) versus B particularly from Figure 8) where higher equilibrium populations occur for the largest values of T and γ with the lowest values of λ. BChl 4 and 6 exhibit more exotic behavior as function of parameters noting there are other alternative input and output channels for the excitation under alternative initial conditions. The results observed in Figure 8 coincide with the energy paths depicted in Figure 6 and they are summarized in the second column of Table 1 reporting the ranges obtained in the analysis for each BChl considering all parametric combinations performed. Note the largest range of variability for BChl 3.

Coherence Depending of T, λ k and γ k
BChls in FMO complexes exhibit an initial increase and sustained quantum coherence, so we are interested to define quantitatively such an aspect. A quantum state is said coherent if it exhibits the form of a pure state: |ψ = ∑ i α i |φ i . Still, in terms of the density matrix, ρ = |ψ ψ| = ∑ i,j ρ ij |φ i φ j |, such state is said non-coherent if ρ ij = ρ ii δ ij with ∑ i ρ ii = 1 and ρ ii = 0 for at least a couple of values of i. It means that a state exhibit coherence if at least one of the ρ ij = 0, i = j. Then, a coherence measure should fulfill certain requirements due to the quantum nature of such a property.
If C(ρ) is a coherence measure, then it should fulfill [59]: (a) Non-negativity, C(ρ) ≥ 0 and C(ρ) = 0 ⇐⇒ ρ is non-coherent. (b) If Λ i is a completely positive trace preserving map (CPTP) lowering the coherence, then C(Λ i [ρ]) ≤ C[ρ]. c) C(∑ i p i ρ i ) ≤ ∑ i p i C(ρ i ), which means C is convex, the mixture of states has less coherence than the average of coherence. There are several candidates for such a measure. If the set of the non-coherent states is Γ, the first one is a measure based on the minimum distance from the any non-coherent state by using the l 1 -norm [60]: for σ diag = ∑ i ρ ii |φ i φ i | having the minimum. C l 1 (ρ) ≤ d − 1 for qudits is bounded by the maximally coherent state |ψ mc = 1 √ d ∑ i |φ i . Another measure, the relative entropy of coherence [59], defined by the relative S(ρ σ) and the von Newmann entropies S(ρ) [61]: measuring the minimum amount of noise to destroy the coherence remaining in the system. The minimum is reached again by σ = ρ diag . It is bounded 0 ≤ C re (ρ) ≤ C l 1 (ρ) ≤ d − 1 when S(ρ) is defined using ln(ρ) instead log 2 (ρ) [62]. Despite C re (ρ) fulfills the entire previous properties for a desired coherence measurement [59], in our analysis, we use C l 1 (ρ) for simplicity (despite last property is not fulfilled by it [59], but alternatively this quantity is nearly related with one of the entanglement measure used for such a purpose also, as we will be seen). For our parametric analysis, Figure 9A is a graphical representation of the parametric dependence of the system's maximum coherence. For the original data obtained from the simulations, the entire general range for C l 1 max (ρ) in the region became [1.22, 2.33]. As before, departing from the group of parametric simulations, we are calculated for each one the maximum value of C l 1 (ρ) reached and then approximated through the method of maximum likelihood (with ∆ 2 Then, some contours for constant values of C l 1 max (ρ) in agreement with the color scale below (note the maximum mathematically possible value is d − 1 = 6). Again, the highest coherence values are found at small values of T and λ (elasticity of the system under the bath interaction), but with larger values of γ (the reciprocal of average coherent time scale). Just as in Figure 8, green, yellow and red contour lines correspond to constant values of T, γ and λ, respectively. The black screens around report specific contour line cases for a fix characteristic value of each parameter as it has been analyzed in similar studies. One must note that while the gradient of colors on the contour surfaces are related to their maximum coherence value, the colors of the lines on the black screens make reference to the plane they belong to. In the case of the contour lines belonging to the constant value of λ = 65 cm −1 (projected on the bottom black screen) one can gather that there is a minimum found between the two brown lines laying on the yellow surface contour. In contrast to the monotone behaviour that prevails in both constant parameters T and λ, the case of constant parameter γ = 50 cm −1 shown in the back black screen has a quadratic like behaviour.
Other convenient definition involves to the coherence time. Quantum regime is reached fast departing from the initial excitation ρ(0) = |1 1|. Thus, we define T ch as the time when the maximum coherence as measured by C l 1 (ρ) (or C re (ρ)) is reached. Otherwise, we define the time of decoherence as an emblematic time during which the main excitation is sustained but it begins to decay considerably. As instance, here we will consider first the half-life time of decoherence T d 50% as the time when the initial excitation on the BChl 1 decays to one half. In the current work, we additionally use C l 1 (ρ) to quantify its correlation with the T d 50% behavior through the relation between both of the last times. Thus, the difference ∆T ch = T d 50% − T ch becomes particularly interesting in the quantum behavior of FMO because it reports the time for the coherence holding. We analyze such parametric behavior in the current section reported in the Figure 10. There, the maximum coherence C l 1 (ρ) for each simulation performed is plotted versus ∆T ch as function of T, λ and γ through the symbols on the right charts. Note that for lower temperatures the maximum of coherence increases as could be expected while the decoherence is really fast together with the decay of initial excitation near from the maximum. This behavior is almost independent from γ and λ despite higher values of γ increases the maximum of coherence while the opposite is true for λ. By rising the temperature ∆T ch increases but the behavior is almost preserved for other parameters. There, for larger T, the maximum of coherence increases with λ and decreases with γ, the inverse of the memory time of the environment. In other issue, T d 50% is shown as good predictor for the decoherence time, while ∆T ch represents the delay to maintain the coherence after the first peak in the maximum. We will analyze the form of the coherence C l 1 (ρ) as depending from the entanglement below.  In other issue, for the original data obtained from the simulations in the region, for T d 50% the entire general range became [0.05, 3.98] ps. Figure 9B represents similarly the corresponding analysis and representation for the parametric dependence of the half-life coherence time T d 50% (with ∆ 2 T d 50% = 4.4 × 10 −2 ) in agreement with the color scale below (in ps). Coherence is maintained in the system during an extremely short period boosted by the initial excitation, reaching its maximum almost immediately and then decreasing quickly. The behavior being found exhibits more exotic issues due to the multi-factorial dependence. The right black screen representing a constant parameter of T = 77 • K on the right shows extreme half-life times at the center of closed curves. The highest half-life coherence time is found at small values of T and γ, and it is few affected by the λ values with exception of folding found at the lowest λ values exhibiting the non-linear behavior in that region depending of the combination of low values of T and lambda.

Entanglement Generation Depending of T, λ k and γ k
Each FMO dimmer has components interacting among them, transmitting the initial excitation in the antennas (BChl 8 for the model with seven BChls, BChl 1 or 6 for the original model with seven BChls) to the remaining BChls through multipartite interactions. This dynamic generates complex entanglement among BChls [63]. Despite the limitation to quantify it, we can still address this phenomenon by using the partial trace criteria to get the concurrence as a measure of entanglement between different parts [61]. First, we trace the whole system except for two BChl states k < l, and then all subsystems except the k-th BChl to get, respectively: Then, taking the concurrence definitions for each kind of expression: the first for a two system density matrix, and the second for one system related with the remainder [58,61], we get: which can be respectively interpreted as the entanglement (a) among systems k and l [43], and (b) between the system k and the remaining BChls [58]. Each concurrence ranges from zero (separable) to one (maximally entangled). Figure 11 is a compact graphical comparative representation of the entanglement dynamics in the excitation process through the previous concurrences present between each BChl and the rest of the system, as well as each one with each other at the fixed γ = 50 cm −1 and (A) λ = 35 cm −1 or (B) λ = 65 cm −1 as a function of t and T. These values are obtained from the application of (36) with the density matrices (ρ kk ) generated by the HEOM calculations previously depicted. The graphs show the time evolution and temperature dependence of the entanglement through a color gradient in the scale 0 (red) to 1 (blue) indicated in the color bar above. Each circular sector of the plot represents the concurrence experienced by each BChl molecule in the FMO complex. The central circular region represents C {k} for each BChl indicated by numbers around. While, the following ring sectors represent C {kl} pairing the numbers around with those for each track on the right side. In each region, time t increases in counterclockwise direction from 0 ps to 2 ps repetitively inside each track (indicated by the blue arrow), while radially outwards each track has a temperature gradient from 77 • K to 347 • K (indicated by the brown arrow).
It is important to note that at the beginning (as well as final states in the equilibrium) of the transfer process, all BChls are separable states. It is only during a brief period of time that entanglement appears. For that reason, only the first 2 ps are presented in Figure 11. For C {k} in the center, we note that the main BChl entangled with the remainder BChls is BChl 1 followed by BChl 2, and then by BChl 3 and BChl 4 with a delay. Other BChls have a lower global entanglement. It is notorious that the highest entanglement readings are usually between BChl 1 and the rest of the system, but particularly with BChl 2 and 3 as seen in the corresponding rings of both Figure 11A,B. Of further interest is the increasing in entanglement between BChls 3 and 4, 4 and 5, and 4 and 7 as time increases. This is due to the fast excitation is being transferred to these pairs coupled pigments (see Figure 6). Note how the increasing of λ slows and reduces the entanglement process particularly for larger temperatures as seen in Figure 11B. BChl 4 is an exception of such behavior particularly notable for the largest temperatures. Results shown in Figure 11 can be also traced back to Figure 6, where we can relate the entanglement effects with the energy transfer paths.  Table 1 also reports the maximum concurrence that each BChl presents with the rest of the system C {k} at the lowest and highest temperature conditions during the evolution process (columns three to four for γ = 50 cm −1 , λ = 35 cm −1 and six to seven for γ = 50 cm −1 , λ = 65 cm −1 ). Columns five and eight of the table enumerates for each BChl with which of other BChl pigments presents the highest value of concurrence C {kl} (in parenthesis). Note the dominant entanglement with BChl 1 for λ = 35 cm −1 .
Another important characterization in the parametric behavior of the FMO is given by the relation between the coherence C l 1 (ρ) and each concurrence C {kl} : C l 1 (ρ) = ∑ k<l C {kl} . In such terms, coherence could be seen partially as sustained by the entanglement between pairs. Thus, we perform an analysis to advise how the coherence is sustained by certain pair through the excitation time in the Figure 12.
For several temperatures depicted in color in agreement with the color-bar besides, the evolution of C l 1 (ρ) (the first 2.5 ps) shows its decay transitioning from the entanglement between pairs (1, 2)-solid line-and (3, 4)-dashed line-as main providers of concurrence. Lower temperatures first slow such a transition except for higher ones, while a larger λ (see panels A and B) accelerates such a transition. As it was previously commented in the previous subsection, higher concurrences are present for lower temperatures. Thus, during the quantum regime, the evolution is dominated by two processes of entanglement in the system involving mainly to those pairs of BChls, despite clearly all of them share the entanglement process in a lower level in the multipartite interaction. Parameters as T and λ boost or slow each one of those main entanglement interactions, while the level of coherence becomes mainly dependent from T and any delay appears as function of those parameters as is observed in Figure 12, remarking the dominance for the process established on the left-side of scheme on Figure 6.

Efficiency Depending of T, λ k and γ k
The so called efficiency η of quantum energy transport inside of each monomer in the FMO is defined as the asymptotic sum of populations in BChls 3 and 4 due they are the responsible to transfer the excitation to the RC: Together, the time of transference is the time T t 95% that system uses to funnel the excitation to the exiting system at the 95% of its final equilibrium value (as instance on the BChl 3 or 4 connected to the RC). In the current work, we use the joint population of BChl 4 and BChl 3 to fix this time because they are the main populations remaining among the whole BChls while they are the responsible to transfer the excitation to the RC.
Following the previous parametric analysis performed through HEOM calculations, we analyze η in Figure 13A as function of parameters T, λ and γ (∆ 2 η = 2.1 × 10 −4 ). While for the original data obtained from the simulations for η in the analysis region, its general range became [0.51, 0.99] ps, the contour surfaces exhibit approximated constant η values by fitting in agreement with the color scale below. As before, on them, green, yellow and red contour lines correspond to constant values of T, γ and λ, respectively. The black screens representing once more characteristic constant parameter values showing the contour lines of the corresponding efficiency contours. In the case of the bottom black screen for λ = 65 cm −1 , we appreciate an inflection point in the middle of the graph where efficiency appears to increase and decrease in crossing directions in the γ − T plane for λ = 65 cm −1 . Efficiency values increases if T decreses in the γ = 50 cm −1 plane, but it is not the same behavior for lower values of γ where higher efficiencies are sustained at larger temperatures. In a related analysis, for the source data obtained from the simulations, the range of T t 95% in the region was [3.11, 20.82] ps, denoting a wide range for the parametric behavior. For the approximation by fitting the data (∆ 2 Figure 13B provides complementary information about the process comparing the transference times T t 95% through the parametric values. Once more, high values for λ appear to have little effect in transfer time, comparable with the behaviour observed for this same parameter in half-life coherence ( Figure 9). Instead, T t 95% values increases with T and γ. On the other hand, transfer times show a better performance (lower T t 95% ) at lower values of T while are combined with small values for γ. Parametric characterization of FMO complexes in the quantum regime becomes important to understand the multivariable functioning of those structures no depends monotonically on T. In fact, FMO complexes are in the chemical mesoscale, thus exhibiting lots of variants in their structure. Chlorophylls and bacteriochlorophylls, pigments of photosynthesis are not exclusive of green or purple sulfur bacteria, instead variants are present in many other live structures as a result of genetics, adaptation, and evolution to their environment, so the comprehension about their operation and characterization is still important to gain knowledge and control for further based bio-applications. Still, those complexes exhibit differences both in their main and scaffolding structure. They are ruled by genetic laws and evolution, not only as a rigid and unique chemical object in nature. Such variants could be modeled and characterized to pretend being artificially designed and produced for concrete applications. Still, quantum behavior will be responsible or not of the main photosynthetic light-harvesting in FMO complexes, their quantum operation properties could provide a precise identification of their structure not affordable for observational techniques. In the final section, we tackle such issues related with the parametric characterization based on alternative studies on another FMO complex being present in near related bacteria in order to comprise the level of differentiation.

Characterization of Genetic Differences in the FMO Complexes with other Species and Strains
This section concerns with variations advisable in the quantum characterization in the FMO complex as compared with other species or strains reported in the literature. Despite, those studies are diverse and still introductory in many cases, so only partial analysis are available. Still, Prosthecochloris aestuarii has been widely studied as Chlorobaculum tepidum, thus there are similar data regarding it. A brief parametric analysis of the quantum regime for the two different species (Chlorobaculum tepidum and Prosthecochloris aestuarii) whose Hamiltonians have been obtained from [1,13,37] (see Appendix A), will be now discussed. The effects of previously mentioned parameters on the quantum behavior will be analyzed according to the results obtained in the modeling. Then, a genetic differentiation between GSB species and its relation with efficiency is discussed. In additional terms, we analyze other alternative initial conditions and pathways in the simulation to compare their differences. The section concludes with the introduction of the concept of genetic engineering and its possible effects on energy transfer efficiency of the complex.

Alternative and Extended Studies on FMO for other Species and Strains
It is important to consider that GSB, as all living organisms do, possess regulatory mechanisms to better adapt to the changing circumstances of their environment. For instance, the bacteriochlorophyll synthesis is strongly regulated by light intensity, meaning that pigments and chlorosomes can be multiplied under light-limiting conditions or their biosynthesis can be halted under rich light conditions, effectively acting in accordance to the situation [15]. Regulation of light harvesting is thus strongly related to the effect the environment has on the genetic expression of the organism. Antenna proteins can be regulated on multiple levels from mRNA (messenger RNA) transcription or even be submitted to protein degradation [64]. Although the environmental effects on gene expression are of interest for the understanding of photosynthesis in vivo, these go beyond the reach of this study. In this work, we focus on the parametric analysis of some functional differences within FMO complexes resulting from the genetic variations between different species and mutations. For this, we refer to the experimental outcomes in literature [19,20,37,65].
In order to set a comparison between the quantum regime behavior of Chlorobaculum tepidum and Prosthecochloris aestuarii, a set of two different Hamiltonians of Chlorobaculum tepidum (Equations (A1) and (A2)) and Prosthecochloris aestuarii (Equations (A3) and (A4)) have been considered. In the Figures 14 and 15, two parallel sets of simulations has been performed for Chlorobaculum tepidum and Prosthecochloris aestuarii, respectively. We have performed simulations using their respective Hamiltonians as well for N = 7 as N = 8 cases considering γ = 50 cm −1 for λ = 35 cm −1 and λ = 65 cm −1 at the room temperature of T = 293 • K. Those values in the present time have been only barely estimated in the literature. They are only measured for some BChls, and for one or other specie, being settled in the same range. Thus, we use the most viable values to have an insight about the comparative behavior of their quantum regime considering synthetically some quantum features previously discussed.

FRET
Due to the analysis for N = 7 should to consider the initial pathway (ρ(0) = |1 1| as in our main analysis, or ρ(0) = |6 6|), we develop those two single main initial site-like conditions to see the differences for such parametric values. For N = 8, instead we will consider ρ(0) = |8 8| as initial condition, we use the so called Förster resonance energy transfer (FRET) initial condition [27]. This condition better reflects the initial excitation coming from a light source with a wide spectrum of frequencies as it happens the natural habitat of FMO complexes. It assumes that the excitation is transferred from the baseplate to the FMO first populating the FMO excitonic states | k , then exciting mainly the closest BChl i in the antenna: In our case, i = 8 will be considered in order to analyze how this condition reflects the combination of two pathways being analyzed by separate for N = 7 (note in the studies with N = 7, the cases for i = 1 or i = 6 commonly are considered separately for this kind of initial conditions [51], despite here we are considered the initial site-like conditions). The construction of such initial states is direct by calculating the eigenstates of H S in each case. Note that other initial conditions are advisable for incoherent light sources by considering their absorption distribution spectrum I(ω) as in [27].
In each one of Figures 14 and 15, left column corresponds to λ = 35 cm −1 and right for λ = 65 cm −1 (γ = 50 cm −1 , T = 293 • K). First and second rows presents the case N = 7 BChls with ρ(0) = |i i|, i = 1, 6, respectively, using the located site initial condition, while third row correspond to ρ(0) = ρ 8 FRET initial condition. All simulations were extended further away than t = 10 ps (see the insets in each plot) almost into the thermal equilibrium, but the main plots are focused on [0, 6] ps to appreciate the initial details in the quantum regime.
In a first glance, some aspects are outstanding. In terms of the two main pathways in N = 7, a very similar behavior is observed for both inclusively for the final populations of ρ 33 and ρ 44 . The case of λ = 65 cm −1 induces a more slow relaxation of the excitation phenomenon. The same aspect is true for both species despite equilibrium populations for P. aestuarii appears larger for λ = 65 cm −1 . For the case N = 8 the same features apply. For the ρ 8 FRET initial condition, it is outstanding that despite the lower population occupancy for ρ 11 , ρ 22 , the equilibrium populations for ρ 33 and ρ 44 remains in a comparable level for the other cases analyzed with N = 7 involving only one of pathways (but depicting a little decreasing for ρ 33 ). For both bacteria, the similitude is notable, despite the initial populations imposed by their respective ρ 8 FRET are different). Note how the ρ 8 FRET initial condition boosts predominantly the pathway involving the BChls 1 and 2, which were selected for our main analysis in the previous sections. In the following section we develop a more deep comparative analysis based on the quantum features studied previously for C. tepidum based on the current simulations in Figures 14 and 15.

Comparative Analysis between Chlorobaculum tepidum and Prosthecochloris aestuarii
We address in this section the analysis of main quantum features predicted for each bacteria during the excitation. Such features are not completely appreciated in the evolution of Figures 14  and 15. We are selected the coherence versus the efficiency to compare their evolution in the process. To compare the relation between the coherence and the gain in efficiency, both as a derived effect of quantum behavior, we are performed the Figure 16 whose information will be complemented with the Table 2, and making reference to Figures 14 and 15. Each plot Figure 16A-C contains the joint evolution between C l 1 (ρ) and ρ 33 + ρ 44 until their convergence to η for ρ(0) = |1 1|, ρ(0) = |6 6|, and ρ(0) = ρ 8 FRET , respectively. Plot for C. tepidum (blue) and P. aestuarii (red) are presented for each λ = 35 cm −1 (clear) and λ = 65 cm −1 (dark) for the room temperature T = 293 • K and γ = 50 cm −1 .
In Table 2, we include the corresponding values of C l 1 (ρ) max , T d 50% , T t 95% and η. As could be expected, the behavior in each plot is barely similar. Note the slowing previously adverted when compare λ = 65 cm −1 with λ = 35 cm −1 (dark) despite the concordance between efficiencies is notable in all cases because in their relatives markovian models final populations depend only on the energy spectrum and temperature. The effect on η is more advisable for the group of λ that the group of bacteria in the most cases. Despite, note the sensible differences in T d 50% for the case ρ 8 FRET despite without impact on η. In contrast, ρ(0) = |6 6| exhibits notable low values for the transference. While, the bacteria type exhibits main differences between them for the case ρ 8 FRET . Despite, for the general differences expected, in a coarse grain terrain, behavior is surprisingly consistent as between bacteria as between pathways as it has been remarked by [20] in terms of the robustness of the light-harvesting process, despite intermediate process become clearly different. It is remarkable the different behavior for the case ρ 8 FRET with a really sudden increasing and decreasing of coherence without a sensible loss in the efficiency. Thus, data reported in Table 2 show for C l 1 (ρ) max main differences for ρ 8 FRET by combining both pathways. Decoherence times T d 50% and T t 95% are notably consistent only depending from λ. Note the slowing for ρ 8 FRET and λ = 65 cm −1 . An outstanding discrepancy is observed for T d 50% for C. tepidum in ρ(0) = |1 1|. Efficiencies reported are very similar at the end inclusively following certain corresponding trends when two λ values are compared. Despite the comparison with respect to λ, which could be uncertainly different for each bacteria, there are sensible and non-trivial differences in the decoherence and transference times between them which constitutes a quantum fingerprint in their structure despite such differences are complex to follow in their biochemical and genetic architecture.

Genetic and Structural Differences between Species and Strains and Their Manipulation
The FMO complex is characteristically present in GSB and it has actually been key for further identification of species and strains [16,18,66]. The gene codifying the FMO complex (fmoA gene) presents small variations within the distinct species and strains of green sulfur bacteria, resulting in a growing phylogenetic tree [16,18]. Regardless of this fingerprint characteristic, these genetic variations have consequences in the final protein structure of the FMO complex. In the cases of C. tepidum and P. aestuarii, a comparison between their crystal structures (with PDB IDs 3BSD and 3EOJ, respectively) using a basic local alignment search tool (BLAST) reveals a 78% of identical amino acids and a total of 87% of similarities, due to the fact that some amino acids that differ between both structures have similar behaviours and characteristics as their substitutes. These structural differences for the fmoA gene can be observed in Figure 17, where the FMO complex of both organisms are being overlapped in Figure 17A, and their differences are highlighted in blue in accordance to their respective sequences shown in Figure 17C. The green arrows and pink serpentine lines are a linear representation of the protein β-sheet and α-helix structures, respectively. Figure 17B, shows an overlap of the 24 BChls (8 BChls per monomer) of both organisms exhibiting tiny differences derived from punctual changes in the orientation by the insertion of alternative amino acids, then reflected in more extended parts of the structures. Slight orientation shifts are observed between the molecules as an effect of the different protein structures, which surely modifies slightly their dipolar momenta as well as their relative orientations, thus changing the respective Hamiltonians H S . Such differences are then responsible from the deviations in the quantum fingerprint observed in Figure 15.
Living organisms can be exploited for chromophore assembly given they possess the genetic blueprints and molecular machinery for their biosynthesis [67]. Genetic engineering and biotechnology provide a new approach towards designing efficient energy transport by manipulating these organisms but their correlation with quantum properties could set a bridge between the design and the functioning effects. For such a goal, the selection and structural arrangement of chromophores, and inter-chromophoric couplings play a crucial role. Some challenges for the customization of energy transfer mechanisms are the structural control of the individual light-harvesting pigments and the achievement of beyond-Förster energy transport [67]. The manipulation of such structures in living organisms is possible by the use of programmable genes that can modify the positioning of the binding sites, thus enabling the creation of defined chromophore networks and ultimately controlling the energy transfer [67] or other notable novel applications in quantum processing taking advantages from biological-based structures as quantum walks [30,32] and incoherent quantum evolution [10]. The possibility to track the discrete differences in the molecular structure of FMO versions into a mapping on the continuous quantum features as those in Figure 15 is a future challenge. showing residues belonging to β-sheets (green arrows) and α-helices (pink curves). Punctual differences in the amino acid sequence (blue) in (C) match the sites illustrated in (A). Protein Data Bank files 3EOJ (P. aestuarii) and 3BSD (C. tepidum). Generated using iCn3D.

Conclusions
Either in the coherent or incoherent regime, the study of the energy dynamics within the FMO complex yields an important insight into other photosynthetic proteins and into the process of photosynthesis as a whole. For quantum theory, FMO complex has been the simplest biochemical nano-structure to be studied, despite sufficiently complex in physics to be faithfully described. It has let the possibility to observe and to model its quantum behavior. Quantum mechanical approaches, either experimental or theoretical, have let us understand the properties shaping its mesoscopic and macroscopic properties [13]. Its more notable feature, its quantum functioning at room temperature, has attracted the interest in an extended part of the scientific community [23].
In the current work, the use of HEOM method has let to approach its non-markovian behavior to reproduce closely its dynamic [46]. Trough of computer simulations, scientists try to reproduce the quantum elements involved there by exploiting the most extreme master equations comprising the interaction among extended systems. This complex system has boosted the development of quantum concepts far away from those simplest just arising in microscopic systems. Those computer simulations representing the most inclusive models and challenge laws are the departing point to understand how to afford complex systems from the more basic quantum laws. We are contributed to map the quantum regime of FMO complex dynamics in terms of unique features related during the coherent regime through the one of the main pathways in the complex unveiling interesting non-monotonic behavior on them.
Thus, under the functioning of such nano-structures until now described through approximated chemical laws, we are discovered how quantum mechanics operates conducting lots of phenomena shaping macroscopic notable features. In the FMO, the presence of quantum entanglement is reflected as long-lived coherence in those light-harvesting complexes representing an important connection between the quantum realm and higher systems with complex but despite natural environmental conditions for the human being [43,52].
Quantum mechanics has provided the necessary basis for the continuous development of approximation models explaining those dynamics. Through this approach, our parametric analysis of the FMO complex allows the identification of key individual parameters and their effects on the dynamics, quantum features and efficiency of energy transfer in the coherent regime, opening the road to understand and to try its modeling. This information results in the determination of optimal operation conditions tentatively controllable.
In addition to this analysis, the comparison between FMO complexes of different species (and their strains in the near future), with respect to those parameters, contributes to a quantum characterization of the different versions of the complex and the species and strains as quantum systems. It lets in the futuristic science a better comprehension about its possible genetic manipulation or applicability as technological resources [19,64]. It establishes a relationship between the observed quantum behavior and the genetics dictating the structure on which those phenomena are observed [20]. In other hand, it provides us the possibility to try the quantum behavior as fingerprints in their complex structure and characterization (as well as it becomes true in the past for the atoms through their emission or absorption spectrum). More extensive genetic comparisons among species and strains may provide information on optimal structural organization of the complex [12].
Although great efforts have been made to describe the energy dynamics within the FMO, there is still much to be unveiled about the actual functioning of this complex in vivo exhibiting non-trivial and non-linear dependence of some of their operation parameters. There is still not enough information about the initial and final conditions of the system since this LHC acts as a quantum bridge between the chlorosome baseplate and the reaction center [1]. The exact way in which excitation is transmitted from the baseplate to the FMO and from the FMO to the reaction center is still unclear. Solvent interaction between these junctions may also affect excitation energy transfer between them.
Much progress should be made in the following years by means of more faithful quantum models and disruptive experimental techniques by considering its inner interactions, but so far this description of the FMO complex can already inspire the possibilities to be used for quantum processing, quantum cryptography, together with controllable more efficient solar energy harvesting [9]. Comprehension of energy dynamics within light-harvesting complexes is key not only to comprehend one of the main factors for the viability of life, but to harness its power towards multiple novel applications based on the quantum realm [10,55,67]. Joint dipole-dipole interactions among BChls inside each monomer in the FMO complex considering seven or eight BChls have been obtained through spectroscopy studies complemented with theoretical analysis. We report those matrices in units of cm −1 used in multiple analysis for the behavior of each monomer of FMO complexes. For N = 7 BChls in C. tepidum (from Tables 1 and 4 in [13]) with a constant diagonal offset of 12,210 cm −1 to set the lowest site energy to zero, H S is: Despite our main analysis is performed for C. tepidum, other specie under an intensive research is P. aestuarii. Thus, for N = 7 BChls in P. aestuarii (from Tables 2 and 4 in [13]) with a constant diagonal offset of 12,230 cm −1 to set the lowest site energy to zero, H S becomes: While for the H S in P. aestuarii with N = 8 BChls [1], we have developed a similar but more limited parametric analysis of its quantum properties [52]. This N = 8 Hamiltonian has been obtained using the charge density coupling method [13,68] assuming a standard protonation pattern. Values presented consider a relative dielectric constant border between baseplate and RC complex border = 80, representing the experimental conditions of a fully solubilized FMO complex. In our analysis, we are also applied the HEOM method on P. aestuarii to get a brief comparison of the quantum regime features as for C. tepidum in our main analysis.

Appendix B. Linblad Equation Derivation in Brief
Lindblad master equation is obtained following different assumptions in order to remove the details of the bath system but leaving sufficient information to predict the effect of it on the main system. The interaction between the system and the bath is assumed bilinear in the form ofH S−B = ∑ N 2 −1 i=1 γ iSi ⊗B i . In the Lindblad model, the integrals in (17) are performed by assuming that the operators of the bath have correlations satisfying Tr(B i (t)ρ B B j (t )) = δ(t − t )Tr(B i (t)ρ B B j (t)), which conducts to the non-diagonal form of Lindblad equation (h ij = 2γ i γ j Tr(B iρBBj )):ρ which by diagonalizing h ij = ∑ α,β T iα J α δ αβ T † βj = ∑ α T iα J α T † αj and returning to the Schrödinger picture can be expressed as:ρ where L α = J α ∑ iSi T iα , represents the so called Lindblad operators. Nonetheless, the relation of L α with the physical arrangement is not always direct because it is based in generic operators S j or L α stating a basis.

Appendix C. Redfield Equation Derivation in Brief
In the Redfield approach, physical information of the system-bath interaction is retained. In fact, by defining C ij = Tr(B iBjρB ) and C ji = Tr(B jBiρB ) = C * ij , and dropping the time dependence by assuming the stability of the bath in (17). If t 0 → −∞ and then performing t → τ = t − t , (13) becomes:ρ which is reduced to its classical form (after to return to the Scrödinger picture):ρ which physically describes the coupling with the environment.

Appendix D. Translating Master Equations into Matrix Differential Linear Equations
We will develop the translation from the HEOM master Equation (28) into the superoperator version through the application of the rule (29). For other quantum master equations, we simply define the supervector ρ = (ρ 11 , ρ 12 , ..., ρ 1N , ρ 21 , ..., ρ NN ) (note other versions can be advised varying the order of the elements of ρ) and we go directly on the application of (29). In the HEOM method, due to the coupling among ρ and the auxiliary matrices ρ n through ρ n k± , we proceed similarly by moving each ρ n into its corresponding ρ n . Then, the HEOM is a system of linear differential equations for the supervector set Ω ≡ { ρ n } labeled with the vectors n of deep D as was defined in Section 3.3.3, each one of size N 2 .
Then, the application of the rule (29) transforms each term in (28) into the matrix product between a matrix of size N 2 × N 2 (superoperators) and each one of the supervectors in Ω: , for s = 1, ..., N, i = 1, 2, 3 are constant matrices. The numerical easiest process to integrate the equation is to consider simply that ρ n (t + δt) ≈ ρ n (t) +ρ n (t)δt (better integration methods could be implemented to avoid tiny values of δt). The main problem with this approach to HEOM is the polynomial growing number of equations, for D = 3 and arbitrary N, there are 1 + N + N(N+1) 2 + N(N+1)(N+2) 6 equations. They are 120 for N = 7 and 165 for N = 8 involving matrices of size 49 × 49 and 64 × 64, respectively, on a fine grained time period of tens of picoseconds.