\section{Introduction}

A large Hall angle approaching $\pi/2$ is one of the key parameters characterizing a high-mobility two-dimensional electron system (2DES) subjected to a perpendicular magnetic field. The electron density $n_\mathrm{e}$ in a 2DES residing at a GaAs/AlGaAs heterointerface is typically $\sim$5$\times$10$^{15}$ m$^{-2}$ or less. There the Hall resistivity $\rho_{xy} \sim B/(n_\mathrm{e}e)$ far exceeding the diagonal resistivity $\rho_{xx}$ can be achieved with a relatively small magnetic field ($\sim$0.1 T), resulting in a large Hall angle $\delta = \arctan (\rho_{xy}/\rho_{xx})$. (We assume throughout the paper that the 2DES is isotropic. Therefore we have $\rho_{xx} = \rho_{yy}$, $\rho_{yx} = -\rho_{xy}$ and the corresponding relations for the conductivity). Let us consider a rectangular sample of a 2DES with the source and drain electrodes attached at opposite ends. The spatial distributions of the electric field ${\bf E}$ and the current density ${\bf j}$ induced in such a sample by the application of a source-drain bias have been calculated employing  classical electromagnetism \cite{Wick54,Wakabayashi78,Rendell81,Neudecker87}. At zero magnetic field, both ${\bf E}$ and ${\bf j}$ are parallel to the longitudinal direction (the direction from the source to the drain). A magnetic field drastically alters the distributions of ${\bf E}$ and ${\bf j}$, primarily by introducing the angle $\delta$ between the two vectors. Rather counterintuitively, ${\bf E}$ becomes nearly perpendicular ($\sim \delta \simeq \pi/2$) to the longitudinal direction, barring the regions close to the electrodes. Small areas with highly concentrated ${\bf j}$ and the equipotential lines, dubbed ``hot spots'', emerge at the two diagonally opposite corners, one facing the source and the other facing the drain electrode. The predicted spatial distributions have been essentially verified by various experimental techniques \cite{Fontein92,Yacoby91,Ahlswede01,Kawano10JJAP}.

Replacing the source-drain bias with a temperature difference, we expect analogous spatial distributions in the temperature gradient $\boldsymbol{\nabla} T$ and the thermal flux density ${\bf j}_\mathrm{Q}$. Good analogy is expected when the electron and the lattice systems are virtually decoupled and major part of ${\bf j}_\mathrm{Q}$ is carried by electrons, as is the case in the GaAs/AlGaAs 2DES at very low temperatures (below $\sim$200 mK). The distribution of $\boldsymbol{\nabla} T$, especially the reorientation of $\boldsymbol{\nabla} T$ by a magnetic field, can profoundly affect the thermoelectric voltage generated by the temperature difference. In the analysis of the thermoelectric voltages of a 2DES, however, $\boldsymbol{\nabla} T$ is usually assumed to simply remain parallel to the longitudinal direction even in the presence of the magnetic field. Strangely, the redistribution of $\boldsymbol{\nabla} T$ has not attracted due attention it deserves. (The effect of the redistribution of $\boldsymbol{\nabla} T$ on the thermoelectric voltages was previously considered for bulk BiSb and InSb alloys \cite{Ertl63,Okumura98,Heremans00,Heremans01}.)

The purpose of the present study is to interpret experimentally measured thermoelectric voltages in terms of $\boldsymbol{\nabla} T$ calculated considering the effect of the magnetic field. Measurements were performed on a 2DES having a Hall-bar geometry. We measure voltages between several pairs of voltage probes and examine how the voltage varies with varying location on the Hall bar, focusing on the amplitude of the quantum oscillations due to the Landau quantization. An important difference between ${\bf j}$ and ${\bf j}_\mathrm{Q}$ in Hall-bar measurements is that the latter is allowed to flow into the voltage contacts, while the former is prohibited. Therefore, the arms of the voltage probes can host non-vanishing $\boldsymbol{\nabla} T$ and thus can contribute to the thermoelectric voltages \cite{Fletcher86}. We numerically calculate the spatial distribution of the electron temperature $T$, taking into account the thermal diffusion through the arms into the contacts as well as the power transferred to the lattice via electron-phonon interaction. The lattice temperature is assumed to be uniform and kept at the lowest temperature in the system. We find that the dependence of the measured thermoelectric voltages on the location and on the magnetic field can be qualitatively explained by the calculated spatial distribution of $T$. The calculations highlight crucial roles played by the voltage probes in reducing the temperature and in generating the thermoelectric voltage, when placed in a magnetic field.
The geometry of the voltage-probe arms designed to minimize the thermal disturbance turns out to be effective only at around zero magnetic field.

Thermoelectric voltages are sensitive to the energy dependence of the conductivity and also are a measure of the entropy of the system \cite{Harman67,Callen85,Pottier10}, and thus have been extensively applied to the studies of transport properties and scattering mechanisms. Numbers of studies have been performed on the thermoelectric properties of a 2DES (see \cite{Gallagher92,Fletcher99} for reviews). Recent theoretical suggestion of the possibility to probe, through the entropy, non-Abelian quasiparticles in $\nu = 5/2$ fractional quantum Hall state \cite{Yang09,Barlas12} revitalized the interest in the thermoelectric voltages of a 2DES subjected to a magnetic field \cite{Chickering10,Chickering13,Kobayakawa13}. The present study suggests precautions to be taken in analysing the thermoelectric voltages measured in a magnetic field.

The paper is organized as follows. In section \ref{meas}, we describe experimental details and the results of the thermoelectric-voltage measurements. Section \ref{calcT} is devoted to the calculation of the spatial distribution of $T$. We start by briefly describing simple one-dimensional and rectangular models neglecting the electron-phonon interaction. Analytic solutions given for these models help us grasp the essence of the role played by a magnetic field. Numerical calculations follow, with which we examine the effect of the voltage probes and of the electron-phonon interaction. Thermoelectric voltages resulting from the distribution of $T$ thus calculated are presented in section \ref{thermvol}. We use thermoelectric coefficients deduced relying on the generalized Mott's formula. The calculated thermoelectric voltages are compared with the measurements presented in section \ref{meas}. We also evaluate separately the contribution of the main bar and of the voltage-probe arms to the thermoelectric voltages. Discussion on quantitative disparities between measured and calculated thermoelectric voltages, as well as on the characteristics and the limitations of the measurements in the present study, is given in section \ref{discussion}. Section \ref{conclusion} concludes the paper.

\section{Measurements \label{meas}}
\begin{figure*}[tb]
\begin{center}
\includegraphics[width=11cm,clip]{figure1.jpg}
\caption{Schematic diagram of the Hall-bar device used in the present study. The main (horizontal) Hall bar (from 1 to 7, length of the effective part $L$, width $2W$) is fitted with six voltage probes (4, 5, 6, 8, 9, 10) to measure the thermoelectric voltages. The arms of the voltage probes are designed to be thin ($2W_\mathrm{V}$) and long ($L_\mathrm{V}$), in order to minimize the thermal disturbance.
The secondary (vertical) Hall bar (from 3 to 11) serves as a heater and also contains two voltage probes (2 and 12) to monitor the electron temperature by the amplitude of the SdH oscillations. The light-green rectangle encompasses the main Hall bar, i.e., the main bar from the heater to the heat sink (7) and the arms of the voltage probes, for which the spatial distribution of the electron temperature is calculated in section \ref{calcTFEM} and plotted in figure \ref{TAephcalc}. Upper right inset: optical micrograph of the device. The dimensions of the Hall bar, including the locations of (the centres of) the arms, are listed to the right of the figure.}
\label{HallBarFig}
\end{center}
\end{figure*}

Thermoelectric-voltage measurements were performed on a Hall bar device depicted in figure \ref{HallBarFig}, fabricated from a GaAs/AlGaAs 2DES wafer having the electron density $n_\mathrm{e} = 3.7 \times 10^{15}$ m$^{-2}$ and the mobility $\mu = 80$ m$^2$V$^{-1}$s$^{-1}$. The device is composed of two crossing Hall bars \cite{Fujita10E}. The main (horizontal) Hall bar contains six voltage probes, allowing us to measure the transverse voltage $V_{yx}$ at three different locations along the main bar and the longitudinal voltage $V_{xx}$ both on the top and the bottom edges. The secondary (vertical) Hall bar is utilized as a heater to introduce the temperature gradient towards the other end (heat sink) of the main Hall bar. Joule heating by an ac (frequency $f$) heating current $I_\mathrm{h}$ raises the electron temperature in the secondary Hall bar, and the temperature is monitored by the amplitude of the Shubnikov-de Hass (SdH) oscillations. Since the Joule heating varies as $\propto {I_\mathrm{h}}^2$, the resulting thermoelectric voltage can be measured by detecting the component of $V_{yx}$ and $V_{xx}$ having the frequency $2f$, employing standard low-frequency ($f = 13$ Hz) ac lock-in technique. The measurements were carried out in a dilution refrigerator (Oxford, Kelvinox TLD) equipped with a superconducting magnet, with the sample immersed in the mixing chamber held at $T_\mathrm{bath} = 20$ mK\@.

Employing the current $I_h = 200$ nA, we obtain the electron temperature $T_\mathrm{H} = 330$ mK at the heater section. The current is chosen to be much larger than that used in measuring the resistivity (typically 0.5--10 nA), but kept small enough to avoid heating the lattice via the electron-phonon interaction. The gradient is thus introduced only into the electron temperature, leaving the lattice temperature, assumed to be at $T_\mathrm{bath}$, intact. This enables us to selectively measure the diffusion contribution \cite{Maximov04,Chickering09,Fujita10E,Goswami11} and eliminate the phonon-drag contribution. Note that the latter can, under certain circumstances (in a standard experimental method employing an external heater, which introduces the temperature gradient also into the lattice temperature), become very large and dominate the thermoelectric voltage in a 2DES embedded in a thick ($\sim$500 $\mu$m) semiconductor wafer \cite{Fletcher86,Fletcher99}.
Although the phonon-drag contribution is expected to be relatively small in the temperature range of the present study, its elimination by the present experimental technique is advantageous in avoiding possible complication of the interpretation.
Ohmic contacts at the end of the main and the secondary Hall bars and of voltage-probe arms are composed of diffused AuGeNi alloy. The large area (200$\times$200 $\mu$m$^2$) of the contact at the heat sink is devised to ensure good thermal contact to the bath. By contrast, the arms of the voltage probes are designed to be thin and long and are terminated by the contact with a small area (21$\times$21 $\mu$m$^2$), in an attempt to minimize the thermal disturbance to the main bar.
In what follows, however, we assume for simplicity that all the Ohmic contacts, including the small contacts for the voltage probes, are at the same temperature $T_\mathrm{L} = T_\mathrm{bath}$.  The appropriateness of this assumption will be examined in section \ref{discMandC}.

\begin{figure}[tb]
\begin{center}
\includegraphics[bb=40 20 455 606,width=8.3cm,clip]{figure2.png}
\caption{Magnetic-field dependence of the thermoelectric voltages $V_{yx}$ measured between the three pairs of probes, depicted by arrows in the top inset, located across the main Hall bar. (The direction of the arrows indicates the sign of $V_{yx}$.)}
\label{MeasVyx}
\end{center}
\end{figure}

The transverse thermoelectric voltages $V_{yx}$ measured with three pairs of voltage probes are plotted in figure \ref{MeasVyx} as a function of the magnetic field $B$. Quantum oscillations due to the Landau quantization (the equivalent of the SdH oscillations in the resistivity) commences at $B \simeq 0.2$ T\@. The amplitude increases with the increase of the magnetic field up to $B \simeq 0.5$ T, and then levels off and shows slight decrease with further increase of $B$. Above $\sim$0.55 T, the lineshape characteristic of the off-diagonal (Nernst) thermoelectric voltage in the quantum Hall systems  \cite{Jonson84} becomes apparent: namely, plateaus with $V_{yx} \sim 0$ at around integer fillings (only even integer fillings are resolved in this magnetic-field range) and the saw-tooth like lineshape with sign reversal in between (see $S_{yx}$ in figure S2 of the Supplementary data). The traces are antisymmetric with the reversal of the magnetic field. An important feature we want to stress here is that the oscillation amplitude becomes smaller for the voltage-probe pairs located farther away from the heater section.

\begin{figure}[tb]
\begin{center}
\includegraphics[bb=40 165 455 947,width=8.3cm,clip]{figure3.png}
\caption{Magnetic-field dependence of the thermoelectric voltages $V_{xx}$ measured between the two pairs of probes, depicted by arrows in the top inset, located along the main Hall bar.  (The direction of the arrows indicates the sign of $V_{xx}$.) The bottom figure illustrates that the top (4-5) and the bottom (10-9) pairs interchange their roles by inverting the magnetic field.}
\label{MeasVxx}
\end{center}
\end{figure}

Figure \ref{MeasVxx} shows the longitudinal thermoelectric voltage $V_{xx}$ measured with pairs of voltage probes on the top (4-5) and the bottom (10-9) edges. Again, the quantum oscillations are observed above $\sim$0.2 T\@. Interestingly, the lineshape here also exhibits the basic traits of the off-diagonal thermoelectric voltage mentioned above, instead of showing the behaviour expected for the diagonal (Seebeck) thermoelectric voltage, characterized by a dip without sign reversal between two adjacent plateaus $V_{xx} \sim 0$ at integer fillings \cite{Jonson84} (see $S_{xx}$ in figure S2). The variation of the oscillation amplitude with $B$ is similar to that in $V_{yx}$. At $B > 0$, the top voltage-probe pair shows much smaller oscillation amplitude with the sign reversed compared to the bottom pair. By inverting the magnetic field ($B < 0$), however, the top and the bottom pairs switch their roles: the top (bottom) pair gains (loses) the oscillation amplitude and behaves like the mirror image of the trace for the bottom (top) pair in $B > 0$. In the following sections, we will interpret the behaviours of $V_{yx}$ and $V_{xx}$ described above in terms of the calculated spatial distribution of the electron temperature.

\section{Calculation of the electron-temperature maps \label{calcT}}

In a GaAs/AlGaAs 2DES, the electron-phonon interaction is weak at low temperatures relevant to the present study,  as we will see in more detail below. We therefore start by neglecting the electron-phonon interaction. This allows us to obtain analytic solutions for the spatial distribution of the electron temperature $T$ for a one-dimensional (1D) model and for a rectangular sample, as will be briefly delineated in sections \ref{calcT1D} and \ref{calcTRect}, respectively. Numerical calculations resuming the electron-phonon interaction, performed on the Hall bar geometry used in the present study, will be presented in section \ref{calcTFEM}.

Neglecting the electron-phonon interaction, the continuity equation is written as,
\begin{equation}
\boldsymbol{\nabla}  \cdot {\bf{j}}_\mathrm{Q} = 0, \label{conteqwo}
\end{equation}
with ${\bf{j}}_\mathrm{Q} = - \hat{\kappa} \boldsymbol{\nabla} T$ the thermal flux density under the temperature gradient $\boldsymbol{\nabla} T$ and the  thermal conductivity tensor $\hat{\kappa}$. (We use ``hat'' to denote a tensor throughout the paper.) Applying the Wiedemann-Franz law (neglecting the thermal flux carried by phonons), $\hat{\kappa} =  L_0 T \hat{\sigma}$, with $L_0 = \pi ^2 k_\mathrm{B}^2 / (3 e^2) = 2.44 \times 10^{ - 8}$ W$\Omega$/K$^2$ the Lorenz number, $k_\mathrm{B}$ the Boltzmann constant and $\hat{\sigma}$ the conductivity tensor. Defining $\psi \equiv (T^2-{T_\mathrm{L}}^2)/2$, we have $\boldsymbol{\nabla} \psi \equiv T \boldsymbol{\nabla} T$ and
\begin{equation}
{\bf{j}}_\mathrm{Q} = -L_0 \hat{\sigma} \boldsymbol{\nabla} \psi \label{jQ}.
\end{equation}
Noting that the conductivity is virtually independent of the temperature below $\sim$1 K (the temperature range where the mobility is limited by impurity scattering) in a high-mobility GaAs/AlGaAs 2DES \cite{Davies98B}, and assuming a uniform 2DES, we can neglect the spatial derivative of $\hat{\sigma}$ even when placed in the spatially varying temperature. With further assumption that the 2DES is isotropic,
we have $\boldsymbol{\nabla} \cdot (\hat{\sigma} \boldsymbol{\nabla} \psi) = \sigma_{xx} \boldsymbol{\nabla} ^2 \psi$. Thus, the equation to be solved becomes simply a Laplace's equation,
\begin{equation}
\boldsymbol{\nabla} ^2 \psi = 0. \label{Laplaceeq}
\end{equation}

\subsection{One-dimensional model \label{calcT1D}}
First, we consider a simplistic 1D model. In this model, we assume that the electron temperature $T$ is uniform across the  Hall bar ($y$ direction) and $T$ depends only on $x$ (the coordinate along the Hall bar, see figure \ref{HallBarFig} for the $x$-$y$ coordinate). We solve equation (\ref{Laplaceeq}) with the following boundary conditions: $T = T_\mathrm{L}$ at $x=0$ (the end of the Hall bar facing the heat sink) and $T = T_\mathrm{H}$ at $x = L$ (the other end, facing the heater). The solution $T(x)$ can readily be found:
\begin{equation}
T(x) = \left[ \left( {T_\mathrm{H}}^2 - {T_\mathrm{L}}^2 \right) \frac{x}{L} + T_\mathrm{L}^2 \right]^{1/2}. \label{Tx1D}
\end{equation}
From equation (\ref{jQ}), we obtain the thermal flux density,
\begin{equation}
j_{\mathrm{Q}x} = -L_0 \sigma_{xx} \frac{{T_\mathrm{H}}^2 - {T_\mathrm{L}}^2}{2 L}. \label{jQ1D}
\end{equation}

\subsection{Two-dimensional analytic solution for a rectangular sample \label{calcTRect}}
\begin{figure}[tb]
\begin{center}
\includegraphics[width=8.3cm,clip]{figure4.png}
\caption{Boundary conditions for a rectangular sample (left) and the Hall bar geometry employed in the present study (right). The temperatures of the low-temperature and the high-temperature ends of the main Hall bar are fixed to $T_\mathrm{L}$ and $T_\mathrm{H}$, respectively. The temperature of the end of the voltage arms is also fixed to $T_\mathrm{L}$. The thermal flux is not permitted to cross the side edges. Refer to figure \ref{HallBarFig} for the dimensions of the sample. (Note that the right figure is not drawn to scale.)}
\label{boundary}
\end{center}
\end{figure}
Next, we take a look at the rectangular sample corresponding to the main part of the Hall bar with the voltage probes removed. The left panel of figure \ref{boundary} depicts the boundary conditions in this model.
By further defining ${\bm{\tau}} \equiv -\boldsymbol{\nabla} \psi$, we can see that solving equation (\ref{Laplaceeq}) with these boundary conditions
is mathematically equivalent to the process obtaining the electric field, reported by Rendell and Girvin \cite{Rendell81}, if we replace the electric field, the current density and the potential by ${\bm{\tau}}$, ${\bf{j}}_\mathrm{Q}/L_0$ and $\psi$, respectively. We can thus make use of the analytic solution presented in the paper \cite{Rendell81} (we follow the coordinate system used in \cite{Hirayama13}) to have
\numparts
\begin{eqnarray}
\tau_x =  - \tau_0 e^\gamma \cos \vartheta \label{taux}, \\
\tau_y = \tau_0 e^\gamma \sin \vartheta, \label{tauy}
\end{eqnarray}
\endnumparts
and
\numparts
\begin{eqnarray}
j_{\mathrm{Q}x} =  - \tau_0 L_0 \frac{\sigma_{xx}}{\cos \delta} e^\gamma \cos \left( \vartheta  - \delta \right), \label{jQx} \\
j_{\mathrm{Q}y} = \tau_0 L_0 \frac{\sigma_{xx}}{\cos \delta} e^\gamma \sin \left( \vartheta  - \delta \right). \label{jQy}
\end{eqnarray}
\endnumparts
The parameters $\gamma$, $\vartheta$ and $\tau_0$ are given, using the Hall angle $\delta = \arctan (\sigma_{yx} / \sigma_{xx})$, as follows:
\begin{equation}
\gamma  =  - 4 \delta \sum\limits_{n = 1}^\infty \frac{\sinh \left[ {\left( {2n - 1} \right)\pi y/L} \right] \cos \left[ {\left( {2n - 1} \right)\pi x/L} \right]}{{(2n - 1) \pi \cosh \left[ {\left( {2n - 1} \right)\pi \alpha/2} \right]}}, \label{gamma}
\end{equation}
\begin{equation}
\vartheta  = 4\delta \sum\limits_{n = 1}^\infty \frac{{\cosh \left[ {\left( {2n - 1} \right)\pi y/L} \right]\sin \left[ {\left( {2n - 1} \right)\pi x/L} \right]}}{{(2n - 1) \pi \cosh \left[ {\left( {2n - 1} \right)\pi \alpha/2}  \right]}} \label{vartheta}
\end{equation}
and
\begin{equation}
\tau _0 \left( \delta ,\alpha \right) = \frac{{T_\mathrm{H}}^2 - {T_\mathrm{L}}^2}{2 I\left( \delta ,\alpha \right)L}, \label{tau0}
\end{equation}
with
\begin{eqnarray}
I\left( {\delta ,\alpha } \right) \equiv \nonumber \\
 \int_0^1 \!\!\! {\cos \left\{ {4\delta \sum\limits_{n = 1}^\infty  {\frac{{\sin \left[ {\left( {2n - 1} \right)\pi \xi } \right]}}{{\left( {2n - 1} \right)\pi }}} {\mathop{\rm sech}\nolimits} \left[ {\left( {2n - 1} \right)\frac{{\alpha \pi }}{2}} \right]} \right\} d\xi }, \nonumber \\ \label{Iinteg}
\end{eqnarray}
where $\alpha \equiv 2W / L$ represents the aspect ratio. The boundary conditions can readily be confirmed by noting that $\vartheta = 0$ at $x = 0$, $L$ and $\vartheta = 4 \delta \Sigma_{n=1}^\infty \{ \sin[(2n-1)\pi x/L]/[(2n-1)\pi]  \} = \delta$ at $y = -W$, $W$\@.
The electron temperature can be obtained by integrating the ${\bm{\tau}}$ given above:
\begin{equation}
T(x,y) = \left[ {T_\mathrm{L}}^2 - 2 \int_0^x \tau_x(x^\prime,y) dx^\prime \right]^{1/2}. \label{TRect}
\end{equation}
At $B = 0$, $\delta = \gamma = \vartheta =0$ and $I(0,\alpha) = 1$, and thus the temperature is given by equation (\ref{Tx1D}) regardless of the value of $y$.

We can find an approximate formula describing the total thermal flux
\begin{equation}
J_\mathrm{Q} = \int_{-W}^W j_{\mathrm{Q}x} dy \label{JQtot}
\end{equation}
flowing down the rectangle. Noting that $J_\mathrm{Q}$ is conserved along the $x$-direction, we evaluate equation (\ref{JQtot}) at $x = L/2$, where  $\gamma = 0$ and, for not too large $\alpha$, $\vartheta \simeq -4 \delta \Sigma_{n=1}^\infty \{ (-1)^n/[(2n-1)\pi] \} = \delta$. The thermal flux density $j_{\mathrm{Q}x} \simeq -\tau_0 L_0 \sigma_{xx} / \cos \delta$ becomes independent of $y$ and we have
\begin{equation}
J_\mathrm{Q} = j_{\mathrm{Q}x} \cdot 2W \simeq - L_0 \sigma_{xx} \frac{\alpha}{I(\delta,\alpha) \cos \delta} \frac{{T_\mathrm{H}}^2 - {T_\mathrm{L}}^2}{2}. \label{JQRect}
\end{equation}
Since both $I(\delta,\alpha)$ and $\cos \delta$ decrease with increasing $\delta$, $J_\mathrm{Q}$ generally increases with $B$.
At $B = 0$, equation (\ref{JQRect}) can be simply rewritten as
\begin{equation}
J_\mathrm{Q} = j_{\mathrm{Q}x} \cdot 2W \simeq - L_0 \sigma_{xx} \alpha \frac{{T_\mathrm{H}}^2 - {T_\mathrm{L}}^2}{2}. \label{JQRectB0}
\end{equation}
and thus $j_{\mathrm{Q}x}$ coincide with that obtained by the 1D model, equation (\ref{jQ1D}). At a high magnetic field with $\delta$ approaching $\pi / 2$, on the other hand, we can show that $I(\delta,\alpha) \simeq \alpha$ for $\alpha < \sim 0.5$, resulting in
\begin{equation}
J_\mathrm{Q} \simeq - L_0 \frac{\sigma_{xx}}{\cos \delta} \frac{{T_\mathrm{H}}^2 - {T_\mathrm{L}}^2}{2} = L_0 \frac{\sin \delta}{\rho_{yx}} \frac{{T_\mathrm{H}}^2 - {T_\mathrm{L}}^2}{2}. \label{JQRectHighB}
\end{equation}
Note that $J_\mathrm{Q}$ does not depend on the aspect ratio $\alpha$ in equation (\ref{JQRectHighB}).

\subsection{Numerical calculation with electron-phonon interaction for a Hall bar sample \label{calcTFEM}}
\begin{figure*}[tb]
\begin{center}
\includegraphics[bb=15 190 810 690,width=17cm,clip]{figure5.png}
\caption{Spatial distribution of the electron temperature $T$ at various magnetic fields ranging from 0.0 to 1.0 T, calculated taking both the thermal diffusion into the voltage contacts and the electron-phonon interaction into consideration.}
\label{TAephcalc}
\end{center}
\end{figure*}
\begin{figure}[tb]
\begin{center}
\includegraphics[bb=30 0 570 600,width=8.3cm,clip]{figure6.png}
\caption{Spatial distribution of the electron temperature $T$ at $B = 0.0$ T\@. In (a) and (b), the voltage probes are neglected, and $T$ is calculated without (a) and with (b) the electron-phonon interaction ($T_\mathrm{(a)}$ and $T_\mathrm{(b)}$). In (c) and (d), the thermal diffusion into the voltage contacts is considered, and $T$ is calculated without (c) and with (d) the electron-phonon interaction ($T_\mathrm{(c)}$ and $T_\mathrm{(d)}$).}
\label{TB00calc}
\end{center}
\end{figure}
\begin{figure}[tb]
\begin{center}
\includegraphics[width=8.3cm,clip]{figure7.png}
\caption{Calculated electron temperature $T$ along the center line ($y = 0$, $0 \leq x \leq L$) of the main Hall bar at $B = 0.0$ T\@. Black dotted and thin solid lines, overlapping each other, represent the 1D model (equation (\ref{Tx1D})) and the numerical FEM calculation neglecting the electron-phonon interaction and using the rectangular boundary condition ($T_\mathrm{(a)}$), respectively. Thick orange dashed and green solid lines show the numerical FEM calculations for the Hall bar geometry without and with the electron-phonon interaction ($T_\mathrm{(c)}$ and $T_\mathrm{(d)}$), respectively. Gray shades indicate the positions of the voltage probes.}
\label{CSB00}
\end{center}
\end{figure}

\begin{figure}[tb]
\begin{center}
\includegraphics[bb=30 0 570 600,width=8.3cm,clip]{figure8.png}
\caption{Spatial distribution of the electron temperature $T$ at $B = 0.5$ T\@. In (a) and (b), the voltage probes are neglected, and $T$ is calculated without (a) and with (b) the electron-phonon interaction ($T_\mathrm{(a)}$ and $T_\mathrm{(b)}$). In (c) and (d), the thermal diffusion into the voltage contacts is considered, and $T$ is calculated without (c) and with (d) the electron-phonon interaction ($T_\mathrm{(c)}$ and $T_\mathrm{(d)}$). }
\label{TB05calc}
\end{center}
\end{figure}
\begin{figure}[tb]
\begin{center}
\includegraphics[width=8.3cm,clip]{figure9.png}
\caption{Calculated electron temperature $T$ along the center line ($y = 0$, $0 \leq x \leq L$) and close to the bottom ($y = -W + w$, $0 \leq x \leq L$) and the top ($y = W - w$, $0 \leq x \leq L$) edges ($w = 1$ $\mu$m) of the main Hall bar at $B = 0.5$ T, plotted by green, red and blue lines respectively. Dotted and thin solid lines, overlapping each other, represent the analytic solution (equation (\ref{TRect})) and the numerical FEM calculation neglecting the electron-phonon interaction and using the rectangular boundary condition ($T_\mathrm{(a)}$), respectively. Thick dashed and solid lines show the numerical FEM calculations for the Hall bar geometry without and with the electron-phonon interaction ($T_\mathrm{(c)}$ and $T_\mathrm{(d)}$), respectively. Gray shades indicate the positions of the voltage probes.}
\label{CSB05}
\end{center}
\end{figure}

\begin{figure*}[tb]
\begin{center}
\includegraphics[bb=15 260 810 450,width=17cm,clip]{figure10.png}
\caption{Spatial distribution of the temperature decrement due to the thermal diffusion into the voltage contacts, $\Delta T_\mathrm{VP} = T_\mathrm{(a)} - T_\mathrm{(c)}$, where $T_\mathrm{(a)}$ and $T_\mathrm{(c)}$ are the temperatures calculated without and with the thermal diffusion, respectively. The electron-phonon interaction is neglected in both  $T_\mathrm{(a)}$ and $T_\mathrm{(c)}$. }
\label{DTAcalc}
\end{center}
\end{figure*}

\begin{figure*}[tb]
\begin{center}
\includegraphics[bb=15 190 810 690,width=17cm,clip]{figure11.png}
\caption{Spatial distribution of the temperature decrement due to the electron-phonon interaction, $\Delta T_\mathrm{e\mathchar`- ph} = T_\mathrm{(c)} - T_\mathrm{(d)}$, where $T_\mathrm{(c)}$ and $T_\mathrm{(d)}$ are the temperatures calculated without and with the electron-phonon interaction, respectively. The thermal diffusion into the voltage contacts is considered in both $T_\mathrm{(c)}$ and $T_\mathrm{(d)}$. }
\label{DTephcalc}
\end{center}
\end{figure*}
In order to examine the effect of the electron-phonon interaction as well as that of the voltage probes, we have to resort to numerical calculations. Incorporating the electron-phonon interaction, the continuity equation is altered from equation (\ref{conteqwo}) to
\begin{equation}
\boldsymbol{\nabla}  \cdot {\bf{j}}_Q + P_\mathrm{e\mathchar`- ph}(T)  = 0, \label{conteqw}
\end{equation}
where
\begin{equation}
P_\mathrm{e\mathchar`- ph}(T) = P_\mathrm{def}(T) + P_\mathrm{pz}(T)
\end{equation}
represents the power (per area) transferred from the 2DES to the lattice via the electron-phonon interaction and is composed of two components: deformation-potential coupling $P_\mathrm{def} (T)$ and piezo-electric coupling $P_\mathrm{pz} (T)$ contributions \cite{Price82,KajiokaCO13}. They are approximately given by
\begin{equation}
P_\mathrm{def} (T) = \frac{P_\mathrm{D} \left( T^7 - {T_\mathrm{L}}^7 \right)}{{n_\mathrm{e}}^{1/2}} , \label{Pdef}
\end{equation}
and
\begin{equation}
P_\mathrm{pz} (T) = \frac{P_\mathrm{P} \left( T^5 - {T_\mathrm{L}}^5 \right)}{{n_\mathrm{e}}^{1/2}}. \label{Ppz}
\end{equation}
Using the material parameters for GaAs, the coefficients are $P_\mathrm{D} \simeq 5.0 \times10^5$ W$\cdot$K$^{-7}$m$^{-3}$ and $P_\mathrm{P} \simeq 1.8 \times10^6$ W$\cdot$K$^{-5}$m$^{-3}$. Detailed derivation of equations (\ref{Pdef}) and (\ref{Ppz}) is presented in the Appendix. Replacing equation (\ref{jQ}) into equation (\ref{conteqw}), and again assuming the uniform and isotropic $\hat{\sigma}$,
we have
\begin{equation}
-L_0 \sigma_{xx} \boldsymbol{\nabla}^2 \psi + P_\mathrm{e\mathchar`- ph}(T) = 0. \label{epheq}
\end{equation}
We solve equation (\ref{epheq}) in the geometry of the main Hall bar used in the present study (enclosed in the light-green rectangle in figure \ref{HallBarFig}) using the boundary conditions illustrated in the right panel of figure \ref{boundary} with $T_\mathrm{L} = 20$ mK and  $T_\mathrm{H} = 330$ mK.

For the convenience of the numerical calculations, we rewrite equation (\ref{epheq}) in a dimensionless formula, using $L$ and $T_\mathrm{L}$ as units of the length and the temperature, respectively. We denote dimensionless quantities with the tilde: $\widetilde{T} \equiv T / T_\mathrm{L}$, $(\widetilde{x}, \widetilde{y}) \equiv (x/L, y/L)$, $\widetilde{\boldsymbol{\nabla}} = L \boldsymbol{\nabla}$, $\widetilde{\psi}  = ({\widetilde{T}}^2 - 1)/2$ and $\widetilde{\boldsymbol{\nabla}} \widetilde{\psi} = \widetilde{T} \widetilde{\boldsymbol{\nabla}} \widetilde{T} = (L/{T_\mathrm{L}}^2) \boldsymbol{\nabla} \psi$. By further introducing dimensionless coefficients,
\begin{equation}
\gamma _\mathrm{D} = \frac{P_\mathrm{D} L^2 {T_\mathrm{L}}^5}{L_0 \sigma_{xx} {n_e}^{1/2}} \label{gammaD}
\end{equation}
and
\begin{equation}
\gamma _\mathrm{P} = \frac{P_\mathrm{P} L^2 {T_\mathrm{L}}^3}{L_0 \sigma_{xx} {n_e}^{1/2}}, \label{gammaP}
\end{equation}
equation (\ref{epheq}) is rewritten as
\begin{eqnarray}
\widetilde{\boldsymbol{\nabla}} ^2 \widetilde{\psi}  - \gamma_\mathrm{D} \left[ \left( {2 \widetilde{\psi}  + 1} \right)^{7/2} - 1 \right] - \gamma_\mathrm{P} \left[ \left( {2 \widetilde{\psi}  + 1} \right)^{5/2} -1 \right] \nonumber \\ \hspace{60mm} = 0. \label{psieq}
\end{eqnarray}
We numerically solve equation (\ref{psieq}) by the finite element method (FEM) using FreeFem++ \cite{Hecht12}.
Since we are primarily interested in the amplitude of the quantum oscillations in the thermoelectric voltages, we used the upper envelop of the SdH oscillations ($\sigma_\mathrm{UE}$ presented in figure S3 in the Supplementary data) as the $\sigma_{xx}$ to be substituted into equations (\ref{gammaD}) and (\ref{gammaP}).

The resulting maps of the electron temperature
\begin{equation}
T = T_\mathrm{L} (2 \widetilde{\psi} + 1)^{1/2} \label{TFEM}
\end{equation}
are presented in figure \ref{TAephcalc} for the magnetic field $B =0.0$, 0.1, .... 1.0 T\@. At $B = 0.0$ T, the temperature exhibits monotonic decrease from the heater to the heat sink. Once a magnetic field is applied, however, the temperature distribution becomes heavily distorted. This is mainly attributable to the large Hall angle $\delta$
(the plot of $\tan \delta$ vs.\ $B$ used in the present calculation is shown in figure S3 in the Supplementary data). The large $\delta$ along with the boundary condition that the thermal flux $\bf j_\mathrm{Q}$ should not cross the top and the bottom edges of the Hall bar require the temperature gradient $\boldsymbol{\nabla} T$ to have a large angle with the $x$ axis in the main Hall bar (see equation (\ref{jQ})). The spatial distribution of the temperature does not vary noticeably above $\sim$0.3 T, where $\delta$ becomes practically saturated.

To elucidate how the voltage probes and the electron-phonon interaction have affected the temperature maps shown in figure \ref{TAephcalc}, we repeated the numerical FEM calculations eliminating these effects. We calculated the electron temperature $T$ in the following cases: (a) a rectangular sample without the voltage probes and without the electron-phonon interaction, (b) a rectangular sample without the voltage probes but with the electron-phonon interaction, (c) a Hall bar sample with the voltage probes but without electron-phonon interaction, in addition to (d) a Hall bar sample with the voltage probes and with the electron-phonon interaction (shown in figure \ref{TAephcalc}). The results of the calculations $T_\mathrm{(a)}$, $T_\mathrm{(b)}$, $T_\mathrm{(c)}$ and $T_\mathrm{(d)}$ at $B = 0.0$ T and at $B = 0.5$ T are presented in figures \ref{TB00calc}
and \ref{TB05calc},
respectively, and some of their selected cross sections are plotted in figures \ref{CSB00} and \ref{CSB05}, respectively. In order to explicitly visualize the effect of the voltage probes and the electron-phonon interaction, we plot the temperature decrement owing to these effects, $\Delta T_\mathrm{VP} = T_\mathrm{(a)} - T_\mathrm{(c)}$ and $\Delta T_\mathrm{e\mathchar`- ph} = T_\mathrm{(c)} - T_\mathrm{(d)}$, in figures \ref{DTAcalc} and \ref{DTephcalc}, respectively.

At $B = 0.0$ T, $T_\mathrm{(a)}$ reproduces $T(x)$ in the 1D model (equation (\ref{Tx1D})) independent of $y$ (see figures \ref{TB00calc} (a) and \ref{CSB00}). The effects of the voltage probes and the electron-phonon interaction are small and are not readily discernible in figure \ref{TB00calc}. Isothermal lines are seen to bend, albeit very faintly, near the voltage probes in $T_\mathrm{(c)}$ and $T_\mathrm{(d)}$. The cross section at $y = 0$ shown in figure \ref{CSB00} and the temperature decrement plotted in figures \ref{DTAcalc} and \ref{DTephcalc} reveal that voltage probes slightly reduces the temperature ($T_\mathrm{(c)}$), while further inclusion of the electron-phonon interaction barely alters the temperature ($T_\mathrm{(d)}$). The very small disturbance of the temperature distribution by the voltage probes attests to the success in the design of the thin and long arms. The almost negligible effect of the electron-phonon interaction is as expected at low temperatures from the $T$ dependence $P_\mathrm{def} \sim T^7$ and $P_\mathrm{pz} \sim T^5$.

A magnetic field considerably enhances the disturbances, as clearly visualized in figure \ref{CSB05} for an example at $B = 0.5$ T\@. First, we can confirm that $T_\mathrm{(a)}$ replicates the analytic solution $T(x,y)$ for a rectangular sample, equation (\ref{TRect}). A substantial drop of the temperature is induced by the voltage contacts ($T_\mathrm{(c)}$), with the decrement $\Delta T_\mathrm{VP}$ becoming more pronounced as we go closer to the bottom (higher-temperature) edge (see also figure \ref{DTAcalc}). Slight reduction $\Delta T_\mathrm{e\mathchar`- ph}$ is further made by the electron-phonon interaction ($T_\mathrm{(d)}$). Figure \ref{DTephcalc} illustrates that $\Delta T_\mathrm{e\mathchar`- ph}$ is amplified by the increase of the magnetic field.

By applying the analytic rectangular model described in section \ref{calcTRect} to the arms of the voltage probe, we can qualitatively understand why the voltage probes cool down the main bar more effectively in the magnetic field. Without the magnetic field, the thermal flux through the arm $J_\mathrm{Q}$ can be made arbitrarily small by reducing the aspect ratio $\alpha$, as indicated by equation (\ref{JQRectB0}). As mentioned earlier, $J_\mathrm{Q}$ increases with $B$, and the strategy of reducing $\alpha$ ceases to be effective at the magnetic field with $\delta \sim \pi/2$, where $J_\mathrm{Q}$ becomes insensitive to $\alpha$ as implied in equation (\ref{JQRectHighB}).
The increase of $\Delta T_\mathrm{e\mathchar`- ph}$ with the increase of the magnetic field is attributable to the decrease in $\sigma_{xx} = \sigma_\mathrm{UE}$, which enhances the coefficients $\gamma_\mathrm{D}$ and $\gamma_\mathrm{P}$ in equations (\ref{gammaD}) and (\ref{gammaP}). Comparison of figure \ref{DTAcalc} and figure \ref{DTephcalc} reveals that the effect of the voltage probes generally outweighs that of the electron-phonon interaction in reducing the electron temperature in our experimental conditions (note the difference in the colour code scales in the two figures).

\section{Simulating thermoelectric voltages from the calculated spatial distribution of the electron temperature \label{thermvol}}
\begin{figure}[tb]
\begin{center}
\includegraphics[width=8.3cm,clip]{figure12.png}
\caption{Top: $\Delta \widetilde{\psi}_{i,j}$ given by equation (\ref{Dpsiyx}) at $B =0.0$, 0.1, ... 1.0 T calculated using the spatial distribution of $T$
presented in figure \ref{TAephcalc}. Red right-pointing triangles, green squares and blue left-pointing triangles represent $(i,j) = (4,10)$, $(5,9)$ and $(6,8)$, respectively, and lines are the interpolating spline curves. Bottom: Magnetic-field dependence of the thermoelectric voltages between contacts $i$ and $j$, $V_{yx} = \phi_{i,j}$, calculated using equation (\ref{phiijc}) with the (interpolated) $\Delta \widetilde{\psi}_{i,j}$ shown in the top panel and $s_{yx} = S_{yx} / T$ deduced from experimentally obtained $\rho_{xx}(B)$ and $\rho_{yx}(B)$.}
\label{SimVyx}
\end{center}
\end{figure}

\begin{figure}[tb]
\begin{center}
\includegraphics[width=8.3cm,clip]{figure13.png}
\caption{Top: $\Delta \widetilde{\psi}_{i,j}$ given by equation (\ref{Dpsixx}) at $B =0.0$, 0.1, ... 1.0 T calculated using the spatial distribution of $T$
presented in figure \ref{TAephcalc}. Dark-green upright triangles and purple downward triangles represent $(i,j) = (4,5)$ and $(10,9)$, respectively, and lines are the interpolating spline curves. Bottom: Magnetic-field dependence of the thermoelectric voltages between contacts $i$ and $j$, $V_{xx} = \phi_{i,j}$, calculated using equation (\ref{phiijc}) with the (interpolated) $\Delta \widetilde{\psi}_{i,j}$ shown in the top panel and $s_{yx} = S_{yx} / T$ deduced from experimentally obtained $\rho_{xx}(B)$ and $\rho_{yx}(B)$.}
\label{SimVxx}
\end{center}
\end{figure}

In this section, we evaluate the thermoelectric voltage resulting from the gradient $\boldsymbol{\nabla} T$ of the temperature calculated in section \ref{calcT}. We compare the calculated thermoelectric voltages with those measured experimentally, presented in section \ref{meas}, focusing on the amplitude of the quantum oscillations.

Local electric field ${\bf E}$ induced by $\boldsymbol{\nabla} T$ is given by
\begin{equation}
{\bf E} = \hat{S} \boldsymbol{\nabla} T, \label{ESnablaT}
\end{equation}
where $\hat{S}$ represents the thermopower tensor. Since $k_\mathrm{B} T \ll E_\mathrm{F}$ in the temperature range studied in the present paper, we can make use of the generalized Mott's formula (applicable to the 2DES subjected to a perpendicular magnetic field)  \cite{Jonson84},
\begin{equation}
\hat{S} = -e L_0 T \hat{\sigma}^{-1} \left. \frac{\partial \hat{\sigma}}{\partial E} \right|_{E = E_\mathrm{F}} = T \hat{s} \label{Mott}.
\end{equation}
Here we introduced $\hat{s} \equiv \hat{S} / T$  for convenience, whose temperature dependence can be neglected under the assumption of the temperature-independent $\hat{\sigma}$. With this, we have
\begin{equation}
{\bf E} = \hat{s} \boldsymbol{\nabla} \psi. \label{Esnablapsi}
\end{equation}
We can readily see that the temperature distribution in the 1D model, equation (\ref{Tx1D}), results in $\boldsymbol{\nabla} \psi = \partial \psi / \partial x =$ constant. Therefore ${\bf E}$, and accordingly $V_{yx}$, does not vary with $x$ in this model. The diminishing oscillation amplitude with decreasing $x$ exhibited by the quantum oscillations of $V_{yx}$ in figure \ref{MeasVyx} is thus clearly irreconcilable with the purely 1D model.

Using the resistivity tensor $\hat{\rho} = \hat{\sigma}^{-1}$, equation (\ref{Mott}) can be rewritten as,
\numparts
\begin{eqnarray}
s_{xx} = s_{yy} = \frac{S_{xx}}{T} = e L_0 \left. \frac{d}{dE} \ln \sqrt{{\rho _{xx}}^2 + {\rho _{yx}}^2} \right|_{E = E_F},  \label{Srhoxx} \\
s_{yx} = -s_{xy} = \frac{S_{yx}}{T} = e L_0 \left. \frac{d}{dE} \arctan \frac{\rho_{yx}}{\rho_{xx}} \right|_{E = E_F}, \label{Srhoyx}
\end{eqnarray}
\endnumparts
where we made use of the relations $\rho_{xx} = \rho_{yy}$ and $\rho_{xy} = -\rho_{yx}$ fulfilled in an isotropic 2DES\@.
Noting that the quantum oscillations are resulting from the location of the Fermi energy with respect to the Landau levels, we can replace the energy derivative by the derivative with respect to the magnetic field,
\[ \frac{d}{dE} \to  - \frac{B}{E_\mathrm{F}} \frac{d}{dB}, \]
for the assessment of oscillation amplitude. This allows us to deduce the components of $\hat{S}$ from experimentally obtained $\rho_{xx}(B)$ and $\rho_{yx}(B)$ (shown in figure S1 in the Supplementary data), using equations (\ref{Srhoxx}) and (\ref{Srhoyx}). The resulting $S_{xx}$ and $S_{yx}$ are also presented in the Supplementary data (figure S2).

The thermoelectric voltage $\phi_{i,j}$ between the contacts $i$ and $j$, to be compared with the experimentally measured value, is obtained by integrating $-{\bf E}$ along the path $C_{j \rightarrow i}$ connecting the contacts,
\begin{eqnarray}
\phi_{i,j}  =  - \int_{C_{j \rightarrow i}}^{} {\!\!\!\! \bf E}  \cdot d{\bf l }
 = -\left( {\int_{C_{j \rightarrow i}}^{} {\!\!\!\! E_x} dx + \int_{C_{j \rightarrow i}}^{} {\!\!\!\! E_y} dy} \right) \nonumber \\
 = -s_{xx} \int_{C_{j \rightarrow i}}^{} {\!\!\!\! \boldsymbol{\nabla} \psi}  \cdot d{\bf l }
-s_{yx} \left( \int_{C_{j \rightarrow i}}^{} {\!\! \frac{\partial \psi}{\partial x}} dy - \int_{C_{j \rightarrow i}}^{} {\!\! \frac{\partial \psi}{\partial y}} dx \right). \nonumber \\ \label{phiij}
\end{eqnarray}
In our boundary conditions, we assumed that the temperatures of all the voltage contacts are the same as the temperature of the mixing chamber $T_\mathrm{L} = T_\mathrm{bath}$, and thus the first term in equation (\ref{phiij}) vanishes for all the combinations of $i$ and $j$. This indicates that under the present boundary conditions and the assumptions of the temperature-independent isotropic $\hat{\sigma}$, primarily resulting from our experimental conditions that the sample is immersed in the low-temperature $^3$He/$^4$He mixture, both $V_{xx}$ and $V_{yx}$ probe only the off-diagonal (Nernst) component $S_{yx}$ of the thermopower tensor and are insensitive to the diagonal (Seebeck) component $S_{xx}$. Note that the Nernst effect manifests itself in a quite counterintuitive manner owing to the deflection of the temperature gradient described in section \ref{calcT}. We therefore only have to calculate the integration in the second term, which can be rewritten as
\begin{equation}
\phi_{i,j} = -s_{yx} {T_\mathrm{L}}^2 \Delta \widetilde{\psi}_{i,j}, \label{phiijc}
\end{equation}
using a dimensionless integral
\begin{equation}
\Delta \widetilde{\psi}_{i,j}  =   \int_{C_{j \rightarrow i}}^{} {\!\! \frac{\partial \widetilde{\psi}}{\partial \widetilde{x}}} d\widetilde{y} - \int_{C_{j \rightarrow i}}^{} {\!\! \frac{\partial \widetilde{\psi}}{\partial \widetilde{y}}} d\widetilde{x}. \label{Deltapsiij}
\end{equation}
To be more specific, we calculated the path through the centre of the arms,
\begin{equation}
\Delta \widetilde{\psi}_{i,j} = \Delta \widetilde{\psi}_{yx}(\widetilde{W}+\widetilde{L}_\mathrm{V}), \label{Dpsiyx}
\end{equation}
corresponding to $V_{yx}$, and the path passing through the centres of the arms before and after going along the small distance $w = 1$ $\mu$m from the edge of the main bar (see the insets to figure \ref{Integxx}),
\begin{equation}
\Delta \widetilde{\psi}_{i,j} = \Delta \widetilde{\psi}_{xx}(2\widetilde{L}_\mathrm{V}+2\widetilde{w}+\widetilde{L}_3), \label{Dpsixx}
\end{equation}
to be compared with $V_{xx}$.
In equations (\ref{Dpsiyx}) and (\ref{Dpsixx}), we used the integrations along the path, $\Delta \widetilde{\psi}_{yx}(\widetilde{y})$ and $\Delta \widetilde{\psi}_{xx}(\widetilde{\eta})$, defined as (see also figures \ref{HallBarFig}, \ref{Integyx} and \ref{Integxx}),
\begin{eqnarray}
\Delta \widetilde{\psi}_{yx}(\widetilde{y}) \equiv \int_{-\widetilde{W}-\widetilde{L}_\mathrm{V}}^{\widetilde{y}}\left. \frac{\partial \widetilde{\psi}(\widetilde{y}^\prime)}{\partial \widetilde{x}}   \right|_{x = x_n} {\!\!\!\!\! d\widetilde{y}}^\prime \hspace{4mm} (n = 1, 2, 3) \label{psiyx} \\
\hspace{25mm} (-W-L_\mathrm{V} < y < W+L_\mathrm{V}), \nonumber
\end{eqnarray}
where $n =1$, 2 and 3 are for the contact pairs 6-8, 5-9 and 4-10, respectively, and
\begin{eqnarray}
\Delta \widetilde{\psi}_{xx}(\widetilde{\eta}) \equiv \nonumber \\
\left\{ {\begin{array}{*{20}{l}}
\displaystyle \int_{\widetilde{y}_0}^{\widetilde{\eta}_1}  \left. \frac{\partial \widetilde{\psi}}{\partial \widetilde{x}} \right|_{x=x_2} {\!\!\!\!\! d\widetilde{y}}
\hspace{15mm} (0 \leq \eta < L_\mathrm{V}+w) \\
\displaystyle \int_{\widetilde{y}_0}^{\widetilde{y}_1}  \left. \frac{\partial \widetilde{\psi}}{\partial \widetilde{x}} \right|_{x=x_2} {\!\!\!\!\! d\widetilde{y}} -
\int_{\widetilde{x}_2}^{\widetilde{x}_2 + \widetilde{\eta}_2} \left. \frac{\partial \widetilde{\psi}}{\partial \widetilde{y}} \right|_{y = y_1} {\!\!\!\!\! d\widetilde{x}} \\
\hspace{25mm} (L_\mathrm{V}+w \leq \eta < L_\mathrm{V}+w+L_3) \\
\displaystyle \int_{\widetilde{y}_0}^{\widetilde{y}_1}  \left. \frac{\partial \widetilde{\psi}}{\partial \widetilde{x}} \right|_{x=x_2} {\!\!\!\!\! d\widetilde{y}}  -
\int_{\widetilde{x}_2}^{\widetilde{x}_3} \left. \frac{\partial \widetilde{\psi}}{\partial \widetilde{y}} \right|_{y = y_1} {\!\!\!\!\! d\widetilde{x}}  +
\int_{\widetilde{y}_1}^{\widetilde{\eta}_3}  \left. \frac{\partial \widetilde{\psi}}{\partial \widetilde{x}} \right|_{x=x_3} {\!\!\!\!\! d\widetilde{y}}, \\
\hspace{20mm} (L_\mathrm{V}+w+L_3  \leq \eta < 2 L_\mathrm{V}+2 w+L_3)
\end{array}} \right. \label{psixx}
\end{eqnarray}
with $x_1 \equiv L_1$, $x_2 \equiv L_1 + L_2$, $x_3 \equiv L_1 + L_2 + L_3$, $y_0 \equiv \pm (W+L_\mathrm{V})$, $y_1 \equiv \pm (W - w)$,
$\eta_1 \equiv y_0 \mp \eta$, $\eta_2 \equiv \eta - (L_\mathrm{V} + w)$, $\eta_3 \equiv y_1 \pm (\eta - L_\mathrm{V} -  w - L_3)$, where the upper (lower) sign in $y_0$, $y_1$, $\eta_1$ and $\eta_3$ is for the contact pair 4-5 (9-10).

In the top panels of figures \ref{SimVyx} and \ref{SimVxx}, we plot $\Delta \widetilde{\psi}_{i,j}$ calculated with equations (\ref{Dpsiyx}) and (\ref{Dpsixx}), respectively, at $B = 0.0$, 0.1,..., 1.0 T for the contact pairs $(i,j)$ noted in the figure. The values of $\Delta \widetilde{\psi}_{i,j}$ for $(i,j) = (4,10), (5,9), (6,8), (4,5)$ and $-\Delta \widetilde{\psi}_{10,9}$ exhibit similar $B$-dependence: after initial steep rise from $B = 0.0$ T, they take a mild peak at $B \sim 0.4$ T, followed by a gentle decline. The behaviour basically tracks the $B$-dependence of the Hall angle $\delta$ (see figure S3 in the Supplementary data). This can be qualitatively understood by noting that $\delta$ plays a key role in determining the orientation of
$\boldsymbol{\nabla} T$. The angle between $\boldsymbol{\nabla} T$ and the side edges generally increases with $\delta$, leading to the enhancement of $|\partial \widetilde{\psi} / \partial \widetilde{x}|$ in the arms and $|\partial \widetilde{\psi} / \partial \widetilde{y}|$ in the main bar, both resulting in the increase of $|\Delta \widetilde{\psi}_{i,j}|$. We note in passing that $\delta$ shows non-monotonic $B$-dependence because we employed $\sigma_\mathrm{UE}$ instead of semiclassical $\sigma_{xx} = n_\mathrm{e} e \mu / (1 + \mu^2 B^2)$ in calculating $\delta$.

We can see the relations $\Delta \widetilde{\psi}_{6,8} < \Delta \widetilde{\psi}_{5,9} < \Delta \widetilde{\psi}_{4,10}$ and  $\Delta \widetilde{\psi}_{10,9} < 0 < \Delta \widetilde{\psi}_{4,5}$ with $|\Delta \widetilde{\psi}_{10,9}| \gg |\Delta \widetilde{\psi}_{4,5}|$ in figures \ref{SimVyx}  and \ref{SimVxx}, respectively. Using an interpolating spline curve connecting the data points, we calculate $\phi_{i,j}$ with equation  (\ref{phiijc}) to simulate the experimentally observed $V_{yx}$ and $V_{xx}$. The calculated $\phi_{i,j}$, plotted in the bottom panels of figures \ref{SimVyx} and \ref{SimVxx}, qualitatively reproduce essential features in the experimental traces shown in figures \ref{MeasVyx} and \ref{MeasVxx}, respectively. The variation of  the oscillation amplitude with the magnetic field is in agreement with the experiments in both figures. Diminishing amplitude with the distance
from the heater seen in $V_{yx}$ and the relation of the sign (phase) and the amplitude between the top and the bottom voltage contacts in $V_{xx}$ are also well reproduced. All of these traits inherit the behaviour of the  $\Delta \widetilde{\psi}_{i,j}$ mentioned above. By inverting the magnetic field, the sign of $\delta$, hence the orientation of $\boldsymbol{\nabla} T$ with respect to ${\bf j}_\mathrm{Q}$, is also inverted. The spatial distribution of the temperature can thus be obtained by simply mirroring those shown in figure \ref{TAephcalc} with respect to the $x$ axis ($y = 0$). Therefore the voltage between the top contacts (4-5) and the bottom contacts (10-9) simply change places, namely, the top contacts at $B < 0$ take the place of the bottom contacts at $B > 0$, and vice versa. This is basically in agreement with what we observe in the experiments (figure \ref{MeasVxx}).
\begin{figure}[tb]
\begin{center}
\includegraphics[width=8.3cm,clip]{figure14.png}
\caption{$\Delta \widetilde{\psi}_{yx}(\widetilde{y})$ given by equation (\ref{psiyx}) for contact pairs 6-8, 5-9 and 4-10 at $B =0.0$ T (top) and $B = 0.5$ T (bottom). }
\label{Integyx}
\end{center}
\end{figure}
\begin{figure}[tb]
\begin{center}
\includegraphics[width=8.3cm,clip]{figure15.png}
\caption{$\Delta \widetilde{\psi}_{xx}(\widetilde{\eta})$ given by equation (\ref{psixx}) for contact pairs 4-5 (top) and 10-9 (bottom) at $B = 0.5$ T\@. Insets: Schematic drawings depicting the integration path.}
\label{Integxx}
\end{center}
\end{figure}

To see which part of the device has made the dominant contribution to the thermoelectric voltages, we examine how the integral in equation (\ref{Deltapsiij}) varies during the course of the integration. In figure \ref{Integyx}, we plot $\Delta \widetilde{\psi}_{yx}(\widetilde{y})$ given by equation (\ref{psiyx}) as a function of $\widetilde{y}$ at $B = 0.0$ T and 0.5 T\@. Without the magnetic field, the contribution of the voltage arms is negligibly small and $\Delta \widetilde{\psi}_{i,j}$ almost exclusively derives from the main bar. Note, however, that the $\Delta \widetilde{\psi}_{i,j}$ do not yield a thermoelectric voltage, since $s_{yx} = 0$ at $B = 0.0$ T\@. In the magnetic field, by contrast, the voltage arms accounts for most of the integral. $\Delta \widetilde{\psi}_{i,j}$ $(>0)$ is primarily determined by the gain at the bottom voltage arm, which far exceeds the loss at the top voltage arm. The integral over the main bar is negligibly small since $\boldsymbol{\nabla} \psi$ is nearly parallel to $y$ direction there. Again by applying the analytic rectangular model to the arms, we can readily find that the integral over an arm increases with the temperature difference between the contact (at $T_\mathrm{L}$) and the edge of the main bar to which the arm is attached. Note that the gradient $|\partial \psi / \partial x|$ (corresponding to $\tau_y$ in equation (\ref{tauy})) is enhanced by the temperature difference through $\tau_0$ in equation (\ref{tau0}). The bottom arm dominates $\Delta \widetilde{\psi}_{i,j}$, since it is attached to the higher-temperature edge of the main bar.  The decrease of $\Delta \widetilde{\psi}_{i,j}$ with the distance from the heater traces the decrease in the temperature of the bottom edge, which, in turn, mainly results from the thermal diffusion through the arms into the contacts (see the thick solid red line in figure \ref{CSB05}).

Figure \ref{Integxx} shows $\Delta \widetilde{\psi}_{xx}(\widetilde{\eta})$ given by equation (\ref{psixx}) as a function of $\widetilde{\eta}$ at $B = 0.5$ T for the top and the bottom contact pairs. Again, the voltage arms dominate the integrals, although the main bar (the second term in equation (\ref{psixx}) containing $|\partial \psi / \partial y|$) also makes a discernible contribution.
For the top pairs, the voltage arms are attached to the locations of the main bar having nearly the same temperature (see the thick solid blue line in figure \ref{CSB05}). Therefore the contributions of the two arms cancel each other, and the relatively small integral from the main bar survives, leading to a small value of $\Delta \widetilde{\psi}_{i,j} > 0$. For the bottom pairs, contribution of the arms is much larger than the top pairs. Near the bottom edge of the main bar, the locations for the two voltage arms have a relatively large temperature difference, as can be perceived from the thick solid red line in figure \ref{CSB05}. Thus, the loss due to the arm of the voltage probe 10 far exceeds the combined gain by the voltage probe 9 and the main bar, resulting in a large negative value of $\Delta \widetilde{\psi}_{i,j}$.

\section{Discussion \label{discussion}}
\subsection{Possible sources of the discrepancy between measured and calculated thermoelectric voltages \label{discMandC}}
So far, we have shown that our calculations qualitatively reproduce the essential features of our measurements. Quantitatively, however, the calculated thermoelectric voltages are larger than the measured values roughly by a factor of 2.5 to 10. We presume that the discrepancy is mainly caused by our use of oversimplified boundary conditions to make the calculations tractable. First, we simply assumed that the higher-temperature end of the main Hall bar ($x = L$ and $-W/2 < y < W/2$) has a constant temperature $T_\mathrm{H}$ inferred from the SdH amplitude measured with the voltage probes located on the other side of the secondary (heater) Hall bar (see figure \ref{HallBarFig}). We have neglected the spatial variation of the temperature within the heater section, which can be especially pronounced when the magnetic field is applied and a hot spot is generated \cite{Hirayama13}. Second, we set all the Ohmic contacts at the temperature $T_\mathrm{L} = T_\mathrm{bath}$ of the mixing chamber in which the sample is immersed. The contacts can, however, have the temperature higher than the bath owing to the Kapitza resistance \cite{Lounasmaa74} between the helium and the metallic (AuGeNi) film constituting the contacts. The effect will be more apparent for the voltage contacts having a smaller area hence a higher thermal resistance, especially for those attached to the higher temperature edge of the main Hall bar (contacts 10, 9, 8) and thus carrying higher thermal flux. The possible increase of the temperature in these contacts can, in principle, reduce the resulting thermoelectric voltages. It is difficult, however, to quantitatively assess the decrement mainly owing to the difficulty in determining the area of the interface between the helium and the granular metallic film.
Note, however, that the temperature difference between voltage contacts, if present, allows for the longitudinal component, $S_{xx}$, to contribute to the thermoelectric voltages by resuming the first term in equation (\ref{phiij}). As pointed out earlier, the footprint of $S_{xx}$ is barely discernible in the lineshape of the experimentally observed oscillations shown in figures \ref{MeasVyx} and \ref{MeasVxx}, suggesting the absence of appreciable temperature difference between the measured contact pairs. We therefore surmise that the temperatures of the voltage contacts are sufficiently close to $T_\mathrm{L}$.

We can also find possible sources of the quantitative disagreement between the calculations and the measurements apart from the boundary conditions. We tacitly assumed that the 2DES in the voltage arms have the same properties as those of the main bar. It is possible, however, that the quality of the 2DES has declined during the process of fabricating the narrow (3 $\mu$m) arms. This can alter the temperature distribution by reducing the amplitude of SdH oscillations hence $\sigma_{xx} = \sigma_\mathrm{UE}$, which determines $\delta$, $\gamma_\mathrm{D}$ and $\gamma_\mathrm{P}$. Probably more importantly, the amplitude of the oscillations in $s_{yx}$ can also be reduced, resulting in the reduction of the oscillation amplitude of the thermoelectric voltage generated at the arms. Neglecting the quantum Hall edge states can also be a source, especially in the higher magnetic field regime. All the effects described here are rather difficult to be unambiguously quantified and cannot readily be incorporated in our model. Despite the quantitative discrepancy, we believe that our simple model captures the essence of the spatial distribution of the temperatures and the resulting thermoelectric voltages in a Hall bar placed in a magnetic field.

\subsection{Characteristics and limitations of the present measurement configurations \label{discCharLim}}
As we have shown in the previous sections, thermoelectric voltages can be a good tool to probe the spatial distribution of $T$ with the aid of the known thermoelectric coefficients. To this end, we used thermoelectric coefficients deduced from the measured resistivity employing the generalized Mott's formula. In many cases, however, the purpose of measuring thermoelectric voltages is to extract unknown thermoelectric coefficients. Our present settings are obviously inappropriate for this purpose. As mentioned in section \ref{thermvol}, the measurement is insensitive to $S_{xx}$. Furthermore, voltages generated at the voltage-probe arms exceed those generated at the main bar, the presupposed arena for the measurements. The situation cannot be improved by making the arms shorter, since the effect of the shortening the integration path in equation (\ref{Deltapsiij}) is cancelled out by the increase in the gradient (see equation (\ref{tau0})). These drawbacks mainly stem from our experimental configuration in which the voltage contact pads are thermally in contact with the surrounding helium bath. In a standard experimental setup for thermoelectric measurements, a sample is placed in a vacuum and is cooled down by thermally connecting one end (the heat sink) to the refrigerant, and the temperature gradient is introduced by heating the other end \cite{Fletcher99,White79}. By selecting a material with a high thermal resistance for the wiring to the voltage contacts \cite{Fletcher86}, the contribution of the voltage arms can be reduced. It will still be difficult, however, to completely eliminate or precisely determine the temperature gradient at the voltage arms and the resulting thermoelectric voltages. Regardless of the contribution of the arms, the redistribution of the temperature within the main bar by the magnetic field should be taken into consideration in interpreting the measured thermoelectric voltages.  It is also worth pointing out that the electron temperature of GaAs/AlGaAs 2DES is quite difficult to cool down to the lowest temperature in a dilution refrigerator in this standard configuration (see, e.g., \cite{Samkharadze11}).

The magnetic field does not affect the thermal flux carried by phonons. Therefore, the effect of the magnetic field on the spatial distribution of the temperature becomes less significant in the measurement
employing an external heater to introduce temperature gradient also to the substrate, provided that the temperature is high enough for the thermal flux to be principally carried by phonons, and also for the electron-phonon coupling to be strong enough so that the electron system has the same temperature as the lattice system. In such circumstances, the thermoelectric power is dominated by the phonon-drag contribution \cite{Fletcher86}.

\section{Conclusions \label{conclusion}}
We have examined the effect of a magnetic field on the spatial distribution of the electron temperature and the resulting thermoelectric voltages generated in a 2DES residing in a Hall bar device fabricated from a GaAs/AlGaAs wafer. At the temperatures investigated in the present study ($T < T_\mathrm{H} =330$ mK), electrons and phonons are virtually decoupled, and the relations between the gradient $\boldsymbol{\nabla} T$  of the electron temperature (or more precisely, ${\bm{\tau}} = -\boldsymbol{\nabla} \psi$), the thermal flux density ${\bf{j}}_\mathrm{Q}$ and the source-drain temperature difference (or, difference in $\psi$) basically duplicate those between the electric field, the electric current density and the source-drain bias in the magnetic field. The temperature difference engendered between the top and the bottom edges of the main bar, resulting from the deflection of $\boldsymbol{\nabla} T$ from the source-drain direction by a large angle $\delta \sim \pi/2$, is most directly detected by the difference in the thermoelectric voltages $V_{xx}$ measured along the top and the bottom edges (figure \ref{MeasVxx}). The top and the bottom edges exchange their roles by inverting the magnetic field, owing to the sign reversal in $\delta$. A notable feature unique to the thermoelectric-voltage measurements is the significant roles played by the voltage arms. A magnetic field lets the thermal flux flow through the arms regardless of the thickness and the length, disabling the design of the low aspect ratio to suppress the thermal flux. The resultant drop in the temperature, as well as the substantial thermoelectric voltage generated within the arms themselves, is responsible for the decrease in the measured thermoelectric voltages $V_{yx}$ with the distance from the heater (figure \ref{MeasVyx}). The effects of the arms also enhance the difference in $V_{xx}$ between the top and the bottom edges. The redistribution and reorientation of the temperature by a magnetic field demonstrated in the present study pose a caveat to be borne in mind in interpreting the thermoelectric-voltage of a 2DES measured in a magnetic field at very low temperatures.

\ack
We acknowledge T.\ Okamoto for calling our attention to the problem with spatially varying $V_{yx}$. This work was supported by JSPS KAKENHI Grant Numbers JP26400311, JP17K05491.


\section{Introduction}

Powerful deep neural networks (DNNs) have been created and investigated for high-level computer vision tasks such as image classification~\cite{krizhevsky2012imagenet,he2016deep}, object detection~\cite{ren2015faster,redmon2017yolo9000}, and semantic segmentation~\cite{longfully} as well as for low-level computer vision tasks such as image denoising~\cite{Vincent:2010vu,burger2012image,Zhang:2017eh,Lefkimmiatis:2017ka}.
Initially, it was challenging for DNNs to outperform powerful classical denoisers such as BM3D~\cite{burger2012image}.
However, recent works with DNNs proposed and demonstrated that it is possible for DNNs to outperform classical denoisers for synthetic Gaussian noise~\cite{Zhang:2017eh} as well as for real noise~\cite{Brooks:2019vh}. All aforementioned DNN based denoisers were trained in a supervised way with noiseless ground truth images.

Collecting high-quality noiseless images for training and evaluating DNN denoisers is challenging. Plotz and Roth collected high-quality benchmark data for denoising by averaging 19 independent noise realizations per one image~\cite{plotz2017benchmarking}. It is a painstaking process to take tens of photos for one high-resolution image of static objects and to perform post-processing to compensate for lighting changes. It seems even more challenging, time-consuming, and even expensive to collect noiseless high-quality ground truth data for slowly moving objects ($e.g.$, animals, humans), for medical imaging, and for airborne hyper-spectral imaging. Even though it is possible to take 19 pictures, it may be inevitable that such a ground truth image contains mild noise if each picture is relatively noisy. For example, an average of 19 pictures contaminated with Gaussian noise of $\sigma=25$ yields a picture with noise level of $\sigma=5.74$ assuming temporal independence of noise.
Thus, it seems desirable to have methods to deal with imperfect ground truth data ($i.e.$, mild noise in ground truth) and/or to train DNN denoisers without clean ground truth.

Recently, there have been several works on unsupervised training of DNN denoisers with noisy images only.
Deep image prior (DIP) exploited the structure of a generator network and minimized the mean-squared error (MSE)
between the output of a DNN and a given noisy image to denoise~\cite{ulyanov2018deep}. While DIP does not train the DNN, it requires to compute MSE minimization for each noisy image, which is slower than other DNN denoisers.
Noise2Noise was proposed to train DNNs for image restoration with the set of two (independent) noise realizations per image in an unsupervised way for various noise models including Gaussian distribution and a wide range of applications including compressive sensing MR recovery~\cite{noise2noise}.
Stein’s unbiased risk estimator (SURE) based training method for Gaussian denoisers was proposed to train DNNs with a set of a single noise realization per image, but it was limited to Gaussian denoising~\cite{soltanayev2018training}.

Both Noise2Noise and SURE based method outperformed classical denoising methods such as BM3D for synthetic Gaussian noise and they often yielded comparable performance to DNN denoisers that were trained with clean ground truth.
Noise2Noise has demonstrated powerful performance in various image denoising and restoration tasks for zero-mean contaminations~\cite{noise2noise}.
However, it required two independent noise realizations per image empirically, while there was no theoretical explanation on
the relationship between two noise realizations.
Thus, it is not clear if Noise2Noise can be used for the case with a single noise realization or for the case with imperfect ground truth data.
Moreover, assuming slowly moving objects or slowly varying light conditions, there is a trade-off between
low noise in ground truth and identical underlying true image over multiple realizations.

Even though SURE based training method is limited to Gaussian noise~\cite{soltanayev2018training}, it has a potential to be extended to more general noise models such as mixed Poisson-Gaussian model~\cite{LeMontagner:2014jd}, exponential family~\cite{Eldar:2009fga} or non-parametric model~\cite{Raphan:2007wj}. It is also extended to unsupervised learning in inverse problems~\cite{Metzler:2019gsure}.
SURE could also be more robust to the noise-blur trade-off than Noise2Noise by using a single noise realization per image for slowly moving objects.
However, SURE based method is also limited since it does not take advantage of having multiple noise realizations per image
when they are available ($e.g.$, two uncorrelated noise realizations per image for Noise2Noise, imperfect ground truth image with mild noise).
In this paper, we address the following questions: 1) can Noise2Noise deal with correlated noise realizations and imperfect ground truth? 2) can SURE be extended to take advantage of having two uncorrelated noise realizations per image as in Noise2Noise? 3) can the extended SURE handle correlated noisy images and well utilize imperfect ground truth?

Here, we propose eSURE to training deep Gaussian denoisers with correlated pairs of noise realizations per image and applied it to
the case with two uncorrelated realizations per image to achieve better performance than the original SURE based method and comparable results to Noise2Noise.
Then, we further investigated the case of training Gaussian denoisers with imperfect ground truth ($i.e.$, mild noise in ground truth) by adding synthetic noise.

For the case of adding noise to imperfect ground truth to yield correlated pairs of images, our proposed eSURE based training method outperformed
conventional SURE based method as well as Noise2Noise.

Here is the summary on the contributions of this paper:
\begin{itemize}
\item Analyzing Noise2Noise theoretically and empirically for correlated pairs of noise realizations per image.
\item Extending SURE to take advantage of having a pair of noise realizations per image for training unlike the conventional SURE that can take a single noise realization per image.
\item Investigating eSURE that exploits two independent or correlated noise realizations per image as well as that can utilize imperfect ground truth with mild noise.
\end{itemize}
We will show that there is a clear theoretical link between Noise2Noise and eSURE in a limited case of Gaussian denoising.
In fact, it turned out that Noise2Noise is a special case of eSURE for independent Gaussian noise in theory.
However, while the performance of Noise2Noise was degraded with correlated pairs of noisy images, the performance of the proposed eSURE remains the same.

This paper is organized as follows: Section~\ref{theory} briefly reviews SURE, Monte-Carlo SURE (MC-SURE) and SURE based denoiser training.
Then, Section~\ref{method} revisits theoretical derivation of Noise2Noise training methods,
 proposes an extended version of MC-SURE (called eSURE) and shows a clear link between our eSURE and Noise2Noise.
 These are followed by experimental results to validate the effectiveness of our proposed unsupervised DNN training method in Section~\ref{results}
 for the case of having two independent realizations per image and for the case with imperfect ground truth.
Finally, Section~\ref{conclusion} draws a conclusion of this work.

\section{Background} \label{theory}
\subsection{Stein`s Unbiased Risk Estimator (SURE)}
Typically, Gaussian contaminated signal (or image) is modeled as a linear equation:
\begin{equation}
\mathbf{y} = \mathbf{x} + \mathbf{n}
\label{eq:signal}
\end{equation}
where $\mathbf{x} \in \mathbb{R}^N$ is an unknown signal,
$\mathbf{y} \in \mathbb{R}^N$ is a known measurement,
$\mathbf{n} \in \mathbb{R}^N$ is an \textit{i.i.d.} Gaussian noise
such that $\mathbf{n} \sim \mathcal{N} (\mathbf{0}, \sigma^2 \mathbf{I})$, and
$\mathbf{I}$ is an identity matrix. We denote $\mathbf{n} \sim \mathcal{N} (\mathbf{0}, \sigma^2 \mathbf{I})$
as $\mathbf{n} \sim \mathcal{N}_{0, \sigma^2}$.

In general, given an estmator $\mathbf{h}( \mathbf{y} )$ of $\mathbf{x}$, the SURE has the following form:
\begin{eqnarray} \label{eq:sure}
\eta (\mathbf{h}( \mathbf{y} )) =  \frac{\| \mathbf{y} - \mathbf{h} ( \mathbf{y} ) \|^2}{N}  - \sigma^2
+ \frac{2 \sigma^2}{N} \sum_{i=1}^N
\frac{\partial \mathbf{h}_i ( \mathbf{y} )}{\partial \mathbf{y}_i }
\end{eqnarray}
Assuming $\mathbf{x}$ to be deterministic signal (or image),
the following theorem for (\ref{eq:sure})  holds.
\begin{thm}\cite{Stein:1981vf, Blu:2007fz}
	The random variable $\eta (\mathbf{h}( \mathbf{y} ))$ is an unbiased estimator of
	\[
	\mathrm{MSE}(\mathbf{h}( \mathbf{y} )) = \frac{1}{N} \| \mathbf{x} - \mathbf{h}( \mathbf{y} ) \|^2
	\]
	or
	\begin{equation}
	\mathbb{E}_{\mathbf{n} \sim \mathcal{N}_{0, \sigma^2}}
	\left\{ \frac{\| \mathbf{x} - \mathbf{h}( \mathbf{y} ) \|^2}{N}  \right\}
	= \mathbb{E}_{\mathbf{n} \sim \mathcal{N}_{0, \sigma^2}}
	\left\{  \eta (\mathbf{h}( \mathbf{y} )) \right\}
	\label{eq:expectsure}
	\end{equation}
	\label{thm:sure}
	\vspace{-1em}
\end{thm}
where $\mathbb{E}_{\mathbf{n} \sim \mathcal{N}_{0, \sigma^2}}\{ \cdot \}$ is the expectation operator in terms of the random vector $\mathbf{n}$.

Although (\ref{eq:sure}) looks appealing in terms of optimizing parameters of an estimator $\mathbf{h}( \mathbf{y})$, the analytical solution for the last divergence term in  (\ref{eq:sure}) is limited only to some special cases such as the estimator $\mathbf{h}( \mathbf{y} )$ to be
non-local mean or linear filters \cite{van2009sure,van2011nonlocal}. Thus, in order to utilize (\ref{eq:sure}),
one needs to find at least an approximate solution of the divergence term for more general cases.

\subsection{Monte-Carlo SURE (MC-SURE)}
A fast Monte-Carlo approximation of the divergence term has been developed by Ramani \textit{et al.} in \cite{Ramani:2008ij}. This method
yielded accurate unbiased estimate of MSE for many denoising methods $\mathbf{h} ( \mathbf{y} )$.
\begin{thm} \cite{Ramani:2008ij} Let $\mathbf{\tilde{n}} \sim \mathcal{N}_{0, 1} \in \mathbb{R}^N$
	be independent of $\mathbf{n}$ or $\mathbf{y}$. Then,
	\begin{equation}
	\sum_{i=1}^K \frac{\partial \mathbf{h}_i ( \mathbf{y} )}{\partial \mathbf{y}_i } =
	\lim_{\epsilon \to 0} \mathbb{E}_{\mathbf{\tilde{n}}} \left\{ \mathbf{\tilde{n}}^\mathrm{t}
	\left(  \frac{\mathbf{h} ( \mathbf{y} + \epsilon \mathbf{\tilde{n}} ) - \mathbf{h} ( \mathbf{y} )}{\epsilon} \right)\right\}
	\end{equation}
	provided that $\mathbf{h} ( \mathbf{y} )$ admits a well-defined second-order Taylor expansion.
	If not, this is still valid in the weak sense provided that $\mathbf{h} ( \mathbf{y} )$ is tempered.
	\label{thm:mcsure}
\end{thm}

Consequently, by applying Theorem~\ref{thm:mcsure} to the divergence term in (\ref{eq:sure}), the divergence approximation of the denoiser $\mathbf{h} ( \mathbf{y} )$
will be:
\begin{equation}
\frac{1}{N} \sum_{i=1}^N \frac{\partial \mathbf{h}_i ( \mathbf{y} )}{\partial \mathbf{y}_i } \approx
\frac{1}{\epsilon N} \mathbf{\tilde{n}}^\mathrm{T}
\left(  \mathbf{h} ( \mathbf{y} + \epsilon \mathbf{\tilde{n}} ) - \mathbf{h} ( \mathbf{y} ) \right),
\label{eq:approx}
\end{equation}
where $\mathbf{\tilde{n}}^\mathrm{T}$ is a transposed $i.i.d$ Gaussian vector ($\mathbf{\tilde{n}} \sim \mathcal{N}_{0, 1}$) and $\epsilon$ is a fixed small positive value.

\subsection{SURE based deep denoiser training}
Recently, SURE was used as a surrogate metric to minimize the MSE between the output of the DNNs and the ground truth for unsupervised training of DNN based Gaussian  denoisers~\cite{soltanayev2018training}. Specifically, MC-SURE allows DNN to learn large-scale weights by minimizing MC-SURE with no noiseless ground truth images for Gaussian denoising. The equation (2) with (4) was reformulated for the DNN $\mathbf{h_\theta} ( \cdot )$ as follows:
\begin{eqnarray}
\eta (\boldsymbol{h_\theta}( \mathbf{y} )) = \frac{1}{M} \sum_{j=1}^M \bigg\{
\| \mathbf{y}^{(j)} - \mathbf{h_\theta} ( \mathbf{y}^{(j)} ) \|^2 -N \sigma^2 \label{eq:riskdenoiser_emp_batch_prop_mcsure} \\
+ \frac{2 \sigma^2}{\epsilon} (\mathbf{\tilde{n}}^{(j)})^\mathrm{t}
\left(  \boldsymbol{h_\theta} ( \mathbf{y}^{(j)} + \epsilon \mathbf{\tilde{n}}^{(j)}) - \boldsymbol{h_\theta} ( \mathbf{y}^{(j)}) \right) \bigg\},
\nonumber
\end{eqnarray}
where $\boldsymbol{\theta}$ is the set of DNN denoiser parameters, $M$ is the size of mini-batch, $\epsilon$ is a small fixed positive constant,
and $\mathbf{\tilde{n}}^{(j)}$ is a single realization from standard normal distribution for each training data $j$.

This approach has been demonstrated to yield state-of-the-art performance in denoising task with synthetic Gaussian noise,
yielding comparable to or slightly worse qualitative and quantitative results than
MSE-trained DNNs.

 \section{Methods} \label{method}
In this section, we first re-visit Noise2Noise method~\cite{noise2noise} and re-derive the method of Noise2Noise in a different approach from~\cite{noise2noise}.
Then, we propose to extend the original SURE and MC-SURE to deal with a pair of correlated noisy images instead of a single noisy image
and to use it for training deep learning based denoisers
with pairs of correlated Gaussian noise realizations per image.
We also show that Noise2Noise is a special case of our proposed extended SURE for Gaussian denoising.
Lastly, we will demonstrate that our proposed method is more robust to correlated noise realization pairs in training data set compared to a Noise2Noise method.
Our eSURE is especially useful for the case of using imperfect ground truth images with mild noise.

\subsection{Revisiting Noise2Noise}

The Noise2Noise method has been proposed to train DNNs for image processing only with noisy images where two noise realizations per image were required~\cite{noise2noise}.
Its theoretical justification required zero-mean noise, but there was no clear assumption on independence or uncorrelated property of two realizations.
However, two independent noise realizations per image were used empirically.

Assuming that the triplet $(\mathbf{x}, \mathbf{y}, \mathbf{z})$ follows a joint distribution and the expectation of two random vectors $\mathbf{y} - \mathbf{x}$, $\mathbf{z} - \mathbf{x}$ are both zero vectors, the MSE for infinite data is as follows:
\begin{eqnarray}
\mathbb{E}_{(\mathbf{x}, \mathbf{y})}
\left\{ \| \mathbf{x} - \boldsymbol{h_\theta}( \mathbf{y} ) \|^2 \right\}  &=&
\mathbb{E}_{\mathbf{x}}  \left[ \mathbb{E}_{(\mathbf{y}, \mathbf{z}) | \mathbf{x} }
\left\{ \| \mathbf{x} - \mathbf{z} + \mathbf{z} - \boldsymbol{h_\theta}( \mathbf{y} ) \|^2 | \mathbf{x} \right\}  \right]  \label{eq:n2nderiv} \\
&=& \mathbb{E}_{\mathbf{x}}  \left[ \mathbb{E}_{(\mathbf{y}, \mathbf{z}) | \mathbf{x}}
\left\{ \| \mathbf{z} - \boldsymbol{h_\theta}( \mathbf{y} ) \|^2 + 2  (\mathbf{z} - \mathbf{x})^T \boldsymbol{h_\theta}( \mathbf{y} ) | \mathbf{x} \right\}  \right] + const.
\nonumber
\end{eqnarray}
Therefore, for a fixed $\mathbf{x}$, if $\mathbf{y}$ and ($\mathbf{z} - \mathbf{x}$) are uncorrelated or independent such that ($\mathbf{z} - \mathbf{x}$) has zero mean vector,
then (\ref{eq:n2nderiv}) is equivalent to the following Noise2Noise loss function in terms of $\boldsymbol\theta$:
\begin{equation}
\mathbb{E}_{(\mathbf{x}, \mathbf{y}, \mathbf{z})} \| \mathbf{z} - \boldsymbol{h_\theta}( \mathbf{y} ) \|^2.
\label{eq:n2n}
\end{equation}
Consequently, the optimal network parameters $\boldsymbol{\theta}$ of a denoiser using (\ref{eq:n2n}) will yield the same solution as the MSE based training with clean ground truth. Noise2Noise achieved outstanding performance in various image restoration tasks including Gaussian noise removal as far as there is the set of two noisy image pairs per one ground truth image~\cite{noise2noise}.

Therefore, this analysis on Noise2Noise can now predict that if there are imperfect ground truth images $\tilde{\mathbf{x}}$ with a mild noise, then denoiser training with
$\tilde{\mathbf{x}}$ plus additional synthetic noise may not be able to yield good performance comparable to the case using two independent noise realizations or using
perfect ground truth data with additional synthetic noise possibly due to non-negligible non-zero term of
\(
\mathbb{E}_{(\mathbf{x} , \mathbf{y}, \mathbf{z})} \left\{ (\mathbf{z} - \mathbf{x})^T \boldsymbol{h_\theta}( \mathbf{y} ) \right\}
\)
in (\ref{eq:n2nderiv}).

\subsection{Extended SURE and MC-SURE}

The original SURE in (\ref{eq:sure}) works well with a single noise realization per image, but it can not take advantage of having multiple noise realizations per image just like Noise2Noise.
Thus, we propose to extend the original SURE to be able to handle pairs of noisy images per ground truth image. The extended SURE (eSURE) can be formulated in the following way:
\begin{thm}
	Let $\mathbf{y_1} \sim \mathcal{N}(\mathbf{x}, \sigma_{\mathbf{y_1}}^2\mathbf{I})$, $\mathbf{z} \sim \mathcal{N}(\mathbf{0}, \sigma_\mathbf{z}^2\mathbf{I})$, and
	$\mathbf{y_2} \triangleq (\mathbf{y_1} + \mathbf{z}) \sim \mathcal{N}(\mathbf{x}, (\sigma_{\mathbf{y_1}}^2 + \sigma_{\mathbf{z}}^2)\mathbf{I})$. Then, the random variable $\gamma (\boldsymbol{h_{\theta}}( \mathbf{y_2}), \mathbf{y_1})$ is an unbiased estimator of MSE:
	\[
	\mathbb{E}_{\mathbf{y_2}}
	\left\{ \frac{1}{N}\| \mathbf{x} - \boldsymbol{h_{\theta}}( \mathbf{y_2} ) \|^2  \right\}
	= \mathbb{E}_{\mathbf{y_2} }
	\left\{  \gamma (\boldsymbol{h_{\theta}}( \mathbf{y_2} ), \mathbf{y_1}) \right\}
	\label{eq:expectnewsure}
	\]
	where $\mathbf{y_1}$ and $\mathbf{z}$ are independent (or uncorrelated) and
	\begin{eqnarray}
	\gamma (\boldsymbol{h_{\theta}}( \mathbf{y_2} ), \mathbf{y_1}) =  \frac{1}{N}\| \mathbf{y_1}
	- \boldsymbol{h_{\theta}}( \mathbf{y_2} ) ) \|^2  - \sigma_{\mathbf{y_1}}^2
	+ \frac{2 \sigma_{\mathbf{y_1}}^2}{N} \sum_{i=1}^N
	\frac{\partial \mathbf{h}_i ( \mathbf{y_2} )}{\partial (\mathbf{y_2})_i }.
	\label{eq:newsure}
	\end{eqnarray}
	\label{thm:newsure}
\end{thm}

Theorem~\ref{thm:newsure} is developed for the general case where imperfect ground truth images with mild Gaussian noise are available and one needs to train the DNN for denoising images contaminated with larger noise level. Moreover, one can train DNN denoisers using the following Corollary of Theorem~\ref{thm:newsure}:
\begin{cor}
	Given a noisy realization pairs of a clean image $(\mathbf{y_3}, \mathbf{y_4})$ from the same distribution $\mathcal{N}(\mathbf{x}, \sigma_\mathbf{y}^2\mathbf{I})$, we calculate less noisy image $\mathbf{w} = \frac{1}{2}(\mathbf{y_3} + \mathbf{y_4}) \sim \mathcal{N}(\mathbf{x}, \frac{1}{2}\sigma_\mathbf{y}^2\mathbf{I})$. Then, we add $i.i.d.$ Gaussian noise
	$\mathbf{z} \sim \mathcal{N}(\mathbf{0}, \frac{1}{2}\sigma_\mathbf{y}^2\mathbf{I})$ to $\mathbf{w}$, so that $\mathbf{v}= (\mathbf{w} + \mathbf{z}) \sim \mathcal{N}(\mathbf{x}, \sigma_\mathbf{y}^2\mathbf{I})$. Finally, by applying Theorem~\ref{thm:newsure} and replacing divergence term with its Monte-Carlo approximation (\ref{eq:approx}), one can minimize extended MC-SURE
	with respect to $\boldsymbol{\theta}$:
	\begin{eqnarray}
	\gamma (\boldsymbol{h_{\theta}}( \mathbf{v} ), \mathbf{w}) = \frac{1}{N}
	\| \mathbf{w} - \boldsymbol{h_\theta} ( \mathbf{v}) \|^2  \label{eq:newsure_dnn2}
	- \frac{1}{2}\sigma_{\mathbf{y}}^2 + \frac{\sigma_{\mathbf{y}}^2}{\epsilon N} (\mathbf{\tilde{n}})^\mathrm{t}\left(  \boldsymbol{h_\theta} ( \mathbf{v} + \epsilon \mathbf{\tilde{n}}) - \boldsymbol{h_\theta} ( \mathbf{v}) \right).
	\end{eqnarray}

\end{cor}
For a training dataset of noisy $M$ pairs $\{ (\boldsymbol{y_3}^{(1)}, \boldsymbol{y_4}^{(1)}), \cdots, (\boldsymbol{y_3}^{(M)}, \boldsymbol{y_4}^{(M)}) \}$, we can generate  $\{ (\boldsymbol{w}^{(j)}, \boldsymbol{v}^{(j)}) \}, j \in [1, M]$ to train deep learning based denoisers with proposed eSURE method. In simulations, our proposed method will show better performance compared to the original MC-SURE based training approach~\cite{soltanayev2018training} for both grayscale and color image denoising. The proof of Theorem~\ref{thm:newsure} and other details can be found in the supplementary material.

\subsection{Link between eSURE and Noise2Noise}
The eSURE framework that we proposed can be applied to a pair of uncorrelated Gaussian noisy images
 ($\mathbf{y} \sim \mathcal{N}(\mathbf{x}, \sigma_{\mathbf{y}}^2)$ and $\mathbf{z} \sim \mathcal{N}(\mathbf{x}, \sigma_{\mathbf{z}}^2)$. In that case, the divergence term vanishes leaving us the following expression:
\begin{eqnarray}
\mathbb{E}_{\mathbf{z, y} }
\left\{  \gamma (\boldsymbol{h_{\theta}}( \mathbf{y} ), \mathbf{z}) \right\}
= \mathbb{E}_{\mathbf{z, y} }
\left\{ \frac{\| \mathbf{z} - \boldsymbol{h_{\theta}}( \mathbf{y} ) ) \|^2}{N} \right\} - \sigma_{\mathbf{z}}^2
\label{eq:link}
\end{eqnarray}
From the above expression, one clearly sees that the first term corresponds to the cost function of Noise2Noise for $i.i.d.$ Gaussian noise denoising case \cite{noise2noise}, while the second term $\sigma^2_{\mathbf{z}}$ is a constant.

Minimization of (\ref{eq:link}) with respect to a set of denoiser parameters $\boldsymbol{\theta}$ should give us the same solution for both Noise2Noise and our eSURE. Although it is not easy to notice the relationship between two different approaches from the first sight, it turns out that Noise2Noise is a special case of proposed extended MC-SURE based training method for $i.i.d.$ Gaussian denoising. A complete derivation can be found in the supplementary material.

\section{Results}

\subsection{Experimental setup}
\label{results}
We have conducted two experiments to evaluate our proposed methods. In the first experiment, we experimentally show that eSURE
efficiently utilizes given two uncorrelated realizations per image to outperform SURE and is
is a general case of Noise2Noise for $i.i.d.$ Gaussian noise.
Second experiment aimed to investigate the effect of noise correlation for Noise2Noise and eSURE with imperfect ground truth.
Our proposed method was compared with
BM3D \cite{Dabov2007fh}, DnCNN trained on MC-SURE \cite{soltanayev2018training}, and DnCNN trained with MSE using noiseless ground truth data.

We used DnCNN \cite{Vincent:2010vu, Zhang:2017eh} as a deep denoising network for grayscale and RGB color images. DnCNN consists of 20 layers of CNN with batch normalization followed by ReLU as a non-linear function. For benchmark test images, we have chosen Berkeley's BSD-68\cite{MartinFTM01}
datasets, and widely used standard test images so called Set 12~\cite{Dabov2007fh}. All experiments were implemented on Tensorflow framework~\cite{Abadi2016} and run on NVidia Titan X GPU.

It is worth to note that $\epsilon$ in (\ref{eq:sure}) and (\ref{eq:newsure}) should be carefully chosen for stable training and high performance. As mentioned in  \cite{soltanayev2018training} and \cite{deledalle2014stein}, $\epsilon$ should be directly proportional to the noise standard deviation $\sigma$. Therefore, $\epsilon$ was fine tuned, so that for our proposed eSURE it was set to be $\epsilon = 1.6 \times 10^{-4} \times \sigma$.

\subsection{The case with two uncorrelated noise realizations per one image} \label{link}
Given an uncorrelated two noisy realizations for each image in the training set, we trained DnCNNs with MC-SURE, Noise2Noise, and our proposed eSURE, respectively. More precisely, by following procedures described in DnCNN paper~\cite{Zhang:2017eh}, we generated 128$\times$2,919 patches with 50$\times$50 size from BSD-400~\cite{MartinFTM01} and produced two independent noisy patches (noise level range $\sigma \in [0-55]$) per clean patch. Using the aforementioned Corollary, we trained DnCNN with eSURE. It is worth to note that MC-SURE requires only a single noisy image per one ground truth in the training set compared to eSURE and Noise2Noise. Therefore, for a fair comparison, we concatenated both noisy datasets to train DnCNN-SURE with twice more data and denoted it as DnCNN-SURE*.
DnCNN denoisers were trained for blind denoising with the noise level range of $\sigma \in [0-55]$ using Adam optimizer \cite{kinga2015method}. An initial learning rate was set to 10\textsuperscript{-3}, which was dropped to 10\textsuperscript{-4} after 40 epochs and the network was further trained for 10 more epochs.

The performance of our approach along with the state-of- the-art methods was tabulated in Table~\ref{table:blind_gray_part1}. Our network training approach demonstrates almost identical quantitative results with Noise2Noise trained DnCNN (DnCNN-N2N) in both test sets.
These results are consistent with our theoretical understandings such as 1) eSURE efficiently utilized two uncorrelated realizations compared to SURE*,
2) (\ref{eq:link}) holds for uncorrelated noisy training set and Noise2Noise is a special case of the extended MC-SURE. Moreover, quantitative analysis on BSD68 test set reveals that our eSURE is consistently better than conventional BM3D for about 0.5dB and outperforms DnCNN-SURE for about 0.15 dB in lower and higher noise cases. The performance gap between the proposed method, BM3D and DnCNN-SURE are still similar for Set12 test set. In addition, we can observe that minimizing MC-SURE with twice more dataset (DnCNN-SURE*) provides a little improvement, but not enough to reach DnCNN-eSURE and DnCNN-N2N.
\begin{table}[!h]
	\caption{PSNR results of blind denoisers on BSD68 and Set12 datasets.}
	\label{table:blind_gray_part1}
	\centering
	\resizebox{\columnwidth}{!}{
		\begin{tabular}{c|ccccc|c}
			\toprule
			\multicolumn{7}{c}{BSD-68}\\
			\midrule
			Methods 		& BM3D & DnCNN-SURE & DnCNN-SURE* & DnCNN-N2N & DnCNN-eSURE & DnCNN-MSE\\
			\midrule
			$\sigma = 25$ &28.56 & 28.92  & 29.00  & \bf{29.08} & \bf{29.08} & 29.20\\
			$\sigma = 50$ &25.62  & 26.00  & 26.07 & 26.13       &\bf{26.15}  & 26.22\\
			\midrule
			\multicolumn{7}{c}{Set 12}\\
			\midrule
			$\sigma = 25$ &29.97  & 30.04  &30.13 & 30.30 & \bf{30.31} & 30.42 \\
			$\sigma = 50$ &26.67  & 26.87  &26.97 & \bf{27.07}   &\bf{27.07}   & 27.16\\
			\midrule
			Noisy dataset	& - & 1  & 2  & 2 & 2 & $\infty$ \\
			\bottomrule
	\end{tabular}}
\end{table}

In terms of visual comparison, our proposed eSURE method effectively removed noise from an image, while preserving texture and edges. In Figure \ref{fig-dncnn_sigma50}, conventional BM3D yielded blurry results. A similar trend is observed for DnCNN-SURE where details of the denoised test image from BSD68 were not fully recovered (see Figure \ref{fig-dncnn_sigma50}). Also one may be able to observe visually similar denoising performance for DnCNN-N2N and our method.

\begin{figure}[ht]
	\centering
	\subfigure[Ground Truth]    		       		 {\includegraphics[width=27mm]{./results/gray_true/clean_sigma50.png}}		\subfigure[BM3D / 25.40 dB] 			 		{\includegraphics[width=27mm]{./results/gray_true/bm3d_sigma50.png}}
	\subfigure[SURE / 25.83 dB] 			  		{\includegraphics[width=27mm]{./results/gray_true/sure_sigma50.png}}
	\subfigure[N2N / \bf{26.01} dB]    				{\includegraphics[width=27mm]{./results/gray_true/n2n_sigma50.png}}
	\subfigure[Ours / \bf{26.01} dB]         	    {\includegraphics[width=27mm]{./results/gray_true/newsure_sigma50.png}}
	\caption{Denoised test (BSD68) results of BM3D, DnCNN trained with various methods for $\sigma=50$.}
	\label{fig-dncnn_sigma50}
\end{figure}

To sum up, it was experimentally shown that for two uncorrelated noise realizations of noisy training set, DnCNN-eSURE and DnCNN-N2N yield almost identical performance. Also, since  eSURE uses less noisy data for training (see Corollary), it better approximates MSE and accordingly outperforms MC-SURE yielding results
that are closer to MSE trained DnCNN.

\subsection{The case with two correlated noise realizations per one image: imperfect ground truth}

In many practical applications, collecting noiseless ground truth images is challenging, expensive, or even infeasible because of camera/object motion and long exposure time. Thus, we may have a limited number of photo shots of a particular scene and by averaging them, we have a ground truth where small amount of noise still remained. Adding synthetic noise to imperfect ground truth to produce noisy images for DNN training produces a dataset where the noise in noisy image may be correlated with the noise in the ground truth. In order to investigate how correlated noise affects our eSURE, Noise2Noise, and original SURE, we have conducted 2 experiments: DnCNN denoisers were trained for a fixed noise level for denoising (e.g. $\sigma_{noisy}=25, 50$) on grayscale BSD-400 and for a blind noise ($\sigma_{noisy} \in [\sigma_{gt}-55]$) on color BSD-432~\cite{MartinFTM01}.

We simulated the case with imperfect ground truth images (or slightly noisy ground truth), by adding synthetic Gaussian noise with $\sigma_{gt}$ to the noiseless clean oracle images. Consequently, noisy training images  were generated by adding $i.i.d.$ Gaussian noise on the top of the imperfect ground truth dataset. Following the same procedures in Section~\ref{link}, 128$\times$2,919 patches with 50$\times$50 size for grayscale and 128$\times$2,019 patches for RGB case were generated. DnCNN denoisers were trained using Adam optimizer \cite{kinga2015method}. The initial learning rate was set to 10\textsuperscript{-3}, which was dropped to 10\textsuperscript{-4} after 40 epochs and the network further trained for 10 more epochs.

\begin{table}[!h]
	\caption{Results of denoising methods on BSD68 and Set 12 datasets (Performance in dB).}
	\label{table:fixed_gray_part2}
	\centering

	\begin{tabular}{l|ccc|cccc}
		\toprule
		\multicolumn{8}{c}{BSD-68} \\
		\midrule
		$\sigma_{noisy}$  & \multicolumn{3}{c}{25} \vline & \multicolumn{4}{c}{50}   \\
		\midrule
		$\sigma_{gt}$ 		 &  1 & 5 & 10 &  1 & 5 & 10 & 20									\\
		\midrule
		BM3D 					  & \multicolumn{3}{c}{28.56} \vline & \multicolumn{4}{c}{25.62} \\
		DnCNN-SURE		  & 29.05 &	29.01 &	29.02 &	25.95 &	25.97 &	25.90 &	25.92 \\
		DnCNN-N2N		   & 29.23 & 29.15 & 28.37 & \bf{26.28} & 26.24 & 25.91	&24.69 \\
		DnCNN-eSURE 	&\bf{29.23}	&\bf{29.23}	&\bf{29.21}	&26.27	&\bf{26.24}	&\bf{26.27}	&\bf{26.25} \\
		\midrule
		DnCNN-MSE		& \multicolumn{3}{c}{29.23} \vline & \multicolumn{4}{c}{26.28} \\
		\midrule
		\multicolumn{8}{c}{Set 12} \\
		\midrule
		BM3D 					 & \multicolumn{3}{c}{29.97} \vline & \multicolumn{4}{c}{26.67} \\
		DnCNN-SURE		 & 30.23 &	30.19 &	30.19 &	26.77 &	26.85 &	26.73 &	26.74 \\
		DnCNN-N2N		 & 30.41 &	30.39 &	29.46 &	\bf{27.28} &	27.20 &	27.08 &	25.39 \\
		DnCNN-eSURE    &\bf{30.47} &	\bf{30.48} &	\bf{30.44} & 27.27&	\bf{27.27} &	\bf{27.25}&	\bf{27.23} \\
		\midrule
		DnCNN-MSE 		& \multicolumn{3}{c}{30.47} \vline & \multicolumn{4}{c}{27.28} \\
		\bottomrule
	\end{tabular
\end{table}

\begin{table}[!h]
	\caption{Results of denoising methods on RGB color BSD68 dataset(Performance in dB).}
	\label{table:blind_rgb_part2}
	\centering

	\begin{tabular}{l|ccc|cccc}
		\toprule
		\multicolumn{8}{c}{RGB BSD-68} \\
		\midrule
		$\sigma_{noisy}$  & \multicolumn{3}{c}{25} \vline & \multicolumn{4}{c}{50}   \\
		\midrule
		$\sigma_{gt}$ 		 &  1 & 5 & 10 & 1 & 5 & 10 & 20									\\
		\midrule
		CBM3D 					  & \multicolumn{3}{c} {30.70} \vline & \multicolumn{4}{c}{27.38} \\
		CDnCNN-SURE		  &30.97	    &30.98			  &30.99   			&27.63	   			 &27.68			   &27.64		 	 &27.63 \\
		CDnCNN-N2N		   &31.18  		  &31.08     		&29.83     		  &27.89			   &27.87			 &27.61   		    &25.62 \\
		CDnCNN-eSURE 	 &\bf{31.20} &\bf{31.18}	 &\bf{31.19}   	  &\bf{27.94}		&\bf{27.91}	    &\bf{27.90}     &\bf{27.78} \\
		\midrule
		DnCNN-MSE			& \multicolumn{3}{c}{31.20} \vline & \multicolumn{4}{c}{27.93} \\
		\bottomrule
	\end{tabular
\end{table}

\begin{figure}[!t]
	\centering
	\subfigure[Ground Truth]    		       		 {\includegraphics[width=46mm]{./results/rgb_gt5_25/gt.png}}
	\subfigure[Noisy / 20.68 dB]        			 {\includegraphics[width=46mm]{./results/rgb_gt5_25/noisy.png}}
	\subfigure[BM3D / 30.82 dB] 			 	   {\includegraphics[width=46mm]{./results/rgb_gt5_25/cbm3d.png}}
	\subfigure[SURE / 30.90 dB] 			  		{\includegraphics[width=46mm]{./results/rgb_gt5_25/sure.png}}
	\subfigure[N2N /29.89 dB]    				     {\includegraphics[width=46mm]{./results/rgb_gt5_25/n2n.png}}
	\subfigure[Ours / \bf{31.17} dB]         	    {\includegraphics[width=46mm]{./results/rgb_gt5_25/esure.png}}
	\caption{CBM3D and CDnCNN results on test image from BSD-68 with noise $\sigma=25$. CDnCNN trained with imperfect ground truth with $\sigma_{gt} = 10$ for blind noise denoising task using various approaches. }
	\label{fig-dncnn_rgb_sigma50}

\end{figure}

Table~\ref{table:fixed_gray_part2} shows the performance of denoising methods trained for a fixed noise given a noisy ground-truth images with $\sigma_{gt} = \{1, 5, 10, 20\}$. The higher the $\sigma_{gt}$ is, the more correlated noise is in a training set. We notice that at low level of ground-truth noise ($\sigma_{gt}=1$), both Noise2Noise (DnCNN-N2N) and eSURE (DnCNN-eSURE) yield the best PSNR results and even comparable to the DnCNN-MSE. However, as noise correlation gets severe ($\sigma_{gt}=5, 10$), Noise2Noise fails to achieve high performance that is consistent with our theoretical derivation.
In contrast, the proposed DnCNN eSURE produces the best quantitative results in a stable manner. Although, DnCNN-SURE is not susceptible to noise correlation, it still yielded worse performance than
DnCNN-eSURE.

The experimental results for RGB color image denoising case are tabulated in Table~\ref{table:blind_rgb_part2}. In this case, we observe the same performance degradation pattern of DnCNN-N2N as the noise in ground truth image increases (more results in the supplementary material). The visual assessment of methods also demonstrates that eSURE trained DnCNN was able to provide high-quality images with preserved texture and color. Moreover, we can see from Figure~\ref{fig-dncnn_rgb_sigma50} that denoised output image of CBM3D is highly smoothed out and DnCNN-N2N recovered image still have some noise, while our eSURE denoised image shows sharp edges with almost no noise.

\section{Conclusion} \label{conclusion}
We have investigated properties of Noise2Noise and proposed eSURE that extended the original SURE to handle correlated pairs of noisy images efficiently
for training DNN denoisers. For two uncorrelated noisy realizations per image, eSURE yielded better performance than SURE that implies efficient utilization of
two uncorrelated noisy realizations as compared to SURE and SURE*. Our eSURE also yielded comparable performance to Noise2Noise that is consistent with our theoretical analysis.
For two correlated noisy realizations per image or imperfect ground truth, eSURE still yielded the best performance among all compared methods such as
BM3D, SURE, and Noise2Noise. However, Noise2Noise did not yield good performance with correlated noisy realizations as predicted based on our theoretical analysis.

	\section*{Acknowledgments}

	This work was supported partly by
	Basic Science Research Program through the National Research Foundation of Korea(NRF)
	funded by the Ministry of Education(NRF-2017R1D1A1B05035810),
	the Technology Innovation Program or Industrial Strategic Technology Development Program
	(10077533, Development of robotic manipulation algorithm for grasping/assembling
	with the machine learning using visual and tactile sensing information)
	funded by the Ministry of Trade, Industry \& Energy (MOTIE, Korea), and a grant of the Korea Health Technology R\&D Project
	through the Korea Health Industry Development Institute (KHIDI),
	funded by the Ministry of Health \& Welfare, Republic of Korea (grant number: HI18C0316).

\small
\bibliographystyle{unsrt}


\section{}

\end{document}