# Imaging with thermal noise induced currents††thanks: Submitted to the
editors May 8 2023. This work was partially funded by the National Science
Foundation grants DMS-2008610 and DMS-2136198.

Trent DeGiovanni Mathematics Department, University of Utah, Salt Lake City,
UT 84112<EMAIL_ADDRESS><EMAIL_ADDRESS>Fernando Guevara
Vasquez22footnotemark: 2 China Mauck Formerly: Mathematics Department,
University of Utah, Salt Lake City, UT 84112. Currently: STV Incorporated, 200
W Monroe St #1650, Chicago, IL 60606.

###### Abstract

We use thermal noise induced currents to image the real and imaginary parts of
the conductivity of a body. Covariances of the thermal noise currents measured
at a few electrodes are shown to be related to a deterministic problem. We use
the covariances obtained while selectively heating the body to recover the
real power density in the body under known boundary conditions and at a known
frequency. The resulting inverse problem is related to acousto-electric
tomography, but where the conductivity is complex and only the real power is
measured. We study the local solvability of this problem by determining where
its linearization is elliptic. Numerical experiments illustrating this inverse
problem are included.

###### keywords:

Conductivity imaging, Hybrid inverse problems, Thermal noise, Johnson-Nyquist
noise

Imaging with thermal noise induced currentsT. DeGiovanni, F. Guevara Vasquez,
and C. Mauck

35R30, 35J25, 35Q61

##  1 Introduction

In an electrical conductor, the excitement of charge carriers due to heat
produces random currents. This phenomenon is called Johnson-Nyquist noise and
was first observed in the early 20th century [21, 34]. Given a single
component with impedance $Z(\omega)$ (in Ohms) at an angular frequency
$\omega$ (in $2\pi$ Hz), the variance
$\left\langle{|J(\omega)|^{2}}\right\rangle$ of the random currents (in
$\textrm{A}^{2}$) is

(1) $\left\langle{|J(\omega)|^{2}}\right\rangle=\frac{2\kappa
T}{\pi}\frac{\text{Re}\,(Z(\omega))}{|Z(\omega)|^{2}}\Delta\omega.$

Here $\kappa\approx 1.36\times 10^{-23}\textrm{J}\cdot\textrm{K}^{-1}$ is
Boltzmann’s constant, $T$ is temperature (in Kelvin) and $\Delta\omega$ is a
bandwidth of interest around $\omega$. Although this noise in a nuisance in
electrical circuits, we show one way to use it to image the conductive
properties of a body. Johnson-Nyquist noise and its generalization to the
Maxwell equations (see, e.g., [37]) are examples of a more general physical
principle called the fluctuation dissipation theorem, which relates the
variance of fluctuations of a linear system about an equilibrium to the
dissipative properties of the system, see e.g. [23, 36, 38].

To image the conductivity of a body, we propose heating the body while
simultaneously measuring the variance of the thermal noise induced currents
using electrodes that are connected to the ground. For instance, this can be
done on a two-dimensional conductive body as illustrated in Fig. 1, with
electrodes on its boundary that are connected to the ground (zero voltage or
potential). The electrical measurements are made while the body is heated at a
known spatial location via an external source, e.g. a laser, and the process
is repeated at different locations to scan the body. In our approach, we also
need to subtract measurements of thermal noise induced currents at a known and
constant background temperature.

Figure 1: A two-dimensional conductive body is attached to the ground via
electrodes which are used to measure the thermal noise currents resulting from
heating the body at particular locations, e.g. the location depicted in red.

Our main contribution is to show that such thermal noise current measurements
are equivalent to measuring the real power dissipated inside the conductive
body, i.e.

(2) $\sigma^{\prime}(x)|\nabla u(x)|^{2},$

where $\sigma^{\prime}(x)$ is the real part of the conductivity and $u$ solves
an appropriate (deterministic) auxiliary problem which depends on electrodes
configuration and the conductivity $\sigma(x)$ which can be complex (see
Section 3).

###  1.1 Related work

Recovering $\sigma^{\prime}(x)$ from functionals of the form (2) is well-
studied for the case of real $\sigma$ in ultrasound modulated electrical
impedance tomography or acousto-electric tomography [2], where the internal
functional (2) is measured by locally perturbing a conductive body using
ultrasound waves, while making electrical measurements on the body’s surface.
Various reconstruction approaches have since been studied for this problem
[14, 24, 25, 9, 7] as well as its well-posedness [20, 26, 4]. The problem of
recovering an anisotropic real conductivity has also been studied [8, 30, 31,
32]. A similar problem using microwaves instead of ultrasound is discussed in
[3]. Optical tomography can also be modulated by ultrasound, allowing
measurements of a functional similar to (2), see [10, 35]. Hybrid inverse
problems (including acousto-electric tomography) have been studied by
formulating them as an overdetermined system of non-linear partial
differential equations and then studying their local uniqueness properties by
linearizing, see [26, 6]. For reviews on hybrid inverse problems, see [1, 5].

The $\sigma$ complex case is considered in [11], but this analysis only
applies if the fields are known (as in elastography). Complex $\sigma$ were
also considered in the case of the Maxwell equations in [12]. However, to our
knowledge, there is no study of the functional (2) where $u$ depends on a
complex $\sigma$, but only its real part $\sigma^{\prime}$ appears explicitly
in the measurements.

###  1.2 Possible applications

The biggest challenge to the applicability of the method that we present here
is that the thermally induced random currents are very small and this could
introduce signal-to-noise issues. We envision two possible applications.

The first possible application would be to Atomic Force Microscopy (AFM). In
this imaging modality, a height-map of a sample is obtained by measuring the
deflections of a cantilever as its tip scans the sample. A heated cantilever
tip can be used to heat the sample locally without touching it, see e.g. [22].
Moreover, electrical measurements of thermal noise induced currents can be
done simultaneously with the AFM scan. An advantage of this approach is that
one can measure height and conductivity of the sample without touching the
sample, possibly making the cantilever tip last longer. We mention that
conductivity variations can be measured in AFM by creating a voltage
difference between the sample and the cantilever (assuming it is conductive).
This method is known as Conductive Atomic Force Microscopy or CAFM, see e.g.
the review in [28].

The second possible application is to monitoring of laser welding, see e.g.
[17]. If two sheets of metal are been welded together, their temperature is
raised significantly near the weld and the sheets also become electrically
connected. We believe that by measuring thermal noise currents, one can
monitor whether the weld was effective. We give an order of magnitude of the
signals and background that would need to be measured in this situation in
Section 3.3.

###  1.3 Contents

We start in Section 2 by deriving a quasi-static model from the Maxwell
equations with a current source modeling the random currents. In Section 3 we
show how variances of the random currents are related to a deterministic
problem. This is done for two different kinds of boundary conditions. Moreover
we give rough magnitude estimates for the currents that would need to be
measured to implement our approach. In Section 4, we analytically and
numerically analyze the linearized real problem (Section 4.1) and linearized
complex problem (Section 4.2). This analysis is based on [6] where the
ellipticity, in the Douglis-Nirenberg sense [18], is established for the real
linearized problem. We give a condition in Lemma 4.1 for the linearized
problem with complex conductivity to be elliptic in the Douglis-Nirenberg
sense [18]. Still, it remains unclear if boundary conditions exist such that
the fields associated with the auxiliary problem satisfy this condition. Then
in Section 5 we present a simple numerical reconstruction approach based on a
finite difference discretization of the problem (Section 5.1). We solve the
inverse problem using data that either comes directly from the internal
functional (2) or from simulated realizations of random currents. In addition,
we show reconstructions in the case that conductivity is real (Section 5.2) or
complex (Section 5.3). Finally, we summarize our results in Section 6.

##  2 The quasi-static model

In an isotropic medium, thermal fluctuations induce fluctuation of charge
carriers near an equilibrium. For non-magnetic media, the thermal fluctuation
currents can be modeled by a random external electric current $j_{e}$
($\textrm{A}/\textrm{m}^{2})$ in the Maxwell equations [37], namely

(3) $\displaystyle\nabla\times H$ $\displaystyle=i\omega\varepsilon E+j_{e}$
$\displaystyle\nabla\times E$ $\displaystyle=-i\omega\mu H.$

Here $E$ and $H$ are the electric and magnetic fields, and the angular
frequency is $\omega$. The convention for time harmonic fields here is that
$\mathcal{E}(x,t)=\text{Re}\,\left[E(x,\omega)\exp[i\omega t]\right]$. The
electric permittivity is $\varepsilon$ and may be written as
$\varepsilon=\varepsilon^{\prime}-i\sigma^{\prime}/\omega$, where
$\varepsilon^{\prime}\equiv\text{Re}\,\varepsilon$ and $\sigma^{\prime}$ is
the real conductivity. The magnetic permeability $\mu$ is assumed real and
equal to that of the vacuum. Note that if $\mu$ had an imaginary part (i.e.
non-zero magnetic losses), then an analogous “random magnetic current” needs
to be added to the Maxwell equations. The fluctuation dissipation theorem (see
e.g. [38, Chapter 1] and the particular application to the Maxwell equations
in [37]) states that the random current field $j_{e}$ has zero mean
$\left\langle{j_{e}}\right\rangle=0$ and its covariance (at a fixed frequency)
is

(4)
$\left\langle{j_{e}(x)j_{e}^{*}(x^{\prime})}\right\rangle=\frac{\kappa}{\pi}T(x)\sigma^{\prime}(x)\delta(x-x^{\prime})\mathbb{I}\Delta\omega,$

where $T(x)$ is the temperature in Kelvin at a point $x$, $\mathbb{I}$ is the
identity matrix and $\Delta\omega$ is a frequency band of interest around
$\omega$. We emphasize that (4) depends only on the temperature and the real
part of the conductivity. The real part of the electrical permittivity (which
is associated with lossless behavior) does not directly appear in (4). While
(4) gives the entire covariance matrix, only its diagonal entries are needed
for our reconstructions. In general (4) should use the energy of a quantum
oscillator [37] instead of $\kappa T$, namely

(5) $\Theta(T,\omega)=\frac{\hslash\omega}{2}\coth\frac{\hslash\omega}{2\kappa
T},$

where $\hslash\approx 1.05\times 10^{-34}\textrm{J}\cdot\textrm{s}$ is
Planck’s constant. Here we assume we work with relatively small frequencies so
that $\kappa T\gg\hslash\omega$ and we can make the approximation
$\Theta(T,\omega)\approx\kappa T$, see also [27]. In particular, this
approximation is valid at room temperature and frequencies of the order of
$1$kHz or $1$MHz.

Instead of working with the Maxwell equations, we use a quasi-static
approximation that is used in electrical impedance tomography, see e.g. [15,
13]. In this approximation, it is convenient to define the complex
conductivity $\sigma$ by

(6) $\sigma(x)=\sigma^{\prime}(x)+i\omega\varepsilon^{\prime}.$

If we assume that $\omega\mu|\sigma|L^{2}\ll 1$, where $L$ is the
characteristic length of the problem, then one can use the approximation
$\nabla\times E\approx 0$. In other words, we may assume that the electric
field comes from a potential $E=-\nabla\phi$. By taking divergence on both
sides of the first equation in (3) we get

(7) $\nabla\cdot[\sigma\nabla\phi]=\nabla\cdot j_{e}.$

######  Remark 1.

As noted in [15], the quasi-static approximation holds for conductivities
consistent with human tissues (see, e.g., [13]). For example if we take
$L=10~{}\mathrm{cm}$,
$\sigma^{\prime}=2~{}\mathrm{cm}^{-1}\mathrm{k}\Omega^{-1}$, $\omega=2\pi
10~{}\mathrm{kHz}$ and $\epsilon^{\prime}=1\mu\mathrm{F}/\mathrm{m}$, we get
$\omega\mu|\sigma|L^{2}\approx 1.7\times 10^{-4}\ll 1$.

##  3 From the stochastic to the deterministic problem

Let $\Omega$ be a smooth simply connected open domain of $\mathbb{R}^{3}$ and
let $\sigma\in C^{1}(\overline{\Omega})$ with its real part satisfying
$\sigma^{\prime}>c$ for some positive constant $c$. We assume a potential
$\phi$ satisfies

(8) $\displaystyle\nabla\cdot[\sigma\nabla\phi]$ $\displaystyle=\nabla\cdot
j_{e},~{}\text{in}~{}\Omega,$ $\displaystyle\phi$
$\displaystyle=0,~{}\text{on}~{}\partial\Omega.$

Here $j_{e}$ is the random current term with
$\left\langle{j_{e}}\right\rangle=0$ and
$\left\langle{|j_{e}|^{2}}\right\rangle$ given in (4). We assume that $j_{e}$
is $C^{1}(\overline{\Omega})$ and hence $\phi$ is $C^{2}(\overline{\Omega})$
[19, Ch. 6.3].

We assume we measure currents flowing out of the domain $\Omega$ at $n$
“electrodes” by the complex vector with $n$ entries

(9) $J=\int_{\partial\Omega}dS(x)\begin{bmatrix}e_{1}(x)\\\ \vdots\\\
e_{n}(x)\end{bmatrix}\sigma(x)\nabla\phi(x)\cdot\nu(x),$

where $\nu(x)$ is the unit outward pointing normal to $\partial\Omega$ at some
$x\in\partial\Omega$, and the function $e_{i}(x)$ are possibly complex
$C^{1}(\partial\Omega)$ “electrode functions” defined on $\partial\Omega$. For
example they could be a continuously differentiable approximation of the
characteristic function of electrodes at the boundary.

In the following result we prove that the $n\times n$ covariance matrix
$\left\langle{JJ^{*}}\right\rangle$ of such measurements can be related to
solutions to deterministic auxiliary problems.

######  Theorem 3.1.

The covariance of the vector of measurements $J$ are given by

(10)
$[\left\langle{JJ^{*}}\right\rangle]_{ij}=\frac{\kappa}{\pi}\int_{\Omega}dy\,\text{Re}\,(\sigma(y))T(y)\nabla
u_{i}(y)\cdot\nabla\overline{u_{j}(y)}\Delta\omega,$

where the functions $u_{j}$ are solutions to the Dirichlet problems

(11) $\displaystyle\nabla\cdot[\sigma\nabla u_{i}]$
$\displaystyle=0,~{}\text{in}~{}\Omega$ $\displaystyle u_{i}$
$\displaystyle=e_{i},~{}\text{on}~{}\partial\Omega.$

######  Proof 3.2.

First, note that we can write the solution to (8) as

(12) $\phi(x)=\int_{\Omega}dy\,G(x,y)\nabla_{y}\cdot j_{e}(y),$

where $G(x,y)$ is the Green function $G(x,y)$ satisfying the equation

(13) $\displaystyle\nabla_{x}\cdot[\sigma(x)\nabla_{x}G(x,y)]$
$\displaystyle=\delta(x-y),~{}x,y\in\Omega$ $\displaystyle G(x,y)$
$\displaystyle=0,~{}x\in\partial\Omega~{}\text{or}~{}y\in\partial\Omega.$

To double check (12), it is clear that $\phi(x)=0$ for $x\in\partial\Omega$
because $G(x,y)=0$ for $x\in\partial\Omega$. Also:

$\displaystyle\nabla_{x}\cdot[\sigma(x)\nabla_{x}\phi]$
$\displaystyle=\int_{\Omega}dy\,\nabla_{x}\cdot[\sigma(x)\nabla_{x}G(x,y)]\nabla_{y}\cdot
j_{e}(y)$ $\displaystyle=\int_{\Omega}dy\,\delta(x-y)\nabla_{y}\cdot j_{e}(y)$
$\displaystyle=\nabla_{x}\cdot j_{e}(x).$

Now it is helpful to use integration by parts to get

(14) $\displaystyle\phi(x)$
$\displaystyle=\int_{\Omega}dy\,G(x,y)\nabla_{y}\cdot j_{e}(y)$
$\displaystyle=\int_{\partial\Omega}dS(y)\,G(x,y)j_{e}(y)\cdot\nu(y)-\int_{\Omega}dy\,\nabla_{y}G(x,y)\cdot
j_{e}(y)$ $\displaystyle=-\int_{\Omega}dy\,\nabla_{y}G(x,y)\cdot j_{e}(y).$

Using (14) and (4) allows us to compute

(15) $\displaystyle\left\langle{JJ^{*}}\right\rangle_{ij}$
$\displaystyle=\left\langle{\int_{\partial\Omega}dS(x)\int_{\partial\Omega}dS(x^{\prime})e_{i}(x)\sigma(x)\nabla_{x}\phi(x)\cdot\nu(x)\overline{e_{j}(x^{\prime})}\overline{\sigma(x^{\prime})}\nabla_{x^{\prime}}\overline{\phi(x^{\prime})}\cdot\nu(x^{\prime})}\right\rangle$
$\displaystyle=\int_{\partial\Omega}dS(x)\int_{\partial\Omega}dS(x^{\prime})\int_{\Omega}dy\frac{\kappa}{\pi}\text{Re}\,(\sigma(y))T(y)\Delta\omega
e_{i}(x)\sigma(x)\overline{e_{j}(x^{\prime})}\overline{\sigma(x^{\prime})}\nu(x)^{T}$
$\displaystyle\hskip
130.88284pt\nabla_{x}\nabla_{y}G(x,y)\nabla_{x^{\prime}}\nabla_{y}\overline{G(x^{\prime},y)}\nu(x^{\prime})$
$\displaystyle=\frac{\kappa}{\pi}\int_{\Omega}dy\,\text{Re}\,(\sigma(y))T(y)\Delta\omega\left[{\int_{\partial\Omega}dS(x)e_{i}(x)\sigma(x)\nabla_{x}\nabla_{y}G(x,y)\nu(x)}\right]^{T}$
$\displaystyle\hskip
128.0374pt\overline{\left[{\int_{\partial\Omega}dS(x^{\prime})e_{j}(x^{\prime})\sigma(x^{\prime})\nabla_{x^{\prime}}\nabla_{y}G(x^{\prime},y)\nu(x^{\prime})}\right]}$
$\displaystyle=\frac{\kappa}{\pi}\int_{\Omega}dy\,\text{Re}\,(\sigma(y))T(y)\Delta\omega\nabla_{y}u_{i}(y)\cdot\nabla_{y}\overline{u_{j}(y)},$

where $u_{i}$ solves the problem (11) and we use that
$\nabla_{x}\nabla_{y}G(x,y)$ is symmetric. The last equality follows by doing
integration by parts twice:

(16)
$\displaystyle\int_{\partial\Omega}dS(x)e_{i}(x)\sigma(x)\nabla_{x}G(x,y)\cdot\nu(x)$
$\displaystyle=\int_{\Omega}dx\nabla_{x}\cdot[\sigma(x)\nabla_{x}G(x,y)u_{i}(x)]$
$\displaystyle=\int_{\Omega}dx\nabla_{x}\cdot[\sigma(x)\nabla_{x}G(x,y)]u_{i}(x)$
$\displaystyle+\int_{\Omega}dx\sigma(x)\nabla_{x}G(x,y)\cdot\nabla_{x}u_{i}(x)$
$\displaystyle=u_{i}(y)+\int_{\Omega}dx\nabla_{x}\cdot[\sigma(x)G(x,y)\nabla_{x}u_{i}(x)]$
$\displaystyle-\int_{\Omega}dxG(x,y)\nabla_{x}\cdot[\sigma(x)\nabla_{x}u_{i}(x)]$
$\displaystyle=u_{i}(y)+\int_{\partial\Omega}dS(x)G(x,y)\sigma(x)\nabla_{x}u_{i}(x)\cdot\nu(x)$
$\displaystyle=u_{i}(y).$

It may be possible to loosen the regularity assumptions on $j_{e}$ and $e_{i}$
and derive a similar result to Theorem 3.1. Since the scope of this work is
focused on establishing the relation between the stochastic and deterministic
problems, we leave this for future work.

###  3.1 Boundary conditions modeling electrodes with insulating gaps

The setup using Dirichlet boundary conditions (8) assumes that
$\phi|_{\partial\Omega}=0$, which would likely be hard to realize in practice
because we expect to have a few electrodes connected to the ground with
insulating gaps between them. This corresponds to a boundary condition of
mixed type: homogeneous Dirichlet on the electrodes and homogeneous Neumann
(zero flux) on the gaps between the electrodes. To be more precise, let
$\Gamma=\operatorname{supp}e_{1}\cup\ldots\cup\operatorname{supp}e_{n}$ then
we replace (8) with

(17) $\displaystyle\nabla\cdot[\sigma\nabla\phi]$ $\displaystyle=\nabla\cdot
j_{e},~{}\text{in}~{}\Omega,$ $\displaystyle\sigma\nabla\phi\cdot\nu$
$\displaystyle=0,~{}\text{on}~{}\partial\Omega-\Gamma,$ $\displaystyle\phi$
$\displaystyle=0,~{}\text{on}~{}\Gamma.$

Then Theorem 3.1 holds in the same fashion, but we assume that the $u_{i}$ are
solutions to the following mixed boundary problem replacing (11) with

(18) $\displaystyle\nabla\cdot[\sigma\nabla u_{i}]$
$\displaystyle=0,~{}\text{in}~{}\Omega,$ $\displaystyle\sigma\nabla
u_{i}\cdot\nu$ $\displaystyle=0,~{}\text{on}~{}\partial\Omega-\Gamma,$
$\displaystyle u_{i}$ $\displaystyle=e_{i},~{}\text{on}~{}\Gamma.$

The proof follows by noting that integration by parts now yields

$\phi(x)=\int_{\partial\Omega-\Gamma}dS(y)G(x,y)j_{e}(y)\cdot\nu(y)-\int_{\Omega}dy\nabla_{y}G(x,y)\cdot
j_{e}(y),$

resulting in four terms when $\phi(x)$ is substituted in
$\left\langle{JJ^{*}}\right\rangle_{ij}$. One of the terms is similar to the
case of the Dirichlet boundary conditions, and all of the others contain
integrals over the zero flux part of the boundary. By invoking the zero flux
boundary conditions, these terms can be easily seen to disappear, leaving us
with the same formula for $\left\langle{JJ^{*}}\right\rangle_{ij}.$

###  3.2 Differential temperature measurements

Utilizing Theorem 3.1, we can now relate the differential temperature
measurements as described in Section 1 to measurements of the internal
functional (2). Concretely, we take a set of measurements of the covariance of
the currents in a body at temperatures $T=T_{0}$ and $T=T_{0}+\delta T(x)$,
where $\delta T(x)$ is a prescribed heating pattern. Then the differential
temperature measurements give

(19) $\left[{\left\langle{J_{T_{0}+\delta T}J_{T_{0}+\delta
T}^{*}}\right\rangle-\left\langle{J_{T_{0}}J_{T_{0}}^{*}}\right\rangle}\right]_{ii}=\frac{\kappa}{\pi}\int_{\Omega}dx\,\delta
T(x)\Delta\omega\text{Re}\,(\sigma(x))|\nabla u_{i}(x)|^{2},$

considering only the diagonal elements of the covariance matrix. As previously
noted, only measurements of the diagonal elements are used for our
reproduction approach. By taking a sufficiently rich set of heating patterns
we get an estimate for the internal functional

(20) $H_{ii}(x)=\sigma^{\prime}(x)|\nabla
u_{i}(x)|^{2},~{}\text{for}~{}x\in\Omega,~{}i,j\in\\{1,\ldots,n\\}.$

For real conductivities ($\sigma=\sigma^{\prime}$), the internal functional
(20) corresponds to the power dissipated inside the domain.

######  Remark 3.3.

In our numerical experiments (see Section 5) we use, for convenience, heating
patterns $\delta T(x)$ that are approximate Dirac delta distributions. Other
patterns such as cosines and sines could also be used. Using spatially
extended patterns may be advantageous in terms of signal to noise ratio.

###  3.3 Rough estimation of thermal noise induced currents

The thermal noise induced currents are very small and may limit the
application of this approach. To get an idea of the magnitude of the signals
that need to be measured to obtain $H_{ii}$ in (19), we need to distinguish
between current measurements with the background temperature $T_{0}$ and with
perturbed temperature $T_{0}+\delta T$. We make rough estimates of these
currents in two situations: the first is consistent with the numerical
experiments and the second one is consistent with laser welding.

#### Conductivities used in the numerical experiments

For the background temperature measurements, recall that Boltzmann’s constant
is on the order of $10^{-23}$ $\textrm{J}\cdot\textrm{K}^{-1}$. For our
numerical experiments we chose $\Delta\omega=10$ kHz, $\Delta z=0.1$ cm and a
domain with area 10 cm2. If the conductivity is about $10^{-3}$ cm-1
$\Omega^{-1}$ and $T_{0}=300$ K, then accounting for the $1/\pi$ factor, the
variance of the random currents is on the order of $10^{-20}$ A2. To reach
this estimate we assumed the squared gradient of the auxiliary fields is
constant and equal to $10^{-2}$ cm-2. For the differential measurements we may
further assume a $\delta T=10$ K on area of $(0.2)^{2}$ cm2. This gives a
current variance of the order $10^{-25}$ A2 and a signal to noise ratio of
$10^{-5}$.

#### Conductivities consistent with welding

The conductivity of gold is much higher than what we used in the numerical
experiments and is on the order of $4.5\times 10^{7}$m-1 $\Omega^{-1}$. For
instance consider a sheet of gold of dimensions 1 cm $\times$ 1 cm $\times$ 1
mm and a bandwidth and central frequencies on the order of 100 Hz. For this
choice of frequencies, the quasi-static approximation (Section 2) is not well
satisfied. Nevertheless, if $T_{0}=300$ K the variance of the random currents
is on the order of $10^{-14}\textrm{A}^{2}$. For the differential measurements
we may further assume a $\Delta T=1300$ K (which is close to the melting point
of gold) on an area of $(0.1)^{2}$ mm2. This gives a current variance of the
order $10^{-17}$ A2 and a signal to noise ratio of $10^{-3}$.

###  3.4 The inverse problem for real conductivities

The inverse problem for a real conductivity $\sigma=\sigma^{\prime}$ consists
of the measurement equation (20) and the auxiliary problem (11). To be more
precise, we seek to recover $u_{i}$ and $\sigma$ given $H_{ii}$ and $e_{i}$
from the real non-linear system of partial differential equations, for
$i=1,...,n$,

(21) $\displaystyle\nabla\cdot\left[{\sigma\nabla u_{i}}\right]$
$\displaystyle=0,\quad x\in\Omega,$ $\displaystyle u_{i}-e_{i}$
$\displaystyle=0,\quad x\in\partial\Omega,$ $\displaystyle
H_{ii}-\sigma|\nabla u_{i}|^{2}$ $\displaystyle=0,\quad x\in\Omega.$

We call the model associated with measurements $H_{ii}$ given by the
expectation in Theorem 3.1 the deterministic model and the model associated
with measurements given by realizations of randomly induced currents the
stochastic model.

###  3.5 The inverse problem for complex conductivities

We write the complex problem by separating the real and complex parts of (20)
and (11). To avoid confusion with the complex number $i$ we use $j$ to denote
experiments for the complex conductivity problem. We use a single prime (resp.
double prime) to denote the real part (resp. imaginary part) of a complex
quantity, e.g. $\sigma=\sigma^{\prime}+i\sigma^{\prime\prime}$,
$u_{j}=u_{j}^{\prime}+iu_{j}^{\prime\prime}$, and
$e_{j}=e_{j}^{\prime}+ie_{j}^{\prime\prime}$. Then the problem is to find
$\sigma^{\prime}$, $\sigma^{\prime\prime}$, $u_{j}^{\prime}$ and
$u_{j}^{\prime\prime}$ given $H_{jj}$, $e_{j}^{\prime}$, and
$e_{j}^{\prime\prime}$ from the non-linear system of partial differential
equations, for $j=1,..,n$,

(22) $\displaystyle\nabla\cdot\left[{\sigma^{\prime}\nabla
u_{j}^{\prime}}\right]-\nabla\cdot\left[{\sigma^{\prime\prime}\nabla
u_{j}^{\prime\prime}}\right]$ $\displaystyle=0,\quad x\in\Omega,$
$\displaystyle\nabla\cdot\left[{\sigma^{\prime}\nabla
u_{j}^{\prime\prime}}\right]+\nabla\cdot\left[{\sigma^{\prime\prime}\nabla
u_{j}^{\prime}}\right]$ $\displaystyle=0,\quad x\in\Omega,$ $\displaystyle
u_{j}^{\prime}-e_{j}^{\prime}$ $\displaystyle=0,\quad x\in\partial\Omega,$
$\displaystyle u_{j}^{\prime\prime}-e_{j}^{\prime\prime}$
$\displaystyle=0,\quad x\in\partial\Omega,$ $\displaystyle
H_{jj}-\sigma^{\prime}(|\nabla u_{j}^{\prime}|^{2}+|\nabla
u_{j}^{\prime\prime}|^{2})$ $\displaystyle=0,\quad x\in\Omega.$

An equivalent formulation of (22) can be found using the conjugates of $u_{j}$
and the $e_{j}$ instead of their real and imaginary components separately.
Both the system (21) and (22) can be modified to instead use the experimental
boundary conditions (18).

######  Remark 3.4.

The non-linear system of equations that would be obtained by allowing the
conductivity to be complex in ultrasound modulated EIT (see e.g. [5]) is
similar to (22) with two real measurement equations per boundary condition
instead of a single one, i.e. for $x\in\Omega$:

$\displaystyle H_{jj}^{\prime}-\sigma^{\prime}(|\nabla
u_{j}^{\prime}|^{2}+|\nabla u_{j}^{\prime\prime}|^{2})$
$\displaystyle=0,~{}\text{and}$ $\displaystyle
H_{jj}^{\prime\prime}-\sigma^{\prime\prime}(|\nabla
u_{j}^{\prime}|^{2}+|\nabla u_{j}^{\prime\prime}|^{2})$ $\displaystyle=0.$

We did not consider this problem because the form of the measurements we
consider (22) is a direct result of using thermal induced random currents, see
Theorem 3.1.

##  4 Linearized problem

Before attempting to reconstruct conductivities numerically, we analyze the
linearizations of the real (21) and complex (22) conductivity problems. Our
goal is to find sufficient conditions for injectivity of the linearized
problems, or in other words, if they admit a unique solution. Our analysis is
based on [6], which includes a proof that the linearized real conductivity
problem is elliptic in the sense of Douglis-Nirenberg under certain boundary
conditions [18]. This was established in [26] for ultrasound modulated EIT and
generalized to other hybrid inverse problem in [6].

The linearization of the real conductivity problem (21) around the solution
$(u_{i},\sigma)$ in the variables $(\delta u_{i},\delta\sigma)$ for $i=1,..,n$
is given by

(23) $\displaystyle\nabla\cdot\left[{\sigma\nabla\delta
u_{i}}\right]+\nabla\cdot\left[{\delta\sigma\nabla u_{i}}\right]$
$\displaystyle=0,\quad x\in\Omega,$ $\displaystyle\delta u_{i}$
$\displaystyle=0,\quad x\in\partial\Omega,$ $\displaystyle\delta
H_{ii}-\delta\sigma|\nabla u_{i}|^{2}-2\sigma\nabla\delta u_{i}\cdot\nabla
u_{i}$ $\displaystyle=0,\quad x\in\Omega.$

The linearization of the complex conductivity problem (22) around the solution
$(u_{j}^{\prime},u_{j}^{\prime\prime},\sigma^{\prime},\sigma^{\prime\prime})$
in the variables $(\delta u_{j}^{\prime},\delta
u_{j}^{\prime\prime},\delta\sigma^{\prime},\delta\sigma^{\prime\prime})$ for
$j=1,..,n$ is given by

(24) $\displaystyle\nabla\cdot\left[{\sigma^{\prime}\nabla\delta
u_{j}^{\prime}}\right]+\nabla\cdot\left[{\delta\sigma^{\prime}\nabla
u_{j}^{\prime}}\right]-\nabla\cdot\left[{\sigma^{\prime\prime}\nabla\delta
u_{j}^{\prime\prime}}\right]-\nabla\cdot\left[{\delta\sigma^{\prime\prime}\nabla
u_{j}^{\prime\prime}}\right]$ $\displaystyle=0,\quad x\in\Omega,$
$\displaystyle\nabla\cdot\left[{\sigma^{\prime}\nabla\delta
u_{j}^{\prime\prime}}\right]+\nabla\cdot\left[{\delta\sigma^{\prime}\nabla
u_{j}^{\prime\prime}}\right]+\nabla\cdot\left[{\sigma^{\prime\prime}\nabla\delta
u_{j}^{\prime}}\right]+\nabla\cdot\left[{\delta\sigma^{\prime\prime}\nabla
u_{j}^{\prime}}\right]$ $\displaystyle=0,\quad x\in\Omega,$
$\displaystyle\delta u_{j}^{\prime}$ $\displaystyle=0,\quad
x\in\partial\Omega,$ $\displaystyle\delta u_{j}^{\prime\prime}$
$\displaystyle=0,\quad x\in\partial\Omega,$ $\displaystyle\delta
H_{jj}-\delta\sigma^{\prime}(|\nabla u_{j}^{\prime}|^{2}+|\nabla
u_{j}^{\prime\prime}|^{2})-2\sigma^{\prime}\nabla\delta
u_{j}^{\prime}\cdot\nabla u_{j}^{\prime}-2\sigma^{\prime}\nabla\delta
u_{j}^{\prime}\cdot\nabla u_{j}^{\prime}$ $\displaystyle=0,\quad x\in\Omega.$

In [6], it is established (23) is elliptic using two boundary conditions if
the gradients of the associated fields are nowhere orthogonal or parallel. We
do not attempt to analyze how this condition might be satisfied in the case of
the mixed boundary conditions (18). We note that this establishes that (23) is
not elliptic in the case of one experiment with mixed boundary conditions.
Instead, we attempt to analyze the problem numerically by estimating the
conditioning of the symbol of the linearized problem. For the case of complex
conductivity (24), we give a sufficient condition in Lemma 4.1 for
ellipticity; however, we do not give boundary conditions that guarantee this
is satisfied, nor do we prove that such boundary conditions exist. We note for
elliptic linear systems, it is possible to obtain stability estimates by
augmenting the system with boundary conditions satisfying the Lopatinskii
condition, following [6].

To establish if (23) and (24) are elliptic, we first compute the principal
symbol of their associated matrix-valued differential operators
$\mathcal{A}(x,D)$ for $x\in\Omega$, where
$D=(\partial_{x_{1}},...,\partial_{x_{n}})$. Since these are linearized
systems, the entries $\mathcal{A}_{ij}(x,D)$ are polynomials in $D$ for each
$x\in\Omega$. We associate each row of $\mathcal{A}$ with an integer $s_{i}$
and each column with an integer $t_{j}$, chosen such that the maximum degree
of each polynomial $\mathcal{A}_{ij}(x,D)$ is $s_{i}+t_{j}$. The principal
component $\mathcal{A}_{0}(x,D)$ is obtained from $\mathcal{A}(x,D)$ by
keeping only the terms in $\mathcal{A}_{ij}(x,D)$ with order exactly
$s_{i}+t_{j}$. If the principal symbol $\mathcal{A}_{0}(x,\xi)$ is injective
for all $\xi\neq 0$, then the problem is elliptic in the Douglis-Nirenberg
sense at $x\in\Omega$.

###  4.1 Injectivity of the linearized real problem

Letting $F_{i}=\nabla u_{i}$ the principal symbol of the real problem (23) is
the $2n\times(n+1)$ matrix

(25) $\mathcal{A}_{0}(x,\xi)=\begin{bmatrix}|F_{1}|^{2}&2\sigma F_{1}\cdot
i\xi&\cdots&0\\\ F_{1}\cdot i\xi&-\sigma|\xi|^{2}&\cdots&0\\\
\vdots&\vdots&\ddots&\vdots\\\ |F_{n}|^{2}&0&\cdots&2\sigma F_{n}\cdot i\xi\\\
F_{n}\cdot i\xi&0&\cdots&-\sigma|\xi|^{2}\end{bmatrix},$

for $i=1,...,n$ and where $\sigma=\sigma^{\prime}$. The system is in Douglis-
Nirenberg form where the row weights $s_{i}$ are given by the $2n$ vector
$(0,1,0,1,\ldots,0,1)$ and the column weights $t_{j}$ are given by the $n+1$
vector $(0,1,1,\ldots,1)$. As noted previously, this symbol is shown to be
injective in two dimensions using two boundary conditions such that $F_{1}$
and $F_{2}$ are nowhere orthogonal, or parallel [6].

We consider the discretized problem on the square $[0,10]^{2}$ using a uniform
$200\times 200$ grid. The conductivity used can be seen in Fig. 4 (a). We
numerically solve (18) to calculate the fields $u_{i}$ for $i=1,..,n$. The
Dirichlet boundary conditions are defined on the set

(26)
$\displaystyle\Gamma=\big{(}([0,4.5]\cup[5.5,10])\times\\{0,10\\}\big{)}\cup\big{(}\\{0,10\\}\times([0,4.5]\cup[5.5,10])\big{)}.$

For $x\in\Gamma$ the boundary conditions are of the form

(27) $\displaystyle g_{n}$
$\displaystyle=5\sin(\theta)\left({\frac{r}{10}}\right)^{n},$ $\displaystyle
h_{n}$ $\displaystyle=5\cos(\theta)\left({\frac{r}{10}}\right)^{n},$

where $(r,\theta)$ is the polar representation of the nodes in $\Gamma$. We
enforce the no flux boundary condition on the gaps, i.e.
$x\in\partial\Omega-\Gamma$.

To establish the ellipticity of the operator, it is sufficient to show that it
is injective for all $\xi$ such that $|\xi|=1$. We check this condition
numerically for $\xi\in\Xi$, where $\Xi$ is a set of 100 vectors uniformly
spaced on the unit circle (since $\xi$ is two-dimensional in our simulations).
At each point $x$ in the grid we use to discretize $\Omega$, we compute the
maximum condition number of the symbol along directions $\xi\in\Xi$, i.e.

(28)
$\max_{\xi\in\Xi}\left[{\frac{\sigma_{\max}(\mathcal{A}_{0}(x,\xi))}{\sigma_{\min}(\mathcal{A}_{0}(x,\xi))}}\right],$

where $\sigma_{\min}(A)$ (resp. $\sigma_{\max}(A)$) is the smallest (resp.
largest) singular value of a matrix $A$. The maximum condition number (28) of
the symbol (25) can be seen in Fig. 2. These numerical results are in line
with the previously established theory: the maximum condition number is higher
under one experiment than under two experiments. Even under one experiment,
the conditioning is still reasonable for most applications. Indeed, the
conductivity can still be “reasonably” recovered using one experiment but with
worsened numerical artifacts.

|
---|---
(a) One boundary condition | (b) Two boundary conditions
Figure 2: Maximum condition number (28) of the symbol of the linearized
problem (23) on the square domain $[0,10]^{2}$ with $\log 10$ scaling. The
left image (a) is the conditioning of the symbol with one boundary condition
given by $g_{1}$ in (27). The right image (b) is the conditioning of the
symbol with two boundary conditions given by $g_{1}$ and $h_{1}$ in (27).

###  4.2 Injectivity of the linearized complex problem

Letting $F_{j}^{\prime}=\nabla u_{j}^{\prime}$, $F_{j}^{\prime\prime}=\nabla
u_{j}^{\prime\prime}$, and $F_{j}=F_{j}^{\prime}+iF_{j}^{\prime\prime}$, the
symbol for the complex linear system (24) is the $3n\times(2+2n)$ matrix

(29)
$\widetilde{\mathcal{A}}_{0}(x,\xi)=\begin{bmatrix}|F_{1}|^{2}&0&2\sigma^{\prime}F_{1}^{\prime}\cdot
i\xi&2\sigma^{\prime}F_{1}^{\prime\prime}\cdot i\xi&\cdots&0&0\\\
F_{1}^{\prime}\cdot i\xi&-F_{1}^{\prime\prime}\cdot
i\xi&-\sigma^{\prime}|\xi|^{2}&\sigma^{\prime\prime}|\xi|^{2}&\cdots&0&0\\\
F_{1}^{\prime\prime}\cdot i\xi&F_{1}^{\prime}\cdot
i\xi&-\sigma^{\prime\prime}|\xi|^{2}&-\sigma^{\prime}|\xi|^{2}&\cdots&0&0\\\
\vdots&\vdots&\vdots&\ddots&\ddots&\vdots&\vdots\\\
|F_{n}|^{2}&0&\cdots&0&\cdots&2\sigma^{\prime}F_{n}^{\prime}\cdot
i\xi&2\sigma^{\prime}F_{n}^{\prime\prime}\cdot i\xi\\\ F_{n}^{\prime}\cdot
i\xi&-F_{n}^{\prime\prime}\cdot
i\xi&\cdots&0&\cdots&-\sigma^{\prime}|\xi|^{2}&\sigma^{\prime\prime}|\xi|^{2}\\\
F_{n}^{\prime\prime}\cdot i\xi&F_{n}^{\prime}\cdot
i\xi&\cdots&0&\cdots&-\sigma^{\prime\prime}|\xi|^{2}&-\sigma^{\prime}|\xi|^{2}\end{bmatrix},$

for $j=1,...,n$. The system is in Douglis-Nirenberg form where the row weights
$s_{i}$ are given by the $3n$ vector $(0,1,1,0,1,1,\ldots,0,1,1)$ and the
column weights $t_{j}$ are given by the $2+2n$ vector $(0,0,1,1,\ldots,1,1)$.

######  Lemma 4.1.

The symbol of the system (24), $\widetilde{\mathcal{A}}_{0}(x,\xi)$, with
$n\geq 2$ is injective at $x\in\Omega$ if there exists two experiments $i$ and
$j$ such that for all $\xi\neq 0$,

(30) $|F_{i}|^{2}|F_{j}\cdot\xi|^{2}\neq|F_{j}|^{2}|F_{i}\cdot\xi|^{2},$

and $\sigma^{\prime}(x),\sigma^{\prime\prime}(x)>0$.

######  Proof 4.2.

If $n=2$, we write the system such that the measurement equations are the
first two rows

$\begin{bmatrix}|F_{1}|^{2}&0&2\sigma^{\prime}F_{1}^{\prime}\cdot
i\xi&2\sigma^{\prime}F_{1}^{\prime\prime}\cdot i\xi&0&0\\\
|F_{2}|^{2}&0&0&0&2\sigma^{\prime}F_{2}^{\prime}\cdot
i\xi&2\sigma^{\prime}F_{2}^{\prime\prime}\cdot i\xi\\\ F_{1}^{\prime}\cdot
i\xi&-F_{1}^{\prime\prime}\cdot
i\xi&-\sigma^{\prime}|\xi|^{2}&\sigma^{\prime\prime}|\xi|^{2}&0&0\\\
F_{1}^{\prime\prime}\cdot i\xi&F_{1}^{\prime}\cdot
i\xi&-\sigma^{\prime\prime}|\xi|^{2}&-\sigma^{\prime}|\xi|^{2}&0&0\\\
F_{2}^{\prime}\cdot i\xi&-F_{2}^{\prime\prime}\cdot
i\xi&0&0&-\sigma^{\prime}|\xi|^{2}&\sigma^{\prime\prime}|\xi|^{2}\\\
F_{2}^{\prime\prime}\cdot i\xi&F_{2}^{\prime}\cdot
i\xi&0&0&-\sigma^{\prime\prime}|\xi|^{2}&-\sigma^{\prime}|\xi|^{2}\end{bmatrix}.$

We consider this as a block matrix with the top left block being $2\times 2$,
and then the bottom right block being $4\times 4$. The bottom right matrix is
block diagonal with invertible diagonal $2\times 2$ blocks, and we use this to
compute the Schur complement

$\begin{bmatrix}|F_{1}|^{2}+\frac{2(\sigma^{\prime})^{2}}{(|\xi|\sigma^{\prime})^{2}+(|\xi|\sigma^{\prime\prime})^{2}}\left({(F_{1}^{\prime}\cdot
i\xi)^{2}+(F_{1}^{\prime\prime}\cdot
i\xi)^{2}}\right)&\frac{2\sigma^{\prime\prime}\sigma^{\prime}}{(|\xi|\sigma^{\prime})^{2}+(|\xi|\sigma^{\prime\prime})^{2}}\left({(F_{1}^{\prime}\cdot
i\xi)^{2}+(F_{1}^{\prime\prime}\cdot i\xi)^{2}}\right)\\\
|F_{2}|^{2}+\frac{2(\sigma^{\prime})^{2}}{(|\xi|\sigma^{\prime})^{2}+(|\xi|\sigma^{\prime\prime})^{2}}\left({(F_{2}^{\prime}\cdot
i\xi)^{2}+(F_{2}^{\prime\prime}\cdot
i\xi)^{2}}\right)&\frac{2\sigma^{\prime\prime}\sigma^{\prime}}{(|\xi|\sigma^{\prime})^{2}+(|\xi|\sigma^{\prime\prime})^{2}}\left({(F_{2}^{\prime}\cdot
i\xi)^{2}+(F_{2}^{\prime\prime}\cdot i\xi)^{2}}\right)\end{bmatrix}.$

The determinant of the Schur complement is then

(31)
$\frac{2\sigma^{\prime\prime}\sigma^{\prime}}{(\sigma^{\prime})^{2}+(\sigma^{\prime\prime})^{2}}\left({\frac{|F_{1}|^{2}}{|\xi|^{2}}\left[{(F_{2}^{\prime}\cdot
i\xi)^{2}+(F_{2}^{\prime\prime}\cdot
i\xi)^{2}}\right]-\frac{|F_{2}|^{2}}{|\xi|^{2}}\left[{(F_{1}^{\prime}\cdot
i\xi)^{2}+(F_{1}^{\prime\prime}\cdot i\xi)^{2}}\right]}\right),$

which gives the desired result.

For $n>2$, without lost of generality we let $i=1$ and $j=2$ in (30). We
proceed by using the same approach considering the Schur complement of the
largest square sub-matrix obtained by deleting the measurement terms of all
experiments $i>2$, that is, the $(2+2n)\times(2+2n)$ matrix,

$\begin{bmatrix}|F_{1}|^{2}&0&2\sigma^{\prime}F_{1}^{\prime}\cdot
i\xi&2\sigma^{\prime}F_{1}^{\prime\prime}\cdot i\xi&0&0&\cdots&0&0\\\
|F_{2}|^{2}&0&0&0&2\sigma^{\prime}F_{2}^{\prime}\cdot
i\xi&2\sigma^{\prime}F_{2}^{\prime\prime}\cdot i\xi&\cdots&0&0\\\
F_{1}^{\prime}\cdot i\xi&-F_{1}^{\prime\prime}\cdot
i\xi&-\sigma^{\prime}|\xi|^{2}&\sigma^{\prime\prime}|\xi|^{2}&0&0&\cdots&0&0\\\
F_{1}^{\prime\prime}\cdot i\xi&F_{1}^{\prime}\cdot
i\xi&-\sigma^{\prime\prime}|\xi|^{2}&-\sigma^{\prime}|\xi|^{2}&0&0&\cdots&0&0\\\
F_{2}^{\prime}\cdot i\xi&-F_{2}^{\prime\prime}\cdot
i\xi&0&0&-\sigma^{\prime}|\xi|^{2}&\sigma^{\prime\prime}|\xi|^{2}&\cdots&0&0\\\
F_{2}^{\prime\prime}\cdot i\xi&F_{2}^{\prime}\cdot
i\xi&0&0&-\sigma^{\prime\prime}|\xi|^{2}&-\sigma^{\prime}|\xi|^{2}&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\ddots&\ddots&\vdots&\vdots\\\
F_{n}^{\prime}\cdot i\xi&-F_{n}^{\prime\prime}\cdot
i\xi&0&0&0&0&\cdots&-\sigma^{\prime}|\xi|^{2}&\sigma^{\prime\prime}|\xi|^{2}\\\
F_{n}^{\prime\prime}\cdot i\xi&F_{n}^{\prime}\cdot
i\xi&0&0&0&0&\cdots&-\sigma^{\prime\prime}|\xi|^{2}&-\sigma^{\prime}|\xi|^{2}\end{bmatrix}.$

We consider the top left $2\times 2$ matrix as a block, and the bottom right
$2n\times 2n$ matrix as a block. The bottom right matrix is still block
diagonal with invertible $2\times 2$ blocks. The determinant of the Schur
complement is given by (31), and the determinant of the sub-matrix is given by

(32)
$2\sigma^{\prime\prime}\sigma^{\prime}|\xi|^{4n-2}|\sigma|^{2(n-1)}\left({|F_{1}|^{2}|F_{2}\cdot\xi|^{2}-|F_{2}|^{2}|F_{1}\cdot\xi|^{2}}\right),$

so by supposition, the sub-matrix is invertible. Since one of the largest
square sub-matrices of the symbol is invertible, the symbol is injective.

We do not present a method for finding boundary conditions such that the
condition in Lemma 4.1 is satisfied, nor do we know if such boundary
conditions exist for all possible
$\sigma^{\prime}(x),\sigma^{\prime\prime}(x)>0$. If $F_{i}$ and $F_{j}$ are
unit length then the condition here resembles that in [6], giving that the
fields cannot be orthogonal or parallel. We note the following corollary,
which is limited to the case where $n=2$.

######  Corollary 4.3.

The symbol of the system (24), $\widetilde{\mathcal{A}}_{0}(x,\xi)$, with
$n=2$ is not injective at $x\in\Omega$ if $\sigma^{\prime}(x)=0$ or
$\sigma^{\prime\prime}(x)=0$.

######  Proof 4.4.

This follows by observing that the determinant of the Schur complement (31) is
zero if $\sigma^{\prime}(x)=0$ or $\sigma^{\prime\prime}(x)=0$.

Showing that a system is not elliptic for all $n$ using determinants is
challenging because we need to ensure that the determinants of all maximal
square sub-matrices are zero. However this approach is sufficient for the case
$n=2$ and we expect Corollary 4.3 to hold with more experiments.

Numerically, we have found that the complex reconstruction is challenging
under many combinations of boundary conditions. This can be expected from
numerically computing the maximum condition number (28) as we illustrate in
the numerical experiment appearing in Fig. 3 and that we describe next.

|  |
---|---|---
(a) Two measurements | (b) Three measurements | (c) Four measurements
Figure 3: The maximum condition number (28) of the symbol of the complex
linearized problem (24) on the discretized square domain with $\log 10$
scaling. The images from left to right are given by two, three, and four
boundary conditions in the form of Eq. 33.

The ground truth conductivity can be seen in Fig. 6 with the real part in (a)
and the imaginary part in (b). We only consider the linearized complex problem
(24) with Dirichlet boundary conditions. The boundary conditions are of the
form

(33) $\displaystyle\tilde{g}_{n}$ $\displaystyle=g_{n}+\frac{i}{2}h_{n},$
$\displaystyle\tilde{h}_{n}$ $\displaystyle=h_{n}+\frac{i}{2}g_{n},$

with $h$ and $g$ defined in (27). The scaling of the imaginary part by $1/2$
is to match the imaginary part of the background conductivity. In Fig. 3, we
can see the numerical condition of the symbol for two, three, and four
boundary conditions. We begin this experiment with $n=2$ since, with one
measurement, the system is underdetermined. The boundary conditions for $n=2$
are $\tilde{g}_{1},\tilde{h}_{1}$, for $n=3$ are
$\tilde{g}_{1},\tilde{g}_{2},\tilde{h}_{1}$, and for $n=4$ are
$\tilde{g}_{1},\tilde{g}_{2},\tilde{h}_{1},\tilde{h}_{2}.$ We use the same
domain and grid as in the real case.

The maximum condition number (28) improves significantly by moving from two to
three measurements, but the improvement from three to four is modest. The
areas with high condition number for two boundary conditions (Fig. 3 (a))
suggest the problem is not elliptic. With more boundary conditions (Fig. 3 (b)
and (c)), the areas where the conditioning is high match up with the
reconstruction artifacts for the complex case (Fig. 6).

##  5 Numerical reconstructions

The following numerical reconstructions use values consistent with the quasi-
static approximation, i.e., values such that $\omega\mu|\sigma|L^{2}\ll 1$. In
particular we let $L=10$ cm and $\sigma^{\prime}\in[1/3,2]$ cm-1
$k\Omega^{-1}$. In the case of non-zero complex conductivity we let
$\omega=2\pi 10$ kHz and
$\sigma^{\prime\prime}=\omega\epsilon^{\prime}\in[1/2,1]$ cm-1 $k\Omega^{-1}$.
Our choice of parameters is near those in human tissues and satisfies the
quasi-static approximation, see e.g. [13]. The examples we consider assume a
thin plate that is homogeneous in the $z$ direction with thickness $\Delta
z=0.1$cm. If we consider $\Omega\subset\mathbb{R}^{2}$, then multiplying the
measurements by $\Delta z$ corresponds to the results in Section 3.

###  5.1 Discrete model

We discretize the system (21) on a square domain $\Omega=[0,10]^{2}$ using a
uniform grid with $n^{2}$ nodes. We denote by $N$ the set of nodes indexed
with their integer coordinates $(i,j)$ and the set of edges by $E\subset
N\times N$. The nodes are partitioned into interior nodes $I$ and boundary
nodes $B$, which are the nodes that are on the boundary $\partial\Omega$. We
use the forward difference operator $D\in\mathbb{R}^{|N|\times|E|}$ defined
such that

$D=\begin{bmatrix}D_{1}\\\ D_{2}\end{bmatrix}.$

Here $D_{1}$ (resp. $D_{2}$) is the horizontal (resp. vertical) first order
difference operator. Given a function $\psi$ defined on the nodes $N$, the
horizontal and vertical difference operators are defined by

$\displaystyle(D_{1}\psi)(i,j)$
$\displaystyle=\frac{\psi(i+1,j)-\psi(i,j)}{\Delta x},$
$\displaystyle(D_{2}\psi)(i,j)$
$\displaystyle=\frac{\psi(i,j+1)-\psi(i,j)}{\Delta y},$

where $\Delta x$ and $\Delta y$ are the horizontal and vertical discretization
steps respectively.

If we use finite differences to discretize (21), we note that the gradient
components in the $x$ and $y$ directions are defined on horizontal and
vertical edges. Thus the norm of the discretized gradient is not defined at
any particular edge. To obtain the gradient at a single spatial location in
the discretized problem, we interpolate the gradient approximated values from
their respective edges to the nodes and compute the gradient norm at the
nodes. Thus it also makes sense to interpret the internal functional $H_{ii}$
as a nodal based quantity and to completely determine the conductivity by
interpolating a node based quantity. These interpolations between edges and
nodes are achieved with the following matrices

$\displaystyle N_{1}:$ $\displaystyle\text{horizontal edges}\to\text{nodes},$
$\displaystyle N_{2}:$ $\displaystyle\text{vertical edges}\to\text{nodes},$
$\displaystyle E_{1}:$ $\displaystyle\text{nodes}\to\text{horizontal edges},$
$\displaystyle E_{2}:$ $\displaystyle\text{nodes}\to\text{vertical edges},$
$\displaystyle E_{1,2}:$ $\displaystyle\text{nodes}\to\text{all edges}.$

To define these matrices we let $\phi$ be the matrix-valued function

$\phi(A)=(\text{Diag}(|A^{T}|\mathbf{1}))^{-1}|A^{T}|,$

where $|\cdot|$ is the entry-wise absolute value, $\mathbf{1}$ is an
appropriately sized vector of ones, and $\text{Diag}(v)$ denotes the matrix
with the vector $v$ on its diagonal. The matrix $\phi(A)$ preserves constant
vectors, more precisely, if $c$ is an appropriate sized constant vector
$\phi(A)c=c$. The interpolation operators are then defined as

$N_{1}=\phi(D_{1}^{T}),\quad N_{2}=\phi(D_{2}^{T}),\quad
E_{1}=\phi(D_{1}),\quad E_{2}=\phi(D_{2}),\quad
E_{1,2}=\begin{bmatrix}E_{1}\\\ E_{2}\end{bmatrix}.$

Given $H_{ii}$ (defined at the nodes) and $e_{i}$ (defined at the boundary
nodes) the discrete inverse problem for real conductivity is then to find $s$
and $u_{i}$ (defined at the nodes) such that for $i=1,...,N$,

(34) $\displaystyle D^{T}[E_{1,2}s\odot(Du_{i})]_{I}=0,$ $\displaystyle
u_{i}|_{B}-e_{i}=0,$ $\displaystyle
H_{ii}-\left[{N_{1}(E_{1}s\odot|D_{1}u_{i}|^{2})+N_{2}(E_{2}s\odot|D_{2}u_{i}|^{2})}\right]_{I}=0,$

where $\odot$ is the Hadamard or componentwise product. We note that in the
first equation of (34), we have a graph Laplacian with edge weights given by
$E_{1,2}s$, see e.g. [16]. This system is modified slightly under the
assumption that the conductivity is known in a small neighborhood of the
boundary. The modified system is solved using Gauss-Newton iteration. However,
the interpolation process introduces a null space into the Jacobian. We use
Tikhonov regularization with a parameter $\gamma$ to prevent this null space
from interfering when solving for the Gauss-Newton step. The parameter
corresponds to adding a penalty term of $\gamma\|w\|^{2}$, when solving the
least squares problem for finding a Gauss-Newton step $w$. An Armijo line
search is used as globalization strategy (see e.g. [33]).

Given $H_{jj}$ and $e_{j}$, the complex inverse problem is to recover
$s^{\prime}$, $s^{\prime\prime}$, $u_{j}^{\prime}$, and $u_{j}^{\prime\prime}$
for $j=1,...,N$,

(35) $\displaystyle
D^{T}[E_{1,2}s^{\prime}\odot(Du_{j}^{\prime})]_{I}-D^{T}[E_{1,2}s^{\prime\prime}\odot(Du_{j}^{\prime\prime})]_{I}=0,$
$\displaystyle
D^{T}[E_{1,2}s^{\prime}\odot(Du_{j}^{\prime\prime})]_{I}+D^{T}[E_{1,2}s^{\prime\prime}\odot(Du_{j}^{\prime})]_{I}=0,$
$\displaystyle u_{j}^{\prime}|_{B}-e_{j}^{\prime}=0,$ $\displaystyle
u_{j}^{\prime\prime}|_{B}-e_{j}^{\prime\prime}=0,$ $\displaystyle
H_{jj}-\left[{N_{1}(E_{1}s^{\prime}\odot|D_{1}u_{j}^{\prime}|^{2})+N_{2}(E_{2}s\odot|D_{2}u_{j}^{\prime}|^{2})}\right]_{I}$
$\displaystyle+\left[{N_{1}(E_{1}s^{\prime}\odot|D_{1}u_{j}^{\prime\prime}|^{2})+N_{2}(E_{2}s^{\prime}\odot|D_{2}u_{j}^{\prime\prime}|^{2})}\right]_{I}=0.$

A similar Gauss-Newton procedure was used to solve (35).

Heating Patterns: Recall that to obtain the measurements $H_{ii}$; we need to
locally heat a region of the conducting plate. This requires numerically
approximating (19) at both the background temperature $T_{0}$ and when the
plate is locally heated according to $\delta T(x)$. For our measurements we
let $\delta T(x)$ be the Gaussian heating pattern,

$g(x,a)=(2\pi a)^{-1}\exp\left({-|x|^{2}(2a)^{-1}}\right)$

where $|\cdot|$ denotes the 2-norm. The heating pattern $g(x,a)$ can be
considered an approximate Dirac since it integrates in $x$ to one. This is
also similar to the heating pattern from a laser covering an area of roughly
$\pi a$.

Deterministic measurements: Continuous measurements from the deterministic
model using this heating pattern can be written as

(36) $H_{ii}(x)=\left\langle{\sigma|\nabla
u_{i}|^{2},T_{0}+g(\cdot-x,a)}\right\rangle_{L_{2}(\Omega)}-\left\langle{\sigma|\nabla
u_{i}|^{2},T_{0}}\right\rangle_{L_{2}(\Omega)}.$

We approximate $H_{ii}(x)$ by evaluating the heating pattern at each node $x$
in the discrete model. Then to approximate the inner products, we use a
uniform fine grid with $\tilde{n}^{2}$ nodes such that $\tilde{n}>n$. The
number of fine grid nodes $\tilde{n}$ is chosen such that there are at least
four fine grid nodes per effective area of the heating pattern, i.e. $\pi a$.

Stochastic measurements: The measurements from the simulated random current
model are given by approximating

(37) $H_{ii}(x)=\left\langle{\left\langle{\sigma|\nabla
u_{i}|^{2},T_{0}+g(\cdot-x,a)}\right\rangle_{L_{2}(\Omega)}}\right\rangle-\left\langle{\left\langle{\sigma|\nabla
u_{i}|^{2},T_{0}}\right\rangle_{L_{2}(\Omega)}}\right\rangle,$

where the outer angular brackets denote ensemble averaging. The average is
approximated empirically with $M$ realizations of random currents where the
inner product for each realization is approximated using a uniform fine grid.
The realizations of the background temperature measurements (the rightmost
term in (37)) are not recalculated for each heating pattern. For each
realization, the currents at every fine grid edge midpoint $e$ are sampled
from a mean zero random normal distribution with standard deviation
$\sqrt{\kappa(T_{0}+g(e,a))s(e)/\pi}$ (heated) or $\sqrt{\kappa
T_{0}s(e)/\pi}$ (unheated) as determined by (4).

We believe the simulation method used for realizations of random currents
creates a more challenging problem than experimental data. In practice,
measurements could be taken over a time interval and then averaged over time.
This gives temporal structure to the data that is not reflected by our
simulations which ignore the temporal correlation of the random currents.

###  5.2 Real conductivity

First, we consider problems of purely real conductivity
($\sigma=\sigma^{\prime}$) in both the case of measurements from simulated
random currents (stochastic model) and measurements using their variances
(deterministic model). For both problems we consider the Dirichlet boundary
conditions (11) using $e_{1}=(x_{1}+x_{2})/10$ and $e_{2}=(1+x_{1}-x_{2})/10$.
For the more challenging problem of the experimental boundary conditions (18),
we only consider the deterministic model of measurements.

We show in Fig. 4 numerical reconstructions for a purely real conductivity in
the deterministic and stochastic models. The deterministic measurements
$H_{ii}$ in (36) are taken with $T_{0}=300$ and $a=0.01$ on a $60\times 60$
coarse grid using a $120\times 120$ fine grid to approximate the integrals.
The same conditions were used for the stochastic model, where in addition we
took $T_{0}=0.01$ using 1000 realizations of random currents. We chose a
particularly low background temperature $T_{0}$ to get clean enough data with
the number of realizations we chose. As can be expected from Section 3.3, the
signal to noise ration for the differential temperature measurements worsens
for large $T_{0}$. We mention that our choice of discretization ensures that
each approximate Dirac heating pattern ($a=0.01$) covers a minimum of four
grid points.

The conductivity is assumed to be known for nodes that are $0.5$cm, or less,
away from the boundary. The Tikhonov regularization parameter is
$\gamma=5^{-4}$, and iterations are run until the 2-norm of the step is less
than $0.1$. The initial guess is the solution to (34) with a constant
conductivity.

|  |
---|---|---
(a) | (b) | (c)
Figure 4: Reconstructions of purely real conductivity values (cm-1
$k\Omega^{-1}$) on a $10$cm$\times$10cm square domain. The ground truth
conductivity in (a) is evaluated on the fine grid. Both the reconstructions
using the deterministic model (b) and a stochastic model (c) are evaluated on
the coarse grid.

The numerical example in Fig. 5 uses data from experimental boundary
conditions in (18).The set $\Gamma$ defined in (26), corresponds to the
electrode functions. On $\Gamma$, we use electrode functions $g_{1}$ and
$h_{1}$ (27) for the boundary conditions. This set has gaps of $1$cm at the
center of each side of the square with no flux conditions. The no flux
conditions are enforced by using centered approximations to the nodes on
$\partial\Omega-\Gamma$ (see, e.g., [29, sec 2.12]). The ground truth
conductivity is given in Fig. 5 (a). The reconstructions are evaluated on a
$100\times 100$ coarse grid, and a $200\times 200$ fine grid is used to
evaluate the measurements (36). A minimum of twelve fine grid points are in
the effective area of each approximate Dirac heating pattern.

The Gauss-Newton iteration is regularized with $\gamma=3^{-3}$, and iterations
are run until the 2-norm of the step size is less than $0.1$. The
reconstructions in Fig. 5 are close to the original conductivity, although
there are some numerical artifacts due to the gaps between the electrodes.

|  |
---|---|---
(a) | (b) | (c)
Figure 5: Conductivity values (cm-1 $k\Omega^{-1}$) for a numerical
reconstruction (c) of the same conductivity as Fig. 4 using the experimental
boundary conditions in Eq. 18. The data $H_{ii}=\text{Re}\,(\sigma(x))|\nabla
u_{i}(x)|^{2}$ used for the reconstructions is shown in (a) for $g_{1}$ and in
(b) for $h_{1}$, with units cm-3 $k\Omega^{-1}$.

###  5.3 Complex conductivity

An example of a complex conductivity reconstruction using the deterministic
model of measurements can be seen in Fig. 6. Four experiments are used with
the Dirichlet boundary conditions
$\tilde{g}_{1},\tilde{g}_{2},\tilde{h}_{1},\tilde{h}_{2}$ given in (33). The
reconstructions are evaluated on a $100\times 100$ coarse grid, and a
$200\times 200$ fine grid is used to evaluate the measurements (36). The
Gauss-Newton iteration now uses $\gamma=1^{-4}$, and iterations are run until
the 2-norm of the step size is less than $0.1$. A constant complex
conductivity and the corresponding solutions to (35) are used for the initial
guess.

The numerical artifacts in the complex reproduction are consistent with the
areas in Fig. 3 (c), where the maximum condition number (28) of the symbol is
largest. Intuitively, reconstructing the imaginary conductivity may be more
challenging as it does not explicitly appear in the measurements. In the
discrete system (35) it appears only when coupled with a gradient of the real
or imaginary auxiliary field. When the gradient of the real part is large it
may overwhelm the contribution of the complex conductivity. This is the exact
behavior we see in our numerical experiments. When the real conductivity is
high, the gradient of the real field is large, and the reconstructed complex
conductivity (Fig. 6 (d)) is lower than the true value.

|
---|---
(a) True $\sigma^{\prime}$ | (b) True $\epsilon\omega$
|
(c) Reconstructed $\sigma^{\prime}$ | (d) Reconstructed $\epsilon\omega$
Figure 6: Real and imaginary conductivity values (cm-1 $k\Omega^{-1}$) for a
numerical reconstruction (c) & (d) of a complex conductivity on a
$10$cm$\times$10cm square domain. The ground truth (a) & (b) is evaluated on
the fine grid and the reconstruction is evaluated on the coarse grid.

##  6 Summary and perspectives

We propose a new hybrid inverse problem for recovering the conductivity of a
body using thermal noise. The fluctuation dissipation theorem for
electrodynamic media allows us to relate the variance of thermal noise
currents taken with different temperature patterns to the real part of the
conductivity of a body. By taking a sufficiently rich set of measurements we
can estimate an internal functional that depends on this real conductivity and
the solution to an associated auxiliary problem. We show this relation holds
for both Dirichlet boundary conditions and mixed Neumann/Dirichlet boundary
conditions, the latter of which is a more realistic description of an
experimental setup where such measurements might be used. For purely real
conductivities, these are power density measurements. This problem of
recovering a real conductivity from power density measurements also appears in
acousto-electric tomography.

Before attempting numerical reconstructions we try and determine if the
linearized problems are elliptic in the Douglis-Nirenberg sense. The
linearized real problem has previously been shown to be elliptic, given the
auxiliary fields are nowhere orthogonal or parallel [6]. For the real problem
with mixed boundary conditions, we make no effort to show our boundary
conditions satisfy this condition. Instead, we numerically evaluate the worst
conditioning of the principal symbol at each point in space. The numerical
results are consistent with the previous theoretical work in [6]; using more
boundary conditions improves the numerical conditioning of the symbol. For the
complex symbol, we give a sufficient condition on the auxiliary fields for the
problem to be elliptic, without proving that the auxiliary fields can be
generated. We perform a similar numerical evaluation of the conditioning of
the symbol under a number of different experiments. This evidence indicates
that the complex conductivity problem with two boundary conditions is not
elliptic. This numerical approach to classification may find use in similar
problems, especially in problems with complicated boundary conditions or
principal symbols. The clear limitation of this method is that it does not
inform the choice of boundary conditions. We note that the conditioning in
these problems would not normally be seen as high for other applications.

Finally, we present a simple discrete model for numerical reconstructions. Our
numerical reconstructions are consistent with the linearization study. We also
present results using simulated random thermal currents for the case of a
purely real conductivity. These simulations ignore temporal correlations,
making the problem more challenging. In this case, we can get an accurate, if
noisy, reconstruction of the conductivity at a low temperature.

This method of thermal noise imaging may find applications in e.g. Atomic
Force Microscopy, laser weld monitoring. A challenge in using our approach is
that the relative size of the measurements due to the background temperature
and the heating pattern (see Section 3.3) results in currents that may hard to
measure reliably in practice.

Our results are for a fixed frequency $\omega$, and removing this limitation
may allow for more accurate reconstructions in the complex case. Considering
multiple frequencies means that $\epsilon$ will be treated as a variable
instead of $\epsilon\omega$. This is an important difference because it
changes the structure of the symbol and requires a separate analysis.
Additionally, the relative scale of the variables of interest, $\omega$ and
$\epsilon$, changes which may introduce other challenges to the
reconstructions.

## Acknowledgments

The authors would like to thank Maxence Cassier for useful comments on a draft
of this manuscript.

## References

  * [1] H. Ammari, An introduction to mathematics of emerging biomedical imaging, vol. 62 of Mathématiques & Applications (Berlin) [Mathematics & Applications], Springer, Berlin, 2008.
  * [2] H. Ammari, E. Bonnetier, Y. Capdeboscq, M. Tanter, and M. Fink, Electrical impedance tomography by elastic deformation, SIAM J. Appl. Math., 68 (2008), pp. 1557–1573, https://doi.org/10.1137/070686408.
  * [3] H. Ammari, Y. Capdeboscq, F. de Gournay, A. Rozanova-Pierrat, and F. Triki, Microwave imaging by elastic deformation, SIAM J. Appl. Math., 71 (2011), pp. 2112–2130, https://doi.org/10.1137/110828241.
  * [4] G. Bal, Cauchy problem for ultrasound-modulated EIT, Anal. PDE, 6 (2013), pp. 751–775, https://doi.org/10.2140/apde.2013.6.751.
  * [5] G. Bal, Hybrid inverse problems and internal functionals, in Inverse problems and applications: inside out. II, vol. 60 of Math. Sci. Res. Inst. Publ., Cambridge Univ. Press, Cambridge, 2013, pp. 325–368.
  * [6] G. Bal, Hybrid inverse problems and redundant systems of partial differential equations, in Inverse problems and applications, vol. 615 of Contemp. Math., Amer. Math. Soc., Providence, RI, 2014, pp. 15–47, https://doi.org/10.1090/conm/615/12289.
  * [7] G. Bal, E. Bonnetier, F. Monard, and F. Triki, Inverse diffusion from knowledge of power densities, Inverse Probl. Imaging, 7 (2013), pp. 353–375, https://doi.org/10.3934/ipi.2013.7.353.
  * [8] G. Bal, C. Guo, and F. Monard, Linearized internal functionals for anisotropic conductivities, Inverse Probl. Imaging, 8 (2014), pp. 1–22, https://doi.org/10.3934/ipi.2014.8.1.
  * [9] G. Bal, W. Naetar, O. Scherzer, and J. Schotland, The Levenberg-Marquardt iteration for numerical inversion of the power density operator, J. Inverse Ill-Posed Probl., 21 (2013), pp. 265–280, https://doi.org/10.1515/jip-2012-0091.
  * [10] G. Bal and J. C. Schotland, Inverse Scattering and Acousto-Optics Imaging, Phys. Rev. Letters, 104 (2010), p. 043902, https://doi.org/10.1103/PhysRevLett.104.043902.
  * [11] G. Bal and G. Uhlmann, Reconstruction of coefficients in scalar second-order elliptic equations from knowledge of their solutions, Comm. Pure Appl. Math., 66 (2013), pp. 1629–1652, https://doi.org/10.1002/cpa.21453.
  * [12] G. Bal and T. Zhou, Hybrid inverse problems for a system of Maxwell’s equations, Inverse Problems, 30 (2014), p. 055013.
  * [13] L. Borcea, Electrical impedance tomography, Inverse problems, 18 (2002), p. R99.
  * [14] Y. Capdeboscq, J. Fehrenbach, F. de Gournay, and O. Kavian, Imaging by modification: numerical reconstruction of local conductivities from corresponding power density measurements, SIAM J. Imaging Sci., 2 (2009), pp. 1003–1030, https://doi.org/10.1137/080723521.
  * [15] M. Cheney, D. Isaacson, and J. C. Newell, Electrical impedance tomography, SIAM Rev., 41 (1999), pp. 85–101, https://doi.org/10.1137/S0036144598333613.
  * [16] F. R. K. Chung, Spectral graph theory, vol. 92 of CBMS Regional Conference Series in Mathematics, Published for the Conference Board of the Mathematical Sciences, Washington, DC; by the American Mathematical Society, Providence, RI, 1997.
  * [17] H. E. Cline and T. R. Anthony, Heat treating and melting material with a scanning laser or electron beam, Journal of Applied Physics, 48 (1977), pp. 3895–3900, https://doi.org/10.1063/1.324261.
  * [18] A. Douglis and L. Nirenberg, Interior estimates for elliptic systems of partial differential equations, Communications on Pure and Applied Mathematics, 8 (1955), pp. 503–538.
  * [19] L. C. Evans, Partial differential equations, vol. 19, of Graduate Studies in Mathematics, American Mathematical Society, second ed., 2010\.
  * [20] B. Gebauer and O. Scherzer, Impedance-acoustic tomography, SIAM J. Appl. Math., 69 (2008), pp. 565–576, https://doi.org/10.1137/080715123.
  * [21] J. B. Johnson, Thermal agitation of electricity in conductors, Physical review, 32 (1928), p. 97.
  * [22] W. P. King, B. Bhatia, J. R. Felts, H. J. Kim, B. Kwon, B. Lee, S. Somnath, and M. Rosenberger, Heated atomic force microscope cantilevers and their applications, Annual Review of Heat Transfer, 16 (2013), pp. 287–326, https://doi.org/10.1615/AnnualRevHeatTransfer.v16.100.
  * [23] R. Kubo, The fluctuation-dissipation theorem, Reports on Progress in Physics, 29 (1966), p. 255, https://doi.org/10.1088/0034-4885/29/1/306, https://dx.doi.org/10.1088/0034-4885/29/1/306.
  * [24] P. Kuchment and L. Kunyansky, Synthetic focusing in ultrasound modulated tomography, Inverse Probl. Imaging, 4 (2010), pp. 665–673, https://doi.org/10.3934/ipi.2010.4.665.
  * [25] P. Kuchment and L. Kunyansky, 2D and 3D reconstructions in acousto-electric tomography, Inverse Problems, 27 (2011), pp. 055013, 21, https://doi.org/10.1088/0266-5611/27/5/055013.
  * [26] P. Kuchment and D. Steinhauer, Stabilizing inverse problems by internal data, Inverse Problems, 28 (2012), p. 084007.
  * [27] L. D. Landau and E. M. Lifshitz, Course of theoretical physics. Vol. 5: Statistical physics, Translated from the Russian by J. B. Sykes and M. J. Kearsley. Second revised and enlarged edition, Pergamon Press, Oxford-Edinburgh-New York, 1968.
  * [28] M. Lanza, Conductive Atomic Force Microscopy: Applications in Nanomaterials, John Wiley & Sons, 2017.
  * [29] R. J. LeVeque, Finite difference methods for ordinary and partial differential equations, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2007, https://doi.org/10.1137/1.9780898717839, https://doi.org/10.1137/1.9780898717839.  Steady-state and time-dependent problems.
  * [30] F. Monard and G. Bal, Inverse anisotropic diffusion from power density measurements in two dimensions, Inverse Problems, 28 (2012), pp. 084001, 20, https://doi.org/10.1088/0266-5611/28/8/084001.
  * [31] F. Monard and G. Bal, Inverse diffusion problems with redundant internal information, Inverse Probl. Imaging, 6 (2012), pp. 289–313, https://doi.org/10.3934/ipi.2012.6.289.
  * [32] F. Monard and G. Bal, Inverse anisotropic conductivity from power densities in dimension $n\geq 3$, Comm. Partial Differential Equations, 38 (2013), pp. 1183–1207, https://doi.org/10.1080/03605302.2013.787089.
  * [33] J. Nocedal and S. J. Wright, Numerical optimization, Springer Series in Operations Research and Financial Engineering, Springer, New York, second ed., 2006.
  * [34] H. Nyquist, Thermal agitation of electric charge in conductors, Physical review, 32 (1928), p. 110.
  * [35] S. Powell, S. R. Arridge, and T. S. Leung, Gradient-based quantitative image reconstruction in ultrasound-modulated optical tomography: First harmonic measurement type in a linearised diffusion formulation, IEEE Transactions on Medical Imaging, 35 (2016), pp. 456–467, https://doi.org/10.1109/TMI.2015.2478742.
  * [36] S. M. Rytov, Y. A. Kravtsov, and V. I. Tatarskiĭ, Principles of statistical radiophysics. 2, Springer-Verlag, Berlin, 1988, https://doi.org/10.1007/978-3-642-61351-7, https://doi.org/10.1007/978-3-642-61351-7.  Correlation theory of random processes, Translated from the second Russian edition by Alexander P. Repyev [A. P. Rep’ev].
  * [37] S. M. Rytov, Y. A. Kravtsov, and V. I. Tatarskiĭ, Principles of statistical radiophysics. 3, Springer-Verlag, Berlin, 1989.  Elements of random fields, Translated from the second Russian edition by Alexander P. Repyev [A. P. Rep’ev].
  * [38] R. Zwanzig, Nonequilibrium statistical mechanics, Oxford University Press, New York, 2001.