# Simultaneous false discovery bounds for invariant causal prediction

Jinzhou Li

###### Abstract

Invariant causal prediction (ICP, Peters et al., (2016)) provides a novel way
to identify causal predictors of a response by utilizing heterogeneous data
from different environments. One advantage of ICP is that it guarantees to
make no false causal discoveries with high probability. Such a guarantee,
however, can be too conservative in some applications, resulting in few or no
discoveries. To address this, we propose simultaneous false discovery bounds
for ICP, which provides users with extra flexibility in exploring causal
predictors and can extract more informative results. These additional
inferences come for free, in the sense that they do not require additional
assumptions, and the same information obtained by the original ICP is
retained. We demonstrate the practical usage of our method through simulations
and a real dataset.

##  1 Introduction

Discovering causal predictors of the response of interest is usually the
primary goal of scientific research. Based solely on observational data, the
causal relationship might be unidentifiable. In such cases, heterogeneous data
from different environments can be helpful for further identifiability and
thus is a valuable source for causal discovery. Under the multi-environments
setting, Peters et al., (2016) proposed a novel method, called invariant
causal prediction (ICP), to identify causal predictors of a response by
exploiting the invariance of the conditional distribution of the response
given its causal predictors across environments. Compared to other causal
discovery algorithms, one main advantage of ICP is that it provides a
statistical confidence guarantee for its output $\widehat{S}^{\text{ICP}}$:

$\displaystyle P\left(\widehat{S}^{\text{ICP}}\subseteq S^{*}\right)\geq
1-\alpha,$

where $\alpha\in(0,1)$ is a nominal level and $S^{*}$ denotes the set of true
causal predictors. That is, ICP ensures that all its discoveries are true
causal predictors with a probability larger than $1-\alpha$. This is also
known as the familywise error rate control guarantee.

The original ICP approach has been generalized to non-linear setting (Heinze-
Deml et al.,, 2018), sequential data (Pfister et al.,, 2019), and
transformation models (Kook et al.,, 2023). All these methods inherit the
familywise error rate control guarantee. However, controlling the familywise
error rate can be too conservative in many applications, especially when
causal predictors are highly correlated with non-causal ones. In such cases,
ICP may result in few or no causal discoveries. To address this, Heinze-Deml
et al., (2018) proposed the so-called defining sets, which are the smallest
sets guaranteeing to contain at least one true causal predictor. They are
especially useful in the case where ICP returns an empty set.

In this paper, we propose a method to obtain simultaneous false discovery
bound (Genovese and Wasserman,, 2006; Goeman and Solari,, 2011; Goeman et
al.,, 2021) for ICP. Specifically, let
$\mathcal{N}^{*}\subseteq[m]=\\{1,\dots,m\\}$ be the index set of non-causal
predictors among $m$ predictors, the simultaneous false discovery bound is a
function $t_{\alpha}:2^{[m]}\rightarrow\\{0,\dots,m\\}$, where $2^{m}$ denotes
the power set of $[m]$, such that

$\displaystyle P(|R\cap\mathcal{N}^{*}|\leq t_{\alpha}(R)\text{ for any
}R\subseteq[m])\geq 1-\alpha.$

As a result, one can freely check any set $R$ of interest, and a high
probability upper bound of the false discoveries can be obtained by
$t_{\alpha}(R)$. This contains the familywise error rate control and the
defining sets of Heinze-Deml et al., (2018) as special cases. In particular,
by searching for the largest set $R$ whose corresponding $t_{\alpha}(R)$ is
$0$, we obtain a set with familywise error rate control; by searching for the
smallest sets whose corresponding $t_{\alpha}(R)\leq|R|-1$, we obtain all
defining sets.

Extra causal information can be extracted using the simultaneous false
discovery bound. As a quick illustration, for a simulated dataset with nine
predictors shown in Section 4.1, no information is obtained by the original
ICP approach as it returns zero discoveries. However, more causal information
can be extracted by using our simultaneous false discovery bounds. For
example, we can know that with a probability larger than $0.95$, there is at
least one causal predictor in $\\{X_{3},X_{7}\\}$, at least two causal
predictors in $\\{X_{3},X_{5},X_{7},X_{8}\\}$, and at least four causal
predictors in all nine predictors. More simulations and a real data
application can be found in Section 4.1.

To obtain the simultaneous false discovery bound for ICP, we proceed by
considering a multiple testing problem (see (6)) that directly tests causal
predictors. We show that the original ICP approach is equivalent to directly
comparing certain p-values (see (8)) to the significance level. Then, by
generalizing the null hypothesis (6) and p-value (7) to (10) and (11),
respectively, we propose a simultaneous false discovery upper bound (12). This
upper bound can be seen as a special case of the closed testing method by
Goeman and Solari, (2011), but in our case, no closed testing adjustment is
needed due to the specific form of the p-values. In comparison to the original
ICP approach, our method provides additional causal information at the cost of
testing more hypotheses, without introducing any additional statistical
assumptions. We also discuss the idea of controlling the false discovery rate
for ICP.

##  2 A brief recap of invariant causal prediction

For the sake of simplicity, we use linear models (Peters et al.,, 2016) for
illustration. Consider the following linear structural equation model:

$Y^{e}\leftarrow\beta_{1}X^{e}_{1}+\cdots\beta_{m}X^{e}_{m}+\epsilon^{e}=\beta_{S^{*}}^{T}X^{e}_{S^{*}}+\epsilon^{e},$
(1)

where $S^{*}=\\{i\in[m]:\beta_{i}\neq 0\\}$ is the set of causal predictors of
$Y$, $e\in\mathcal{E}$ denotes one environment, $\epsilon^{e}\sim
F_{\epsilon}$, $\epsilon^{e}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X^{e}_{S^{*}}$, and $X^{e}_{i}$ can have
arbitrary distribution for different $e$. The goal is to estimate the set of
causal predictors $S^{*}$.

The key idea of ICP is to test whether the conditional distribution of $Y^{e}$
is invariant across environments given predictors $X^{e}_{S}$ for some
$S\subseteq[m]$. This leads to the following null hypothesis:

$\displaystyle H_{0,S}(\mathcal{E}):$ $\displaystyle\text{there exists a
distribution $F_{\epsilon}$ and a vector }\gamma\in\mathbb{R}^{m}\text{ with
support }S\text{ such that }$ (2) $\displaystyle\text{for all
}e\in\mathcal{E},Y^{e}=\gamma^{T}X^{e}+\epsilon^{e},\text{where
}\epsilon^{e}\sim F_{\epsilon}\text{ and
}\epsilon^{e}\mathchoice{\mathrel{\hbox
to0.0pt{$\displaystyle\perp$\hss}\mkern
2.0mu{\displaystyle\perp}}}{\mathrel{\hbox
to0.0pt{$\textstyle\perp$\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptstyle\perp$\hss}\mkern
2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox
to0.0pt{$\scriptscriptstyle\perp$\hss}\mkern
2.0mu{\scriptscriptstyle\perp}}}X^{e}_{S}.$

Under the assumptions that there is no latent variable and $\mathcal{E}$ does
not contain the environment where $Y$ is intervened, it is clear that
$H_{0,S^{*}}(\mathcal{E})$ is true. However, there may be other sets than
$S^{*}$ for which $H_{0,S}(\mathcal{E})$ is true, causing an identifiability
issue. To this end, Peters et al., (2016) defined the so-called identifiable
causal predictors under $\mathcal{E}$ as follows:

$\displaystyle S(\mathcal{E})=\bigcap_{S:H_{0,S}(\mathcal{E})\text{ is true
}}S.$ (3)

Note that $S(\mathcal{E})\subseteq S^{*}$ as $H_{0,S^{*}}(\mathcal{E})$ is
true.

ICP estimates the causal predictors by using the sample version of
$S(\mathcal{E})$ in two steps: (i) For every $S\subseteq[m]$, test
$H_{0,S}(\mathcal{E})$ at level $\alpha\in(0,1)$. (ii) Obtain the selected
causal predictors by

$\widehat{S}^{\text{ICP}}(\mathcal{E})=\underset{S:H_{0,S}(\mathcal{E})\text{
is not rejected }}{\bigcap}S.$ (4)

The familywise error rate guarantee of $\widehat{S}^{\text{ICP}}(\mathcal{E})$
holds at level $\alpha$ because

$\displaystyle P\left(\widehat{S}^{\text{ICP}}(\mathcal{E})\subseteq
S^{*}\right)=P\left(\underset{S:H_{0,S}(\mathcal{E})\text{ is not rejected
}}{\bigcap}S\subseteq S^{*}\right)\geq P(H_{0,S^{*}}(\mathcal{E})\text{ is not
rejected })\geq 1-\alpha.$ (5)

##  3 Simultaneous false discovery bounds for ICP

###  3.1 A multiple testing formulation of ICP

Instead of directly estimating $S(\mathcal{E})$ as in (4), another natural way
is to form a multiple testing problem:

$\displaystyle H^{*}_{0,i}:i\in\mathcal{N}^{*},\quad i\in[m],$ (6)

where $\mathcal{N}^{*}=[m]\setminus S^{*}$. To this end, we propose the
following p-value for $H^{*}_{0,i}$:

$\displaystyle p^{*}_{i}=\max_{S\subseteq[m]\setminus\\{i\\}}p_{S},$ (7)

where $p_{S}$ is a valid p-value for testing $H_{0,S}(\mathcal{E})$ (see (2)),
that is, $P_{H_{0,S}(\mathcal{E})}(p_{S}\leq c)\leq c$ for any $c\in(0,1)$.
Note that it can be obtained in the first step of ICP. All p-values in this
paper depend on the environment set $\mathcal{E}$, but we omit this dependence
for simplicity of notation. The validity of $p^{*}_{i}$ is shown in the
following proposition.

######  Proposition 3.1.

For any $c\in(0,1)$, $P_{H^{*}_{0,i}}(p^{*}_{i}\leq c)\leq c$.

To connect to the original ICP approach, we need to introduce a slightly more
complicated p-value. For a given $\tau\in(0,1)$, let

$\displaystyle\tilde{p}^{*}_{i}(\tau)=\max\\{p^{*}_{i},\Pi_{S\subseteq[m]}\mathbbm{1}_{p_{S}\leq\tau}\\}.$
(8)

Note that $\tilde{p}^{*}_{i}(\tau)$ is also valid because
$\tilde{p}^{*}_{i}(\tau)\geq p^{*}_{i}$. Then, as the following proposition
shows, the discovery set obtained by using ICP is equivalent to the one
obtained by directly comparing $\tilde{p}^{*}_{i}(\alpha)$ to the significance
level $\alpha$.

######  Proposition 3.2.

For $\alpha\in(0,1)$, let
$\tilde{S}(\mathcal{E})=\\{i:\tilde{p}^{*}_{i}(\alpha)\leq\alpha\\}$. Then
$\tilde{S}(\mathcal{E})=\widehat{S}^{\text{ICP}}(\mathcal{E})$.

Hence, $\tilde{S}(\mathcal{E})$ possesses the familywise error rate control
guarantee (see (5)). In fact, the discovery set

$\displaystyle\widehat{S}(\mathcal{E})=\\{i:p^{*}_{i}\leq\alpha\\}$ (9)

obtained by directly comparing $p^{*}_{i}$ to the significance level also
controls the familywise error rate, as shown in the following proposition.

######  Proposition 3.3.

For $\alpha\in(0,1)$, we have $P\left(\widehat{S}(\mathcal{E})\subseteq
S^{*}\right)\geq 1-\alpha.$

Therefore, both $\widehat{S}(\mathcal{E})$ and
$\widehat{S}^{\text{ICP}}(\mathcal{E})$ (equivalent to
$\tilde{S}(\mathcal{E})$) control the familywise error rate, and using
$\widehat{S}(\mathcal{E})$ seems better because
$p^{*}_{i}\leq\tilde{p}^{*}_{i}(\alpha)$. However, the scenario where
$p^{*}_{i}<\tilde{p}^{*}_{i}(\alpha)$ only occurs when $p_{S}\leq\alpha$ for
all $S\subseteq[m]$. In such cases, we have
$\widehat{S}^{\text{ICP}}(\mathcal{E})=\emptyset$ and
$\widehat{S}(\mathcal{E})=[m]$, which appears not meaningful. In practice,
$\widehat{S}(\mathcal{E})$ and $\widehat{S}^{\text{ICP}}(\mathcal{E})$ are
generally equivalent.

###  3.2 Simultaneous false discovery bounds for ICP

Based on the multiple testing formulation (6), we propose simultaneous false
bounds by considering all intersection hypotheses. Specifically, for a given
$S\subseteq[m]$, by generalizing (6) and (7), we consider null hypothesis

$\displaystyle H^{*}_{0,S}:i\in\mathcal{N}^{*}\text{ for all }i\in S$ (10)

and p-value

$\displaystyle p^{*}_{S}=\max_{I\subseteq[m]\setminus S}p_{I}.$ (11)

By using a similar argument as in Proposition 3.1, one can see that
$p^{*}_{S}$ is valid.

######  Proposition 3.4.

For any $c\in(0,1)$, $P_{H^{*}_{0,S}}(p^{*}_{S}\leq c)\leq c$.

For any set $R\subseteq[m]$, let $t_{\alpha}(R)$ be the size of the largest
subset of $R$ whose corresponding p-value is larger than $\alpha$. That is,

$\displaystyle t_{\alpha}(R)=\max\\{|I|:I\subseteq R\text{ and
}p^{*}_{I}>\alpha\\}.$ (12)

Then, $t_{\alpha}(R)$ is a simultaneous false discovery upper bound, as shown
in the following Theorem.

######  Theorem 3.1.

Let $\alpha\in(0,1)$, then

$\displaystyle P(|R\cap\mathcal{N}^{*}|\leq t_{\alpha}(R)\text{ for any
}R\subseteq[m])\geq 1-\alpha.$ (13)

Based on the simultaneous guarantee (13), one can freely check any set $R$ of
interest, and a high probability upper bound (12) on the false discoveries can
be calculated. Equivalently, $|R|-t_{\alpha}(R)$ is a simultaneous true
discovery lower bound:

$\displaystyle P(|R\cap S^{*}|\geq|R|-t_{\alpha}(R)\text{ for any
}R\subseteq[m])\geq 1-\alpha.$ (14)

The simultaneous upper bound (12) is the same as the upper bound proposed by
Goeman and Solari, (2011) for closed testing (Marcus et al.,, 1976). In fact,
our method can be seen as a special case of closed testing. Compared to the
standard closed testing procedure, the difference is that no closed testing
adjustments are needed in our case. This is due to the specific form of
$p^{*}_{S}$. In particular, if $p^{*}_{S}\leq\alpha$, we must have
$p^{*}_{S^{\prime}}\leq\alpha$ for all $S\subseteq S^{\prime}$. Thus, if a
hypothesis is locally rejected, it must be rejected by closed testing, so
there is no need to implement closed testing adjustments.

At last, we mention that calculating the false discovery bound (12) can be
computationally intractable for large $|R|$ (or large $m$). In fact,
computation is also a fundamental issue for the original ICP approach, which
requires testing $2^{|m|}$ hypotheses. Some tricks are suggested to deal with
this issue, including first reducing the size of the potential causal
predictors by using a pre-screening method such as Lasso. See Peters et al.,
(2016) for more discussions.

###  3.3 A discussion about false discovery rate control for ICP

With a multiple testing formulation (6) and valid p-values (7), an alternative
strategy to obtain a less conservative result than
$\widehat{S}^{\text{ICP}}(\mathcal{E})$ is to control the commonly used false
discovery rate (Benjamini and Hochberg,, 1995; Benjamini and Yekutieli,, 2001;
Lehmann and Romano,, 2022). Due to Proposition 3.3, however, directly
comparing $p^{*}_{1},\dots,p^{*}_{m}$ to the significance level already
controls the familywise error rate, so there is no point in further applying
false discovery rate control procedures on these p-values. This observation is
interesting because, typically, multiple testing corrections (e.g., Bonferroni
correction) on raw p-values are required to control the familywise error rate.
Thus, it is natural to ask: Can we obtain smaller valid p-values than
$p^{*}_{1},\dots,p^{*}_{m}$? If so, one may apply false discovery rate control
procedures to obtain more discoveries than ICP. However, in the following, we
argue that this approach does not seem promising.

Specifically, we show that if some $p^{\prime}_{i}$ satisfying

$\displaystyle{\mathbb{P}}(p^{\prime}_{i}\leq
p^{*}_{i})=1\quad\text{and}\quad{\mathbb{P}}_{H^{*}_{0,i}}(p^{\prime}_{i}<p^{*}_{i})>0,$
(15)

then it is not valid in general. In particular, suppose that we have an ideal
test for $H_{0,S}(\mathcal{E})$ that yields p-values $p_{S}\sim U(0,1)$ for
true nulls and $p_{S}=0$ for false nulls. Consider a setting where the only
true null hypothesis is $H_{0,S^{*}}(\mathcal{E})$. This is the desired
setting where all causal predictors are identifiable, and it can happen, for
example, when the environments are rich enough. The following proposition
shows that $p^{\prime}_{i}$ satisfying (15) is not valid in such settings.

######  Proposition 3.5.

Let $S^{*}$ be the set of true causal predictors. Consider the desired case
where only $H_{0,S^{*}}(\mathcal{E})$ is true, $p_{S^{*}}\sim U(0,1)$ and
$p_{S}=0$ for any $S\neq S^{*}$. Then, for $p^{\prime}_{i}$ satisfying (15),
there exists some $c^{*}\in(0,1)$ such that
$P_{H^{*}_{0,i}}(p^{\prime}_{i}\leq c^{*})>c^{*}$.

This result excludes the possibility of constructing $p_{i}^{\prime}$ based on
some analytic function, such that $p_{i}^{\prime}$ is strictly smaller than
$p_{i}^{*}$ in probability $1$ (recall that $p_{i}^{*}$ is obtained by a
maximum function (7)). Consequently, considering false discovery rate control
for ICP does not appear to be a promising idea.

##  4 Simulations and a real data application

###  4.1 Numerical simulations

We implement simulations to verify the simultaneous guarantee (13) empirically
and to illustrate the extra information one may obtain by using the
simultaneous false discovery bound compared to the original ICP approach. All
simulations were carried out in R, and the code is available at
https://github.com/Jinzhou-Li/ICPsimultaneousBounds.

We generate samples based on the following linear structural equation model:

$X^{e}\leftarrow BX^{e}+\epsilon^{e},$

where $X^{e}=(X^{e}_{1},\dots,X^{e}_{10})^{T}$, $B$ is a strict lower-
triangular matrix, $\epsilon^{e}\sim N_{10}(b^{e},\text{diag}(\sigma))$ and
$e\in\\{1,2,3,4,5\\}$. That is, we consider mean-shift intervention and five
environments. We treat $X_{6}$ as the response in our simulations.

We implement $500$ simulations. In each simulation, we randomly sample the
non-zero elements of matrix $B$ with a non-zero probability $0.8$. Their
values are then uniformly sampled from $1$ to $2$. Each entry of $\sigma$ is
uniformly sampled from $0.5$ to $1.5$. For the mean-shift intervention, we set
$b^{1}=0$. For other $b_{i}$, we randomly select $3$ non-zero entries
excluding $X_{6}$. That is, we consider $3$ intervened variables. The non-zero
entries of $b_{i}$ are sampled from $N(10,5^{2})$. We generate $100$ samples
for each environment and use a significance level $\alpha=0.05$.

For each generated dataset, we apply ICP and record the number of discoveries.
Based on the same p-values used by ICP, we calculate the simultaneous false
discovery upper bound $t_{\alpha}(R)$ (see (12)) and the true discovery lower
bound $|R|-t_{\alpha}$ (see (14)) for all $2^{9}=512$ sets.

Figure 1: The false discovery upper bounds and true discovery lower bounds for
all $512$ sets, as well as the size of each set. The first plot is the
averaged result over $500$ simulations, and the second and third plots are the
result based on two simulated datasets.

Over these $500$ simulations, the empirical probability of
$P(|R\cap\mathcal{N}^{*}|\leq t_{\alpha}(R)\text{ for any }R\subseteq[m])$ is
$0.988$, which is larger than $1-\alpha=0.95$, so the simultaneous guarantee
(13) empirically holds. The average number of causal predictors is $4.144$,
and the average number of discoveries of ICP is $0.52$. The first plot in
Figure 1 shows the average false discovery upper bounds and true discovery
lower bounds for all $512$ sets, as well as the size of each set. Note that
the x-axis is only the index for each set, and we didn’t present which sets
these indices refer to for simplicity. Compared to ICP, extra information can
be obtained by looking at these bounds. For example, by looking at the true
discovery lower bound of the set containing all predictors (with set index
$512$), we know that, on average, there are about three true causal predictors
in all nine predictors.

To better illustrate the usage of the simultaneous bounds in practice, we look
at the results of two simulated datasets rather than the averaged results.
These results are shown in the second and third plots of Figure 1,
respectively.

For the first simulated dataset, the true causal predictors are
$X_{1},X_{2},X_{3}$ and $X_{4}$. ICP discovers $X_{2}$ and $X_{3}$, meaning
that with a probability larger than $0.95$, $X_{2}$ and $X_{3}$ are causal
predictors. That is all the information we can get by using ICP. By looking at
the simultaneous bounds (see the second plot in Figure 1), however, extra
information can be obtained. For example, we can see that the sets with
indices $131$ and $257$ has true discovery lower bounds $3$ and $4$,
respectively. These two sets are $\\{X_{1},X_{2},X_{3},X_{4}\\}$ and
$\\{X_{1},X_{2},X_{3},X_{4},X_{5}\\}$. Hence, we know that with probability
larger than $0.95$, there are at least three causal predictors in
$\\{X_{1},X_{2},X_{3},X_{4}\\}$, and at least four causal predictors in
$\\{X_{1},X_{2},X_{3},X_{4},X_{5}\\}$. Note that the same information that
both $X_{2}$ and $X_{3}$ are true causal predictors can also be obtained,
because the true discovery lower bound for $\\{X_{2},X_{3}\\}$ (with set index
$13$) is two.

For the second simulated dataset, the true causal predictors are
$\\{X_{1},X_{2},X_{3},X_{4},X_{5}\\}$. ICP returns zero discoveries. In
particular, the corresponding p-values (see (7)) for predictors are
$1,1,1,1,\\\ 0.848,0.329,1,1,1$. Hence, one can not get much information about
the causal predictors by using ICP in this case. By looking at the
simultaneous bounds (see the third plot in Figure 1), more information can be
obtained. For example, we know that with a probability larger than $0.95$,
there are at least one causal predictor in $\\{X_{3},X_{7}\\}$ (with set index
$23$), at least two causal predictors in $\\{X_{3},X_{5},X_{7},X_{8}\\}$ (with
set index $164$), at least three causal predictors in
$\\{X_{1},X_{3},X_{5},X_{7},X_{8},X_{9}\\}$ (with set index $406$), and at
least four causal predictors in
$\\{X_{1},X_{2},X_{3},X_{4},X_{5},X_{7},X_{8},X_{9}\\}$ (with set index
$503$).

###  4.2 Real data application

We look at a real dataset used in Stock et al., (2003) about the educational
attainment of teenagers (Rouse,, 1995). This dataset contains information of
$4739$ students from approximately 1100 US high schools, including gender,
ethnicity, achievement test score, whether the father or mother is a college
graduate, whether the family owns their home, whether the school is in an
urban area, county unemployment rate, state hourly wage in manufacturing,
average college tuition, whether the family income above 25000 per year, and
region.

We follow Peters et al., (2016) to obtain two datasets of two environments,
with $2231$ and $2508$ samples, respectively. The response variable $Y$ is a
binary variable indicating whether the student attained a BA degree or higher.
We use dummy variables to encode the factors in the dataset, resulting in $13$
predictors in total. The goal is to find the causal predictors of the response
$Y$.

Figure 2: The false discovery upper bounds, true discovery lower bounds, and
the size of set for all $8192$ sets based on the real data.

We first apply ICP with a significance level $\alpha=0.05$, returning one
variable ‘score’ as the causal predictor. It seems reasonable that test score
has a causal influence on whether the student would attain a BA degree or
higher. We then calculate the simultaneous bounds for all $2^{13}=8192$ sets,
and the results are shown in Figure 2. More information can be obtained by
looking at these bounds. For example, by looking at the set with index $112$,
we know that with a probability larger than $0.95$, there are at least two
causal predictors in {‘score’, ’fcollege$\\_$no’, ’mcollege$\\_$no’}. It seems
plausible that the true discovery lower bound is $2$ for this set, because the
two variables of whether father or mother went to college are highly
correlated, so it is difficult to distinguish between whether both are causal
predictors or only one of them is.

## Acknowledgement

The author thanks Jelle Goeman and Nicolai Meinshausen for helpful discussions
and suggestions on the draft of this paper. The author gratefully acknowledges
support by the SNSF Grant P500PT-210978.

## References

  * Benjamini and Hochberg, (1995) Benjamini, Y. and Hochberg, Y. (1995).  Controlling the false discovery rate: a practical and powerful approach to multiple testing.  Journal of the Royal statistical society: series B (Methodological), 57(1):289–300.
  * Benjamini and Yekutieli, (2001) Benjamini, Y. and Yekutieli, D. (2001).  The control of the false discovery rate in multiple testing under dependency.  Annals of statistics, pages 1165–1188.
  * Genovese and Wasserman, (2006) Genovese, C. R. and Wasserman, L. (2006).  Exceedance control of the false discovery proportion.  Journal of the American Statistical Association, 101(476):1408–1417.
  * Goeman et al., (2021) Goeman, J. J., Hemerik, J., and Solari, A. (2021).  Only closed testing procedures are admissible for controlling false discovery proportions.  The Annals of Statistics, 49(2):1218–1238.
  * Goeman and Solari, (2011) Goeman, J. J. and Solari, A. (2011).  Multiple testing for exploratory research.  Statistical Science, 26(4):584–597.
  * Heinze-Deml et al., (2018) Heinze-Deml, C., Peters, J., and Meinshausen, N. (2018).  Invariant causal prediction for nonlinear models.  Journal of Causal Inference, 6(2).
  * Kook et al., (2023) Kook, L., Saengkyongam, S., Lundborg, A. R., Hothorn, T., and Peters, J. (2023).  Model-based causal feature selection for general response types.  arXiv preprint arXiv:2309.12833.
  * Lehmann and Romano, (2022) Lehmann, E. and Romano, J. P. (2022).  Multiple testing and simultaneous inference.  In Testing Statistical Hypotheses, pages 405–491. Springer.
  * Marcus et al., (1976) Marcus, R., Eric, P., and Gabriel, K. R. (1976).  On closed testing procedures with special reference to ordered analysis of variance.  Biometrika, 63(3):655–660.
  * Peters et al., (2016) Peters, J., Bühlmann, P., and Meinshausen, N. (2016).  Causal inference by using invariant prediction: identification and confidence intervals.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 78(5):947–1012.
  * Pfister et al., (2019) Pfister, N., Bühlmann, P., and Peters, J. (2019).  Invariant causal prediction for sequential data.  Journal of the American Statistical Association, 114(527):1264–1276.
  * Rouse, (1995) Rouse, C. E. (1995).  Democratization or diversion? the effect of community colleges on educational attainment.  Journal of Business & Economic Statistics, 13(2):217–224.
  * Stock et al., (2003) Stock, J. H., Watson, M. W., et al. (2003).  Introduction to econometrics, volume 104.  Addison Wesley Boston.

##  Appendix A Supplementary material

###  A.1 Proof of Proposition 1

######  Proposition A.1.

For any $c\in(0,1)$, $P_{H^{*}_{0,i}}(p^{*}_{i}\leq c)\leq c$.

###### Proof.

When $H^{*}_{0,i}$ is true, we have $S^{*}\subseteq[m]\setminus\\{i\\}$ as
$i\not\in S^{*}$, so $p^{*}_{i}=\max_{S\subseteq[m]\setminus\\{i\\}}p_{S}\geq
p_{S^{*}}$. Therefore, for any $c\in(0,1)$,

$\displaystyle P_{H^{*}_{0,i}}(p^{*}_{i}\leq c)\leq
P_{H^{*}_{0,i}}(p_{S^{*}}\leq c)\leq c,$

where the last inequality holds because $H_{0,S^{*}}(\mathcal{E})$ is true. ∎

###  A.2 Proof of Proposition 2

######  Proposition A.2.

For $\alpha\in(0,1)$, let
$\tilde{S}(\mathcal{E})=\\{i:\tilde{p}^{*}_{i}(\alpha)\leq\alpha\\}$. Then
$\tilde{S}(\mathcal{E})=\widehat{S}^{\text{ICP}}(\mathcal{E})$.

To prove Proposition 2, we first introduce the following lemma, which gives
the relation between
$\widehat{S}(\mathcal{E})=\\{i:p^{*}_{i}(\alpha)\leq\alpha\\}$ and
$\widehat{S}^{\text{ICP}}(\mathcal{E})$.

######  Lemma A.1.

Assume that there exists $R\subseteq[m]$ such that $H_{0,R}(\mathcal{E})$ is
not rejected. Let
$\widehat{S}(\mathcal{E})=\\{i:p^{*}_{i}(\alpha)\leq\alpha\\}$. Then,
$\widehat{S}(\mathcal{E})=\widehat{S}^{\text{ICP}}(\mathcal{E})$.

We first prove Lemma A.1.

###### Proof.

If $\widehat{S}(\mathcal{E})=\emptyset$, that is,
$p^{*}_{i}=\max_{S\subseteq[m]\setminus\\{i\\}}p_{S}>\alpha$ for all
$i\in[m]$, then we must have
$\widehat{S}^{\text{ICP}}(\mathcal{E})=\emptyset$. Otherwise, if there exists
some $j\in\widehat{S}^{\text{ICP}}(\mathcal{E})$, then for any $R$ whose
corresponding $p_{R}>\alpha$, we must have $j\in R$ by the definition of
$\widehat{S}^{\text{ICP}}(\mathcal{E})$. But
$p^{*}_{j}=\max_{S\subseteq[m]\setminus\\{j\\}}p_{S}>\alpha$ implies that
there exists some $R^{\prime}$ not containing $j$ and $p_{R^{\prime}}>\alpha$,
which leads to a contradiction.

If $\widehat{S}^{\text{ICP}}(\mathcal{E})=\emptyset$, by the assumption in the
proposition, there must exist some $R_{1}$ and $R_{2}$ such that
$p_{R_{1}}>\alpha$, $p_{R_{2}}>\alpha$ and $R_{1}\cap R_{2}=\emptyset$. Then
we must have $\widehat{S}(\mathcal{E})=\emptyset$. Otherwise, if there exists
some $j\in\widehat{S}(\mathcal{E})$, that is,
$p^{*}_{j}=\max_{S\subseteq[m]\setminus\\{j\\}}p_{S}\leq\alpha$, then in order
for $p_{R_{1}}>\alpha$ and $p_{R_{2}}>\alpha$, both $R_{1}$ and $R_{2}$ must
contain $j$, which contradicts the fact that $R_{1}\cap R_{2}=\emptyset$.

Now assume both $\widehat{S}(\mathcal{E})$ and
$\widehat{S}^{\text{ICP}}(\mathcal{E})$ are non-empty.

If $i\in\widehat{S}^{\text{ICP}}(\mathcal{E})$, that is, for any $S$ such that
$H_{0,S}(\mathcal{E})$ is not rejected (that is, $p_{S}>\alpha$), we have
$i\in S$. So for any $R\subseteq[m]\setminus\\{i\\}$, $H_{0,R}(\mathcal{E})$
must be rejected, so
$p^{*}_{i}=\max_{R\subseteq[m]\setminus\\{i\\}}p_{R}\leq\alpha$, which implies
that $i\in\widehat{S}(\mathcal{E})$.

If $i\in\widehat{S}(\mathcal{E})$, we have
$p^{*}_{i}=\max_{S\subseteq[m]\setminus\\{i\\}}p_{S}\leq\alpha$, which implies
that for any $S$ not containing $i$, $H_{0,S}(\mathcal{E})$ must be rejected.
So for any $R$ such that $H_{0,R}(\mathcal{E})$ is not rejected, we must have
$i\in R$, which implies that $i\in\widehat{S}^{\text{ICP}}(\mathcal{E})$. ∎

Now we prove Proposition 2.

###### Proof.

If there exists some $R\subseteq[m]$ such that $H_{0,R}(\mathcal{E})$ is not
rejected, that is, $p_{R}>\alpha$. Then
$\tilde{p}^{*}_{i}(\alpha)=\max\\{p^{*}_{i},\Pi_{S\subseteq[m]}\mathbbm{1}_{p_{S}\leq\alpha}\\}=\max\\{p^{*}_{i},0\\}=p^{*}_{i}$,
so
$\tilde{S}(\mathcal{E})=\widehat{S}(\mathcal{E})=\widehat{S}^{\text{ICP}}(\mathcal{E})$
by Lemma A.1.

Now we consider the remaining case that $H_{0,S}(\mathcal{E})$ is rejected for
all $S\subseteq[m]$, that is, $p_{S}\leq\alpha$ for all $S\subseteq[m]$. Then,
we have $\widehat{S}^{\text{ICP}}(\mathcal{E})=\emptyset$ by definition, and
$\tilde{p}^{*}_{i}(\alpha)=\max\\{p^{*}_{i},\Pi_{S\subseteq[m]}\mathbbm{1}_{p_{S}\leq\alpha}\\}=\max\\{p^{*}_{i},1\\}=1>\alpha$
for any $i\in[m]$. Hence
$\tilde{S}(\mathcal{E})=\emptyset=\widehat{S}^{\text{ICP}}(\mathcal{E})$. ∎

###  A.3 Proof of Proposition 3

######  Proposition A.3.

For $\alpha\in(0,1)$, we have $P\left(\widehat{S}(\mathcal{E})\subseteq
S^{*}\right)\geq 1-\alpha.$

###### Proof.

When $H_{0,S^{*}}(\mathcal{E})$ is not rejected, for any $i\not\in S^{*}$, we
have $p^{*}_{i}=\max_{S\subseteq[m]\setminus\\{i\\}}p_{S}\geq
p_{S^{*}}>\alpha$, so $i\not\in\widehat{S}(\mathcal{E})$. Hence, we have
$\widehat{S}(\mathcal{E})\subseteq S^{*}$. Therefore,

$P\left(\widehat{S}(\mathcal{E})\subseteq S^{*}\right)\geq
P\left(H_{0,S^{*}}(\mathcal{E})\text{ is not rejected}\right)\geq 1-\alpha.$

∎

###  A.4 Proof of Proposition 4

######  Proposition A.4.

For any $c\in(0,1)$, $P_{H^{*}_{0,S}}(p^{*}_{S}\leq c)\leq c$.

###### Proof.

When $H^{*}_{0,S}$ is true, we have $S^{*}\subseteq[m]\setminus S$ as $S\cap
S^{*}=\emptyset$, so $p^{*}_{S}=\max_{I\subseteq[m]\setminus S}p_{I}\geq
p_{S^{*}}$. Therefore, for any $c\in(0,1)$,

$\displaystyle P_{H^{*}_{0,S}}(p^{*}_{S}\leq c)\leq
P_{H^{*}_{0,S}}(p_{S^{*}}\leq c)\leq c.$

∎

###  A.5 Proof of Theorem 1

######  Theorem A.1.

Let $\alpha\in(0,1)$, then

$\displaystyle P(|R\cap\mathcal{N}^{*}|\leq t_{\alpha}(R)\text{ for any
}R\subseteq[m])\geq 1-\alpha.$

###### Proof.

For any $R\subseteq[m]$,

$\displaystyle
p^{*}_{R\cap\mathcal{N}^{*}}=\max_{I\subseteq[m]\setminus(R\cap\mathcal{N}^{*})}p_{I}\geq\max_{I\subseteq[m]\setminus\mathcal{N}^{*}}p_{I}=p^{*}_{\mathcal{N}^{*}}.$

So if $p^{*}_{\mathcal{N}^{*}}>\alpha$, we have
$p^{*}_{R\cap\mathcal{N}^{*}}>\alpha$, which implies that
$|R\cap\mathcal{N}^{*}|\leq t_{\alpha}(R)$ by definition. Hence

$P(|R\cap\mathcal{N}^{*}|\leq t_{\alpha}(R)\text{ for any }R\subseteq[m])\geq
P(p^{*}_{\mathcal{N}^{*}}>\alpha).$

In addition, we have

$\displaystyle
P(p^{*}_{\mathcal{N}^{*}}\leq\alpha)=P(\max_{I\subseteq[m]\setminus\mathcal{N}^{*}}p_{I}\leq\alpha)\leq
P(p_{S^{*}}\leq\alpha)\leq\alpha,$

or equivalently

$P(p^{*}_{\mathcal{N}^{*}}>\alpha)\geq 1-\alpha,$

which completes the proof. ∎

###  A.6 Proof of Proposition 5

######  Proposition A.5.

Let $S^{*}$ be the set of true causal predictors. Consider the desired case
where only $H_{0,S^{*}}(\mathcal{E})$ is true, $p_{S^{*}}\sim U(0,1)$ and
$p_{S}=0$ for any $S\neq S^{*}$. Then, for $p^{\prime}_{i}$ satisfying (15),
there exists some $c^{*}\in(0,1)$ such that
$P_{H^{*}_{0,i}}(p^{\prime}_{i}\leq c^{*})>c^{*}$.

###### Proof.

For any $i$ such that $H^{*}_{0,i}$ is true, we have
$p^{*}_{i}=\max_{S\subseteq[m]\setminus\\{i\\}}p_{S}=p_{S^{*}}$.

Since $P_{H^{*}_{0,i}}(p^{\prime}_{i}<p^{*}_{i})>0$, there must exist some
rational number $c^{*}\in(0,1)$ such that

$\displaystyle P_{H^{*}_{0,i}}(p^{\prime}_{i}\leq c^{*}<p^{*}_{i})>0.$ (16)

Then,

$\displaystyle P_{H^{*}_{0,i}}(p^{\prime}_{i}\leq c^{*})$ $\displaystyle=$
$\displaystyle P_{H^{*}_{0,i}}(p^{\prime}_{i}\leq
c^{*},p^{\prime}_{i}<p^{*}_{i})+P_{H^{*}_{0,i}}(p^{\prime}_{i}\leq
c^{*},p^{\prime}_{i}=p^{*}_{i})$ $\displaystyle=$ $\displaystyle
P_{H^{*}_{0,i}}(p^{\prime}_{i}\leq c^{*},p^{\prime}_{i}<p^{*}_{i}\leq
c^{*})+P_{H^{*}_{0,i}}(p^{\prime}_{i}\leq
c^{*}<p^{*}_{i})+P_{H^{*}_{0,i}}(p^{\prime}_{i}\leq
c^{*},p^{\prime}_{i}=p^{*}_{i})$ $\displaystyle>$ $\displaystyle
P_{H^{*}_{0,i}}(p^{*}_{i}\leq
c^{*},p^{\prime}_{i}<p^{*}_{i})+P_{H^{*}_{0,i}}(p^{*}_{i}\leq
c^{*},p^{\prime}_{i}=p^{*}_{i})$ $\displaystyle=$ $\displaystyle
P_{H^{*}_{0,i}}(p^{*}_{i}\leq c^{*})$ $\displaystyle=$ $\displaystyle
P_{H^{*}_{0,i}}(p_{S^{*}}\leq c^{*})$ $\displaystyle=$ $\displaystyle c^{*},$

which completes the proof. ∎


# Electron-phonon coupling and non-equilibrium thermal conduction in ultrafast
heating systems

Chuang Zhang<EMAIL_ADDRESS>Department of Physics, Hangzhou Dianzi
University, Hangzhou 310018, China  Rulei Guo  Department of Mechanical
Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-8656,
Japan  Meng Lian  School of Physics, Institute for Quantum Science and
Engineering and Wuhan National High Magnetic Field Center, Huazhong University
of Science and Technology, Wuhan 430074, China  Junichiro Shiomi
Corresponding author<EMAIL_ADDRESS>Department of Mechanical
Engineering, The University of Tokyo, 7-3-1 Hongo, Bunkyo, Tokyo 113-8656,
Japan  Institute of Engineering Innovation, The University of Tokyo, 7-3-1
Hongo, Bunkyo, Tokyo 113-8656, Japan

###### Abstract

The electron-phonon coupling in ultrafast heating systems is studied within
the framework of Boltzmann transport equation (BTE) with coupled electron and
phonon transport. A discrete unified gas kinetic scheme is developed to solve
the BTE, in which the electron/phonon advection, scattering and electron-
phonon interactions are coupled together within one time step by solving the
BTE again at the cell interface. Numerical results show that the present
scheme can correctly predict the electron-phonon coupling constant, and is in
excellent agreement with typical two-temperature model (TTM) and experimental
results in existing literatures and our performed time-domain
thermoreflectance technique. It can also capture the ballistic or thermal wave
effects when the characteristic length/time is comparable to or smaller than
the mean free path/relaxation time where the TTM fails. Finally, the electron-
phonon coupling in transient thermal grating geometry and Au/Pt bilayer metals
with interfacial thermal resistance is simulated and discussed. For the
former, heat flow from phonon to electron is predicted in both the ballistic
and diffusive regimes. For the latter, the reflected signal increases in the
early tens of picoseconds and then decreases with time after the heat source
is removed.

##  I INTRODUCTION

Ultrafast laser heating plays a more and more important role in industrial
manufacturing of micro/nano electronic devices, medical detection and the
exploration of frontier basic science, which involves the interactions between
various energy-carrying (quasi)particles Chen (2005); Kittel  _et al._ (1996).
One of the key factors in multiscale energy transfer and conversion is the
coupling between electron and phonon (lattice vibration) in solid materials
Fann  _et al._ (1992a); Karna  _et al._ (2023); Block  _et al._ (2019); Allen
(1987); Caruso and Novko (2022). The electron-phonon coupling at the
microscopic level can explain many macroscopic thermal and electrical
transport phenomena, including thermoelectric conversion Mao  _et al._ (2021);
Deng  _et al._ (2021), electrothermal power consumption and transfer in
semiconductor chips Pop  _et al._ (2004); Pop (2010) and so on.

Over the past decades, many macroscopic phenomenological heat conduction
models are developed to describe the electron-phonon coupling Carpene (2006);
Schoenlein  _et al._ (1987); Groeneveld  _et al._ (1992); Caruso and Novko
(2022); Waldecker  _et al._ (2016); Lu  _et al._ (2018); Zhou  _et al._
(2015). One of the most widely used theoretical models is the two-temperature
model (TTM) Schoenlein  _et al._ (1987); Sun  _et al._ (1994); Qiu and Tien
(1994); Lin  _et al._ (2008); Rethfeld  _et al._ (2002); Caruso and Novko
(2022) proposed by Anisimov $et~{}al.$ Anisimov  _et al._ (1974). In this
empirical model, the electron temperature $T_{e}$ and phonon temperature
$T_{p}$ are introduced individually and their interactions are represented by
a single phenomenological electron-phonon coupling constant $G$. Although it
is simple, it is widely used in ultrafast pump-prope experiments for
calculating the electron-phonon coupling constant in metals Groeneveld  _et
al._ (1995); Brorson  _et al._ (1987); Fujimoto  _et al._ (1984); Guo and Xu
(2014); Qiu and Tien (1992); Qiu  _et al._ (1994). The Fourier’s law is used
to express the evolution process of electron and phonon in the spatial and
temporal spaces so that the TTM is a parabolic two-step heat conduction
equation with infinite heat propagation speed Joseph and Preziosi (1989). In
order to remove this non-physical assumption, a hyperbolic two-temperature
model is developed by introducing a delay term of the derivative of heat flow
with respect to time for electron or phonon transport Qiu and Tien (1994);
Sobolev (2021), which is like an extension of Cattaneo equation Cattaneo
(1948). In addition, the non-thermal lattice model Waldecker  _et al._ (2016)
or multitemperature model Lu  _et al._ (2018); Sobolev (2021) is developed, in
which phonons with different modes or branches are not in local thermal
equilibrium and many lattice temperatures and electron-phonon coupling
constants are introduced to describe the complex physical interactions.

Although above phenomenological heat conduction models made great success in
describing ultrafast energy exchange, they are not suitable to capture the
highly non-equilibrium situations, for example when the system characteristic
length is comparable to or smaller than electron/phonon mean free path where
the diffusive transport is broken Pattamatta and Madnia (2009a); Caruso and
Novko (2022). To study the multiscale energy transfer in ultrafast heating
systems, the time-dependent Boltzmann transport equation (BTE) Chen  _et al._
(2006); Pattamatta and Madnia (2009b); Miao and Wang (2021); Tong and Bernardi
(2021); Wang  _et al._ (2016) with coupled electron and phonon thermal
transport becomes an optimal compromise between efficiency and accuracy, which
can capture the diffusive and ballistic thermal transport simultaneously.
Instead of directly describe the evolution of macroscopic temperatures, the
macroscopic fields are obtained by taking the moment of distribution function
within the framework of BTE, so that the key of BTE is tracing the evolution
of phonon/electron distribution function in seven-dimensional phase space
(time, position, momentum). Two kinds of numerical methods are usually used to
solve the BTE with coupled electron and phonon thermal transport. One is the
statistics Monte Carlo method Jacoboni and Reggiani (1983); Muthukunnil Joseph
and Cao (2022), and the other is the discrete ordinate method Adams and Larsen
(2002); Miao and Wang (2021) which discretizes the whole three-dimensional
phase spaces into a lot of small pieces. Both of them made progress in highly
non-equilibrium heat transfer problems. However the advection and scattering
are treated separately in single time step at the numerical discrete level, so
that they have large numerical dissipations in the diffusive regime.

In this work, non-equilibrium thermal conduction in ultrafast heating systems
is studied by the BTE accounting for electron-phonon coupling, and a discrete
unified gas kinetic scheme (DUGKS) Guo and Xu (2021); Zhang and Guo (2019) is
developed in which the electron/phonon advection, scattering and electron-
phonon interactions are coupled together within one time step by solving the
BTE again at the cell interface. Numerical results show that the present
results are in excellent agreement with the experiments and TTM results in the
diffusive regime, and can capture the ballistic or thermal wave effects.

The remainder of this article is organized as follows. The BTE for coupled
electron and phonon thermal transport and the DUGKS solver are introduced in
Sec. II and Sec. III, respectively. Results and discussions are shown in Sec.
IV. Finally, a conclusion is made in Sec. V.

##  II Boltzmann transport equation

The Boltzmann transport equation for coupled electron and phonon thermal
transport is Chen  _et al._ (2006); Pattamatta and Madnia (2009b); Wang  _et
al._ (2016)

$\displaystyle\frac{\partial f_{e}}{\partial t}+\bm{v}_{e}\cdot\nabla f_{e}$
$\displaystyle=Q_{e-e}+Q_{e-p}+w_{s}S_{e},$ (1) $\displaystyle\frac{\partial
f_{p}}{\partial t}+\bm{v}_{p}\cdot\nabla f_{p}$
$\displaystyle=Q_{p-p}+Q_{p-e},$ (2)

where the subscripts $e$ and $p$ represent electron and phonon, respectively.
Note that this subscript representation will be used throughout the rest of
this article. $f$ is the distribution function, $\bm{v}$ is the group
velocity, $S_{e}$ is the external heat source at the mesoscopic level, $w_{s}$
is the associated weight satisfying

$\displaystyle\int w_{s}d\bm{K}=1,$ (3)

where $\int d\bm{K}$ represents the integral over the whole first Brillouin
zone. In actual pump-probe experiments, the external heat source originates
from the complex and ultrafast photon-electron-phonon coupling process, which
is not discussed in detail in this work Karna  _et al._ (2023); Block  _et
al._ (2019).

The first term on the left-hand side of Eqs.(1,2) represents temporal
evolution of distribution function and the second term represents
electron/phonon advection. $Q_{e-e}$ and $Q_{p-p}$ represent the electron-
electron scattering and phonon-phonon scattering, respectively. $Q_{e-p}$ and
$Q_{p-e}$ are the electron scattering term and phonon scattering term
characterizing electron-phonon interaction, respectively. The actual
scattering processes between (quasi) particles are much complicated so that
they are reduced in relaxation form in this work,

$\displaystyle Q_{e-e}$ $\displaystyle=\frac{f_{e}^{eq}-f_{e}}{\tau_{e}},$ (4)
$\displaystyle Q_{p-p}$ $\displaystyle=\frac{f_{p}^{eq}-f_{p}}{\tau_{p}},$ (5)

where $\tau$ is the associated relaxation time, and $f^{eq}$ is the
equilibrium state satisfying the Fermi-Dirac and Bose-Einstein distribution
for electron and phonon Chen (2005); Kittel  _et al._ (1996), respectively,

$\displaystyle f_{e}^{eq}(T_{e})$
$\displaystyle=\frac{1}{\exp\left(\frac{\varepsilon-\mu}{k_{B}T_{e}}\right)+1},$
(6) $\displaystyle f_{p}^{eq}(T_{p})$
$\displaystyle=\frac{1}{\exp\left(\frac{\hbar\omega}{k_{B}T_{p}}\right)-1},$
(7)

where $T$ is the temperature, $\varepsilon$ is the electron energy level,
$\mu$ is the chemical potential, $k_{B}$ is the Boltzmann constant, $\hbar$ is
the planck constant reduced by $2\pi$, $\omega$ is the phonon angular
frequency.

Equations (1) and (2) can be written as in terms of energy density,

$\displaystyle\frac{\partial u_{e}}{\partial t}+\bm{v}_{e}\cdot\nabla u_{e}$
$\displaystyle=\frac{u_{e}^{eq}-u_{e}}{\tau_{e}}-w_{g}G(T_{e}-T_{p})$
$\displaystyle+w_{s}S,$ (8) $\displaystyle\frac{\partial u_{p}}{\partial
t}+\bm{v}_{p}\cdot\nabla u_{p}$
$\displaystyle=\frac{u_{p}^{eq}-u_{p}}{\tau_{p}}+w_{g}G(T_{e}-T_{p}),$ (9)

where

$\displaystyle u_{e}-u_{e}^{eq}(T_{\text{ref}})$
$\displaystyle=\left(f_{e}-f_{e}^{eq}(T_{\text{ref}})\right)\left(\varepsilon-
E_{f}\right)D_{e},$ (10) $\displaystyle u_{p}-u_{p}^{eq}(T_{\text{ref}})$
$\displaystyle=\left(f_{p}-f_{p}^{eq}(T_{\text{ref}})\right)\hbar\omega
D_{p},$ (11) $\displaystyle S$ $\displaystyle=S_{e}\left(\varepsilon-
E_{f}\right)D_{e},$ (12)

where $E_{f}$ is the Fermi energy, $D$ is the density of states,
$T_{\text{ref}}$ is the reference temperature. We assume $\mu\approx E_{f}$
when $T\ll E_{f}/k_{B}$. The electron-phonon coupling process in Eqs. (8,9) is
simplified as

$\displaystyle Q_{p-e}=-Q_{e-p}=w_{g}G(T_{e}-T_{p})$ (13)

by invoking an electron-phonon coupling parameter $G$, which represents the
inelastic scattering that relaxes the electron and phonon to thermal
equilibrium. $w_{g}$ is the associated weights satisfying $\int
w_{g}d\bm{K}=1$ and the interaction strength between phonons with various
modes and electrons with various energy level may be different. The energy
conservation is satisfied during the scattering process so that

$\displaystyle\int\frac{u_{e}^{eq}-u_{e}}{\tau_{e}}d\bm{K}$ $\displaystyle=0,$
(14) $\displaystyle\int\frac{u_{p}^{eq}-u_{p}}{\tau_{p}}d\bm{K}$
$\displaystyle=0.$ (15)

The macroscopic variables including the energy $U$ and heat flux $\bm{q}$ are
obtained by taking the moment of distribution function,

$\displaystyle U_{e}$ $\displaystyle=\int u_{e}d\bm{K},$ (16) $\displaystyle
U_{p}$ $\displaystyle=\int u_{p}d\bm{K},$ (17) $\displaystyle\bm{q}_{e}$
$\displaystyle=\int\bm{v}_{e}u_{e}d\bm{K},$ (18) $\displaystyle\bm{q}_{p}$
$\displaystyle=\int\bm{v}_{p}u_{p}d\bm{K}.$ (19)

The temperature is calculated by following constrictions with Newton iteration
method,

$\displaystyle\int u_{e}^{eq}(T_{e})d\bm{K}$ $\displaystyle=\int
u_{e}d\bm{K},$ (20) $\displaystyle\int u_{p}^{eq}(T_{p})d\bm{K}$
$\displaystyle=\int u_{p}d\bm{K}.$ (21)

Actually in non-equilibrium systems the temperature calculated by above
formulas is not the thermodynamic temperature but more like the symbol of
local energy density Chen (2005).

To simplify the computation, the frequency-independent assumption is used so
that $|\bm{v}|$ and $\tau$ are constants for a given temperature. The three-
dimensional materials are considered and the first Brillouin zone is assumed
isotropic just like a spherical surface with fixed radius. $w_{g}$ and $w_{s}$
are both set to be $1/(4\pi)$. The temperature-dependent specific heat is
introduced in Eqs. (10,11) Pattamatta and Madnia (2009b),

$\displaystyle u_{e}-u_{e}^{eq}(T_{\text{ref}})$
$\displaystyle\approx\frac{1}{4\pi}\int_{T_{\text{ref}}}^{T_{e}}C_{e}dT,$ (22)
$\displaystyle u_{p}-u_{p}^{eq}(T_{\text{ref}})$
$\displaystyle\approx\frac{1}{4\pi}\int_{T_{\text{ref}}}^{T_{p}}C_{p}dT,$ (23)

where $C_{e}=\partial\left(f_{e}^{eq}\left(\varepsilon-
E_{f}\right)D_{e}\right)/\partial T_{e}$,
$C_{p}=\partial\left(f_{p}^{eq}\hbar\omega D_{p}\right)/\partial T_{p}$.

##  III Discrete unified gas kinetic scheme

To solve Eqs. (8,9), the discrete unified gas kinetic scheme (DUGKS) Guo and
Xu (2021); Guo  _et al._ (2023) is used, which has made great success in
multiscale particle transport Guo and Xu (2021). The whole time space, spatial
space and first Brillouin zone are discretized into many small pieces, and
Eqs. (8,9) in integral form over a control volume $i$ from time $t_{m}$ to
$t_{m+1}=t_{m}+\Delta t$ can be written as follows,

$\displaystyle u_{e,i,n}^{m+1}-u_{e,i,n}^{m}+\frac{\Delta t}{V_{i}}\sum_{j\in
N(i)}\left(\bm{v}_{e}\cdot\mathbf{n}_{ij}u_{e,ij,n}^{m+1/2}A_{ij}\right)$
$\displaystyle=\frac{\Delta t}{2}\left(H_{i,n}^{m+1}+H_{i,n}^{m}\right),$ (24)
$\displaystyle u_{p,i,n}^{m+1}-u_{p,i,n}^{m}+\frac{\Delta t}{V_{i}}\sum_{j\in
N(i)}\left(\bm{v}_{p}\cdot\mathbf{n}_{ij}u_{p,ij,n}^{m+1/2}A_{ij}\right)$
$\displaystyle=\frac{\Delta t}{2}\left(F_{i,n}^{m+1}+F_{i,n}^{m}\right),$ (25)

where the trapezoidal quadrature is used for the time integration of the
scattering, electron-phonon coupling and heat source terms, while the mid-
point rule is used for the flux term.
$H=(u_{e}^{eq}-u_{e})/\tau_{e}-w_{g}G(T_{e}-T_{p})+w_{s}S$,
$F=(u_{p}^{eq}-u_{p})/\tau_{p}+w_{g}G(T_{e}-T_{p})$, $n$ represents the index
of discretized first Brillouin zone, $V_{i}$ is the volume of the cell $i$,
$N(i)$ denotes the sets of neighbor cells of cell $i$, $ij$ denotes the
interface between cell $i$ and cell $j$, $A_{ij}$ is the area of the interface
$ij$, $\mathbf{n}_{ij}$ is the normal unit vector of the interface $ij$
directing from cell $i$ to cell $j$, $\Delta t$ is the time step and $m$ is an
index of time step. Reformulate above two equations (24,25) as

$\displaystyle\tilde{I}_{e,i,n}^{m+1}$
$\displaystyle=\tilde{I}_{e,i,n}^{+,m}-\frac{\Delta t}{V_{i}}\sum_{j\in
N(i)}\left(\bm{v}_{e}\cdot\mathbf{n}_{ij}u_{e,ij,n}^{m+1/2}A_{ij}\right),$
(26) $\displaystyle\tilde{I}_{p,i,n}^{m+1}$
$\displaystyle=\tilde{I}_{p,i,n}^{+,m}-\frac{\Delta t}{V_{i}}\sum_{j\in
N(i)}\left(\bm{v}_{p}\cdot\mathbf{n}_{ij}u_{p,ij,n}^{m+1/2}A_{ij}\right),$
(27)

where

$\displaystyle\tilde{I}_{e}$ $\displaystyle=u_{e}-\Delta tH/2,$ (28)
$\displaystyle\tilde{I}_{e}^{+}$ $\displaystyle=u_{e}+\Delta tH/2,$ (29)
$\displaystyle\tilde{I}_{p}$ $\displaystyle=u_{p}-\Delta tF/2,$ (30)
$\displaystyle\tilde{I}_{p}^{+}$ $\displaystyle=u_{p}+\Delta tF/2.$ (31)

In order to obtain the distribution function at the cell interface at the mid-
point time step ($u_{e,ij,n}^{m+1/2}$, $u_{p,ij,n}^{m+1/2}$), integrating Eqs.
(8,9) from time $t_{m}$ to $t_{m+1/2}=t_{m}+\Delta t/2$ along the
characteristic line with the end point $\bm{x}_{ij}$ locating at the center of
the cell interface $ij$ between cell $i$ and cell $j$,

$\displaystyle
u_{e}^{m+1/2}(\bm{x}_{ij})-u_{e}^{m}(\bm{x}_{ij}-\bm{v}_{e}\Delta t/2)$
$\displaystyle=\Delta
t/4\left(H^{m+1/2}(\bm{x}_{ij})+H^{m}(\bm{x}_{ij}-\bm{v}_{e}\Delta
t/2)\right),$ (32) $\displaystyle
u_{p}^{m+1/2}(\bm{x}_{ij})-u_{p}^{m}(\bm{x}_{ij}-\bm{v}_{p}\Delta t/2)$
$\displaystyle=\Delta
t/4\left(F^{m+1/2}(\bm{x}_{ij})+F^{m}(\bm{x}_{ij}-\bm{v}_{p}\Delta
t/2)\right).$ (33)

Reformulate above two equations as

$\displaystyle\bar{I}_{e}^{m+1/2}(\bm{x}_{ij})$
$\displaystyle=\bar{I}_{e}^{+,m}(\bm{x}_{e,ij}^{\prime}),$ (34)
$\displaystyle\bar{I}_{p}^{m+1/2}(\bm{x}_{ij})$
$\displaystyle=\bar{I}_{p}^{+,m}(\bm{x}_{p,ij}^{\prime}),$ (35)

where

$\displaystyle\bar{I}_{e}$ $\displaystyle=u_{e}-\Delta tH/4,$ (36)
$\displaystyle\bar{I}_{e}^{+}$ $\displaystyle=u_{e}+\Delta tH/4,$ (37)
$\displaystyle\bar{I}_{p}$ $\displaystyle=u_{p}-\Delta tF/4,$ (38)
$\displaystyle\bar{I}_{p}^{+}$ $\displaystyle=u_{p}+\Delta tF/4,$ (39)
$\displaystyle\bm{x}_{e,ij}^{\prime}$
$\displaystyle=\bm{x}_{ij}-\bm{v}_{e}\Delta t/2,$ (40)
$\displaystyle\bm{x}_{p,ij}^{\prime}$
$\displaystyle=\bm{x}_{ij}-\bm{v}_{p}\Delta t/2.$ (41)

$\bar{I}_{e}^{+,m}(\bm{x}_{e,ij}^{\prime})$ and
$\bar{I}_{p}^{+,m}(\bm{x}_{p,ij}^{\prime})$ are reconstructed by numerical
interpolation,

$\displaystyle\bar{I}_{e}^{+,m}(\bm{x}_{e,ij}^{\prime})=\bar{I}_{e}^{+,m}(\bm{x}_{c})+(\bm{x}_{e,ij}^{\prime}-\bm{x}_{c})\bm{\sigma}_{e,c},$
(42)
$\displaystyle\bar{I}_{p}^{+,m}(\bm{x}_{p,ij}^{\prime})=\bar{I}_{p}^{+,m}(\bm{x}_{c})+(\bm{x}_{p,ij}^{\prime}-\bm{x}_{c})\bm{\sigma}_{p,c},$
(43)

where $\bm{\sigma}_{c}$ is the spatial gradient of the distribution function
$\bar{I}^{+}$ at the cell $c$. If $\bm{v}\cdot\bm{n}_{ij}>0$, $c=i$; else
$c=j$. The first-order upwind scheme, van Leer limiter or least square method
can be adopted to calculate the spatial gradient to ensure the numerical
stability.

Once $\bar{I}_{e}^{m+1/2}$ and $\bar{I}_{p}^{m+1/2}$ are obtained, taking an
integral of Eqs. (36,38) over the whole first Brillouin zone leads to

$\displaystyle\frac{\Delta t}{4}S+\int\bar{I}_{e}d\bm{K}$
$\displaystyle=U_{e}+\frac{\Delta t}{4}G(T_{e}-T_{p}),$ (44)
$\displaystyle\int\bar{I}_{p}d\bm{K}$ $\displaystyle=U_{p}-\frac{\Delta
t}{4}G(T_{e}-T_{p}).$ (45)

Combining Eqs. (22,23,14,15) and above two equations, the electron and phonon
temperatures at the cell interface can be calculated by Newton iteration
method Zhang and Guo (2019). Then the original distribution function
$u_{e}^{m+1/2}$ and $u_{p}^{m+1/2}$ at the cell interface can be obtained
based on Eqs. (36,38). Similar treatment is conducted for the evolutions of
the distribution function and macroscopic variables at the cell center.
$\tilde{I}_{e}$ and $\tilde{I}_{p}$ can be updated based on Eqs. (26,27), then
taking an integral of Eqs. (28,30) over the whole first Brillouin zone leads
to

$\displaystyle\frac{\Delta t}{2}S+\int\tilde{I}_{e}d\bm{K}$
$\displaystyle=U_{e}+\frac{\Delta t}{2}G(T_{e}-T_{p}),$ (46)
$\displaystyle\int\tilde{I}_{p}d\bm{K}$ $\displaystyle=U_{p}-\frac{\Delta
t}{2}G(T_{e}-T_{p}).$ (47)

Combining Eqs. (22,23,14,15) and above two equations, the electron and phonon
temperatures at the cell center can be calculated by Newton iteration method.
Then the original distribution function $u_{e}^{m+1}$ and $u_{p}^{m+1}$ at the
cell interface can be updated.

Above procedures are the main evolution processes of the phonon and electron
distribution function in the DUGKS. Actually the numerical evolutions process
of DUGKS is not limited to frequency-independent and isotropic assumptions.
Related extension will be conducted in the future. The key difference between
the DUGKS and typical Monte Carlo method or discrete ordinate method is the
reconstruction of the distribution function at the cell interface. Instead of
direct numerical interpolation, the BTE is solved again at the cell interface
(Eqs. 32,33) so that the physical evolution process is included adaptively and
the electron/phonon advection, scattering and electron-phonon interactions are
coupled together within one time step Guo  _et al._ (2023); Guo and Xu (2021).

##  IV RESULTS AND DISCUSSIONS

In this section, numerical simulations are conducted and the present results
are compared to the experimental data or results predicted by macroscopic two-
temperature model (TTM, see appendix A). The time step of DUGKS is $\Delta
t=\text{CFL}\times\Delta x/v_{max}$, where $0<\text{CFL}<1$ is the
Courant–Friedrichs–Lewy number, $\Delta x$ is the minimum cell size and
$v_{max}$ is the maximum group velocity of electron and phonon. In the
following simulations, $\text{CFL}=0.4$ or $0.8$. The first-order upwind
scheme is used in the ballistic regime, and van Leer limiter is used in the
diffusive regime.

The remainder of this section is organized as follows. The first three
subsections are mainly to verify the effectiveness of the present scheme. In
subsection IV.1, the DUGKS results are compared with the experimental data in
existing references and the TTM. All input parameters and experimental data
are obtained from the existing references. In subsection IV.2, a time-domain
thermoreflectance (TDTR) experimental platform is used to measure the
reflected signals, and the experimental data is fitted by the DUGKS to get a
reasonable electron-phonon coupling constant. In subsection IV.3, the steady
quasi-1D cross-plane heat conduction is simulated and results show that the
present scheme could predict the thermal transport accurately even when the
cell size is much larger than the mean free path. In subsection IV.4, the
dynamics of electron and phonon in transient thermal grating geometry is
studied by DUGKS, and the thermal behaviors and differences in the ballistic
and diffusive regimes are discussed. In the final subsection IV.5, the heat
conduction in bilayer metals is studied accounted for interfacial thermal
resistance.

###  IV.1 Ultrafast laser heating for electron-phonon coupling constant

Figure 1: Ultrafast laser heating of single-layer metals.

(a) $L=100$ nm, front surface (b) $L=100$ nm, rear surface

(c) $L=200$ nm, front surface (d) $L=200$ nm, rear surface

Figure 2: Comparison of the normalized electron temperature predicted by the
present DUGKS results, two-temperature model (TTM) and experiments Brorson
_et al._ (1987); Qiu and Tien (1993) for single-layer Au film, where
$T_{e}^{*}=(T_{e}-T_{\text{ref}})/(T_{max}-T_{\text{ref}})$ and $T_{max}$ is
the maximum electron temperature.

(a) $P_{input}=17.6$ J$\cdot$m-2 (b) $P_{input}=70.6$ J$\cdot$m-2

Figure 3: Comparison of the normalized electron temperature predicted by the
present DUGKS results, two-temperature model (TTM) and experiments Guo and Xu
(2014) for single-layer Au film with fixed thickness $L=~{}1\mu$m and various
laser energy input $P_{input}$, where
$T_{e}^{*}=(T_{e}-T_{\text{ref}})/(T_{max}-T_{\text{ref}})$ and $T_{max}$ is
the maximum electron temperature.

The quasi-1D thermal transport of the single-layer Au films with different
thicknesses $L$ was studied, as shown in Fig. 1. The initial temperature in
the whole domain is $T_{\text{ref}}=300$ K and the ultrafast laser heating is
implemented on the front surface $x=0$,

$\displaystyle
S=\frac{P_{input}(1-R_{r})}{t_{s}d_{pump}}\exp{\left(-\frac{x}{d_{pump}}-\frac{(4\ln{2})t^{2}}{t_{s}^{2}}\right)},$
(48)

where $R_{r}$ is the optical reflectivity, $d_{pump}$ is the optical
penetration depth, $P_{input}$ is the total energy carried by a laser pulse
divided by the laser spot cross section, $t_{s}=t_{pump}+t_{h}$, $t_{h}$ is
the delay in the electron thermalization time after pulse absorption Hopkins
_et al._ (2011); Fann  _et al._ (1992b); Sun  _et al._ (1993); Hohlfeld  _et
al._ (2000) after pulse absorption, $t_{pump}$ is the full-width-at-half-
maximum (FWHM) duration of the laser pulse.

Thermal physical parameters of electrons in Au metals are listed below:
specific heat is $C_{e}=\gamma T_{e}$, $\gamma=71$ J$\cdot$m-3$\cdot$K-2,
group velocity is $|\bm{v}_{e}|=1.36\times 10^{6}$ m$\cdot$s-1, $E_{f}=5.51$
eV, thermal conductivity is Pattamatta and Madnia (2009b); Anisimov and
Rethfeld (1997)

$\displaystyle\kappa_{e}=\chi\frac{(\vartheta_{e}^{2}+0.16)^{5/4}(\vartheta_{e}^{2}+0.44)\vartheta_{e}}{(\vartheta_{e}^{2}+0.092)^{51/2}(\vartheta_{e}^{2}+\eta\vartheta_{p})},$
(49)

where $\chi=353$ W$\cdot$m-1$\cdot$K-1, $\eta=0.16$,
$\vartheta_{e}=k_{B}T_{e}/E_{f}$, $\vartheta_{p}=k_{B}T_{p}/E_{f}$. The
relaxation time is calculated by the kinetic relation
$\tau_{e}=3\kappa_{e}/(C_{e}|\bm{v}_{e}|^{2})$. Thermal physical parameters of
phonons in Au metals are listed below: specific heat is $C_{p}=2.5\times
10^{6}$ J$\cdot$m-3$\cdot$K-1, group velocity is $|\bm{v}_{p}|=2142.86$
m$\cdot$s-1, thermal conductivity is $\kappa_{p}=2.6$ W$\cdot$m-1$\cdot$K-1,
relaxation time is $\tau_{p}=3\kappa_{p}/(C_{p}|\bm{v}_{p}|^{2})$. The phonon
and electron mean free path are about $1.5$ nm and $33$ nm, respectively.

The transient heat conduction in single-layer Au film with different thickness
$L$ is simulated by DUGKS, and the dynamics of electron temperature at the
front $x=0$ and rear $x=L$ surfaces are plotted. In these simulations
Mozafarifard  _et al._ (2023); Qiu and Tien (1993), $t_{pump}=96$ fs,
$t_{h}=0$ fs, $d_{pump}=15.3$ nm, $R_{reflection}=0.93$, $P_{input}=10$
J$\cdot$m-2, $G=2.6\times 10^{16}$ W$\cdot$m-3$\cdot$K-1. The computational
domain is discretized with $20-80$ uniform cells, and $40$ discrete points in
the $|\bm{v}|\cos\theta$ direction is used to ensure the accuracy of the
numerical quadrature, where $\theta\in[0,\pi]$. The specular reflection
boundary conditions (Eq. (62)) are adopted for both the front and rear
surfaces. The evolutions of normalized electron temperature
$T_{e}^{*}=(T_{e}-T_{\text{ref}})/(T_{max}-T_{\text{ref}})$ are shown in Fig.
2, where $T_{max}$ is the maximum electron temperature and we assume that the
change in the reflected signal is linear to the change in electron temperature
Qiu and Tien (1993); Mozafarifard  _et al._ (2023). It can be found that the
present DUGKS results remain broadly consistent with the TTM and previous
experiments Brorson  _et al._ (1987); Qiu and Tien (1993).

The deviations between different models and experiments result from following
reasons. 1) Uncertainty of experimental parameters Guo and Xu (2014); Qiu and
Tien (1993); Sivan and Spector (2020); Sun  _et al._ (1993), for example, it
is very difficult to precisely determine the practical FWHM or optical
reflectivity of laser heating, which significantly influence the maximum
electron temperature rise. 2) The TTM assumes diffusive electron/phonon
transport with infinite heat propagation speed, which may be invalid when the
system characteristic length/time is comparable to or smaller than the mean
free path/relaxation time.

We also made a comparison with the experimental data measured by Guo
$et~{}al.$ Guo and Xu (2014) with fixed thickness $L=1~{}\mu$m and various
laser heating pump $P_{input}$. The electron ballistic depth length
$d_{ballistic}$ is introduced and the external heat source is

$\displaystyle S=$ $\displaystyle
P_{input}\frac{0.94(1-R_{r})}{t_{s}d_{length}\left(1-\exp(-L/d_{length})\right)}$
$\displaystyle\times\exp{\left(-\frac{x}{d_{length}}-(4\ln{2})\frac{t^{2}}{t_{s}^{2}}\right)},$
(50)

where $t_{s}=280$ fs, $d_{length}=d_{pump}+d_{ballistic}$, $d_{pump}=12.44$
nm, $d_{ballistic}=200$ nm, $R=0.970$. Thermal physical parameters of
electrons and phonon are the same as above except $\kappa_{p}=0.311$
W$\cdot$m-1$\cdot$K-1 and $G=1.5\times 10^{16}$ W$\cdot$m-3$\cdot$K-1.

The evolutions of normalized electron temperature $T_{e}^{*}$ are shown in
Fig. 3 with different laser heating input $P_{input}$, where the variations of
measured reflection signals are approximately linear to the electron
temperature increment Guo and Xu (2014). When $L=1~{}\mu$m, both electron and
phonon suffer a diffusive transport process, and the profiles show that the
DUGKS results are in excellent agreement with the TTM and experiments. When
$P_{input}=17.6$ J$\cdot$m-2, the maximum electron temperature predicted by
the DUGKS and TTM are $370$ K and $369.0$ K, respectively. When
$P_{input}=70.6$ J$\cdot$m-2, the maximum electron temperature predicted by
the DUGKS and TTM are $533$ K and $527.4$ K, respectively. Those results are
in consistent with the data shown in Ref. Guo and Xu (2014).

###  IV.2 Time-domain thermoreflectance experiments

(a) $t_{h}=0$ (b) $t_{h}=260$ fs

Figure 4: Comparison of the normalized reflected signal $\Delta R/R$ between
the TDTR experiments and DUGKS with various electron-phonon coupling constant
$G$ and delay in the electron thermalization time after pulse absorption
$t_{h}$ Fann  _et al._ (1992b); Sun  _et al._ (1993).

The transient reflected signal for Au thin film with thickness $18.2$ nm at
environment temperature $T_{\text{ref}}=288.15$ K was measured using a time-
domain thermoreflectance (TDTR) experimental platform, which is a widely used
thermal property measurement method based on pump-probe technology. In this
method, a modulated pulsed pump beam heats the surface periodically and a
delayed pulsed probe beam detects the variation in thermoreflectance signal,
which is a function of the temperature of electron and phonon. The signal
picked up by a photodiode and a lock-in amplifier is fitted with an analytical
heat conduction solution, such as two-temperature model and Boltzmann’s
transport equation. A Ti:Sapphire femtosecond laser (Chameleon II from
Coherent Inc.) is utilized to generate pulsed beam with a pulse width of $140$
fs. The wavelength of pump beam and probe beam are $400$ nm and $800$ nm,
respectively. These two beams are focused on the surface of sample by a $10$x
objective lens with diameters of $39.2~{}\mu$m and $10.9~{}\mu$m,
respectively. The delay time between probe beam and pump beam is controlled by
a motorized linear stage (PRO165SL(E)-600 from Aerotech Inc.) and four-fold
light path, which provides a minimum delay time step of $0.01$ ps.

The thermal physical parameters of electrons and phonons are the same as those
used in Fig. 2 as well as the heat source (48) except $t_{pump}=140$ fs and
$P_{input}(1-R_{r})=0.0112$ J$\cdot$m-2. The heating spot radius is much
larger than the film thickness so that the heat conduction can be approximated
as a quasi-1D transient heat conduction problem. The normalized reflected
signal $\Delta R/R$ predicted by DUGKS with fitted electron-phonon coupling
constant $G$ and the delay in the electron thermalization time after pulse
absorption $t_{h}$ Hopkins  _et al._ (2011); Fann  _et al._ (1992b); Ibrahim
_et al._ (2004); Chen  _et al._ (2006) are shown in Fig. 4, where the
relationship between the reflected signals and electron/phonon temperature is
introduced in Appendix C. It can be found that when $t_{h}=260$ fs and
$G\approx 2.0\times 10^{16}$ W$\cdot$m-3$\cdot$K-1, the numerical results are
in excellent agreement with the four group TDTR measured data.

There are many uncertainties compared to existing references: 1) different
experimental samples and external heat source, 2) different empirical formulas
and coefficients of the thermal physical parameters, for example the
temperature-dependent electron specific heat. Hence it can be found that the
measured electron-phonon coupling constant of Au in various references are
different. The electron-phonon coupling coefficient for Au thin film measured
by our performed TDTR is basically consistent with previous experimental data
Guo and Xu (2014); Brorson  _et al._ (1987); Qiu and Tien (1993); Sivan and
Spector (2020); Sun  _et al._ (1993) ranging from $1.5\times 10^{16}$ to
$4.0\times 10^{16}$ W$\cdot$m-3$\cdot$K-1.

###  IV.3 Steady cross-plane heat conduction

Figure 5: Spatial distributions of the electron and phonon temperatures with
different thickness $L$, where $T^{*}=(T-T_{R})/(T_{L}-T_{R})$, the red line
is the phonon temperature $T_{p}$ and the dark green is the electron
temperature $T_{e}$.

(a) $L=100$ nm

(b) $L=10~{}\mu$m

Figure 6: Spatial distributions of the electron and phonon temperatures with
different thickness $L$, where $T^{*}=(T-T_{R})/(T_{L}-T_{R})$, the red line
is the phonon temperature $T_{p}$ and the dark green is the electron
temperature $T_{e}$. ‘M=40’ and ‘M=100’ represent the discretized cell number.

The quasi-1D cross-plane heat conduction in Au metal film with different
thickness $L$ is studied. The wall temperature at each side of the film is
fixed at $T_{L}=T_{0}+\Delta T$ and $T_{R}=T_{0}-\Delta T$, respectively. The
thermal physical parameters of electrons and phonons are the same as those
used in Fig. 1 as well as the numerical discretizations. Spatial distributions
of the electron and phonon temperatures with different thickness $L$ are shown
in Fig. 5, where isothermal boundary conditions (Eq. (61)) are used for two
walls. It can be found that the temperature slip appears near the wall
boundaries, and the temperature slip of electron is more obvious because the
electron mean free path (33 nm) is larger than that of phonon (1.5 nm) in Au
metals. When the thickness increases, the electron-phonon coupling becomes
frequent so that the deviations between phonon and electron temperature
decreases.

We also make a test of grid independence in Fig. 6. It can be found that the
numerical results predicted by coarse grid (‘$M=40$’) are in excellent
agreement with that predicted by fine grid (‘$M=100$’). When the thickness
$L=10~{}\mu$m is much larger than mean free path, the temperature slip
disappears in the diffusive regime. Furthermore, the present scheme could
capture the temperature distributions accurately even the coarse cell size
($L/M=250$ nm) is much larger than the mean free path.

###  IV.4 Transient thermal grating

Table 1: Physical parameters of BTE with coupled electron and phonon transport in TTG geometry at $300$ K and $25$ K Pattamatta and Madnia (2009b); Freedman  _et al._ (2019); Block  _et al._ (2019, 2023). Au |  $300$ K |  $25$ K
---|---|---
$C_{p}$ (J$\cdot$m-3$\cdot$K-1) | 2.50E6 | 5.10E5
$\kappa_{p}$ (W$\cdot$m-1$\cdot$K-1) | 2.75 | 23.0
$|\bm{v}_{p}|$ (m$\cdot$s-1) | 2142.857 | 2142.857
$|\tau_{p}|$ (s) | 7.187E-13 | 2.946E-11
$\lambda_{p}=|\bm{v}_{p}||\tau_{p}|$ (m) | 1.540E-9 | 6.314E-8
$\alpha_{p}=\kappa_{p}/C_{p}$ (m2$\cdot$s-1) | 1.10E-6 | 4.51E-5
$C_{e}$ (J$\cdot$m-3$\cdot$K-1) | 2.0E4 | 1.70E3
$\kappa_{e}$ (W$\cdot$m-1$\cdot$K-1) | 320.0 | 977.0
$|\bm{v}_{e}|$ (m$\cdot$s-1) | 1.36E6 | 1.36E6
$|\tau_{e}|$ (s) | 2.595E-14 | 9.322E-13
$\lambda_{e}=|\bm{v}_{e}||\tau_{e}|$ (m) | 3.53E-8 | 1.268E-6
$\alpha_{e}=\kappa_{e}/C_{e}$ (m2$\cdot$s-1) | 0.016 | 0.575
$G$ (W$\cdot$m-3$\cdot$K-1) | 3.0E16 | 5.0E15

(a) $L=10$ nm (b) $L=100$ nm (c) $L=1~{}\mu$m

(d) $L=10$ nm (e) $L=100$ nm (f) $L=1~{}\mu$m

Figure 7: Time-dependent electron and phonon temperatures at $x=L/4$ with
different $L$ when $T_{\text{ref}}=300$ K, where $\Delta T=1$ K,
$T_{e}^{*}=(T_{e}-T_{\text{ref}})/\Delta T$,
$T_{p}^{*}=(T_{p}-T_{\text{ref}})/\Delta T$. Reference data originates from
Ref. Maznev  _et al._ (2011).

(a) $L=100$ nm (b) $L=1~{}\mu$m (c) $L=10~{}\mu$m

(d) $L=100$ nm (e) $L=1~{}\mu$m (f) $L=10~{}\mu$m

Figure 8: Time-dependent electron and phonon temperatures at $x=L/4$ with
different $L$ when $T_{\text{ref}}=25$ K, where $\Delta T=1$ K,
$T_{e}^{*}=(T_{e}-T_{\text{ref}})/\Delta T$,
$T_{p}^{*}=(T_{p}-T_{\text{ref}})/\Delta T$.

Transient thermal grating (TTG) Maznev  _et al._ (2011); Sivan and Spector
(2020) in a Au metal is analyzed. In quasi-1D geometry, the crossed pump
lasers produce a sinusoidal interference pattern on the surface of a sample
with period $L$, which results in a spatially sinusoidal temperature profile,

$\displaystyle T_{e}(x,t=0)$ $\displaystyle=T_{\text{ref}}+\Delta T\sin(2\pi
x/L),$ (51) $\displaystyle T_{p}(x,t=0)$ $\displaystyle=T_{\text{ref}},$ (52)

where $\Delta T$ is the temperature difference. We assume that the grating
initially heats only the electronic subsystem. The computational domain is
discretized with $20-50$ uniform cells, and $8-100$ discrete points in the
$|\bm{v}|\cos\theta$ direction is used. The periodic boundary conditions are
adopted for the domain. Detailed physical parameters of phonon and electron in
Au metal are listed in Table. 1, where $\alpha=\kappa/C$ is the thermal
diffusivity and $\lambda=|\bm{v}|\tau$ is the mean free path.

Numerical results of time-dependent electron and phonon temperatures at
$x=L/4$ with different spatial grating length $L$ are shown in Fig. 7, where
$T_{\text{ref}}=300$ K and $\Delta T=1$ K. It can be found that when
$L=1~{}\mu$m, the DUGKS results at room temperature are the same as those
obtained from Ref. Maznev  _et al._ (2011) or predicted by the TTM, namely,
the heat conduction is in the diffusive regime. Actually the BTE could recover
the TTM in the diffusive limit, see Appendix D. When the grating period length
decreases and becomes comparable to the electron mean free path, the results
predicted by the TTM deviate from the DUGKS results. The electron temperature
wave appears in the ballistic regime, which breaks the diffusion equation
exactly.

Another interesting thing is that initially the electron temperature is higher
than phonon temperature, so energy/heat is released from the electron to the
phonon, but after a while the phonon temperature is higher than the electron
temperature, namely, the energy is transferred back from the phonon to the
electron, as shown in Fig. 7(c,f). Similar heat flow from phonon to electron
$T_{e}-T_{p}<0$ also exist when $L=10$ nm or $100$ nm based on the DUGKS
results. This phenomenon $T_{e}-T_{p}<0$ has also been measured experimentally
Block  _et al._ (2023) or studied theoretically based on TTM Maznev  _et al._
(2011) at micron scale.

The underlying physical mechanisms of $T_{e}-T_{p}<0$ can be summarized as
below. Actually the heat conduction in electron or phonon subsystems is
composed of three parts: advection $\bm{v}\cdot\nabla u$, scattering
$(u^{eq}-u)/\tau$ and electron-phonon coupling $G(T_{e}-T_{p})$. The purpose
of advection and scattering is to make the temperature tend to be a constant
at different spatial locations, i.e., $T(\bm{x})\rightarrow T_{\text{ref}}$.
While the electron-phonon coupling aims to establish a local equilibrium
between phonon and electron subsystems for a given spatial position, i.e.,
$|T_{e}-T_{p}|\rightarrow 0$. In the diffusive regime, the phonon-phonon and
electron-electron scattering are sufficient so that the heat dissipations
depend on the thermal diffusivity $\alpha$ and electron-phonon coupling $G$.
In the ballistic regime, the phonon-phonon and electron-electron scattering
are rare so that the heat dissipations depend on the phonon/electron advection
speed $|\bm{v}|$ and electron-phonon coupling $G$. The thermal diffusivity and
group velocity of electron are much larger than those of phonon, which
indicates that the phonon subsystem is more like a energy accumulator relative
to the electron subsystem in both diffusive and ballistic regimes. Hence the
electron temperature varies between phonon temperature $T_{p}$ and reference
temperature $T_{\text{ref}}$.

When the reference temperature decreases from $300$ K to $25$ K, both thermal
diffusivity and mean free path increases. When $L\leq 1~{}\mu$m, electron
ballistic transport is obvious and temperature wave appears. The TTM results
deviate significantly from the DUGKS results, because the diffusive transport
is no longer valid due to the small grating period size. It can also be found
that the anomalous thermal phenomena $T_{e}-T_{p}<0$ exist. When
$L=10~{}\mu$m, the TTM results are in consistent with the DUGKS data, and the
thermal transport is in the diffusive regime.

Although the TTM is valid in the diffusive regime, the numerical results in
Fig. 8(c,f) are quite different from those in Fig. 7(c,f). In Fig. 7(c,f), the
electron-phonon coupling dominates heat dissipation at the initial stage so
that the electron temperature decreases significantly, and then the electron
and phonon temperatures decay slowly due to the small thermal diffusivity. In
Fig. 8(c,f), the thermal diffusivity of electron and phonon are much larger
than that at room temperature, so that it can be found that both electron and
phonon temperature tend to a constant in a shorter time.

###  IV.5 Thermal transport in a metallic bilayer

Figure 9: Schematic of the physical evolution process of electron/phonon transport in bilayer metals with interfacial thermal resistance.  Table 2: Thermal physical parameters of electrons and phonons in Pt metals Wang  _et al._ (2016) at $300$ K. $C_{p}$ (J$\cdot$m-3$\cdot$K-1) | 2.67E6
---|---
$\kappa_{p}$ (W$\cdot$m-1$\cdot$K-1) | 5.80
$|\bm{v}_{p}|$ (m$\cdot$s-1) | 1.91E3
$C_{e}$ (J$\cdot$m-3$\cdot$K-1) | $\gamma T_{e}$
$\gamma$ J$\cdot$m-3$\cdot$K-2 | $748.1$
$\kappa_{e}$ (W$\cdot$m-1$\cdot$K-1) | $\kappa_{0}T_{e}/T_{p}$
$\kappa_{0}$ (W$\cdot$m-1$\cdot$K-1) | $65.80$
$|\bm{v}_{e}|$ (m$\cdot$s-1) | 0.46E6
$G$ (W$\cdot$m-3$\cdot$K-1) | 108.80E16

(a)

(b)

Figure 10: Transient (a) temperature evolution and (b) reflected signal in
Pt/Au bilayer metals, where $\Delta T=T-T_{\text{ref}}$ and $\Delta R=a\Delta
T_{e}+b\Delta T_{p}$. For Pt, $a/b=0.25$.

(a) Au/Pt

(b) Au/Pt

Figure 11: Transient (a) temperature evolution and (b) reflected signal in
Au/Pt bilayer metals, where $\Delta T=T-T_{\text{ref}}$ and $\Delta R=a\Delta
T_{e}+b\Delta T_{p}$. For Au, $a/b=0.02$.

(a) $t=20\Delta t$ (b) $t=60\Delta t$

(c) $t=1000\Delta t$ (d) $t=10000\Delta t$

Figure 12: Spatial distributions of electron and phonon temperature in Au/Pt
at different moment, where $\Delta t=0.0294$ ps.

The heat conduction in the last four subsections only consider a single-layer
metal film. In the practical experiments or applications, the heat conduction
in bilayer or multilayer dissimilar materials is more common Giri  _et al._
(2015); Hopkins  _et al._ (2009); Karna  _et al._ (2023); Choi  _et al._
(2014). In this subsection, the thermal transport in a metallic bilayer (Au/Pt
or Pt/Au Choi  _et al._ (2014)) with ultrafast laser heating is studied by the
DUGKS and the interfacial thermal resistance between two dissimilar materials
are accounted Chen  _et al._ (2022); Wang  _et al._ (2016); Miao and Wang
(2023).

The diffuse mismatch model Chen  _et al._ (2022); Wang  _et al._ (2016); Miao
and Wang (2023), in which electron/phonon scattering at the metal-metal
interface is assumed to be completely diffuse, is used to describe the
interfacial thermal resistance. The transmittance $t$ and reflectance $r$ on
each side of the interface satisfy the following constrictions due to energy
conservation,

$\displaystyle r_{ij}+t_{ij}$ $\displaystyle=1,$ (53)

where $t_{ij}$ represents the transmittance from medium $i$ to medium $j$
across the interface, and $r_{ij}$ represents the reflectance in the medium
$i$ reflected back from the interface. In addition, at the thermal equilibrium
state (the temperature of the left and right side of interface is the same),
the net heat flux across the interface is zero due to the principle of
detailed balance so that

$\displaystyle t_{ij}C_{i}v_{i}=t_{ji}C_{j}v_{j},$ (54)

where $C_{i}$ and $v_{i}$ are the specific heat and group velocity of medium
$i$. We only consider the interaction between electron (phonon) in medium $i$
and electron (phonon) in medium $j$. The schematic of particle transport in
bilayer metals is shown in Fig. 9.

We consider a combination of a 23 nm Pt film and a 58 nm Au film with initial
environment temperature $T_{\text{ref}}=300$ K. The thermal physical
parameters of electrons and phonons in Au metals are shown in previous Tab. 1
with $\gamma=71.0$ J$\cdot$m-3$\cdot$K-2, and the thermal physical parameters
of electrons and phonons in Pt metals are shown in Tab. 2 Wang  _et al._
(2016).

Firstly, we simulate the transient heat conduction in Au/Pt bilayer metal, and
the external heat source is shown in Eq. (48) with $t_{pump}=400$ fs and
$P_{input}(1-R_{r})=1.60$ J$\cdot$m-2 Choi  _et al._ (2014). The specular
boundary conditions are adopted for the heating top surface regardless phonon
and electron. For the unheated bottom surface, the specular boundary condition
(Eq. (62)) is used for electron transport and the isothermal boundary
condition (Eq. (61)) with environment temperature $300$ K is used for phonon
transport, which indicates that the net electron heat flux at the bottom
surface is zero and the phonon heat flux is nonzero. The uniform cell is used
with cell size $1$ nm, and the CFL number is $0.80$. $40$ discrete points in
the $|\bm{v}|\cos\theta$ direction is used. Similarly, we also simulate the
transient heat conduction in Pt/Au bilayer metal.

Numerical results are shown in Fig. 10. It can be found that in Pt/Au metals,
the phonon temperature in Pt is always higher than the phonon temperature in
Au. The reflected signal increases firstly and then decreases gradually with
time. For Au/Pt metals the thermal behaviors are completely different as shown
in Fig. 11. It can be found that the phonon temperature in the unheated metal
is higher than the phonon temperature in the heated metal in the early tens of
picoseconds. The reflected signal is also non-monotonic with time after the
laser heating is almost removed. It increases with time in the early tens of
picoseconds and finally decreases with time.

In order to understand this anomalous heat conduction phenomenon, the spatial
distributions of electron and phonon temperature in Au/Pt bilayer metal at
different moment are plotted in Fig. 12. At the early stage $t=20\Delta t$,
the electron temperature in Au is highest. There are obvious temperature slip
or interfacial thermal resistance in the Au/Pt interface for both electron and
phonon transport. In addition, the phonon temperature in the Pt side is higher
than that in the Au side. That’s because the electron-phonon coupling of Pt is
much stronger than that of Au. When $t=60\Delta t$, the electron temperature
in Pt side near the Au/Pt interface is higher than the electron temperature in
Au side. When $t=1000\Delta t$, the phonon temperature in Au metal increases
and the lattice energy come from two part. One part is the electron-phonon
coupling in Au because the electron temperature is still higher than the
phonon temperature. The other is the phonon heat flows back from the Pt to Au
side, which leads to that the deviations of the phonon temperature on the each
side of Au/Pt interface decreases. When $t=10000\Delta t$, the electron-
phonon, electron-electron, phonon-phonon scattering are sufficient so that all
temperatures decreases with time gradually as shown in Fig. 11(a).

##  V Conclusion

Non-equilibrium thermal conduction in ultrafast heating systems is studied by
the BTE accounting for electron-phonon coupling. A discrete unified gas
kinetic scheme is developed to directly solve the BTE, in which the
electron/phonon advection, scattering and electron-phonon interactions are
coupled together within one time step by solving the BTE again at the cell
interface. Numerical results show that the present scheme not only correctly
predicts the heat conduction in the diffusive regime with coarse cell size,
but also captures the ballistic or thermal wave effects when the
characteristic length is comparable to or smaller than the mean free path
where the TTM fails. For ultrafast laser heating problem, the present results
are in excellent agreement with the experimental results in existing
literatures and our performed TDTR.

In transient thermal grating geometry, heat flow from phonon to electron is
predicted in both the ballistic and diffusive regimes. It results from the
competition of the thermal diffusivity and electron-phonon coupling in the
diffusive regime, and in the ballistic regime it results from the competition
of the phonon/electron advection and electron-phonon coupling. Furthermore, in
Au/Pt bilayer metals with interfacial thermal resistance, the predicted
reflected signal increases in the early tens of picoseconds and then decreases
with time after the heat source is removed. The stronger electron-phonon
coupling strength in the unheated Pt side results in the higher phonon
temperature in Pt, so that phonon heat flows back from Pt to Au side.

## Acknowledgments

This work is supported by the China Postdoctoral Science Foundation
(2021M701565). The authors are grateful to Prof. Yonatan Sivan in Ben-Gurion
University and Dr. Wuli Miao in Tsinghua University for communications on
electron-phonon coupling.

##  Appendix A Two-temperature model

Two-temperature model (TTM) SINGH (2010); Qiu and Tien (1994); Lin  _et al._
(2008); Rethfeld  _et al._ (2002) with coupled electron and phonon thermal
transport is

$\displaystyle C_{e}\frac{\partial T_{e}}{\partial t}$
$\displaystyle=\nabla\cdot(\kappa_{e}\nabla T_{e})-G(T_{e}-T_{p})+S,$ (55)
$\displaystyle C_{p}\frac{\partial T_{p}}{\partial t}$
$\displaystyle=\nabla\cdot(\kappa_{p}\nabla T_{p})+G(T_{e}-T_{p}).$ (56)

In the discretized time space, the diffusion term and heat source term are
dealt with the semi-implicit scheme and the electron-phonon coupling term is
implemented with explicit scheme,

$\displaystyle C_{e}^{m}\frac{\delta T_{e}^{m}}{\Delta
t}-0.5\nabla\cdot(\kappa_{e}^{m}\nabla\delta T_{e}^{m})$
$\displaystyle=\nabla\cdot(\kappa_{e}^{m}\nabla
T_{e}^{m})-G(T_{e}^{m}-T_{p}^{m})+\frac{S^{m}+S^{m+1}}{2},$ (57)
$\displaystyle C_{p}^{m}\frac{\delta T_{p}^{m}}{\Delta
t}-0.5\nabla\cdot(\kappa_{p}^{m}\nabla\delta T_{p}^{m})$
$\displaystyle=\nabla\cdot(\kappa_{p}^{m}\nabla
T_{p}^{m})+G(T_{e}^{m}-T_{p}^{m}),$ (58)

where $\delta T^{m}=T^{m+1}-T^{m}$ is the temperature increment for each time
step $m$. Actually the electron-phonon coupling term can also be implemented
with semi-implicit scheme,

$\displaystyle C_{e}^{m}\frac{\delta T_{e}^{m}}{\Delta
t}-0.5\nabla\cdot(\kappa_{e}^{m}\nabla\delta T_{e}^{m})+0.5G(\delta
T_{e}^{m}-\delta T_{p}^{m})$ $\displaystyle=\nabla\cdot(\kappa_{e}^{m}\nabla
T_{e}^{m})-G(T_{e}^{m}-T_{p}^{m})+\frac{S^{m}+S^{m+1}}{2},$ (59)
$\displaystyle C_{p}^{m}\frac{\delta T_{p}^{m}}{\Delta
t}-0.5\nabla\cdot(\kappa_{p}^{m}\nabla\delta T_{p}^{m})-0.5G(\delta
T_{e}^{m}-\delta T_{p}^{m})$ $\displaystyle=\nabla\cdot(\kappa_{p}^{m}\nabla
T_{p}^{m})+G(T_{e}^{m}-T_{p}^{m}).$ (60)

But the latter method will generate a larger coefficient matrix with coupled
electron and phonon temperature increment $G(\delta T_{e}^{m}-\delta
T_{p}^{m})$ compared to the former. We adopt the first strategy in this work
for the simplicity. To simulate the present quasi-1D numerical cases, $10-100$
uniform cells are used and the central scheme is used to discrete the
diffusion term in the spatial space.

##  Appendix B Boundary conditions

Boundary conditions are one of the key parts in the numerical simulations. Two
kinds of boundary conditions are used in this work.

  1. 1.

Thermalization/isothermal boundary assumes that the incident particles
(phonons/electron) are all absorbed by the boundary $\bm{x}_{b}$, and the
particles emitted from the boundary are the equilibrium state with the
boundary temperature $T_{b}$, i.e.,

$u(\bm{x}_{b})=u^{eq}(T_{b}),\quad\bm{v}\cdot\mathbf{n}_{b}>0,$ (61)

where $\mathbf{n}_{b}$ is the normal unit vector of the boundary pointing to
the computational domain.

  2. 2.

The specular reflecting boundary condition is

$\displaystyle
u(\bm{v})=u(\bm{v}^{\prime}),\quad\bm{v}^{\prime}\cdot\mathbf{n}_{b}<0,$ (62)

where $\bm{v}^{\prime}=\bm{v}-2\mathbf{n}_{b}(\bm{v}\cdot\mathbf{n}_{b})$ is
the incident direction.

##  Appendix C Temperature dependent thermoreflectance signals

The reflectivity signal $R$ Block  _et al._ (2019); Hopkins  _et al._ (2011)
for a thin Au film with thickness $h$ and refractive index
($n_{2}=\sqrt{\epsilon}$) under perpendicular incidence from air ($n_{1}=1$)
and supported on a sapphire substrate ($n_{3}=1.7$) is $R=|r|^{2}$, where

$\displaystyle r$
$\displaystyle=\frac{r_{12}+r_{23}\cdot\exp(2i\beta)}{1+r_{12}r_{23}\cdot\exp(2i\beta)},$
(63) $\displaystyle\beta$ $\displaystyle=2\pi n_{2}h/\lambda_{0},$ (64)
$\displaystyle r_{jk}$ $\displaystyle=\frac{n_{j}-n_{k}}{n_{j}-n_{k}},$ (65)

where $i$ is the imaginary number. $\lambda_{0}$ is the wave length of probe
pulse. The permittivity or dielectric function $\epsilon$ of Au thin film is
related to the electron temperature $T_{e}$ and phonon temperature $T_{p}$,
i.e.,

$\displaystyle\epsilon$
$\displaystyle=\epsilon_{\infty}-\frac{\omega_{p}^{2}(T_{p})}{\omega_{0}(\omega_{0}+i\gamma_{re}(T_{e},T_{p}))},$
(66) $\displaystyle\gamma_{re}(T_{e},T_{p})$
$\displaystyle=A_{ee}T_{e}^{2}+B_{ep}T_{p},$ (67)

where $\omega_{p}\approx 1.37\times 10^{16}$ rad/s is the plasma frequency
Kittel  _et al._ (1996); Guo  _et al._ (2012), $\epsilon_{\infty}=9.50$,
$A_{ee}=1.77\times 10^{7}$ K-2 $\cdot$s-1, $B_{ep}=1.45\times 10^{11}$ K-1
$\cdot$s-1. $\omega_{0}$ is the angular frequency of probe pulse. The
nonlinear thermoreflectance model is complicated, hence in many references the
change in the reflected signal is usually assumed to be linear with the
increase in temperature.

##  Appendix D Dimensional analysis

Make a dimensional analysis of the BTE (8,9) without external heat source,

$\displaystyle\frac{\partial u_{e}^{*}}{\partial
t^{*}}+|\bm{v}_{e}^{*}|\bm{s}\cdot\nabla_{\bm{x}^{*}}u_{e}^{*}$
$\displaystyle=\frac{u_{e}^{eq,*}-u_{e}^{*}}{\tau_{e}^{*}}-\frac{(T_{e}^{*}-T_{p}^{*})/(4\pi)}{1/G^{*}},$
(68) $\displaystyle\frac{\partial u_{p}^{*}}{\partial
t^{*}}+|\bm{v}_{p}^{*}|\bm{s}\cdot\nabla_{\bm{x}^{*}}u_{p}^{*}$
$\displaystyle=\frac{u_{p}^{eq,*}-u_{p}^{*}}{\tau_{p}^{*}}+\frac{(T_{e}^{*}-T_{p}^{*})/(4\pi)}{1/G^{*}},$
(69)

where $\bm{s}$ is the unit directional vector and the thermal physical
parameters of electrons are regarded as the reference variables
$v_{\text{ref}}=|\bm{v}_{e}|$, $C_{\text{ref}}=C_{e}$, so that

$\displaystyle t^{*}$ $\displaystyle=\frac{t}{t_{\text{ref}}},\quad$
$\displaystyle G^{*}$
$\displaystyle=\frac{t_{\text{ref}}}{C_{\text{ref}}/G},\quad$
$\displaystyle\bm{x}^{*}$ $\displaystyle=\frac{\bm{x}}{L_{\text{ref}}},$ (70)
$\displaystyle u_{p}^{*}$
$\displaystyle=\frac{u_{p}}{C_{\text{ref}}T_{\text{ref}}},\quad$
$\displaystyle\tau_{p}^{*}$
$\displaystyle=\frac{\tau_{p}}{t_{\text{ref}}},\quad$
$\displaystyle\bm{v}_{p}^{*}$
$\displaystyle=\frac{\bm{v}_{p}}{v_{\text{ref}}},$ (71) $\displaystyle
u_{e}^{*}$ $\displaystyle=\frac{u_{e}}{C_{\text{ref}}T_{\text{ref}}},\quad$
$\displaystyle\tau_{e}^{*}$
$\displaystyle=\frac{\tau_{e}}{t_{\text{ref}}},\quad$
$\displaystyle\bm{v}_{e}^{*}$
$\displaystyle=\frac{\bm{v}_{e}}{v_{\text{ref}}},$ (72) $\displaystyle
T_{p}^{*}$ $\displaystyle=\frac{T_{p}}{T_{\text{ref}}},\quad$ $\displaystyle
T_{e}^{*}$ $\displaystyle=\frac{T_{e}}{T_{\text{ref}}},\quad$ $\displaystyle
t_{\text{ref}}$ $\displaystyle=\frac{L_{\text{ref}}}{v_{\text{ref}}},$ (73)

where $L_{\text{ref}}$ is the system characteristic size. When
$t^{*}\gg\tau_{e}^{*}$ and $x^{*}\gg|\bm{v}_{e}^{*}|\tau_{e}^{*}$, diffusive
electron transport happens. When $t^{*}\gg\tau_{p}^{*}$ and
$x^{*}\gg|\bm{v}_{p}^{*}|\tau_{p}^{*}$, diffusive phonon transport happens.
When $t^{*}\gg 1/G^{*}$, $x^{*}\gg|\bm{v}_{e}^{*}|/G^{*}$ and
$x^{*}\gg|\bm{v}_{p}^{*}|/G^{*}$, electron-phonon coupling is sufficient so
that the deviations between electron and phonon temperatures tend to zero.

When both electron and phonon suffer diffusive transport processes, according
to the first-order Chapman-Enskpg expansion the distribution function can be
approximated as

$\displaystyle u_{e}^{*}$ $\displaystyle\approx
u_{e}^{eq,*}-\tau_{e}^{*}\left(\frac{\partial u_{e}^{eq,*}}{\partial
t^{*}}+|\bm{v}_{e}^{*}|\bm{s}\cdot\nabla_{\bm{x}^{*}}u_{e}^{eq,*}+\frac{(T_{e}^{*}-T_{p}^{*})/(4\pi)}{1/G^{*}}\right),$
(74) $\displaystyle u_{p}^{*}$ $\displaystyle\approx
u_{p}^{eq,*}-\tau_{p}^{*}\left(\frac{\partial u_{p}^{eq,*}}{\partial
t^{*}}+|\bm{v}_{p}^{*}|\bm{s}\cdot\nabla_{\bm{x}^{*}}u_{p}^{eq,*}-\frac{(T_{e}^{*}-T_{p}^{*})/(4\pi)}{1/G^{*}}\right).$
(75)

Combined above six equations and taking an integral of the BTE (68,69) over
the whole first Brillouin zone, we can get

$\displaystyle\frac{\partial}{\partial
t^{*}}\left<u_{e}^{eq,*}-\tau_{e}^{*}\left(\frac{\partial
u_{e}^{eq,*}}{\partial
t^{*}}+\frac{(T_{e}^{*}-T_{p}^{*})/(4\pi)}{1/G^{*}}\right)\right>+\nabla_{\bm{x}^{*}}\cdot\left<|\bm{v}_{e}^{*}||\bm{v}_{e}^{*}|\bm{s}\bm{s}\tau_{e}^{*}u_{e}^{eq,*}\right>$
$\displaystyle=-\frac{T_{e}^{*}-T_{p}^{*}}{1/G^{*}},$ (76)
$\displaystyle\frac{\partial}{\partial
t^{*}}\left<u_{p}^{eq,*}-\tau_{p}^{*}\left(\frac{\partial
u_{p}^{eq,*}}{\partial
t^{*}}+\frac{(T_{e}^{*}-T_{p}^{*})/(4\pi)}{1/G^{*}}\right)\right>+\nabla_{\bm{x}^{*}}\cdot\left<|\bm{v}_{p}^{*}||\bm{v}_{p}^{*}|\bm{s}\bm{s}\tau_{p}^{*}u_{p}^{eq,*}\right>$
$\displaystyle=\frac{T_{e}^{*}-T_{p}^{*}}{1/G^{*}},$ (77)

where $<>$ represents the integral over the whole first Brillouin zone. Then
the associated macroscopic heat conduction is

$\displaystyle C_{e}\frac{\partial T_{e}}{\partial
t}-\tau_{e}\frac{\partial^{2}U_{e}}{\partial
t^{2}}-\tau_{e}\frac{\partial(GT_{e}-GT_{p})}{\partial t}$
$\displaystyle=\nabla\cdot(\kappa_{e}\nabla T_{e})-G(T_{e}-T_{p}),$ (78)
$\displaystyle C_{p}\frac{\partial T_{p}}{\partial
t}-\tau_{p}\frac{\partial^{2}U_{p}}{\partial
t^{2}}+\tau_{p}\frac{\partial(GT_{e}-GT_{p})}{\partial t}$
$\displaystyle=\nabla\cdot(\kappa_{e}\nabla T_{p})+G(T_{e}-T_{p}),$ (79)

where $\kappa=\int C|\bm{v}|^{2}\tau/3d\bm{K}$ is the bulk thermal
conductivity. When $\sqrt{\tau_{e}^{*}}\ll t^{*}$, $\sqrt{\tau_{p}^{*}}\ll
t^{*}$, $\tau_{e}^{*}/(t^{*}/G^{*})\ll 1$, $\tau_{p}^{*}/(t^{*}/G^{*})\ll 1$,
the first order term of relaxation time in Eqs. (78,79) could be removed and
above two equations recover the typical two-temperature model in the diffusive
limit.

## References

  * Chen (2005) G. Chen,  _Nanoscale energy transport and conversion: A parallel treatment of electrons, molecules, phonons, and photons_ (Oxford University Press, 2005).
  * Kittel  _et al._ (1996) C. Kittel, P. McEuen, and P. McEuen,  _Introduction to solid state physics_ , Vol. 8 (Wiley New York, 1996).
  * Fann  _et al._ (1992a) W. S. Fann, R. Storz, H. W. K. Tom, and J. Bokor, Phys. Rev. Lett. 68, 2834 (1992a).
  * Karna  _et al._ (2023) P. Karna, M. S. B. Hoque, S. Thakur, P. E. Hopkins, and A. Giri, Nano Lett. 23, 491 (2023).
  * Block  _et al._ (2019) A. Block, M. Liebel, R. Yu, M. Spector, Y. Sivan, F. J. G. de Abajo, and N. F. van Hulst, Sci. Adv. 5, eaav8965 (2019).
  * Allen (1987) P. B. Allen, Phys. Rev. Lett. 59, 1460 (1987).
  * Caruso and Novko (2022) F. Caruso and D. Novko, Advances in Physics: X 7, 2095925 (2022).
  * Mao  _et al._ (2021) J. Mao, G. Chen, and Z. Ren, Nat. Mater. 20, 454 (2021).
  * Deng  _et al._ (2021) C. Deng, Y. Huang, M. An, and N. Yang, Mater. Today Phys. 16, 100305 (2021).
  * Pop  _et al._ (2004) E. Pop, R. W. Dutton, and K. E. Goodson, J. Appl. Phys. 96, 4998 (2004).
  * Pop (2010) E. Pop, Nano Res. 3, 147 (2010).
  * Carpene (2006) E. Carpene, Phys. Rev. B 74, 024301 (2006).
  * Schoenlein  _et al._ (1987) R. W. Schoenlein, W. Z. Lin, J. G. Fujimoto, and G. L. Eesley, Phys. Rev. Lett. 58, 1680 (1987).
  * Groeneveld  _et al._ (1992) R. H. M. Groeneveld, R. Sprik, and A. Lagendijk, Phys. Rev. B 45, 5079 (1992).
  * Waldecker  _et al._ (2016) L. Waldecker, R. Bertoni, R. Ernstorfer, and J. Vorberger, Phys. Rev. X 6, 021003 (2016).
  * Lu  _et al._ (2018) Z. Lu, A. Vallabhaneni, B. Cao, and X. Ruan, Phys. Rev. B 98, 134309 (2018).
  * Zhou  _et al._ (2015) J. Zhou, N. Li, and R. Yang,  88, 156 (2015).
  * Sun  _et al._ (1994) C.-K. Sun, F. Vallée, L. H. Acioli, E. P. Ippen, and J. G. Fujimoto, Phys. Rev. B 50, 15337 (1994).
  * Qiu and Tien (1994) T. Qiu and C. Tien, Int. J. Heat Mass Transfer 37, 2789 (1994).
  * Lin  _et al._ (2008) Z. Lin, L. V. Zhigilei, and V. Celli, Phys. Rev. B 77, 075133 (2008).
  * Rethfeld  _et al._ (2002) B. Rethfeld, A. Kaiser, M. Vicanek, and G. Simon, Phys. Rev. B 65, 214303 (2002).
  * Anisimov  _et al._ (1974) S. I. Anisimov, B. L. Kapeliovich, and T. L. Perelman, Sov. Phys.-JETP 39, 375 (1974).
  * Groeneveld  _et al._ (1995) R. H. M. Groeneveld, R. Sprik, and A. Lagendijk, Phys. Rev. B 51, 11433 (1995).
  * Brorson  _et al._ (1987) S. D. Brorson, J. G. Fujimoto, and E. P. Ippen, Phys. Rev. Lett. 59, 1962 (1987).
  * Fujimoto  _et al._ (1984) J. G. Fujimoto, J. M. Liu, E. P. Ippen, and N. Bloembergen, Phys. Rev. Lett. 53, 1837 (1984).
  * Guo and Xu (2014) L. Guo and X. Xu, Journal of Heat Transfer 136 (2014).
  * Qiu and Tien (1992) T. Qiu and C. Tien, Int. J. Heat Mass Transfer 35, 719 (1992).
  * Qiu  _et al._ (1994) T. Qiu, T. Juhasz, C. Suarez, W. Bron, and C. Tien, Int. J. Heat Mass Transfer 37, 2799 (1994).
  * Joseph and Preziosi (1989) D. D. Joseph and L. Preziosi, Rev. Mod. Phys. 61, 41 (1989).
  * Sobolev (2021) S. L. Sobolev, Nanoscale and Microscale Thermophysical Engineering 25, 153 (2021).
  * Cattaneo (1948) C. Cattaneo, Atti Sem. Mat. Fis. Univ. Modena 3, 83 (1948).
  * Pattamatta and Madnia (2009a) A. Pattamatta and C. K. Madnia, Numerical Heat Transfer, Part A: Applications 55, 611 (2009a).
  * Chen  _et al._ (2006) J. Chen, D. Tzou, and J. Beraun, Int. J. Heat Mass Transfer 49, 307 (2006).
  * Pattamatta and Madnia (2009b) A. Pattamatta and C. K. Madnia, J. Heat Transf. 131 (2009b).
  * Miao and Wang (2021) W. Miao and M. Wang, Phys. Rev. B 103, 125412 (2021).
  * Tong and Bernardi (2021) X. Tong and M. Bernardi, Phys. Rev. Research 3, 023072 (2021).
  * Wang  _et al._ (2016) Y. Wang, Z. Lu, A. K. Roy, and X. Ruan, Journal of Applied Physics 119, 065103 (2016).
  * Jacoboni and Reggiani (1983) C. Jacoboni and L. Reggiani, Rev. Mod. Phys. 55, 645 (1983).
  * Muthukunnil Joseph and Cao (2022) A. Muthukunnil Joseph and B.-Y. Cao, Int. J. Therm. Sci 181, 107742 (2022).
  * Adams and Larsen (2002) M. L. Adams and E. W. Larsen, Prog. Nucl. Energ. 40, 3 (2002).
  * Guo and Xu (2021) Z. Guo and K. Xu, Adva. Aerodyn. 3, 6 (2021).
  * Zhang and Guo (2019) C. Zhang and Z. Guo, Int. J. Heat Mass Transfer 134, 1127 (2019).
  * Guo  _et al._ (2023) Z. Guo, J. Li, and K. Xu, Phys. Rev. E 107, 025301 (2023).
  * Qiu and Tien (1993) T. Q. Qiu and C. L. Tien, J. Heat Transf. 115, 835 (1993).
  * Hopkins  _et al._ (2011) P. E. Hopkins, L. M. Phinney, and J. R. Serrano, Journal of Heat Transfer 133 (2011).
  * Fann  _et al._ (1992b) W. S. Fann, R. Storz, H. W. K. Tom, and J. Bokor, Phys. Rev. B 46, 13592 (1992b).
  * Sun  _et al._ (1993) C.-K. Sun, F. Vallée, L. Acioli, E. P. Ippen, and J. G. Fujimoto, Phys. Rev. B 48, 12365 (1993).
  * Hohlfeld  _et al._ (2000) J. Hohlfeld, S.-S. Wellershoff, J. Güdde, U. Conrad, V. Jähnke, and E. Matthias, Chemical Physics 251, 237 (2000).
  * Anisimov and Rethfeld (1997) S. I. Anisimov and B. Rethfeld, in  _Nonresonant Laser-Matter Interaction (NLMI-9)_ , Vol. 3093 (SPIE, 1997) pp. 192–203.
  * Mozafarifard  _et al._ (2023) M. Mozafarifard, Y. Liao, Q. Nian, and Y. Wang, Int. J. Heat Mass Transfer 202, 123759 (2023).
  * Sivan and Spector (2020) Y. Sivan and M. Spector, ACS Photonics 7, 1271 (2020).
  * Ibrahim  _et al._ (2004) W. M. Ibrahim, H. E. Elsayed-Ali, C. E. Bonner, and M. Shinn, Int. J. Heat Mass Transfer 47, 2261 (2004).
  * Freedman  _et al._ (2019) J. P. Freedman, R. F. Davis, and J. A. Malen, Phys. Rev. B 99, 054308 (2019).
  * Block  _et al._ (2023) A. Block, R. Yu, I.-W. Un, S. Varghese, M. Liebel, N. F. van Hulst, S. Fan, K.-J. Tielrooij, and Y. Sivan, ACS Photonics 10, 1150 (2023).
  * Maznev  _et al._ (2011) A. A. Maznev, J. A. Johnson, and K. A. Nelson, Journal of Applied Physics 109 (2011).
  * Giri  _et al._ (2015) A. Giri, J. T. Gaskins, B. F. Donovan, C. Szwejkowski, R. J. Warzoha, M. A. Rodriguez, J. Ihlefeld, and P. E. Hopkins, Journal of Applied Physics 117, 105105 (2015).
  * Hopkins  _et al._ (2009) P. E. Hopkins, J. L. Kassebaum, and P. M. Norris, Journal of Applied Physics 105 (2009), 10.1063/1.3068476.
  * Choi  _et al._ (2014) G.-M. Choi, R. B. Wilson, and D. G. Cahill, Phys. Rev. B 89, 064307 (2014).
  * Chen  _et al._ (2022) J. Chen, X. Xu, J. Zhou, and B. Li, Rev. Mod. Phys. 94, 025002 (2022).
  * Miao and Wang (2023) W. Miao and M. Wang, Int. J. Heat Mass Transfer 200, 123538 (2023).
  * SINGH (2010) N. SINGH, International Journal of Modern Physics B 24, 1141 (2010).
  * Guo  _et al._ (2012) L. Guo, S. L. Hodson, T. S. Fisher, and X. Xu, Journal of Heat Transfer 134 (2012).