# Photonic heat transport through a Josephson junction in a resistive
environment

A. Levy Yeyati  Departamento de Física Teórica de la Materia Condensada,
Condensed Matter Physics Center (IFIMAC), and Instituto Nicolás Cabrera,
Universidad Autónoma de Madrid, 28049 Madrid, Spain  D. Subero  PICO Group,
QTF Centre of Excellence, Department of Applied Physics, Aalto University
School of Science, P.O. Box 13500, 0076 Aalto, Finland  J. P. Pekola  PICO
Group, QTF Centre of Excellence, Department of Applied Physics, Aalto
University School of Science, P.O. Box 13500, 0076 Aalto, Finland  R. Sánchez
Departamento de Física Teórica de la Materia Condensada, Condensed Matter
Physics Center (IFIMAC), and Instituto Nicolás Cabrera, Universidad Autónoma
de Madrid, 28049 Madrid, Spain

###### Abstract

Motivated by recent experiments (Subero et. al. Nature Comm. 14, 7924 (2023))
we analyze photonic heat transport through a Josephson junction in a
dissipative environment. For this purpose we derive the general expressions
for the heat current in terms of non-equilibrium Green functions for the
junction coupled in series or in parallel with two environmental impedances at
different temperatures. We show that even on the insulating side of the Schmid
transition the heat current is sensitive to the Josephson coupling exhibiting
an opposite behavior for the series and parallel connection and in qualitative
agreement with experiments. We also predict that this device should exhibit
heat rectification properties and provide simple expressions to account for
them in terms of the system parameters.

Introduction.— The physics of Josephson junctions (JJs) has attracted great
interest for many decades [1]. Even the apparently simple case of a single JJ
in a resistive environment provides a rich playground to study the interplay
of quantum tunneling and dissipation [2], which continues to be explored until
today [3, 4]. For this system a transition between a superconducting and an
insulating phase was predicted to appear as a function of the environmental
resistance ($R$), regardless of the Josephson coupling energy ($E_{J}$). This
is the so-called Schmid-Bulgadaev (SB) transition [5, 6], which has been the
object of intense theoretical [7, 8, 9, 10, 11] and experimental debate in
recent years [12, 13, *Murani2021, 15].

In spite of this intense activity little is known regarding this system
properties beyond dc charge transport. In this respect, heat transport [16]
can provide additional insights. In a recent experiment, Subero et al. [17]
explored heat transport through a JJ in the supposedly insulating regime
($R>R_{Q}=h/4e^{2}$). Even when charge transport could be well described by
dynamical Coulomb blockade theory, it was found that heat transport was
sensitive to the $E_{J}$ value (tunable through a magnetic flux in a SQUID
configuration) indicating inductive response at high frequencies. The
experimental results were fitted using a phenomenological theory including an
inductive term to the junction effective impedance.

Figure 1: Schematics of the configurations considered, consisting in two
resistors, $R_{L}$ and $R_{R}$, connected in (a) parallel and (b) series to a
Josephson junction (defined by a Josephson, $E_{J}$, and a charging energy
$E_{C}$), and holding different temperatures, $T_{L}$ and $T_{R}$, that lead
to a heat current $J$.

In the present work we analyze this problem from a microscopic perspective.
For this purpose we develop a non-equilibrium Green functions (GFs) approach
which takes into account interaction effects due to finite $E_{J}/E_{c}$
values, where $E_{c}$ is the junction charging energy. We find that even for
$R>R_{Q}$, heat transport is sensitive to the Josephson coupling, in
qualitative agreement with experiments. Moreover, we find that the heat
current response is radically different depending on whether the junction is
connected in parallel or in series to the source and drain leads. We finally
demonstrate that asymmetric devices exhibit heat rectification properties.

In Fig. 1 we schematically represent the two situations to be considered: in
panel 1(a) the JJ is connected in parallel with the left and right resistors
($R_{L}$ and $R_{R}$) which are held at temperatures $T_{L}$ and $T_{R}$
respectively, whereas in panel 1(b) these circuit elements are connected in
series, corresponding to the experimental configuration in Ref. [17]. The aim
of the theory is to determine the resulting heat current $J$ as a function of
the model parameters in these two situations, assuming that temperatures are
low enough to neglect excited quasiparticles in the superconducting leads.

Parallel configuration.— Let us start with the parallel configuration [18].
This is modeled by the following Hamiltonian

$\displaystyle H$ $\displaystyle=$ $\displaystyle
E_{c}N^{2}-E_{J}\cos\varphi-\sum_{j,k}\left[\lambda_{j,k}\left(a_{jk}+a^{\dagger}_{jk}\right)\varphi-\frac{\lambda^{2}_{jk}}{\hbar\omega_{jk}}\varphi^{2}\right]$
(1) $\displaystyle+\sum_{j,k}\omega_{jk}a^{\dagger}_{jk}a_{jk}\;,$

where $j\equiv L,R$ and the dissipative terms satisfy

$\frac{\pi}{2}\sum_{k}\lambda_{jk}^{2}\delta(\omega-\omega_{jk})=\frac{\hbar^{2}\omega}{4e^{2}}\mbox{Re}\left[\frac{1}{Z_{j}(\omega)}\right]\;,$

$Z_{j}(\omega)$ being the leads impedances, which in our case correspond to
pure resistances $R_{j}$. In this configuration only one superconducting phase
variable, $\varphi$, conjugate to the charge which is transferred through the
junction, $N$, plays a role.

Following Ref. [19], the heat current through lead $j$ can then be computed as

$\displaystyle J_{j}$ $\displaystyle=$
$\displaystyle\langle\dot{H}_{j}\rangle=-i\sum_{k}\lambda_{jk}\omega_{Lk}\left[\langle
a_{jk}(t)\varphi(t)\rangle-\mbox{h.c.}\right]$ (2) $\displaystyle=$
$\displaystyle
2\sum_{k}\lambda_{jk}\omega_{jk}\mbox{Re}G^{+-}_{a_{jk},\varphi}(t,t)\;,$

where $G^{\alpha\beta}_{A,B}(t,t^{\prime})=-i\langle T_{\cal
C}B(t^{\prime})A(t)\rangle$ denotes the Keldysh $\alpha\beta$ GF associated
with operators $A,B$. Using the equation of motion method and transforming it
into frequency representation one can express $G^{+-}_{a_{jk},N}$ as

$G^{+-}_{a_{jk},\varphi}(\omega)=\lambda_{jk}\left[D^{r}(\omega)g^{+-}_{a_{jk}}(\omega)+D^{+-}(\omega)g^{a}_{a_{jk}}(\omega)\right]\;,$
(3)

where the $g^{+-,a}_{a_{jk}}$ correspond to the uncoupled leads and, for
simplicity, we denote $D^{\alpha\beta,r,a}\equiv
G^{\alpha\beta,r,a}_{\varphi\varphi}$ (the superscripts $r,a$ indicate the
retarded, advanced components). In this way we can express the heat current as

$J_{j}=2\sum_{k}\lambda_{jk}^{2}\omega_{jk}\int\frac{d\omega}{2\pi}\mbox{Re}\left[D^{r}(\omega)g^{+-}_{a_{jk}}(\omega)+D^{+-}(\omega)g^{a}_{a_{jk}}(\omega)\right]\;.$
(4)

We can further use $g^{+-}_{a_{jk}}(\omega)=-2\pi
in_{j}(\omega_{jk})\delta(\omega-\omega_{jk})$ and
$g^{a}_{a_{jk}}(\omega)=1/(\omega-\omega_{jk}-i0^{+})$, where $n_{j}(\omega)$
denotes the Bose function for the modes in the $j$-lead. Then, replacing in
Eq. (4) and taking the real part we obtain

$\displaystyle
J_{j}{=}{\int}\frac{d\omega}{2\pi}\frac{\hbar^{2}\omega^{2}}{e^{2}}\mbox{Im}\left[D^{r}(\omega)\coth\left(\frac{\omega\beta_{j}}{2}\right){-}\frac{D^{K}(\omega)}{2}\right]\mbox{Re}\left[\frac{1}{Z_{j}(\omega)}\right]\;,$
(5)

where we have introduced the Keldysh GFs $D^{K}=D^{+-}+D^{-+}$ in the
triangular representation [20, *keldysh_diagram_1965] and
$\beta_{j}=1/k_{B}T_{j}$. Further simplification is obtained if, as in our
case, the lead admittances $1/Z_{j}(\omega)\equiv 1/R_{j}$ are frequency
independent. Then, imposing heat current conservation $J_{R}=-J_{L}=J$ one
gets [19]

$\displaystyle J=-\frac{\hbar^{2}}{2\pi e^{2}}\int
d\omega\frac{\omega^{2}\mbox{Im}D^{r}(\omega)}{R_{L}+R_{R}}\left[n_{L}(\omega)-n_{R}(\omega)\right]\;,$
(6)

which is the heat transport analog of Meir-Wingreen formula for mesoscopic
charge transport [22], allowing us to define an effective heat transmission
coefficient for the parallel configuration

$\tau^{\rm eff}_{\rm
parallel}(\omega)=-\frac{4\eta_{L}\eta_{R}\omega\mbox{Im}D^{r}(\omega)}{\eta_{L}+\eta_{R}}\;,$
(7)

where $\eta_{j}=R_{Q}/(2\pi R_{j})\equiv\alpha_{j}/2\pi$. One can further
check that for the perfect transmission case one obtains $J=J_{\rm
max}\equiv{\pi k_{B}^{2}(T_{L}^{2}-T_{R}^{2})/(12\hbar)}$ as it corresponds to
a ballistic channel [23].

Perturbation theory.— The problem is thus mapped into the evaluation of the
Keldysh GFs $D^{r,a,K}(\omega)$. In the present work we focus on the regime
$R_{L,R}\gtrsim R_{Q}$ and $E_{J}\ll E_{c}$ for which perturbation theory in
$E_{J}$ should be valid. To zero-th order in $E_{J}$ the GFs are given by

$\displaystyle\begin{aligned} D_{0}^{r,a}(\omega)&=\frac{1}{m\omega^{2}\pm
i\omega\left(\eta_{L}+\eta_{R}\right)-m\omega_{c}^{2}},\\\
D_{0}^{K}(\omega)&=\frac{-2i\omega\left(\eta_{L}(1+2n_{L}(\omega))+\eta_{R}(1+2n_{R}(\omega)\right)}{\left(m(\omega^{2}-\omega_{c}^{2})\right)^{2}+\left(\omega(\eta_{L}+\eta_{R})\right)^{2}}\;,\end{aligned}$
(8)

where $m=\hbar/E_{C}$ and $\omega_{c}$ is a small frequency cutoff which
warrants convergence in the Fourier transforms. At zero-th order in $E_{J}$
and $\omega_{c}\rightarrow 0$, from (7) we obtain

$\tau^{0}_{\rm
parallel}(\omega)=\frac{4\eta_{L}\eta_{R}}{(m\omega)^{2}+(\eta_{L}+\eta_{R})^{2}}\;,$
(9)

which reaches a maximum value $4\eta_{L}\eta_{R}/(\eta_{L}+\eta_{R})^{2}$ at
zero frequency.

Higher order corrections can be introduced through the self-energies
$\Sigma^{r,a,K}(\omega)$ associated to the $E_{J}$ term in the Hamiltonian,
which would allow us to evaluate the needed GFs as

$\displaystyle D^{r,a}(\omega)$ $\displaystyle=$
$\displaystyle\left[\left(D_{0}^{r,a}(\omega)\right)^{-1}-\Sigma^{r,a}(\omega)\right]^{-1}$
$\displaystyle D^{K}(\omega)$ $\displaystyle=$
$\displaystyle\left[1+D^{r}(\omega)\Sigma^{r}(\omega)\right]D^{K}_{0}(\omega)\left[1+\Sigma^{a}(\omega)D^{a}_{0}(\omega)\right]$
(10) $\displaystyle+D^{r}(\omega)\Sigma^{K}(\omega)D^{a}(\omega).$

To lowest order in $E_{J}$ we have [24]

$\Sigma^{r,a(1)}(\omega)=E_{J}\exp\left[-\frac{i}{4}D_{0}^{K}(t=0)\right]\rightarrow
0\;,$ (11)

since $\mbox{Im}D_{0}^{K}(t=0)\rightarrow\infty$ for $\omega_{c}\rightarrow
0$. One also obtains $\Sigma^{K(1)}(\omega)=0$ as it corresponds to an
effective static potential. To obtain the frequency-dependent corrections it
is then necessary to go to higher order. To second order in $E_{J}$ one
obtains [24]

$\displaystyle\Sigma^{r(2)}(t)$ $\displaystyle=$ $\displaystyle
E_{J}^{2}\sin\frac{D_{0}^{r}(t)}{2}\exp\left[\frac{i}{2}\left(D_{0}^{K}(t)-D^{K}_{0}(0)\right)\right]-B\delta(t)$
$\displaystyle\Sigma^{K(2)}(t)$ $\displaystyle=$ $\displaystyle-
iE_{J}^{2}\cos\frac{D_{0}^{r}(t)}{2}\cos\frac{D_{0}^{a}(t)}{2}\exp\left[\frac{i}{2}\left(D_{0}^{K}(t)-D^{K}_{0}(0)\right)\right]\;,$

where $B$ is a constant such that $\Sigma^{r(2)}(\omega=0)=0$.

From these expressions one can obtain numerical results for the self-energies
[25]. In addition, fully analytical results for the Keldysh self-energies can
be obtained in the limit $(R_{L}+R_{R})/R_{Q}\rightarrow 1$ and
$T_{j}\rightarrow 0$ [25] (see also comments below).

Replacing these results into Eq. (10), we obtain the results in Fig. 2(a) for
the heat transmission coefficient, which exhibits a strong dependence on
$E_{J}$, in spite of the full suppression of the dc Josephson effect that
occurs in the insulating phase [26]. As shown in the inset of Fig. 2(a), the
total heat current in this configuration decreases with increasing $E_{J}$,
and, as illustrated by the dashed lines in this inset, the sensitivity of the
total heat current with $E_{J}$ exhibits a strong dependence on the lead
resistances.

Figure 2: (a) Effective heat transmission coefficient in the parallel
configuration for increasing $E_{J}$ values (from 0 to 0.4) and $T_{L}=0.2$,
$T_{R}=0.1$, $\omega_{c}=2\times 10^{-3}$ (in units of $E_{c}$) and
$\eta_{L,R}=1/4\pi$. The full lines in the insets correspond to the total heat
current $J$, normalized to the maximum heat current for a ballistic channel
$J_{\rm max}$, as a function of $E_{J}$, while the dashed lines illustrate its
behavior with $\eta_{L}=\eta_{R}$ varying between $1/\pi$ to $1/5\pi$. (b)
$J/J_{\rm max}$ for $E_{J}=0.2$ as a function of mean temperature
$\bar{T}=(T_{L}+T_{R})/2$ in the linear response regime for the same
$\eta_{L}=\eta_{R}$ range.

One gets more insight by analyzing the low frequency renormalization of the
heat transmission coefficient $\tau_{\rm parallel}(\omega)$. In fact, at low
frequencies $\Sigma^{r(2)}(\omega)\sim-i\delta\eta\omega+\delta m\omega^{2}$,
where $\delta\eta$ and $\delta m$ provide renormalization of the damping and
mass model parameters. An analytical expression for the temperature and
resistance scaling of these effective parameters is presented in [25]. At low
temperatures ($T_{L,R}\ll\gamma$) and for $\eta_{L,R}\sim 1/2\pi$ we obtain

$\displaystyle\delta\eta$
$\displaystyle\sim\left(\frac{E_{J}}{\gamma}\right)^{2}\left(\frac{\eta_{L}T_{L}+\eta_{R}T_{R}}{m^{2}\gamma^{2}}\right)^{-2}\prod_{j}\left(\frac{T_{j}}{\gamma}\right)^{\eta_{j}/\pi
m^{2}\gamma^{2}}\,,$ (13)

where $\gamma=(\eta_{L}+\eta_{R})/m$. Notice that for $T_{L}=T_{R}=T$ and
$\alpha_{L}=\alpha_{R}=2\alpha$, the SB transition at $\alpha=1$ is encoded in
the scaling law $\sim T^{2/\alpha-2}$, which results for the renormalization
of the damping parameter [27, 24].

In terms of these renormalized parameters $\tau_{\rm parallel}(\omega)$ can be
written as

$\tau^{\rm eff}_{\rm
parallel}(\omega)=\frac{4\eta_{L}\eta_{R}}{\tilde{m}^{2}\omega^{2}+\tilde{\eta}^{2}}\frac{\tilde{\eta}}{\eta_{L}+\eta_{R}}\;,$
(14)

where $\tilde{\eta}=\eta_{L}+\eta_{R}+\delta\eta$ and $\tilde{m}=m-\delta m$.
This expression accounts for the qualitative behavior of the heat transmission
coefficient that can be observed in Fig. 2(a), i.e. a suppression and
broadening of the $\omega\sim 0$ peak for increasing $E_{J}$.

On the other hand, lower temperatures are required to observe signatures of
the SB transition in the heat transport properties. In Fig. 2(b) we show the
heat current $J$ as a function of $\bar{T}=(T_{L}+T_{R})/2$ in the linear
regime ($\delta T=T_{L}-T_{R}\rightarrow 0$) for different values of
$\eta_{L}=\eta_{R}$, ranging from $1/\pi$ to $1/5\pi$, i.e. across the
transition which for the parallel case should occur at $\eta_{L,R}=1/4\pi$.
While for $\bar{T}\gtrsim\gamma$, $J$ decreases with increasing resistance due
to the narrowing of the $\tau^{\rm eff}_{\rm parallel}(\omega)$ zero-frequency
peak, the tendency is reversed at low temperatures ($\bar{T}\ll\gamma$) due to
the divergence of $\delta\eta$ for $\alpha>1$. The transition between these
two opposite behaviors occurs around $\bar{T}=T^{*}\sim 1/(2\pi e^{3/2})$,
where the damping renormalization is weakly dependent on the environment
resistance [25]. This behavior is similar to the non-monotonous dependence of
the mobility found for quantum Brownian motion in a periodic potential [28].

Figure 3: Same as Fig. 2(a) but for the series configuration.

Series configuration.— In the actual setup in Ref. [17] left and right leads
are placed in a loop for the heat current measurements, which actually
corresponds to a series connection, as depicted in Fig. 1(b). One should
describe this situation in terms of two phase variables $\varphi_{1}$ and
$\varphi_{2}$, each one coupled to the left and right reservoirs respectively,
while the charge $N$ through the junction is conjugate to
$\varphi_{1}-\varphi_{2}$. Thus, the bare GFs adopt a matrix form in the
$\varphi_{1,2}$ space and can be obtained from the following equations

$\displaystyle\left(\begin{array}[]{cc}m(\omega^{2}-\omega_{c}^{2})\pm
i\omega\eta_{L}&-m(\omega^{2}-\omega_{c}^{2})\\\
-m(\omega^{2}-\omega_{c}^{2})&m(\omega^{2}-\omega_{c}^{2})\pm
i\omega\eta_{R}\end{array}\right)\hat{D}^{r,a}_{0}$ $\displaystyle=$
$\displaystyle\hat{I},\;$ (17)

and $\hat{D}^{K}_{0}=\hat{D}^{r}_{0}\hat{d}^{K}_{0}\hat{D}^{a}_{0}$, where

$\hat{d}^{K}_{0}=-2i\omega\;\mbox{diag}\left(\eta_{L}\coth(\beta_{L}\omega/2);\eta_{R}\coth(\beta_{R}\omega/2)\right)\;.$

To include the effect of finite $E_{J}$ one should notice that the junction
term in the Hamiltonian is now
$-E_{J}\cos\left(\varphi_{1}-\varphi_{2}\right)$, thus the corresponding self-
energies in Eq. (LABEL:sigK) have to be expressed in terms of
$\tilde{D}_{0}(t)=\mbox{Tr}\left[\hat{D}_{0}(t)(\hat{I}-\sigma_{x})\right]$,
where $\sigma_{x}$ is a Pauli matrix in $\varphi_{1,2}$ space. The resulting
self-energies $\tilde{\Sigma}^{r,K}(\omega)$ exhibit a similar behavior as in
the parallel configuration (see SM [25]).

On the other hand, the Dyson equations for the dressed GFs (10) acquire a
matrix form in the $\varphi_{1,2}$ space and the corresponding self-energies
are obtained from $\tilde{\Sigma}^{r,a,K}$ as
$\hat{\Sigma}^{r,a,K}=\tilde{\Sigma}^{r,a,K}(\hat{I}-\sigma_{x})$.

The heat current in the series connection cannot be expressed in the compact
form (6). It is, in contrast, necessary to reproduce the steps leading to Eq.
(5), in order to derive the corresponding expression for the heat current in
the series connection

$\displaystyle J_{j}$ $\displaystyle=$ $\displaystyle\int
d\omega\omega^{2}\eta_{j}\mbox{Im}\left[D_{l_{j}l_{j}}^{r}(\omega)\left(1{+}2n_{j}(\omega)\right)-\frac{D_{l_{j}l_{j}}^{K}(\omega)}{2}\right],$
(18)

where $l_{j}=1,2$ for $j=L,R$ respectively. This expression can be decomposed
[25] as $J_{j}=(-1)^{l_{j}}J^{\rm el}+J^{\rm in}_{j}$, where $J^{\rm el}=\int
d\omega\hbar\omega\tau^{\rm eff}_{\rm
series}(\omega)\left[n_{L}(\omega)-n_{R}(\omega)\right]$ corresponds to the
elastic contribution with

$\tau^{\rm eff}_{\rm
series}(\omega)=4\eta_{L}\eta_{R}|D^{r}_{12}(\omega)|^{2}\;,$ (19)

while the inelastic contribution is given by

$J^{\rm in}_{j}=\int
d\omega\omega^{2}\eta_{j}\mbox{Im}\left[\hat{D}^{r}\left\\{\left(\hat{\Sigma}^{r}{-}\hat{\Sigma}^{a}\right)(1{+}2n_{j})-\hat{\Sigma}^{K}\right\\}\hat{D}^{a}\right]_{l_{j}l_{j}}\;.$
(20)

In the non-interacting case (i.e. $E_{J}=0$) the inelastic term vanishes and
the heat current can again be written in the usual Landauer form with a
transmission coefficient at $\omega_{c}\rightarrow 0$

$\tau^{0}_{\rm
series}(\omega)=\frac{4m^{2}\omega^{2}\eta_{L}\eta_{R}}{\left((\eta_{L}+\eta_{R})m\omega\right)^{2}+\left(\eta_{L}\eta_{R}\right)^{2}}\;,$
(21)

which exhibits a zero frequency dip and tends to
$4\eta_{L}\eta_{R}/(\eta_{L}+\eta_{R})^{2}$ at large frequencies, as
illustrated in Fig. 3.

The behavior of the effective transmission coefficient for finite $E_{J}$ is
illustrated in Fig. 3 for the same parameters as in Fig. 2(a). We observe that
at finite $E_{J}$ the zero frequency dip is progressively suppressed. On the
other hand, the inelastic term provides an additional contribution to the heat
current [25]. Thus, in contrast to the parallel configuration, there is a
slight increase of the heat current with $E_{J}$, which is confirmed by the
integrated results in the inset of Fig. 3. As can be observed, this increase
is more pronounced as the leads resistances approach $R_{Q}$. This behavior is
in qualitative agreement with the experimental results in Ref. [17].

Rectification properties.— Another interesting property of a JJ in a resistive
environment is its possible behavior as a heat rectifier. The JJ anharmonicity
has been used to define thermal diodes based on weakly coupled qubits [29],
see also Ref. [19]. Though heat rectifiers using JJs have been measured in the
absence of environmental effects (but additionally coupled to a phonon bath)
[30] or in superconducting qubits [31, 32] (see also Ref. [33] for a review),
their rectifying properties in the insulating side of the SB transition have
not been considered. For simplicity we consider here only the parallel
configuration with an asymmetry provided by $R_{L}\neq R_{R}$. We can then
define a “forward" ($J^{+}$) and a “reverse" ($J^{-}$) heat current given by

$J^{\pm}=\int\frac{d\omega}{\pi}\hbar\omega\tau^{\rm
eff(\pm)}(\omega)\left[n^{\pm}_{L}(\omega)-n^{\pm}_{R}(\omega)\right]\;,$ (22)

which correspond to positive and negative temperature bias $\Delta
T=T_{L}-T_{R}$. As is customary, we also define a heat rectification ratio
${\cal R}=|J^{+}/J^{-}|$ which, as can be seen from Eq. (22), would deviate
from unity provided $\delta\tau^{\rm eff}=\tau^{\rm eff(+)}-\tau^{\rm
eff(-)}\neq 0$.

In Fig. 4 we show the behavior of $\delta\tau^{\rm eff}(\omega)$ for
increasing values of $E_{J}$ in an asymmetric configuration with
$\alpha_{L}=2$ and $\alpha_{R}=0.5$. We observe that, for finite $E_{J}$, heat
transport in the forward direction is favored at low frequencies, while the
opposite occurs at higher frequencies. There is not, however, a compensation
in the total heat current and thus one obtains that ${\cal R}>1$ as
illustrated in the inset of Fig. 4. Rectification ratios of the order of
$10\%$ or larger are obtained for the set of parameters in this figure. These
ratios can increase further by decreasing the mean temperature
$\bar{T}=(T_{L}+T_{R})/2$ (see inset in Fig. 4).

Figure 4: Asymmetry in the heat transmission coefficient $\delta\tau^{\rm
eff}(\omega)=\tau^{\rm eff(+)}(\omega)-\tau^{\rm eff(-)}(\omega)$ for the
parallel configuration with $\eta_{L}=1/\pi$, $\eta_{R}=1/4\pi$, $\Delta
T=0.1$, ${\bar{T}}=(T_{L}+T_{R})/2=0.15$, and increasing values of $E_{J}$
from 0 to 0.4. The corresponding heat rectification ratio ${\cal
R}=|J^{+}/J^{-}|$ as a function of $E_{J}$ is shown as a full line in the
inset. The dashed lines illustrate its variation with $\bar{T}$ from 0.175 to
0.1.

Further insight on these properties can be obtained from the analytical
approximation to the system self-energies. As shown in [25] the rectification
ratio in the parallel configuration can be expressed as

${\cal
R}-1\sim(mE_{J})^{2}\frac{\eta_{R}-\eta_{L}}{\eta_{L}+\eta_{R}}\left(\frac{2}{\alpha}-2\right)\bar{T}^{\frac{2}{\alpha}-2}\frac{\delta
T}{\bar{T}}$ (23)

in the limit $T_{L,R}\ll\gamma$ and $\alpha\gtrsim 1$ upon a small temperature
difference $\delta T\ll\bar{T}$. This expression accounts for the main
features of the numerical results, i.e. the quadratic increase with $E_{J}$
and the increase of ${\cal R}$ as the mean temperature is decreased. Eq. (23)
also indicates that, as far as $\alpha>1$, the heat current is larger when the
colder reservoir corresponds to the one with larger resistance. This behavior
is similar to the case of a two-level system with a nonseparable coupling to
the heat baths [29, 34].

Conclusions.— In summary, we have presented non-equilibrium Green functions
calculations for photonic heat transport through a Josephson junction in a
resistive environment. Our results are in qualitative agreement with the
experimental ones from [17] and suggest how signatures of the SB transition
can be detected in heat transport measurements at sufficiently low
temperatures [35]. We further demonstrate rectification properties for this
device that can be tested in future experiments.

Note added: While completing this manuscript we become aware of a related
study [36] for the same system in the scaling regime, yielding results which
are complementary to the ones in the present work.

The authors wish to thank M. Houzet and P. Joyez for enlightening discussions.
ALY also acknowledges Prof. U. Weiss for useful comments on the validity of
the perturbative approach adopted in this work. We acknowledge funding from
the Spanish Ministerio de Ciencia e Innovación via grants No.
PID2020-117671GB-I00, PID2019-110125GB-I00 and PID2022-142911NB-I00, and
through the “María de Maeztu” Programme for Units of Excellence in R&D
CEX2018-000805-M. This work was partially funded by the Research Council of
Finland Centre of Excellence program grant 336810 and grant 349601 (THEPOW).

## References

  * Tafuri [2019] F. Tafuri, ed.,  _Fundamentals and Frontiers of the Josephson Effect_ (Springer International Publishing, Cham, Switzerland, 2019).
  * Caldeira and Leggett [1983] A. O. Caldeira and A. J. Leggett, Quantum tunnelling in a dissipative system, Ann. Phys. 149, 374–456 (1983).
  * Leppäkangas  _et al._ [2013] J. Leppäkangas, G. Johansson, M. Marthaler, and M. Fogelström, Nonclassical photon pair production in a voltage-biased josephson junction, Phys. Rev. Lett. 110, 267004 (2013).
  * Jebari  _et al._ [2018] S. Jebari, F. Blanchet, A. Grimm, D. Hazra, R. Albert, P. Joyez, D. Vion, D. Estève, F. Portier, and M. Hofheinz, Near-quantum-limited amplification from inelastic Cooper-pair tunnelling, Nature Electronics 1, 223–227 (2018).
  * Schmid [1983] A. Schmid, Diffusion and localization in a dissipative quantum system, Phys. Rev. Lett. 51, 1506–1509 (1983).
  * Bulgadaev [1984] S. A. Bulgadaev, Phase diagram of a dissipative quantum system, JETP Lett. 39, 264–267 (1984).
  * Masuki  _et al._ [2022] K. Masuki, H. Sudo, M. Oshikawa, and Y. Ashida, Absence versus Presence of Dissipative Quantum Phase Transition in Josephson Junctions, Phys. Rev. Lett. 129, 087001 (2022).
  * Sépulcre  _et al._ [2023] T. Sépulcre, S. Florens, and I. Snyman, Comment on “Absence versus Presence of Dissipative Quantum Phase Transition in Josephson Junctions”, Phys. Rev. Lett. 131, 199701 (2023).
  * Altimiras  _et al._ [2023] C. Altimiras, D. Esteve, Ç. Girit, H. le Sueur, and P. Joyez, Absence of a dissipative quantum phase transition in Josephson junctions: Theory , arXiv:2312.14754 (2023).
  * Kashuba and Riwar [2023] O. Kashuba and R.-P. Riwar, Dissipation and phase transitions in superconducting circuits coupled to normal metal resistor , arXiv:2312.15951 (2023).
  * Giacomelli and Ciuti [2023] L. Giacomelli and C. Ciuti, Emergent quantum phase transition of a Josephson junction coupled to a high-impedance multimode resonator, arXiv:2307.06383 (2023).
  * Murani  _et al._ [2020] A. Murani, N. Bourlet, H. le Sueur, F. Portier, C. Altimiras, D. Esteve, H. Grabert, J. Stockburger, J. Ankerhold, and P. Joyez, Absence of a Dissipative Quantum Phase Transition in Josephson Junctions, Phys. Rev. X 10, 021003 (2020).
  * Hakonen and Sonin [2021] P. J. Hakonen and E. B. Sonin, Comment on “Absence of a Dissipative Quantum Phase Transition in Josephson Junctions”, Phys. Rev. X 11, 018001 (2021).
  * Murani  _et al._ [2021] A. Murani, N. Bourlet, H. le Sueur, F. Portier, C. Altimiras, D. Esteve, H. Grabert, J. Stockburger, J. Ankerhold, and P. Joyez, Reply to Comment on “Absence of a Dissipative Quantum Phase Transition in Josephson Junctions", Phys. Rev. X 11, 018002 (2021).
  * Kuzmin  _et al._ [2023] R. Kuzmin, N. Mehta, N. Grabon, R. A. Mencia, A. Burshtein, M. Goldstein, and V. E. Manucharyan, Observation of the Schmid-Bulgadaev dissipative quantum phase transition, arXiv:2304.05806 (2023).
  * Pekola and Karimi [2021] J. P. Pekola and B. Karimi, Colloquium: Quantum heat transport in condensed matter systems, Rev. Mod. Phys. 93, 041001 (2021).
  * Subero  _et al._ [2023] D. Subero, O. Maillet, D. S. Golubev, G. Thomas, J. T. Peltonen, B. Karimi, M. Marín-Suárez, A. L. Yeyati, R. Sánchez, S. Park, and J. P. Pekola, Bolometric detection of coherent Josephson coupling in a highly dissipative environment, Nature Comm. 14, 7924 (2023).
  * Thomas  _et al._ [2019] G. Thomas, J. P. Pekola, and D. S. Golubev, Photonic heat transport across a Josephson junction, Phys. Rev. B 100, 094508 (2019).
  * Ojanen and Jauho [2008] T. Ojanen and A.-P. Jauho, Mesoscopic photon heat transistor, Phys. Rev. Lett. 100, 155902 (2008).
  * Keldysh [1964] L. V. Keldysh, Diagram technique for nonequilibrium processes, Zh. Eksp. Teor. Fiz 47, 1515 (1964).
  * Keldysh [1965] L. V. Keldysh, Sov. Phys. JETP 20, 1018 (1965).
  * Meir and Wingreen [1992] Y. Meir and N. S. Wingreen, Landauer formula for the current through an interacting electron region, Phys. Rev. Lett. 68, 2512–2515 (1992).
  * Pendry [1983] J. B. Pendry, Quantum limits to the flow of information and entropy, J. Phys. A: Math. Gen. 16, 2161–2171 (1983).
  * Eckern and Pelzer [1987] U. Eckern and F. Pelzer, Quantum Dynamics of a Dissipative Object in a Periodic Potential, Europhys. Lett. 3, 131 (1987).
  * [25] See supplemental information.
  * Ingold and Nazarov [1992] G.-L. Ingold and Y. V. Nazarov, Charge tunneling rates in ultrasmall junctions, in  _Single Charge Tunneling_, edited by H. Grabert and M. Devoret (Springer, 1992) Chap. 2, p. 21.
  * Aslangul  _et al._ [1987] C. Aslangul, N. Pottier, and D. Saint-James, Quantum Brownian motion in a periodic potential: a pedestrian approach, J. Phys. France 48, 1093–1110 (1987).
  * Fisher and Zwerger [1985] M. P. A. Fisher and W. Zwerger, Quantum Brownian motion in a periodic potential, Phys. Rev. B 32, 6190–6206 (1985).
  * Segal and Nitzan [2005] D. Segal and A. Nitzan, Spin-boson thermal rectifier, Phys. Rev. Lett. 94, 034301 (2005).
  * Martínez-Pérez and Giazotto [2013] M. J. Martínez-Pérez and F. Giazotto, Efficient phase-tunable Josephson thermal rectifier, Appl. Phys. Lett. 102, 182602 (2013).
  * Senior  _et al._ [2020] J. Senior, A. Gubaydullin, B. Karimi, J. T. Peltonen, J. Ankerhold, and J. P. Pekola, Heat rectification via a superconducting artificial atom, Commun. Phys. 3, 40 (2020).
  * Upadhyay  _et al._ [2024] R. Upadhyay, D. S. Golubev, Y.-C. Chang, G. Thomas, A. Guthrie, J. T. Peltonen, and J. P. Pekola, Microwave quantum diode, Nat. Commun. 15, 630 (2024).
  * Arrachea [2023] L. Arrachea, Energy dynamics, heat production and heat–work conversion with qubits: toward the development of quantum machines, Rep. Prog. Phys. 86, 036501 (2023).
  * [34] Note that the temperature behaviour of the rectification is strongly system and parameter dependent.
  * [35] We notice that our calculations assume a purely Ohmic environment. Departures from this idealized case might lead to deviations from our predictions in an actual experiment [37].
  * Yamamoto  _et al._ [2024] T. Yamamoto, L. Glazman, and M. Houzet, Thermal transport across a josephson junction in a dissipative environment (2024).
  * Daviet and Dupuis [2023] R. Daviet and N. Dupuis, Nature of the Schmid transition in a resistively shunted Josephson junction, Phys. Rev. B 108, 184514 (2023).

## Supplementary information for “Photonic heat transport through a Josephson
junction in a resistive environment"

##  I Details on numerical calculations

###  I.1 Self-energies

From Eq. (8) in the main text, the zero order GFs in time representation are
given by

$\displaystyle D_{0}^{r}(t)$ $\displaystyle=$
$\displaystyle-\frac{\theta(t)}{m\omega_{1}}e^{-\gamma
t/2}\sinh\left(\omega_{1}t\right)$ $\displaystyle D_{0}^{K}(t)$
$\displaystyle=$ $\displaystyle S_{1}(t)+S_{2}(t)$ (S 1)

where $\gamma=(\eta_{L}+\eta_{R})/m$,
$\omega_{1}=\sqrt{\gamma^{2}/4-\omega_{c}^{2}}$ and

$\displaystyle S_{1}(t)$ $\displaystyle=$
$\displaystyle-\sum_{\mu=L,R}\frac{\eta_{\mu}}{2m\omega_{1}\gamma}\left(\coth\left(\frac{\beta_{\mu}\lambda_{1}}{2}\right)e^{i\lambda_{1}|t|}-\coth\left(\frac{\beta_{\mu}\lambda_{2}}{2}\right)e^{i\lambda_{2}|t|}\right)$
$\displaystyle S_{2}(t)$ $\displaystyle=$ $\displaystyle
4i\sum_{\mu=L,R}\frac{\eta_{\mu}}{m\beta_{\mu}}\sum_{n>0}\frac{\nu^{\mu}_{n}e^{-\nu^{\mu}_{n}|t|}}{\left(\omega_{c}^{2}+(\nu^{\mu}_{n})^{2}\right)^{2}-\gamma^{2}(\nu^{\mu}_{n})^{2}}\;,$
(S 2)

where $\lambda_{1,2}=i(\gamma/2\pm\omega_{1})$, and
$\nu_{n}^{\nu}=2n/\beta_{\nu}$ are the Matsubara frequencies.

By introducing these analytical expressions for the propagators in time
representation into Eqs. (12) in the main text and performing their Fourier
transform numerically, we obtain the results shown in Fig. 2 for the different
self-energies in frequency space. For comparison we show the corresponding
results for the series configuration.

Figure S 1: Self-energies in frequency representation obtained numerically
from Eqs. (12) in the main text. Energies in units of $E_{c}$. The
corresponding parameters are indicated in the figure title. The dashed lines
correspond to the self-energies in the series configuration for the same set
of parameters. Figure S 2: Inelastic contribution to the heat current in the
series configuration for the same parameters as in Fig. 3 in the main text.

In Fig. S 2 we show the inelastic contribution to the heat current in the
series configurations for the same parameters as in Fig. 3 in the main text.
As can be observed, this contribution is positive so that it tends to increase
the total heat current when added to the elastic contribution.

##  II Analytical approximations

Here we give details on the approximations used to obtain the analytic
behaviour of the self-energies $\Sigma^{r(2)}(\omega)$ and
$\Sigma^{K(2)}(\omega)$ in the parallel configuration.

###  II.1 Retarded self-energy

We start with the time-dependent retarded selfenergy as expressed in Eq. (12)
in the main text:

$\Sigma^{r(2)}(t)=E_{J}^{2}\sin\frac{D_{0}^{r}(t)}{2}\exp\left[\frac{i}{2}\left(D_{0}^{K}(t)-D^{K}_{0}(0)\right)\right]-B\delta(t),$
(S 3)

of which we want to obtain the Fourier transform. For this, we consider the
first term:

$\tilde{\Sigma}(t)=E_{J}^{2}\sin\frac{D_{0}^{r}(t)}{2}\exp\left[\frac{i}{2}\left(D_{0}^{K}(t)-D^{K}_{0}(0)\right)\right].$
(S 4)

The second term, proportional to $B$ will be cancelled out with
$\tilde{\Sigma}(0)$. We can simplify it by writing
$D_{0}^{r}(t)=D_{0}^{a}(-t)=-\theta(t)(1-\exp[-\gamma t])/m\gamma$, and using
the approximation

${\rm
Re}\left[\frac{i}{2}(D_{0}^{K}(t){-}D_{0}^{K}(0))\right]\approx-\sum_{\mu}\frac{\eta_{\mu}}{\pi
m^{2}\gamma^{2}}\ln\left[\frac{\gamma\beta_{\mu}}{\pi}\sinh\left(\frac{\pi|t|}{\beta_{\mu}}\right)\right]$
(S 5)

for $|\gamma t|\gg 1$ and $k_{\text{B}}T\ll\gamma$ [24]. Doing so, we get:

$\tilde{\Sigma}(t)=\theta(t)E_{J}^{2}\sin\left(\frac{e^{-\gamma
t}-1}{2m\gamma}\right)\prod_{\mu}\left[\frac{\gamma\beta_{\mu}}{\pi}\sinh\left(\frac{\pi|t|}{\beta_{\mu}}\right)\right]^{-\eta_{\mu}/\pi
m^{2}\gamma^{2}},$ (S 6)

where $\eta_{\mu}=R_{Q}/2\pi R_{\mu}$, $m=\hbar/E_{C}$ and
$\gamma=\sum_{\mu}\eta_{\mu}/m$, as defined in the main text. We need to solve
its Fourier transform:

$\tilde{\Sigma}(\omega)={\cal
A}\gamma\int_{0}^{\infty}dt\sin\left(\frac{e^{-\gamma
t}{-}1}{2m\gamma}\right)\prod_{\mu}\left[2\sinh\left(\frac{\pi|t|}{\beta_{\mu}}\right)\right]^{-\frac{\eta_{\mu}}{\pi
m^{2}\gamma^{2}}}e^{i\omega t},$ (S 7)

where we have defined the prefactor

${\cal
A}\equiv\frac{E_{J}^{2}}{\gamma}\prod_{\mu}\left(\frac{\gamma\beta_{\mu}}{2\pi}\right)^{-\eta_{\mu}/\pi
m^{2}\gamma^{2}}.$ (S 8)

To reach the low frequency limit $T_{j}\gg\omega$, we approximate
$\sinh(\pi|t|/\beta_{\mu})\approx e^{\pi|t|/\beta_{\mu}}/2$. Then, performing
the change of variables $u=e^{-\gamma t}$ we get:

$\tilde{\Sigma}(\omega)={\cal
A}\int_{0}^{1}du\sin\left(\frac{u{-}1}{2m\gamma}\right)u^{c-1},$ (S 9)

with

$c\equiv\frac{1}{m^{2}\gamma^{3}}\left(\frac{\eta_{L}}{\beta_{L}}+\frac{\eta_{R}}{\beta_{R}}\right)-\frac{i\omega}{\gamma}.$
(S 10)

The integral in Eq. (S 9) has an analytical solution if Re$(c)>0$:

$\int_{0}^{1}du\sin[a(u{-}1)]u^{c-1}=g_{1}(a,c)+g_{2}(a,c),$ (S 11)

with

$\displaystyle\begin{aligned}
g_{1}(a,c)&=\frac{\cos(a)}{(1{+}c)/a}\leavevmode\nobreak\
{{}_{1}F_{2}}\left(\frac{1{+}c}{2};\frac{3}{2},\frac{3{+}c}{2};-\frac{a^{2}}{4}\right)\\\
g_{2}(a,c)&=-\frac{\sin(a)}{c}\leavevmode\nobreak\
{{}_{1}F_{2}\left(\frac{c}{2};\frac{1}{2},\frac{2{+}c}{2};-\frac{a^{2}}{4}\right)},\end{aligned}$
(S 12)

and where ${}_{1}F_{2}(a_{1};b_{1},b_{2};z)$ is a “generalized hypergeometric
function". Formally, this then gives an analytical solution:

$\displaystyle\tilde{\Sigma}(\omega)={\cal
A}\left[g_{1}\left(\frac{1}{2m\gamma},c\right)+g_{2}\left(\frac{1}{2m\gamma},c\right)\right],$
(S 13)

but it is not very informative until we plot its real and imaginary parts,
which is done in Fig. S 3.

Figure S 3: Self-energies for different temperatures $T_{\mu}=(T_{L},T_{R})$
and damping parameters $\eta_{\mu}=(\eta_{L},\eta_{R})$. (a) and (b):
asymmetric case with $\eta_{L}=1/2\pi$ and $\eta_{R}=1/4\pi$. (c) and (d):
symmetric case with $\eta_{L}=\eta_{R}=1/4\pi$. Dashed lines show the
numerical results. In all panels, $E_{J}=0.2$, $m=1$.  Figure S 4:  Comparison
of numeric and analytic selfenergies in the low impedance regime,
$R_{\mu}=0.5R_{Q}$, for increasing temperatures. Other parameters are as in
Fig. S 3.

In a linear expansion, the constant term is real, giving the value of $B$ in
Eq. (S 3), and the linear term is pure imaginary, so the linear expansion
simplifies to:

$\tilde{\Sigma}(\omega)=i{\cal A}{\rm
Im}\left[\left.(g^{\prime}_{1}+g^{\prime}_{2})\right|_{\omega=0}\right]\omega.$
(S 14)

###  II.2 Keldysh self-energy

Now we consider the time-dependent Keldysh self-energy in Eq. (12) in the main
text:

$\displaystyle\begin{aligned}
\Sigma^{K(2)}(t)=\frac{E_{J}^{2}}{i}\cos\left(\frac{D_{0}^{r}(t)}{2}\right)\cos\left(\frac{D_{0}^{a}(t)}{2}\right)\exp\left[\frac{i}{2}\left(D_{0}^{K}(t){-}D^{K}_{0}(0)\right)\right].\end{aligned}$
(S 15)

We separate the time integral in the Fourier transform in two terms, with
$D_{0}^{r}(t<0)=0$ and $D_{0}^{a}(t>0)=0$:

$\displaystyle\Sigma^{K(2)}(\omega)=$ $\displaystyle-i{\cal
A}\gamma\left\\{\int_{-\infty}^{0}dt\cos\left(\frac{e^{\gamma
t}{-}1}{2m\gamma}\right)\prod_{\mu}\left[2\sinh\left(\frac{-\pi
t}{\beta_{\mu}}\right)\right]^{\frac{-\eta_{\mu}}{\pi
m^{2}\gamma^{2}}}e^{i\omega t}\right.$
$\displaystyle\left.+\int_{0}^{\infty}dt\cos\left(\frac{e^{-\gamma
t}{-}1}{2m\gamma}\right)\prod_{\mu}\left[2\sinh\left(\frac{\pi
t}{\beta_{\mu}}\right)\right]^{\frac{-\eta_{\mu}}{\pi
m^{2}\gamma^{2}}}e^{i\omega t}\right\\},$ (S 16)

where we have again used the approximation in Eq. (S 5) [24]. Doing a change
$t\rightarrow-t$ in the first integral of Eq. (S 16), we can write

$\Sigma^{K(2)}(t)=-i{\cal A}\gamma\left({\cal I}^{+}+{\cal I}^{-}\right),$ (S
17)

with

${\cal I}^{\pm}=\int_{0}^{\infty}dt\cos\left(\frac{e^{-\gamma
t}{-}1}{2m\gamma}\right)\prod_{\mu}\left[2\sinh\left(\frac{\pi
t}{\beta_{\mu}}\right)\right]^{\frac{-\eta_{\mu}}{\pi m^{2}\gamma^{2}}}e^{\pm
i\omega t}.$ (S 18)

Following the same steps as for the retarded case ($u=e^{-\gamma t}$), we
arrive to:

$\displaystyle\Sigma^{K(2)}=-i{\cal A}$
$\displaystyle\int_{0}^{1}du\cos\left(\frac{u{-}1}{2m\gamma}\right)\left(u^{c-1}+u^{c^{*}-1}\right)$
(S 19)

In this case, we can use the integral:

$\int_{0}^{1}du\cos[a(u{-}1)]u^{c-1}=g_{3}(a,c)+g_{4}(a,c),$ (S 20)

with

$\displaystyle\begin{aligned}
g_{3}(a,c)&=\frac{\cos(a)}{c}\leavevmode\nobreak\
{{}_{1}F_{2}\left(\frac{c}{2};\frac{1}{2},\frac{2{+}c}{2};-\frac{a^{2}}{4}\right)}\\\
g_{4}(a,c)&=a\frac{\sin(a)}{1{+}c}\leavevmode\nobreak\
{{}_{1}F_{2}\left(\frac{1{+}c}{2};\frac{3}{2},\frac{3{+}c}{2};-\frac{a^{2}}{4}\right)},\end{aligned}$
(S 21)

provided Re$(c)>0$. We note that
$g_{3}(1/2m\gamma,c^{*})=g_{3}^{*}(1/2m\gamma,c)$ and
$g_{4}(1/2m\gamma,c^{*})=g_{4}^{*}(1/2m\gamma,c)$. Hence, the Keldysh self-
energy is purely imaginary:

$\displaystyle\Sigma^{K(2)}=-i2{\cal A}{\rm
Re}\left[g_{3}\left(\frac{1}{2m\gamma},c\right)+g_{4}\left(\frac{1}{2m\gamma},c\right)\right].$
(S 22)

The result is plotted in orange lines in Fig. S 3.

###  II.3 Frequency expansion

We can expand on the integrand of Eqs. (S 9) and (S 19) for low frequencies,
getting:

$\Sigma^{\alpha(2)}\approx\Sigma_{0}^{\alpha(2)}+\omega\Sigma_{1}^{\alpha(2)}+\omega^{2}\Sigma_{2}^{\alpha(2)},$
(S 23)

with

$\displaystyle\Sigma_{k}^{r(2)}$
$\displaystyle=-\frac{1}{k!}\left(\frac{i}{\gamma}\right)^{k}{\cal A}{\cal
I}_{k}(c_{0})$ (S 24) $\displaystyle\Sigma_{0}^{K(2)}$ $\displaystyle=-2i{\cal
A}\int_{0}^{1}du\cos\left(\frac{u{-}1}{2m\gamma}\right)u^{c_{0}-1}$ (S 25)
$\displaystyle\Sigma_{1}^{K(2)}$ $\displaystyle=0$ (S 26)
$\displaystyle\Sigma_{2}^{K(2)}$ $\displaystyle=\frac{i}{\gamma^{2}}{\cal
A}\int_{0}^{1}du\cos\left(\frac{u{-}1}{2m\gamma}\right)u^{c_{0}-1}\ln^{2}u,$
(S 27)

where

$c_{0}\equiv
c(\omega=0)=\frac{1}{m^{2}\gamma^{3}}\left(\frac{\eta_{L}}{\beta_{L}}+\frac{\eta_{R}}{\beta_{R}}\right)$
(S 28)

and

${\cal
I}_{k}(c_{0})=\int_{0}^{1}du\sin\left(\frac{u{-}1}{2m\gamma}\right)u^{c_{0}-1}\ln^{k}u.$
(S 29)

These integrals again have analytical solutions in terms of hypergeometric
functions of higher order. For this, we define:

${\cal
G}_{pqs}^{(n)}={}_{n}F_{n+1}\left(\frac{p+c_{0}}{2}1_{n};\frac{q}{2},\frac{s+c_{0}}{2}1_{n};-\frac{1}{16m^{2}\gamma^{2}}\right),$
(S 30)

where $1_{n}$ is a list of $n$ units.With these, we write:

$\displaystyle{\cal I}_{0}(c_{0})$
$\displaystyle=\frac{\cos\left(\frac{1}{2m\gamma}\right)}{2m\gamma(1{+}c_{0})}{\cal
G}_{1,3,3}^{(1)}-\frac{\sin\left(\frac{1}{2m\gamma}\right)}{c_{0}}{\cal
G}_{0,1,2}^{(1)}$ $\displaystyle{\cal I}_{1}(c_{0})$
$\displaystyle=-\frac{2m\gamma(1{+}2c_{0})\cos\left(\frac{1}{2m\gamma}\right)\sin\left(\frac{1}{2m\gamma}\right)}{c_{0}^{2}(1{+}c_{0})^{2}}$
$\displaystyle+\frac{\cos\left(\frac{1}{2m\gamma}\right)}{12m^{2}\gamma^{2}(1{+}c_{0})(3{+}c_{0})}\left[\frac{{\cal
G}_{3,5,5}^{(1)}}{1{+}c_{0}}+\frac{{\cal G}_{3,5,5}^{(2)}}{3{+}c_{0}}\right]$
(S 31) $\displaystyle-\frac{\sin\left(\frac{1}{2m\gamma}\right)}{2m\gamma
c_{0}^{2}(2{+}c_{0})^{2}}\left[(2{+}c_{0}){\cal G}_{2,3,4}^{(1)}+c_{0}{\cal
G}_{2,3,4}^{(2)}\right]$ $\displaystyle{\cal I}_{2}(c_{0})$
$\displaystyle=2\left[\frac{\cos\left(\frac{1}{2m\gamma}\right)}{2m\gamma(1{+}c_{0})^{3}}{\cal
G}_{2,3,3}^{(3)}-\frac{\sin\left(\frac{1}{2m\gamma}\right)}{c_{0}^{3}}{\cal
G}_{0,1,2}^{(3)}\right].$

For the Keldysh self-energy, it is sufficient to look at the value at the
origin, $\Sigma^{K(2)}(\omega)=\Sigma_{0}^{K(2)}+{\cal O}(\omega^{2})$:

$\displaystyle\Sigma_{0}^{K(2)}=-2i{\cal
A}\left[\frac{\cos\left(\frac{1}{2m\gamma}\right)}{c_{0}}{\cal
G}_{0,1,2}^{(1)}+\frac{\sin\left(\frac{1}{2m\gamma}\right)}{2m\gamma(1{+}c_{0})}{\cal
G}_{1,3,3}^{(1)}\right].$ (S 32)

####  II.3.1 Low temperature scaling

The linear coefficient of the retarded self-energy corresponds to the damping
parameter correction, $\delta\eta$, defined in the main text. Taking the low
temperature limit in the above expressions it is straightforward to show that
for $T_{L,R}\ll\gamma$ and $m\gamma\sim 1/(2\pi)$

$\delta\eta\sim{\cal
A}\frac{\gamma}{c_{0}^{2}}\equiv\left(\frac{E_{J}}{\gamma}\right)^{2}\frac{\left(\frac{T_{L}}{\gamma}\right)^{\frac{\eta_{L}}{\pi
m^{2}\gamma^{2}}}\left(\frac{T_{R}}{\gamma}\right)^{\frac{\eta_{R}}{\pi
m^{2}\gamma^{2}}}}{\left(\frac{\eta_{L}T_{L}+\eta_{R}T_{R}}{m^{2}\gamma^{2}}\right)^{2}}\;.$
(S 33)

In the limit $T_{L}=T_{R}=T$ the scaling reduces to $\delta\eta\sim
T^{2/\alpha-2}$, where $\alpha=2\pi(\eta_{L}+\eta_{R})$ is the inverse of the
total resistance in units of $R_{Q}$, which signals the SB transition at
$\alpha=1$ [24].

Based on this scaling it is possible to determine the transition point $T^{*}$
for heat transport indicated in Fig. 3(b) of the main text. For that purpose
remind that at low frequency the heat transmission coefficient scales as
$\eta/\bar{\eta}$ for a symmetric case with $\eta_{L}=\eta_{R}=\eta$ (see Eq.
(14) in main text). Thus, in order to get the same heat current for two cases
with different resistances at the same mean temperature we need

$\displaystyle\frac{\bar{\eta}_{1}}{\eta_{1}}=\frac{\bar{\eta}_{2}}{\eta_{2}}\longrightarrow\frac{T^{2/\alpha_{1}-2}}{\gamma_{1}^{2/\alpha_{1}+1}}=\frac{T^{2/\alpha_{2}-2}}{\gamma_{2}^{2/\alpha_{2}+1}}\;.$
(S 34)

Setting $\alpha_{1}=1$ in the limit $\alpha_{2}\rightarrow 1$ we obtain
$T=T^{*}=1/(2\pi e^{3/2})\simeq 0.036$, in agreement with the numerical
results in Fig. 3(b) of the main text.

##  III Rectification

Figure S 5: Rectification effect as a function of the temperature difference
using the approximation of Eq. (S 36), for different temperatures. Thin solid
lines correspond to the linear expansion of Eq. (S 38). Parameters:
$E_{J}=0.2$, $m=1$, $\eta_{L}=1/\pi$, $\eta_{R}=1/4\pi$.

The term $\Sigma_{1}^{r(2)}$ modifies the damping parameter
$\Sigma_{1}^{r(2)}=-i\delta\eta$, which affects the transmission probability.
If the resistances are different, the response is different depending on their
temperature. Linearizing Eq. (14) in the main text at $\omega=0$:

$\tau_{\rm
eff}^{\pm}\approx\frac{4\eta_{L}\eta_{R}}{(\eta_{L}+\eta_{R})^{2}}\left(1-\frac{\delta\eta^{\pm}}{\eta_{L}+\eta_{R}}\right).$
(S 35)

A rectification effect appears when $\delta\tau_{\rm parallel}^{\rm
eff}=\tau_{\rm eff}^{+}-\tau_{\rm eff}^{-}\neq 0$:

$\delta\tau_{\rm parallel}^{\rm
eff}\approx\frac{4\eta_{L}\eta_{R}}{(\eta_{L}+\eta_{R})^{3}}\left(\delta\eta^{-}-\delta\eta^{+}\right).$
(S 36)

We hence need

$\delta\eta^{-}-\delta\eta^{+}=i\left[\Sigma_{1,L}^{r(2)}-\Sigma_{1,R}^{r(2)}\right]$
(S 37)

to be finite, where $\Sigma_{1,\mu}^{r(2)}$ is defined as the linear
coefficient $\Sigma_{1}^{r(2)}$ when terminal $\mu$ is at temperature
$T+\delta T$ and the other one is at $T$. They are obtained using Eqs. (S 24)
and (II.3). The result is plotted in Fig. S 5 as a function of the temperature
difference $\delta T$.

We gain further insight by assuming a small temperature difference $\delta T$,
so we can linearize use the scaling of Eq. (S 33), getting

$\delta\eta^{-}-\delta\eta^{+}\approx(mE_{J})^{2}\frac{\eta_{R}-\eta_{L}}{\eta_{L}+\eta_{R}}\left(\frac{2}{\alpha}-2\right)T^{\frac{2}{\alpha}-2}\frac{\delta
T}{T},$ (S 38)

see Fig. S 5. This expression gives the scaling of the rectificaion
coefficient:

${\cal
R}-1\sim(mE_{J})^{2}\frac{\eta_{R}-\eta_{L}}{\eta_{L}+\eta_{R}}\left(\frac{2}{\alpha}-2\right)T^{\frac{2}{\alpha}-2}\frac{\delta
T}{T}$ (S 39)

given in the main text.


# Object level footprint uncertainty quantification in infrastructure based
sensing

Arpan Kusari  Member, IEEE  Asma Almutairi  Mark E.Gilbert and David J.LeBlanc
Member, IEEE This work was supported by a grant as part of the Ford Motor
Company and University of Michigan collaboration on Center for Smart Vehicles
in Smart World (C-SVSW).A. Kusari, A. Almutairi, M.E. Gilbert and D.J. LeBlanc
are with University of Michigan Transportation Research Institute, University
of Michigan, Ann Arbor, MI 48109 USA (e-mail: {kusari, asmaalm, mgilbe,
leblanc}@umich.edu).

###### Abstract

We examine the problem of estimating footprint uncertainty of objects imaged
using the infrastructure based camera sensing. A closed form relationship is
established between the ground coordinates and the sources of the camera
errors. Using the error propagation equation, the covariance of a given ground
coordinate can be measured as a function of the camera errors. The uncertainty
of the footprint of the bounding box can then be given as the function of all
the extreme points of the object footprint. In order to calculate the
uncertainty of a ground point, the typical error sizes of the error sources
are required. We present a method of estimating the typical error sizes from
an experiment using a static, high-precision LiDAR as the ground truth.
Finally, we present a simulated case study of uncertainty quantification from
infrastructure based camera in CARLA to provide a sense of how the uncertainty
changes across a left turn maneuver.

object detection, footprint, uncertainty

##  1 Introduction

Autonomous driving has long being proposed as a measure to enhance road
safety, promote better traffic efficiency and address issues of large scale
transportation mobility. Autonomous vehicles (AV) rely on a suite of different
sensors, including LiDAR, radar, camera, ultrasonic sensors etc to provide
360$\degree$ perception around the vehicle. The sensing data is then fused in
real time to navigate through the dynamic environment. Different AV companies
are working hard towards bringing a reliable AV and related technologies [9].
However, a cursory glance at AV related current news articles show that
accidents are being continuously reported for AVs on test runs across
different cities in US including fatalities [6]. While there are different
reasons being proposed for the causes of the accidents - inclement weather
conditions, poor lighting conditions, software and hardware faults and
occluded dynamic objects, etc., it is becoming exceedingly clear that relying
solely on onboard sensors would not provide the intended coverage required by
these technologies to operate safely and effectively.

As an alternative, different kinds of Internet of Vehicles communication
mechanisms have being proposed to augment onboard sensing where the
information about other dynamic objects (positional attributes or high level
intent) is shared with the ego vehicle to enhance coverage and provide
robustness to aforementioned limitations. Chief among these is infrastructure
(IX) based sensing where the data from the infrastructure based sensors is
packaged in a sparse fashion and sent to the onboard sensing unit [5]. The
sensors placed on the IX node have to be designed to be low-cost from an
operations perspective where the impetus is to provide continuous and
unfettered monitoring. These IX nodes are then connected to an edge processing
device which processes the sensing information on the fly and transmits to the
vehicles on the road though a basic safety message (BSM). On the ego vehicle’s
end, it not only needs to receive the information about the dynamic objects
that IX node sends but also have a measure of trust on the quality of
information in case of field-of-view (FOV) overlap. Our aim through this
proposed approach is to provide a methodology to denote a quantitative measure
of uncertainty for object level information from IX based sensing. Here we
propose camera based sensing only but our analysis can be extended easily for
other sensors.

Object detection is the field of computer vision dedicated towards locating
and classifying objects in an image or video along with a confidence measure
in the classification [2]. It could be thought of as an amalgamation of image
localization and image classification. Traditional machine learning (ML) based
image processing methods have been applied to extract handcrafted features
such as Histogram of oriented gradients (HOG) [13], Speeded-up robust features
(SURF) [3] etc. Deep learning (DL) based object detection methods have
provided the ability to automatically extract features from a large number of
labeled samples such as Region-based Convolutional Neural Network (R-CNN) and
its variants [4], You only look once (YOLO) [11], Single Shot Detector (SSD)
[8] etc with a much larger jump in classification and localization accuracy
than the traditional ML methods. The DL based object detection methods, thus
perform a classification task (label of the object) along with the regression
task (position in terms of bounding box).

Uncertainty quantification (UQ) in DL based object detection has seen a lot of
research in the last few years [1]. Uncertainty can be neatly divided into two
categories - aleatoric (data-based) and epistemic (model-based). The aleatoric
uncertainty is the irreducible uncertainty present in data that gives rise to
uncertainty in predictions. The objective of UQ in DL based object detection
is to estimate the epistemic uncertainty present in the model while also
estimating the aleatoric uncertainty from the noise in the ground truth
bounding box. However, all these UQ papers do not inherently model the
uncertainty in the physical process of capturing the image through the camera.
Thus, while it can report on the confidence level of the predictions based on
the neural network model, it cannot provide the uncertainty in the ground
coordinates as a function of the bounding box. In this research, we aim to
close this gap by determining the uncertainty in the absolute position of the
object given the bounding box position in the earth fixed frame. We combine
the uncertainties due to position and orientation of the camera with the
object level uncertainty.

Our paper is laid out as follows: Section 2 provides the overview and the
mathematical closed form derivatives for computing the object level
uncertainties, Section 3 provides the typical error sizes of the error sources
based on ground truth experiments and Section 4 provides a simulated case
study for a left turn across intersection with two different vantage points.

##  2 Uncertainty quantification in camera images

###  2.1 Computing the location of a point on the ground plane from a camera
image

Define a camera-fixed coordinate system with $x$ along the boresight and $y$
and $z$ axes parallel to the image plane, as shown in Figure 1. A pinhole
camera model is used, so that the image of a point $p(x,y,z)$ appears at row
$r=f\cdot z\textfractionsolidus x$ and column $c=f\cdot y\textfractionsolidus
x$ in the image, where $f$ is the camera focal length in pixels and $c$ and
$r$ are measured from the center of the image plane. Point $F$ in the figure
is the focal point and can serve as the origin of the camera frame.

Figure 1: Camera coordinate system

The camera is placed on a rigid structure a height $h$ above a ground plane,
as shown in Figure 2. Define an earth-fixed coordinate system $(X,Y,Z)$ with
an origin at $F$ (a height $h$ above a ground plane) with $Z$ pointed downward
and $X$ and $Y$ parallel to the ground plane. The camera’s orientation can be
described as initially identical to the earth-fixed $(X,Y,Z)$ coordinates, but
then the camera is first panned (rotated around its $z$ axis) an angle
$\alpha$, followed by a pitching down of the camera frame (rotated about the
camera $y$ axis) of angle -$\theta$.

Figure 2: Earth-fixed coordinates $(X,Y,Z)$, camera coordinates $(x,y,z)$, a
rotated camera, and a the image of a point on the ground plane

For a given point on the image plane, a ray is defined that begins at
coordinate frame origin $F$, passes through the image plane at a point on the
image plane $(c,r)$, and continues until it intersects the ground plane at
$P$. Denote the distance $D$ as the distance along this ray from $F$ to the
point of intersection with the ground plane, so that the location of $P$ in
camera coordinates is:

$P_{c}=D\begin{bmatrix}f\\\ c\\\ r\\\ \end{bmatrix}1/L$

where $L=\sqrt{f^{2}+c^{2}+r^{2}}$.

Sources of error | Error parameters | Error units | X sensitivity to error | Y sensitivity to error
---|---|---|---|---
Intrinsic calibration
Focal length | $df$ | pixels | $\frac{\partial X}{\partial f}$ | $\frac{\partial Y}{\partial f}$
Column error | $dc$ | pixels | $\frac{\partial X}{\partial c}$ | $\frac{\partial Y}{\partial c}$
Row error | $dr$ | pixels | $\frac{\partial X}{\partial r}$ | $\frac{\partial Y}{\partial r}$
Extrinsic Calibration
Camera location | dX | m | 1 | 0
Camera location | dY | m | 0 | 1
Camera location | dh | m | $\frac{\partial X}{\partial h}$ | $\frac{\partial Y}{\partial h}$
Camera azimuth | $d\alpha$ | rad | $\frac{\partial X}{\partial\alpha}$ | $\frac{\partial Y}{\partial\alpha}$
Camera pitch | $d\theta$ | rad | $\frac{\partial X}{\partial\theta}$ | $\frac{\partial Y}{\partial\theta}$
Camera roll | NA | NA | NA | NA
Ground plane error
Ground to camera height error | $dh=dZ$ | m | $\frac{\partial X}{\partial h}$ | $\frac{\partial Y}{\partial h}$
Imaging errors
Column error | dc | pixels | $\frac{\partial X}{\partial c}$ | $\frac{\partial Y}{\partial c}$
Row error | dr | pixels | $\frac{\partial X}{\partial r}$ | $\frac{\partial Y}{\partial r}$
Resolution
Column error | dc | pixels | $\frac{\partial X}{\partial c}$ | $\frac{\partial Y}{\partial c}$
Row error | dr | pixels | $\frac{\partial X}{\partial r}$ | $\frac{\partial Y}{\partial r}$
Table 1: Error sensitivities

In earth-fixed coordinates, this vector from $F$ to $P$ is simply
$P_{e}=(X,Y,Z)^{T}$ where $Z=h$. The distance $D$ is simply
$D=\sqrt{X^{2}+Y^{2}+Z^{2}}$. We can now compute the location in earth
coordinates using the transformation matrix associated with the two rotations
described earlier – a pan of angle $\alpha$ and a pitch downward of $\theta$:

$P_{e}=T_{c\rightarrow e}P_{c}$ (1)

where $T_{c\rightarrow e}=\begin{bmatrix}cos\theta\cdot
cos\alpha&-sin\alpha&-sin\theta\cdot cos\alpha\\\ cos\theta\cdot
sin\alpha&cos\alpha&-sin\theta\cdot sin\alpha\\\
sin\theta&0&cos\theta\end{bmatrix}$

The third element of this equation 1 then becomes

$Z=h=\dfrac{(sin\theta\cdot f+cos\theta\cdot r)D}{L}$ (2)

So that

$D=\dfrac{hL}{sin\theta\cdot f+cos\theta\cdot r}$

Then X and Y can simply be computed from equation 1

$\displaystyle X$ $\displaystyle=\dfrac{(cos\theta\cdot cos\alpha\cdot
f-sin\alpha\cdot c-sin\theta\cdot cos\alpha\cdot r)D}{L}$ (3) $\displaystyle
Y$ $\displaystyle=\dfrac{(cos\theta\cdot sin\alpha\cdot f+cos\alpha\cdot
c-sin\theta\cdot sin\alpha\cdot r)D}{L}$

This provides the closed-form relationship between the point on the image
plane and the point on the ground plane. We can now compute $(X,Y)$ given the
point on the image plane. (Note we can also compute the image plane projection
of point $P$, given $X$ and $Y$, but we have not shown that straightforward
result.)

###  2.2 Computing sensitivity to errors by closed-form partial derivatives

To determine the influence of an error in one of the input parameters, denoted
$par={h,\alpha,\theta,f,c,r}$ on the accuracy in estimating the location of
point $P$ in earth-fixed coordinates, we can simply compute $(X,Y)$ with
“perfect” information and then perturb one or more of the input parameters and
compute $(X,Y)$ using the perturbation. The difference is the error that is
introduced by the perturbation (error) in the input parameter. Note that we
can also solve for changes in $(c,r)$ caused by a motion (perturbation) in the
known location $(X,Y)$ by writing out $(c,r)$ as a function of $(X,Y)$, by
going back to equation 1 and solving for the image parameters as a function of
$(X,Y)$, as mentioned earlier. Then we could perturb $(X,Y)$ and find how the
image plane location changes.

Errors in camera calibration, knowledge of the ground plane, and so on will
introduce errors in the computed location $(X,Y)$ of the imaged point. Table 1
below shows selected error sources and the parameters whose values are
influenced by the errors. For instance, an incorrect calibration of the
camera’s earth-fixed location $(X,Y)$ will introduce an error $(\Delta
X,\Delta Y)$.

The partial derivatives referred in the table above are given below:

$\displaystyle\begin{split}\frac{\partial X}{\partial
f}={}&\frac{cos\theta\cdot cos\alpha\cdot h}{sin\theta\cdot f+cos\theta\cdot
r}-\\\ &\frac{(cos\theta\cdot cos\alpha\cdot f-sin\alpha\cdot c-sin\theta\cdot
cos\alpha\cdot r)\cdot h\cdot sin\theta}{(sin\theta\cdot f+cos\theta\cdot
r)^{2}}\end{split}$ (4) $\displaystyle\begin{split}\frac{\partial Y}{\partial
f}={}&\frac{cos\theta\cdot sin\alpha\cdot h}{sin\theta\cdot f+cos\theta\cdot
r}-\\\ &\frac{(cos\theta\cdot sin\alpha\cdot f+cos\alpha\cdot c-sin\theta\cdot
sin\alpha\cdot r)\cdot h\cdot sin\theta}{(sin\theta\cdot f+cos\theta\cdot
r)^{2}}\end{split}$ (5) $\displaystyle\frac{\partial X}{\partial c}$
$\displaystyle=\frac{-sin\alpha\cdot h}{sin\theta\cdot f+cos\theta\cdot r}$
(6) $\displaystyle\frac{\partial Y}{\partial c}$
$\displaystyle=\frac{cos\alpha\cdot h}{sin\theta\cdot f+cos\theta\cdot r}$ (7)
$\displaystyle\begin{split}\frac{\partial X}{\partial
r}={}&\frac{-sin\theta\cdot cos\alpha\cdot h}{sin\theta\cdot f+cos\theta\cdot
r}-\\\ &\frac{(cos\theta\cdot cos\alpha\cdot f-sin\alpha\cdot c-sin\theta\cdot
cos\alpha\cdot r)\cdot h\cdot cos\theta}{(sin\theta\cdot f+cos\theta\cdot
r)^{2}}\end{split}$ (8) $\displaystyle\begin{split}\frac{\partial Y}{\partial
r}={}&\frac{-sin\theta\cdot sin\alpha\cdot h}{sin\theta\cdot f+cos\theta\cdot
r}-\\\ &\frac{(cos\theta\cdot sin\alpha\cdot f+cos\alpha\cdot c-sin\theta\cdot
sin\alpha\cdot r)\cdot h\cdot cos\theta}{(sin\theta\cdot f+cos\theta\cdot
r)^{2}}\end{split}$ (9) $\displaystyle\frac{\partial X}{\partial h}$
$\displaystyle=\frac{(cos\theta\cdot cos\alpha\cdot f-sin\alpha\cdot
c-sin\theta\cdot cos\alpha\cdot r)}{sin\theta\cdot f+cos\theta\cdot r}$ (10)
$\displaystyle\frac{\partial Y}{\partial h}$
$\displaystyle=\frac{(cos\theta\cdot sin\alpha\cdot f+cos\alpha\cdot
c-sin\theta\cdot sin\alpha\cdot r)}{sin\theta\cdot f+cos\theta\cdot r}$ (11)
$\displaystyle\frac{\partial X}{\partial\alpha}$
$\displaystyle=\frac{(-cos\theta\cdot sin\alpha\cdot f-cos\alpha\cdot
c+sin\theta\cdot sin\alpha\cdot r)\cdot h}{sin\theta\cdot f+cos\theta\cdot r}$
(12) $\displaystyle\frac{\partial Y}{\partial\alpha}$
$\displaystyle=\frac{(cos\theta\cdot cos\alpha\cdot f-sin\alpha\cdot
c-sin\theta\cdot cos\alpha\cdot r)\cdot h}{sin\theta\cdot f+cos\theta\cdot r}$
(13) $\displaystyle\begin{split}\frac{\partial
X}{\partial\theta}={}&\frac{(-sin\theta\cdot cos\alpha\cdot f-cos\theta\cdot
cos\alpha\cdot r)\cdot h}{sin\theta\cdot f+cos\theta\cdot r}-\\\
&\frac{(cos\theta\cdot cos\alpha\cdot f-sin\alpha\cdot c-sin\theta\cdot
cos\alpha\cdot r)\cdot h\cdot(cos\theta\cdot f-sin\theta\cdot
r)}{(sin\theta\cdot f+cos\theta\cdot r)^{2}}\end{split}$ (14)
$\displaystyle\begin{split}\frac{\partial
Y}{\partial\theta}={}&\frac{(-sin\theta\cdot sin\alpha\cdot f-cos\theta\cdot
sin\alpha\cdot r)\cdot h}{sin\theta\cdot f+cos\theta\cdot r}-\\\
&\frac{(cos\theta\cdot sin\alpha\cdot f+cos\alpha\cdot c-sin\theta\cdot
sin\alpha\cdot r)\cdot h\cdot(cos\theta\cdot f-sin\theta\cdot
r)}{(sin\theta\cdot f+cos\theta\cdot r)^{2}}\end{split}$ (15)

The sensitivity of the computed location to the error is computed using the
partial derivative of X and Y to the error sources as given below:

$Cov(X,Y)=J\cdot Cov(par))\cdot J^{T}$ (17)

where $J$ is the Jacobian comprised of the partial derivatives of X and Y with
respect to the parameters and $Cov(par)$ represents the covariance matrix of
the parameters.

###  2.3 Estimating uncertainty of bounding box

In the previous section, we determined the uncertainty of a given point in the
image. However, for image based object detection, the output usually is an
axis aligned bounding box encompassing the object along with a class label
associated with the object and the corresponding confidence the detector has
in classifying the object. Therefore, we need a pooled uncertainty estimate
for the bounding box. We can compute the uncertainties of the extreme points
of the bounding box but its not entirely not apparent how the uncertainty of
the bounding box can be computed.

Since, we are concerned with the uncertainty of a given imaged point in terms
of ground coordinates, it translates very neatly to estimating the footprint
uncertainty in case of an object. Intuitively, this is also important from an
operations perspective where the uncertainty in the X-Y ground plane is
important for safe navigation. A thing to note is that unless the center of
the bounding box coincides with the center of the camera coordinate system,
the four extreme points will have different uncertainties with the points away
from the center having largest uncertainties. The 3D bounding box of a
simulated vehicle is shown in Figure 3.

There are a few different ways of reporting the uncertainty of a bounding box
footprint:

  * •

Largest point uncertainty - we can report the largest uncertainty of the
object footprint since that provides the upper bound on the positional
uncertainty.

  * •

Footprint center uncertainty - the footprint center is the centroid of the
convex quadrilateral that forms the footprint. Typically, at least one of the
edges are obscured which can complicate the estimation of uncertainty.

Figure 3: 3D bounding box of a simulated vehicle along with the footprint of
the vehicle

##  3 Typical size of error parameters

Given the methodology for examining the effect of the observed parameters have
on the determined ground coordinates, we would like to understand the expected
level of error in each of the observations. Typically, the errors are
determined using a known calibration pattern such as a checkerboard pattern
consisting of alternating white and black squares of equal size. The seminal
work of Zhang [14] provided a simple algorithm to estimate the intrinsic
calibration from a checkerboard pattern given different poses. In the paper,
the author shows that the method is stable and is able to estimate the
parameters with a small amount of variation. But under addition of synthetic
random noise in the model points, the error rate in the parameters increases
almost linearly while under systematic spherical and cylindrical non-
planarity, the scale factors increase linearly as well. We derive inspiration
from this study in trying to understand the typical sizes of the error sources
mentioned in Table 1.

We set up multiple data collections using two of each type of camera listed
below -

  * •

Basler Ace 2 camera with 3.5mm focal length wide angle low distortion lens

  * •

Bad Wolf HD-SDI Mini Cube camera with 2.8mm focal length CS-mount lens
(referred to as BW-Cube henceforth)

  * •

Bad Wolf HD-SDI Mini Bullet camera with 3.6mm focal length S-mount lens
(referred to as BW-Bullet henceforth)

and Ouster OS1 128 beam LiDAR with complete overlap of field of views, mounted
on top of a car as a surrogate for IX, as shown in Figure 4 (a). Our aim is to
test the variability of each source for these different kinds of cameras with
the LiDAR serving as the ground truth. A checkerboard pattern of 14x9 inner
corners precisely printed on a rigid white board is used to perform
calibration. A continuous data collection is performed using the checkerboards
from which the individual frames are extracted. The example of all the valid
targets are given in Figure 4 (b) with respect to the camera coordinate
system.

Below, we present an exhaustive analysis on using this data to estimate the
error sizes of the aforementioned parameters. All the error values quoted are
the standard deviation values.

Figure 4: (a) Ouster LiDAR and Basler Camera rigidly mounted for the
experiment (b) position of the checkerboard patterns with respect to the
camera

###  3.1 Error sizes for intrinsic calibration

Table 2: Variability of intrinsic parameters for the different cameras

Camera type |  $\sigma$fx |  $\sigma$fy |  $\sigma$f |  $\sigma$c |  $\sigma$r |  $\sigma$k1 |  $\sigma$k2 |  $\sigma$p1 |  $\sigma$p2 |  $\sigma$k3
---|---|---|---|---|---|---|---|---|---|---
| (in pixels) | (in pixels) | (in pixels) | (in pixels) | (in pixels) | (in units) | (in units) | (in units) | (in units) | (in units)
Basler1 | 0.1992 | 0.1923 | 0.2768 | 0.1713 | 0.1314 | $5E-4$ | 0.0019 | $3E-5$ | $5E-5$ | 0.0002
Basler2 | 0.1479 | 0.1470 | 0.2085 | 0.1486 | 0.1465 | 0.0004 | 0.0011 | $2E-5$ | $3E-5$ | 0.0009
BW-Cube1 | 0.3057 | 0.3180 | 0.4411 | 0.3031 | 0.2437 | 0.0004 | 0.0005 | $7E-5$ | $8E-5$ | 0.0005
BW-Cube2 | 0.670 | 0.680 | 0.9546 | 0.3075 | 0.3656 | 0.0001 | 0.0001 | $1E-5$ | $3E-5$ | $0.0003$
BW-Bullet1 | $0.4225$ | $0.4255$ | $0.5996$ | $0.4115$ | $0.3405$ | $0.001$ | $0.0014$ | $6E-5$ | $1E-5$ | $0.0009$
BW-Bullet2 | $0.332$ | $0.325$ | $0.4645$ | $0.416$ | $0.299$ | $0.0001$ | $4E-5$ | $5E-5$ | $7E-5$ | $0.0001$

Table 3: Variability of extrinsic parameters for the different cameras

Camera type |  $\sigma$X |  $\sigma$Y |  $\sigma$h | $\sigma\alpha$ | $\sigma\theta$
---|---|---|---|---|---
| (m) | (m) | (m) | ($10^{-3}$ deg) | ($10^{-3}$ deg)
Basler1 | 0.1061 | 0.0861 | 0.1936 | 0.1524 | 0.1480
Basler2 | 0.1106 | 0.1077 | 0.2483 | 0.1488 | 0.1354
BW-Cube1 | 0.4912 | 0.1910 | 0.8484 | 0.8032 | 0.7049
BW-Cube2 | 0.1396 | 0.1103 | 0.1876 | 0.3311 | 0.2784
BW-Bullet1 | 0.0802 | 0.0632 | 0.1744 | 0.3543 | 0.2558
BW-Bullet2 | 0.0658 | 0.0588 | 0.1546 | 0.3123 | 0.2239

The intrinsic parameters represent the optical center and focal length of the
camera along with geometric characteristics of the camera which provides a
framework of how the camera captures a given image. Typically, the intrinsic
parameters are a function of the camera system and as such, provided by the
manufacturer as part of the technical specifications. It can also be very
easily determined using the method given by Zhang which comes standard as part
of the OpenCV camera calibration implementation. As part of the process, the
variability of the intrinsic parameters are outputted as well. Since we
capture a lot of frames with the targets, we restrict the frames which have a
mean reprojection error lesser than a user-defined threshold.

We perform the intrinsic calibration for the cameras and present the
variability of the intrinsic parameters of each camera in Table 3. From the
table, it is quite evident that the uncertainty of parameters is higher in the
BW-Cube and BW-Bullet cameras as opposed to the Basler cameras especially in
the estimates of focal lengths and principal points. In certain cases, the
uncertainty reaches over half a pixel which can have significant effect on the
object uncertainty. We also tried to understand the difference between
individual cameras based on the uncertainty of the parameters. We found that
BW camera lenses had widely varying focal lengths as opposed to the Basler
lenses which is also not surprising given the cost differences between them.

While the previous results demonstrate the differences in uncertainty of the
cameras, we recognize that the variability of the intrinsic parameters is
directly influenced by the number and positions of checkerboard patterns. One
of the main reasons behind this is the uncertainty in the position and
orientation of the model plane from the ‘n’ images presented during the
process. Therefore, in order to understand what the spatial distribution of
the checkerboards are with respect to the image space, we project the convex
hull of each checkerboard onto the image plane. The image space is gridded
into blocks of 5 pixels each and the number of the projected checkerboards on
each grid cell is counted. The resulting figure of histograms for each kind of
camera is given in left hand side of the Figure 7 (Fig. 7(a) Basler, Fig. 7(c)
BW-Cube, Fig. 7(e) BW-Bullet). The distribution shows the concentration of the
checkerboards is somewhat non-symmetric in the latter cases and there are
portions of the image space (specifically at the corners) where there are no
checkerboards, which also have largest amount of distortions.

The two kinds of errors in lenses are radial and tangential. While they
account for the most of the systematic errors in lenses, we are also
interested in knowing if there are non-symmetric errors in lenses that these
two error types do not model. Therefore, we devise a non-parametric approach
towards understanding the total error budget in each grid cell (of 5x5
pixels). We estimate the reprojection error of each checkerboard corner in
each image and project the corner into the image plane. Then we grid the image
space and calculate the mean reprojection error in each grid cell. The results
are provided in the right hand of Figure 7. The maximum mean reprojection
error for the BW cameras is about two times greater than the Basler camera.
Fig. 7(f) shows that there is some non-symmetric error in the right corner of
the image space which needs to be identified.

Figure 5: Distributions of checkerboards in image space (left) and mean
reprojection errors in image space (right) for Basler (first row), BW-Cube
(second row) and BW-Bullet (third row).

###  3.2 Error sizes for extrinsic calibration

Figure 6: Distribution of position and orientation of checkerboards for Basler
(first row), BW-Cube (second row) and BW-Bullet (third row).

The extrinsic parameters provide the position and orientation of the camera
with reference with the ground truth coordinate system (in our case, the
LiDAR). Extrinsic parameter estimation of camera and LiDAR can be done in
variety of ways as shown in the literature. We utilize the well known
Perspective-n-Point (PnP) algorithm to estimate the pose of a calibrated
camera given a set of n 3D world points and their corresponding 2D projections
in the image space. We specifically utilize Efficient PnP [7] which solves the
general problem for $n\geq 4$ combining with random sample consensus (RANSAC)
algorithm to weed out outliers in point correspondences. Similar to previous
analysis, we perform the extrinsic parameter estimation for different cameras
in Table 3 and the corresponding figure is given in Figure 6. The results show
that there is a significant difference in the vertical and angular uncertainty
of extrinsic parameters from Basler cameras as compared to the BW cameras.

###  3.3 Error size for ground plane error

The estimation of the ground $X$ and $Y$ coordinates hinge on estimating the
value of distance $D$, which in turn is estimated based on knowledge of the
height of the camera from the ground plane. Even if the height of the camera
is very accurately measured with the help of surveying grade equipment, the
undulation of the ground plane and distance from the camera can lead to
substantial errors. This can be decomposed neatly in terms of distance $D$
from the camera and the angle of incidence of the local terrain normal.

The distance from the surface to the camera can lead to the lengthening of the
surface area captured by a pixel. Thus, the farther the camera is, the more
spread out the surface area which contributes to a single pixel. The effect is
shown in Figure 7.

Figure 7: Response of the light projected into the camera as a function of
distance

We assume that the surface is perfectly flat in our formulation. However, in
reality, the surface normal has a substantial effect on the light intensity
received by the camera and is usually given as a function of cosine of the
angle with respect to the surface normal as shown in Figure 8.

Figure 8: Response of the light projected into the camera as a function of
surface normal

Since we do not have inherent knowledge of the surface normal, we assume the
largest error in the estimation of surface normal.

###  3.4 Error sizes for imaging and resolution errors

Estimating imaging and resolution errors need greater precision than afforded
by our measuring apparatus and so we turn to existing literature to determine
these errors. Specifically, digital image correlation (DIC) has been adopted
by the solid mechanics community to measure strain fields [12]. DIC utilizes a
set of points in the reference image and obtain the corresponding set in the
deformed image. [10] provides the imaging noise as $0.01~{}pixels$ and
resolution error as $0.001~{}pixels$. We understand that their camera
resolution is far higher than that of our application. So we relax the
uncertainties and use $\Delta r=\Delta c=0.1~{}pixels$ for imaging errors and
$\Delta r=\Delta c=0.01~{}pixels$ for resolution errors.

##  4 Case study

Figure 9: Left turn maneuver with camera 1 with the trajectory samples and the
XY-uncertainties in image frame for each point under different types of
cameras - (a) Basler camera (b) BW-Bullet camera (c) BW-Cube camera Figure 10:
Left turn maneuver with camera 2 with the trajectory samples and the XY-
uncertainties in image frame for each point under different types of cameras -
(a) Basler camera (b) BW-Bullet camera (c) BW-Cube camera

We show two instances of infrastructure mounted camera in the CARLA simulator,
where the first camera captures the diagonal between two intersections while
other captures the vehicle across from the first camera. We run a left turn
maneuver for both cameras and capture the ground truth bounding boxes in world
coordinate-frame, camera coordinate-frame and image coordinate-frame.

Figures 9 and 10 provide the result of two left turn maneuvers from the two
infrastructure mounted cameras. We plot the X and Y-uncertainties in image
frame as ellipses for the three kinds of cameras chosen previously. The first
panels show the uncertainty ellipse over the entire trajectory in image frame
for a Basler camera, the second panels show the uncertainty ellipse for a BW-
Bullet camera and the third panels show the uncertainty ellipse for a BW-Cube
camera. From the figures, we can observe the transition of uncertainties over
the trajectories, with the uncertainty in X dominating in the beginning and
equaling the uncertainty in Y in the middle of the intersection. The
uncertainty in Y dominates from that point onward.

We can also clearly see the differences in uncertainties between the cameras
with the uncertainty in BW-Cube cameras being almost equivalent to the lane
width. This shows that for infrastructure mounted cameras, the uncertainty in
object footprint is an important consideration.

##  5 Conclusion

In this work, a generic approach is proposed to determine object footprint
uncertainty in infrastructure mounted cameras. A closed form relationship
relating the object footprint (X and Y) to the parameters are established and
the covariance matrix of the object footprint uncertainty is derived along
with the individual closed-form partial derivatives. We also looked at the
typical sizes of the error parameters, deriving them from established
calibration techniques and using them to estimate uncertainties for simulated
left turn with two infrastructure mounted cameras. The analysis of the
uncertainty is broad enough that it can be easily applied to determine the
analytical uncertainty for deep learning based models.

## References

  * [1] Moloud Abdar, Farhad Pourpanah, Sadiq Hussain, Dana Rezazadegan, Li Liu, Mohammad Ghavamzadeh, Paul Fieguth, Xiaochun Cao, Abbas Khosravi, U Rajendra Acharya, et al.  A review of uncertainty quantification in deep learning: Techniques, applications and challenges.  Information Fusion, 76:243–297, 2021.
  * [2] Abhishek Balasubramaniam and Sudeep Pasricha.  Object detection in autonomous vehicles: Status and open challenges.  arXiv preprint arXiv:2201.07706, 2022.
  * [3] Herbert Bay, Andreas Ess, Tinne Tuytelaars, and Luc Van Gool.  Speeded-up robust features (surf).  Computer vision and image understanding, 110(3):346–359, 2008.
  * [4] Ross Girshick.  Fast r-cnn.  In Proceedings of the IEEE international conference on computer vision, pages 1440–1448, 2015.
  * [5] Myungsoo Jun and Anthony J Markel.  Infrastructure-based sensors augmenting efficient autonomous vehicle operations.  Technical report, National Renewable Energy Lab.(NREL), Golden, CO (United States), 2017.
  * [6] Philip Koopman.  How safe is safe enough?: Measuring and predicting autonomous vehicle safety.  2022\.
  * [7] Vincent Lepetit, Francesc Moreno-Noguer, and Pascal Fua.  Ep n p: An accurate o (n) solution to the p n p problem.  International journal of computer vision, 81:155–166, 2009.
  * [8] Wei Liu, Dragomir Anguelov, Dumitru Erhan, Christian Szegedy, Scott Reed, Cheng-Yang Fu, and Alexander C Berg.  Ssd: Single shot multibox detector.  In Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11–14, 2016, Proceedings, Part I 14, pages 21–37. Springer, 2016.
  * [9] Piergiuseppe Mallozzi, Patrizio Pelliccione, Alessia Knauss, Christian Berger, and Nassar Mohammadiha.  Autonomous vehicles: state of the art, future trends, and challenges.  Automotive systems and software engineering: State of the art and future trends, pages 347–367, 2019.
  * [10] Bing Pan, Liping Yu, and Dafang Wu.  Accurate ex situ deformation measurement using an ultra-stable two-dimensional digital image correlation system.  Applied optics, 53(19):4216–4227, 2014.
  * [11] Joseph Redmon, Santosh Divvala, Ross Girshick, and Ali Farhadi.  You only look once: Unified, real-time object detection.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 779–788, 2016.
  * [12] Michael A Sutton, Jean Jose Orteu, and Hubert Schreier.  Image correlation for shape, motion and deformation measurements: basic concepts, theory and applications.  Springer Science & Business Media, 2009.
  * [13] Carlo Tomasi.  Histograms of oriented gradients.  Computer Vision Sampler, pages 1–6, 2012.
  * [14] Zhengyou Zhang.  A flexible new technique for camera calibration.  IEEE Transactions on pattern analysis and machine intelligence, 22(11):1330–1334, 2000.