Universal scaling of quench-induced correlations in a one-dimensional channel at finite temperature

It has been shown that a quantum quench of interactions in a one-dimensional fermion system at zero temperature induces a universal power law $\propto t^{-2}$ in its long-time dynamics. In this paper we demonstrate that this behaviour is robust even in the presence of thermal effects. The system is initially prepared in a thermal state, then at a given time the bath is disconnected and the interaction strength is suddenly quenched. The corresponding effects on the long times dynamics of the non-equilibrium fermionic spectral function are considered. We show that the non-universal power laws, present at zero temperature, acquire an exponential decay due to thermal effects and are washed out at long times, while the universal behaviour $\propto t^{-2}$ is always present. To verify our findings, we argue that these features are also visible in transport properties at finite temperature. The long-time dynamics of the current injected from a biased probe exhibits the same universal power law relaxation, in sharp contrast with the non-quenched case which features a fast exponential decay of the current towards its steady value, and thus represents a fingerprint of quench-induced dynamics. Finally, we show that a proper tuning of the probe temperature, compared to that of the one-dimensional channel, can enhance the visibility of the universal power-law behaviour.


Introduction
Among all open problems in the field of quantum many-body systems, of special interest is the study of real-time dynamics far from equilibrium. In this context, recent state-of-the-art experiments performed on cold atoms [1,2,3,4,5,6,7] have shown the possibility to modulate in time various parameters and to detect transport properties as a useful tool to investigate the system response [8,9,10,11]. Due to their high degree of tunability, cold atoms represent an ideal platform to study quantum quenches [12,13,14], which consist in the rapid variation over time of one of the system parameters in a controlled way. Such protocol naturally settles an out-of-equilibrium state with highly non-trivial time evolution. Non-equilibrium physics of one-dimensional (1D) systems and time-resolved dynamics has been also recently investigated in pioneering experiments in solid-state implementation, see for example [15,16,17,18,19].
When dealing with out-of-equilibrium systems, two important questions arise: do such systems eventually settle to a steady state? And if so, what characterises their relaxation dynamics? Several theoretical studies [20,21,22,23] and experiments have addressed such topics, providing important results which strongly depend on the system considered. For example, some cold atom systems have shown a faster-than-expected relaxation towards a state well described by a grand-canonical ensemble [5]. Quantum quench experiments performed splitting a 1D gas of cold atoms [4], initially prepared in a thermal state, have shown the system to settle towards a quasi-thermal state, which can be well-described in terms of an effective temperature lower than the initial one. On the contrary, other cold-atom implementations have reported relaxation towards a steady state which cannot be described as a thermal one [2,24]. In this respect, integrable systems are particularly interesting, since it has been conjectured that they approach a steady state that can be characterised by the so-called generalised Gibbs ensemble [25,26], whose associated density matrix is in general very different from the thermal one. However, a complete characterisation of their relaxation dynamics is yet to be found.
In a recent paper [45], the relaxation dynamics of a Luttinger liquid subjected to a sudden quench of the interaction strength has been studied at zero temperature. It has been shown that the latter generates entanglement between counter-propagating excitations [45,57], resulting in finite bosonic cross-correlators among the 1D channel. This, in turn, affects the time evolution of fermionic observables such as the charge and the energy currents locally injected from a biased probe. Indeed, in the long-time limit they display a universal power-law relaxation dynamics ∝ t −2 , not sensitive to the details of the quench.
In this paper, we investigate whether this universal power law is still present when the system is prepared, before the quench, in a thermal state. We therefore study the behaviour of the fermionic non-equilibrium spectral function and its relaxation dynamics toward the steady state at finite temperature. We demonstrate that at long times its universal quench decay ∝ t −2 is stable against thermal effects, while the non-universal interaction-induced power laws, present at T = 0, acquire a fast exponential decay dictated by the finite initial temperature. The initial preparation in a thermal state can then be useful to highlight and enhance the visibility of universal features in comparison to the zero temperature case. We also consider the impact of this quench protocol with finite temperature preparation on transport. The behaviour of the spectral function, indeed, determines the relaxation dynamics of the charge current injected from a biased probe, which displays a universal scaling law ∝ t −2 , again robust against thermal effects. This is in sharp contrast with a conventional non-quenched situation, in which a finite temperature induces a fast, non-universal exponential decay towards the steady current value. We also consider the case of different initial temperatures for the 1D channel and the probe, showing that a proper tuning of the latter can enhance significantly the visibility of the ∝ t −2 universal decay in the current dynamics.
The outline of the paper is the following. In Sec. 2 the model is described. In Sec. 3 the fermionic Green functions, central to the evaluation of the transport properties, are introduced. Their basic building blocks, the bosonic two-point correlation functions are also defined and evaluated, and their time dynamics is analysed in detail. In Sec. 4 we focus on the non-equilibrium spectral function and its relaxation dynamics towards the steady state. In Sec. 5 we consider a possible transport setup, where a probe is assumed to be weakly tunnelcoupled to 1D channel, and we evaluate the resulting charge current and its time-evolution after the quench. Section 6 contains the conclusions.

Model
We consider an interacting 1D channel of spinless fermions with short-range repulsive interactions. At times t < 0, they are governed by the Hamiltonian (in all the paper we set = k B = 1) with v the Fermi velocity, ψ r (x) the fermion field of the r-channel (r = L/R for left-/rightbranches), and ϑ R/L = +/−. The interaction term is given by where n r (x) =: ψ † r (x)ψ r (x) : is the particle density on the r-channel and g (2) i , g (4) i are the coupling strengths of the intra-and inter-channel interactions, respectively [27,28]. For simplicity, hereafter we will consider the Galilean invariant case g This fact has no relevant consequence, as the more general case g (2) i = g (4) i leads to the same qualitative results discussed here. The non-interacting part of H i can be written in terms of free bosonic fields φ r (x), related to the fermionic field operator by the bosonization identity [27,28] where q F is the Fermi wave vector and a is the usual short-length cut-off properly introduced in bosonization techniques [28]. Here, we have safely omitted Klein factors. Different regularisation schemes are possible and have been considered in the literature, see for example the constructive bosonization approach described in [58,59]. The Hamiltonian in Eq. (1) can then be diagonalised introducing new bosonic operators φ i,± (x), connected to φ r (x) by the canonical transformation with the Bogoliubov coefficients 2B a parameter encoding the repulsive interaction strength, with 0 < K i ≤ 1 and K i = 1 for non-interacting fermions. The Hamiltonian H i then reads where v i = v/K i is the renormalised propagation velocity of the collective excitations described by the chiral fields For t < 0 the channel is prepared in a thermal state described by the equilibrium density matrix ρ eq = Z −1 e −H i /T , where Z = Tr{e −H i /T } is the partition function. At t = 0, the bath is disconnected and the interaction strength is suddenly quenched Here, v f = v/K f and K f = [1 + (g f /πv)] −1/2 describes the interaction strength in the postquench state. Note that also these fields obey chiral properties, and they are connected to φ i,η (x, t) by the canonical transformation [44,51] and to the free bosonic fields by Note that the non-quenched case ( Finally, the different Bogoliubov coefficients are related among themselves via B ± = A ± θ + + A ∓ θ − .

Fermionic and bosonic correlation functions
In the following we will discuss the local lesser fermionic Green function of the 1D channel at the generic position z. Here, the brackets . . . eq denotes a quantum average performed on the initial thermal density matrix ρ eq and the last bosonic expression is obtained using the identity of Eq. (3). The first step to evaluate the Green function is to compute the time evolution of the free bosonic fields φ r (z, t). It crucially depends on whether the operators are evaluated before or after the quench. Indeed, using Eqs. (4) and (12) one can write having exploited the proper chirality properties of Eq. (7) and Eq. (9). Note that, since space translational invariance is not broken by the quench protocol, the Green function will not depend on the generic position z. On the other hand, it will feature four different time regimes, depending on the sign of t andt, due to the breaking of time translational invariance. One can thus write where ϑ(t) is the unit step function. If both time arguments are negative (n), the Green function is not affected by the subsequent quench and it thus reduces to the standard equilibrium result Here, we have introduced the equal time bosonic correlator whose expectation value on the initial thermal state is [28] C with Γ(z) the Euler Gamma function and ω i = v i a −1 a cut-off frequency. From now on, we will safely set a −1 = q F . It is worth to notice that G nn r (t,t) depends only on the time The opposite regime (t < 0,t > 0) can be easily obtained by exploiting G np r (t,t) = −G pn r (t, t) * . Finally, the fourth regime, when both t andt are positive, is the most interesting one and can be expressed as Indeed, as we will see, it is the only one responsible for the presence of the universal features which characterise the post-quench relaxation dynamics. Moreover, it is the only one controlling the transport properties after the quench. It is therefore worth to analyse it in details, highlighting the physical origin of its peculiar quench-induced features. To this end, it is useful to rewrite it in terms of correlation functions between the final chiral bosonic fields. This reads with Here α, β = ± and the time difference should satisfy τ ≡ t −t < t. Exploiting the chirality of the fields as well as Eqs. (10) and (18), these correlation functions can be decomposed as with D (0) α,β (t, τ ) the zero-temperature contribution and ∆D α,β (t, τ ) corrections due to the finite-temperature of the initial state. In particular, one obtains with defined in Eq. (11). As a general feature, we note that auto-correlators (α = β) only depend on τ and not on t. In contrast, cross-correlators (α = −β) feature a full dependence on both t and τ . A direct comparison between the quenched and the non-quenched cases turns out to be particularly useful. In the latter one has θ − = 0 and cross-correlators vanish, This is a consequence of the fact that chiral fields φ f,η (x, t) are completely decoupled from each other in the Hamiltonian H f . Conversely, we observe that a quantum quench always leads to finite cross-correlations, see Eqs. (26) and (27).
To study the relaxation dynamics it is useful to derive asymptotic expressions for D α,−α (t, τ ). For the zero-temperature contribution of D (0) α,−α (t, τ ) the relevant time scale is τ and for t τ one finds via [45] Inspecting the temperature-dependent term ∆D α,−α (t, τ ), an additional time scale ( T ) −1 emerges and two different regimes of the gamma functions, present in Eq. (27), have to be considered where Γ(n, z) is the n-th order polygamma function. Inserting Eq. (30) in the expression for ∆D α,−α (t, τ ) of Eq. (27) and considering the reasonable temperature regime T ω i , one obtains the leading terms The validity of this asymptotic expansion is confirmed in Fig. 1(a), where we have numerically evaluated ∆D α,−α (t, τ ) as a function of time t for different temperatures and at fixed τ . At short time t ( T ) −1 the thermal component of the cross-correlator saturates to a timeindependent value, while a power law decay ∝ t −2 is evident for t ( T ) −1 , in accordance with Eq. (31).
Combining Eqs. (29) and (31) we obtain the asymptotic expression for the full crosscorrelator As one can see there are two different regimes, both featuring a power-law decay ∝ t −2 but with different prefactors. For long time t ( T ) −1 the cross-correlator is positive and proportional to the temperature. By contrast, when finite temperature effects have not yet kicked in, i.e. at shorter time t ( T ) −1 , the prefactor of the t −2 power-law decay has a negative prefactor and is temperature-independent. This behaviour emerges clearly in Fig.  1(b) where we plotted |D α,−α (t, τ )| for fixed τ . The transition between the two regimes spans an order of magnitude.

Non-equilibrium spectral function
To fully characterise the effects of the quench, we now focus on the local (lesser) nonequilibrium spectral function. This is a key quantity to inspect the presence of universal features in the relaxation dynamics and where the role of finite temperature plays a nontrivial role. Moreover, it is a important ingredient to evaluate transport properties. The local non-equilibrium spectral function is defined as [48] Our task is to investigate its time evolution after the quench, i.e. for t > 0. Since the integration range over τ extends to +∞, the calculation of the above expression requires to distinguish between two different regimes of the Green function: Equation (34) is exact and suitable for numerical investigation, using Eqs. (19,20). This quantity admits a finite steady stateĀ r (ω) ≡ lim t→∞ A < r (ω, t) given bȳ with It is instructive to derive some helpful analytical expansions. In particular, we can observe that the integrand e iωτ G pp r (t, t − τ ) contributes to A pp r (ω, t) only in two distinct regions: near τ ∼ 0, where the Green function approaches its poles -see Eqs. (20,24) -and close to the boundary of the integration domain τ ∼ t. As a consequence, in the long-time limit t ω −1 one has [45] A pp the former term stemming from an expansion of the integrand for τ → 0 and the latter from an expansion for τ → t. In particular, using Eq. (31), one finds where Considering A r (ω, t), it is useful to introduce z = t − τ and expand the integrand retaining only the leading term in z/t. The result reads with Here ω = 0.1ω f , K i = 0.9, K f = 0.6, time units ω −1 f and temperature units ω i . We have included the zero temperature case (green curve) to underline the difference with the finite temperature case.
Using Eq. (30), two asymptotic regimes can be identified a non-universal exponent encoding the strength of the interaction quench which, for reasonable values of the interaction quench, can take values [45] 1 < ξ ≤ 2. The first regime, which is present as long as t ( T ) −1 , is a typical interaction-dependent power law. On the other hand, for t ( T ) −1 , a fast exponential decay appears which basically kills A r (ω, t) in the long-time limit.
The same arguments can be used to derive the asymptotic behaviour of A pn r (ω, t) present in Eq. (34). In this case, the integrand contributes only for τ ∼ t and τ ∼ t where the contribution stemming from the expansion for τ → t is In the latter equation we have introduced the function As one can see, A pn r (ω, t) qualitatively behaves like A r (ω, t): at long times t ( T ) −1 it becomes negligibly small due to the presence of a fast non-universal exponential decay induced by the finite temperature.
In the end, the transient dynamics of the non-equilibrium spectral function is dominated by where ξ, given in Eq. (43), attains values 1 < ξ ≤ 2 for reasonable quenches. At short times (t ( T ) −1 ) a non-universal power law behaviour with exponent ξ is present and could dominate the relaxation dynamics, masking the universal features. For sufficiently long times t ( T ) −1 the situation changes and the non-equilibrium spectral function presents a very clear t −2 decay, directly linked to the interaction quench. This is an interesting feature, since in a standard equilibrium case (non-quenched system) one would expect a exponential behaviour for fermionic correlation functions at finite temperature. It is also in contrast with the zero temperature case discussed in Ref. [45] and shown with green curves in Fig. 2 for reference. At T = 0, indeed, the non-equilibrium spectral function always shows non-universal power laws ∝ t −ξ , masking the ∝ t −2 behaviour. On the contrary, the initial preparation in a thermal state allows to identify of ∝ t −2 features in the relaxation dynamics: provided that t ( T ) −1 this universal power-law decay, a fingerprint of the quench-induced entanglement, remarkably emerges in the non-equilibrium spectral function at long enough times.
All features described so far can be seen in Fig. 2, where the real and the imaginary part of the transient spectral function δA r (ω, t) are evaluated numerically as a function of t for different temperatures. Notice that we have also inserted the zero temperature case (green line) to better clarify the differences at finite temperature. The universal decay ∝ t −2 clearly emerges in the real part once the exponential decay of the non-universal contributions sets in, around t ∼ 5( T ) −1 . A similar behaviour is observed looking at the imaginary part of the non-equilibrium spectral function even if, in this case, the power-law decay is faster and ∝ t −3 .

Transport properties
In this section we focus on a possible setup to observe the intrinsic features present in the fermionic spectral function induced by the quantum quench. To this end we focus on the relaxation dynamics of the average current injected into the channel from a weakly tunnelcoupled 1D non-interacting probe. It is worth to underline that the probe is a tool to inspect the intrinsic properties of the fermionic channel out-of-equilibrium and, thus, it is supposed to be as non-invasive as possible. A perturbative approach in the weak tunnel coupling is therefore fully justified in evaluating the related transport properties [60,61]. Fermions of the probe are described by the Hamiltonian with χ(x) the fermionic field. The probe is initially (t < 0) prepared in equilibrium at a given temperature T p (not necessarily equal to T ) with associated density matrix ρ (p) eq = e −Hp/Tp /Z p with Z p = Tr[e −Hp/Tp ]. It is also subjected to a bias voltage V taken with respect to the Fermi level of the channel.
A localised injection at x = x 0 is switched on at t = 0 + -immediately after the quenchand is described by the tunneling Hamiltonian [44,45,63,64,65] with λ the tunneling amplitude and ϑ(t) the step function. We note that after the tunneling process both the channel and the probe will evolve dynamically and can be considered as a closed system. We are interested into the total injected current I(V, t) = I + (V, t) is the chiral current associated to the η channel, with q the fermion charge. Here, where ρ tot (0) = ρ eq ⊗ ρ (p) eq and its time evolution is computed in the interaction picture with respect to H t . The chiral particle density reads At the lowest order in the tunneling amplitude, the total current can be written as [45] where G < r (t, t − τ ) is the lesser Green function of the fermionic channel defined in Eq. (13) and is the greater Green function of the probe with ω p = v/a the associated energy cut-off. Recalling that the lesser Green function can be expressed as G < r (t, t − τ ) = +∞ −∞ e −iωτ A < r (ω, t)dω by Fourier transform, one immediately recognise the direct link with the non-equilibrium spectral function of the channel discussed in the previous Section. Plugging Eq. (21), with D α,β (t, τ ) given in Eqs. (24)(25)(26)(27), and Eq. (55) into Eq. (54) one can eventually write the total current as where ν ± and γ are defined in Eqs. (36) and (39), respectively. Equation (56) contains both the initial temperatures of the probe T p and of the channel T . Useful analytical expansions in the long-time limit can be obtained following the same procedure described in the previous Section. In particular, in the long-time limit the current decomposes into two terms the former stemming from an expansion of the integrand for τ → 0 and the latter from an expansion for τ → t [45]. I (1) (V, t) is responsible for the universal features and reads where with κ(t) given in Eq. (39). The steady-state current is Looking at the term I (2) (V, t), it is useful to introduce z = t − τ and expand the integrand of Eq. (56), retaining only the leading term in z/t. One thus obtains where is an oscillating function of time with period ∝ (qV ) −1 . Exploiting the asymptotic expression for the gamma functions in Eq. (30), one finds Here, are power-law corrections to the exponential decays in time. As one can see I (2) (V, t) contains all the non-universal contributions.
To discuss the analytical results obtained so far, we begin by recalling the T = T p → 0 limit. In this case Eq. (58) with κ(t) = −1/2 and Eq. (62a) contribute to I(V, t): the current decay shows then a competition between the t −2 power law and a non-universal power law ∝ t −ξ−1 . Since for reasonable quenches ξ + 1 2, the two power laws are similar: although the universal decay is leading, distinguishing its contribution in I(V, t) can be cumbersome. Let us now consider the case of a thermal preparation for the channel and the probe starting with the case T = T p > 0, considering quenches which increase the interaction strength, i.e. K f < K i (implying > 1). Here, two different time regimes must be distinguished: • For t ( T ) −1 , the same competition between the universal ∝ t −2 and the non-universal ∝ t −ξ−1 power laws as in the T = 0 case is present; contribute to the current. In both cases, the non-universal power law is now exponentially suppressed and the decay of I(V, t) is entirely dominated by a universal power law t −2 with a prefactor ∝ T .
Thus we can conclude that a thermal preparation with a common temperature T for the system strongly enhances the visibility of the universal decay, setting a time scale ( T ) −1 beyond which the latter becomes the only relevant contribution to I(V, t).
The situation is even better if the case T p > T > 0 is considered. Now, three regimes arise Thus, by suitably choosing a larger temperature T p for the initial preparation of the probe with respect to the initial temperature T of the channel, the universal power law is the only leading term beyond the time scale T −1 p and two regimes with universal decay but different prefactors are now visible (see the last two points above). We close this discussion by commenting the non quenched case (γ = 0, = 1) with T p > T > 0: without quench the universal power is not present -see Eq. (58). For t T −1 p the non-universal power law ∝ t −ξ−1 is found, while for t T −1 p the current decays exponentially as normally expected. Thus, the presence of a universal decay is a unique fingerprint of the quantum quench, stable even in the case of a preparation at a finite temperature.
All these features are confirmed by the numerical evaluation of Eq. (56). In Fig. 3 (main panel) the quantity |δI(V, t)| = |I(V, t)−Ī(V )| is plotted as a function of t, on a bi-logarithmic scale, for different initial temperatures of the whole setup, with T = T p . The two regimes t ( T ) −1 and t ( T ) −1 are clearly visible: the former features oscillations as a result of the strong competition between the universal power-law decay and the non-universal one which, for the quench considered here, has a non-universal exponent −2.3. By contrast, the latter regime shows a very clean decay ∝ t −2 , resulting from the exponential suppression of the non-universal oscillating terms. The crosses in Fig. 3 (main panel) show |δI(V, t)| in the non-quenched case with K i = K f : the fast exponential decay for t ( T ) −1 is evident (note that, to improve the readability of the Figure, the data sampling rate is kept very low). In both Panels, the quench is from K i = 0.9 to K f = 0.6.
A more detailed plot of the non-quenched current is shown with a logarithmic scale in the Inset, where the oscillations of the non-universal decay can be clearly seen enveloped by an exponential decay, which is in good agreement with Eqs. (62c,62d) -see the dashed black line.
The t ( T ) −1 regime is the most interesting one in order to assess the different time evolution of quenched and non-quenched channels at finite temperature with T = T p . It is therefore worth to discuss how the prefactor F(V, T ) of the universal power law, introduced in Eqs. (58) and (60), depends on the temperature T and on the bias V in this regime. In Fig. 4(a), we plot |F(V, T )| focusing on a reasonable range of the parameters. Interestingly, this function turns out to be non-monotonic with respect to both T and V . This fact leads to a non-trivial feature: considering temperatures lower than 10 −2 ω i , a decrease of the voltage bias can actually determine an increase by orders of magnitude of the prefactor |F(V, T )|. As an example of this behaviour, in Fig. 4(b) we plotted the current |δI(V, t)| as a function of time for two different voltage biases but with the same temperature T = 10 −3 ω i . In the t ( T ) −1 regime, the dotted line (qV = 0.004ω f ) lies more than one order of magnitude above the solid one (qV = 0.1ω f ). Therefore, a proper tuning of the voltage bias can be useful to magnify and to detect the peculiar features induced by the presence of the quantum quench.
In Fig. (5), we show |δI(V, t)| for the case T p T > 0. The two regimes T −1 p t ( T ) −1 with |δI(V, t)| ∝ 1/(2t 2 ) and t ( T ) −1 with |δI(V, t)| ∝ T /(ω i t 2 ) are clearly evident, as well as the different prefactors that characterise them. This confirms that preparing the probe in a state with a higher temperature strongly enhances the visibility of the universal decay of the current. Plot of |δI(V, t)| (units 4q|λ| 2 ω −1 f a −2 ) as a function of time with K i = 0.9 and K f = 0.6. The temperatures of the channel and of the probe are T = 2 · 10 −4 ω i and T p = 2 · 10 −1 ω i , respectively. The black lines shows decays ∝ t −2 . Here qV = 0.1ω f , time is in units ω −1 f .

Conclusions
We have investigated the relaxation dynamics of a 1D channel of spinless interacting fermions, initially prepared in a thermal state at T > 0, subject to a sudden quench of the interaction strength. The channel is modeled as a Luttinger liquid. The out-of-equilibrium fermionic spectral function is analytically evaluated and studied in details. One of the main results is that it shows a universal t −2 decay towards the steady state in presence of finite temperature, showing the peculiar role associated to the quantum quench and its link with universal powerlaw decay. Indeed, at finite temperature, non-universal factors acquire an exponential decay, greatly enhancing the visibility of the decay t −2 in the long time limit. This is in sharp contrast with the expectation in absence of an interaction quench, where fermionic correlation functions would acquire an exponential behaviour due to the presence of a finite temperature. These features are also reflected in observable quantities such as charge current. In particular, we have analysed the charge current flowing to the channel from a locally tunnel-coupled biased probe, prepared before the quench in a thermal state with a temperature T p . A weak tunneling is switched on immediately after an interaction quench, and the evolution of the closed system composed by channel and probe is analysed to lowest order in the tunneling amplitude. The current decays towards its steady value with the same universal power law ∝ t −2 which characterises the universal relaxation of the non-equilibrium spectral function. Not only the persistence of this universal decay confirms the robustness of the results previously obtained in the T → 0 limit [45], but finite temperature can be a tool to enhance the visibility of this universal scaling induced by the quench. Finally we have shown that, by acting on the initial temperature of a tunnel-coupled probe, the visibility of the ∝ t −2 decay can be enhanced even further. For a typical channel of cold atoms [4], one can estimate temperatures of about T ∼ 100 nK, for which the thermal time-scale is about T −1 ∼ 500 µs. Typical experiments in such a system can explore time domains well in excess of the hundred of ms, so it seems feasible, at least in principle, to observe the striking differences between quenched and non-quenched systems in the presence of thermal effects.
A. C. and T. L. S. acknowledge support from the National Research Fund Luxembourg (AT-TRACT 7556175). M. C. acknowledges support from the CNR-CONICET cooperation programme "Energy conversion in quantum, nanoscale, hybrid devices". A. C., F.M. G., F. C. and M. S. gratefully acknowledge support from UNIGE -Genova University.