{'fraction_non_alphanumeric': 0.04137552405669794, 'fraction_numerical': 0.01357556398482731, 'mean_word_length': 4.18006721820062, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 17, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-linear transport in an out-of-equilibrium single-site Bose Hubbard model: scaling, rectification and time dynamics

Recent experiments in hybrid-quantum systems facilitate the potential realization of one of the most fundamental interacting Hamiltonian-Reservoir system, namely, the single-site Bose-Hubbard model coupled to two reservoirs at different temperatures. Using Redfield equations in Born-Markov approximation, we compute non-equilibrium average particle number, energy and currents beyond linear response regime, both time-dynamics and steady state and investigate its dependence on various tunable parameters analytically. We find interesting scaling laws in high temperature regimes that are independent of choice of bath spectral functions. We also demonstrate that the system shows very interesting particle and energy current rectification properties which can be controlled via the relative strength of interaction and temperatures, as well as via the degree of asymmetry in system-bath coupling. Specifically, we find inversion of direction of thermal rectification as a function of the relative strength of the interaction strength and the temperatures. We also show that, in the limit of low-temperature and high interaction strength our results are consistent with the non-equilibrium spin-Boson model. Our results are experimentally relevant not only to hybrid quantum systems, but also in other areas such as molecular junctions.


A. Introduction
Far-from-equilibrium systems beyond the paradigm of linear response has been a subject of growing theoretical and experimental interest [1][2][3][4][5][6][7][8]. Such studies are important from fundamental perspective and from the point of view of device applications. Especially, understanding the role of interactions is of great interest. Important fundamental problems to explore include, role of far-from-equilibrium physics in physical quantities of interest and the intricate interplay between system and bath degrees of freedom when interactions are involved. Some of the main consequences of interactions range from fundamental physics (such as the phenomenon of Kondo effect [9,10]) to device applications (such as diodes/rectifiers [11]). Even though the role of interaction in device applications in fermionic systems is well studied [12][13][14], the corresponding bosonic set-ups are much less explored. However, with the advent of the relatively new field of photonics, bosonic set-ups for device applications such as a potential optical diode are gaining popularity [15][16][17]. With recent cutting-edge technology in hybrid quantum systems, it is now possible to design interacting bosonic Hamiltonians and reservoirs [8,[18][19][20][21][22][23]. There is progress not only in fabricating Hamiltonian-Reservoir systems, but also, great advance in measuring physical quantities of interest such as photon number, photon statistics [24] and photon current [25] .
One of the most basic interacting bosonic Hamiltonian system that one can think of is a single bosonic site with Bose-Hubbard interaction, hereafter, called the singlesite Bose Hubbard (SSBH) model. In optics, such interaction is often called the Kerr interaction strength. In the field of hybrid quantum systems a Hamiltonian with such an interaction can be potentially experimentally realized in more than one way. The Jaynes-Cummings model, an experimentally realized light-matter system, can be tuned to the dispersive regime where it behaves like a SSBH model. The photon-spin interaction in this limit can be integrated out "pertubatively" to generate a non-linear Bose-Hubbard like interaction between photons [19,26]. Alternatively there are interesting potential realizations of Bose-Hubbard interactions between photons which involve 4-level atoms in an optical cavity [27,28]. These realizations offer large tunability of parameters. Specifically, while the former realization involving Jaynes-Cummings model has a small interaction strength compared to the linear term, the later realizations have very large interaction strength.
Another area of applicability for such interacting bosonic Hamiltonians, which is perhaps more suited to non-equilibrium measurements, are the fields of molecular thermoelectrics and nanophononics. These fields typically deal with molecular junctions connecting two reservoirs, experimentally, which are often two large chemical compounds [29][30][31][32][33][34]. Role of interactions in phononic transport through such systems is of interest both experimentally and theoretically [35].
There has been a large amount of work on such a SSBH model with a finite interaction strength coupled to a single bath [36][37][38][39][40][41]. On the other hand, non-equilibrium spin boson (NESB) model, which corresponds to the limit of very large interaction strength, has been well studied in out-of-equilibrium set-ups [10,[42][43][44][45][46] and is of growing experimental significance. The conductance of an anharmonic junction with quartic anharmonicity has been studied recently [47,48]. However, there is essentially no investigation of the SSBH model with finite interaction strength in a far-from-equilibrium setup via connection with multiple reservoirs. In this paper, we investigate the SSBH model weakly coupled to two bosonic baths at different temperatures beyond the linear response regime.
The organization of our paper is as follows: In section B, we describe our setup. In section C, we give the Redfield quantum master equation for our setup and in section D, we find the non-equilibrium steady state (NESS) properties. In section E, we look at scaling behaviour of average particle number and energy at NESS and in section F, we discuss the scaling behaviour and rectification of NESS particle and energy currents. In section G, we present results for time dynamics of various physical quantitites which has become of growing recent interest. Finally, in section H, we summarize our main results along with an outlook.

B. Model and Setup
We consider a single site with Bose-Hubbard interaction connected to two bosonic baths in non-equilibrium. The baths are taken to be quadratic and the systembath couplings are taken to be bilinear. Thus our set-up is given by the full system+bath Hamiltonian whereâ correspond to bosonic annihilation operators and B r to those of th bath ( = {1, 2}). Note thatN = a †â is a conserved quantity with respect to the system Hamiltonian H S . The energy spectrum of the system Hamiltonian can be easily written down and it has the non-linear form E = Ω 0N + χN 2 . The baths are quadratic and have infinite degrees of freedom. ε is a dimensionless parameter that controls system bath coupling. We assume that, initially, there is no coupling between the system and the baths, and the two baths are at thermal equilibrium with their own inverse temperature β 1 and β 2 and chemical potential µ 1 and µ 2 . Thus the initial bath correlation functions satisfy the thermal properties : where n (ω) = [e β (ω−µ ) − 1] −1 is the bosonic distribution function. We also introduce the bath spectral functions: We will assume, i.e, the two baths have same density of states, but the system-bath coupling is different in general. Γ has dimensions such that J (ω) has dimensions of energy. The quantitative nature of the results will depend on our choice of spectral function. We assume a general spectral function which is commonly used in bosonic systems: where ω c gives the cut-off frequency. The cut-off frequency is considered very large so that the system energy levels near the edge of the bath spectrum correspond to extremely high energies, which do not really contribute to the system properties at the chosen set of temperatures and chemical potentials. This is satisfied when ω c Ω 0 , χ, β −1 1 , β −1 2 , µ 1 , µ 2 . Also, we are concerned with a photonic or phononic system, so we will set µ 1 = µ 2 = 0 finally.

C. The Quantum Master Equation (QME)
We want to investigate this set-up without any restrictions on the interaction strength χ. So, we adopt the method of QME under Born-Markov approximation, which is only valid when system-bath coupling is weak, and gives results only to leading order in system-bath coupling. The microscopically derived Redfield QME for our system coupled to single bath was first written down in [37]. It is straighforward to generalize to two baths. Let us define ρ ≡ T r B (ρ f ull ) with ρ f ull being the full density matrix of system+bath and T r B (..) implying trace taken over bath degrees of freedom. The weak-coupling Redfield equation for the density matrix of the system ρ of our set-up is given by : wherễ ω = Ω 0 + χ(2N + 1), and h.c. stands for Hermitian conjugate. Note that the Redfield equation assumes ε 1 and keeps terms only upto O(ε 2 ), i.e, only upto quadratic in system-bath coupling (the Born approximation). In deriving the above equation, we have also done the Markov approximation which assumes that the observation time is much larger than the time scale of relaxation of bath correlation functions. To ensure weak system-bath coupling, we choose ε = 0.1 and κ r to be of the same order as Ω 0 .
In equilibrium, i.e, when β 1 = β 2 = β, µ 1 = µ 2 = µ, it can be checked by direct substitution that the thermal state is the steady state of the Eq. (6), where Z = T r e −β(H S −µN is the equilibrium partition function.

D. The Non-equilibrium Steady State (NESS)
To calculate various physical observables in NESS beyond linear response, we need to find the NESS density matrix. Also, in non-equilibrium, since we have no guess for direct substitution, we need to find the ρ directly from Eq 6. Since both the number operatorN and the system Hamiltonian H S are diagonal in occupation number basis, the NESS transport properties as well as average occupation and energy of the system can be found from the steady state diagonal elements of ρ in this basis. The occupation number basis satisfieŝ The evolution equation for the diagonal elements ρ n = n|ρ|n is given by where n C ( ) n = (n + 1)J (ω n )n (ω n ), D ( ) n = nJ (ω n−1 )(n (ω n−1 ) + 1) ω n = Ω 0 + χ(2n + 1) .
In the steady state we set the LHS of Eq. 10) to zero. This leads to a difference equation. Noting that C ( ) = D 0 = 0, we obtain, by recursion, the solution ρ n D n = ρ n−1 C n−1 (12) The constant ρ 0 is fixed from the normalization condition n ρ n = 1, i.e, trace of density matrix is unity. Thus, whereZ is a normalization constant and is analogous to the partition function in equilibrium systems. It can be easily checked by putting n 1 (ω) = n 2 (ω) = n(ω) in Eq. (13) that ρ n indeed gives Eq. (8) in equilibrium. Also, in the equilibrium case, Eq. (12) corresponds to the detailed balance condition. In the non-equilibrium case, while this still looks like the detailed balance condition, we note that in general it is not possible to define an effective temperature. Also note that ρ n is independent of the choice of spectral function of the bath.
The explicit expression for population ρ n of bosons in NESS, given by Eq. 13, is the central result that allows us to go beyond linear response in this interacting bosonic problem. The population, in itself, is a physically measurable quantity and, as we show below, can be used to compute various other physical observables.
In Fig. 1, we show the plots of population density ρ n of the system under asymmetric system-bath coupling. Since system-bath coupling is asymmetric, the out-ofequilibrium distribution changes under interchange of hot and cold baths. For high interaction strength, i.e, for χ Ω 0 , T 1 , T 2 , only the lowest two levels have nonnegligible probability. In this regime, we can truncate the energy spectrum in just two levels. Then ρ 0 and ρ 1 become with ω 0 = Ω 0 + χ. The above results are exactly the same as obtained for a non-equilibrium spin-boson model (NESB) by using the Redfield equation [42]. Thus for χ Ω 0 , T 1 , T 2 , the system becomes identical to the NESB. Since NESB is already a well explored problem, in the following, we will be mainly interested in the physics beyond this regime. We also note that for χ = 0 the system reduces to a harmonic oscillator and in this case the population ρ χ=0 n is given in terms of an effective temperature which is consistent with the finding in Ref. [49].

E. Average occupation and energy
First we will look at the average occupation and energy of the system. These quantities are measurable in current state-of-the-art experiments in quantum light-matter hybrid systems [19,24]. The expressions for these are given by with E n = Ω 0 n + χn 2 .
In NESB limit, these average quantities can be trivially found from Eq. 15. They become N ≈ ρ 1 , and Ĥ S ≈ ω 0 ρ 1 . We are interested in going beyond the NESB regime. So let us look at the regime of high temperatures T 1 , T 2 χ, Ω 0 . First, we look at the normalization constant defined in Eq. 14. which can be written asZ Note that f (ω p ) > 1, and hence log f (ω p ) > 0. It follows that there is a energy level cut-off n * beyond which the energy levels have negligible contribution toZ. For high enough temperatures, we can assume, β ω p 1 ∀p < n * . Under this condition, we can expand f (ω p ) to obtain (after some amount of algebra), The contribution of terms n > n * is small and so their precise form is irrelevant. Hence we get Thus, for T 1 , T 2 Ω 0 , χ the normalization constant has the same form as the equilibrium partition function with the effective temperatureT . This is consistent with the effective temperature for harmonic oscillator (χ = 0) given in Eq. 16. For high temperatures, β ef f = 1/T . It is also interesting to note that for symmetric systembath coupling, i.e, Γ 1 = Γ 2 , the effective temperature is just the mean temperature,T = T m = (T 1 +T 2 )/2. However, the description in terms of an effective temperature is not possible for low temperatures, except when χ = 0.
The high temperature scaling of the normalization constant can now be easily found by noting that for high temperatures, the summation can be converted into an integral. So we havẽ The second step requires the conditionT χ Ω 2 0 , under which we see that the normalization constant scales as T /χ.
The above trick can be used to find high temperature scaling of average of any operator which is diagonal in the eigenbasis of the system Hamiltonian. The average of any operatorĥ which is diagonal in the eigenbasis of the system has the form which, for T 1 , T 2 Ω 0 , χ, exactly following above arguments, becomes Using this, we readily obtain the high temperature be- Thus N should show a data collapse for various χ and vary as a function of T /χ, whereas, Ĥ S /χ varies as T /χ and should show a data collapse for χ Ω 0 . The dotted plots correspond to the NESB model. For χ Ω0 (Ω0 = 1), there is data collapse over the entire range of temperature. Even for small χ(= 0.1), there is substantial deviation from linear (χ = 0) behaviour. For high temperatures, N scales as T1/ω0, while, Ĥ S scales as T1.
Here the system-bath coupling is taken symmetric : Γ1 = Γ2 = 1. All energy variables are measured in units of Ω0.
To check the above high temperature discussion and the connection with equilibrium behaviour, we now first define which quantifies the degree of deviation from equilibrium, r = 1 corresponding to equilibrium. To ensure that the system is far from equilibrium, we keep r fixed at r < 1. We note that, if r is kept fixed, all the NESS results become a function of only one temperature, say, T 1 . We choose T 1 /ω 0 , where ω 0 = Ω 0 + χ, as the scaling variable since it can be used both in the highly interacting regime (where χ ≈ ω 0 ) and the linear regime (where χ = 0). Also in the NESB regime, with r fixed, all NESS quantities vary as a function of T 1 /ω 0 , as can be checked from Eq. 15. In Fig. 2, we plot N and Ĥ S /ω 0 as a function of T 1 /ω 0 . Note that Eq. 17, and not any simplified expression, was used to calculate N and Ĥ S in the plots. The plots clearly show data collapse over the entire temperature regime for χ Ω 0 . Also, important to note, is the substantial effect of small interaction strengths at high temperatures. The high temperature scaling behaviour for small interaction strengths is same as that for large interaction strengths, but there is no data collapse. This is because, the conditionT χ Ω 2 0 means that, for small χ, there are large subleading terms. The low temperature behaviour matches with NESB. Since all these observations are valid for any choice of r, it fol- lows that they are all valid in equilibrium (i.e, r = 1) also. This is consistent with previous results in [37]. The crucial point in above discussion was to find that high temperature non-equilibrium results can be described via an effective temperatureT . In the following, we will see that even transport properties at high temperatures can be described in terms ofT .

F. Currents
Now we look at the average transport properties of the system, in particular we compute particle and energy currents. To calculate current we look at the evolution equations of the expectation values N and Ĥ S ofN and H S . SinceN andĤ S are diagonal in the eigenbasis of the system Hamiltonoan we can directly obtain the evolution of their expectation values from Eq. (10). This gives Collecting all terms depending on each bath separately, the above equations can be written like continuity equations of the forms d N /dt = I 1 −I 2 , d Ĥ S /dt = J 1 −J 2 , where I (J ) is the particle (energy) current flowing into the system from th bath. In steady state, I 1 = I 2 = I and J 1 = J 2 = J. The steady state expressions for cur-rents are : with I(ω n−1 ) = ε 2 Γ 1 Γ 2 J (ω n−1 )(n 1 (ω n−1 ) − n 2 (ω n−1 )) Γ 1 n 1 (ω n−1 ) + Γ 2 n 2 (ω n−1 ) .
The second steps of Eq. 29 have been arrived at from the first steps after some simplification using the property of NESS density matrix given in Eq. 12. Note that energy and particle currents are not independent. But, in general, there is no way of directly finding one current given the other and they can have quite different behaviour. Fig. 3 shows variation of energy and particle currents with interaction strength χ for Ohmic baths ((s = 1 in Eq. 5)), for both forward (∆T > 0) and backward (∆T < 0) biases. The mean temperature T m = (T 1 + T 2 )/2 is kept fixed in the plots, and systembath coupling is asymmetric (Γ 1 = Γ 2 ). The plots immediately show us a number of physical aspects of the system.
Firstly, we note that the particle current decreases with increase in interaction strength χ. This is expected because of increasing repulsive interaction in the system. On the other hand, energy current shows non-monotonic behaviour with χ. This is plausible because, while, with increasing χ, system allows less number of particles to pass, higher energy particles have a larger probability to pass through the system.
Secondly, we see that there is rectification of both energy and particle currents, since the particle and energy currents for forward and backward biases do not match. This is to be expected because the expressions for currents in Eq. 29 are not antisymmetric under interchange of hot and cold baths (i.e, n 1 ↔ n 2 ) in general. It is only so under special conditions. Two of such special conditions where there is no rectification are when χ = 0, i.e, when the system is linear, and when Γ 1 = Γ 2 (for any χ). These can be easily checked from the expressions for currents (Eq. 29). Hence, in general, there will be rectification effects in both particle and energy currents for χ = 0 and Γ 1 = Γ 2 . This is the generic behaviour in non-linear (interacting) systems.
Thirdly, as discussed before, for χ Ω 0 , T 1 , T 2 , the system behaves as NESB, and currents match with the NESB results. But, the rectification in the NESB limit is less than that for smaller interaction strengths. Thus rectification behaviour is non-monotonic as a function of χ. Our findings therefore suggest that a careful engineering of the system Hamiltonian is required to get maximum rectification from a given system. Finally, and most interestingly, for small interaction strength, the rectification of energy current occurs in the opposite direction to rectification of particle current. Also, there is a non-zero value of χ where the forward and backward energy currents match, and hence there is no rectification. At this point, the system rectifies particle current but not the energy current. Beyond this value of χ, energy and particle rectification occur in the same direction (see Fig. 3, bottom panel).
In what follows, we investigate the behaviour of particle and energy currents and rectification in more detail along with their scaling behaviour.

Scaling behaviour of currents
As we have seen with average system properties, transport properties also behave differently for different relative values of temperatures and interaction strength. In the NESB regime, χ Ω 0 , T 1 , T 2 , the currents are given by with ω 0 = Ω 0 + χ. This is identical to the expression for current previously derived for NESB [42]. Note that in NESB regime, energy current is proportional to particle current. This is because, in this limit, transport is allowed through transfer of exactly one particle through the system, and that particle has energy Ω 0 + χ. This is not valid beyond the NESB regime.
Now, let us look at the high temperature regime, T 1 , T 2 χ, Ω 0 , where NESB results are not valid. We note that the expression for the currents in Eq. 29 has the same form as Eq. 23. So, the high temperature trick in Eq. 24 can be readily applied to obtain with A = ε 2 Γ 1 Γ 2 /(Γ 1 + Γ 2 ) and ∆T = T 1 − T 2 . In the second step above, we have expanded the Bose distributions in I(ω n ) (Eq. 30) for high temperatures to obtain I(ω n ) ≈ AJ (ω n )∆T /T . Using the general form of spectral function given in Eq. 5, and after some algebra, we obtain For the choice of spectral function in Eq. 5, we can relate particle and energy currents via the function K(s). We now look at the properties of the function K(s). First, we look at the regimeT χ Ω 0 . In this regime, where F(z, s) = A √ π χ T ∞ 0 dy y − 1 2 e −y × √ zy + 1 1 + 2 √ zy s .

This then gives
We immediately make the following observations. Firstly, for T 1 , T 2 χ Ω 0 , K s /(χ s−1 ∆T ) varies as a function ofT /χ. Thus, I/ χ s−1 ∆T and J/ χ s ∆T for various values of χ should show a data collapse when plotted withT /χ. Secondly, we see K(1) ≈ A∆T . It follows that for Ohmic baths (s = 1 in Eq. 5) at high temperatures, particle current is independent of the interaction strength χ and one gets a linear-response-like relation I ≈ A∆T , even for large temperature bias. This is consistent with high temperature result for a harmonic oscillator (χ = 0).
However, the energy current J ∼ T /χ∆T and shows the effect of interaction. On the other hand for constant bath (s = 0), the energy current J always satisfies the linear-response-like relation, whereas the particle current is suppressed by a factor χ/T . Thirdly, we note that the quantity J/(χI) varies as a function ofT /χ and scales as T /χ for any s, In the other regime, T 1 , T 2 Ω 0 χ, the scaling behaviour is the same, but the data collapse is difficult to see because of large sub-leading terms.
Note that, I/∆T and J/∆T actually give the beyond linear response analog of particle and energy conductance. For linear response,T = T in the RHS of Eq. 34, T being the equilibrium temperature. Then, the above discussion gives the high temperature scaling of conductance.
In Fig. 4, we show the scaling behaviour of I/∆T and J/∆T as a function of T 1 /ω 0 . As in Sec. ( E), we have chosen T 1 /ω 0 as the relevant scaling variable and have kept r fixed at r < 1 to ensure beyond linear response regime. r → 1 gives the linear response conductance. The plots (a) and (b) of Fig. 4 are then the beyondlinear-response equivalent of temperature scaling of conductance. Note that Eq. 29 (and not any simplified expression) was used to calculate currents. The plots show the scaling behaviour discussed above. The plots for χ Ω 0 show data collapse over the entire range of temperatures. Also, the low temperature behaviour is given by NESB.
In Fig. 4(a), we plot I/∆T for fixed r for Ohmic baths (s = 1 in Eq. 5). We see that, for NESB, this quantity behaves non-monotonically with temperature, while for our non-equilibrium SSBH model, this quantity monotonically increases with temperature. In fact, the deviation from NESB result starts precisely at the point where the NESB result reaches a maximum. Similar behaviour is observed for J/(∆T ω 0 ) in Fig. 4(b). We conclude that, the non-monotonic behaviour of conductance of NESB, both at and beyond linear response, comes as a result of truncation of the energy spectrum and is not observed when all energy levels are considered. This demonstrates a major difference between the high temperature behaviour of non-equilibrium SSBH model and the NESB model.
Also, it follows from Fig. 4(a), that, for Ohmic baths (s = 1 in Eq. 5), at high temperatures the particle current (for fixed ∆T ) becomes independent of both interaction strength and the effective temperature. Fig. 4(b) shows a considerable deviation from linear behaviour even for small interaction strengths. Fig. 4(c) shows that, consistent with our previous discussion, irrespective of the choice of bath spectral function, the quantity J/(Iω 0 ), for χ Ω 0 , shows a data collapse and goes as ∼ T 1 /ω 0 for T 1 ω 0 (with r kept fixed).
Having discussed the scaling behaviour of currents in detail, we now look into another important property of interacting systems, the rectification of current.

Rectification
Rectification of current is a generic behaviour of nonlinear (interacting) systems in non-equilibrium. As we have seen, in a non-equilibrium set-up, two kinds of currents through the system can be defined, the particle current and the energy current, and their rectification behaviour can also be quite different. To our knowledge, there has been no previous work where both particle and energy current rectification for a bosonic non-linear system has been investigated. Also, note that rectification can only be observed beyond linear response regime.
Since rectification occurs only for asymmetric systembath coupling, we use the following definition to describe the degree of asymmetry where 0 ≤ γ ≤ 1 is dimensionless. Given a value of asymmetry parameter γ, we define a measure of rectification as R I and R J are the particle and energy current rectifications. This measure of rectification is as used in [42,43]. Note that, by this definition, rectification is positive if higher current flows when the cold bath is more strongly coupled to the system. Also, R I and R J are zero when γ = 0, 1. In our following discussion of rectification, we will primarily confine ourselves to the Ohmic baths (s = 1 in Eq. 5). The variation of R I and R J with γ is shown in Fig. 5. Irrespective of the strength of interaction strength, we notice that the maximum rectification occurs when γ ≈ 0.6. The figure also shows, both particle and energy rectifications behave non-monotonically with χ. For small χ, R J is negative while R I is positive, hence, direction of energy rectification is opposite to particle rectification. For large χ, the rectification is same as that obtained from NESB. In the NESB regime, particle and energy rectifications are same, because, particle current is proportional to energy current. All these observations are consistent with our discussion of Fig. 3.
To concisely investigate the rectification behaviour of the system as a function of the interaction strength and the temperatures, we again resort to the scaling variable T 1 /ω 0 with r fixed at r < 1. The plots of R I and R J as a function of T 1 /ω 0 for fixed r for Ohmic baths are shown in Fig. 6. We readily make the following observations : Firstly, for χ < Ω 0 , rectification is small. For χ Ω 0 , there is data collapse as expected from the scaling of currents. Secondly, where NESB matches the non-equilibrium SSBH model, there is small rectification. This can be understood from expression for currents in NESB regime given in Eq. 31. NESB result holds when χ Ω 0 , T 1 , T 2 . Therefore, in this regime ω 0 T 1 , T 2 . So the bose distributions in Eq. 31 are exponentially small. Hence, 1 + n(ω 0 ) ≈ 1. With this approximation, the expressions of currents in Eq. 31 become antisymmetric under interchange of hot and cold baths. Thus the NESB regime of the SSBH model gives very small rectification. Maximum particle rectification R I is reached when ω 0 < T 1 which is outside this regime.
Thirdly, after the maximum, R I gradually approaches zero with increase in T 1 /ω 0 . This is expected because, as we have seen before, at high temperatures, for Ohmic baths, the particle current eventually has the form, I ≈ A∆T , which is antisymmetric under interchange of hot and cold baths. On the other hand, the corresponding NESB rectification (the dotted line in Fig. 6) continues to increase with T 1 /ω 0 , until it saturates to a high value. Thus, as an effect of having all energy levels, and not truncating at two levels, the particle rectification is suppressed at high temperatures.
Fourthly, and most interestingly, energy rectification R J also has a peak for ω 0 < T 1 . However, R J = 0 at ω 0 ≈ T 1 . At this point, the particle rectification is still positive. So the system behaves as a particle recti-fier and not as an energy or heat rectifier. Beyond that point, R J changes sign. Thus the direction of rectification is reversed. With further increase in T 1 /ω 0 , energy rectification continues to grow in the reversed direction. Therefore, at high temperatures, the heat rectification occurs in the opposite direction to particle rectification, and continues to grow as temperatures of both hot and cold baths are increased.
Note that all the above observations are for Ohmic baths (s = 1 in Eq. 5). Finally, we discuss the case of non-Ohmic baths. For superohmic baths (i.e, for s > 1), both particle and energy currents show a reversal of direction of rectification. The reversal of energy rectification for superohmic bath occurs at a lower value of T 1 /ω 0 than that for the Ohmic bath. For subohmic baths (i.e, for s < 1), the reversal of energy rectification occurs at higher values of T 1 /ω 0 than that for the Ohmic bath. The particle rectification does not show reversal for subohmic baths. For the constant bath, s = 0, neither particle nor energy rectification shows reversal. However, in all cases, both particle and energy rectifications vary nonmonotonically with the interaction strength.

G. Time dynamics
Until now, we have discussed the properties of the NESS of the out-of-equilibrium SSBH model. In this section, we look at the transient time dynamics of the various physical quantities we have so far calculated in NESS.
To do this we revert to the equation for time evolution of ρ n , Eq. 10. The equation can be re-written and solved in the form whereρ is an infinite dimensional column vector containing diagonal elements of the density matrix and M is a infinite dimensional non-Hermitian square matrix containing the entries of Eq. 10. M has the form where C n and D n are as defined in Eq. 11. Note that the matrix M has the form of a Markov matrix. The sum of each column is zero (D 0 = 0 by definition). This corresponds to the fact that the trace of the density matrix is preserved, i.e, n ρ n = 1. To calculate the time dynamics, we choose an initial state with no particles in the system, i.e, initially, ρ 0 (0) = 1, and ρ n (0) = 0 ∀n = 0. The Eq. 39 is used to numerically obtain the time evolution. Even though the matrices involved are infinite dimensional, for given interaction χ and temperatures T 1 and T 2 , only a finite number of levels, determined by the ratio of the temperatures and the interaction, effectively contribute. Thus, starting from a finite matrix size, a convergence is reached as the matrix size is increased. Smaller interaction and higher temperatures require larger matrix sizes. A subtle point to note is that, if the matrix M is truncated at any finite size, say p, then the constraint that the sum of each column should be zero is not satisfied for the pth column, unless C p = 0. Consequently, the matrix M can be truncated at size p only if C p D p . Since we are using Redfield equation under Born-Markov approximation to obtain transient time dynamics, we need to be careful in choosing the observation times. This is because, as mentioned in Sec. C, Markov approximation is valid only when observation times are much larger than the time for decay of bath correlation functions. We have made some estimates (similar to that in Ref. [50]), which indicate that the bath relaxation times are indeed much smaller than the observational transient times of our interest.
Time dynamics of physical quantities N (t) , H s (t) , I(t), J(t) for Ohmic baths (s = 1 in Eq. 5) are shown in Fig. 7, where I 1 (t) and J 1 (t) are respectively particle and energy currents from the left bath into the system. (Before reaching steady state, the currents from left and right baths are not the same.) Similar time dynamics, but for bosonic system of two sites without interactions, was calculated in Ref. [50]. We observe that, unlike the two site non-interacting case in Ref. [50], here, the currents show no oscillations. More interestingly, we also observe that the time to reach steady state, called t ss hereafter, decreases with increase in the interaction. In the following, we investigate the dependence of t ss on set-up parameters more carefully.
We note from the solution ofρ(t) in Eq. 39 that time can be scaled as ε 2 t. It follows that t ss ∝ ε −2 . Thus, time to reach steady state increases as system-bath coupling becomes weaker.
It is also clear from Eq. 39 that the steady state is given by the eigenvector of the matrix M corresponding to zero eigenvalue. The fact that a unique steady state is reached in long time then implies that all other eigenvalues of the M have positive real part. Interestingly, we have found in our numerical computation that the eigenvalues of M are all real (and hence no oscillations in time). The smallest non-zero eigenvalue then gives a measure of t ss . So, we define λ 1 being the smallest non-zero eigenvalue of M . Even though t ss cannot be calculated for all set-up parameters analytically, in two limiting cases, analytical results can be obtained. The first case corresponds to the NESB regime χ Ω 0 , T 1 , T 2 . In this case, C 1 D 1 , only two levels effectively contribute and matrix M has the form The non-zero eigenvalue of M is The time to reach steady state therefore also depends on the temperatures of the baths. However, for large ω 0 = Ω 0 + χ, the Bose distributions are expotentially small. Thus, for χ Ω 0 , T 1 , T 2 and for general bath of the form given in Eq. 5, which is independent of the temperatures of the baths.
We also see that, for constant baths (s = 0 in Eq. 5), t ss is independent of interaction strength. For other baths, t ss decreases with increase of interaction strength as a power law. The second case where t ss can be analytically, calculated corresponds to the linear system, χ = 0. In this case, C n does not decay with n (Eq. 11 for χ = 0) and hence the matrix M cannot be truncated at any finite size. So, the above method of finding t ss fails. However, since, in this case, we have a non-interacting system, we can find t ss directly from the evolution equation for N (t) . The evolution equation for N (t) can be obtained from Eq. 27 by setting χ = 0. The resulting equation can be written and solved in the form with N ss = 2 =1 J (Ω 0 )n (Ω 0 ) / 2 =1 J (Ω 0 ) being the NESS occupation. From above equation, it is clear that for χ = 0, the t ss is given by This is again independent of the temperatures of the baths. Except for these two limiting cases, t ss needs to be found numerically using the definition Eq. 41. Fig. 8 shows log-log plot of numerically obtained t ss as a function of interaction strength χ for Ohmic baths for two different choices of temperatures of hot and cold baths. It is seen that except for the limiting cases, t ss depends on the temperatures of the baths. The limiting cases show the behaviour discussed above.

H. Conclusion
We have investigated a system consisting of a single bosonic site with Bose-Hubbard interaction weakly coupled to two bosonic baths at different temperatures. We have used the Redfield QME method to obtain analytical results beyond linear response regime. Below, we summarize our main findings and their potential applications.
We have found an analytical result for the population density (diagonal elements of the density matrix in the eigenbasis of system Hamiltonian) of the system in NESS (Eq. 13). This has further allowed us to find various physical observables like average occupation and energy, as well as the particle and energy currents in NESS. We have then analytically found interesting scaling behaviour of the physical observables. Our main finding in this respect is that the high temperature behaviour of the system can be described in terms of an effective temperature (Eq. 20). Then, it follows that, with the ratio r of the temperatures of cold and hot baths fixed, there occurs a data collapse for various strengths of interaction strength χ, when various physical observables of the system are plotted in terms of the scaling variable T 1 /ω 0 (Figs. 2,4,6). The scaling behaviours hold for a general choice of bath spectral functions of the form given in Eq. 5. We have also found very interesting rectification behaviour of the system. The most interesting finding is that, for Ohmic and subohmic baths (0 < s ≤ 1 in Eq. 5), the energy current shows a reversal in direction of rectification (Fig. 6). It follows that there is a nonzero strength of interaction strength, (χ ≈ T 1 − Ω 0 , for Ohmic baths), where energy or heat rectification is zero.
At this point, the system behaves as a particle rectifier but not as a heat rectifier. For superohmic baths (s > 1 in Eq. 5), both particle and energy currents show reversal in direction of rectification. For constant baths (s = 0 in Eq. 5), there is no change in direction of rectification for both particle and energy currents. Therefore careful engineering of baths can lead to various interesting rectification behaviour of the system. Such phenomena can be potentially used to create quantum devices, such as, optical diodes. Reversal of direction of thermal rectification of a quantum system has also been previously theoretically seen for a Heisenberg spin chain out of equilibrium [51]. Further, we have computed non-unitary timedynamics of all the above mentioned physical quantities. We found that, except for constant baths, the time to reach steady state is shorter for higher interactions and higher system-bath couplings. For constant baths, with increase in interaction, the time to reach steady state approaches a constant independent of the strength of interaction. For large χ, the time to reach steady state goes as t ss ∼ ε 2 χ s −1 for a general choice of bath spectral functions of the form given in Eq. 5. All our results also are consistent with the linear, χ = 0 (harmonic oscillator) case, as well as the NESB, χ Ω 0 , T 1 , T 2 , case.
Our results are experimentally relevant in quantum hybrid systems, where a single site Bose-Hubbard model can be realized, as well as, in molecular junction systems, where our set-up describes a model for anharmonic junctions. Future work includes going to strong systembath coupling regime, as well as to generalize to nonequilibrium interacting systems of more than one site, such as, the non-equlibrium Bose Hubbard chain [52,53] and Jaynes-Cummings Hubbard model [54].