{'fraction_non_alphanumeric': 0.026130961508971016, 'fraction_numerical': 0.0068394418034043855, 'mean_word_length': 4.7253731343283585, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A compact theoretical model for opto-electronic devices based on quantum dot arrays

We present a theoretical model describing the photo-electrical response of a system composed of quantum dots embedded in a dielectric matrix. The model is based on the non-coherent rate equations and the Transfer Hamiltonian formalism. The self charge was included. Within this methodology, the response of the system only depends on fundamental material parameters and its geometry. Transport through several quantum dot configurations was simulated obtaining current-voltage curves in dark and illuminating conditions for three different scenarios: single one quantum dot and two quantum dots in parallel and serial configurations. Despite the simplicity of the model, it has been used to reproduce successfully experimental results.


Introduction
Silicon quantum dots (Si Qd) embedded in insulator matrices have opened a new branch of possibilities in the electronics based on quantum confinement and Coulomb blockade effects [1,2]. Compared to the standard bulk Si technology, devices based in these structures have increased the fabrication and conceptual complexity. Qd can be used in single-electron transistors [3], new memory concepts [4] and photon or electroluminescent devices [5]. Concerning opto electronic devices, the possibility to modify the energy band gap as a function of the Qd radius [6,7] can be exploited in order to build light absorbers for photovoltaic applications [8][9][10]. This property can be used to solve the fundamental problem existing with silicon photovoltaic allowing to increase the photo-response of the device in a wide range of the solar spectrum. The spatial confinement in Qds also helps to overcome the Shockley-Queisser limit [11,12], i.e. the maximum theoretical efficiency of opto-electronic devices. Besides, discrete energy states appear inside the wide band gap of the insulator matrix, making possible tunable band gap devices [13]. Nowadays, to increase the conversion efficiency of silicon photovoltaics tandem cells are used [14,15], i.e., a stack of absorber layers with different bandgaps cover a wider range of the solar spectrum than a single bandgap absorber layer. These structures conform the third generation of solar cells that can increase the energy conversion efficiency limit up to 60% [16]. Therefore, a whole new generation devices are based on arrays of Qds that form the stack layers, thus, is needed a detailed explanation of the opto electronic response of these arrays. The movement of the electrons under illumination conditions in the Qds are governed by three basic mechanisms: i) promotion to an excited state assisted by a photon (photon absorption), ii) relaxation transitions from excited states to ground state (stimulated photon emission) and iii) tunneling processes to the neighbor Qds or electrodes, since the phonon scattering is suppressed by the phonon bottleneck effect [17]. Therefore, in order to generate a net photo-current, the photo-generated electron has to be extracted before recombining, thus, the charge transport from the location of photogeneration to the leads becomes a crucial point that governs the response of the system. Indeed, to obtain a good optoelectronic device not only the photo-generation has to be effective, also the charge transport extraction mechanisms must be efficient to enable electronic current flow. These transport mechanisms depend on tunneling processes thus being strongly dependent on the insulator material and also on the geometrical arrangement of the Qds. A compact model to explain the elastic electrical transport in the self-consistent field regime [18] in these kind of systems was reported by the authors previously [19][20][21]. It is based in the Transfer Hamiltonian formalism [22] in order to describe the currents between the different parts of the system. Assuming that the Qds are independent, the non equilibrium distribution function can be obtained solving a rate equation [23]. Moreover, the effect of the trapped charge in the Qd has been taken into account solving the transport equations and the Poisson equation self consistently, within the self-consistent field regime, explaining a first order the Coulomb blockade effects [18]. The aim of this work is to extend the previous developed transport formalism [19] and incorporate the optical processes in order to create a complete model to describe and simulate opto-electronic devices from a small number of basic material parameters. The rate equations have been modified in order to reflect the stimulated generation/recombination carrier processes and the transition rates have been derived from the classical electron light interaction taking into account all the possible transitions due to the optical processes. We present here a complete and compact methodology that can derive the electronic and optical properties and the response of a device based on Qds embedded in an insulator matrix placed between two electrodes from basic material parameters. Some cases based on Si Qd embedded in SiO 2 matrix have been studied in detail, the single Qd in two different scenarios, symmetric and asymmetric configuration. Different rates between the optical and current terms have been studied obtaining the key features in order to obtain a net photo-current without external bias voltage applied. Once this case has been studied, two more cases were analyzed, the parallel and serial case. Both cases were studied under external illumination and applied voltage. Finally, in order to validate the presented methodology it was compared to experimental results extracted from Prins et al. [24] obtaining successfully agreement between the simulations and the experiments.

Theoretical model
The general framework has been presented in previous works [19,20]. It describes transport between N Qds embedded in a dielectric matrix coupled to two electrodes (L lead and R lead) as a network of multi-tunnel-junctions. It is based on the non coherent rate equations [23,25] and the Transfer Hamiltonian approach [22]. In that framework, each Qd is treated independently and the non-equilibrium distribution function of each Qd can be obtained. This model provides intuitive results and it has successfully been validated to other methodologies like Non Equilibrium Green Function Formalism (NEGFF) with positive results [20]. From the Transfer Hamiltonian approach, the currents between the different parts of the system are described by the density of states of each part, the different occupation probabilities and the transmission probability trough the insulator barrier. The charge evolution inside each dot can be written and solved in the steady state extracting the non equilibrium distribution function of each Qd. Moreover, the influence of the applied bias voltage has been introduced taking into account the capacitive coupling between the different elements of the system. The effect of the accumulated charge in each dot is also taken into account within the self-consistent field regime. In order to describe the system as realistic as possible all the variables that appear in the model have been expressed as a function of basic material constant. For the transmission coefficients, a WKB expressions have been used in order to describe the tunnel probability under low and high electric fields. The capacity values have been described using expressions between a infinite plane and a sphere and between two spheres to calculate the Qd-lead and Qd-Qd capacities respectively. Therefore, the complete system can be fully determined by the material constants and some geometrical parameters. A complete explanation of the used expressions are presented in Ref. [21]. The previous model is used to evaluate the electrical response of the system under external bias voltage, our model considers that the electronic structure inside the Qd is not affected by external influences (light interaction or external electrical polarization), therefore, in the previous equation the equilibrium properties of the Qd are used: effective masses, confinement potentials and energy levels. We use it as a base to develop the photo-response of an array of Qds but some modifications are needed. We introduce the stimulated generation/recombination assuming that the optical and transport processes are independent as proposed by [26,27]. Therefore, the rate equations for each energy levels are modified as Neighboring Qds contribution where the subscript i refers to the i th Qd, j and k refer to the j th and k th energy level of the corresponding Qd. n i j is the non-equilibrium distribution function of the j th level in the i th Qd. q is the electronic charge andh is the reduced Planck constant. ρ L and ρ R are the density of states (DOS) of the leads evaluated at the energy of the energy level j, and ρ i j is the degeneracy of the j th energy level of the i th Qd. f Lj and f Rj are the distribution functions of the leads that we consider in equilibrium. Therefore, they are described by the equilibrium distribution function, Fermi Dirac, with the corresponding electrochemical potential, µ L − µ R = qV , where V is the applied bias voltage. |T Lj | and |T Rj | are the transmission coefficients between the Qd and the leads. If the j th energy level belongs to the conduction band we will use the barrier for electrons otherwise we will use the barrier for holes. Concerning the neighboring contribution, three kinds of processes are taken into account; 1) the conduction-conduction band transmission, 2) the conduction-valence band transmission and 3) the valence-valence band transmission depending on the energy levels involved in the tunnel process [21]. The neighboring Qd contributions are summed over all the Qds, except for that under consideration. The last terms correspond to the generation and recombination terms due to the light interaction [26].
These terms involve all the energy levels of the i th Qd as a function of the optical transition values R jk . The first light term is the generation rate while the second one is the recombination. The coefficient R jk is the transition probability between the levels j and k, it is explained in the next section. As expected, the generation/recombination terms depend on the occupancy of the energy levels involved in the optical transitions [28][29][30]. Solving Eq. 1 in the steady state condition, the non-equilibrium distribution function for each energy level of each Qd can be obtained. Therefore, the charge stored in the i th Qd is N i = j 2ρ i j n i j , where we take into account the spin factor 2. Once the charge has been calculated, the Poisson equation is solved and the self-consistent solution of the local potential and the charge is imposed [31].

Optical properties
As we have seen in Eq. 1 the transport properties depend on the DOS of the Qds which it describes the electrical and optical properties of the Qds. Therefore, a correct description of the electronic structure of the embedded Qd is needed. Following our strategy to keep the model as general as possible, we propose a simplified model to describe the Qds as a finite spherical potential well [32]. The confinement of the carriers is described by the potentials V e (r) and V h (r) for electron and holes, respectively. The height of the potential well is the difference between the conduction/valence band energy level of the dielectric matrix and the ones that form the Qd. The width of the well is 2R, where R is the Qd radius. We use the effective mass approximation (EMA) [33] to describe the electrons by replacing the free electron mass m 0 by an affective one, denoted m * , which takes into account the periodicity and the interaction with the atomic structure. Moreover, the carriers effective masses m * are assumed to have different values in the Qd and in the dielectric matrix. In this case, the coupling of different bands and Coulomb repulsion inside the Qd are neglected. Therefore, this model can be described by simple material parameters of the Qd and the insulator matrix. From now on, we write all the calculations for the conduction band, the valence band calculations are straightforward. The single particle radial Schrodinger equation is written as Here, m * (r) is the position-dependent electron's effective mass, l is the angular momentum quantum number, V (r) is the finite confining potential, E l is the electron energy eigenvalues for the quantum number l and R l (r) is the radial wave function. The angular dependence of the wave function, Y (θ, φ), is given by the spherical harmonics. The general solution of the wave function is where j l (k 0 r) is the spherical Bessel function (l th order), [34]. V 0 is the band offset between the Qd and the surrounding material, and m * Qd and m * Ox are the effective masses in the Qd and in the oxide, respectively. A l and B l are the normalization constants. The boundary conditions of the wave function and its first derivate require From these conditions, the eigenvalues for each l binding states are obtained. We remark that the degeneracy of each binding state is (2l + 1), which implies that the number of binding states increases dramatically when the Qd radius increases. Up to here, the internal structure of the Qd is described.
In the presence of a photon field, the effective electron-photon interaction can be introduced via the minimal coupling in the Hamiltonian [34]. Assuming that the photon can be treated as a perturbation, the interaction can be calculated starting from the Fermi's Golden Rule in the dipole approximation. This is the usual treatment of the light interaction and a formal derivation can be found in many text books [29,30,34]. For simplicity, we treat only the optical transition rates between different binding states for z-polarized incident light. The transition rate between the l and l' state can be written as where α is the fine structure constant, r the Qd dielectric permitivity, V Qd is the Qd volume, I(hν) is the flux of the incident light and ν is the photon frequency. z l l is the dipole matrix element, it reads as where we have used z = rcosθ. The selection rules for such transitions are determined by the angular part of the integral. It can be readily shown that the allowed transitions must fulfill the condition |l − l | = 1 and m = m. In addition, the oscillator strength which is related to the dipole transition can be expressed as P l l = 2m refers to the energy difference of the two levels. From Eq. 5, the absorption coefficient is given by In order to obtain the total absorption coefficient, a summation is performed over all possible transitions taking into account the occupation of the states. To obtain a smooth absorption spectrum, we replace the δ function in Eq. 7 with a Lorentzian function with a half-width Γ. Γ is related to the thermal broadening.

Results and discussion
We have done a systematic study of the electrical and optical response of a system composed by two electrodes, left and right, and a central region in which the Qds have been placed. The scheme of the system is shown in Fig. 1. The Qds were embedded in an insulator matrix. In our model, the system is described by the geometrical properties of the Qds, radius and relative position respect to the leads or other Qds and the material parameters that describe the insulator and the Qds. In table I, we show the list of parameters used to describe Si Qds embedded and in SiO 2 insulator [21]. The effective masses inside the Si Qd (m * Qd,CB and m * Qd,V B ), in the oxide (m * ECB , m * EV B and m * HV B ), the confinement potentials φ 1,ECB and φ 1,HV B for electrons and holes respectively, the Si bulk gap E gSi = E shif t,CB + E shif t,V B , and the SiO 2 dielectric permittivity ε r . We use T=300K. The study of the energy states of the Qd has been done assuming the finite spherical potential well using the Si and SiO 2 effective masses and the confinement potentials assuming bulk values. The binding states appear from the numeric solution of Eq. 4. Since the eigenstates and the wave functions have been obtained; the optical properties of the Qd, the absorption spectra, can be calculated. The Table 1: Parameters used in the simulation in order to describe Si Qds embedded in SiO 2 insulator matrix. m * ECB , m * EV B and m * HV B refer to the effective masses used to describe the different tunneling processes: conduction-conduction band transmission, conduction-valence band transmission and valence-valence band transmission, respectively. E shif t,CB and E shif t,V B are the energy shift of the binding states due to the Si bulk band gap. absorption spectra gives important information about the electron transitions between different energy levels produced by the interaction with an external photon and its probability. We present, in Fig. 2, the evolution of the simulated spectra as a function of the Qd radius using Eq. 7 for a z-polarized incident light. When the Qd becomes larger, the number of energy states in the Qd increases obtaining more possible transitions. Moreover, the value of E g decreases [7] thus, a redshift of the spectra is observed. For high photon energies the absorption decrease since we only consider the transitions between the discrete energy states of the Qd. Therefore, we neglect the continuum conduction and valence band states and all the possible transitions from/to these bands. Once the absorption spectra is calculated, the photo-current response of the different systems can be obtained. The systems are composed by two electrodes, right and left leads and a central region in which the Qds are placed. For simplicity, we assume an incident monochromatic light with irradiance I(hν) = 1W m −2 . Using the rate equation formalism described previously, the current can be obtained under external bias voltage conditions, illuminated or in both cases.

One single Qd
The first system under study is a single Si Qd under illumination connected to two electrodes with a constant external bias voltage applied. From Eq. 1 many scenarios appear as a function of the value  In the inset a scheme of the system is presented, a Qd of R = 1.06nm is placed between the two electrodes. a) Symmetric system d = d = 1.78nm. For V = 0 case the current is zero while for V = 0 the symmetry of the system is broken and net current appears. b) Asymmetric system d = 1.47nm and d = 2.09nm. A net current is obtained even at V = 0. The current peaks refers the position of the maximum optical transition probabilities.
of the optical and transport terms. We study the particular case of V = 0, if the transport terms are grater than the optical ones the non-equilibrium distribution function of the energy levels will follow the Fermi Dirac distribution function since the optical terms are a small perturbation. When the optical term increases, an equal generation of electron and holes appears. Both carriers tend to diffuse to the electrodes creating currents and, depending on the probability of these transitions, the Qd could be charged. Therefore, depending on the relation between the the transmission and transition probabilities the low and high injection level can be achieved.
When an external bias voltage is applied the distribution functions of the two leads differ and a net current appears. Overlapped to this, the optical processes appear filling/emptying levels and adding new conductive channels to the transport enhancing the total current. The results of the single Qd are shown in Fig. 3. We present two scenarios: a) the Qd symmetrically connected to the leads and (b) in asymmetric configuration. An interesting result appears for V=0. In the symmetrically coupled system (a), the current is zero since the incoming hole currents for each side equals the outgoing electron currents. Therefore, the net current is zero because the electron and hole currents compensate each other. This result derives intuitively from the rate equation, as it can be viewed as a sum of intensities in each energy level equals to zero. Therefore, in order to generate a net photo-current the symmetry of the Qd respect to the leads must be broken, hence, different coupling to the leads are needed as show in the asymmetric case. When an external bias voltage is applied, the transmissions coefficient between the Qd and the two leads change thus the system becomes asymmetric and a net current appears. The current peaks are related to the maximum transition probabilities for an incident photon reflecting the absorption spectra. When the voltage increases the current tends to be independent of hν. This effect is the result of competition between two processes, the pure light current term and the external bias voltage term. For small voltages, the optical terms dominate and the optical transition peaks are observed, but when the voltage increases the tunneling currents become the most important terms and the current appears as a photon energy independent. In the same way, in Fig. 4 we show the obtained current voltage curve I(V) under external illumination for the same previous system as a function of the energy of the incident photon. In Fig. 4 a-b) we show the current for the symmetric and asymmetric structure respectively. As it was demonstrated in previous works [19], the symmetry of the dark I(V) curves depends on the symmetry of the coupling to the leads, and current plateaus appear when conduction channels are opened [35]. For the illuminated case, the main differences appear for small voltages. The I(V) curve looses the step like behavior since the occupancy of the energy levels are not sequential with the applied bias voltage due to optical transitions. Moreover, the optical transitions tend to fill the conduction states while the valence states become emptied increasing the electron/hole currents. As a consequence, more conduction channels are opened and the total current is bigger than in the dark case. The next step is the study of a system composed by two Qd. We study two configurations that give us the main properties of bigger systems. In previous works [20], the parallel and serial configurations revealed the basic mechanisms behind the transport in large Qd arrays, therefore, these systems are the best starting point in order to extend the methodology to larger systems. Figure 5: The radii Qds are R 1 = 1.0nm and R 2 = 0.8nm respectively. The distances between the first Qd (Qd 1 ) and the leads are d L1 = 2.5nm and for the second one (Qd 2 ) d L2 = 2.5nm. The distance between the Qds is d 12 = 3nm. All the distances are measured from the center of the Qds. a) The total I(V) curve (in absolute value) for the parallel system in dark case and under different illumination conditions (incident photon energies). b) Photo-current as a function of the photon incident energy for different applied bias voltage.

Parallel case
We start with the case of two Qd in a parallel configuration. We use two Qd of different radius in order to obtain different optical transitions in each Qd for a given incident photon. Now, the rate equations for each energy level of the Qds have three contributions; the leads, the neighbor Qd and the optical contributions. In Fig. 5a) we show the I(V) curve under different incident photon energy. These photons energies are related to the maximum absorption peaks, for R = 0.8nm and R = 1.0nm Qds. In a previous study, we have demonstrated that in parallel configuration the total dark current trough the system is the sum of the individual Qd terms [20]. At low voltages, the optical terms dominate in the rate equation and the obtained I(V) curve differs from the dark curve. As shown before (for one single Qd), this trend is a competition between the optical transition probabilities and the transmission probabilities. For low voltages, the barriers are not transparent enough. Therefore, electron/holes are photo-generated in the Qd changing the distribution function of the Qd creating net fluxes incoming/outgoing from the Qd to the leads. The electrons tend to move from the Qd to leads and the holes follow the opposite direction, obtaining a net current if the currents are different. When the voltage increases, the barriers are bended and the tunneling probability increases as well. Then, the optical term becomes smaller than the electrical terms and the dark trend is recovered. When the energy of the incident photon increases, the optical transitions involve the higher energy levels (valence and conduction states for electron and holes respectively) which have the maximum transmission probabilities increasing the current.
Moreover, in this case, the electrons can also move to the other Qd. Thus, a current between Qds appears and breaks the symmetry respect to the leads obtaining a net current for the case V=0V (Fig. 5b). As in the previous case, when the polarization voltage increases the photo-current peaks tend to disappear. More peaks appear in the current term that take into account the photo-current generated in each Qd since the two Qds have different radii and the optical transitions occur at different photon energy. The third type of arrangement is two Qd in a serial configuration. In previous works, this kind of arrangement demonstrates the possibility of filtering the current trough the position of the energy levels [20,36]. Since overlapping between the energy levels of the Qds is necessary, only a few conduction channels are opened and the current is strongly dependent on the electrostatic coupling of the Qds obtaining NDR in the current-voltage curve [19,37]. We proceed in a similar way as in the previous case, the electrical response under an external bias voltage and under illumination with different photon energy. In this configuration, each Qd is connected to one electrode (left and right) and the neighbor Qd. First of all, we present the obtained dark I(V) and illuminated curve in Fig. 6. In the dark case, the NDR and current resonant peaks are obtained. The explanation related to these effects can be found elsewhere [37]. When the light is applied, the current increases, since the occupation of the higher energy levels increase, making favorable the tunneling processes. For certain photon energies, the current saturates in a voltage region as the transitions involves the higher/lowest energy states of the conduction/valence band. Overlapped to these optical transitions, we have the electronic transport due to the applied voltage. In order to obtain electronic transport between the Qds, the condition of the overlapping of the energy levels has to be fulfilled. Therefore the NDR and the current peaks still remain in the I(V) curve. In addition, we present the photo-current generated as a function of the energy of the incident photon in Fig. 6b). The values of the photon energy have been chosen in order to maximize the transitions probabilities in each Qd obtaining different photo generation rates in each dot. The cases under external bias voltages are also presented. We obtain a similar trend as in the previous cases, the current reflects the absorption spectra of the systems obtaining current peaks when the transitions probabilities are maximal in the cases when the photon energy equals the difference between the energy levels involved in the transition. As a difference as in the previous cases, the current saturates when the voltage is increased as a consequence of the no overlapping between the energy levels of the Qds. For these reason, we study this kind of system; a small Qd (with few energy levels) connected in series with a bigger one (with large number of energy levels). Therefore, the small Qd dominates the behavior of the photo-current since it controls the number of conduction channels. In this configuration, the photo-current retains the current peaks when the external voltage increases since the small Qd acts as a current filter.

Comparison with experiments
To conclude, we present the comparison between the presented theoretical model and a real system. The experimental results have been taken from Prins et al. [24]. They fabricated a system based in two electrodes separated by 4nm. A single layer of PeSe Qd of 4nm in size was deposited on top of the electrodes. Therefore, they obtained a system of Qds placed in parallel configuration. This kind of structure is the same as we have studied in the previous section. In order to simulate this new system, we have changed the material parameters to describe the PeSe Qds. The value of the effective masses and confinement potentials were taken from Cheng et al. [38]. We used m * Qd,CB = 0.07 and m * Qd,V B = 0.068 for electron and hole effective masses respectively. The confinement potentials were φ 1,ECB = 3.1eV and φ 1,HV B = 2.43eV . The dielectric constant value was ε r = 23ε 0 [39]. As an example, we created a system composed by three Qds placed in parallel configuration. With this small system, the main trends and the physics behind the device can be studied in a simple form. The Qd radius is generated randomly using a normal distribution with mean radius < R >= 2nm and σ = 0.1nm. The Qds were placed in parallel configuration but the distance between them and the leads was chosen as a free parameter breaking the symmetry of the system.
In Fig. 7 we present the comparison between the experimental and theoretical I(V) curves in dark and light conditions. The simulated curves follow the trends of the experimental ones in both cases, in dark and with light. The simulated curve still retains the step behavior but when the system becomes bigger and there is a random size distribution of Qds the current contributions of each Qd is summed and the plateaus disappear. Moreover, the theoretical curve reflects the asymmetry in the I(V). We also compute the normalized photo current in Fig. 8 and its comparison to Prins et al. Besides, we represent the total simulated absorption coefficient of the system. In both cases, the photo-current follows the same trend as the absorption coefficient but there are differences between the simulated and experimental curves. The differences appear since the absorption coefficients are not equal. Our model treats the PbSe Qd as a finite spherical potential well and uses the EMA approximation and it can be viewed as zero order approximation but other interactions effects must be included. Otherwise, the difference arises from the radius Qd distribution. If the Qd radius increases the energy gap decrease and the absorption spectra will tend to redshift, shifting the current to larger wavelengths. However, the methodology presented here can be improved using density of states and absorption coefficients that derive from powerful techniques as DFT.

Conclusions
The next generation of optoelectronic devices based on the singular properties of the Qds embedded in an insulator matrix poses new requirements on theoretical models in order to describe their response. We propose a theoretical framework based on the non-coherent rate equations combined with the transfer Hamiltonian approach to explain the transport properties of these kinds of systems in the selfconsistent field regime. Within this approach, the rate equations were modified to take into account the optical processes, stimulated generation and recombination. Moreover, the whole methodology presented here is only dependent on basic material parameters and the geometry of the device. This approach has been used to study the electrical response of three systems: one single Qd, two Qds in parallel configuration and two Qd in serial configuration. An important result emerges from the geometrical disposition of the Qds respect to the leads, obtaining that the generated photo-current is zero if the Qd is symmetrically coupled to the leads for the V=0 case. Moreover, the methodology has been used to compare with experimental results. The simulated I(V) curves in dark and under illumination conditions can reproduce the experimental observations, making possible the use of this model as a guideline to design new opto-electronic devices.

Acknowledgments
S. Illera is supported by the FI programme of the Generalitat de Catalunya. A. Cirera acknowledges support from ICREA academia program. The authors thankfully acknowledge the computer resources, technical expertise and assistance provided by the Barcelona Supercomputing Center -Centro Nacional de Supercomputación.