Éric Larose, ISTerre, Grenoble, France.

# Ultrasonic monitoring of stress and cracks of the 1/3 scale mock-up of
nuclear reactor concrete containment structure

Qi Xue11affiliationmark:  Eric Larose11affiliationmark:  Ludovic
Moreau11affiliationmark:  Romain Thery1,21,2affiliationmark:  Odile
Abraham22affiliationmark:  Jean-Marie Henault33affiliationmark:
11affiliationmark:  Univ. Grenoble Alpes, CNRS, ISTerre, CS 40700, 38058
Grenoble, France
22affiliationmark:  Uni. Eiffel, GERS-GeoEND, allée des ponts et chaussées -
CS5004, 44344 Bouguenais Cedex, FRANCE
33affiliationmark:  EDF R&D, Chatou, France  eric.larose@univ-grenoble-
alpes.fr

###### Abstract

To evaluate the stress level and damage of a reinforced concrete containment
wall (similar to those used in nuclear power plants) and its reaction to
pressure variations, we implemented successive ultrasonic experiments on the
exterior surface of the containment wall in the gusset area for three
consecutive years (2017, 2018 and 2019). During each experiment, the pressure
inside the containment wall increased gradually from 0 MPa to 0.43 MPa and
then decreased back to 0 Mpa. From the analysis of the ultrasonic coda waves
obtained in the multiple scattering regime (80-220 kHz), we performed Coda
Wave Interferometry to calculate the apparent velocity changes in the
structure (denoted by $dV/V_{a}$) and Coda Wave Decorrelation (DC)
measurements to produce 3D cartographies of stress and crack distribution.
From three source-receiver pairs, located at the top, middle and bottom of the
experimental region, we observe that coda waves dilate, shrink and remain
almost unchanged, respectively. This corresponds to the decreasing, increasing
and invariant pressure inside the concrete. The comparison of three years’
results demonstrates that the variation of $dV/V_{a}$ and DC under the same
pressure test increases through the years, which indicates the progressive
deterioration and aging of the concrete. From a large collection of source-
receiver pairs at different times, the spatial-temporal variations of
$dV/V_{a}$ and DC are then used to produce a map of the structural velocity
and scattering changes, respectively. We observe a decreasing velocity on the
top part and an increasing in the middle one, which is in line with the
$dV/V_{a}$ analysis. The reconstructed scattering changes (or structural
changes) highlight the active region during the inflation-deflation procedure,
corresponding to the opening and closing (and sometimes the development) of
cracks. The larger magnitude in 2019 than in 2017 indicates the increasing
damage in the concrete.

###### keywords:

Coda wave interferometry, diffuse ultrasound, cracks in concrete, velocity
change, stretching and correlation

##  1 Introduction

Concrete containments are generally used to protect the nuclear reactor of a
nuclear power plant in case of any leakage of the nuclear fuel or waste, or
from possible external aggression. The humidity, temperature and/or pressure
changes induce the deformation of the concrete shell and, together with
additional chemical processes, progressively alter the mechanical integrity of
the structure. Small cracks (millimeter scale) are generated and progressively
develop during the lifetime of the structure. The assessment of the state of
damage and the characterization of cracks are vital issues for the security of
the nuclear power plant in case of emergency, and for the estimation of the
residual life time. In France, a leak test of each concrete containment is
enforced by the law every ten
years111https://www.legifrance.gouv.fr/jorf/id/JORFTEXT000000469544, during
which the pressure inside the containment is increased from 0 MPa to about
0.43 MPa and then decreased back to 0 MPa.

In order to study the aging of the concrete wall and follow the evolution of
the leak rate owning to drying creep and loss of pre-stress in a realistic
environment, Électricité de France (EDF) launched a research program on a 1/3
scale mock-up, called VeRCoRs masson2013objectives ; mathieu2018temperature ;
henault2018characterize . This mock-up is a double walled concrete containment
building, as a 1:3 scale specimen representative for 24 (about 40 %) of the
French reactors. The present work explores the applicability and performance
of a method based on diffuse ultrasound analysis to image and quantify stress
changes and cracks within the structure at VeRCoRs during a standard internal
pressure test.

To this end, we take advantage of the high sensitivity of multiply scattered
waves at the ultrasonic scale, for frequencies of few hundreds of kHz. The
high heterogeneity of the concrete results in the loss of ballistic and
coherent properties of the waveforms after a propagation over a few
wavelengths. The late-arriving ultrasonic signals, also named coda waves, are
known to be highly repeatable and sensitive to small perturbations of the
medium. At the seismic scale, the sensitivity of the coda wave has been
extensively exploited to detect small changes in the Earth’s crust, e.g.,
fault zones wegler2007fault ; brenguier2008postseismic  or volcanic areas
gret2005monitoring ; sens2006passive ; brenguier2008towards .

Coda wave interferometry (CWI), as introduced in poupinet1984monitoring  and
further developed in snieder2002coda ; snieder2006theory , is a technique to
detect slight velocity perturbations in a highly heterogeneous medium, by
comparing coda waves recorded before and after a perturbation. A quantitative
way to compare two waveforms is to compute their correlation coefficient (CC).
One way of applying CWI consists in stretching or shrinking the waveform to
match two coda signals, i.e., to maximize the CC in a given time window, with
a dilation coefficient. This coefficient, associated to the local velocity
change lobkis2003coda , is generally referred to as the apparent relative
velocity change $dV/V_{a}$. Alternatively, it can be useful to calculate the
decorrelation coefficient (DC) instead of the CC, because the DC can be linked
to small structural changes planes2015imaging . The $dV/V_{a}$ and DC
variables are now widely used to detect ultrasonic velocity changes in a
concrete structure induced by thermal or stress variations
larose2006observation ; larose2009monitoring ; schurr2011damage , and/or to
detect micrometer scale cracks hilloulin2014small . The reader can refer to
planes2013review  for a review.

The relation between the $dV/V_{a}$ and the local velocity change is proposed
in pacheco2005time  by comparing wave propagation with random walk. The
relation between DC and the local scattering change (for example, induced by
cracks) follows a similar formula proposed later larose2010locating . The
reader may refer to planes2014decorrelation  for a detailed explanation and
proof. These two formulas have the same form: an integral equation with the
same sensitivity kernel. The sensitivity kernel is a function of the source,
receiver, space and time, that statistically estimates the density of time
spent by the diffuse wave at a given region of the medium. It can be estimated
either from i) the explicit solution of diffusive or radiative-transfer
approximation of ultrasonic waves for simple geometries (e.g., infinite plate
or cuboid) planes2014decorrelation , ii) numerical methods (e.g., finite
difference or finite element) for any shape xue2019locating  or iii) directly
solving the wave equation kanu2015numerical . By combining information from
several sources and receivers, we are able to obtain thousands of $dV/V_{a}$
and DC values between different source-receiver pairs in different times,
which can be used to image the velocity and scattering changes through solving
the integral equation. Please refer to planes2015imaging  for a numerical
simulation, zhang2016diffuse  and zhan2020three  for laboratory experiments
and zhang2018three  for a field experiment.

We implemented similar diffuse ultrasonic experiments on the outer surface of
the gusset area of VeRCoRs (see Fig. 1(a)) during the internal pressure test
over three consecutive years, in order to detect cracks and the effect of
damage on stress changes during each test. From the analysis of the $dV/V_{a}$
and DC for three different source-receiver pairs between the initial state and
other states during the pressure test, we aim at observing that different
parts of the wall are expected to have different responses to the same
inflation-deflation procedure. We also aim at quantifying the evolution of
$dV/V_{a}$ and DC versus aging and damage over the years. An imaging inversion
procedure is then applied to obtain 3D cartographies of relative velocity
changes and density of probability of cracks.

The comparison of different years’ results may quantify the deterioration of
the concrete: the sensitivity of $dV/V_{a}$ to internal stress is expected to
increase with time (thus with damage) as observed elsewhere planes2013review
and the decorrelation is expected to show some non-reversible (thus
cumulative) behavior.

##  2 Experimental setup

(a)

(b)

(c)

Figure 1: Configuration of the experiment. (a) Sketch of the concrete
structure and the experimental region (green). We focus on the lower corner
(the gusset) of the exterior surface of VeRCoRs, a 1/3 replica of a nuclear
containment wall. (b) source and receiver positions in 2017 and 2018. (c)
source and receiver positions in 2019 with new transducer positions in green.
The horizontal solid lines represent the change of thickness of the wall. The
thickness of the upper part is 40 cm and the middle part (gusset) 60 cm. The
bottom layer is the base slab. The dashed horizontal lines are positions of
the construction joint. Dotted lines are positions of cracks on the surface of
the wall observed in 2018.

VeRCoRs is a 1/3 scale experimental mock-up of a nuclear reactor containment
(see Fig. 1(a)) built by EDF near Paris222https://fr.xing-events.com/EDF-
vercors-project.html?page=1426249. Construction of the mock-up began in 2013
and was completed in March 2016. The construction procedure is layer by layer.
The horizontal dashed lines in Fig. 1(b) and (c) indicate positions of the
joint in our experimental region. After the construction, two auxiliary air-
conditioning units started to impose relative realistic humidity and
temperature inside and outside the containment to simulate regular operations
of a nuclear reactor. The 1/3 scale was chosen to accelerate aging, ten times
faster than in a realistic environment hopefully.

To monitor the properties of the concrete, the pressure inside the structure
increases gradually from 0 MPa to 0.43 MPa, which is followed by a static
period, during which the pressure decreases a little bit due to small leakage
from cracks, and then the whole structure is deflated to its initial state
(see Fig. 2(a)). This test was repeated several times, twice in 2017, and once
in 2018 and 2019, which corresponds to the standard leak test for real nuclear
power plant containment every ten years. We foresee a large change of the
concrete in the first test in 2017, because it was the first pressure test and
a lot of micro and macro-cracks were created. Therefore, in this paper we only
present the second test in 2017, to take into account the so-called Kaiser
effect that has been observed and commented previously in concrete
zhang2012study . Modifications observed in the following are thus not related
to the fact that the pressure level exceeds the first one. Our experiment was
implemented on the outer surface of the gusset region (a surface of about
2m$\times$2m). Fig. 1(a) provides a sketch of the structure and the
experimental region.

Several transducers are glued directly to the concrete to generate and receive
ultrasonic signals. Fig. 1(b) shows the positions of sources and receivers in
2017 and 2018, and Fig. 1(c) shows positions in 2019. Since the base of the
concrete structure (base slab) is extremely thick, the signals recorded by the
bottom source-receiver pair (S12-R12) remain almost unchanged during the
experiment in 2017 and 2018. We decided to remove them and to add several new
transducers to improve the resolution of the upper part in 2019. In addition,
we did not receive reasonable signals from R1 in 2018 which may be caused by
the failure or detachment of the transducer.

During the experiment, a chirp signal with frequency ranging from 80 to 220
kHz is generated, amplified and sent into the concrete from each source one
after another, and then recorded by all the receivers in order to reconstruct
all the source-receiver impulse responses, through deconvolution of recorded
waveforms with the chirp source signal. The whole acquisition process takes
about 3 minutes, within which the change of concrete can be neglected. The
same process is repeated every hour to obtain measurements at different
loading states.

##  3 Apparent velocity change and decorrelation analysis

Figure 2: (a) pressure variation during the inflation-deflation procedure. We
took the measurement every an hour, during which the pressure has a change of
+ or - 0.02 MPa except in steady states. (b) DC between the original state and
other states; (c) $dV/V_{a}$ variation with respect to pressure. Blue, green
and red lines in (b) and (c) correspond to the top, middle and bottom source-
receiver pairs (S2-R3, S6-R6 and S10-R10 in Fig. 1(b)), respectively. Solid,
dashed and dotted lines are results from 2017, 2018 and 2019, respectively.

Let us consider ultrasounds generated and recorded by the same source-receiver
pair in two different states: $\phi_{1}(t)$ and $\phi_{2}(t)$. The apparent
velocity change $dV/V_{a}$ on the time interval $[t_{1},t_{2}]$ is defined by

$dV/V_{a}=\arg\max_{\alpha}\frac{\int_{t_{1}}^{t_{2}}\phi_{1}((1+\alpha)t)\phi_{2}(t)dt}{\sqrt{\int_{t_{1}}^{t_{2}}\phi_{1}^{2}((1+\alpha)t)dt\int_{t_{1}}^{t_{2}}\phi_{2}^{2}(t)dt}},$

i.e., the amount of stretching that maximizes the cross-correlation (CC)
between the signals, and the DC is one minus CC, i.e.,

$\text{DC}=1-\max_{\alpha}\frac{\int_{t_{1}}^{t_{2}}\phi_{1}((1+\alpha)t)\phi_{2}(t)dt}{\sqrt{\int_{t_{1}}^{t_{2}}\phi_{1}^{2}((1+\alpha)t)dt\int_{t_{1}}^{t_{2}}\phi_{2}^{2}(t)dt}}.$

As an example, if the velocity decreases in the second state and if we focus
on a peak at $t_{0}$ in the wave field $\phi_{1}$, then the same peak should
arrive after $t_{0}$ in $\phi_{2}$. Therefore, to match the peak in these two
waveforms, $\alpha$ should be negative. That is to say, a dilation of the time
axis (a negative $dV/V_{a}$) indicates a decreasing velocity.

To quantitatively study the response of concrete to pressure changes in
different regions and to compare the state of the concrete over different
years, we compute $dV/V_{a}$ and DC for three source-receiver pairs between
the original state and other states during the inflation-deflation procedure:
the top pair (S2-R3 in Fig. 1), the middle pair (S6-R6) and the bottom pair
(S10-R10).

Fig. 2(c) shows the value of $dV/V_{a}$ for source-receiver pairs located at
the top (in blue), middle (in green) and bottom (in red) zones of the
structure, during inflation cycles from 2017 to 2019. When the the structure
is inflated, stress builds up in concrete, as shown in Fig. B1 and B2 in
appendix B, where we plot tangential stress distribution calculated with
linear elastic simulations. The stress distribution shows traction in the
upper part of the structure. This is consistent with what is observed in Fig.
2(c), which exhibits a decrease of velocity for the source-receiver pair
located at the top region. However, the apparent velocity change remains small
in the middle and marginal at the bottom, where the simulation shows a
constant stress or compression.

We have also checked the uncertainties of the apparent velocity changes.
Appendix A describes how the uncertainties are calculated, following the
formula introduced by weaver2011precision . These uncertainties are shown as
colored areas in Fig. A1, for the source-receiver pair located at the upper
region of the structure. One may see that they remain small despite strong
apparent velocity changes.

From 2017 to 2019 one notices an increase in the slope of $dV/V_{a}$, which
suggests that the sensitivity of concrete to the same pressure change becomes
larger. That is to say, the deterioration of the concrete is more pronounced
planes2013review .Another interesting phenomenon is that in 2017 the line
during inflation almost overlaps with the line during deflation which means
that expansion and contraction of the concrete are fully reversible, because
the concrete properties have not evolved significantly between the first ever
pressure test in 2017 and the one performed soon after presented here, which
is no more the case in 2019. This might be a good criterion to quantify the
damaging impact of the pressure tests, and to test the quality of the
concrete. Similar conclusions can be drawn for the middle and bottom pairs,
though a much lower resolution: the middle region experiences a contraction
and the pressure in the bottom region keeps almost unchanged. We leave the
precision study of the $dV/V_{a}$ due to changes in wave speed in the
Appendix.

In 2017, the concrete retrieves its initial state after the whole inflation-
deflation procedure, as indicated by DC (solid lines) in Fig. 2(b), which is
not the case in 2018 and 2019. We conclude from lines with the same color that
the change of DC becomes more irreversible as time goes by, which is caused by
the cumulative deterioration of the concrete (including micro-cracks and
cracks). The large value of DC between the initial and final states implies
that the damage induced by pressure changes is permanent: some cracks open and
close reversibly while others develop further and the crack surface is altered
such that they do not close exactly back to their initial state. What’s also
interesting is the bottom pair: we do not observe significant pressure changes
in the bottom region from $dV/V_{a}$, but they do suffer dramatic scattering
(or structural) changes (see Fig. 2(b)).

##  4 3D cartography of velocity and structural changes

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: (a)-(c) average relative velocity change between consecutive
measurements (one hour apart) at a depth of 5 cm during the second inflation
stage in 2017, 2018 and 2019, respectively; (d)-(f) results during the first
deflation stage. The white crosses and squares indicate positions of sources
and receivers, respectively. Black lines are the same as those in Fig. 1.

(a)

(b)

(c)

Figure 4: (a)-(c) average structural (scattering) change between two adjacent
measurements at the depth of 5 cm during the second inflation stage in 2017,
2018 and 2019, respectively. The white crosses and squares indicate positions
of sources and receivers, respectively. Black lines are the same as those in
Fig. 1.

The relation between the relative velocity change $\tfrac{dV}{V}(x)$ and
$dV/V_{a}$ is given by

$\frac{dV/V_{a}}{t}=\int\tfrac{dV}{V}(x)K(s,r,x,t)dx,$

and the formula connecting the scattering change and DC has a similar form

$\text{DC}=\int\tfrac{V\sigma(x)}{2}K(s,r,x,t)dx,$

where $V$ is the average velocity of the ultrasound and $\sigma$ is the local
scattering cross section density. The reader may refer to
planes2014decorrelation  for more details. These two formula share a common
sensitivity kernel $K(s,r,x,t)$, which is computed from the diffusive
approximation of the coda wave. The sensitivity kernel reflects the
probability that a wave emitted at $s$ at time $0$ and collected at $r$ at
time $t$ has passed through the position $x$. Here $t$ is in practice the
center of the time window where we evaluate $dV/V_{a}$ and DC experimentally.
Refer to xue2019locating  for a detailed discussion about computing the
sensitivity kernel and solving the integral equations.

We compute the velocity change between every two consecutive states and take
the average value from all consecutive pairs in the second inflation stage
(from 0.21 MPa to 0.41 MPa) and the first deflation stage (from 0.40 Mpa to
0.21 Mpa) to generate Fig. 3. Each step has, in average, a change of + or -
0.02 MPa. The decreasing (increasing, respectively) velocity is observed in
the inflation stage in the top region (in the lower middle region,
respectively), and the opposite phenomenon occurs during the deflation stage.
The magnitude in 2019 is much larger than in 2017, which indicates again the
cumulative damage of the concrete over years. We only present velocity changes
at a depth of 5 cm from the surface. We obtain similar results for other
depth, except a decrease of magnitude for deeper layers.

Fig. 4 shows the average structural changes, or more precisely the scattering
cross section density $\sigma(x,y)$ of the changes, between two consecutive
states during the second inflation stage. The units ($m^{2}/m^{3}$) can be
understood as a density of micro-cracks per volume unit. We observe an
increase of scattering change from 2017 to 2019. Our result shows that the
most active region is along the top construction joint. We also observe the
developing of the crack near dotted lines from 2017 to 2019.

##  5 Conclusion and discussion

In this work, we use diffuse ultrasound to study the pressure and structural
changes in the gusset area of a 1/3 replica of a nuclear power plant
containment wall. Changes are induced by the increase and decrease of the
inner pressure. We list the main results below:

  1. 1.

We are able to image the ultrasonic relative velocity change with a magnitude
as small as $10^{-4}$ thanks to the high sensitivity of diffuse ultrasound.

  2. 2.

The stretching of diffuse ultrasound and reconstructed velocity change
indicate that the concrete near the base is in compression and in traction in
the upper region in the inflation process, in agreement with mechanical
models.

  3. 3.

The increase of the slope of $dV/V_{a}$ versus internal pressure change from
2017 to 2019, in other words the sensitivity of $dV/V_{a}$ to internal
pressure, suggests that the concrete is progressively damaged and weakened.
That is to say, the concrete deteriorates over time.

  4. 4.

Although we do not detect significant pressure change in the base (bottom)
region, it could actually also suffer severe cumulative changes. The large DC
value between the initial and final states suggests that the concrete does not
restore to its initial state during the test, i.e., the damage is permanent
and irreversible, at least on the time scale of the test.

  5. 5.

The reconstructed structural changes (scattering cross section) shows active
regions during the inflation-deflation procedure that could correspond to
cracking inside the concrete.

##  6 Acknowledgements

This work was partially funded by the ANR ENDE (PIA 11-RSNR-0009 ) project,
and by the IDEX Université Grenoble Alpes innovation Grant (CND-US).

## Appendix A

In weaver2011precision , authors proposed to evaluate the uncertainties of
$dV/V_{a}$ due to the change in wave speed through the computation of its root
mean square, such that

$\text{rms}(dV/V_{a})=\sqrt{\frac{6\sqrt{\frac{\pi}{2}}T}{\omega_{c}^{2}(t_{2}^{3}-t_{1}^{3})}}\frac{\sqrt{1-\text{CC}^{2}}}{2\text{CC}}$
(1)

For our experiment, the parameters are as follows:

$\displaystyle\omega_{c}$ $\displaystyle=$ $\displaystyle 2\pi\times 1.5\times
10^{5}\text{ rad/sec},$ $\displaystyle T$ $\displaystyle=$
$\displaystyle\frac{\ln 10}{(2.2\times 10^{5}-0.8\times 10^{5})\pi},$
$\displaystyle t_{1}$ $\displaystyle=$ $\displaystyle 1.6\times 10^{-4}\text{
sec},$ $\displaystyle t_{2}$ $\displaystyle=$ $\displaystyle 4.8\times
10^{-4}\text{ sec}.$

Application of this formula to data recorded in 2017, 2018 and 2019 at the
upper source-receiver pair gives the uncertainties represented as colored
region in Fig. A1.

Figure A1: Apparent velocity change calculated in 2017 (solid line), 2018
(dashed line) and 2019 (dotted line) from source-receiver pair at the upper
region, and the corresponding uncertainties evaluated from Eq. (1) (colored
region).

## Appendix B

Here we evaluate numerically the state of stress inside the concrete around
the gusset. We perform elastic linear numerical simulations using the ASTER
code CodeAster , without taking into account the possible role of cracks. For
clarity, we only represent the difference of stress between the initial state
(additional internal pressure 0 MPa) and the final state (internal pressure
0.43 MPa). Fig. B1 and B2 represent a 2D vertical cut and a vertical profile
of this tangential stress change. Positive values correspond to traction,
where cracks are expected to develop and open, and where $dV/V_{a}$ is
expected to be negative. Negative values of stress correspond to compression
areas, where $dV/V_{a}$ are expected to be positive. Experimental results of
$dV/V_{a}$ are perfectly in agreement with this mechanical simulations,
demonstrating the ability of $dV/V_{a}$ obtained from ultrasonic coda waves to
cartography the state of stress inside the structure.

Figure B1: Numerical simulation of stress around the gusset (vertical cut). We
represent here the difference in tangential stress from the initial state (no
internal pressure) to the plateau at 0.4 MPa. The colored areas indicate
traction (form 0 to 0.3 MPa approximately). Black color indicates compression
or neutral areas. Figure B2: Tangential stress difference between the initial
(0 MPa) and final (0.43 MPa) state along a vertical line (y-axis) at a 5 cm
depth inside the concrete (from the outside surface). Negative values hold for
compression, positive values for traction.

## References

  * (1) Masson B and Alliard PM.  Objectives and design of the new experimental program vercors based on a 1/3 scaled PWR containment building.  _TINCE 2013_ 2013; .
  * (2) Mathieu J, Charpin L, Sémété P et al.  Temperature and humidity-driven ageing of the vercors mock-up.  In _Computational Modelling of Concrete Structures: Proceedings of the Conference on Computational Modelling of Concrete and Concrete Structures (EURO-C 2018)_ , volume 215.
  * (3) Henault JM, Laviron P, Desforges S et al.  How to characterize the airtightness of containment structure; overview of monitoring techniques tested on VeRCoRs mock up.  In _Congres on Technological Innovations in Nuclear Civil Engineering (TINCE), France, Paris-Saclay_.
  * (4) Wegler U and Sens-Schönfelder C.  Fault zone monitoring with passive image interferometry.  _Geophysical Journal International_ 2007; 168(3): 1029–1033.
  * (5) Brenguier F, Campillo M, Hadziioannou C et al.  Postseismic relaxation along the san andreas fault at parkfield from continuous seismological observations.  _science_ 2008; 321(5895): 1478–1481.
  * (6) Grêt A, Snieder R, Aster RC et al.  Monitoring rapid temporal change in a volcano with coda wave interferometry.  _Geophysical Research Letters_ 2005; 32(6).
  * (7) Sens-Schönfelder C and Wegler U.  Passive image interferometry and seasonal variations of seismic velocities at merapi volcano, indonesia.  _Geophysical research letters_ 2006; 33(21).
  * (8) Brenguier F, Shapiro NM, Campillo M et al.  Towards forecasting volcanic eruptions using seismic noise.  _Nature Geoscience_ 2008; 1(2): 126–130.
  * (9) Poupinet G, Ellsworth W and Frechet J.  Monitoring velocity variations in the crust using earthquake doublets: An application to the calaveras fault, california.  _Journal of Geophysical Research: Solid Earth_ 1984; 89(B7): 5719–5731.
  * (10) Snieder R, Grêt A, Douma H et al.  Coda wave interferometry for estimating nonlinear behavior in seismic velocity.  _Science_ 2002; 295(5563): 2253–2255.
  * (11) Snieder R.  The theory of coda wave interferometry.  _Pure and Applied geophysics_ 2006; 163(2-3): 455–473.
  * (12) Lobkis OI and Weaver RL.  Coda-wave interferometry in finite solids: Recovery of P-to-S conversion rates in an elastodynamic billiard.  _Physical review letters_ 2003; 90(25): 254302.
  * (13) Planès T, Larose E, Rossetto V et al.  Imaging multiple local changes in heterogeneous media with diffuse waves.  _The Journal of the Acoustical Society of America_ 2015; 137(2): 660–667.
  * (14) Larose E, de Rosny J, Margerin L et al.  Observation of multiple scattering of khz vibrations in a concrete structure and application to monitoring weak changes.  _Physical Review E_ 2006; 73(1): 016609.
  * (15) Larose E and Hall S.  Monitoring stress related velocity variation in concrete with a 2$\times$ 10- 5 relative resolution using diffuse ultrasound.  _The Journal of the Acoustical Society of America_ 2009; 125(4): 1853–1856.
  * (16) Schurr DP, Kim JY, Sabra KG et al.  Damage detection in concrete using coda wave interferometry.  _NDT & E International_ 2011; 44(8): 728–735.
  * (17) Hilloulin B, Zhang Y, Abraham O et al.  Small crack detection in cementitious materials using nonlinear coda wave modulation.  _NDT & E International_ 2014; 68: 98–104.
  * (18) Planès T and Larose E.  A review of ultrasonic coda wave interferometry in concrete.  _Cement and Concrete Research_ 2013; 53: 248–255.
  * (19) Pacheco C and Snieder R.  Time-lapse travel time change of multiply scattered acoustic waves.  _The Journal of the Acoustical Society of America_ 2005; 118(3): 1300–1310.
  * (20) Larose E, Planes T, Rossetto V et al.  Locating a small change in a multiple scattering environment.  _Applied Physics Letters_ 2010; 96(20): 204101.
  * (21) Planès T, Larose E, Margerin L et al.  Decorrelation and phase-shift of coda waves induced by local changes: multiple scattering approach and numerical validation.  _Waves in Random and Complex Media_ 2014; 24(2): 99–125.
  * (22) Xue Q, Larose E and Moreau L.  Locating structural changes in a multiple scattering domain with an irregular shape.  _The Journal of the Acoustical Society of America_ 2019; 146(1): 595–602.
  * (23) Kanu C and Snieder R.  Numerical computation of the sensitivity kernel for monitoring weak changes with multiply scattered acoustic waves.  _Geophysical Supplements to the Monthly Notices of the Royal Astronomical Society_ 2015; 203(3): 1923–1936.
  * (24) Zhang Y, Planes T, Larose E et al.  Diffuse ultrasound monitoring of stress and damage development on a 15-ton concrete beam.  _The Journal of the Acoustical Society of America_ 2016; 139(4): 1691–1701.
  * (25) Zhan H, Jiang H and Jiang R.  Three-dimensional images generated from diffuse ultrasound wave: detections of multiple cracks in concrete structures.  _Structural Health Monitoring_ 2020; 19(1): 12–25.
  * (26) Zhang Y, Larose E, Moreau L et al.  Three-dimensional in-situ imaging of cracks in concrete using diffuse ultrasound.  _Structural Health Monitoring_ 2018; 17(2): 279–284.
  * (27) Zhang Y, Abraham O, Grondin F et al.  Study of stress-induced velocity variation in concrete under direct tensile force and monitoring of the damage level by using thermally-compensated coda wave interferometry.  _Ultrasonics_ 2012; 52(8): 1038–1045.
  * (28) Weaver RL, Hadziioannou C, Larose E et al.  On the precision of noise correlation interferometry.  _Geophysical Journal International_ 2011; 185(3): 1384–1392.
  * (29) Electricité de France.  Finite element $\bf\it code\\_aster$, structures and thermomechanics analysis for studies and research.  Open source on www.code-aster.org, 1989–2021.


# Manipulation of non-linear heat currents in the dissipative Anderson-
Holstein model

Bitan De  Bhaskaran Muralidharan<EMAIL_ADDRESS>Department of Electrical
Engineering, Indian Institute of Technology Bombay, Powai, Mumbai-400076,
India

(August 27, 2024)

###### Abstract

The anomalous behavior of electron induced phonon transport is investigated
using an Anderson-Holstein based dissipative quantum dot setup under two
relevant bias situations: (a) a voltage bias in the absence of an electronic
temperature gradient and (b) an electronic temperature gradient at zero
voltage. It is shown that the direction of phonon transport in the non-linear
regime is different in the two cases since the first case facilitates the
accumulation of phonons in the dot and the second case leads to the absorption
of phonons in the dot. In the linear regime, both the phonon and electronic
transport get decoupled and Onsager’s symmetry is verified. We explain the
observed cumulative effects of voltage and electronic temperature gradients on
the non-linear phonon currents by introducing a new transport coefficient that
we term as the electron induced phonon thermal conductivity. It is
demonstrated that under suitable operating conditions in Case (a) the dot can
pump in phonons into the hotter phonon reservoirs and in Case (b) the dot can
extract phonons out of the colder phonon reservoirs. Finally, we elaborate on
how the non-linear electronic heat current can be stimulated and controlled by
engineering the temperature of the phonon reservoirs even under vanishing
effective electron flow.

##  I Introduction

The physics of electron-phonon interaction has gained significant interest
Giustino (2017) with its fingerprints being captured in diverse physical
systems like superconductorsLanzara et al. (2001); Walborn et al. (2006); Lee
et al. (2006); Hong et al. (2009); Bianchi et al. (2010), cavity-coupled
mesoscopic conductorsHartke et al. (2018); Frey et al. (2012), optomechanical
systemsAres et al. (2016); Park et al. (2004) and electronic devicesLaird et
al. (2012); Javey et al. (2003), to name a few. So far, the research in this
area is primarily carried out along two predominant directions: (a)
exploration of novel experimental techniques to detect such interactions Yu et
al. (2004a); Rakhmilevitch et al. (2014); Ward et al. (2010); Ioffe et al.
(2008); Benyamini et al. (2014) and (b) theoretical studies on the effect of
electron-phonon coupling on quantum phenomena in the nanoscale Härtle and
Thoss (2011); Leijnse et al. (2010); Galperin et al. (2004, 2007); De and
Muralidharan (2016); Agarwalla et al. (2015); Arrachea et al. (2014); Sothmann
et al. (2014); Siddiqui et al. (2006); Mitra et al. (2004). Typically in
microscopic thermoelectric setups, the former direction concerns materials
synthesis, device fabrication and the integration at the systems level Kim et
al. (2014); Reddy et al. (2007). The principal goal in this case, is to
control the flow of waste heat to optimize the efficiency of energy
conversionRoßnagel et al. (2016); Steele et al. (2009); Lassagne et al.
(2009); Woodside and McEuen (2002). On the other hand, the latter direction
deals with the dynamics of charge and heat transport from a quantum transport
perspective. From the application point of view, this approach aims to propose
new devices to realize novel paradigms like phonon computation Jiang et al.
(2015); Wang and Li (2007, 2008); Li et al. (2012); Joulain et al. (2016) .
Many notable works have established that in the zero-dimensional systems such
as quantum dots,Mahan and Sofo (1996); Kubala and König (2006); Kubala et al.
(2008) the interplay of electrons and phonons gives rise to some unique
features manifesting as transport signatures Park et al. (2000); Leroy et al.
(2004); LeRoy et al. (2005); Yu et al. (2004b); Sapmaz et al. (2006). Here,
typical to quantum transport set ups, electronic and phonon currents are
driven by applying a bias voltage and/or a thermal bias across the contacts.
While voltage induced phonon transport has been the subject of many
theoretical worksSchmidt and Komnik (2009); Urban et al. (2010); Egger and
Gogolin (2008); Seoane Souto et al. (2014); Ness (2014); Wilner et al. (2014),
understanding the non-equilibrium phonon transport under an electronic
temperature bias has not been explored so far, which will be the objective of
this paper.

Figure 1: Non-equilibrium charge current(shown in red bold arrow) is driven
through the central dot in two ways: either by applying (a) voltage bias or
(b) by applying electronic temperature gradient across the contacts $L$ and
$R$. Dot phonon current (shown by orange bold arrows) is set up by creating a
thermal gradient across the reservoirs $H$ and $C$. In both (a) and (b) the
dot temperature is monitored by phonon thermometer bath. (c) The state
transition diagram in the electronic Fock space. (d) The state transition
diagram in the many body electron-phonon Fock space. Direct electron
tunnelings (depicted by thin black double arrows) occur between states
$|n,q\rangle$ and $|n\pm 1,q\rangle$, while phonon assisted electron
tunnelings (depicted by gray thin double arrows) happen between states
$|n,q\rangle$ and $|n\pm 1,q^{\prime}\rangle$. Reservoir assisted phonon
transitions (represented by black bold arrows) take place between the states
$|n,q\rangle$ and $|n,q\pm 1\rangle$.

Our objective in this work is to compare the nature of the electronically
stimulated phonon transport in two cases via: Case (a): a voltage bias in the
absence of an electronic temperature bias and Case (b): an electronic
temperature bias at zero voltage. It is shown by the simulation framework that
the polarity of phonon current in the nonlinear regime is opposite in the two
cases since the first case facilitates the accumulation of phonons in the dot
and the second case leads to the absorption of phonons in the dot.
In Case (a), surplus phonons heat up the dot and are suitably extracted by the
phonon reservoirs kept in equilibrium. In Case (b), the dot cools down for
being phonon deficient and the reservoirs pump in phonons into it. However in
the linear regime, phonon and electron transport get decoupled and Onsager’s
reciprocity is verified.
In our previous work De and Muralidharan (2018), we explored the non-linear
phonon Peltier effect by noticing the voltage dependence of the non-linear
phonon current. In the current work, we propose a new parameter called,
electron-assisted phonon thermal conductivity to signify the mutual dependence
of phonon currents on the voltage and electronic temperature gradient.
We also notice that the polarity of phonon transport can be controlled by
differing the temperature of phonon reservoirs from the temperature of
electronic contacts. In Case (a), phonons flow from the reservoirs to the dot
at low voltage when the reservoirs are hotter than the contacts. At large
voltage, the polarity of phonon flow is reversed and the phonons are pumped
into the hotter reservoirs. An opposite situation is observed in Case (b)
given the reservoirs are colder than the contacts. In this case, phonons relax
from the dot to the reservoirs at low electronic temperature gradient. When
the temperature gradient is large, the direction of phonon flow is reversed
and phonons are extracted from the colder reservoirs. In both cases, the
direction of phonon current is verified by estimating the dot temperature with
a phonon thermometer bath weakly coupled to the dot. Finally, we elaborate on
how the non-linear electronic heat current can be stimulated and controlled by
engineering the temperature of the phonon reservoirs even even under vanishing
effective electron flow.
This paper is organized as follows: Section II introduces the device model and
formulates the transport equations. In the Sec. III A, we focus on the
comparative study of the characteristics of the phonon currents in the two
cases under consideration. Section III B introduces the electron assisted
phonon thermal conductivity coefficient and compares it with the conventional
electronic thermal conductivity. In the Sec. III C, we describe the
engineering of the temperature of the phonon reservoirs to modulate the flow
of the phonon and electronic heat currents. In the Sec. IV, we summarize our
results and conclude.

##  II Physics and Formulation

The schematics of the set ups to be studied here are depicted in Fig.1 and
Fig.1. Each of them typically comprises a single quantum dot weakly coupled to
the electronic contacts $\alpha_{1}$ ($\alpha_{1}\in L,R$), and the phonon
reservoirs $\alpha_{2}$ ($\alpha_{2}\in H,C$). The dot is described by the
dissipative Anderson-Holestein model which has both electronic and phonon
degrees of freedom interacting with each other. The electronic current is
driven through the device in two ways: Case (i): by applying a voltage bias
(shown in Fig.1(a)) across the contacts $L$ and $R$ and Case (ii) by applying
an electronic temperature gradient (shown in Fig.1) across the contacts $L$
and $R$ . The phonon current is set up by applying a thermal bias across the
phonon reservoirs $H$ and $C$. In addition, the electron current gives rise to
a phonon current even in the absence of a temperature bias across the
reservoirs $H$ and $C$, provided the electron-phonon coupling is finite. The
direction of phonon flow depends on the dot temperature ($T_{M}$), which is
estimated by a phonon thermometer bath weakly coupled to the dot.

####  II.0.1 Model Hamiltonian

The composite system Hamiltonian $\hat{H}$ can be written as
$\hat{H}=\hat{H}_{D}+\hat{H}_{\alpha_{1}}+\hat{H}_{\alpha_{2}}+\hat{H}_{\alpha_{1}D}+\hat{H}_{\alpha_{2}D}$.
Here $\hat{H}_{D}$, $\hat{H}_{\alpha_{1}}$, $\hat{H}_{\alpha_{2}}$ denote the
respective Hamiltonians of the dot, the contacts and the reservoirs, while
$\hat{H}_{\alpha_{1}D}$ and $\hat{H}_{\alpha_{2}D}$ represent the dot-to-
contact electronic tunneling processes and the dot-to-reservoir phonon
relaxation processes respectively. The dot Hamiltonian is further divided as,

$\begin{gathered}\hat{H}_{D}=\hat{H}_{el}+\hat{H}_{ph}+\hat{H}_{el-
ph}.\end{gathered}$ (1)

The electronic part ($\hat{H}_{el}$) consists of a single spin degenerate
energy level with an on-site energy $\epsilon$ and a Coulomb interaction
energy $U$ for double occupancy. The phonon part ($\hat{H}_{ph}$) comprises of
a single phonon mode $\nu$ with angular frequency $\omega_{\nu}$. Inside the
dot, the electrons and phonons interact through a dimensionless parameter
$\lambda_{\nu}$. Hence, $H_{el}$, $H_{ph}$,$H_{el-ph}$ are given by,

$\begin{gathered}\hat{H}_{el}=\bigg{(}\sum\limits_{\sigma}\epsilon
d_{\sigma}^{\dagger}d_{\sigma}+Ud_{\uparrow}^{\dagger}d_{\uparrow}d_{\downarrow}^{\dagger}d_{\downarrow}\bigg{)},\\\
\hat{H}_{ph}=\hbar\omega_{\nu}\hat{b}_{\nu}^{\dagger}\hat{b}_{\nu},\\\
\hat{H}_{el-
ph}=\sum\limits_{\sigma}\lambda_{\nu}\hbar\omega_{\nu}d_{\sigma}^{\dagger}d_{\sigma}(\hat{b}_{\nu}^{\dagger}+\hat{b}_{\nu}).\end{gathered}$
(2)

In the above expressions, $d_{\sigma}^{\dagger}$($d_{\sigma}$) and
$b_{\nu}^{\dagger}$($b_{\nu}$) are the respective creation(annihilation)
operators of the electrons and phonons. The electronic contacts $\alpha_{1}$
characterize a macroscopic body of non-interacting electrons in the momentum
eigen state $\alpha_{1}k$ and spin $\sigma^{\prime}$ with energies
$\epsilon_{\alpha_{1}k\sigma^{\prime}}$. Similarly, the phonon reservoirs
$\alpha_{2}$ comprise of numerous phonon modes $\alpha_{2}r$ with an angular
frequency $\omega_{\alpha_{2}r}$. They are defined by the following
Hamiltonians,

$\begin{gathered}\hat{H}_{\alpha_{1}}=\sum\limits_{\alpha_{1}\in
L,R}\sum\limits_{\alpha_{1}k\sigma^{\prime}}\hat{c}_{\alpha_{1}k\sigma^{\prime}}^{\dagger}\hat{c}_{\alpha_{1}k\sigma^{\prime}},\\\
\hat{H}_{\alpha_{2}}=\sum\limits_{\alpha_{2}\in
H,C}\sum\limits_{\alpha_{2}}\hat{B}_{\alpha_{2}r}^{\dagger}\hat{B}_{\alpha_{2}r},\end{gathered}$
(3)

where
$\hat{c}_{\alpha_{1}k\sigma^{\prime}}^{\dagger}(\hat{c}_{\alpha_{1}k\sigma^{\prime}})$
creates (annihilates) an electron with momentum $k$ and spin $\sigma^{\prime}$
in the contact $\alpha_{1}$, and
$\hat{B}_{\alpha_{2}r}^{\dagger}(\hat{B}_{\alpha_{2}r})$ creates (annihilates)
a phonon of angular frequency $\omega_{\alpha_{2}r}$ in the reservoir
$\alpha_{2}$. If the electrons in the dot are coupled to the electrons in the
contacts $\alpha_{1}$ through an energy
$\tau_{\alpha_{1}k\sigma^{\prime}\sigma}^{el}$ and the phonons in the dot are
coupled to the phonons in the reservoirs $\alpha_{2}$ through an energy
$\tau_{\alpha_{2}r\nu}^{ph}$ then the coupling Hamiltonians
$\hat{H}_{\alpha_{1}D}$ and $\hat{H}_{\alpha_{2}D}$ are given as follows

$\begin{gathered}\hat{H}_{\alpha_{1}D}=\sum\limits_{k\sigma^{\prime}\sigma}[\tau_{\alpha_{1}k\sigma^{\prime}\sigma}^{el}\hat{c}_{\alpha_{1}k\sigma^{\prime}}^{\dagger}\hat{d}_{\sigma}+H.c],\\\
\hat{H}_{\alpha_{2}D}=\sum\limits_{\nu,\alpha_{2}}\tau_{\alpha_{2}r\nu}^{ph}(\hat{B}_{\alpha
2r}^{\dagger}+\hat{B}_{\alpha
2r})(\hat{b}_{\nu}^{\dagger}+\hat{b}_{\nu}).\end{gathered}$ (4)

In this work, we assume that the phonon transition processes are mode
independent and the electron tunneling processes are momentum conserved and
spin independent. With this approximation we can rewrite
$\tau_{\alpha_{1}k\sigma^{\prime}\sigma}^{el}$ and
$\tau_{\alpha_{2}r\nu}^{ph}$ simply as $\tau^{el}$ and $\tau^{ph}$
respectively.
The dot Hamiltonian is diagonalized by the polaron transformation (such that
$\hat{H}_{D}\rightarrow\tilde{\hat{H}}_{D}=e^{S}\hat{H}_{D}e^{-S}$, where
$S=\sum_{\nu}\lambda_{\nu}[\hat{b}_{\nu}-\hat{b}_{\nu}^{\dagger}$]). It
transforms the dot fermionic operators
$\hat{d}_{\sigma}^{\dagger}$($\hat{d}_{\sigma}$) and leads to the
renormalization of $\epsilon$ and $U$, such that
$\tilde{\epsilon}=\epsilon-\lambda_{\nu}^{2}\hbar\omega_{\nu}$ and
$\tilde{U}=U-2\lambda_{\nu}^{2}\hbar\omega_{\nu}$. The eigen states of
$\tilde{\hat{H}}_{D}$ are represented by $|n,q\rangle$, where $n$ and $q$ are
the electron and phonon number of the eigen states. The energy eigen values of
those states become $E_{(n,q)}=\tilde{E}_{n}+q\hbar\omega_{\nu}$, where,
$\tilde{E}_{0}=0,\tilde{E}_{1}=\tilde{\epsilon},\tilde{E}_{2}=2\tilde{\epsilon}+\tilde{U}$,
corresponding to the $n=0,1$ and $2$ electron number spaces respectively. This
transformation also renormalizes the electronic tunneling energy $\tau^{el}$,
such that
$\tilde{\tau}^{el}=\tau^{el}\exp[-\lambda_{\nu}(\hat{b_{\nu}}-\hat{b_{\nu}}^{\dagger})]$.
However, the phonon coupling energy $\tau^{ph}$ is left unaltered due to the
polaron transformation since the operator $S$ commutes with the phonon
operators $b_{\nu}^{\dagger}$($b_{\nu}$) of the dot .
With the derived expressions of $\tilde{\tau}_{el}$ and $\tau_{ph}$, we can
evaluate the dot-to-contact electron tunneling rate $\gamma_{\alpha_{1}}$ and
the dot-to-phonon relaxation rate $\beta_{\alpha_{2}}$ using the Fermi’s
golden rule. They are represented as:
$\gamma_{\alpha_{1}}=\frac{2\pi}{\hbar}\sum_{\alpha_{1}}\absolutevalue{\tilde{\tau}_{\alpha_{1}}^{el}}^{2}\rho_{\alpha_{1}\sigma^{\prime}}$
and
$\beta_{\alpha_{2}}=\frac{2\pi}{\hbar}\sum_{\alpha_{2}}\absolutevalue{\tau_{\alpha_{2}}^{ph}}^{2}D_{\alpha_{2}}$,
where $\rho_{\alpha_{1}\sigma^{\prime}}$ and $D_{\alpha_{2}}$ are the constant
electron and phonon density of states associated with the contacts
$\alpha_{1}$ and the reservoirs $\alpha_{2}$ respectively. The study of the
model Hamiltonian enables us to formulate the transport methodology, which we
will analyze in the next subsection.

####  II.0.2 Transport methodology

Before we start the formulation of transport equations, it is imperative to
explain the important approximations we have considered in this work. First,
we set the rate of dot-to-reservoir phonon relaxation processes much lower
than the rate of dot-to-contact electron tunneling processes
($\hbar\gamma_{\alpha_{1}}>>\hbar\beta_{\alpha_{2}}$) to rule out the system
dampingBraig and Flensberg (2003). In this limit, the phonon currents
emanating from the terminal phonon reservoirs remain uncorrelated and hence
can be computed independentlySegal and Nitzan (2005a); Segal (2006); Segal and
Nitzan (2005b). Second, we set
$\hbar\gamma_{\alpha_{1}},\hbar\beta_{\alpha_{2}}<<k_{B}T$, so that the
transport of the electrons and phonons through a single quantum dot weakly
coupled to the electronic contacts and the phonon reservoirs can be described
in the sequential tunneling regime in which charge and phonon currents are
calculated via the quantum Master equationMuralidharan et al. (2006);
Muralidharan and Datta (2007); Beenakker (1991); Muralidharan and Grifoni
(2013, 2012). Last, we neglect the overlap of the adjacent phonon sidebands,
by setting the energy gap between the sidebands larger than the tunneling
induced broadening of energy levels
($\hbar\omega_{\nu}>>\hbar\gamma_{\alpha_{1}}$)Timm (2008). With this
assumption, the Markov approximation is justified and two consecutive electron
tunneling processes are completely uncorrelated since the memory of
macroscopic contacts and reservoirs is negligiblePiovano (2012). Hence, the
diagonal terms of the system density matrix get decoupled from the off-
diagonal terms and the Master equation reduces to the rate equationBraig and
Brouwer (2005); König and Martinek (2003); Braun et al. (2004).
The tunneling rate between two eigen states $|n,q\rangle$ and $|n\pm
1,q^{\prime}\rangle$ is determined by the contact Fermi function of the energy
difference of the two states and it is given by:

$\begin{gathered}R_{(n,q)\rightarrow(n+1,q^{\prime})}=\sum\limits_{\alpha_{1}\in
L,R}\gamma_{\alpha_{1}}|\langle
n,q|\hat{d}_{\sigma}|n+1,q^{\prime}\rangle|^{2}\\\ \times
f_{\alpha_{1}}(E_{n+1,q^{\prime}}-E_{n,q})\end{gathered},$ (5)
$\begin{gathered}R_{(n,q)\rightarrow(n-1,q^{\prime})}=\sum\limits_{\alpha_{1}\in
L,R}\gamma_{\alpha_{1}}|\langle
n,q|\hat{d}_{\sigma}^{\dagger}|n-1,q^{\prime}\rangle|^{2}\\\
\times[1-f_{\alpha_{1}}(E_{n,q}-E_{n-1,q^{\prime}})],\end{gathered}$ (6)

where
$f_{\alpha_{1}}(\zeta)=1/(1+exp(\frac{\zeta-\mu_{\alpha_{1}}}{k_{B}T_{\alpha_{1}}}))$
is the Ferni-Dirac distribution function of the contact $\alpha_{1}$ with
chemical potential $\mu_{\alpha_{1}}$ and temperature $T_{\alpha_{1}}$. The
phonon relaxation process between the reservoirs and the dot lead to the
transition between two eigen states $|n,q\rangle$ and $|n,q\pm 1\rangle$ and
the relaxation rate follows the Boltzmann ratio:

$\begin{split}R_{(n,q)\rightarrow(n,q+1)}=\sum\limits_{\alpha_{2}\in
H,C}\beta_{\alpha_{2}}(q+1)exp\bigg{(}-\frac{\hbar\omega_{\nu}}{k_{B}T_{\alpha_{2}}}\bigg{)},\end{split}$
(7) $\begin{split}R_{(n,q)\rightarrow(n,q-1)}=\sum\limits_{\alpha_{2}\in
H,C}\beta_{\alpha_{2}}(q+1).\end{split}$ (8)

Using the rate equations, the master equation for the probabilities
$P_{(n,q)}$ of the many-body electron-phonon states $|n,q\rangle$ take the
following form:

$\begin{gathered}\frac{dP_{(n,q)}}{dt}=\sum\limits_{q^{\prime}}^{N_{q}}\bigg{[}R_{(n^{\prime},q^{\prime})\rightarrow(n,q)}^{el}P_{(n^{\prime},q^{\prime})}-R_{(n,q)\rightarrow(n^{\prime},q^{\prime})}^{el}P_{(n,q)}\bigg{]}\\\
+\bigg{[}R_{(n^{\prime},q^{\prime})\rightarrow(n,q)}^{ph}P_{(n^{\prime},q^{\prime})}-\\\
R_{(n,q)\rightarrow(n^{\prime},q^{\prime})}^{ph}P_{(n,q)}\bigg{]}\delta(n\pm
1,n^{\prime})\delta(q\pm 1,q^{\prime}).\end{gathered}$ (9)

In the steady state, the derivative in the left hand side of (9) vanishes and
the steady state probabilities $P_{(n,q)}$ can be determined by solving the
algebric equations. We use the probabilities to compute the charge currents,
the electronic heat currents associated with the contacts $\alpha_{1}$ and the
phonon heat currents associated with reservoirs $\alpha_{2}$. They are given
as,

$\begin{gathered}I_{\alpha_{1}}=\sum\limits_{q=0}^{N_{q}}\sum\limits_{q^{\prime}=0}^{N_{q}}-q\bigg{[}R_{(n+1,q^{\prime})\rightarrow(n,q)}^{el_{\alpha_{1}}}P_{(n+1,q^{\prime})}\\\
-R_{(n+1,q^{\prime})\rightarrow(n,q)}^{el_{\alpha_{1}}}P_{(n,q)}\bigg{]},\end{gathered}$
(10)
$\begin{gathered}I_{el_{\alpha_{1}}}^{Q}=\sum\limits_{q=0}^{N_{q}}\sum\limits_{q^{\prime}=0}^{N_{q}}(E_{n+1,q^{\prime}}-E_{n,q}-\mu_{\alpha})\\\
\bigg{[}R_{(n+1,q^{\prime})\rightarrow(n,q)}^{el_{\alpha_{1}}}P_{(n+1,q^{\prime})}\\\
-R_{(n+1,q^{\prime})\rightarrow(n,q)}^{el_{\alpha_{1}}}P_{(n,q)}\bigg{]},\end{gathered}$
(11)
$\begin{gathered}I_{ph_{\alpha_{2}}}^{Q}=\sum\limits_{q=0}^{N_{q}}\sum\limits_{q^{\prime}=0}^{N_{q}}\hbar\omega_{\nu}\bigg{[}R_{(n,q)\rightarrow(n^{\prime},q^{\prime})}^{ph_{\alpha_{2}}}P_{(n,q)}\\\
-R_{(n^{\prime},q^{\prime})\rightarrow(n,q)}^{ph_{\alpha_{2}}}P_{(n^{\prime},q^{\prime})}\bigg{]}\delta(n,n^{\prime})\delta(q\pm
1,q^{\prime}).\end{gathered}$ (12)

Now it is evident from (5) and (6) that the transition between the states
$|n,q\rangle$ and $|n\pm 1,q^{\prime}\rangle$ leads to the net phonon
generation (or absorption) in the dot. In the next subsection we will focus on
the role of $\lambda_{\nu}$ in controlling the phonon generation(or
absorption). Considering the law of charge conservation, in the rest of the
paper we will denote $I_{L}=-I_{R}=I$.

####  II.0.3 Effect of interaction on phonon transport

The effective electron tunneling rate ($\gamma_{\alpha_{1}}^{eff}$) between
the states $|n,q\rangle$ and $|n\pm 1,q^{\prime}\rangle$ is modified by the
Frank-Condon overlapping factor between the two statesKoch and Von Oppen
(2005); Koch et al. (2006); Koch (2006). The effective electron tunneling rate
is given as

$\begin{gathered}\gamma_{\alpha_{1}}^{eff}=\gamma_{\alpha_{1}}|FC_{q,q^{\prime}}|^{2}=\gamma_{\alpha_{1}}\bigg{[}|\langle
n,q|\hat{d}^{\dagger}|n^{\prime},q^{\prime}\rangle|^{2}\delta(n^{\prime},n-1)\\\
+|\langle
n,q|\hat{d}|n^{\prime},q^{\prime}\rangle|^{2}\delta(n^{\prime},n+1)\bigg{]},\end{gathered}$
(13)

where
$|FC_{q,q^{\prime}}|^{2}=exp(-\lambda_{\nu}^{2})\frac{k!}{K!}\lambda^{2(K-k)}[L_{k}^{K-k}(\lambda_{\nu}^{2})]^{2}$
is the Frank-Condon factor between the two states with phonon number $q$ and
$q^{\prime}$ and $L_{k}^{K-k}$ is the associated Laguerre polynomial with
$k=min(q,q^{\prime})$ and $K=max(q,q^{\prime})$. The net phonon generation (or
absorption) in the dot due to the transition between the states $|n,q\rangle$
and $|n\pm 1,q\rangle$, takes place at a rateSiddiqui et al. (2006)

$\begin{gathered}GE_{ph}^{\alpha_{1}}=\sum\limits_{n,q}\sum\limits_{n\pm
1}(q^{\prime}-q)P_{n,q}R_{(n,q)\rightarrow(n\pm
1,q^{\prime})}^{\alpha_{1}},\end{gathered}$ (14)

It is evident that in the limit
$\hbar\gamma_{\alpha_{1}}>>\hbar\beta_{\alpha_{2}}$, the phonon distribution
in the dot is primarily determined by the electron transport. When $q\neq
q^{\prime}$, the phonons generate (or get absorbed) in the dot and the average
phonon number in the dot $\langle N_{ph}\rangle=\sum\limits_{n,q}qP_{n,q}$
deviates from the equilibrium phonon distribution $\langle N_{ph}^{eq}\rangle$
of the phonon reservoirs. The excess phonons are extracted by the reservoirs
at a rateSiddiqui et al. (2006)

$\begin{gathered}\begin{split}RE_{ph}^{\alpha_{2}}=\sum\limits_{q=0}^{N_{q}}\sum\limits_{q^{\prime}=0}^{N_{q}}\hbar\omega_{\nu}\bigg{[}R_{(n,q)\rightarrow(n^{\prime},q^{\prime})}^{ph_{\alpha_{2}}}P_{(n,q)}\\\
-R_{(n^{\prime},q^{\prime})\rightarrow(n,q)}^{ph_{\alpha_{2}}}P_{(n^{\prime},q^{\prime})}\bigg{]}\delta(n,n^{\prime})\delta(q\pm
1,q^{\prime})\\\ RE_{ph}^{\alpha_{2}}=\beta_{\alpha_{2}}\frac{\langle
N_{ph}\rangle-\langle N_{ph}\rangle^{eq}}{1+\langle
N_{ph}\rangle^{eq}}.\end{split}\end{gathered}$ (15)

It is evident from (12) and (15) that the phonon current
$I_{ph}^{Q}=\hbar\omega RE_{ph}^{\alpha_{2}}$. It is noticed from (13), that
when $\lambda_{\nu}=0$, $\gamma_{\alpha_{1}}^{eff}$ vanishes unless
$q=q^{\prime}$. In this case, $\langle N_{ph}\rangle$ equals with $\langle
N_{ph}^{eq}\rangle$ and $I_{ph_{\alpha_{2}}}^{Q}$ vanishes. On the other way,
when $\lambda_{\nu}$ is non-zero, $\langle N_{ph}\rangle$ deviates from
$\langle N_{ph}^{eq}\rangle$ and $I_{ph_{\alpha_{2}}}^{Q}$ becomes finite. The
deviation of $\langle N_{ph}\rangle$ from $\langle N_{ph}^{eq}\rangle$ causes
the variation of dot temperature $T_{M}$ from the equilibrium reservoir
temperature $T_{\alpha_{2}}$. The dot temperature can be computed from the
Boltzmann ratio by with a quasi-equilibrium approximationArrachea et al.
(2014) as:

$\begin{gathered}T_{M}=\frac{\hbar\omega_{\nu}}{k_{B}}\bigg{[}ln\bigg{(}\frac{P_{n,q}}{P_{n,q+1}}\bigg{)}\bigg{]}^{-1}\end{gathered}$
(16)

We will notice in the next section that the electron induced phonon current
shows different characteristics depending on whether the electron flow is
stimulated by a voltage bias or by an electronic temperature bias. From now
on, we will denote $\lambda_{\nu}$, $\omega_{\nu}$, $GE_{ph}^{\alpha_{1}}$ and
$RE_{ph}^{\alpha_{2}}$ simply as $\lambda$, $\omega$, $G_{ph}$ and $R_{ph}$.
Also, unless otherwise mentioned, the electronic contact coupling and
reservoir phonon couplings are assumed to be symmetric (i.e
$\gamma_{\alpha_{1}}=\gamma$, $\gamma_{\alpha_{1}}^{eff}=\gamma^{eff}$ and
$\beta_{\alpha_{2}}=\beta$, where $\alpha_{2}\in H,C$).

Figure 2: Electron induced phonon transport at zero phonon temperature
gradient. Modified setups in two operating limits: (a) Current is driven by a
voltage bias when $\Delta T_{el}=0$. (b) Current is driven set up by an
electronic temperature bias when $V=0$. (c) Variation of $I_{Q}^{ph}$ with
voltage for different $\lambda$ in Case (a). (d) Variation of $I_{Q}^{ph}$
with $\Delta T_{el}$ for different $\lambda$ in Case (b). This figure
establishes that $V$ and $\Delta T_{el}$ decide the direction of $I_{Q}^{ph}$
and $\lambda$ works like a phonon switch.

##  III Results

###  III.1 Electron-phonon coupled transport

In this section we focus on the electron-phonon coupled transport in our
setup. Each operating point is signified by a voltage bias, an electronic
temperature gradient ($\Delta T_{el}=T_{L}-T_{R}$) applied across the contacts
and a phonon temperature gradient ($\Delta T_{ph}=T_{H}-T_{C}$) applied
between the reservoirs. In our analysis, we set
($k_{B}$T$>>\hbar\omega>>\hbar\gamma,\hbar\beta$) the tunnel induced
broadening of the energy levels in the dot small to ensure that the levels do
not overlap. Under this condition, transport is described in the sequential
tunneling limit, where charge and heat current are calculated via rate
equation. Additionally we assume the dot functions as n-type, i.e
$\epsilon>>\mu$.
First we demonstrate the effect of electron current on $I_{Q}^{ph}$ in two
regimes: Case (a) the electron current is driven by a voltage bias when
$\Delta T_{el}$ ($\Delta T_{el}=T_{L}-T_{R}$) is zero and Case (b) the
electron current is set up by an electronic temperature gradient ($\Delta
T_{el}=T_{L}-T_{R}\neq 0$) at zero voltage. Each case operates at zero $\Delta
T_{ph}$ to ensure that only electron current influences the phonon
distribution. To put things simple, reservoirs $H$ and $C$ are merged to a
single reservoir $B$ with temperature $T_{B}$ and relaxation rate $\beta$. The
modified setups in Case (a) and (b) are depicted in Fig. 2 and 2 respectively.
In each operating point, the direction of phonon current is decided by the
gradient between $T_{B}$ and dot temperature $T_{M}$.
The preceding section established that $I_{ph}^{Q}$ vanishes at $\lambda=0$.
As we turn on a finite $\lambda$, the magnitude of $I_{ph}^{Q}$ increases with
voltage and $\Delta T_{el}$. Both Case (a) and (b) observe this via Fig. 2 and
2 respectively. However, we find that the direction of $I_{ph}^{Q}$ is
opposite in these cases. Figure 2 notes that in Case (a), the dot phonons
relax to the reservoir and the magnitude of $I_{Ph}^{Q}$ increases in the
negative direction. On the other way, Fig. 2 shows that in Case (b), the dot
becomes phonon deficient and the reservoir pumps phonons into the dot. This
clearly indicates that $\lambda$ acts like a phonon switch whereas voltage and
$\Delta T_{el}$ control the direction of $I_{Ph}^{Q}$. These results can be
utilized in the design of thermal transistors and rectifiers which aim to
implement digital logic by modulating the heat current.

Figure 3: Study of electron induced phonon generation and absorption at finite
electron-phonon interaction. (a) Variation of $G_{ph}$ with voltage for
different $\lambda$ in Case (a). (b) Positive deviation of $\langle
N_{ph}^{eq}\rangle$ from $\langle N_{ph}^{eq}\rangle$ as voltage increases for
finite $\lambda$ in Case (a). The magnitude of phonon current increases with
voltage when $\lambda$ is non-zero (c) Variation of $G_{ph}$ as a function of
$\Delta T_{el}$ for different $\lambda$ in Case (b). (d) Negative deviation of
$\langle N_{ph}^{eq}\rangle$ from the level of $\langle N_{ph}^{eq}\rangle$ as
$\Delta T_{el}$ increases for non-zero $\lambda$ in Case (b). We infer that in
Case (a) phonons accumulate in the dot and in Case (b) phonons get absorbed
and explain opposite polarity of phonon current in the two cases.

The anomaly in the direction of $I_{Ph}^{Q}$ can be explained by studying the
phonon generation rate ($G_{ph}$) in Case (a) and (b). Figure 3 notes that in
Case (a), $G_{ph}$ is positive and rises with voltage when $\lambda$ is
finite. Hence, phonons accumulate in the dot as $\langle N_{ph}\rangle$
exceeds $\langle N_{ph}^{eq}\rangle$. The voltage response of $\langle
N_{ph}\rangle$ is presented in Fig. 3. In contrast, Fig. 3 observes that in
Case (b), $G_{ph}$ becomes negative and falls with $\Delta T_{el}$ when
$\lambda$ is non-zero. Hence the rise of $\Delta T_{el}$ facilitates phonon
absorption and $\langle N_{ph}\rangle$ falls below the level of $\langle
N_{ph}^{eq}\rangle$ as depicted in Fig. 3. As a result in Case(b), reservoirs
pump phonons into the dot.

Figure 4:  Linear response of $I_{Q}^{ph}$. (a) Variation of $I_{Q}^{ph}$ with
voltage in the linear regime for different $\lambda$ in Case (a). (b)
Variation of $I_{Q}^{ph}$ with $\Delta T_{el}$ in the linear regime for
different $\lambda$ in Case (b). (c) Variation of $P_{n,q}$ with linear
voltage in Case (a) at non-zero $\lambda$ ($\lambda=0.5$). (d) Variation of
$P_{n,q}$ with linear range of electronic temperature bias in Case (b) at non-
zero $\lambda$ ($\lambda=0.5$). This figure notes that in the linear regime
the probabilities of the states $|n,q\rangle$ remain constant and consequently
$I_{Q}^{ph}$ show no variation with voltage or $\Delta T_{el}$.

One should note that the variation of $I_{ph}^{Q}$ with voltage and $\Delta
T_{el}$ takes place in a non-linear fashion. Now it is essential to test the
features of $I_{ph}^{Q}$ in the linear regime. From the fundamental
standpoint, this is imperative since Onsager’s reciprocity should be validated
in the linear regime. Figure 4 depicts that in Case (a), the phonon heat
current does not vary with voltage in the linear limit (q$V<<k_{B}T$).
Similarly, Fig.4 shows that in Case (b), $I_{Ph}^{Q}$ does not change with
$\Delta T_{el}$ in the linear range ($\Delta T_{el}<<T$). To explore the
linear response of $I_{Q}^{ph}$, we test the nature of many body electron-
phonon probabilities $P_{n,q}$ of the states $|n,q\rangle$ in Case (a) and
(b), when $\lambda$ is non-zero. Figure.4 presents that in Case (a), $P_{n,q}$
does not vary with voltage. Similarly, Fig.4 depicts that in Case (b),
$P_{n,q}$ remains constant with $\Delta T_{el}$. Hence, in the linear regime,
the average population of phonons in the dot ($\langle
N_{ph}\rangle=\sum\limits_{n,q}qP_{n,q}$) remains unaltered and $I_{Ph}^{Q}$
remains constant as dictated in (15). This implies that the electron and
phonon transport are uncoupled in the linear regime .

Figure 5: Test of the inter-dependence of $I$ and $\Delta T_{ph}$. Color
variation of $I$ as a function of $\Delta T_{ph}$ and $V$ for (a) $\lambda=0$
and (b) $\lambda=0.5$ when $\Delta T_{el}=0$. Color variation of $I$ as a
function of $\Delta T_{ph}$ and $\Delta T_{el}$ for (c) $\lambda=0$ and (d)
$\lambda=0.5$ at zero voltage. It implies that in our setup $I$ can not be
stimulated by $\Delta T_{ph}$ even in the non-linear regime.

After the discussion of linear response of $I_{Q}^{ph}$, it is essential to
examine the inter-relation of electron current $I$ and the phonon temperature
gradient ($\Delta T_{ph}=T_{H}-T_{C}$). In this course, first we set $\Delta
T_{el}$ at zero and test the variation of $I$ as a function of voltage and
$\Delta T_{ph}$. Figure.5 and 5 depict the color variation of $I$ as a
function of $V$ and $\Delta T_{ph}$ for $\lambda=0$ and $\lambda=0.5$
respectively. It is observed that even in the non-linear regime, $I$ shows no
dependence on $\Delta T_{ph}$. Next, we study the variation of $I$ with
$\Delta T_{el}$ and $\Delta T_{ph}$ at zero voltage. The color variation of
$I$ as a function of $\Delta T_{el}$ and $\Delta T_{ph}$ is depicted in Fig.5
and 5 for $\lambda=0$ and $\lambda=0.5$ respectively. Here also, we notice
that $I$ shows no dependence on $\Delta T_{ph}$. This result is consistent
with the previous reportLü et al. (2016) which demonstrated why $I$ can not be
stimulated by $\Delta T_{ph}$ in quantum dot based devices with a single
phonon mode. In the linear regime, several charge and heat currents
$I$,$I_{Q}^{el}$ and $I_{Q}^{ph}$ are related to voltage and various
temperature gradient $V$, $\Delta T_{el}$ and $\Delta T_{ph}$ by the Onsager’s
matrix, such that

$\begin{bmatrix}{\frac{{I}}{{q}}}\\\ \\\ {\frac{{I_{Q}^{el}}}{{k_{B}T}}}\\\
\\\
{\frac{{I_{Q}^{ph}}}{{k_{B}T}}}\end{bmatrix}=\begin{bmatrix}L_{11}&L_{12}&L_{13}\\\
\\\ L_{21}&L_{22}&L_{23}\\\ \\\
L_{31}&L_{32}&L_{33}\end{bmatrix}\begin{bmatrix}qV\\\ \\\ k_{B}\Delta
T_{el}\\\ \\\ k_{B}\Delta T_{ph}\end{bmatrix}.$

In the present scenario, we find that $L_{13}=L_{31}=0$ and $L_{23}=L_{32}=0$.
It confirms that Onsager’s reciprocity is obeyed. Now the dependence of
$I_{Q}^{ph}$ on the voltage and $\Delta T_{el}$, immediately hints at the
concept of a novel thermal conductivity coefficient which is different fron
the conventional electronic thermal conductivity. We will elaborate on it in
the next subsection.

Figure 6: Phonon Peltier and electron-assisted phonon thermal conductivity
coefficients. (a) Color variation of electronic Peltier coefficient $\pi_{el}$
as a function of $\lambda$ and $V$. (b) Variation of the phonon Peltier
coefficient as a function of voltage for different $\lambda$. (c) Color
variation of electronic thermal conductivity $\kappa_{el}$ as a function of
$\lambda$ and $\Delta T_{el}$. (d) Variation of $\kappa_{el-I_{p}h}$ as a
function of $\Delta T_{el}$ for different $\lambda$. We infer from this figure
that $\pi_{el}$ and $\kappa_{el}$ are independent of $\lambda$ although
$\pi_{ph}$ and $\kappa_{el-I_{p}h}$ are strong functions of $\lambda$.

###  III.2 Non-linear phonon Peltier and electron-induced phonon thermal
conductivity

The conventional electronic Peltier coefficient is defined as
$\pi_{el}=\frac{dI_{Q}^{el}}{dI}$, when there is no electronic temperature
gradient ($\Delta T_{el}=0$). In analogy with the electronic Peltier
coefficient, our previous work De and Muralidharan (2018) has introduced the
phonon Peltier coefficient as $\pi_{ph}=\frac{dI_{Q}^{ph}}{dI}$. In this
subsection, we draw a comparative study between the two Peltier co-efficients.
Both Peltier coefficients are computed in the differential form to avoid the
singularity at the short-circuit point ($V=0$). First, we plot the color
variation of $\pi_{el}$ as a function of voltage and $\lambda$ in Fig.6. It
shows that $\pi_{el}$ varies with voltage but remains completely independent
of $\lambda$. It implies that in the absence of an electronic temperature
gradient, the ratio of the electronic heat current $I_{Q}^{el}$ and the charge
current does not vary $\lambda$. In contrast, Fig.6 notes that $\pi_{ph}$ is a
joint function of voltage and $\lambda$. $\pi_{ph}$ vanishes when $\lambda$ is
zero and it varies linearly with the voltage as we switch on a finite
$\lambda$. Earlier Fig.2 showed that the magnitude of $I_{Q}^{ph}$ increases
with $\lambda$ and Fig.5 showed that the magnitude of $I$ falls with $\lambda$
at fixed voltage. Therefore, the absolute magnitude of $\pi_{ph}$ increases
with $\lambda$ as shown in Fig.6. Since $I_{Q}^{ph}$ varies with voltage only
in the non-linear regime, we mention the phonon Peltier co-efficient
($\pi_{ph}$) as non-linear phonon Peltier coefficient.
The concept of $\pi_{ph}$ was manifested since voltage bias stimulates and
modulates $I_{Q}^{ph}$, when $\Delta T_{el}$ is zero. We have already observed
that, $I_{Q}^{ph}$ is a joint function of $V$ and $\Delta T_{el}$. This
motivates us to propose a new thermal conductivity coefficient $\kappa_{el-I-
ph}$, which is different from the conventional electronic thermal conductivity
$\kappa_{el}$. Both thermal conductivity coefficients $\kappa_{el}$ and
$\kappa_{el-I-ph}$ are defined as,

$\left.\kappa_{el}=-\frac{\partial I_{Q}^{el}}{\partial\Delta
T_{el}}\right|_{\mathrlap{I=0}},$ $\left.\kappa_{el-I-ph}=-\frac{\partial
I_{Q}^{ph}}{\partial\Delta T_{el}}\right|_{\mathrlap{I=0}}.$

The color plot shown in Fig.6, depicts that $\kappa_{el}$ remains constant
with $\lambda$ and shows minor variation with $\Delta T_{el}$. On the other
hand, Fig.6 notes that the absolute magnitude of $\kappa_{el-I-ph}$ increases
with both $\lambda$ and $\Delta T_{el}$. At constant $\lambda$, the open-
circuit voltage of the device increases with the increase of $\Delta T_{el}$.
On the other hand, the magnitude of $I_{Q}^{ph}$ increases with the voltage
bias as shown in Fig.2. Therefore, the absolute magnitude of $\kappa_{el-I-
ph}$ rises with $\Delta T_{el}$ as shown in Fig.6. Therefore, the absolute
magnitude of $\kappa_{el-I-ph}$ rises with $\Delta T_{el}$. Since $I_{Q}^{ph}$
varies with voltage and $\Delta T_{el}$ only in the non-linear regime, we
mention $\kappa_{el-I-ph}$, as the non-linear electron assisted phonon thermal
conductivity.

Figure 7: Control of $I_{Q}^{ph}$ by voltage bias in Case (a), when $\Delta
T_{el}=0$. Two situations are considered: Case (HR) when the reservoir is
hotter than the electronic contacts ($T_{L}=T_{R}=100$ K, $T_{B}=T_{BH}=150$
K) and Case (CR) when the reservoir is colder than the electronic contacts
($T_{L}=T_{R}=150$ K, $T_{B}=T_{BC}=100$ K). (a) The modified device diagrams
in the two cases. (b) Voltage Variation of $I_{Q}^{ph}$ for $\lambda=0$ and
$\lambda\neq 0$ in Case (HR). (c) Voltage variation of $\langle N_{ph}\rangle$
in Case (HR). Inset shows the variation of dot temperature $T_{M}$ with
voltage. (d) Voltage variation of $I_{Q}^{ph}$ for $\lambda=0$ and
$\lambda\neq 0$ in Case (CR). (e) Voltage variation of $\langle N_{ph}\rangle$
with voltage in Case (CR). Inset shows the variation of $T_{M}$ with voltage.
(f) The schematic of $I_{Q}^{ph}$ for Case (HR). This figure shows that in
Case (a), the polarity of $I_{Q}^{PH}$ can be reversed by driving the voltage
in Case (HR).

###  III.3 Engineering of $I_{Q}^{ph}$ and $I_{Q}^{el}$ by varying reservoir
temperature

It is evident from (7),(8) and (12) that $I_{Q}^{ph}$ can be controlled by
differing the temperature ($T_{B}$) of the reservoir from the temperature of
the electronic contacts ($T_{L(R)}$). In this subsection, we demonstrate the
engineering of $I_{Q}^{ph}$ in two situations: Case (HR) where, the reservoir
is hotter than the contacts ($T_{B}=T_{BH}>T_{L(R)}$) and Case (CR) where, the
reservoir is colder than the contacts ($T_{B}=T_{BC}<T_{L(R)}$). First, we
illustrate these cases in the limit of Case (a) where $I_{Q}^{ph}$ is
controlled by driving the voltage bias at zero electronic temperature gradient
($\Delta T_{el}=0$) .

Figure 8: Control of $I_{Q}^{ph}$ by driving $\Delta T_{el}$, at $V=0$. Two
cases are considered: Case (HR) Case (HR) when the reservoir is hotter than
the electronic contacts ($T_{H}=T_{C}=T_{BH}=300K>T_{L}=T_{R}=100K$) and Case
(CR) Cold Bath: when the reservoir is colder than the electronic contacts
($T_{H}=T_{C}=T_{BC}=250K<T_{L}=T_{R}=350K$). (a) The modified device
schematics in two cases. (b) Variation of $I_{Q}^{ph}$ as a function of
$\Delta T_{el}$ for $\lambda=0$ and $\lambda\neq 0$ for Case (HR). (c)
Variation of $\langle N_{ph}\rangle$ with $\Delta T_{el}$ in Case (HR). Inset
shows the variation of $T_{M}$ with $\Delta T_{el}$ in the same case. (d)
Variation of $I_{Q}^{ph}$ with $\Delta T_{el}$ for $\lambda=0$ and
$\lambda\neq 0$. (e) Variation of $\langle N_{ph}\rangle$ with $\Delta
T_{el}$. Inset shows the variation of $T_{M}$ with $\Delta T_{el}$ in the same
case. (f) The schematic of $I_{Q}^{ph}$ in Case (CR). This figure infers that
the schematic of that in Case (b), the polarity odf $I_{Q}^{ph}$ can be
flipped by driving $\Delta T_{el}$ in Case (CR).

We present the setups in Case (HR) and (CR) in Fig.7. Figure 7 plots the
voltage response of $I_{Q}^{ph}$ in Case (HR). In this case, $I_{Q}^{ph}$
flows from the dot to the reservoir in the low voltage range but reverses it’s
polarity at $|V|=V_{0}$, when $\lambda$ is non-zero. This non-trivial
phenomenon can be explained by studying the voltage variation of $\langle
N_{ph}^{eq}\rangle$ depicted in Fig.7. As we increase the voltage, phonons
accumulate in the dot and at $V=|V_{0}|$, the average phonon number in the dot
exceeds the equilibrium phonon number $\langle N_{ph}^{eq}\rangle$ of the
reservoir. Hence, the polarity of $I_{Q}^{ph}$ is flipped. However, Fig.7
captures no such non-trivial polarity reversal of $I_{Q}^{ph}$ in Case (CR).
In this case, Fig.7 notes that $\langle N_{ph}\rangle$ is always higher than
the mark of $\langle N_{ph}^{eq}\rangle$ and $I_{Q}^{ph}$ is directed from the
dot to the cold reservoir. The voltage variation of $\langle N_{ph}\rangle$
presented in Fig.7 and 7 is consistent with the voltage variation of the dot
temperature $T_{M}$ shown in their respective insets. It is interesting to
observe that the electron current emanating from electronic contacts can pump
phonons into the hotter reservoir. This is a type of counter-intuitive heating
of a hot reservoir. Figure 7 depicts the schematic of counter-intuitive
heating in Case (HR).

Figure 9: Study of small voltage response of electronic heat current
$I_{QL(R)}^{el}$ associated with contact $L(R)$, when $\Delta T_{el}=0$. We
consider two cases:(1) Hot Bath: $T_{H}=T_{C}>T_{L}=T_{R}$ and (2) Cold Bath:
$T_{H}=T_{C}<T_{L}=T_{R}$. (a) and (b) plot the respective voltage variation
of $I_{QL}^{el}$ and $I_{QR}^{el}$ for Hot Bath as $\lambda$ varies. Similarly
(c) and (d) show the respective voltage variation of $I_{QL}^{el}$ and
$I_{QR}^{el}$ for Cold Bath as $\lambda$ varies. (e) Present the schematic of
$I_{QL(R)}^{el}$ in two cases. (f) and(g) depict the voltage dependence of
$I_{QL}^{el},I_{QR}^{el}$ for non-zero $\lambda$ for cases (1) and (2)
respectively, when $\gamma_{L}>>\gamma_{R}$. (h) presents the schematic for
two cases given $\gamma_{L}>>\gamma_{R}$.

An opposite response of $I_{Q}^{ph}$ is observed when we analyze Case (HR) and
(CR) in the limit of Case (b) where we drive the electronic temperature
gradient at zero voltage. Figure 8 depicts the modified device schematics in
Case (HR) and (CR). Figure.8 plots the variation of $I_{Q}^{ph}$ with $\Delta
T_{el}$ in the Case (HR) for $\lambda=0$ and $\lambda\neq 0$. In this case, we
note no non-triviality in the direction of $I_{Q}^{ph}$ and reservoir pumps
phonons into the dot for all values of $\Delta T_{el}$. Fig.8 shows that in
this case, $\langle N_{ph}\rangle$ is always lower than the level of $\langle
N_{ph}^{eq}\rangle$ and $I_{Q}^{ph}$ does not flip its polarity. In contrast,
Fig.8 depicts that in Case (CR), $I_{Q}^{ph}$ reverses it’s direction at
$\Delta T_{el}=|T_{0}|$ for non-zero $\lambda$ and cold reservoir pumps
phonons into the dot. As explained earlier, the increase of $\Delta T_{el}$
enhances the absorption of phonons and the average phonon number in the dot
falls below the level of at $\Delta T_{el}=|T_{0}|$. This phenomenon is
explained with the help of Fig.8.The $\Delta T_{el}$ variation of $\langle
N_{ph}\rangle$ presented in Fig.8 and 8 is consistent with the $\Delta T_{el}$
variation of the dot temperature $T_{M}$ shown in their respective insets.
Therefore, this case observes a counter-intuitive cooling is of a cold
reservoir. Figure 8 depicts the schematic of counter-intuitive cooling in the
Case (CR).
Finally, we elaborate on the engineering of the electronic heat current
($I_{Q}^{el}$) in Case (HR) and (CR). In this context we rename $I_{Q}^{el}$
associated with the contacts $L$ and $R$ as $I_{QL}^{el}$ and $I_{QR}^{el}$
respectively. Figure 9 and 9 depict the voltage variation of $I_{QL}^{el}$ and
$I_{QR}^{el}$ for different values of $\lambda$ in Case (HR). We observe that
both $I_{QL}^{el}$ and $I_{QR}^{el}$ vanishes when $\lambda=0$. When $\lambda$
is non-zero, both $I_{QL}^{el}$ and $I_{QR}^{el}$ becomes finite and they
disperse uniformly into the contacts even at the short-circuit point ($V=0$).
At $\lambda\neq 0$, the Hot Bath push phonons into the dot and the excess heat
flows away into the contacts in the form of $I_{Q}^{el}$. Similarly Fig.9 and
9 explain the same phenomenon in Case (CR). In this case, the Cold Bath
extracts phonons out of the dot and $I_{Q}^{el}$ flows from the contacts to
the dot. The whole scenario is pictorially presented in Fig.9. It is
intriguing to notice that $I_{Q}^{el}$ can be stimulated by varying the
temperature of phonon reservoir even at the short circuit point $V=0$, where
the effective electron flow vanishes (i.e $I=0$). Now if the contact $L$ is
strongly coupled to the dot as compared to the contact $R$, (i.e
$\gamma_{L}>>\gamma_{R}$), $I_{Q}^{el}$ is dragged out from (or pushed into)
the dot by the contact $L$. Fig.9 and 9 plot the small voltage response of
$I_{QL}^{el}$ and $I_{QR}^{el}$ for non-zero $\lambda$ in Case (HR) and (CR)
respectively. They show that the hotter (colder) reservoir can selectively
heat up (or cool down) $L$ even when $V=0$ as shown in the schematic presented
in Fig.9. Therefore, the temperature of phonon reservoir plays a major role in
the stimulation and control of the electronic and phonon heat current in the
dot. We believe that this intriguing physics can be implemented in the design
of thermal transistor or thermal logic devices.

##  IV  Conclusion

The anomalous behavior of phonon transport due to finite electron-phonon
interaction was investigated using an Anderson-Holstein based dissipative
quantum dot setup in two relevant cases: (a) electron flow stimulated by a
voltage bias in the absence of an electronic temperature gradient and (b)
electron flow driven by the electronic temperature gradient at zero voltage.
We explained the observed cumulative effects of voltage and electronic
temperature gradients on the non-linear phonon currents, using a new transport
coefficient termed as electron induced phonon thermal conductivity. It was
demonstrated that under suitable operating conditions in Case (a) the dot
pumped in phonons into the hotter phonon reservoirs and in Case (b) the dot
extracted phonons out of the colder phonon reservoirs. Finally, we elaborated
how the non-linear electronic heat current was stimulated and controlled by
engineering the temperature of the phonon reservoirs coupled to the dot.

## References

  * Giustino (2017) F. Giustino, Rev. Mod. Phys. 89, 015003 (2017).
  * Lanzara et al. (2001) A. Lanzara, P. V. Bogdanov, X. J. Zhou, S. A. Kellar, D. L. Feng, E. D. Lu, T. Yoshida, H. Eisaki, A. Fujimori, K. Kishio, et al., Nature 412, 510 (2001).
  * Walborn et al. (2006) S. P. Walborn, P. H. Souto Ribeiro, L. Davidovich, F. Mintert, and A. Buchleitner, Nature 440, 1170 (2006).
  * Lee et al. (2006) J. Lee, K. Fujita, K. McElroy, J. A. Slezak, Y. Wang, M andAiura, H. Bando, M. Ishikado, T. Masui, J.-X. Zhu, A. V. Balatsky, et al., Nature 442, 546 (2006).
  * Hong et al. (2009) I.-P. Hong, C. Brun, F. m. c. Patthey, I. Y. Sklyadneva, X. Zubizarreta, R. Heid, V. M. Silkin, P. M. Echenique, K. P. Bohnen, E. V. Chulkov, et al., Phys. Rev. B 80 (2009).
  * Bianchi et al. (2010) M. Bianchi, E. D. L. Rienks, S. Lizzit, A. Baraldi, R. Balog, L. Hornekær, and P. Hofmann, Phys. Rev. B 81, 041403 (2010).
  * Hartke et al. (2018) T. R. Hartke, Y.-Y. Liu, M. J. Gullans, and J. R. Petta, Phys. Rev. Lett. 120, 097701 (2018).
  * Frey et al. (2012) T. Frey, P. J. Leek, M. Beck, A. Blais, T. Ihn, K. Ensslin, and A. Wallraff, Phys. Rev. Lett. 108, 046807 (2012).
  * Ares et al. (2016) N. Ares, T. Pei, A. Mavalankar, M. Mergenthaler, J. H. Warner, G. A. D. Briggs, and E. A. Laird, Phys. Rev. Lett. 117, 170801 (2016).
  * Park et al. (2004) J.-Y. Park, S. Rosenblatt, Y. Yaish, V. Sazonova, H. Üstünel, S. Braig, T. A. Arias, P. W. Brouwer, and P. L. McEuen, Nano Letters 4, 193 (2004).
  * Laird et al. (2012) E. A. Laird, F. Pei, W. Tang, G. A. Steele, and L. P. Kouwenhoven, Nano Letters 12, 517 (2012).
  * Javey et al. (2003) A. Javey, J. Guo, Q. Wang, M. Lundstrom, and H. Dai, Nature 424, 654 (2003).
  * Yu et al. (2004a) L. H. Yu, Z. K. Keane, J. W. Ciszek, L. Cheng, M. P. Stewart, J. M. Tour, and D. Natelson, Phys. Rev. Lett. 93, 266802 (2004a).
  * Rakhmilevitch et al. (2014) D. Rakhmilevitch, R. Korytár, A. Bagrets, F. Evers, and O. Tal, Phys. Rev. Lett. 113, 236603 (2014).
  * Ward et al. (2010) D. R. Ward, D. A. Corley, J. M. Tour, and D. Natelson, Nature Nanotechnology 6, 33 (2010).
  * Ioffe et al. (2008) Z. Ioffe, T. Shamai, A. Ophir, G. Noy, I. Yutsis, K. Kfir, O. Cheshnovsky, and Y. Selzer, Nature Nanotechnology 3, 727 (2008).
  * Benyamini et al. (2014) A. Benyamini, A.Hamo, S. Kusminskiy, F. V. Oppen, and S. Ilani, Nat .Phys. 10, 151 (2014).
  * Härtle and Thoss (2011) R. Härtle and M. Thoss, Phys. Rev. B 83, 115414 (2011).
  * Leijnse et al. (2010) M. Leijnse, M. R. Wegewijs, and K. Flensberg, Phys. Rev. B 82, 045412 (2010).
  * Galperin et al. (2004) M. Galperin, M. A. Ratner, and A. Nitzan, J. Chem. Phys. 121, 11965 (2004).
  * Galperin et al. (2007) M. Galperin, A. Nitzan, and M. A. Ratner, Phys. Rev. B 75, 155312 (2007).
  * De and Muralidharan (2016) B. De and B. Muralidharan, Phys. Rev. B 94, 165416 (2016).
  * Agarwalla et al. (2015) B. K. Agarwalla, J.-H. Jiang, and D. Segal, Phys. Rev. B 92, 245418 (2015).
  * Arrachea et al. (2014) L. Arrachea, N. Bode, and F. von Oppen, Phys. Rev. B 90, 125450 (2014).
  * Sothmann et al. (2014) B. Sothmann, R. Sánchez, and A. N. Jordan, Nanotechnology 26, 32001 (2014).
  * Siddiqui et al. (2006) L. Siddiqui, A. W. Ghosh, and S. Datta, Phys. Rev. B p. 085433 (2006).
  * Mitra et al. (2004) a. Mitra, I. Aleiner, and a. J. Millis, Phys. Rev. B 69, 245302 (2004).
  * Kim et al. (2014) Y. Kim, W. Jeong, K. Kim, W. Lee, and P. Reddy, Nat. Nanotechnol. 9, 881 (2014).
  * Reddy et al. (2007) P. Reddy, S.-Y. Jang, R. A. Segalman, and A. Majumdar, Science 315, 1568 (2007).
  * Roßnagel et al. (2016) J. Roßnagel, S. T. Dawkins, K. N. Tolazzi, O. Abah, E. Lutz, F. Schmidt-Kaler, and K. Singer, Science 352, 325 (2016).
  * Steele et al. (2009) G. A. Steele, A. K. Hüttel, B. Witkamp, M. Poot, H. B. Meerwaldt, L. P. Kouwenhoven, and H. S. J. van der Zant, Science 325, 1103 (2009).
  * Lassagne et al. (2009) B. Lassagne, Y. Tarakanov, J. Kinaret, D. Garcia-Sanchez, and A. Bachtold, Science 325, 1107 (2009).
  * Woodside and McEuen (2002) M. T. Woodside and P. L. McEuen, Science 296, 1098 (2002).
  * Jiang et al. (2015) J.-H. Jiang, M. Kulkarni, D. Segal, and Y. Imry, Phys. Rev. B 92, 045309 (2015).
  * Wang and Li (2007) L. Wang and B. Li, Phys. Rev. Lett. 99, 177208 (2007).
  * Wang and Li (2008) L. Wang and B. Li, Phys. Rev. Lett. 101, 267203 (2008).
  * Li et al. (2012) N. Li, J. Ren, L. Wang, G. Zhang, P. Hänggi, and B. Li, Rev. Mod. Phys. 84, 1045 (2012).
  * Joulain et al. (2016) K. Joulain, J. Drevillon, Y. Ezzahri, and J. Ordonez-Miranda, Phys. Rev. Lett. 116, 200601 (2016).
  * Mahan and Sofo (1996) G. D. Mahan and J. O. Sofo, Proceedings of the National Academy of Sciences of the United States of America 93, 7436 (1996).
  * Kubala and König (2006) B. Kubala and J. König, Phys. Rev. B 73, 195316 (2006).
  * Kubala et al. (2008) B. Kubala, J. König, and J. Pekola, Phys. Rev. Lett 100, 066801 (2008).
  * Park et al. (2000) H. Park, J. Park, A. Lim, E. Anderson, A. Alivisatos, and P. McEuen, Nature 407, 57 (2000).
  * Leroy et al. (2004) B. J. Leroy, S. G. Lemay, J. Kong, and C. Dekker, Nature 432, 371 (2004).
  * LeRoy et al. (2005) B. J. LeRoy, J. Kong, V. K. Pahilwani, C. Dekker, and S. G. Lemay, Phys. Rev. B 72, 075413 (2005).
  * Yu et al. (2004b) L. H. Yu, Z. K. Keane, J. W. Ciszek, L. Cheng, M. P. Stewart, J. M. Tour, and D. Natelson, Phys. Rev. Lett 93, 266802 (2004b).
  * Sapmaz et al. (2006) S. Sapmaz, P. Jarillo-Herrero, Y. M. Blanter, C. Dekker, and H. S. J. Van Der Zant, Phys. Rev. Lett 96, 026801 (2006).
  * Schmidt and Komnik (2009) T. L. Schmidt and A. Komnik, Phys. Rev. B 80, 041307 (2009).
  * Urban et al. (2010) D. F. Urban, R. Avriller, and A. Levy Yeyati, Phys. Rev. B 82, 121414 (2010).
  * Egger and Gogolin (2008) R. Egger and A. O. Gogolin, Phys. Rev. B 77, 113405 (2008).
  * Seoane Souto et al. (2014) R. Seoane Souto, A. L. Yeyati, A. Martín-Rodero, and R. C. Monreal, Phys. Rev. B 89, 085412 (2014).
  * Ness (2014) H. Ness, Phys. Rev. B 89, 045409 (2014).
  * Wilner et al. (2014) E. Y. Wilner, H. Wang, M. Thoss, and E. Rabani, Phys. Rev. B 89, 205129 (2014).
  * De and Muralidharan (2018) B. De and B. Muralidharan, Scientific Reports 8, 5185 (2018).
  * Braig and Flensberg (2003) S. Braig and K. Flensberg, Phys. Rev. B 68, 205324 (2003).
  * Segal and Nitzan (2005a) D. Segal and A. Nitzan, J. Chem. Phys. 122 (2005a).
  * Segal (2006) D. Segal, Phys. Rev. B 73, 205415 (2006).
  * Segal and Nitzan (2005b) D. Segal and A. Nitzan, Phys. Rev. Lett 94, 034301 (2005b).
  * Muralidharan et al. (2006) B. Muralidharan, A. W. Ghosh, and S. Datta, Phys. Rev. B 73, 155410 (2006).
  * Muralidharan and Datta (2007) B. Muralidharan and S. Datta, Phys. Rev. B 76, 035432 (2007).
  * Beenakker (1991) C. W. J. Beenakker, Phys. Rev. B 44, 1646 (1991).
  * Muralidharan and Grifoni (2013) B. Muralidharan and M. Grifoni, Phys. Rev. B 88, 045402 (2013).
  * Muralidharan and Grifoni (2012) B. Muralidharan and M. Grifoni, Phys. Rev. B 85, 155423 (2012).
  * Timm (2008) C. Timm, Phys. Rev. B 77, 195416 (2008).
  * Piovano (2012) G. Piovano, Effects of Electron-Vibron Coupling in Nano-electromechanical Systems Phd Thesis (2012).
  * Braig and Brouwer (2005) S. Braig and P. W. Brouwer, Phys. Rev. B 71, 195324 (2005).
  * König and Martinek (2003) J. König and J. Martinek, Phys. Rev. Lett. 90, 166602 (2003).
  * Braun et al. (2004) M. Braun, J. König, and J. Martinek, Phys. Rev. B 70, 195345 (2004).
  * Koch and Von Oppen (2005) J. Koch and F. Von Oppen, Physical Review Letters 94, 206804 (2005).
  * Koch et al. (2006) J. Koch, M. Semmelhack, F. Von Oppen, and A. Nitzan, Phys. Rev. B 73, 155306 (2006).
  * Koch (2006) J. Koch, Quantum transport through single-molecule devices Phd Thesis (2006).
  * Lü et al. (2016) J.-T. Lü, J.-S. Wang, P. Hedegård, and M. Brandbyge, Phys. Rev. B 93, 205404 (2016).