Entanglement and interaction in a topological quantum walk

Alberto D. Verga1*, Ricardo Gabriel Elías2

1 Aix-Marseille Université, CPT, Campus de Luminy, case 907, 13288 Marseille,
France

2 Departamento de Física and CEDENNA, Universidad de Santiago de Chile, Avda.
Ecuador 3493, Santiago Chile.

*<EMAIL_ADDRESS>
## Abstract

We study the quantum walk of two interacting particles on a line with an
interface separating two topologically distinct regions. The interaction
induces a localization-delocalization transition of the edge state at the
interface. We characterize the transition through the entanglement between the
two particles.

###### Contents

  1. 1 Introduction
  2. 2 Model
  3. 3 Results
  4. 4 Discussion and conclusions

##  1 Introduction

Quantum walks are not only important in quantum information [1, 2] as
universal algorithms, [3, 4] but also in condensed matter [5, 6], as models of
topological phases [7, 8]. Quantum walks can also provide a simulation tool of
quantum systems [9, 10] or be used as an instrument to probe complex systems
[11, 12, 13]. Various experimental implementations demonstrated their
feasibility using ions or photons, including various walkers and interactions
[14, 15, 16, 17, 18].

Topological nontrivial phases of matter are an active field of investigation
since the discovery of the quantum Hall effect [19]. These condensed phases
can emerge in widely different materials, from topological insulators [20, 21]
to chiral magnetic films [22, 23]. A variety of topological effects arise from
the interaction between spin and momentum, as in the spin-orbit coupling: edge
states in topological insulators and skyrmion lattices in magnetic materials
arise from this coupling. It is interesting to note that the quantized
anomalous Hall effect was experimentally exhibited in a three dimensional
topological insulator [24], doped with magnetic impurities. Yet, quantum walks
are based on the same kind of coupling. Indeed, in a quantum walk a particle
at a site moves to the neighboring sites according to its spin state; in the
simplest one dimensional walk the spin up projection moves to the right and
the spin down projection to the left: after one walk step the particle is in a
superposition of two sites and two spin projections. Therefore, it is a
natural idea to use quantum walks from quantum information theory to explore
topological effects from condensed matter physics [25].

In this paper we investigate the motion of two interacting walkers on a one
dimensional lattice [26, 27, 28, 29]. The lattice is partitioned into two
regions with different topology, in such a way that an edge channel appears at
the interface (in one dimension it reduces to a localized state at the origin)
[30, 31]. One question to be asked is how does the interaction modify the edge
state, or equivalently, what is the interplay between interaction and
topology. We will see that the interaction can indeed create a bound state at
the interface, when without interaction there is no edge state, and inversely,
it can destroy an existing edge state. Another question that naturally arises
is whether the existence of localized states modifies the behavior of the
entanglement entropy. We will show that this is indeed the case, the
entanglement growths at different rates depending on the way particles
propagate.

After a detailed presentation of the model, we describe the phenomenology of
the two particle quantum walk, then, we present our results on its
localization and entanglement properties, followed by a brief discussion and a
conclusion.

##  2 Model

The Hilbert space of a walker is defined by two quantum numbers, its position
$x=-L/2,\ldots,L/2$ that can take $L$ integer values, and its spin $s=0,1$ (up
and down, respectively). The Hilbert space of the two particles is the
Kronecker product of the one-particle Hilbert spaces
$\mathcal{H}=\mathcal{H}_{1}\otimes\mathcal{H}_{2}$,

$\ket{x_{1}s_{1}x_{2}s_{2}}=\ket{x_{1}s_{1}}\otimes\ket{x_{2}s_{2}}\in\mathcal{H}\,.$
(1)

The dimension of the two particle space is
$\mathrm{dim}\,\mathcal{H}=(2L)^{2}$.

Figure 1: Topological index (charge) as a function of
$(\theta_{\\!+},\theta_{\\!-})$; dots are for typical numerical parameters
used: $\theta_{\\!-}=\pi/4$ (dashed line),
$\theta_{\\!+}=-\pi/3,\pi/16,3\pi/8$ (for $x>0$, red points), sweeping regions
of charge $c=-1,0,1$, respectively, and $\theta_{\\!+}=-\pi/16,9\pi/16$ (for
$x<0$, blue points), corresponding to the charges $c=0,1$. The solid lines
separates parameter regions with different charge.

The spin state can be modified at each site by the unitary operator,

$R(\theta)=1_{x}\otimes\exp(\mathrm{i}\sigma_{y}\theta)$ (2)

where $1_{x}$ is the unit in the position Hilbert subspace, and the rotation
angle $\theta\in(-\pi,\pi)$, could in principle depend on the site
($\sigma_{y}$ is a Pauli matrix). Each particle moves from site $x$ to a
neighboring site $x\pm 1$, depending on its spin polarization; the operators
$T_{\pm}$ shift the up or down spins to the right or the left, respectively,

$\displaystyle T_{\\!+}$
$\displaystyle=\sum_{x}\left(\ket{x+1}\bra{x}\otimes\ket{0}\bra{0}+\ket{x}\bra{x}\otimes\ket{1}\bra{1}\right)\,,$
(3) $\displaystyle T_{\\!-}$
$\displaystyle=\sum_{x}\left(\ket{x}\bra{x}\otimes\ket{0}\bra{0}+\ket{x-1}\bra{x}\otimes\ket{1}\bra{1}\right)\,,$
(4)

which are also unitary (we take the lattice spacing as the unit length). One
step of the two particle quantum walk is defined by the unitary “evolution”
operator,

$U=V(U_{0}\otimes U_{0})\,,\quad
U_{0}=T_{\\!-}R(\theta_{\\!-})T_{\\!+}R(\theta_{\\!+})$ (5)

where $U_{0}$ is the one-particle operator and $V$ is the interaction
operator,

$V=1+\sum_{xs_{1}s_{2}}(\mathrm{e}^{\mathrm{i}\phi}-1)\ket{xs_{1}xs_{2}}\bra{xs_{1}xs_{2}}\,,$
(6)

where $1$ is the unit matrix in $\mathcal{H}$ (the $-1$ term cancels with the
first term when $x_{1}=x_{2}$). The interaction acts only when the two
particles share the same site, by adding a phase $\phi$ to the corresponding
amplitude. The quantum walk is then characterized by the set of angles
$\theta_{\\!\pm}$ and the interaction phase $\phi$. One step of the walk is,

$\ket{\psi(t+1)}=U\ket{\psi(t)}$ (7)

(we take the time step as the unit of time, and $\hbar=1$) where
$\ket{\psi(t)}$ is the quantum state of the walk at time $t$ (a vector of
dimension $(2L)^{2}$). The dimension of the matrix $U$ is $(2L)^{2\times 2}$.

The split-step quantum walk defined by $U_{0}$ is known to possess nontrivial
topology, depending on the choice of the pair $(\theta_{\\!+},\theta_{\\!-})$
[32]. In Fig. 1 we represented the topological phases of the split-step walk
in the angles plane, and a set of typical parameters used in the numerical
calculations. To study the interplay of interaction and topology in the
quantum walk we introduce an interface separating two regions with different
values of the topological charge $c$. These regions are defined by different
choices of the rotation angles for the left ($\theta_{L}$) and right
($\theta_{R}=\theta$) parts of the interface:

$\theta_{\\!+}(x)=\theta_{L}+(\theta-\theta_{L})\mathrm{H}(x),\quad\theta_{\\!-}=\pi/4$
(8)

where $\theta_{L}$ can be $-\pi/16$ to set the left region charge to $c=0$ or
$9\pi/16$ to set it to $c=1$ (as shown in Fig. 1 by the blue points), and we
denote these two cases with the label ‘c’ ($\mbox{c}=0,1$); for each
$\theta_{L}$ we vary the right angle $\theta=(-\pi/3,\pi/16,3\pi/8)$, to
account for $c=-1,0,1$ (the three black points of Fig. 1), respectively
($\mathrm{H}$ is the Heaviside function). We label ‘i’ ($\mbox{i}=0,1,2$) the
choice of $\theta$ for the three charges $c=-1,0,1$, respectively. We fixed
$\theta_{\\!-}=\pi/4$. Combinations of the labels ‘c’ (left charge) and ‘i’
(right charge) will be used to identify the parameters used in the
calculations.

Initial states are chosen to be Bell states with $x_{1},x_{2}=0$,

$\displaystyle\ket{\phi_{\pm}}$
$\displaystyle=\frac{1}{\sqrt{2}}\big{(}\ket{0000}\pm\ket{0101}\big{)}$ (9)
$\displaystyle\ket{\psi_{\pm}}$
$\displaystyle=\frac{1}{\sqrt{2}}\big{(}\ket{0001}\pm\ket{0100}\big{)}$ (10)

or a separable state $\ket{\psi_{0}}=\ket{0000}$, where the $+$-sign states
are symmetric and $-$-sign states are antisymmetric. We assign a label $b$ to
these states which takes the values ‘0,1,2,3’, for the Bell states, and ‘4’
for the product state $\ket{\psi_{0}}$. Typical values for the interaction
are,

$\phi=(0,\pi/3,\pi/2,3\pi/4,\pi)$

labeled ‘g’ ($\mbox{g}=0,\ldots,4$), respectively. Therefore, the code

$(\mbox{cbgi}),\;\mbox{c}=0,1,\;\mbox{b}=0,\ldots,4,\;\mbox{g}=0,\ldots,4,\;\mbox{i}=0,1,2,$

specifies the set of parameters used in the numerical computations: for
instance $(0321)$ corresponds to (left charge 0, state $\psi_{\\!-}$, strength
$\phi=\pi/2$, right charge $c=0$).

Figure 2: Quasienergies in the non interacting case, $\phi=0$ (black points).
(top) Distribution on the unit circle for
$(\theta_{\\!+},\theta_{\\!-})=(3\pi/7,2\pi/9)$. (bottom) The same
quasienergies in increasing order, showing the band gaps. Red points are the
eigenvalues of the quasimomentum $P$, whose eigenvectors are common to $U$,
for a system size of $L=21$. The gap arround $E=0$ in the single particle
system, splits into two gaps, as a consequence of the band wrapping
characteristic of levels defined modulo $2\pi$.

In summary, the model consists in two particles moving in a one dimensional
lattice with an interface at the origin; their evolution is ensured by a one
time step operator $U$, composed of a coin (defined with different rotation
angles in the two regions), which couples to the shift operator for each
particle, and a local on-site interaction between the two particles. Although
the interaction is spin independent, the global unitary operator cannot be
splitted into a spin dependent part and a position dependent part. Coin and
shift entangle spin and position in the one particle sector, and the
interaction the positions in the two particle sector; as a consequence, the
interaction can modify differently the evolution of a given initial condition
depending on its symmetry: symmetric and antisymmetric Bell states (9). In
other words, $U=V\,(U_{0}\otimes U_{0})$ entangles spin and position ($U_{0}$)
of each particle and the two particle positions ($V$), resulting in
entanglement of all degrees of freedom labeled by the basis vectors
$\ket{x_{1}s_{1}x_{2}s_{2}}$ of the Hilbert space.

##  3 Results

Figure 3: Band structure in the $(p_{n},E_{n})$ plane, for $\phi=\pi/3$ (other
parameters as in Fig. 2). Note the appearance of bound modes within the gaps
in the presence of interaction. These modes correspond to the the two
particles occupying the same location. (The bands are colored for clarity.)

The spectrum of the evolution operator is given by the solution of the
eigensystem,

$U\ket{E}=\mathrm{e}^{-\mathrm{i}E}\ket{E}\,,$ (11)

with $E$ the quasienergies and $\ket{E}$ the eigenvectors. Since the
interaction depends only on the distance between the two particles, the total
quasimomentum $P=k_{1}+k_{2}$, is conserved; here $(k_{1},k_{2})$ are the
common eigenvalues of the translation operators $T_{\pm}$, for the two
particles. Therefore, the energy eigenvectors $\ket{E}$, can be classified
according to the eigenvalues $p_{n}$ of $P$ (see Figs. 2 and 3), which are the
roots of the unity. In Fig. 2 we plot the quasienergies on the unit circle for
the free system; the spectrum for the interacting system is shown in Fig. 3.
We observe, in the interacting case, the appearance of bound states in the
gaps, which correspond to the binding of the two particles (molecular state).
This effect is similar to the Hadamard coin case investigated by Ahlbrecht et
al. Ref. [27], and to a related model of two interacting fermions studied by
Bisio et al. Ref. [29]. It is worth noting that, due to the wrapping of the
energy bands (energies are defined modulo $2\pi$), the position and width of
the gaps of the composite system do not trivially follow the one particle
system, even for vanishing interaction strength. As a consequence, the
parameter range for the existence of edge states (localized states at the
origin $x=0$), change for the two particles quantum walk.

To study how do the edge states depend upon the quantum walk parameters, we
consider the eigenstates $\ket{E}$ of $U$. Let us start with the one particle
case. The edge states should correspond to the $E=0,\pi$ eigenmodes (Dirac
cones). We are interested in the behavior of the wavefunction
$\psi_{E}(x)=\braket{x}{E}$, far form the interface $x=0$. For a localized
mode, we expect an exponentially decreasing function of $x$,

$\psi_{E}(x)\sim\mathrm{e}^{-x/\xi_{E}}$ (12)

with a characteristic scale $\xi_{E}=\xi_{E}(\theta_{+},\theta_{-})$, the
localization length of the eigenmode of energy $E$. To determine the
localization length it is convenient to use the transfer matrix method: [33]

$\xi_{E}^{-1}=\lim_{x\rightarrow\infty}\frac{1}{x}\log\bigg{|}\mathrm{Tr}\,\prod_{n=1}^{x}M_{n}(E)\bigg{|}$
(13)

where $M_{n}$, the transfer matrix, relates the state at points $(x,x-1)$ to
points $(x+1,x)$. The general case, with two particles, is not amenable
analytically, however, we can consider the special case where the two
particles share the same position, and restrict the evolution operator to the
subspace $x_{1}=x_{2}$:

$U\,\rightarrow\,\mathrm{e}^{\mathrm{i}\phi}U_{0}\,,$

leading to an effective one particle split-step walk with a phase change at
each step. In this case the transfer matrix reduces to a two blocks $4\times
4$ matrix $M$ independent of $x$:

$\begin{pmatrix}\psi_{0}(x+1)\\\ \psi_{1}(x+1)\\\ \psi_{0}(x)\\\
\psi_{1}(x)\end{pmatrix}=M\begin{pmatrix}\psi_{0}(x)\\\ \psi_{1}(x)\\\
\psi_{0}(x-1)\\\ \psi_{1}(x-1)\end{pmatrix}$ (14)

that can be written as

$M=\begin{pmatrix}m&0\\\ 0&m\end{pmatrix}$ (15)

where, using $c_{\pm}=\cos\theta_{\pm}$ and $s_{\pm}=\sin\theta_{\pm}$,

$m=\frac{1}{c_{+}c_{-}}\begin{pmatrix}m_{00}&m_{01}\\\
m_{10}&m_{11}\end{pmatrix}\,,$ (16)

with

$\displaystyle m_{00}$
$\displaystyle=\mathrm{e}^{-\mathrm{i}(E+\phi)}+2\mathrm{e}^{\mathrm{i}\phi}s_{+}s_{-}+\mathrm{e}^{\mathrm{i}(E+2\phi)}s_{+}^{2}$
$\displaystyle m_{01}$
$\displaystyle=\mathrm{e}^{\mathrm{i}(E+\phi)}c_{+}s_{+}+c_{+}s_{-}=m_{10}$
$\displaystyle m_{11}$
$\displaystyle=\mathrm{e}^{\mathrm{i}(E+\phi)}c_{+}^{2}\,.$

The eigenvalues $\lambda_{\pm}$ of $m$ give, from Eq. (13), the explicit
expression of the localization length,

$\xi_{E}^{-1}=\log\big{|}\max\left[\lambda_{+}(E),\lambda_{-}(E)\right]\big{|}=\Lambda(E)\,.$
(17)

Because the determinant $\det(M)=\det(m)=1$, we have
$\lambda_{+}\lambda_{-}=1$. For a related computation see Rakovszky and Asboth
work on the localization in the split-step quantum walk Ref. [34]. For $E=0$,
we obtain,

$\lambda_{\pm}(0)=\frac{\mathrm{e}^{-\mathrm{i}\phi}}{2c_{+}c_{-}}\bigg{[}\mathrm{e}^{2\mathrm{i}\phi}+2\mathrm{e}^{\mathrm{i}\phi}s_{+}s_{-}+1\pm\sqrt{\left(\mathrm{e}^{2\mathrm{i}\phi}+2\mathrm{e}^{\mathrm{i}\phi}s_{+}s_{-}+1\right)^{2}-4\mathrm{e}^{2\mathrm{i}\phi}c_{+}^{2}c_{-}^{2}}\bigg{]}$
(18)

and,

$\lambda_{\pm}(\pi)=\frac{\mathrm{e}^{-\mathrm{i}\phi}}{2c_{+}c_{-}}\bigg{[}\mathrm{e}^{2\mathrm{i}\phi}+2\mathrm{e}^{\mathrm{i}\phi}s_{+}s_{-}-1\mp\sqrt{\left(\mathrm{e}^{2\mathrm{i}\phi}+2\mathrm{e}^{\mathrm{i}\phi}s_{+}s_{-}+1\right)^{2}+4\mathrm{e}^{2\mathrm{i}\phi}c_{+}^{2}c_{-}^{2}}\bigg{]}$
(19)

for $E=\pi$. Graphs of (17) for various parameters are shown in Fig. 4.

Figure 4: Inverse of the localization length $\Lambda(E)=\xi_{E}^{-1}$, for
$E=0$ (first column) and $E=\pi$ (second column), and two values of the
interaction phase $\phi=0,\pi/3$. The main effect of the interaction is to
increase the area (white regions), in the parameter space
$(\theta_{+},\theta_{-})$ with diverging $\xi$ (unlocalized states). Grey
scales from white (zero level) to black ($\Lambda(E)=6$).

Figure 5: Two particle probability distribution (after 60 steps), with initial
state Bell ‘0’ and left charge ‘0’ and right charge ‘1’. In the non-
interacting case there is a localized state at the origin (left, code
$(0000)$); in the presence of interaction, $\phi=\pi/2$, a binding state of
the two particle becomes itinerant (right, code $(0020)$); for a weak
interaction strength the localized state persists. In the trivial interacting
system, a localized state appears $(0031)$, for $\phi=3\pi/4)$, while for the
antisymmetric state $(0330)$ is delocalized, in spite of the interface.

This simple computation shows that the interaction changes the localization
properties of the quantum walk, allowing for instance the existence of
delocalized states in regions of the parameter space where, without
interaction, the states are localized, and vice-versa. These interacting,
delocalized states correspond to the binding states, for which the two
particles occupy the same position (molecular states). This is just the effect
found in the calculation of the localization length, predicting the spreading
over the line of the molecular state.

To illustrate this point, we plotted in Fig. 5, the distribution probability
of the two particles, for an initial symmetric Bell state (state ‘0’), and for
two strengths of the interaction. The origin separates a zero charge $(x<0)$
region to a one charge $(x>0)$ region, creating at zero interaction, a
localized state (Fig. 5, first pannel). For finite interaction, a two particle
bound state appears, which may propagate away from the origin. However, the
probability distribution along the propagation line $x_{1}=x_{2}$ lost its
symmetry with respect to the origin (Fig. 5, second pannel). This is a
consequence of the edge state at the origin, which do not disappear with the
interaction in this range of parameters. Another behavior is observed without
interface for $\phi=3\pi/4$ (Fig. 5, third pannel), a localized molecular
state appears together with enhanced repulsion for an initial symmetric state,
more reminiscent of an antisymmetric state. The opposite situation, a free
localized state becoming extended in the presence of interaction, is found for
exemple in the case of an antisymmetric initial state (last pannel of Fig. 5);
propagation in the $x_{1},x_{2}<0$ region, inhibited for $\phi=0$, is
reestablished for $\phi=3\pi/4$.

In order to delimitate the regions in parameter space possessing localized
states, we computed the probability to stay at the origin for one or two
particles:

$\displaystyle P(t)$
$\displaystyle=\sum_{s_{1}x_{2}s_{2}}|\braket{0s_{1}x_{2}s_{2}}{\psi(t)}|^{2}\,,$
(20) $\displaystyle P_{0}(t)$
$\displaystyle=\sum_{s_{1}s_{2}}|\braket{0s_{1}0s_{2}}{\psi(t)}|^{2}\,,$ (21)

respectively. If these probabilities remain finite, decreasing at a rate
slower than $1/N$, after $N$ time steps, we say that localized states exists.
This is a qualitative definition that seems to be precise enough to detect a
range of parameter values corresponding to localization at the origin. In Fig.
6 we fixed two values of the left rotation angle $\theta_{L}$, and two initial
states $\phi_{\\!+}$ and $\psi_{\\!-}$, and varied the right angle $\theta$
and interaction phase $\phi$. We plot $P(t)$ as a function of $(\theta,\phi)$
for the two Bell states, one symmetric and the other one antisymmetric
(columns), and two types of interfaces, with charge zero or one on the left
region (rows). Note the complex pattern of the localized-delocalized regions
with the emergence of localization in areas of the parameter space, and their
intricate boundary, especially in the case where the left region $x<0$
possesses a $c=1$ topological charge.

Figure 6: Probability of one particle to stay at the origin for different
parameters $(\theta,\phi)$, in a logarithmic scale. The level $1/N$ defines
the frontier of the localized states ($N=65$). The interface separates a left
region in a ‘0’ phase (first row) or ‘1’ phase (second row), and a right
region, which is in phase $-1$ for $\theta<\pi/4$, 1 for $\theta>\pi/4$, and 0
in between. The first column corresponds to an initial ‘0’ Bell state
(symmetric), while the second column corresponds to a ‘3’ Bell state
(antisymmetric). One observes (first row) the emergence of a localized region,
induced by the interaction $\phi$, for the ‘00‘ symmetric case that disappears
for the antisymmetric ‘03’ case, for which almost all parameter space
corresponds to delocalized states. A similar qualitative behavior is found for
an interface with charge ‘1’ on the left (second row).

One important characteristic of the quantum behavior of the system is related
to the growth of the entanglement entropy. As a measure of the entanglement we
use the von Neumann entropy of the position of particle one degree of freedom
$S_{x}$, obtained by tracing out the remaining quantum degrees of freedom
(second particle position, and spin state),

$S_{x}(t)=-\mathrm{Tr}\,\rho_{x}(t)\log\rho_{x}(t),\;\rho_{x}(t)=\mathrm{Tr}_{\bar{x}}\rho_{x,\bar{x}}(t)$
(22)

where $\rho=\rho_{x,\bar{x}}$ is the density matrix of the total system, and
the partial trace $\mathrm{Tr}_{\bar{x}}$, is over the Hilbert space labeled
by the set of quantum numbers $\bar{x}=\\{c_{1},x_{2},c_{2}\\}$, the
complementary set of $\\{x=x_{1}\\}$.

Figure 7: Entanglement entropy of particle one position as a function of time
for case $(0011)$, without topological interface, and case $(0010)$, with a
‘01’ interface. In the first case (left pannel), there is a binding state for
the two particles, and the entropy grows following a $\log$ law, with an
exponent $\alpha\approx 0.54$; in contrast, for the second case (right
pannel), a localized state is present at the origin, and the entropy growth is
only $\log\log$. The insets show the unscaled graph.

To analyse the entropy growth, we assume that for enough long times, there
exists a probability distribution $p_{S}=p_{S}(x,t)$ such that the von Neumann
entropy can be approached by the expression,

$S_{x}(t)\approx-\int_{-\infty}^{\infty}\mathrm{d}x\,p_{S}(x,t)\log
p_{S}(x,t)\,,$ (23)

with the normalization,

$\int_{-\infty}^{\infty}\mathrm{d}x\,p_{S}(x,t)=1\,.$ (24)

In the delocalized state, we consider that the growth of the entanglement can
be described by a self-similar form:

$p_{S}(x,t)=\frac{1}{t^{\alpha}}P_{S}\left(\frac{x}{t^{\alpha}}\right)\,,$
(25)

with the self-similar variable $X=x/t^{\alpha}$, characterized by a single
exponent $\alpha$. This form guarantees that the normalization condition is
automatically satisfied if $P_{S}$ is itself normalized. Substituting into
(23), and using (24), we obtain

$S_{x}(t)=S_{0}+\alpha\log t\,,$ (26)

where

$S_{0}=-\int_{-\infty}^{\infty}\mathrm{d}X\,P_{S}(X)\log P_{S}(X)\,.$

The exponent depends in general on the interaction strength
$\alpha=\alpha(\phi)$, vanishing for $\phi=0$. The self-similar form (25) is
borrowed from the one particle quantum walk, for which the ballistic law
$x$\mathrm{t}$$ applies [35, 2]; this is certainly appropriated when the
probability distribution is dominated by the motion of the bounded state of
the two particles, and then we have an effective one particle walk, or in the
other limit, when there is repulsion (as in the antisymmetric case) and the
two particles are almost independent. The fact that the characteristic
exponent varies with the interaction can be related to the leak of one
particle position probability in correlations with the other degrees of
freedom: the effective motion of one particle is no more ballistic.

In the localized state, for times of the order of $t(x)\sim\mathrm{e}^{x/\xi}$
for $x\gg\xi$ (where $\xi$ is the localization length), the entanglement
probability $p_{S}(x,t)\rightarrow P_{S}(x/\log t)$ tends towards an almost
stationary distribution. This leads to a very slow, double logarithmic growing
law,

$S_{x}(t)\sim\log(\xi\log t)\,,$ (27)

difficult to test numerically.

In Fig. 7 we represented the entanglement entropy for a delocalized two
particle bounded state (left, labeled ‘0011’) and for a localized state
(right, labeled ‘0010’), for the initial Bell state ‘0’, illustrating the two
growing modes.

##  4 Discussion and conclusions

We investigated the localization and entanglement properties of a system of
two interacting particles executing a topological quantum walk. The behavior
of a two particle quantum walk cannot be straightforwardly deduced from the
one particle case, even without interaction [26, 18], new effects arise
related to the symmetry of the initial state (which is preserved by the
interaction) and the wrapping of the energy bands.

The presence of edge states at the interface of two regions with different
topology, and the appearance of interacting molecular states, further modifies
the properties of the system. In fact, the evolution of the system depends on
the type of the interface; an interface between charge $-1$ and $0$ regions,
or $1$ and $0$ are not equivalent. Furthermore, the overlapping of the initial
state with the edge state (for example, putting the particles on the left or
on the right of the interface), will affect the details of the system’s
evolution. However, the basic properties, like the localization and binding
effects, are robust, they do not depend on the details of the initial
condition.

More specifically, we showed that new localization-delocalization transitions
were possible, depending on the interaction phase and the symmetry of the
quantum states. Antisymmetric states favor delocalization of the walk, and the
unbinding of the molecule. Symmetric states, on the contrary, preserve the
edge states and the molecular binding. For the particular value $\phi=\pi$,
strongly localized states are observed. For intermediate values, in the range
$\phi=(\pi/3,\pi/2)$, localized states in the free delocalized region appear,
for symmetric states, but also for antisymmetric states depending on the
topology.

Therefore, the effect of the interaction is that it can induce localized
states, absent for free particles, in the trivial topological case. This
effect is present, for different values of the interaction strength,
independently of the initial state. The molecular states, which can propagate
along the line, compete with the localized state, which tends to stick the
particles at the interface. Moreover, depending on the phase of the
interaction, a delocalization is possible in the presence of the interface.
Indeed, the interaction can also destroy the edge states to allow the
propagation of the information between regions with different topology.
Therefore, the interaction can be used as a control parameter of the
information transport.

Interaction not only contributes to modify the localization properties of the
quantum walk, it is the physical mechanism of entanglement generation. We
measured the entanglement entropy of one particle and found that it steadily
increases over time [12, 36]. However, at variance with other interacting
systems, its growth in not simply proportional to time, but shows two
different regimes according to the localization of the walk. For delocalized
states the entropy growth logarithmically, while for localized states it
growth double logarithmically. We showed some numerical evidence of these
different growth laws.

In conclusion, using a quantum walk we studied various basic mechanisms, such
as the interplay of topology and interaction, or the relation between
entanglement and interaction. These effects also have interesting analogies
with condensed matter systems. Indeed, quantum walks may be used to simulate
material systems [37], or can be related to specific condensed matter systems
via their effective hamiltonian. Simple quantum walks defined by a set of
local unitary operations, can be equivalent to complex nonlocal effective
quantum systems. We observed for instance, the growing of the entanglement
entropy associated with long range correlations induced by a local interaction
rule; or the appearance of localized states without disorder, induced by the
interaction in a background trivial topology.

## Acknowledgements

We thanks Simon Lucchesini who wrote a first version of the numerical code. We
benefited from discussions with Giuseppe Di Molfetta, Laurent Raymond and Luis
Foa Torres. Part of this work was funded by the Université de Toulon.

## References

  * [1] J. Kempe,  _Quantum random walks: An introductory overview_ ,  Contemporary Physics 44(4), 307 (2003),  10.1080/00107151031000110776.
  * [2] S. E. Venegas-Andraca,  _Quantum walks: a comprehensive review_ ,  Quantum Information Processing 11(5), 1015 (2012),  10.1007/s11128-012-0432-5.
  * [3] A. M. Childs,  _Universal computation by quantum walk_ ,  Phys. Rev. Lett. 102, 180501 (2009),  10.1103/PhysRevLett.102.180501.
  * [4] N. B. Lovett, S. Cooper, M. Everitt, M. Trevers and V. Kendon,  _Universal quantum computation using the discrete-time quantum walk_ ,  Phys. Rev. A 81, 042330 (2010),  10.1103/PhysRevA.81.042330.
  * [5] T. Kitagawa,  _Topological phenomena in quantum walks: elementary introduction to the physics of topological phases_ ,  Quantum Information Processing 11(5), 1107 (2012),  10.1007/s11128-012-0425-4.
  * [6] B. Zeng, X. Chen, D.-L. Zhou and X.-G. Wen,  _Quantum Information Meets Quantum Matter – From Quantum Entanglement to Topological Phase in Many-Body Systems_ ,  ArXiv e-prints (2015),  1508.02595.
  * [7] T. Kitagawa, M. S. Rudner, E. Berg and E. Demler,  _Exploring topological phases with quantum walks_ ,  Phys. Rev. A 82, 033429 (2010),  10.1103/PhysRevA.82.033429.
  * [8] J. K. Asbóth,  _Symmetries, topological phases, and bound states in the one-dimensional quantum walk_ ,  Phys. Rev. B 86, 195414 (2012),  10.1103/PhysRevB.86.195414.
  * [9] G. Di Molfetta and A. Pérez,  _Quantum walks as simulators of neutrino oscillations in a vacuum and matter_ ,  New J. Phys. 18(10), 103038 (2016),  10.1088/1367-2630/18/10/103038.
  * [10] P. Arrighi, S. Facchini and M. Forets,  _Quantum walking in curved spacetime_ ,  Quantum Information Processing pp. 1–20 (2016),  10.1007/s11128-016-1335-7.
  * [11] B. L. Douglas and J. B. Wang,  _A classical approach to the graph isomorphism problem using quantum walks_ ,  J. Phys. A: Math. Gen. 41(7), 075303 (2008),  10.1088/1751-8113/41/7/075303.
  * [12] S. D. Berry and J. B. Wang,  _Two-particle quantum walks: Entanglement and graph isomorphism testing_ ,  Phys. Rev. A 83, 042317 (2011),  10.1103/PhysRevA.83.042317.
  * [13] G. D. Paparo and M. A. Martin-Delgado,  _Google in a quantum network_ ,  Sci. Rep. 2, 444 (2012),  10.1038/srep00444.
  * [14] F. Zähringer, G. Kirchmair, R. Gerritsma, E. Solano, R. Blatt and C. F. Roos,  _Realization of a quantum walk with one and two trapped ions_ ,  Phys. Rev. Lett. 104, 100503 (2010),  10.1103/PhysRevLett.104.100503.
  * [15] A. Peruzzo, M. Lobino, J. C. F. Matthews, N. Matsuda, A. Politi, K. Poulios, X.-Q. Zhou, Y. Lahini, N. Ismail, K. Wörhoff, Y. Bromberg, Y. Silberberg _et al._ ,  _Quantum walks of correlated photons_ ,  Science 329(5998), 1500 (2010),  10.1126/science.1193515.
  * [16] A. Schreiber, K. N. Cassemiro, V. Potocek, A. Gabris, I. Jex and C. Silberhorn,  _Decoherence and disorder in quantum walks: From ballistic spread to localization_ ,  Phys. Rev. Lett. 106, 180403 (2011),  10.1103/PhysRevLett.106.180403.
  * [17] A. Schreiber, A. Gábris, P. P. Rohde, K. Laiho, M. Štefaňák, V. Potoček, C. Hamilton, I. Jex and C. Silberhorn,  _A 2d quantum walk simulation of two-particle dynamics_ ,  Science 336(6077), 55 (2012),  10.1126/science.1218448.
  * [18] A. Crespi, R. Osellame, R. Ramponi, V. Giovannetti, R. Fazio, L. Sansoni, F. De Nicola, F. Sciarrino and P. Mataloni,  _Anderson localization of entangled photons in an integrated quantum walk_ ,  Nature Photon. 7(4), 322 (2013),  10.1038/nphoton.2013.26.
  * [19] R. Prange, S. Girvin, A. Chang, F. Duncan, R. Laughlin, A. Pruisken and D. Thouless,  _The Quantum Hall Effect_ ,  Graduate Texts in Contemporary Physics. Springer New York, 2nd ed. edn.,  10.1007/978-1-4612-3350-3 (1990).
  * [20] B. A. Bernevig, T. L. Hughes and S.-C. Zhang,  _Quantum spin Hall effect and topological phase transition in HgTe quantum wells_ ,  Science 314, 1757 (2006),  10.1126/science.1133734.
  * [21] X.-L. Qi and S.-C. Zhang,  _Topological insulators and superconductors_ ,  Rev. Mod. Phys. 83, 1057 (2011),  10.1103/RevModPhys.83.1057.
  * [22] S. Muhlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii and P. Boni,  _Skyrmion Lattice in a Chiral Magnet_ ,  Science 323(5916), 915 (2009),  10.1126/science.1166767.
  * [23] W. Jiang, G. Chen, K. Liu, J. Zang, S. G. E. te Velthuis and A. Hoffmann,  _Skyrmions in magnetic multilayers_ ,  Phys. Rep. 704(Supplement C), 1 (2017),  https://doi.org/10.1016/j.physrep.2017.08.001.
  * [24] C.-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L.-L. Wang, Z.-Q. Ji, Y. Feng _et al._ ,  _Experimental observation of the quantum anomalous hall effect in a magnetic topological insulator_ ,  Science 340(6129), 167 (2013),  10.1126/science.1234414.
  * [25] A. D. Verga,  _Edge states in a two-dimensional quantum walk with disorder_ ,  Eur. Phys. J. B 90(3), 41 (2017),  10.1140/epjb/e2017-70433-1.
  * [26] Y. Omar, N. Paunković, L. Sheridan and S. Bose,  _Quantum walk on a line with two entangled particles_ ,  Phys. Rev. A 74, 042304 (2006),  10.1103/PhysRevA.74.042304.
  * [27] A. Ahlbrecht, A. Alberti, D. Meschede, V. B. Scholz, A. H. Werner and R. F. Werner,  _Molecular binding in interacting quantum walks_ ,  New J. Phys. 14(7), 073050 (2012),  https://doi.org/10.1088/1367-2630/14/7/073050.
  * [28] J. Lockhart and M. Paternostro,  _Discrete-time quantum walks as generators of multipartite entanglement_ ,  ArXiv e-prints (2015),  1509.05816.
  * [29] A. Bisio, G. M. D’Ariano, P. Perinotti and A. Tosini,  _Thirring quantum cellular automaton_ ,  Phys. Rev. A 97, 032132 (2018),  10.1103/PhysRevA.97.032132.
  * [30] H. Obuse and N. Kawakami,  _Topological phases and delocalization of quantum walks in random environments_ ,  Phys. Rev. B 84, 195139 (2011),  10.1103/PhysRevB.84.195139.
  * [31] T. Kitagawa, M. A. Broome, A. Fedrizzi, M. S. Rudner, E. Berg, I. Kassal, A. Aspuru-Guzik, E. Demler and A. G. White,  _Observation of topologically protected bound states in photonic quantum walks_ ,  Nature Commun. 3, 882 (2012),  10.1038/ncomms1872.
  * [32] C. Cedzich, F. A. Grünbaum, C. Stahl, L. Velázquez, A. H. Werner and R. F. Werner,  _Bulk-edge correspondence of one-dimensional quantum walks_ ,  J. Phys. A: Math. Gen. 49(21), 21LT01 (2016),  10.1088/1751-8113/49/21/21LT01.
  * [33] A. Crisanti, G. Paladin and A. Vulpiani,  _Products of Random Matrices in Statistical Physics_ , vol. 104 of _Springer Series in Solid-State Sciences_ ,  Springer-Verlag, Berlin, Heidelberg,  10.1007/978-3-642-84942-8 (1993).
  * [34] T. Rakovszky and J. K. Asboth,  _Localization, delocalization, and topological phase transitions in the one-dimensional split-step quantum walk_ ,  Phys. Rev. A 92, 052311 (2015),  10.1103/PhysRevA.92.052311.
  * [35] A. Nayak and A. Vishwanath,  _Quantum Walk on the Line_ ,  eprint arXiv:quant-ph/0010117 (2000),  quant-ph/0010117.
  * [36] A. Nahum, J. Ruhman, S. Vijay and J. Haah,  _Quantum entanglement growth under random unitary dynamics_ ,  Phys. Rev. X 7, 031016 (2017),  10.1103/PhysRevX.7.031016.
  * [37] A. Aspuru-Guzik and P. Walther,  _Photonic quantum simulators_ ,  Nature Phys. 8, 285 (2012),  10.1038/nphys2253.


# The photothermal nonlinearity in plasmon-assisted photocatalysis

Ieng Wai Un,a∗, Yonatan Dubib and Yonatan Sivana
aSchool of Electrical and Computer Engineering, Ben-Gurion University of the
Negev, Israel
bDepartment of Chemistry, Ben-Gurion University of the Negev, Beer-Sheva,
8410501, Israel.
∗E-mail<EMAIL_ADDRESS>
###### Abstract

Understanding the intricate relation between illumination and temperature in
metallic nano-particles is crucial for elucidating the role of illumination in
various physical processes which rely on plasmonic enhancement but are also
sensitive to temperature. Recent studies have shown that the temperature rise
in optically-thick ensembles of metal nanoparticles under intense illumination
is dominated by the thermal conductivity of the host, rather than by the
optical properties of the metal or the host. Here, we show that the
temperature dependence of the thermal conductivity of the host dominates the
nonlinear photothermal response of these systems. In particular, this
dependence typically causes the temperature rise to become strongly sublinear,
reaching even several tens of percent. We then show that this can explain
experimental observations in several recent plasmon-assisted photocatalysis
experiments. This shows that any claim for dominance of non-thermal electrons
in plasmon-assisted photocatalysis must account first for this photothermal
nonlinear mechanism.

##  1 Introduction

The use of illuminated metallic surfaces to enhance the yield of chemical
reactions (aka plasmon-assisted photocatalysis) was proposed by Nitzan and
Brus already in 1981 [1] and implemented experimentally shortly later [2].
However, only after several decades of slow progress (see e.g., Refs. 3, 4, 5,
6, 7, 8), this line of research has rapidly gained popularity following
several high impact publications (see, e.g. Refs. 9, 10, 11, 12 for some
recent reviews). The growing interest was propelled by claims in some of the
more famous papers on the topic [13, 14, 15, 16, 17] that the reaction rate
increases due to the excitation of high energy non-thermal electrons in the
metal (aka “hot” electrons), which then tunnel out of the metal, and provide
the necessary energy for the reactants to allow them to be converted into the
products more efficiently.

However, these very papers (Refs. 13, 14, 15, 16, 17) were shown to suffer
from technical and conceptual flaws (including improper temperature
measurements, improper data normalization etc., see discussion in Refs. 18,
19, 20, 21, 22). Instead, a purely thermal mechanism was shown to be able to
explain the experimental data quite convincingly [18, 20, 23, 22]. In
particular, a shifted Arrhenius Law for the reaction rate,
$R\sim\exp\left(-\frac{\mathcal{E}_{a}}{k_{B}T({\bf
r})+aI_{\textrm{inc}}}\right)$ whereby the temperature of the system was
corrected for the illumination-induced heating was shown to provide an
excellent fit to the published data, essentially with no fit parameters. This
result was corroborated with the first ever complete calculation of the
steady-state electron non-equilibrium in metals [24], a consequent Fermi
golden-rule argument [25] that pointed to the improbability of nonthermal
electrons to cause the catalysis, and by detailed thermal simulations where
the dynamics of the heat generated from each of the nanoparticles (NPs) in the
system was properly modelled [20]. Similar criticism was raised in an
independent study [21], and characterization which is inline with our approach
was employed by several groups (see, e.g., Refs. 26, 27, 28, 29, 30).

In a consequent paper [31], it was shown that in many typical configurations,
the tedious modelling of the contributions of each of the heated NPs in the
sample can be replaced by an effective medium approximation. This approach
also enabled accounting for the exact reactor geometry, constituent materials
and boundary conditions, thus, enabling a quantitative comparison to the
measured data (see Fig. 1(a)). This series of works was lately extended to
account also for fluid dynamics effects and for redox reactions [32] (Fig.
1(b)).

The bottom line of the thermal modelling was that when attempting to
quantitatively separate thermal and non-thermal effects in plasmon-assisted
photocatalysis experiments, one has to overcome a conceptual difficulty - the
thermocatalysis control experiments must reproduce the exact spatially non-
uniform temperature profile induced by the illumination, otherwise, when
subtracting the thermocatalysis rate from the photocatalysis rate (e.g., as in
Refs. 17, 33), any difference between the temperature distributions in an
inaccurate control and the corresponding photocatalysis experiment is bound to
be incorrectly interpreted as “hot” electron action. This issue is
particularly important because the Arrhenius Law shows that the reaction rate
has an exponential sensitivity to the temperature distribution [20].

Detailed measurements and/or calculations of the temperature distribution in
the studied samples indeed constituted a central role in several recent
demonstrations of non-thermal effects in plasmon-assisted photocatalysis [26,
33, 34, 27, 28]. However, while the simple thermal calculations done so far
were sufficient for relatively simple scenarios, they may not be sufficient to
account for more complicated ones. Those include, in particular, high
intensity illumination which invokes steady-state nonlinear thermo-optic and
photothermal effects [35, 36, 37, 38, 39], which are usually simply ignored
(without any justification). Already in the context of thin metal layers [40,
41, 42, 43] and single NPs [36, 44, 45, 37, 38, 39], this effect was shown to
cause deviations of several tens to hundreds of percent in the permittivity,
and hence in the field and temperature distributions compared with the purely
uniform linear thermal response. Note that such a nonlinear effect is far
greater than conventional nonlinear optical effects.

In this work, we go beyond the study of the linear response, and evaluate the
importance of nonlinear photothermal effects in large random ensembles of
metal NPs, suitable to plasmon-assisted photocatalysis experiments (Fig.
1(a)-(b), but also to many other types of experiments in nonlinear optics (see
e.g., Refs. 46, 47, 48, 49, 35, 50). First, in Section 2 we provide a
qualitative analysis that points to the most important parameter that affects
the overall nonlinear photothermal response of typical plasmon-assisted
photocatalysis systems, namely, the thermal conductivity of the host. We also
identify those parameters which have a negligible effect on the nonlinearity,
namely, the optical parameters. Then, in Section 3 we switch to a rigorous
analysis and describe the methodology we employ to calculate the temperature
distribution in the samples considered. The qualitative analysis is then
applied in Section 4 to two sets of experimental data taken from recent high-
impact plasmon-assisted photocatalysis experiments. The good match of our
analysis with the experimental data indicates that the photothermal
nonlinearity is indeed a significant effect, such that a neglect to account
for it is bound to lead to an overestimate of the role of non-thermal
electrons. Section 5 provides a discussion and outlook.

##  2 A qualitative analysis

In order to achieve a qualitative understanding of the high temperature and/or
intensity response of plasmon-assisted photocatalysis systems (i.e., the
photothermal nonlinearity), we start by considering a simplified
configuration, namely, we assume that the sample consists of metal NPs (with
dielectric permittivity
$\varepsilon_{m}=\varepsilon_{m}^{\prime}+i\varepsilon_{m}^{\prime\prime}$,
thermal conductivity $\kappa_{m}$ and NP number density $n_{p}$) distributed
in a cylinder-shape volume and immersed in a uniform host material (with
thermal conductivity $\kappa_{h}$).

In our previous work [31], we have shown that in the weak illumination limit,
the temperature rise at the top center of such a sample can be approximately
written as

$\displaystyle\Delta
T^{\textrm{top}}\approx\dfrac{I_{\textrm{inc}}\rho_{b}}{2\kappa_{h}}\left(1-e^{-H/\delta_{\textrm{skin}}}\right),$
(1)

where $\rho_{b}$ is the beam radius, $H$ is the sample thickness, and
$\delta_{\textrm{skin}}$ is the penetration (skin) depth of light in the
catalyst sample. The inverse of the penetration (skin) depth (i.e., the
absorption coefficient) is related to the NP number density $n_{p}$ and the
absorption cross-section $\sigma_{\textrm{abs}}$ via

$\displaystyle
1/\delta_{\textrm{skin}}(\omega)=n_{p}\sigma_{\textrm{abs}}(\omega).$ (2)

Figure 1: (Color online) (a) A schematic illustration of a photocatalytic
chamber (left) and of a pellet-based catalyst sample (right). Note the
sparsity of the metal nanoparticles and the near-percolation nature of the
micron-size oxide particles. (b) A schematic illustration of a typical
plasmon-assisted photocatalysis setup based on metal nanoparticles liquid
suspension. In both set-ups, a heat flux boundary condition is applied to all
domain boundaries (marked by dashed lines). (c) The heat source density
(orange solid line) and temperature (orange dashed line) along the optical
axis into the catalyst sample (labeled by the orange solid line in (a)). The
black dotted lines represent the top ($T_{1}$) and the bottom ($T_{2}$)
surface temperature.

We now recall that in plasmon-assisted photocatalysis experiments, the light
penetration depth is usually designed to be smaller than the sample thickness
(i.e., such that $H\ll\delta_{\textrm{skin}}(\omega)$ for all wavelengths in
the illumination) to ensure that all the illumination energy is absorbed, see
Fig. 1(c). Eq. (1) shows that under such conditions, the temperature rise
simply becomes $\Delta
T^{\textrm{top}}\approx\dfrac{I_{\textrm{inc}}\rho_{b}}{2\kappa_{h}}$, so that
the overall temperature rise is weakly-sensitive to the illumination spectrum,
NP shape, size, and density, but exhibits an inverse proportion to the host
thermal conductivity [31].

When the illumination intensity is increased, the illumination-induced heating
of the NPs causes a modification of the optical and thermal properties of the
NPs and their surrounding. This effect gives rise to a nonlinear dependence of
the sample temperature on the illumination intensity. Eq. (1) and its
approximation shows that the photothermal nonlinearity would be dominated by
the temperature dependence of the thermal conductivity of the host,
$\kappa_{h}$. To understand this potentially non-intuitive result, we note
that the increase of the imaginary part of the metal permittivity
($\varepsilon_{m}^{\prime\prime}$) with temperature reduces the quality factor
of the plasmonic resonance of the NPs [37, 38, 39]. Meanwhile, the change of
the real part of the metal permittivity ($\varepsilon_{m}^{\prime}$) causes a
resonance shift of the absorption spectrum. Although the impact of the change
of the metal permittivity on the absorption cross-section varies in a complex
manner with the NP size and the illumination wavelength, see e.g. Refs. 37,
38, 39, and although the sensitivity of the metal permittivity to the rising
temperature is relatively high [51, 52, 48, 41, 42, 37, 43], the metal
permittivity has a relatively small contribution to the overall photothermal
response of the sample when the penetration depth is much thinner than the
sample thickness, as shown in Eq. (1) because the sample absorbs all light
regardless of these changes. Moreover, since the temperature dependence of the
optical properties of the host has a similar effect on the absorption cross-
section as the real part of the metal permittivity [38, 39] (i.e., it causes a
resonance shift), the nonlinear response due to the change of the host
permittivity is also small.

On the other hand, the contribution of the thermal properties to the overall
photothermal nonlinearity is significant; it naturally depends on the volume
fraction of the various materials. Because the metal occupies a small fraction
of the sample volume, one can appreciate that the change of the metal thermal
properties hardly contributes to the nonlinear response of the sample. In
contrast, Eq. (1) shows that the thermal properties of the host matters much
more (be it gas, liquid, a porous composite or even a (dielectric) solid).
Clearly, an illumination-induced increase of the thermal conductivity with the
temperature means that the overall temperature rise in the sample would become
sublinear as a function of the illumination intensity. Judging by the typical
thermoderivative of these properties, the nonlinearity is expected to manifest
itself at a temperature rise of several hundreds of degrees, see Refs. 38, 39;
this estimate is found below to be in good agreement with the experimental
data.

##  3 Rigorous analysis

Now, having understood the expected qualitative behavior of the photothermal
nonlinearity in plasmon-assisted photocatalysis, we turn to describe the
methodology employed to study it rigorously. In order to properly account for
the non-trivial reactor geometry and the multitude of materials, we used a
numerical software package (COMSOL Multiphysics) to solve the time-independent
heat (or Poisson) equation with temperature-dependent parameters.

Two generic configurations were used in plasmon-assisted photocatalysis
experiments. One is based on pellet geometries (Fig. 1(a)) and the other on
NPs in a liquid suspension (Fig. 1(b)). In both cases, we distinguish between
the thermal conductivities of the catalyst sample and of its surrounding to
account for the inhomogeneity of the thermal properties, namely, we solve

$\displaystyle\begin{cases}\nabla\cdot\left[\kappa_{\textrm{cata}}(T({\bf
r}))\nabla T({\bf r})\right]=-p_{\textrm{abs}}({\bf r}),&\textrm{inside the
catalyst,}\\\ \nabla\cdot\left[\kappa_{\textrm{holder}}(T({\bf r}))\nabla
T({\bf r})\right]=0,&\textrm{in the sample holder},\\\
\nabla\cdot\left[\kappa_{\textrm{gas}}(T({\bf r}))\nabla T({\bf
r})\right]=0,&\textrm{elsewhere.}\end{cases}$ (3)

Here, $\kappa_{\textrm{cata}}(T({\bf r}))$, $\kappa_{\textrm{holder}}(T({\bf
r}))$ and $\kappa_{\textrm{gas}}(T({\bf r}))$ are the temperature-dependent
thermal conductivities of the catalyst [20, 31], sample holder, and
surrounding gas, respectively. Finally, $p_{\textrm{abs}}$ is the heat source
density induced by the light absorption in the random metal NP array; it is
well described by the effective medium approximation for the electromagnetic
properties of the catalyst sample [31]. When the skin (penetration) depth is
much smaller than the sample thickness one can neglect the temperature
variation within the skin depth, namely,

$\displaystyle p_{\textrm{abs}}({\bf
r})=\int\dfrac{i_{\textrm{inc}}(\rho,\omega)}{\delta_{\textrm{skin}}(\omega,T_{1})}\exp\left(-\dfrac{z}{\delta_{\textrm{skin}}(\omega,T_{1})}\right)d\omega.$
(4)

Here, $\rho$ is the distance from the propagation optical axis, $z$ is the
distance along the propagation direction of the incident beam from the top
surface of the catalyst sample, and $i_{\textrm{inc}}(\rho,\omega)$ describes
the transverse spatial and spectral profile of the incident beam. $T_{1}$
represents the temperature of the top layer of the catalyst sample and is now
an unknown variable which needs to be determined by solving Eq. (3) and Eq.
(4) self-consistently.

In practice, in Eq. (4) (specifically, when calculating the absorption cross-
section), $T_{1}$ is chosen to be 300 K since the only available data for the
metals used as catalysts in the papers we analyze below (Ru and Cu) are at 300
K. In fact, this is a very good approximation for Eq. (4) because, as pointed
out in Section 2, the temperature dependence of the permittivities has a minor
effect on the overall photothermal response when the skin (penetration) depth
is much smaller than the sample thickness. In addition, this approximation
also allows significant time and computational resource-saving. In contrast,
when the skin depth is comparable to the sample thickness, one needs to
account for the temperature dependence of the permittivities and the
temperature variation within the skin depth (see Supplemental Information
Section LABEL:suppsec:T_dep_derivation for a complete derivation). In this
case, the expression for the heat source density becomes much more complicated
than Eq. (4) (see Eq. (LABEL:suppeq:pabs_closeform)) and one needs to solve
Eq. (3) with Eq. (LABEL:suppeq:pabs_closeform) self-consistently. The earlier
analysis of the temperature dependence of the permittivities [37, 38, 39]
imply that even in this case one should expect a sublinear growth of the
temperature.

In order to limit the simulation domain to a realistic and significant volume,
we set a convection heat flux boundary condition at the outer surfaces of the
simulation domain; this models the heat transfer driven by the temperature
difference between the simulation domain and the distant surrounding (assumed
to be at 20∘C). The associated heat transfer coefficient (denoted by $h$ in
COMSOL Multiphysics) is used as an adjustable parameter to fit the
experimental results in the small intensity limit for which one can neglect
the temperature dependence of the thermal conductivities. In this case, Eq.
(3) becomes a linear Poisson equation, the solution of which is the linear
approximation of the solution of Eq. (3). The obtained value for $h$ has a
negligible effect on the photothermal nonlinearity (see Fig.
LABEL:suppfig:JieLiu_NanoLett_2019_temp_map_T1_T2_h_para). We emphasize that
the frequently employed fixed temperature boundary condition does not
represent correctly the heat transfer from the reactor to the environment.
Instead, it effectively brings the distant cool regions closer to the reactor,
and thus, cools the reactor in a manner that is inconsistent with measured
experimental data.

For simplicity, we apply the non-convection approximation. This was justified
by simplistic estimates [20] which showed that the standard gas flow level
used is not expected to be significant in removing the generated heat, and
further supported by simulations of natural and forced convection done in Ref.
32; this effect anyhow is not expected to affect the nonlinear response.

##  4 Analysis of Experiments

In the following, we apply the model described above to two representative
experiments. Both involved pellet geometries which are typically composed of a
large number ($\sim 10^{12}-10^{14}$) of few nm metal NPs randomly distributed
within a highly sparse 3D powder of micron-size metal oxide particles, see
Fig. 1(a); gases occupy the volume between the NPs. In a typical reaction
chamber, the catalyst sample (typically a few mm in size) is placed on a
(stainless steel) sample holder and the reaction rate is then measured under a
specific illumination and/or under resistive heating.

The heat conduction in such catalyst sample might be expected to be dominated
by the solids (the sample holder and the oxide) since their thermal
conductivities are much larger than that of the gases. However, this is not
the case and instead, as shown previously [31], it is the thermal conductivity
of the gases that dominates. To see that, we first note that since for all
practically useful samples and in all the cases we study below, the heat
generation occurs primarily on the top layer of the catalyst sample, i.e.,
away from most of the sample holder, the sample holder is effective in
reducing the bottom-surface temperature but is less effective in reducing the
top-surface temperature (see Supplemental Information Section
LABEL:suppsec:thermal_simulations for details).

Second, in such samples, the oxide fill factor typically only reaches $\sim
10$% (see Eq. (LABEL:suppeq:volume_fraction)). Indeed, the sample consists of
a highly sparse random array of micron size oxide particles, which are barely
touching each other, see Fig. 1(a); thus, somewhat unintuitivelythe gas serves
as the bottle-neck for the heat conduction in the catalyst sample and the heat
conduction through the oxide is highly inefficient. This conceptual picture
has been verified by a hierarchy of effective medium approaches which was
found to be in good agreement with experimental results (see Ref. 53).
Specifically, that work showed that the Maxwell-Garnett model [54, 55]
provides an excellent approximation for the effective thermal conductivity
$\kappa^{\textrm{eff}}_{\textrm{cata}}$ of the catalyst sample we analyze
below, namely (see Supplemental Information Section
LABEL:suppsec:thermal_simulations for details),

$\displaystyle\kappa^{\textrm{eff}}_{\textrm{cata}}=\kappa_{\textrm{gas}}+\dfrac{3f_{\textrm{oxide}}\kappa_{\textrm{gas}}}{\dfrac{\kappa_{\textrm{oxide}}+2\kappa_{\textrm{gas}}}{\kappa_{\textrm{oxide}}-\kappa_{\textrm{gas}}}-f_{\textrm{oxide}}},$
(5)

where $\kappa_{\textrm{oxide}}$ is the thermal conductivity of the oxide.
Thus, since the fill factor of the oxide is small ($f_{\textrm{oxide}}\ll 1$),
the effective thermal conductivity of the catalyst sample becomes
approximately
$\kappa^{\textrm{eff}}_{\textrm{cata}}\approx(1+3f_{\textrm{oxide}})\kappa_{\textrm{gas}}$
(see Eqs. (LABEL:suppeq:kappa_cata)-(LABEL:suppeq:kappa_cata_approx)), i.e.,
it is very close to that of the gas mixture (see Fig.
LABEL:suppfig:kappa_cata_pellet). Due to the above, in the examples below we
set $\kappa_{cata}=\kappa_{cata}^{eff}$.

###  4.1 Analysis of experiments from Li et al., Nano Letters 2019, 19,
1706-1711 [56]

First, we look at the experimental results of Li et al. [56] who studied
ammonia synthesis using a cesium-promoted, magnesium-oxide supported,
ruthenium (Ru-Cs/MgO) catalyst. The Ru NPs in this study were estimated to be
$\sim 2$nm in diameter. The catalyst sample (3mm height and 6mm diameter) was
put in a reaction chamber equipped with a quartz window (which allows the
catalyst sample to be illuminated at varied intensity and wavelength) and a
temperature controller (which is used to heat up the catalyst sample). The
illumination spot was set to have the same size as the catalyst surface area.
Two thin thermocouples were inserted into the catalyst sample, one to measure
the top-surface temperature (denoted by $T_{1}$) and the other to measure the
bottom-surface temperature (denoted by $T_{2}$). A mixture of N2, H2, and Ar
with a ratio of 1:3 for N2/H2 flowed into the reactor at a total flow rate of
75 sccm. The gaseous product (NH3) was monitored by an online mass
spectrometer. For these parameters, an estimate similar to the one done in
Ref. 20 shows that the gas flow may affect the temperature by no more than
10%, justifying the use of the non-convection approximation in the following
analysis.

Figure 2: (Color online) (a) Temperature dependence of the thermal
conductivities of N2 [57] (blue dotted line), H2 [58] (orange dashed line),
Argon [59] (green dash-dotted line) and the mixture used in Ref. 56 (red solid
line). (b) Light penetration (skin) depth and the absorption cross-section of
the Ru NP (inset) as a function of the illumination wavelength.

The temperature-dependent thermal conductivities of the input gases [57, 58,
59] are shown in Fig. 2(a). Since H2 has a small molecular mass and a small
molecule size, its thermal conductivity is much larger than that of N2 and Ar.
As a result, the thermal conductivity of the gas mixture is $\sim 60\%$ of the
H2 thermal conductivity. The thermal conductivity of the catalyst sample is
related to the volume fraction of the oxide using the Maxwell Garnett equation
[54, 55, 20] (see Supplemental Information Section
LABEL:suppsec:thermal_simulations).

The absorption cross-section of the Ru NPs and the light penetration depth in
the sample is calculated using the permittivity of Ru at 300K from Ref. 60 and
Eq. (2), see Fig. 2(b). For wavelengths $300<\lambda<600$nm, the penetration
(skin) depth (Eq. (2)) is $\delta_{\textrm{skin}}<750\mu$m, indeed much
smaller than the sample thickness; thus, by Ref. 31, we do not expect the
numerical results to be sensitive to the exact parameters in this spectral
regime. However, for $\lambda>700$nm, we find that
$\delta_{\textrm{skin}}\gtrsim$ 1 mm, so that here the penetration (skin)
depth is only $2-3$ times thinner than the sample thickness, see Fig. 2(b);
accordingly, one may expect a slight sensitivity to the various parameters in
this regime. In addition, we assume that the transverse profile of the
illumination intensity is uniform (denoted as $I_{\textrm{inc}}$).

To obtain the temperature distribution, we perform a full simulation including
the catalyst sample, the sample holder, and the reaction chamber with a few
simplifications. Specifically, we simplified the complicated reaction chamber,
by assuming it has a cylindrical shape (1 cm height and 2 cm diameter) and
that the catalyst sample (having the same size as used in the experiment) is
placed on a steel holder, see Fig. 3(a) 111The geometry of the reactor
chamber, of the sample holder and of the catalyst sample were directly
obtained from the authors of Refs. 56, 61..

In Ref. 56, two experiments were performed. In the first, no resistance-
heating was used; the simulation results of the top surface temperature
$T_{1}$ for $h\approx 70$ W/(m${}^{2}\cdot$K) demonstrate an excellent match
to the experimental data, see Fig. 3(b) (the simulation result of the bottom
surface temperature $T_{2}$ is shown in Fig.
LABEL:suppfig:JieLiu_NanoLett_2019_temp_map_T1_T2(a)). In order to demonstrate
the actual level of the photothermal nonlinearity, we also perform a
simulation under the linear approximation, namely, we neglect the temperature
dependence of the thermal conductivity of the gas mixture. Fig. 3(b) shows
that the linear approximation overestimates the nonlinear solution by more
than $30\%$ for a temperature rise of $\sim 300^{\circ}$C; this is an
unusually large nonlinearity.

In the second experiment, the authors studied the effect of the temperature
gradient on the reaction rate and the dependence of the temperature gradient
on the illumination wavelength and intensity. To do that, they measured the
intensity-dependent $T_{1}$, $T_{2}$ and reaction rate using four different
light sources (UV, blue light, white light and NIR) and adjusted the resistive
heating such that the equivalent temperature222The equivalent temperature
$T_{\textrm{e}}$ is defined in Refs. 61, 56 through the relation
$e^{-E_{\textrm{a}}/k_{\textrm{B}}T_{\textrm{e}}}=\dfrac{1}{T_{2}-T_{1}}\int_{T_{1}}^{T_{2}}e^{-E_{\textrm{a}}/k_{\textrm{B}}T}dT$.
remains the same (325∘C) for all four light sources and all intensities [56].
Accordingly, we use the same simulation configuration to calculate the
temperature distribution. Since the description of the resistive heating
apparatus is not available in Ref. 56, we simply set the temperature of the
bottom of the catalyst sample to the reported $T_{2}$ as a constraint so as to
mimic the resistive heating. The simulation results for the top surface
temperature $T_{1}$ again demonstrate an excellent match with the experimental
data for blue, UV, and white light sources333The simulation results for UV and
white light sources are shown in Fig.
LABEL:suppfig:JieLiu_NanoLett_2019_temp_map_T1_T2(b) because they almost
overlap with the results for the blue light source. This is a direct
consequence of the weak sensitivity of the temperature distribution to the
illumination wavelength when the penetration depth is much smaller than the
sample thickness, conforming with the analysis in Ref. 31., see Fig. 3(c) and
Fig. LABEL:suppfig:JieLiu_NanoLett_2019_temp_map_T1_T2(b). We find that the
linear approximation deviates from the nonlinear solution by $\sim 10$% for a
temperature rise of $\sim 100^{\circ}$C. For the NIR light source, the
simulation results fit well the experimental data, see Fig.
LABEL:suppfig:JieLiu_NanoLett_2019_temp_map_T1_T2(b). As mentioned, in this
case, the light penetration depth is closer to the sample thickness (see Fig.
2(b)) so that the temperature distribution might become less insensitive to
the change of permittivity of the Ru NPs with temperature. In that regard, the
use of the Ru permittivity data at $300$K in the simulation is one possible
reason for the small mismatch between the simulation results and the
experimental data. Except for this minor discrepancy, the analysis here
reinforces the conclusion of Ref. 56 that the catalytic effect of the Ru-
Cs/MgO system on the ammonia synthesis reaction is purely thermal.

Figure 3: (Color online) (a) Details of the photocatalytic chamber and the
temperature distribution at the illumination intensity
$I_{\textrm{inc}}=5$W/cm2. (b) $T_{1}$ (blue solid symbols: experimental data
extracted from Ref. 56, blue solid line: COMSOL simulation, blue dotted line:
COMSOL simulation with the linear approximation) as a function of the
illumination intensity for the blue light source without resistive heating.
(c) Same as (b) but under resistive heating. The blue opened symbols represent
the experimental data of $T_{2}$.

An earlier work [61] by the same group employed a similar experimental system
to study the carbon dioxide hydrogenation reaction using a titanium oxide
supported Rhodium catalyst. Although the gas composition was different from
that in Ref. 56, the thermal conductivities of the gas mixtures [62, 58, 59]
were similar since the gas mixtures had a similar fraction of H2 in these two
works. We simulated the temperature for the experiment without resistive
heating at low temperatures (25∘C $<T_{1}<120^{\circ}$C) and found very good
agreement with the experimental data (see Supplemental Information Fig.
LABEL:suppfig:JieLiu_NanoLett_2018_T1_T2). This explains the observed
photothermal nonlinearity ($\sim$ 10%) for a similar temperature rise of $\sim
100^{\circ}$C, see Fig. S8 in Ref. 61. However, a much larger nonlinearity was
reported in the experiments with resistive heating at high temperatures
($250^{\circ}$C $<T_{1}<450^{\circ}$C); in fact, some of the data shows an
unusually large increase in $T_{1}$ for $300^{\circ}$C $\lesssim T_{1}\lesssim
375^{\circ}$C and then a small growth rate in $T_{1}$ for $T_{1}\gtrsim
375^{\circ}$C. Accordingly, it cannot be explained by our (nonlinear) thermal
model. Ruling out possible measurement artifacts in the data, it is natural to
advocate for the possibility of non-thermal electrons contributing to the
catalyzed reaction rate. However, a convincing interpretation of this
experiment might require a dedicated explanation of the origin of the large
nonlinearity.

###  4.2 Analysis of experiments from Zhou et al., Science, 2018, 362, 69 [17]

We now move on to show that the photothermal nonlinearity also explains the
experimental measurements described in Ref. 17. This work employed a similar
setup to study ammonia decomposition on a MgO-Al2O3 supported Cu-Ru catalyst,
however, unlike the papers studied in Section 4.1, the nonlinearity is much
stronger due to the higher temperature rises caused by the much smaller
thermal conductivity of NH3 [63, 64] (see Fig. 4(a)), as well as the higher
incident intensities used (see Fig. 4(c)). More importantly, while the authors
of [56] measured the sample temperature properly, significant concerns
regarding the validity of the temperature measurements in [17] have been
raised (see discussion in Refs. 18, 20, 19, 23, 22). Therefore, in what
follows, we rely only on a calculation of the temperature as well as on an
extraction of it from the reaction rate using the Arrhenius equation. These
approaches were shown to match well with the fitted temperatures, see Refs.
18, 20, 22.

Figure 4: (Color online) (a) Temperature dependence of the thermal
conductivity of NH3 [63, 64]. (b) Light penetration depth and the absorption
cross-section of the Cu-Ru NP (inset) as a function of the wavelength. (c)
Simulation results of the temperature at the center of the top surface as a
function of the illumination intensity (blue solid line) and its linear
approximation (blue dotted line); The black dashed line represents the upper
bound of the temperature range of the NH3 thermal conductivity measured in
Ref. 63. (d) A zoom-in of (c) at the low illumination region. (e) Calculated
reaction rates (the blue solid line and the blue dotted line) as a function of
the illumination intensity. The blue circles in (c), (d) and (e) represent the
numerical fit of the experimental data using the temperature-shifted Arrhenius
law.

We follow the procedure described in Ref. 20 to extract the temperature out of
the reaction rate. Specifically, we first obtain the activation energy 1.18eV
by fitting the experimental data in the dark (Fig. 1d in Ref. 17) to an
Arrhenius curve. Next, due to the failure to measure the light-induced
temperature rise in Ref. 17 (see 18, 19), we distinguish between the actual
temperature of the reactor (denoted as $T(I_{\textrm{inc}})$) and the
experimentally-measured temperature (denoted as $T_{\textrm{M}}$, see Fig.
S11D in Ref. 17). The difference between $T(I_{\textrm{inc}})$ and
$T_{\textrm{M}}$ was found to be
$T(I_{\textrm{inc}})=T_{\textrm{M}}+\tilde{a}I_{\textrm{inc}}+\tilde{b}I_{\textrm{inc}}^{2}$,
where $I_{\textrm{inc}}$ is in W/cm2, $T(I_{\textrm{inc}})$ and
$T_{\textrm{M}}$ are in K. The experimentally-measured temperature
$T_{\textrm{M}}$ vs. $I_{\textrm{inc}}$ is also fitted to a second-order
polynomial, giving
$T_{\textrm{M}}=298+80I_{\textrm{inc}}-3.8I_{\textrm{inc}}^{2}$. Then, we fit
the measured reaction rate data under illumination (Fig. 1d in Ref. 17) to a
temperature-shifted Arrhenius curve [20], i.e.,
$R(I_{\textrm{inc}})=R_{0}\exp\left(-\dfrac{E_{\textrm{a}}}{k_{B}T(I_{\textrm{inc}})}\right)$,
leading to
$T(I_{\textrm{inc}})=T_{\textrm{M}}+180I_{\textrm{inc}}-8I_{\textrm{inc}}^{2}$,
see Figs. 4(c), (d) and (e). This shows that the photothermal nonlinearity
becomes nearly $\sim 50\%$ at the highest intensity used. This also shows that
the temperature rise due to photon absorption becomes a sublinear function of
the illumination intensity at high temperatures (as predicted in Refs. 37, 38,
39 for a single metal nanoparticle) so that unlike the claims in Ref. 65 (see
p. 270 and on), the maximal temperature reached is $\sim 1400^{\circ}$C rather
than $\sim 2700^{\circ}$C. This result shows that claims that our thermal
model predicts temperatures which are unrealistically high are simply
incorrect.

The question remains - what is the reason for this massive sublinearity in the
temperature-rise-versus-intensity? To answer this question, we adapt the
simulation configuration used in Section 4.1 to the experimental setup
described in Ref. 17 and simulate the temperature distribution via the heat
transfer module of COMSOL Multiphysics. The thermal conductivity of NH3 is
taken from the experimental measurement of Refs. 63, 64 for the temperature
range of 52 - 652∘C 444Although the thermal conductivity of NH3 was measured
at the pressures of 12.9, 26.5, and 45.0 kN/m2 (0.127, 0.262 and 0.444 atm) in
Ref. 63, this data can be used for the simulation at 1 atm since the gas
thermal conductivity is very weakly-dependent on the pressure [66]. To justify
this, we compared the data from Ref. 63 with the data measured at 1 atm but in
a lower temperature range ($-33.6$ \- $426.9^{\circ}$C) provided in Ref. 64
and found good agreement between these two sets of data.. The measured data
were also fitted in Ref. 63 by the following cubic polynomial,
$\kappa_{\textrm{NH\textsubscript{3}}}(T)=5.237\times 10^{-4}+5.179\times
10^{-5}T+8.404\times 10^{-8}T^{2}+1.557\times 10^{-11}T^{3}$, where
$\kappa_{\textrm{NH\textsubscript{3}}}$ is in W/(m$\cdot$K) and $T$ is in K,
as shown in Fig. 4(a). This cubic polynomial is used in our simulations to
extrapolate the NH3 thermal conductivity for temperatures higher than 652∘C.

The catalyst sample was illuminated by a pulsed broadband white-light source
(Fianium, WL-SC-400-8, 400-900 nm, pulse duration 4 ps, repetition rate 80
MHz, and a 2 mm diameter beam profile on the catalyst sample surface) without
applying any resistive heating [17]. The penetration (skin) depth to the
sample is calculated using the permittivity data of Cu and Ru at 300K [60]. We
find that the penetration depth is much smaller than the sample thickness, see
Fig. 4(b). In addition to the temperature distribution, we also calculated the
reaction rate based on the Arrhenius equation using the calculated temperature
at the center of the top sample surface. Our simulation results show
satisfactory agreement with the fitted data, although they are independent of
the fit procedure. In fact, the agreement extends up to 1000∘C, i.e., even
beyond the expected bound for the validity of the used values for the thermal
conductivity of the host, see Figs. 4(c) and (d). In particular, the
temperature at the center of the top surface increases monotonically with the
illumination intensity but with a decreasing slope, see Fig. 4(c). In the
original manuscript [17], this sublinear growth of the temperature rise with
the illumination intensity was incorrectly attributed to the temperature
dependence of the thermal conductivity of the oxide support. This claim is,
however, invalid for two reasons. First, this claim would lead to an opposite
trend to that observed since the thermal conductivity of the MgO-Al2O3 support
decreases with temperature [67]. Second, the temperature dependence of the
thermal conductivity of the oxide support, in fact, was shown in Ref. 31 to
have a negligible effect on the overall nonlinearity because of the small
volume fraction of oxides in the catalyst sample (using the sample mass and
the sample volume reported in Ref. 17). Instead, in agreement with the
qualitative analysis and the modelling of the papers studied in Section 4.1,
the sublinear growth of the temperature rise with the illumination intensity
is mainly due to the increase of the gas thermal conductivity with temperature
[63, 64].

One of the possible reasons for the slight discrepancy observed for
$T>1000^{\circ}$C between the simulation results and the fitted data at high
intensities is the inaccuracy induced by the extrapolation to the thermal
conductivity data of NH3. This is because the third-order term ($T^{3}$) in
the cubic polynomial used for the extrapolation has the same order of
magnitude as the second-order term ($T^{2}$) for temperatures higher than
1000∘C, i.e., it is likely that higher-order terms are required for better
accuracy. Another possible but minor reason is the usage of the permittivity
data at 300K for Cu in the simulation555Only limited data is available, e.g.,
in Ref. 68, only the imaginary part of the permittivity for photon energy 1.8
- 6eV at 5 different temperatures (77K, 295K, 575K, 770K, and 920K) are
provided; in Ref. 69, the permittivity data are provided at only 3 different
temperatures (78K, 293K, and 423K).. However, since the light penetration
depth is quite smaller than the sample thickness, the change of the optical
properties of Cu could only have a minor effect on the overall photothermal
nonlinearity, see discussion in Ref. 31. Since comprehensive data of the
temperature dependence of the metal permittivity and of the gas thermal
conductivity hardly exist at such high-temperature regimes, the resolution of
this discrepancy requires further experimental study. Nevertheless, our
simulation result shows that the photothermal nonlinearity plays a non-
negligible role in reducing the growth rate of the rising temperature and,
hence, of the photocatalysis reaction rate.

##  5 Discussion

The analysis presented above was based on the initial modeling of the low-
temperature response. Once the unknown system parameters were determined by
fit to the experimental data, the nonlinear photothermal response observed
experimentally was modeled accurately using the known temperature dependence
of the various material constituents, i.e., using no additional fit
parameters. The success of the photothermal analysis shows that as
qualitatively predicted in Section 2, the nonlinear response in the
temperature rise originates from the temperature dependence of the (effective)
thermal conductivity of the host, and not from the response of the metal
itself.

This result also confirms the error in previous claims on a different source
of nonlinearity for these systems, and shows that claims that the thermal
model leads to unrealistic high temperatures are simply incorrect, see
discussion in Ref. 19. The metal nonlinearity could be of significance only
for very thin plasmonic catalysts or even on the single nanoparticle level
(see discussion in Refs. 37, 38, 39). These systems are, however, usually of
more fundamental rather than practical importance.

Unlike the generally weak effect of “hot” electrons, the photothermal
nonlinearity is very strong. Therefore, this effect must be quantified before
any claim for “hot” electron action can become convincing, and should not be
ignored even at low illumination intensities. In practice, the rather large
uncertainty in the magnitude of the associated nonlinear response coefficients
means that only “hot” electron effects which are clearly greater than this
uncertainty can be deduced. Unfortunately, satisfying these conditions poses a
severe constraint on claims for “hot” electron dominance.

While the effect of the rising temperature on the absorptivity of the metal
NPs may have a negligible effect on the overall temperature distribution, it
may have a significant effect on the thermal emissivity via the Kirchhoff Law
of Radiation, an effect already demonstrated experimentally [70]. This is
relevant for a correct determination of the temperature using thermal imaging
at mid-IR frequencies, i.e., the change of emissivity at those frequencies
would need to be accounted for at high temperatures. To the best of our
knowledge, this was not done so far in the context of plasmon-assisted
photocatalysis, see discussion in Ref. 21.

At the high temperatures at which significant photothermal nonlinearity may be
observed, the NPs themselves may undergo geometrical and morphological changes
and eventually may even melt (at a temperature which may be significantly
lower than the bulk melting temperature). This possibility was discussed in
great detail in Ref. 65 p. 270 and on. Briefly, while this effect may be
possible, it was not probed directly, and it is a-priori expected not to
affect the temperature distribution so much, especially not when
$\delta_{\textrm{skin}}\ll H$ [31]. Therefore, although melting may occur, it
does not affect the results discussed in the current work in a significant
manner.

Lastly, we note that system similar to those studied here were extensively
studied in the past in the context of composites with a high thermo-optical
nonlinearity, see e.g., Refs. 46, 47, 48, 49, 35 and in the context of various
applications such as optical limiting [50] or tunable optical devices [35,
71]. The main difference is that the particle density in those systems was
typically much lower, such that some (even most) of the light was transmitted
through the sample. In addition, the focus in these systems was on (nonlinear)
changes of the optical response (permittivity, transmission etc.) rather than
on the temperature rise (as above).

Thermal effects in these systems were usually ignored, and the optical
response was typically interpreted using a temporally- and spatially-local
response. In that sense, it would be intriguing to study the thermal response
in such systems, to see if (the temporally- and spatially-nonlocal) thermal
effects could explain some of the experimental observations, in particular,
the strong dependence on the spatial [72, 73, 74] and temporal extent of the
illumination [45, 75, 76] and on the host properties [77, 37, 38, 39]. In
these cases, changes to the host permittivity may be more important, because
of macroscopic transmission changes and thermal lensing effects, which were
the main motivation for these studies in the first place.

## Author Contributions

I.W. Un derived the formulation, performed the calculations, and wrote the
paper. Y.D. performed the initial fits to the experimental data and together
with Y.S. supervised the project and contributed to the writing of the paper.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

I.W.U. and Y.S. were supported by Israel Science Foundation (ISF) grant
(340/2020) and by Lower Saxony - Israel cooperation grant no. 76251-99-7/20
(ZN 3637).

## References

  * Nitzan and Brus [1981] A. Nitzan and L. E. Brus, _J. Chem. Phys._ , 1981, 75, 2205.
  * Chen and Osgood [1983] C. J. Chen and R. M. Osgood, _Phys. Rev. Lett._ , 1983, 50, 1705\.
  * Jin _et al._ [2001] R. Jin, Y. Cao, C. A. Mirkin, K. L. Kelly, G. C. Schatz and J. G. Zheng, _Science_ , 2001, 294, 1901–1903.
  * Maillard _et al._ [2003] M. Maillard, P. Huang and L. Brus, _Nano Letters_ , 2003, 3, 1611–1615.
  * Tian and Tatsuma [2005] Y. Tian and T. Tatsuma, _J. Am. Chem. Soc._ , 2005, 127, 7632\.
  * Ueno _et al._ [2008] K. Ueno, S. Juodkazis, T. Shibuya, Y. Yokota, V. Mizeikis, K. Sasaki and H. Misawa, _J. Am. Chem. Soc._ , 2008, 130, 6928–6929.
  * Tsuboi _et al._ [2009] Y. Tsuboi, R. Shimizu, T. Shoji and N. Kitamura, _J. Am. Chem. Soc._ , 2009, 131, 12623–12627.
  * Volpe _et al._ [2012] G. Volpe, M. Noack, S. S. Acímović, C. Reinhardt and R. Quidant, _Nano Lett._ , 2012, 12, 4864–4868.
  * Clavero [2014] C. Clavero, _Nat. Photonics_ , 2014, 8, 95–103.
  * Baffou and Quidant [2014] G. Baffou and R. Quidant, _Chem. Soc. Rev._ , 2014, 43, 3898.
  * Li and Valentine [2016] W. Li and J. Valentine, _Nanophotonics_ , 2016, 6, 177–191.
  * Chen and Wang [2021] K. Chen and H. Wang, _Molecular Systems Design and Engineering_ , 2021, 6, 250–280.
  * Linic _et al._ [2011] S. Linic, P. Christopher and D. Ingram, _Nat. Mater._ , 2011, 10, 911–921.
  * Mukherjee _et al._ [2013] S. Mukherjee, F. Libisch, N. Large, O. Neumann, L. V. Brown, J. Cheng, J. B. Lassiter, E. A. Carter, P. Nordlander and N. J. Halas, _Nano Lett._ , 2013, 13, 240–247.
  * Christopher _et al._ [2012] P. Christopher, H. Xin, A. Marimuthu and S. Linic, _Nat. Mater._ , 2012, 11, 1044–1050.
  * Mukherjee _et al._ [2014] S. Mukherjee, L. Zhou, A. Goodman, N. Large, C. Ayala-Orozco, Y. Zhang, P. Nordlander and N. J. Halas, _J. Am. Chem. Soc._ , 2014, 136, 64–67.
  * Zhou _et al._ [2018] L. Zhou, D. F. Swearer, C. Zhang, H. Robatjazi, H. Zhao, L. Henderson, L. Dong, P. Christopher, E. A. Carter, P. Nordlander and N. J. Halas, _Science_ , 2018, 362, 69.
  * Sivan _et al._ [2019] Y. Sivan, J. Baraban, I. W. Un and Y. Dubi, _Science_ , 2019, 364, eaaw9367.
  * Sivan _et al._ [2020] Y. Sivan, J. Baraban and Y. Dubi, _OSA Continuum_ , 2020, 3, 483–497.
  * Sivan _et al._ [2020] Y. Sivan, I. W. Un and Y. Dubi, _Chem. Sci._ , 2020, 11, 5017–5027.
  * Baffou _et al._ [2020] G. Baffou, I. Bordacchini, A. Baldi and R. Quidant, _Light: Sci. Appl._ , 2020, 9, 108.
  * Dubi _et al._ [2021] Y. Dubi, I. W. Un, J. Baraban and Y. Sivan, _accepted in Nature Catalysis, https://arxiv.org/abs/2105.06382_ , 2021.
  * Sivan and Dubi [2020] Y. Sivan and Y. Dubi, _Appl. Phys. Lett._ , 2020, 117, 130501\.
  * Dubi and Sivan [2019] Y. Dubi and Y. Sivan, _Light: Sci. Appl._ , 2019, 8, 89.
  * Sivan _et al._ [2019] Y. Sivan, I. W. Un and Y. Dubi, _Faraday Discuss._ , 2019, 214, 215–233.
  * Kamarudheen _et al._ [2018] R. Kamarudheen, G. W. Castellanos, L. P. J. Kamp, H. J. H. Clercx and A. Baldi, _ACS Nano_ , 2018, 12, 8447–8455.
  * Xu _et al._ [2020] X. Xu, A. Dutta, J. Khurgin, A. Wei, V. M. Shalaev and A. Boltasseva, _Laser Photonics Rev._ , 2020, 14, 1900376.
  * Yu and Jain [2019] S. Yu and P. K. Jain, _Nat. Commun._ , 2019, 10, 2022.
  * Dai _et al._ [2019] X. Dai, Q. Wei, T. Duong and Y. Sun, _ChemNanoMat_ , 2019, 5, 1000–1007.
  * Kamarudheen _et al._ [2020] R. Kamarudheen, G. Kumari and A. Baldi, _Nature Communications_ , 2020, 11, 3957.
  * Un and Sivan [2020] I. W. Un and Y. Sivan, _Nanoscale_ , 2020, 12, 17821–17832.
  * Un and Sivan [2021] I. W. Un and Y. Sivan, _ACS Photonics_ , 2021, 8, 1183–1190.
  * Li _et al._ [2019] X. Li, H. O. Everitt and J. Liu, _Nano Res._ , 2019, 19, 1706–1711.
  * Pensa _et al._ [2019] E. Pensa, J. Gargiulo, A. Lauri, S. Schlücker, E. Cortés and S. A. Maier, _Nano Lett._ , 2019, 19, 1867–1874.
  * Donner _et al._ [2015] J. S. Donner, J. Morales-Dalmau, I. Aldaa, R. Marty and R. Quidant, _ACS Photonics_ , 2015, 2, 355–360.
  * Chu _et al._ [2014] S.-W. Chu, T.-Y. Su, R. Oketani, Y.-T. Huang, H.-Y. Wu, Y. Yonemaru, M. Yamanaka, H. Lee, G.-Y. Zhuo, M.-Y. Lee, S. Kawata and K. Fujita, _Phys. Rev. Lett._ , 2014, 112, 017402.
  * Sivan and Chu [2017] Y. Sivan and S.-W. Chu, _Nanophotonics_ , 2017, 6, 317–328.
  * Gurwich and Sivan [2017] I. Gurwich and Y. Sivan, _Phys. Rev. E_ , 2017, 96, 012212.
  * Un and Sivan [2020] I. W. Un and Y. Sivan, _Phys. Rev. Mater._ , 2020, 4, 105201.
  * Wilson _et al._ [2012] R. B. Wilson, B. A. Apgar, L. W. Martin and D. G. Cahill, _Opt. Exp._ , 2012, 20, 28829 – 28838.
  * Reddy _et al._ [2017] H. Reddy, U. Guler, K. Chaudhuri, A. Dutta, A. V. Kildishev, V. M. Shalaev and A. Boltasseva, _ACS Photonics_ , 2017, 4, 1083–1091.
  * Reddy _et al._ [2016] H. Reddy, U. Guler, A. V. Kildishev, A. Boltasseva and V. M. Shalaev, _Optical Materials Express_ , 2016, 6, 2776–2802.
  * Shen _et al._ [2016] P.-T. Shen, Y. Sivan, C.-W. Lin, H.-L. Liu, C.-W. Chang and S.-W. Chu, _Opt. Exp._ , 2016, 24, 19254.
  * Osaka _et al._ [2016] Y. Osaka, S. Sugano and S. Hashimoto, _Nanoscale_ , 2016, 8, 18187–18196.
  * Stoll _et al._ [2014] T. Stoll, P. Maioli, A. Crut, N. D. Fatti and F. Vallée, _Eur. Phys. J. B_ , 2014, 87, 260.
  * Smith _et al._ [1997] D. D. Smith, G. Fischer, R. W. Boyd and D. A. Gregory, _J. Opt.Soc. Am. B_ , 1997, 14, 1625.
  * Boyd _et al._ [2014] R. W. Boyd, Z. Shi and I. D. Leon, _Opt. Comm._ , 2014, 326, 74–79.
  * Rashidi-Huyeh and Palpant [2006] M. Rashidi-Huyeh and B. Palpant, _Phys. Rev. B_ , 2006, 74, 075405\.
  * Khurgin and Sun [2013] J. B. Khurgin and G. Sun, _Opt. Exp._ , 2013, 21, 27460–27480.
  * Wang and Blau [2009] J. Wang and W. J. Blau, _J. Opt. A: Pure Appl. Opt._ , 2009, 11, 024001\.
  * Winsemius _et al._ [1975] P. Winsemius, M. Guerrisi and R. Rosei, _Phys. Rev. B_ , 1975, 12, 4570–4572.
  * Winsemius _et al._ [1976] P. Winsemius, F. F. van Kampen, H. P. Lengkeek and C. G. van Went, _J. Phys. F: Met. Phys._ , 1976, 12, 1583.
  * Gusarov and Kovalev [2009] A. V. Gusarov and E. P. Kovalev, _Phys. Rev. B_ , 2009, 80, 024202\.
  * [54] R. B. Bird, W. E. Stewart and E. N. Lightfoot, _Transport Phenomena 2nd edition,(2002)_.
  * Pietrak and Wísniewski [2014] K. Pietrak and T. S. Wísniewski, _J. Power Technol._ , 2014, 95, 14–24.
  * Li _et al._ [2019] X. Li, X. Zhang, H. O. Everitt and J. Liu, _Nano Lett._ , 2019, 19, 1706–1711.
  * Engineering ToolBox [2018] Engineering ToolBox, _Nitrogen - Thermal Conductivity_ , 2018.
  * Engineering ToolBox [2018] Engineering ToolBox, _Hydrogen - Thermal Conductivity_ , 2018.
  * Saxena and Chen [1975] S. C. Saxena and S. H. P. Chen, _Molecular Physics_ , 1975, 29, 1507–1519.
  * Adachi [2012] S. Adachi, _The Handbook on Optical Constants of Metals_ , World Scientific, 2012.
  * Zhang _et al._ [2018] X. Zhang, X. Li, M. E. Reish, D. Zhang, N. Q. Su, Y. Gutiërrez, F. Moreno, W. Yang, H. O. Everitt and J. Liu, _Nano Lett._ , 2018, 18, 1714–1723.
  * Engineering ToolBox [2018] Engineering ToolBox, _Carbon dioxide - Thermal Conductivity_ , 2018.
  * Afshar _et al._ [1981] R. Afshar, S. Murad and S. C. Saxena, _Chemical Engineering Communications_ , 1981, 10, 1–11.
  * Engineering ToolBox [2018] Engineering ToolBox, _Ammonia - Thermal Conductivity at Varying Temperature and Pressure_ , 2018.
  * Aizpurua _et al._ [2019] J. Aizpurua, F. Baletto, J. Baumberg, P. Christopher, B. d. Nijs, P. Deshpande, Y. Diaz Fernandez, L. Fabris, S. Freakley, S. Gawinkowski, A. Govorov, N. Halas, R. Hernandez, B. Jankiewicz, J. Khurgin, M. Kuisma, P. V. Kumar, J. Lischner, J. Liu, A. Marini, R. J. Maurer, N. S. Mueller, M. Parente, J. Y. Park, S. Reich, Y. Sivan, G. Tagliabue, L. Torrente-Murciano, M. Thangamuthu, X. Xiao and A. Zayats, _Faraday Discuss._ , 2019, 214, 245–281.
  * Sears and Salinger [1975] F. W. Sears and G. L. Salinger, _Thermodynamics, kinetic theory, and statistical thermodynamics_ , Addison-Wesley, Boston, MA, 3rd edn, 1975.
  * Hofmeister [2014] A. M. Hofmeister, _Physics and Chemistry of Minerals_ , 2014, 41, 361–371.
  * Pells and Shiga [1969] G. Pells and M. Shiga, _J. Phys. C_ , 1969, 2, 1835.
  * Johnson and Christy [1975] P. B. Johnson and R. W. Christy, _Phys. Rev. B_ , 1975, 11, 1315–1323.
  * Ariceta [2016] T. E. Ariceta, _PhD Thesis_ , 2016.
  * Donner _et al._ [2011] J. S. Donner, G. Baffou, D. McCloskey and R. Quidant, _ACS Nano_ , 2011, 5, 5457–5462.
  * Maznev _et al._ [2011] A. A. Maznev, J. A. Johnson and K. A. Nelson, _Appl. Phys._ , 2011, 109, 073517.
  * Block _et al._ [2019] A. Block, M. Liebel, R. Yu, M. Spector, Y. Sivan, J. G. de Abajo and N. F. van Hulst, _Science Advances_ , 2019, 5, eaav8965.
  * Sivan and Spector [2020] Y. Sivan and M. Spector, _ACS Photonics_ , 2020, 7, 1271–1279.
  * Masia _et al._ [2012] F. Masia, W. Langbein and P. Borri, _Phys. Rev. B_ , 2012, 85, 235403\.
  * Marini _et al._ [2013] A. Marini, M. Conforti, G. D. Valle, H. W. Lee, T. X. Tran, W. Chang, M. A. Schmidt, S. Longhi, P. S. J. Russell and F. Biancalana, _New. J. Phys._ , 2013, 15, 013033.
  * Baffou _et al._ [2010] G. Baffou, R. Quidant and F. J. G. de Abajo, _ACS Nano_ , 2010, 4, 709–716.