Current-induced forces in single-resonance systems.

In this paper, we will consider the tight-binding system depicted in Fig. 1. It consists of a single energy level with energy $E_{d}$ coupled, with hopping constants $V_{L}$ and $V_{R}$, to two semi-infinite tight-binding chains, with site energy $E_{0}$ and hopping $V_{0}$. This Hamiltonian represents a minimum model for a quantum dot coupled to two reservoirs and is widely used in the context of quantum transport. The total Hamiltonian reads

where $E^{(0)}=E_{d}$ and $E^{(n)}=E_{0}$ for $n\neq 0$, and the hopping constants $V^{(-1)}=V_{L}$, $V^{(0)}=V_{R}$, and $V^{(n)}=V_{0}$ for $n\neq\left\{-1,0\right\}$.

For convenience, we will define $\hat{H}_{S}$ as the Hamiltonian of the system which includes the first sites of each semi-infinite chain representing the leads. Then, the three by three matrix $\boldsymbol{H}_{S}$ contains the site energies $E^{(-1)}=E_{0}$, $E^{(0)}=E_{d}$, and $E^{(1)}=E_{0}$, as well the couplings $V_{L}$ and $V_{R}$.

The matrix representation of the eigenvalue equation $\hat{H}\left|\psi\right\rangle=\varepsilon\left|\psi\right\rangle$, which in principle is of infinite dimension, can be reduced to an effective system of finite dimension by using the decimation technique Pastawski and Medina Dagger (2001); Cattena et al. (2014). There, the sites connected to leads are corrected by a self-energy associated to the $\boldsymbol{G}^{r}(\varepsilon)$

The real ($\Delta$) and imaginary ($\Gamma$) parts of $\Sigma$ are

where $\varepsilon-E_{0}\geq 2V_{0}$ for case A, $\varepsilon-E_{0}\in[-2V_{0},2V_{0}]$ for case B, and $\varepsilon-E_{0}\leq-2V_{0}$ for case C. The region where $\Gamma\neq 0$ (case B), the band, corresponds to the energy regions where excitations propagate freely inside the leads. Using the above self-energy we can define the effective Hamiltonian $\hat{H}^{r}_{\rm eff}\left(\varepsilon\right)=\hat{H}_{S}+\hat{\Sigma}^{r}\left(\varepsilon\right)$, where $\hat{\Sigma}^{r}\left(\varepsilon\right)$ is given by $\hat{\Sigma}^{r}\left(\varepsilon\right)=\Sigma_{0}^{r}\left(\varepsilon\right)\left(\left|-1\right\rangle\left\langle-1\right|+\left|1\right\rangle\left\langle 1\right|\right)$. In matrix form $\hat{H}^{r}_{\rm eff}$ is

By means of Eq. 1, $\boldsymbol{H}^{r}_{\rm eff}(\varepsilon)$ can be used to calculate $\boldsymbol{G}^{r}$, while $\boldsymbol{G}^{a}=\left[\boldsymbol{G}^{r}\right]^{\dagger}$. The resulting eigenvalue equation $\left(\varepsilon\boldsymbol{I}-\boldsymbol{H}^{r}_{\rm eff}\left(\varepsilon\right)\right)\left|\psi\right\rangle=0$ is now finite and exact. The price to be paid is that now the coefficients of $\boldsymbol{H}^{r}_{\rm eff}$ are energy-dependent. Notice that, because $\Sigma^{r}\in\mathbb{C}$, the effective Hamiltonian becomes non-Hermitian, and thus, its eigenvalues, or equivalently, the poles of the associated Green’s function $\boldsymbol{G}^{r}$ are not necessarily real. Importantly, these poles can change abruptly with the parameters of the Hamiltonian and, in consequence, also the LDOS [see panels from (c) to (e) in Fig. 2 ] and the dynamics of the system [see panels from (f) to (h) in Fig. 2 ].

Finally, one last detail. Due to technical reasons, we perturbatively couple a third lead to the system, see $\Gamma_{\eta}$ in Eq. 8. This reservoir is modeled by a self-energy in the wideband approximation $\Sigma_{\eta}=-i\Gamma_{\eta}$ with $\Gamma_{\eta}$ positive. The purpose of this extra reservoir is to add a small broadening to the LDOS and to guarantee an occupation to the system. Both effects are particularly important when the eigenenergies of the effective Hamiltonian lay outside the band or when $V_{L}=V_{R}=0$. Despite this, in all calculations, we take the limit $\Gamma_{\eta}\rightarrow 0$, and thus its effects on the poles structure is negligible.

The poles of the system correspond to the singular points of the Green’s function. The specific tight-binding model described in the previous section has two poles $\left\{\varepsilon_{p,+},\varepsilon_{p,-}\right\}$ given by the solutions of the secular equation $\det\left(\varepsilon\boldsymbol{I}-\boldsymbol{H}_{\rm eff}\right)=0$. Only three variables are really necessary to describe such poles, $E_{d}$, $V_{L}$, and $V_{R}$. The effect of $E_{0}$ is only to shift the poles’ position while $V_{0}$ gives the energy scale. For that reason, we define the dimensionless parameters

where $v_{L}=\frac{V_{L}}{V_{0}}$ and $v_{R}=\frac{V_{R}}{V_{0}}$. The real and imaginary part of the dimensionless poles are, respectively,Bustos-Marún et al. (2010); Garmon et al. (2012)

for $v^{2}\neq 1$. The first case of $\mathrm{Re}\left(\epsilon_{p,\pm}\right)$ and $\mathrm{Im}\left(\epsilon_{p,\pm}\right)$ corresponds to the condition $\left(\frac{\epsilon_{d}}{2}\right)^{2}+v^{2}\geq 1$, while the other one corresponds to the condition $\left(\frac{\epsilon_{d}}{2}\right)^{2}+v^{2}\leq 1$

In Fig. 2-(a) and (b) we show the real and imaginary part of the poles, $\epsilon_{p,\pm}$, as a function of the dimensionless parameter $\epsilon_{d}$. In the lower panels of the figure, see panels (c)-(h), we show the typical spectral and dynamical characteristic of three types of poles defining the resonant state, the virtual state, and the localized state.

Resonant states. They correspond to poles with a nonzero imaginary part. See the green dot marked as (1) in Fig. 2-(a). When the system presents such a pole, the LDOS is typically a Lorentzian peak centered at approximately $\mathrm{Re}(\epsilon_{p})$ and whose width is given by $\mathrm{Im}(\epsilon_{p})$. This description is particularly accurate when $\mathrm{Re}(\epsilon_{p})$ is close to the center of the band. However, when $\mathrm{Re}(\epsilon_{p})$ approaches one of the band edges the maximum of the LDOS is shifted toward the closest band edge and its shape deviated from a perfect Lorentzian function. After a short quadratic decay, the dynamics of such systems typically shows an exponential decay. Rufeil-Fiori and Pastawski (2006)

Localized states. They correspond to poles with $\mathrm{Im}(\epsilon_{p})=0$ and where the $\mathrm{Re}(\epsilon_{p})$ lays outside the band ($|\mathrm{Re}(\epsilon_{p})|>2$). For instance, see the green dot marked as (3) in Fig. 2-(b). When the system presents such a pole, its LDOS typically shows a very narrow peak outside the band (a Dirac’s delta in the limit $\Gamma_{\eta}\rightarrow 0$). The survival probability may show a small decay at short times but then it remains constant with time.

Virtual states. They also correspond to poles with $\mathrm{Im}(\epsilon_{p})=0$ and $|\mathrm{Re}(\epsilon_{p})|>2$, see the green dot marked as (2) in Fig. 2-(b). However, they are poles of the nonphysical Riemann sheet of $\Sigma^{r}$, see Refs. Dente et al., 2008 and Bustos-Marún et al., 2010. for more details. Due to this, there is not a peak in the LDOS at $\mathrm{Re}(\epsilon_{p})$. However, the effect of these poles on the LDOS within the band is the same as that of localized states. There is an accumulation of states near the closest band edge, which grows as the pole approaches this band edge. See the examples shown in panels (d) and (e).

The three type of poles described above determine two QDPT namely the resonant-virtual and the virtual-localized transitions. The Resonant-virtual QDPT is given by Bustos-Marún et al. (2010); Garmon et al. (2012)

while the Virtual-localized QDPT is given by Bustos-Marún et al. (2010); Garmon et al. (2012)

A more detailed classification may lead to another QDPT as resonant states can be subdivided into two additional type of poles, see Ref. Dente et al., 2008; Bustos-Marún et al., 2010. However, for the purpose of this work this subclassification is irrelevant and will be ignored.

One final remark is still in order. In more general scenarios of a quantum dot coupled to $m$ identical semi-infinite chains with coupling $V_{i}$, the corresponding poles are still given by Eq. 10, but the quadratic coupling $v^{2}$ should be replaced by

The mean value of the force operator can be calculated in terms of Green’s functions by

where $\boldsymbol{\Lambda}_{\nu}=-\frac{\partial\boldsymbol{H}}{\partial X_{\nu}}$ and $\boldsymbol{\mathcal{G}}^{<}$ is the lesser Green’s function. This function, which for shortness we omitted its energy and time dependence ($\boldsymbol{\mathcal{G}}^{<}\equiv\boldsymbol{\mathcal{G}}^{<}(\varepsilon,\tau)$), is indeed the Wigner transform of what generalizes the density matrix and is called the particle propagator $\mathcal{G}^{<}(x,t;x^{\prime},t^{\prime})=(i/\hbar)\left<\psi^{\dagger}(x^{\prime},t^{\prime})\psi(x,t)\right>$.Haug and Jauho (2008); Bode et al. (2011) The closely related Green function $\mathcal{G}^{>}(x,t;x^{\prime},t^{\prime})=-(i/\hbar)\left<\psi(x,t)\psi^{\dagger}(x^{\prime},t^{\prime})\right>$ is dubbed the hole propagator, since the order of the creation and annihilation operators are reversed. Both functions are directly linked to observables and kinetic properties, such as particle densities and currents. Pastawski (1992); Haug and Jauho (2008)

In complex scenarios like the one treated here, where the time-dependent Hamiltonians are slowly varying functions of some parameters, the function $\boldsymbol{\mathcal{G}}^{<}$ can be calculated by resorting to an adiabatic expansion in terms of the time-dependent parameters. Haug and Jauho (2008); Pastawski (1992) There, $\boldsymbol{\mathcal{G}}^{<}$ is given, up to first order, by, Bode et al. (2011)

where $\boldsymbol{G}^{<}$ is the frozen lesser Green’s function, given by

Here, $\boldsymbol{G}^{r/a}$ are given by Eq. 1 and $\boldsymbol{\Sigma}^{<}$ is the lesser self-energy, $\boldsymbol{\Sigma}^{<}=2i\underset{\alpha}{\sum}f_{\alpha}\boldsymbol{\Gamma}_{\alpha}$. The latter contains the information of the reservoirs, where $f_{\alpha}$ is the Fermi-Dirac distribution function of reservoir $\alpha$ with chemical potential $\mu_{\alpha}$ and temperature $T_{\alpha}$. The function $\Gamma_{\alpha}$ is the imaginary part of the self-energy of the reservoir $\alpha$, and describes the escape rate from the region of the contact in units of $\hbar$. Importantly, in deriving Eq. 17 it was assumed that the self-energies $\boldsymbol{\Sigma}^{<}$ and $\boldsymbol{\Sigma}^{r/a}$ are time-independent. However, as we included the first sites of the leads into the effective Hamiltonian, Eq. 8, a variation of the couplings $V_{L}$ and $V_{R}$ does not change the self-energies.

Then, under the adiabatic approximation, $\mathcal{F}_{\nu}$ is

where the “frozen” force or, from now on, simply the force, is

The terms $\gamma_{\nu\nu^{\prime}}$ are the components of the electronic-friction tensor $\boldsymbol{\gamma}$. The diagonal elements of $\boldsymbol{\gamma}$ provide the mechanical dissipation while the nondiagonal components are Lorentz-like terms which allow energy transfer between modes. This tensor can be linearly decomposed into a symmetric and an antisymmetric contribution $\boldsymbol{\gamma}=\boldsymbol{\gamma}^{s}+\boldsymbol{\gamma}^{a}$. The antisymmetric component is finite only for nonzero voltage biases or temperature gradients between the reservoirs. Thus, close to equilibrium conditions only the symmetric contribution of the electronic-friction tensor needs to be taken into account.

In this work, we are interested only on the leading orders of the expansion of $\mathcal{F}$ in terms of the different nonequilibrium sources $\delta\mu$, $\delta T$, and $\dot{\overrightarrow{X}}$,

The equilibrium contribution of the force,

is conservative and thus can be written as the gradient of a potential, $\overrightarrow{F}^{\mathrm{eq}}=-\nabla U^{\mathrm{eq}}$. As we show in appendix A, the equilibrium potential $U^{\mathrm{eq}}$ is

where $f_{0}$ is the equilibrium Fermi-Dirac distribution function given by $\mu_{0}$ and $T_{0}$, the equilibrium chemical potential and temperature respectively.

Comparing Eqs. 21 and 19, it is obvious that $\left.\frac{\partial\mathcal{F}_{\nu}}{\partial\dot{X}_{\nu^{\prime}}}\right|_{\mathrm{eq}}=\gamma_{\nu,\nu^{\prime}}|_{\mathrm{eq}}$. Therefore, only the equilibrium contribution of the symmetric component of $\boldsymbol{\gamma}$ is relevant for the present work, where $\gamma_{\nu\nu^{\prime}}^{s,\mathrm{eq}}\equiv\gamma_{\nu\nu^{\prime}}^{\mathrm{eq}}$ and

Here, $\boldsymbol{G}^{\lessgtr}_{\mathrm{eq}}$ are the equilibrium greater ($>$) and lesser ($<$) Green’s functions. The above equation can be rewritten, see appendix B, in a simpler form,

where $\boldsymbol{A}$ is the spectral function, $\boldsymbol{A}=i\left(\boldsymbol{G}^{r}-\boldsymbol{G}^{a}\right)$. This equation shows that only states with energies close to the Fermi energy contribute to the equilibrium electronic-friction tensor. Those states are also responsible for the electric current Pastawski and Medina Dagger (2001) and the nonequilibrium contribution of the CIFs, as we will see afterward. Therefore, Eq. 25 highlights a fundamental physical connection between these quantities. In the low temperature limit, Eq. 25 becomes

Stochastic forces, in Eq. 14, can be obtained by using the above equations and the fluctuation-dissipation relation, Bode et al. (2011) given by

where $k_{B}$ is the Boltzmann constant and it is assumed that stochastic forces are described by $\left\langle\xi_{\nu}\left(t\right)\xi_{\nu^{\prime}}\left(t^{\prime}\right)\right\rangle=D_{\nu\nu^{\prime}}^{\mathrm{eq}}\delta\left(t-t^{\prime}\right)$.

At low voltage biases and temperature gradients, the nonequilibrium forces are given by

which in the low-temperature limit (see appendix C) gives

Here, $\delta\mu_{\alpha}=\mu_{\alpha}-\mu_{0}$, $\delta T_{\alpha}=T_{\alpha}-T_{0}$, and the expression is only valid for $\delta T_{\alpha}\leq T_{0}$.

One of the main motivations of this work is to study the general conditions that maximize the performance of different nanoelectromechanical devices and quantum machines. In this respect, one of the central quantities to evaluate is the work $W$ done by CIFs during cyclic motions. If one is dealing with systems where only two or three parameters are required to describe their motion, then it is useful to study the curl of the CIFs, defined as

where $\rho$, $\nu$, and $\nu^{\prime}$ are the indexes of the coordinates and $\nu<\nu^{\prime}$. Now, the work can be obtained from

where the integration is done on the surface $S$ enclosed by the closed curve $C$ that describes the cyclic motion.

Now inserting Eq. 29 into Eq. 30, gives

where we used $\boldsymbol{G}^{r/a}\boldsymbol{\Lambda}_{\nu}\boldsymbol{G}^{r/a}=\partial_{\nu}\left(\boldsymbol{G}^{r/a}\right)$, $\partial_{\nu^{\prime}}\boldsymbol{\Lambda}_{\nu}=\partial_{\nu}\boldsymbol{\Lambda}_{\nu^{\prime}}$, and $\left\{\boldsymbol{\Lambda}_{\nu}..\boldsymbol{\Lambda}_{\nu^{\prime}}..\right\}_{a}=\left(\boldsymbol{\Lambda}_{\nu}..\boldsymbol{\Lambda}_{\nu^{\prime}}..-\boldsymbol{\Lambda}_{\nu^{\prime}}..\boldsymbol{\Lambda}_{\nu}..\right)$. Eq. 32 gives the work done per unit area in the low temperatures limit and for small bias voltages or temperature gradients.

Let us return to the tight-binding model described in Sec. II.1.1. For the analysis of the forces, especially for the equilibrium forces, it is crucial to have a definite occupation of the system, even for localized states or when $V_{L}=V_{R}=0$. As discussed in section II.1.2, when the system is at the localized state, the poles of $\boldsymbol{G}^{r}$, or the eigenenergies of $\boldsymbol{H}_{\mathrm{eff}}$, lay outside the band. At this condition, the imaginary part of the self-energies, $\Gamma$, is zero. The same happens for $V_{L}=V_{R}=0$. This implies the nonphysical condition that localized states do not contribute to the forces, or even worse, that isolated systems do not have forces at all. To solve this problem, we added the third reservoir ($\eta$) to the model, assuming its occupation is given by the equilibrium Fermi function $f_{0}=(f_{L}+f_{R})/2$, and taking the limit $\Gamma_{\eta}\rightarrow 0$ in all calculations.

Given the Hamiltonian of Eq. 8 and the definition of $\boldsymbol{G}^{r/a}$ (Eq. 1) and $\boldsymbol{G}^{<}$ (Eq. 18), we obtain a closed-form expression for the Green functions, (see appendix D). As can be seen from the expressions of the appendix, the third reservoir does not affect the calculations for energies within the band, provided $\Gamma_{\eta}$ is small. However, it allows the calculation of equilibrium $\boldsymbol{G}^{<}$ (and thus equilibrium forces) at all conditions. Then, all occupied states, independently if they are inside or outside the band, contribute to the equilibrium forces. One important point to emphasize is that the nonequilibrium contribution of $\boldsymbol{G}^{<}$ (and thus $\overrightarrow{F}^{(\mathrm{ne})}$) is always zero for states outside the band (localized states), or for $V_{L}=V_{R}=0$. This will have important consequences when studying the maximization of the curl of the forces.

According to Eq. 16, the forces depend on how the elements of the electronic Hamiltonian are affected by the movement of the mechanical degrees of freedom, but this is system-dependent. Therefore, to gain generality, we take as “mechanical” variables the same parameters of the Hamiltonian ($E_{d},V_{L},V_{R}$), which leads to dimensionless forces $F_{\nu}$, see appendix E. One can readily calculate the forces in terms of “physical” variables $q_{i}$ by

where $X_{\nu}=\left\{E_{d},V_{L},V_{R}\right\}$. Notice that while forces may depend on the chosen variables, the work in Eq. 31 is independent of them, as long as the curve $C$ remains the same, see appendix E. Therefore, using our atypical choice of “mechanical” variables in Eq. 32 provides a useful tool to analyze in a general way the performance of nanoelectromechanical devices or quantum machines without having to specify the details of their physical implementation.