\section{Effective conductivity problem}
We start by recalling the notion of effective conductivity in the sense of stochastic homogenization theory for an heterogeneous material made of inclusions in a given matrix of homogeneous conductivity. Let $d\ge1$ denote the space dimension.

\subsection{Stochastic setting}
We shall use a statistical description for the set of inclusions in the material. Restricting to spherical inclusions for notational simplicity, we let
\begin{equation}\label{eq:inclusion}
\Bc(\Pc)\,:=\,\cup_{x\in\Pc}B(x),
\end{equation}
where $B(x)=B+x$ stands for the unit ball centered at $x$ in $\R^d$, and where $\Pc$ is the set of centers of the inclusions. It remains to define statistical ensembles for the latter.
We call {\it point set} any countable subset $\Pc\subset\R^d$ that is {locally finite} in the sense that for any bounded $E\subset\R^d$ the number of points of~$\Pc$ in~$E$ is finite, $\Pc(E):=\sharp\{\Pc \cap E\}<\infty$.
A point set~$\Pc$ can be represented by the associated locally finite measure $\sum_{x \in \Pc} \delta_x$, which acts on the space of compactly supported continuous functions via $f \mapsto\Pc(f):=\sum_{x\in \Pc} f(x)$. We endow the space $\Omega$ of point sets with the smallest $\sigma$-algebra that makes all evaluation maps $\Pc\mapsto f(\Pc)$ measurable.
A {\it random point process} is then defined as a probability measure $\Pm$ on $\Omega$, and we denote by $\expec{\cdot}$ the associated expectation.
We further define stationarity and ergodicity with respect to translations $\Pc+z:=\{x+z:x\in\Pc\}$ of point sets: The point process is said to be {\it stationary} (or {\it statistically translation-invariant}) if for any measurable set $A\subset\Omega$ we have $\pr{A+z}=\pr{A}$ for all $z\in\R^d$, where we use the notation $A+z:=\{\Pc+z:\Pc\in A\}$. The process is said to be {\it ergodic} if any measurable set $A\subset\Omega$ that is translation-invariant, in the sense that $\pr{A\setminus(A+z)}=\pr{(A+z)\setminus A}=0$ for all $z\in\R^d$, satisfies~$\pr{A}=0$ or~$1$.
In the sequel, the set of inclusions in the material is modeled by~\eqref{eq:inclusion} with $\Pc$ sampled according to some stationary and ergodic random point process. Note that the inclusions are allowed to overlap in general.

\subsection{Effective conductivity}
Given a stationary ergodic random point process as defined above, we consider the associated coefficient field
\begin{equation}\label{e.coeff}
A(y):=A_1 +\mathds{1}_{\Bc(\Pc)}(y)(A_2-A_1) ,\qquad y\in\R^d,
\end{equation}
where $A_1,A_2\in\R^{d\times d}$ are two strongly elliptic matrices and where $\Pc$ is sampled according to the point process.
This models a homogeneous material of conductivity $A_1$ that is perturbed by disordered spherical inclusions of another material of conductivity $A_2$.
On large scales, in the sense of the homogenization theory, this two-phase heterogeneous material behaves like a homogeneous material with some effective conductivity $\bar A$ defined  by
\begin{equation}\label{0-1}
\bar A e = \expec{A(\nabla \phi_e+e)},
\end{equation}
where $\phi_e$ is the so-called corrector defined as the unique weak solution of the whole-space equation
\begin{equation}\label{0-2}
-\nabla \cdot A (\nabla \phi_e+e)=0,\qquad \text{in $\R^d$},
\end{equation}
in the following class: $\phi_e$ is almost surely in $H^1_\loc(\R^d)$, satisfies the anchoring condition~$\fint_{B(0)} \phi_e=0$, and its gradient $\nabla \phi_e$ is a stationary random field with vanishing expectation $\expec{\nabla \phi_e}=0$, and finite second moments~$\expec{|\nabla \phi_e|^2}<\infty$; see \cite{MR542557,MR712714,MR1329546}.
Note that $\bar A$ is not explicit in general.

\section{Dilute homogenization}

\subsection{The Clausius--Mossotti formula}
Homogenization was neither born in the mathematical community in the 1970s, nor in the engineering community the decade before: it emerged much earlier, in the second half of the 19th century, in the physics community, in the context of two-phase dispersed media.
Motivated by the works of Poisson~\cite{Poisson} and Faraday~\cite{Faraday}, Mossotti and Clausius were the first to investigate the question of the effective dielectric constant of a homogeneous background material perturbed by sparse spherical inclusions~\cite{Mossotti-36,Mossotti-50,Clausius-79}.
The problem was largely
revisited by Maxwell~\cite{Maxwell-81} for the effective conductivity of two-phase media; we refer to~\cite{Markov-00} for a detailed account of the historical context.
In modern language, these authors argued that the effective conductivity~\eqref{0-1} associated with the two-phase model~\eqref{e.coeff} takes on the following guise, in dimension $d\ge1$, in case of isotropic conductivities $A_1=\alpha\Id$ and $A_2=\beta\Id$,
\begin{equation}\label{0-10}
\bar A=\alpha \Id +\varphi \frac{\alpha d(\beta-\alpha)}{\beta+\alpha(d-1)} \Id+o(\varphi),\qquad\text{as $\varphi\downarrow0$},
\end{equation}
where $\varphi$ stands for the volume fraction of the inclusions, that is, by the ergodic theorem,
\begin{equation}\label{eq:def-varphi}
\varphi\,:=\,\expecm{\mathds1_{\Bc(\Pc)}}\,=\,\lim_{R\uparrow\infty}\tfrac{|\Bc(\Pc)\cap RB|}{|RB|},\qquad\text{for $\Pm$-almost all $\Pc$}.
\end{equation}
In case of disjoint inclusions, note that $\varphi=\lambda|B|$ where $\lambda$ is the intensity of the point process (defined by $\lambda|E|=\expec{\Pc(E)}$ for any Borel set $E\subset\R^d$).
This dilute approximation~\eqref{0-10} for the effective medium is known as the Clausius--Mossotti formula and was soon adapted to various other physical settings, in particular by Lorenz and Lorentz for the effective refractive index of two-phase media in optics~\cite{Lorenz-80,Lorentz-09}, and by Einstein for the effective viscosity of a Stokes fluid with a dilute suspension of rigid particles~\cite{Einstein-05,Einstein-11}. Einstein's result was actually part of his PhD thesis,
where he used it to design a celebrated experiment to measure the Avogadro number; see e.g.~the inspiring historical account in~\cite{Straumann-05}.

We start by describing a heuristic argument for~\eqref{0-10}.
Recalling that~$\expec{\nabla \phi_e}=0$,
the effective conductivity~\eqref{0-1} can be decomposed as
\begin{eqnarray*}
e \cdot \bar A e &=& \expec{e\cdot A(\nabla \phi_e+e)}\\
&=&\expec{e\cdot \big(\alpha+(\beta-\alpha)\mathds 1_{\Bc(\Pc)}\big) (\nabla \phi_e+e)}\\
&=& \alpha|e|^2+ (\beta-\alpha)\,\expec{\mathds 1_{\Bc(\Pc)}\, e\cdot (\nabla \phi_e+e)},
\end{eqnarray*}
and thus, by the ergodic theorem, assuming for simplicity that inclusions are almost surely disjoint,
\begin{equation}\label{0-12}
e \cdot \bar A e
\,=\,\alpha|e|^2+ (\beta-\alpha)\,\lim_{R\uparrow \infty} \tfrac1{|RB|} \sum_{x \in \Pc \cap RB}e \cdot\int_{B(x)} (  \nabla \phi_e+e).
\end{equation}
In the dilute regime~\mbox{$\varphi\ll 1$}, inclusions are typically far from one another and therefore do not `interact' much when solving the corrector equation~\eqref{0-2}.
More precisely, for all~$x\in\Pc$, we may heuristically approximate the corrector in the inclusion $B(x)$ by the solution of a corresponding single-inclusion problem,
\begin{equation}\label{0-13}
\nabla \phi_e|_{B(x)} \simeq \nabla \psi_e(\cdot-x),
\end{equation}
where $\psi_e$ is the unique weak solution in $\dot H^1(\R^d)$ of the whole-space single-inclusion equation
\begin{equation*}
-\nabla \cdot \big(\alpha+ (\beta-\alpha)  \mathds1_{B(0)}\big) (\nabla \psi_e+e)=0,\qquad\text{in $\R^d$}.
\end{equation*}
In the present case of spherical inclusions and isotropic conductivity, this equation is explicitly solvable in form of
\begin{equation}\label{0-11}
\nabla \psi_e(x)=
\begin{cases}
Ke,&\text{for $|x|<1$};\\
\frac{K}{|x|^d}\left(e-d\frac{x\cdot e}{|x|}\frac{x}{|x|}\right)
,&\text{for $|x|>1$};
\end{cases}
\end{equation}
with
\begin{equation*}
K\,=\,\frac{\alpha-\beta}{\beta+\alpha(d-1)}.
\end{equation*}
Inserting this form into~\eqref{0-13} and~\eqref{0-12}, the Clausius--Mossotti formula~\eqref{0-10} heuristically follows.

\subsection{Main result}
The aim of the present note is to prove the Clausius--Mossotti formula~\eqref{0-10} in the most general setting possible and to establish a sharp error bound.
The above heuristic argument indicates that the error mostly comes from `interactions' between inclusions, as it amounts to locally replacing the corrector by solutions of single-inclusion problems, cf.~\eqref{0-13}.
To quantify this error, we need to recall the notion of {\it second-order intensity} $\lambda_2$ of the point process, which we introduced in~\cite{2008.03837}. For that purpose, we first define the minimal lengthscale $\ell\ge0$ of the point process,
\begin{equation}\label{0-14}
\ell\,:=\, \inf_{x,y\in\Pc\atop x\ne y} |x-y|_\infty,\qquad\text{for $\Pm$-almost all $\Pc$,}
\end{equation}
which is deterministic by ergodicity.
In case $\ell>0$ (that is, if the point process is hardcore), the {\it second-order intensity} is defined as
\begin{equation}\label{0-15}
\qquad\lambda_2\,:=\,\sup_{z_1,z_2\in\R^d}\,\expecM{\sum_{x_1,x_2\in\Pc\atop x_1\ne x_2}\ell^{-d}\,\mathds1_{Q_{\ell}(z_1)}(x_1)\,\ell^{-d}\,\mathds1_{Q_{\ell}(z_2)}(x_2)},
\end{equation}
where $Q_r(z):=r Q+z$ is the cube of sidelength $r$ centered at $z$. Note that, by definition~\eqref{0-14}, each cube $Q_\ell(z)$ contains at most one point of $\Pc$ almost surely.
In other words, $\lambda_2$ is the maximum expected number of couples of points that lie in the $\ell$-neighborhood of a given element of $(\R^d)^2$, properly normalized by $\ell$.
Alternatively, recalling that the {\it$2$-point density} is the non-negative function~$f_{2}$ defined
by the following relation,
\begin{equation*}
\qquad\expecM{\sum_{x_1, x_2\in\Pc\atop x_1\ne x_2}\zeta(x_{1},x_{2})} \,=\, \int_{(\R^d)^2}\zeta f_2,\qquad\text{for all $\zeta \in C_c((\R^d)^2)$},
\end{equation*}
the definition~\eqref{0-15} of $2$-point intensity can be reformulated as
\begin{equation}\label{eq:high-intens/Re}
\lambda_2\,=\,\sup_{z_1,z_2\in\R^d}\,\fint_{Q_{\ell}(z_1)\times Q_\ell(z_2)}f_2.
\end{equation}
In case $\ell=0$ (that is, if the point process is not hardcore), this definition is naturally extended to $\lambda_2:=\|f_2\|_{\Ld^\infty((\R^d)^2)}$. For a Poisson point process, due to the tensor structure, the $2$-point intensity is simply the square of the intensity, $\lambda_2=\lambda^2$, but for a general mixing point process it can be anything in the interval $[\lambda^2,\lambda]$ and its smallness describes some form of local independence.

We can now state the main result of this note, which is an adaptation of our recent work on Einstein's formula~\cite[Theorem~1]{2008.03837} to the effective conductivity problem.
It proves the validity of the Clausius--Mossotti formula~\eqref{0-10} for the first time in the setting of general point processes.
The error bound $O(\lambda_2|\!\log\lambda_2|)$ below is new and is sharp in general, cf.~\cite[Theorem~7]{2008.03837}, showing that the approximation~\eqref{0-10} is only valid in the dilute regime if~$\lambda_2|\!\log\lambda_2|=o(\varphi)$.
In the specific case of a Poisson point process, this result follows from~\cite{CRMATH_2022__360_G8_909_0} with the improved error bound $O(\lambda_2)=O(\lambda^2)$ (without logarithmic correction). The particular case of dilute point processes obtained by Bernoulli deletion or by dilation of a given process was already treated in~\cite{MR3458165,Pertinand}.
Note that inclusions here are allowed to overlap and that no upper bound is assumed on the number of points per unit volume.

\begin{theor}\label{th:CM-DG}
Given a stationary and ergodic point process, and given $\alpha,\beta>0$,
the effective conductivity~\eqref{0-1} associated with the two-phase model~\eqref{e.coeff} with isotropic conductivities $A_1=\alpha\Id$ and $A_2=\beta\Id$ satisfies the following quantitative version of~\eqref{0-10},
\begin{equation*}
\Big|\bar A-\Big(\alpha \Id + \varphi\frac{\alpha d(\beta-\alpha)}{\beta+\alpha(d-1)} \Id \Big)\Big| \,\lesssim\,\lambda_2|\!\log\lambda_2|.
\qedhere
\end{equation*}
\end{theor}

We shall prove this result in the following slightly more general form, where conductivities $A_1$ and $A_2$ are no longer assumed to be isotropic and where $A_2$ may itself be heterogeneous.
We consider spherical inclusions for notational convenience, but we emphasize that, as in~\cite{MR3458165,2008.03837}, this is not essential (only the above explicit form of the Clausius--Mossotti formula then needs to be changed).

\begin{theor} \label{th:CM-DG+}
Let $A_1\in\R^{d\times d}$ be a strictly elliptic (non-necessarily symmetric) matrix, and let $A_2$ be a
stationary and ergodic random field of uniformly elliptic (non-necessarily symmetric) matrices.
Given a stationary and ergodic point process that is independent of~$A_2$, consider the associated coefficient field
\begin{equation}\label{eq:model-A}
A(y):=A_1+\mathds{1}_{\Bc(\Pc)}(y)(A_2(y)-A_1),\qquad y\in\R^d,
\end{equation}
where $\Pc$ is sampled according to the point process.
Then, the effective coefficient $\bar A$ associated via~\eqref{0-1} satisfies the following expansion,
\begin{align} \label{0-26+}
&\big|\bar A-(A_1 +\varphi\widehat A_2 )\big| \,\lesssim\,\lambda_2 |\!\log\lambda_2|,
\end{align}
where the first-order effective correction $\widehat A_2$ is given by
\begin{equation}\label{eq:def-corrA2}
\widehat A_2 e \,:=\,\expecM{\fint_{B} (A_2-A_1)(\nabla \psi_e+e)},
\end{equation}
where $\psi_e$ is the unique weak solution in $\dot H^1(\R^d)$ of the whole-space single-inclusion problem
\begin{equation}\label{eq:single-incl-1}
-\nabla \cdot \big(A_1+\mathds 1_{B}(A_2-A_1)\big)(\nabla \psi_e+e)=0,\qquad\text{in $\R^d$.}\qedhere
\end{equation}
\end{theor}

\subsection{Previous contributions}

The asymptotic analysis of the effective conductivity in case of a periodic array of inclusions with a small volume fraction was first addressed by Berdichevski{\u\i}~\cite{Berdichevski-83}; see also~\cite[Section~1.7]{MR1329546}.
The first justification of the  Clausius--Mossotti formula in a random setting is due to Almog in dimension $d=3$, whose results in~\cite{Almog-13,Almog-14} precisely yield~\eqref{0-10} when combined with elementary homogenization theory.
The proof is based on (scalar) potential theory and crucially relies on the facts that the space dimension is $d=3$, that~$A$ is everywhere isotropic, and that the inclusions are spherical and disjoint.
Another contribution is due to Mourrat~\cite{Mourrat-13}, who studied for all~$d\ge 2$ a discrete elliptic equation (instead of a continuum one) with sparse i.i.d.\@ perturbations of the conductivity, proving~\eqref{0-10} in that setting by strongly relying on quantitative stochastic homogenization results of~\cite{MR2789576,MR2932541}. We also highlight the inspiring work~\cite{Anantharaman-LeBris-11,Anantharaman-LeBris-12} by Anantharaman and Le Bris, who obtained related results on sparse i.i.d.\@ perturbations of a periodic array of inclusions; see also~\cite{Anantharaman-these}.
Those different previous results were largely surpassed in~\cite{MR3458165}, where we studied the case of a general stationary and ergodic inclusion process, focusing on a dilute regime obtained by a Bernoulli deletion procedure where each inclusion is preserved independently with low probability, and where we established real analyticity with respect to the Bernoulli parameter --- only assuming that the number of inclusions per unit volume be uniformly bounded. This was recently extended in~\cite{doi:10.1142/S0219199722500274,CRMATH_2022__360_G8_909_0} to prove Gevrey regularity when starting from a Poisson point process (for which the uniform boundedness assumption fails).
In a different vein, partly inspired by Berdichevski{\u\i}'s approach in the periodic setting, Pertinand addressed in~\cite{Pertinand} the case when the dilution of the point process is obtained by dilating a given hardcore process, and he proved the real analyticity with respect to the inverse of the dilation parameter.
Very recently, G\'erard-Varet~\cite{MR4400909} proposed an alternative approach where
he bypasses homogenization theory and directly quantifies in terms of the volume fraction $\varphi$ the distance between the solution of a problem with sparse inclusions and that of an effective problem with conductivity given by the Clausius--Mossotti formula.
Various related contributions concern the validity of Einstein's formula for the effective viscosity of a Stokes fluid with a dilute suspension of rigid particles; see~\cite{MR2982744,MR4098775,MR4102716,MR4260456,2008.03837}.
Our approach in the present note is an adaptation of the recent short proof that we obtained in~\cite{2008.03837} for Einstein's formula: it allows to justify the Clausius--Mossotti formula for the first time in the setting of general point processes and to determine the optimal error estimate. Note that several points of the proof simplify in the present setting. In particular we manage to fully bypass the variational formulation of~\cite{2008.03837}.

While the Clausius--Mossotti formula is universal in the sense that it only depends on the set of inclusions via its volume fraction (and on the shape of the inclusions, here assumed to be spherical),
the next-order correction further depends on the two-point correlation function. The identification of this correction was first discussed in~\cite{MR3458165}, and it has been the object of many recent contributions in the context of Einstein's formula~\cite{BG-72,GVH,GV-20,GVM-20,2008.03837}; we refer in particular to our recent work~\cite{2008.03837} where all higher-order corrections are systematically described in form of a cluster expansion.

\section{Proof of Theorem~\ref{th:CM-DG+}}
We denote by $\Pc=\{x_n\}_n$ the point set sampled according to the underlying random point process.
In order to justify the approximation~\eqref{0-13} of the corrector in terms of single-inclusion problems, we start by singling out clusters of close inclusions:
let $\{K_{q}\}_q$ be the family of connected components of the fattened set $\Bc(\Pc)+ B$, and consider the corresponding clusters
\[\qquad J_{q}\,:=\,\bigcup_{n:B(x_n)\subset K_{q}}B(x_n).\]
As we focus here on spherical inclusions, we note that a point $x_n\in\Pc$ does not belong to a cluster $J_q$ made of at least two inclusions if and only if
\[\rho_n\,:=\,\tfrac12\inf_{m:m\ne n}|x_n-x_m| \,\ge\,2.\]
Let then $\Sc:=\{n:\rho_n \ge2\}$ be the set of indices corresponding to well-separated inclusions
and let~$\Sc'$ be the set of indices $q$ such that the cluster $J_{q}$ is made of at least two inclusions. (In what follows,~$\Sc$ or $\Sc'$ can be empty.)
Next, in order to define suitable neighborhoods of the inclusions, we consider the Voronoi tessellation~$\{V_n\}_{n}$ associated with the point set~$\Pc=\{x_n\}_n$, that is,
\begin{equation}\label{eq:Voronoi}
V_n\,:=\,\Big\{z\in\R^d:|z-x_n|<\inf_{m:m\ne n}|z-x_m|\Big\}.
\end{equation}
We merge Voronoi cells that correspond to inclusions in the same cluster, setting
\[W_q\,:=\,\bigcup_{n:B(x_n)\subset J_{q}}V_n,\]
and we then partition the whole space as
\[\R^d\,=\,\Big(\bigcup_{n\in\Sc}V_n\Big)\cup\Big(\bigcup_{q\in\Sc'}W_q\Big).\]
We shall repeatedly use the following elementary property of Voronoi tessellations, see~\cite[Lemma~2.5]{2008.03837}: for all stationary random fields $\zeta$ with $\expec{|\zeta|}<\infty$, we have
\begin{equation}\label{eq:Voronoi}
\expec{\zeta}\,=\,\expecM{\sum_{n\in\Sc}\frac{\mathds1_{0\in B(x_n)}}{|B|}\int_{V_n}\zeta}+\expecM{\sum_{q\in\Sc'}\frac{\mathds1_{0\in J_q}}{|J_q|}\int_{W_q}\zeta}.
\end{equation}
We shall also use the short-hand notation $\tilde\Bc(\Pc):=\cup_{n\in \Sc} B(x_n)$ for the union of well-separated inclusions.
We split the proof into six steps. Let $e\in\R^d$ with $|e|=1$.

\medskip
\step1 Representation formula for the error: proof that
\begin{equation}\label{eq:error-estim0}
\big|\bar Ae-(A_1e+\varphi\widehat A_2e)\big|
\,\lesssim\,
\expecM{\sum_{n\in\Sc}\mathds1_{B(x_n)}|\nabla(\phi_e-\psi_{e,n})|}
+\expecm{\mathds1_{\Bc(\Pc)\setminus\tilde\Bc(\Pc)}\big(1+|\nabla \phi_e|\big)},
\end{equation}
where $\widehat A_2$ is given by~\eqref{eq:def-corrA2} in the statement, and where for all $n$ we define $\psi_{e,n}$ as the unique weak solution in $\dot H^1(\R^d)$ of the whole-space single-inclusion problem centered at~$x_n$, that is,
\begin{equation*}
-\nabla \cdot \big(A_1+\mathds1_{B(x_n)}(A_2-A_1)\big)(\nabla \psi_{e,n}+e)=0,\qquad\text{in $\R^d$}.
\end{equation*}
Recall that as in the statement we also denote by $\psi_e$ the solution of the corresponding single-inclusion problem with center~$x_n$ replaced by the origin~$0$, cf.~\eqref{eq:single-incl-1}.

\medskip\noindent
By definition of the two-phase coefficient field $A$, cf.~\eqref{eq:model-A}, and of the associated effective conductivity~$\bar A$, cf.~\eqref{0-1},  we find
\begin{eqnarray*}
\bar A e& =& \expec{A(\nabla \phi_e+e)}\\
&=&\expec{A_1 (\nabla \phi_e+e)}+\expec{\mathds1_{\Bc(\Pc)}(A_2-A_1) (\nabla \phi_e+e)}\\
&=& A_1e+\expec{\mathds1_{\Bc(\Pc)}(A_2-A_1) (\nabla \phi_e+e)},
\end{eqnarray*}
where the last identity follows from the fact that $A_1$ is constant and $\expec{\nabla\phi_e}=0$.
Decomposing $\Bc(\Pc)$ into well-separated inclusions and into clusters of close inclusions, and comparing the corrector $\phi_e$ to single-inclusion solutions $\{\psi_{e,n}\}_n$ in well-separated inclusions, we can decompose
\begin{multline}\label{eq:pre-decomp-barA}
\bar A e
\,=\, A_1e+\expecM{\sum_{n\in\Sc}\mathds1_{B(x_n)}(A_2-A_1) (\nabla\psi_{e,n}+e)}\\
+\expecM{\sum_{n\in\Sc}\mathds1_{B(x_n)}(A_2-A_1) \nabla(\phi_e-\psi_{e,n})}
+\expec{\mathds1_{\Bc(\Pc)\setminus\tilde\Bc(\Pc)}(A_2-A_1) (\nabla \phi_e+e)}.
\end{multline}
In order to reformulate the first right-hand side term, we appeal to~\eqref{eq:Voronoi}, to the effect of
\begin{equation*}
\expecM{\sum_{n\in\Sc}\mathds1_{B(x_n)}(A_2-A_1) (\nabla\psi_{e,n}+e)}
\,=\,\expecM{\sum_{n\in\Sc}\frac{\mathds1_{0\in B(x_n)}}{|B|}\int_{B(x_n)}(A_2-A_1) (\nabla\psi_{e,n}+e)},
\end{equation*}
which can be further rewritten as follows, recalling that $A_2$ is stationary and independent of the point process,
\begin{eqnarray*}
\expecM{\sum_{n\in\Sc}\mathds1_{B(x_n)}(A_2-A_1) (\nabla\psi_{e,n}+e)}
&=&\expecM{\sum_{n\in\Sc}\frac{\mathds1_{0\in B(x_n)}}{|B|}}\expecM{\int_{B}(A_2-A_1) (\nabla\psi_{e}+e)}\\
&=&\expecm{\mathds1_{\tilde\Bc(\Pc)}}\expecM{\fint_{B}(A_2-A_1) (\nabla\psi_{e}+e)}.
\end{eqnarray*}
Recognizing the definition~\eqref{eq:def-corrA2} of~$\widehat A_2$, noting that the definition~\eqref{eq:def-varphi} of the volume fraction yields
\[\expecm{\mathds1_{\tilde\Bc(\Pc)}}\,=\,\varphi-\expecm{\mathds1_{\Bc(\Pc)\setminus\tilde\Bc(\Pc)}},\]
and using the energy estimate $\int_{\R^d}|\nabla\psi_e|^2\lesssim1$ for the solution of~\eqref{eq:single-incl-1},
we deduce
\begin{eqnarray*}
\bigg|\expecM{\sum_{n\in\Sc}\mathds1_{B(x_n)}(A_2-A_1) (\nabla\psi_{e,n}+e)}-\varphi\widehat A_2e\bigg|
\,\lesssim\,\expecm{\mathds1_{\Bc(\Pc)\setminus\tilde\Bc(\Pc)}}.
\end{eqnarray*}
Combined with~\eqref{eq:pre-decomp-barA}, this yields the claim~\eqref{eq:error-estim0}.

\medskip
\step2 Approximation of the corrector by local Dirichlet problems.\\
Instead of directly comparing the corrector $\phi_e$ to whole-space single-inclusion solutions $\{\psi_{e,n}\}_n$ as in~\eqref{eq:error-estim0}, we start by comparing to solutions of single-inclusion Dirichlet problems in Voronoi cells.
More precisely, for all $n\in\Sc$, we define $\psi_{e,n}^\circ$ as the unique weak solution in $H^1_0(V_n)$ of the single-inclusion problem
\begin{equation}\label{eq:psi0n}
-\nabla \cdot A (\nabla \psi_{e,n}^\circ+e)=0,\qquad\text{in $V_n$},
\end{equation}
and, for all $q\in\Sc'$, we define $\psi_{e,q}^{\prime\circ}$ as the unique weak solution in $H^1_0(W_q)$ of
\[-\nabla \cdot A (\nabla \psi_{e,q}^{\prime\circ}+e)=0,\qquad\text{in $W_q$}.\]
Implicitly extending $\psi_{e,n}^\circ$ by zero outside $V_n$, and $\psi^{\prime\circ}_{e,q}$ by zero outside $W_q$, we then set
\[\psi_e^\circ\,:=\,\sum_{n\in \Sc} \psi_{e,n}^\circ+\sum_{q \in W_q} \psi_{e,q}^{\prime\circ}.\]
By definition, $\nabla \psi_e^\circ$ is stationary and we claim that it has vanishing expectation and finite second moments,
\begin{equation}\label{eq:prop-nab-psi0}
\expec{\nabla\psi_e^\circ}=0,\qquad\expec{|\nabla \psi_e^\circ|^2} \lesssim 1.
\end{equation}
Indeed, for all $R>0$, applying~\eqref{eq:Voronoi} to $|\nabla\psi_e^\circ|^2\wedge R$, we find
\begin{equation*}
\expec{|\nabla \psi_e^\circ|^2\wedge R}
\,\le\,\expecM{\sum_{n\in \Sc} \frac{\mathds{1}_{0\in B(x_n)}}{|B|} \int_{V_n} |\nabla \psi_e^\circ|^2}+ \expecM{\sum_{q\in \Sc'} \frac{\mathds{1}_{0\in J_{q}}}{|J_{q}|} \int_{W_q} |\nabla \psi_e^\circ|^2},
\end{equation*}
and thus, using energy estimates for $\psi_e^\circ$ in Voronoi cells, and further applying~\eqref{eq:Voronoi} to the constant function $1$,
\begin{equation*}
\expec{|\nabla \psi_e^\circ|^2\wedge R}
\,\lesssim\,\expecM{\sum_{n\in \Sc} \frac{\mathds{1}_{0\in B(x_n)}}{|B|} |V_n|}+ \expecM{\sum_{q\in \Sc'} \frac{\mathds{1}_{0\in J_{q}}}{|J_{q}|} |W_q|}
\,=\,1.
\end{equation*}
By the monotone convergence theorem, this proves the claim $\expec{|\nabla\psi_e^\circ|^2}\lesssim1$.
Next, we can apply~\eqref{eq:Voronoi} to $\nabla\psi_e^\circ$, to the effect of
\begin{eqnarray*}
\expec{\nabla \psi_e^\circ}
&=& \expecM{\sum_{n\in \Sc} \frac{\mathds{1}_{0\in B(x_n)}}{|B|} \int_{V_n} \nabla \psi_e^\circ}+ \expecM{\sum_{q\in \Sc'} \frac{\mathds{1}_{0\in J_{q}}}{|J_{q}|} \int_{W_q} \nabla \psi_e^\circ},
\end{eqnarray*}
where the right-hand side vanishes due to homogeneous Dirichlet boundary conditions, and the claim~\eqref{eq:prop-nab-psi0} follows.

\medskip
\step3 Approximation error estimate: proof that
\begin{equation}\label{0-20}
\expec{ |\nabla (\phi_e-\psi_e^\circ)|^2}+\expecm{ \mathds1_{\Bc(\Pc)} |\nabla (\phi_e-\psi_e^\circ)|}
\,\lesssim\,\expecM{\sum_{n\in\Sc}\frac{\mathds1_{0\in B(x_n)}}{|B|}\rho_n^{-d}}+\expecm{\mathds1_{\Bc(\Pc)\setminus\tilde\Bc(\Pc)}}.
\end{equation}
We start by proving the estimate on $\expec{ |\nabla (\phi_e-\psi_e^\circ)|^2}$.
Using the corrector equation~\eqref{0-2} for~$\phi_e$ in form of
\[\expec{ \nabla (\phi_e-\psi_e^\circ)\cdot A(\nabla\phi_e+e)}=0,\]
we find
\begin{eqnarray*}
\expec{ |\nabla (\phi_e-\psi_e^\circ)|^2}
&\lesssim&\expec{ \nabla (\phi_e-\psi_e^\circ)\cdot A\nabla (\phi_e-\psi_e^\circ)}\\
&=&-\expec{ \nabla (\phi_e-\psi_e^\circ)\cdot A(\nabla\psi_e^\circ+e)},
\end{eqnarray*}
which we can further decompose into
\begin{multline*}
\expec{ |\nabla (\phi_e-\psi_e^\circ)|^2}
\,\lesssim\,
-\expecm{\nabla(\phi_e-\psi_e^\circ)\cdot Ae}\\
+\expecm{(\nabla\psi_e^\circ+e)\cdot A\nabla\psi_e^\circ}
-\expecm{(\nabla\phi_e+e)\cdot A\nabla\psi_e^\circ}.
\end{multline*}
As $A_1$ is constant and as $\expec{\nabla(\phi_e-\psi_e^\circ)}=0$, the first right-hand side term is equal to
\begin{equation*}
\expecm{\nabla(\phi_e-\psi_e^\circ)\cdot Ae}\,=\,\expecm{\nabla(\phi_e-\psi_e^\circ)\cdot \mathds1_{\Bc(\Pc)}(A_2-A_1)e},
\end{equation*}
hence
\begin{multline}\label{eq:pre-decomp-phi-psi0}
\expec{ |\nabla (\phi_e-\psi_e^\circ)|^2}
\,\lesssim\,
-\expecm{\nabla(\phi_e-\psi_e^\circ)\cdot \mathds1_{\Bc(\Pc)}(A_2-A_1)e}\\
+\expecm{(\nabla\psi_e^\circ+e)\cdot A\nabla\psi_e^\circ}
-\expecm{(\nabla\phi_e+e)\cdot A\nabla\psi_e^\circ}.
\end{multline}
If $A$ was symmetric, then the corrector equation would ensure that the last right-hand side term vanishes. Using that $A_1$ is constant, we shall show that, even though this term does not vanish in the general non-symmetric case, it can be localized inside inclusions.
For that purpose, we rewrite the corrector equation as
\[-\nabla\cdot A_1\nabla\phi_e\,=\,\nabla\cdot Ae+\nabla\cdot \mathds1_{\Bc(\Pc)}(A_2-A_1)\nabla\phi_e,\qquad\text{in $\R^d$},\]
and we note that, as the coefficient $A_1$ is constant, it can be replaced by its pointwise transpose $A_1^T$ in the left-hand side,
\[-\nabla\cdot A_1^T\nabla\phi_e\,=\,\nabla\cdot Ae+\nabla\cdot \mathds1_{\Bc(\Pc)}(A_2-A_1)\nabla\phi_e\qquad\text{in $\R^d$}.\]
Testing this equation with $\psi_e^\circ$ then yields
\[\expecm{\nabla\psi_e^\circ\cdot A_1^T\nabla\phi_e}\,=\,-\expecm{\nabla\psi_e^\circ\cdot Ae}-\expecm{\nabla\psi_e^\circ\cdot\mathds1_{\Bc(\Pc)}(A_2-A_1)\nabla\phi_e},\]
or equivalently, adding and subtracting several terms,
\begin{multline*}
\expecm{(\nabla\phi_e+e)\cdot A\nabla\psi_e^\circ}
\,=\,\expecm{e\cdot A\nabla\psi_e^\circ}
-\expecm{\nabla\psi_e^\circ\cdot Ae}\\
+\expecm{\nabla(\phi_e-\psi_e^\circ)\cdot \mathds1_{\Bc(\Pc)}(A_2-A_1)\nabla\psi_e^\circ}
-\expecm{\nabla\psi_e^\circ\cdot\mathds1_{\Bc(\Pc)}(A_2-A_1)\nabla(\phi_e-\psi_e^\circ)}.
\end{multline*}
Inserting this into~\eqref{eq:pre-decomp-phi-psi0}, we get
\begin{multline}\label{eq:pre-decomp-phi-psi1}
\expec{ |\nabla (\phi_e-\psi_e^\circ)|^2}
\,\lesssim\,
\expecm{\nabla\psi_e^\circ\cdot A(\nabla\psi_e^\circ+e)}
-\expecm{\nabla(\phi_e-\psi_e^\circ)\cdot \mathds1_{\Bc(\Pc)}(A_2-A_1)e}\\
+\expecm{\nabla\psi_e^\circ\cdot\mathds1_{\Bc(\Pc)}(A_2-A_1)\nabla(\phi_e-\psi_e^\circ)}\\
-\expecm{\nabla(\phi_e-\psi_e^\circ)\cdot \mathds1_{\Bc(\Pc)}(A_2-A_1)\nabla\psi_e^\circ}.
\end{multline}
We note that the first right-hand side term vanishes: indeed, appealing to~\eqref{eq:Voronoi}, we find
\begin{eqnarray*}
\lefteqn{\expecm{ \nabla\psi_e^\circ\cdot A(\nabla\psi_e^\circ+e)}}\\
&=&\expecM{\sum_{n\in\Sc}\frac{\mathds1_{0\in B(x_n)}}{|B|}\int_{V_n} \nabla\psi_e^\circ\cdot A(\nabla\psi_e^\circ+e)}+\expecM{\sum_{q\in\Sc'}\frac{\mathds1_{0\in J_q}}{|J_q|}\int_{W_q} \nabla\psi_e^\circ\cdot A(\nabla\psi_e^\circ+e)},
\end{eqnarray*}
where for all $n\in\Sc$ and $q\in\Sc'$ the equations for $\psi_e^\circ|_{V_n}=\psi_{e,n}^\circ$ and $\psi^\circ_e|_{W_q}=\psi_{e,q}^{\prime\circ}$ precisely give
\begin{eqnarray*}
\int_{V_n} \nabla\psi_e^\circ\cdot A(\nabla\psi_e^\circ+e)\,=\,0,\qquad\int_{W_q} \nabla\psi_e^\circ\cdot A(\nabla\psi_e^\circ+e)\,=\,0.
\end{eqnarray*}
The estimate~\eqref{eq:pre-decomp-phi-psi1} then leads us to
\[\expec{ |\nabla (\phi_e-\psi_e^\circ)|^2}
\,\lesssim\,\expecm{\mathds1_{\Bc(\Pc)}\big(1+|\nabla\psi_e^\circ|\big)|\nabla (\phi_e-\psi_e^\circ)|}.\]
Appealing to~\eqref{eq:Voronoi}, together with the Cauchy--Schwarz inequality, this yields
\begin{multline}\label{eq:decomp-phi-psi1}
\expec{ |\nabla (\phi_e-\psi_e^\circ)|^2}
\,\lesssim\,\expecM{\sum_{n\in\Sc}\frac{\mathds1_{0\in B(x_n)}}{|B|}\Big(1+\int_{B(x_n)}|\nabla\psi_e^\circ|^2\Big)^\frac12\Big(\int_{B(x_n)} |\nabla (\phi_e-\psi_e^\circ)|^2\Big)^\frac12}\\
+\expecM{\sum_{q\in\Sc'}\frac{\mathds1_{0\in J_q}}{|J_q|}\Big(|J_q|+\int_{J_q}|\nabla\psi_e^\circ|^2\Big)^\frac12\Big(\int_{J_q}|\nabla (\phi_e-\psi_e^\circ)|^2\Big)^\frac12}.
\end{multline}
For all $n\in\Sc$, as $A_1$ is constant, the defining equation~\eqref{eq:psi0n} for $\psi_e^\circ|_{V_n}=\psi_{e,n}^\circ\in H^1_0(V_n)$ yields
\begin{equation*}
\int_{V_n}|\nabla\psi_e^\circ|^2\,\lesssim\,\int_{V_n}\nabla\psi_{e,n}^\circ\cdot A\nabla\psi_{e,n}^\circ
\,=\,-\int_{V_n}\nabla\psi_{e,n}^\circ\cdot Ae
\,=\,-\int_{B(x_n)}\nabla\psi_{e,n}^\circ\cdot (A_2-A_1)e,
\end{equation*}
which implies the energy estimate
\[\int_{V_n}|\nabla\psi_e^\circ|^2\,\lesssim\,1.\]
Similarly, for all $q\in\Sc'$, we find
\begin{equation}\label{eq:energy-Wq-psi0}
\int_{W_q}|\nabla\psi_e^\circ|^2\,\lesssim\,|J_q|.
\end{equation}
Inserting this into~\eqref{eq:decomp-phi-psi1}, we deduce
\begin{multline*}
\expec{ |\nabla (\phi_e-\psi_e^\circ)|^2}
\,\lesssim\,\expecM{\sum_{n\in\Sc}\frac{\mathds1_{0\in B(x_n)}}{|B|}\Big(\int_{B(x_n)} |\nabla (\phi_e-\psi_e^\circ)|^2\Big)^\frac12}\\
+\expecM{\sum_{q\in\Sc'}\frac{\mathds1_{0\in J_q}}{|J_q|}|J_q|^\frac12\Big(\int_{W_q}|\nabla (\phi_e-\psi_e^\circ)|^2\Big)^\frac12}.
\end{multline*}
We now appeal to the following mean-value property, which we shall prove in Step~4 below,
\begin{equation}\label{0-16}
\int_{B(x_n)} |\nabla (\phi_e-\psi_e^\circ)|^2 \,\lesssim\, \rho_n^{-d} \int_{V_n}  |\nabla (\phi_e-\psi_e^\circ)|^2.
\end{equation}
By Young's inequality, we then get for any $K>0$,
\begin{multline*}
\expec{ |\nabla (\phi_e-\psi_e^\circ)|^2}
\,\lesssim\,K\expecM{\sum_{n\in\Sc}\frac{\mathds1_{0\in B(x_n)}}{|B|}\rho_n^{-d}}+K\expecM{\sum_{q\in\Sc'}\frac{\mathds1_{0\in J_q}}{|J_q|}|J_q|}\\
+\tfrac1K\expecM{\sum_{n\in\Sc}\frac{\mathds1_{0\in B(x_n)}}{|B|}\int_{V_n}|\nabla (\phi_e-\psi_e^\circ)|^2}+\tfrac1K\expecM{\sum_{q\in\Sc'}\frac{\mathds1_{0\in J_q}}{|J_q|}\int_{W_q}|\nabla (\phi_e-\psi_e^\circ)|^2}.
\end{multline*}
As~\eqref{eq:Voronoi} implies that the last two right-hand side terms are equal to $\frac1K\expec{ |\nabla (\phi_e-\psi_e^\circ)|^2}$, choosing $K\simeq1$ large enough, and recalling the notation $\Bc(\Pc)\setminus\tilde\Bc(\Pc)=\cup_{q\in\Sc'}J_q$, we finally obtain
\begin{equation}\label{0-20/a}
\expec{ |\nabla (\phi_e-\psi_e^\circ)|^2}\,\lesssim\,\expecM{\sum_{n\in\Sc}\frac{\mathds1_{0\in B(x_n)}}{|B|}\rho_n^{-d}}+\expecm{\mathds1_{\Bc(\Pc)\setminus\tilde\Bc(\Pc)}}.
\end{equation}
To conclude the proof of the claim~\eqref{0-20}, it remains to establish the corresponding estimate on $\expec{\mathds1_{\Bc(\Pc)}|\nabla (\phi_e-\psi_e^\circ)|}$.
For that purpose, we start by appealing to~\eqref{eq:Voronoi} in form of
\begin{multline*}
\expec{\mathds1_{\Bc} |\nabla (\phi_e-\psi^\circ)|}
\\
\,=\,\expecM{\sum_{n\in \Sc} \frac{\mathds{1}_{0 \in B(x_n)}}{|B|} \int_{B(x_n)}  |\nabla (\phi_e-\psi^\circ)| }+ \expecM{\sum_{q\in \Sc'} \frac{\mathds{1}_{0 \in J_{q}}}{|J_{q}|} \int_{J_{q}} |\nabla(\phi_e-\psi^\circ)|}.
\end{multline*}
By the Cauchy--Schwarz inequality and the mean-value property~\eqref{0-16}, we find
\begin{multline*}
\expec{\mathds1_{\Bc} |\nabla (\phi_e-\psi_e^\circ)|}
\,\lesssim\,\expecM{\sum_{n\in \Sc} \frac{\mathds{1}_{0 \in B(x_n)}}{|B|}\rho_n^{-\frac d2}\Big(\int_{V_n}  |\nabla (\phi_e-\psi_e^\circ)|^2\Big)^\frac12}\\
+\expecM{\sum_{q\in\Sc'}\frac{\mathds1_{0\in J_q}}{|J_q|}|J_q|^\frac12\Big(\int_{W_q}|\nabla (\phi_e-\psi_e^\circ)|^2\Big)^\frac12}.
\end{multline*}
By Young's inequality, arguing as above, we obtain
\begin{equation*}
\expec{\mathds1_{\Bc} |\nabla (\phi_e-\psi_e^\circ)|}
\,\lesssim\,\expecM{\sum_{n\in \Sc} \frac{\mathds{1}_{0 \in B(x_n)}}{|B|}\rho_n^{-d}}+\expecm{\mathds1_{\Bc(\Pc)\setminus\tilde\Bc(\Pc)}}+\expec{|\nabla (\phi_e-\psi_e^\circ)|^2},
\end{equation*}
and the claim~\eqref{0-20} now follows from~\eqref{0-20/a}.

\medskip
\step4 Proof of the mean-value property~\eqref{0-16}. \\
We reformulate \eqref{0-16} in the following form: for all $r\ge1$, for all $v\in H^1(rB)$ that satisfies in the weak sense,
\[-\nabla \cdot \big(A_1+\mathds 1_{B}(A_2-A_1)\big) \nabla v=0,\qquad\text{in $rB$},\]
we have
\begin{equation}\label{0-21}
\int_{B} |\nabla v|^2\, \lesssim\, r^{-d} \int_{rB} |\nabla v|^2.
\end{equation}
To prove this result, we decompose $v=v_1+v_2$, where $v_1$ is the unique weak solution in $v+H^1_0(rB)$ of
\[-\nabla \cdot A_1 \nabla v_1=0,\qquad\text{in $rB$},\]
and where $v_2$ is the unique weak solution in $H^1_0(rB)$ of
$$
-\nabla \cdot\big(A_1+\mathds 1_{B}(A_2-A_1)\big)  \nabla v_2 =\nabla \cdot\mathds1_B(A_2-A_1)\nabla v_1,\qquad\text{in $rB$}.
$$
As $A_1$ is constant, we can apply to $v_1$ the mean-value property for harmonic functions, to the effect of
\[\int_{B} |\nabla v_1|^2\, \lesssim\, r^{-d} \int_{rB} |\nabla v_1|^2.\]
As the defining Dirichlet problem for $v_1$ yields the energy estimate
\[\int_{rB}|\nabla v_1|^2\lesssim\int_{rB}|\nabla v|^2,\]
we deduce
\begin{equation}\label{0-22a}
\int_{B} |\nabla v_1|^2\, \lesssim\, r^{-d} \int_{rB} |\nabla v|^2.
\end{equation}
In order to estimate $v_2$, we start from the corresponding energy estimate
\[\int_{rB} |\nabla v_2|^2 \, \lesssim\, \int_{B} |\nabla v_1|^2.\]
Combined with~\eqref{0-22a}, this yields the claim~\eqref{0-21} by the triangle inequality.

\medskip
\step5 Comparison of whole-space and of Dirichlet single-inclusion problems: proof that for all $n\in \Sc$ we have
\begin{equation}\label{0-22}
\int_{B(x_n)} |\nabla (\psi_{e,n}-\psi^\circ_{e})|^2 \,\lesssim\, \rho_n^{-2d}.
\end{equation}
By definition, the difference $\delta \psi_{e,n}:=\psi_{e,n}-\psi^\circ_{e,n}$ satisfies $\delta\psi_{e,n}|_{V_n}\in\psi_{e,n}+H^1_0(V_n)$ and, in the weak sense,
\begin{equation}\label{eq:delta-psin}
-\nabla  \cdot A \nabla \delta \psi_{e,n}=0,\qquad\text{in $V_n$}.
\end{equation}
The mean-value property~\eqref{0-21} on $B_{\rho_n}(x_n)\subset V_n$ then yields
\begin{equation}\label{eq:estim-0st4}
\int_{B(x_n)} |\nabla \delta \psi_{e,n}|^2 \,\lesssim \,  \rho_n^{-d} \int_{B_{\rho_n}(x_n)} |\nabla \delta \psi_{e,n}|^2\,\le\,\rho_n^{-d} \int_{V_n} |\nabla \delta \psi_{e,n}|^2.
\end{equation}
Given a smooth cut-off $\chi_n$ such that
\[\chi_n|_{B_{\frac12\rho_n}(x_n)}\,=\,0,\qquad\chi_n|_{\R^d\setminus B_{\rho_n}(x_n)}\,=\,1,\qquad|\nabla \chi_n|\,\lesssim\, \rho_n^{-1},\]
the Dirichlet problem~\eqref{eq:delta-psin} for $\delta\psi_{e,n}$ implies the energy estimate
\[\int_{V_n}|\nabla\delta \psi_{e,n}|^2\,\lesssim\,\int_{V_n}|\nabla(\chi_n \psi_{e,n})|^2,\]
and thus, as the single-inclusion solution $\psi_{e,n}$ enjoys the same decay properties as in~\eqref{0-11}, we deduce
\[\int_{V_n}  |\nabla \delta \psi_{e,n}|^2\,\lesssim\,  \rho_n^{-d}.\]
Combined with~\eqref{eq:estim-0st4}, this proves the claim~\eqref{0-22}.

\medskip
\step6 Conclusion.\\
Starting from~\eqref{eq:error-estim0} and comparing single-inclusion solutions to their Dirichlet version $\psi_e^\circ$, we have
\begin{multline*}
\big|\bar Ae-(A_1e+\varphi\widehat A_2e)\big|
\,\lesssim\,
\expecm{\mathds1_{\Bc(\Pc)}|\nabla(\phi_e-\psi_e^\circ)|}
+\expecm{\mathds1_{\Bc(\Pc)\setminus\tilde\Bc(\Pc)}\big(1+|\nabla \psi^\circ_e|\big)}\\
+\expecM{\sum_{n\in\Sc}\mathds1_{B(x_n)}|\nabla(\psi_{e,n}-\psi_e^\circ)|}.
\end{multline*}
Using~\eqref{0-20} to estimate the first right-hand side term, appealing to~\eqref{eq:Voronoi} together with the energy estimate~\eqref{eq:energy-Wq-psi0} in form of
\begin{equation*}
\expecm{\mathds1_{\Bc(\Pc)\setminus\tilde\Bc(\Pc)}|\nabla \psi_e^\circ|}\,=\,\expecM{\sum_{q\in\Sc'}\frac{\mathds1_{0\in J_q}}{|J_q|}\int_{J_q}|\nabla \psi_e^\circ|}
\,\lesssim\,\expecM{\sum_{q\in\Sc'}\frac{\mathds1_{0\in J_q}}{|J_q|}|J_q|}\,\le\,\expecm{\mathds1_{\Bc(\Pc)\setminus\tilde\Bc(\Pc)}},
\end{equation*}
and appealing to~\eqref{eq:Voronoi} together with~\eqref{0-22} in form of
\begin{eqnarray*}
\expecM{\sum_{n\in\Sc}\mathds1_{B(x_n)}|\nabla (\psi_{e,n}-\psi_{e}^\circ)|}&=&\expecM{\sum_{n\in\Sc}\frac{\mathds1_{0\in B(x_n)}}{|B|}\int_{B(x_n)}|\nabla (\psi_{e,n}-\psi_{e}^\circ)|}\\
&\lesssim&\expecM{\sum_{n\in\Sc}\frac{\mathds1_{0\in B(x_n)}}{|B|}\rho_n^{-d}},
\end{eqnarray*}
we deduce
\begin{equation}\label{bingo}
\big|\bar A e- (A_1e+\varphi \widehat A_2e)\big|
\,\lesssim\,
\expecM{\sum_{n}\frac{\mathds1_{0\in B(x_n)}}{|B|}\langle\rho_n\rangle^{-d}}
+\expecm{\mathds1_{\Bc(\Pc)\setminus\tilde\Bc(\Pc)}}.
\end{equation}
It remains to evaluate the two right-hand side terms.
For that purpose, we shall appeal to the following observation that we first made in~\cite{2008.03837}: for any non-increasing function \mbox{$g\in\Ld^\infty(\R^+)$} with $g(r)\downarrow0$ as $r\uparrow\infty$, there holds
\begin{equation}\label{0-23}
\expecM{\sum_{n} \mathds1_{0 \in B(x_n)} g(\rho_n)}\,\lesssim\, \lambda_2 \|g\|_{\Ld^\infty(\R^+)}+\int_{0}^\infty |g'(r)|\,\Big((\lambda_2 \,\langle r\rangle^d)\wedge \lambda \Big)\,dr,
\end{equation}
with the short-hand notation $\langle r\rangle:=(1+r^2)^\frac12$.
Applying this to $g(r)=\langle r\rangle^{-d}$, we find
\begin{eqnarray*}
\expecM{\sum_{n} \mathds{1}_{0 \in B(x_n)} \langle\rho_n\rangle^{-d}}&\lesssim&\lambda_2+\int_{0}^\infty\langle r\rangle^{-d-1}\Big((\lambda_2\langle r\rangle^d)\wedge \lambda \Big)\,dr\\
&\lesssim &
\lambda_2 \log(2+\tfrac{\lambda}{\lambda_2}) ~\lesssim~ \lambda_2 |\!\log \lambda_2|.
\end{eqnarray*}
Applying it to $g(r)=\mathds1_{r\le2}$, we further get
\[\expecm{\mathds 1_{\Bc(\Pc)\setminus\tilde \Bc(\Pc)}} \,=\,\expecM{\sum_n\mathds1_{0\in B(x_n)}\mathds1_{\rho_n\le2}}\,\lesssim\,
 \lambda_2.\]
Combined with~\eqref{bingo}, this concludes the proof of~\eqref{0-26+}.

\medskip\noindent
Finally, for completeness, we include a short proof of~\eqref{0-23}. For that purpose, we start by rewriting the left-hand side as
\begin{equation}\label{0-24}
\expecM{\sum_{n}\mathds1_{0\in B(x_n)} g(\rho_n)}\,=\,\int_0^\infty g(r)\,d \Lambda(r),
\end{equation}
where the positive measure $\Lambda$ on $\R^+$ is defined by its distribution function
\[\Lambda([0,r])\,:=\,\expecM{\sum_{n} \mathds1_{0\in B(x_n)} \mathds1_{\rho_n\le r}}
\,=\,\expecM{\sum_{n}\mathds1_{|x_n|<1}\,\mathds1_{\exists m\ne n:\,|x_m-x_n|\le 2r}}.\]
By definition of the minimal length $\ell$ of the point process, cf.~\eqref{0-14}, note that $\Lambda([0,r])=0$ for $r<\frac12\ell$.
Moreover, we can bound
\[\Lambda([0,r])\,\le\, \expecM{\sum_{n} \mathds1_{|x_n|<1}}\,=\, \lambda|B|,\]
and alternatively, for $r\ge\frac12\ell$, by definition of the second-order intensity, cf.~\eqref{eq:high-intens/Re},
\begin{multline*}
\Lambda([0,r])\,\le\, \expecM{\sum_{n\ne m}\mathds1_{|x_n|<1}\,\mathds1_{|x_m-x_n|\le 2r}}
\,=\,\iint_{B\times B_{2r}}f_2(x,x+y)\,dxdy\\
\,=\,(2r)^{-d}\iint_{B_{2r}\times B_{2r}}f_2(x,x+y)\,dxdy
\,\lesssim\,\lambda_2r^{d}.
\end{multline*}
Combining these estimates yields for all $r\ge0$,
\begin{equation}\label{0-25}
\Lambda([0,r])\,\lesssim\,
(\lambda_2 \langle r\rangle^d)\wedge \lambda .
\end{equation}
Under our assumptions on $g$, an integration by parts yields
\[\int_0^\infty g(r)\,d\Lambda(r)\,=\,-g(0)\Lambda(\{0\})+\int_{0}^\infty |g'(r)|\,\Lambda([0,r])\,dr,\]
and the claim~\eqref{0-23} follows in combination with~\eqref{0-24} and~\eqref{0-25}.
\qed

\section*{Acknowledgements}
MD acknowledges financial support from F.R.S.-FNRS, and AG from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant Agreement n$^\circ$~864066).

\bibliographystyle{plain}


\section*{AUTHOR DECLARATIONS}

\subsection*{Acknowledgments}

The authors acknowledge the financial support via the projects C5, A2, and Z3 in the CRC 1375 “NOA–Nonlinear optics down to atomic scales” funded by the Deutsche Forschungsgemeinschaft (DFG) under project-ID 398816777.

\subsection*{Conflict of Interest}

The authors have no conflicts to disclose.

\subsection*{Author Contributions}

\textbf{Francesco Vitale}: conceptualization (equal); investigation – experimental (lead); writing – original draft (lead). \textbf{Daniel Repp}: formal analysis, software – simulations (lead); writing – review \& editing (equal). \textbf{Thomas Siefke}: resources –  sample fabrication (lead); writing – review \& editing (equal). \textbf{Uwe Zeitner}: writing – review \& editing (equal). \textbf{Ulf Peschel}: writing – review \& editing (equal). \textbf{Thomas Pertsch}: funding acquisition (equal); writing – review \& editing (equal). \textbf{Carsten Ronning}: funding acquisition (equal); conceptualization (equal); project administration (lead); supervision (lead); writing – review \& editing (lead).

\section{Dielectric/metal substrate fabrication and EBL patterning}
The dielectric/metal substrates were fabricated by first depositing aluminum films with a thickness of $t \approx$ 68 nm on a fused silica substrate via ion beam deposition (Ionfab 300LC+ by Oxford Instruments). High relative densities and good adhesion were achieved with this method. However, the high normalized energy flux and the generalized temperature during the film growth fostered diffusion processes and, hence, recrystallization \cite{Anders}. Therefore, materials with a low melting point, such as aluminum, tend to become polycrystalline (PC) with very large grains, leading to high surface roughness. Figure \ref{Fig.S1}(a) shows an AFM image of an as-deposited Al layer: the RMS surface roughness was found to be $R_q \approx$ 3.8 nm, calculated from measurements conducted in tapping mode over a $10 \upmu \textup{m} \times 10 \upmu \textup{m}$ area with a resolution of 40 nm. This resulted in a decreased propagation length and increased scattering for the surface plasmon polaritons at the spacer/nanowire interface \cite{Raether}. In order to reduce the roughness, 10 at\% silicon was co-deposited by switching between an aluminum and a silicon target. The deposition time was 25 s at an Ar-ion energy of 1.3 keV and a current of 250 mA for aluminum, and 10 s at an Ar-ion energy of 1.3 keV and a current of 100 mA for silicon respectively. The cycles were repeated 13 times and ended with a Si-free termination layer of thickness $t \approx$ 5-6 nm. On top of this termination layer, a native aluminum oxide layer formed with an estimated thickness of $t \approx$ 2-3 nm.
The co-deposition with silicon led to the reduction of the surface roughness by about one order of magnitude down to $R_q \approx$ 0.2 nm: the corresponding AFM height profile is shown in Figure \ref{Fig.S1}(b).

\begin{figure} [ht!]
\centering
    \includegraphics[width=\textwidth]{figures/S1.png}
        \caption{AFM height profiles for (a) pristine and (b) 10\% Si co-deposited polycrystalline Al substrates. The analysis revealed an average RMS surface roughness $R_q  \approx$ 3.8 nm for the pristine Al layer and  $R_q  \approx$ 0.2 nm for the co-deposited one.
}
        \label{Fig.S1}
\end{figure}

 Figure \ref{Fig.S2} shows the respective fractions of lasing yield for a sample of 30 nanowires on top of Al substrates with and without the Si co-deposition. As a consequence of the smoother surface, the yield of nanowires that showed plasmonic-like lasing at room temperature, significantly increased. This can be ascribed to the larger propagation length and reduced scattering of the surface plasmon polaritons at the spacer/nanowire interface \cite{Raether}, together with an overall better thermal dissipation that increased the optical damage threshold.

\begin{figure} [ht!]
\centering
    \includegraphics[width=\textwidth]{figures/S2.png}
        \caption{Lasing yield for nanowires on top of (a) pristine and (b) Si co-deposited PC aluminum layers. The improvement of the surface roughness led to an increase in the yield of devices exhibiting plasmonic-like lasing action.
}
        \label{Fig.S2}
\end{figure}

Subsequently, the grating structures were realized by character projection electron beam lithography (Vistec SB350 OS) \cite{Kley}. The mask, consisting of a positive tone CAR FEP171 (Fujifilm) was transferred onto the aluminum film by ion beam etching (Ionfab 300LC+, Oxford Instruments). The resist was, then, stripped by reactive ion etching in an oxygen atmosphere (Sentech SI-591), and, ultimately, the samples were laser diced.

\section{Sample preparation and single-nanowire nano-manipulation}

\begin{figure} [ht!]
\centering
    \includegraphics[width=\textwidth]{figures/S3.png}
        \caption{(a) Optical microscope image of a section of the Al$_2$O$_3$/Al substrate, exhibiting an exemplary EBL-patterned metallic grating. (b) Photograph of the nano-manipulation setup: the arrow points to the manipulator arm, equipped with the whisker mounting a W tip. (c) Representative optical image of a single nanowire and the W tip, acquired prior to the nano-manipulation. The same nanowire was successfully placed on top of the grating with a ridge (trench) orientation of $\theta \approx 50\degree$, as shown in (d).
}
        \label{Fig.S3}
\end{figure}

The ZnO nanowires were dispersed onto the EBL-patterned Al$_2$O$_3$/Al substrates, shown in Figure \ref{Fig.S3}(a), via a dry imprint technique consisting in gently rubbing the growth sample surface against that of the target substrate, which was, in turn, imaged under an optical microscope (Zeiss Axiolab 5). The mechanical nano-manipulation was performed in air under the PL microscope, as shown in the photograph of \ref{Fig.S3}(b), in order to avoid the modification of the optical properties induced by the electron beam irradiation, as discussed in the main text of the manuscript. To this end, we made use of a remotely controlled nano-manipulator (Kleindiek MM3A-EM) featuring a W tip with a radius of curvature of about 1 micron installed in between the sample and the 50x NUV objective. The sample was imaged via a white light fiber illuminator (Thorlabs OSL1-EC) on a microscopy camera (PixeLink PL-B776U). After the tip was brought in contact with the substrate, finer displacement steps were achieved via a remotely controlled piezo-controller (Thorlabs BPC303) connected to the piezo-actuators of the sample stage. This allowed for an accurate “nudging” of the single nanowire toward the metallic grating and a successful placement on top of it, as shown in Figures \ref{Fig.S3}(c)(d), in order to avoid abrupt and uncontrolled movements of the W tip, eventually leading to the nanowire breakage.

\section{Further optical measurements of nanowires on gratings}

\begin{figure} [ht!]
\centering
    \includegraphics[width=\textwidth]{figures/S4.png}
        \caption{PL spectra of a ZnO nanowire with $L \approx$  13 $\upmu$m and $\langle d \rangle \approx$ ≈ 230 nm placed on a grating with $p$ = 300 nm and $FF$ = 0.5 at different configurations.
}
        \label{Fig.S4}
\end{figure}

Figure \ref{Fig.S4} shows the PL spectra from another nanowire featuring a length $L \approx$ 13 $\upmu$m and average diameter $\langle d \rangle \approx$ 230 nm, placed on a grating with pitch $p$ = 300 nm and fill factor $FF = 0.5$. Also in this case, we can clearly observe the emergence of additional peaks in the low-energy side of the gain envelope, for an orientation of the nanowire perpendicular to the metallic ridges. The additional peaks disappeared for a subsequent re-orientation of the wire at $\theta \approx 45\degree$ as well as for the planar situation. We attribute the slight shift ($\approx$ 1 nm) of the spectral location of the long-wavelength peak, compared to the other two nanowires discussed in the main text, to the variations related to the cavity imperfections (both nanowire- and grating-related).

\section{Active lasing simulations}

The active lasing simulations were accomplished by choosing the 4-level-2-electron material model provided by Lumerical for the nanowire material. A non-dispersive dielectric with an index of $n$ = 2.4 was chosen as a base material. The nanowire was illuminated by a Total-Field Scattered-Field (TFSF) source, modeling the experimental pump, for an excitation time $\tau_{\textup{exc}}$ = 500 fs at a wavelength $\lambda_{\textup{exc}}$ = 350 nm, in resonance with the pump transition in the material. The pump was, then, switched off, leading to the excitation of temporal Fourier components resonant with the lasing transition, thereby achieving a broad-band seeding of the lasing modes. This resulted in the appearance of a lasing pulse, which experienced gain. The lasing wavelength was identified via the observation of a maximum in the spectral power distribution from a monitor placed perpendicularly to the nanowire axis in a transversal cross-section. The field distribution was given by a second monitor plane, spanned by the nanowire axis and the surface normal to the substrate.


\section{Hot, massive stars and their ionizing fluxes}

Massive stars spend most of their lifetime as hot stars with $T_\mathrm{eff} > 10\,$kK. Consequently, these stars have their flux maximum at ultraviolet (UV) wavelengths, providing an efficient source of flux for driving a stellar wind and ionizing their surrounding environment. Photometric and spectroscopic observations in the UV are particularly crucial to determine the parameters of hot, massive stars and their outflows. As the extreme ultraviolet (EUV) is largely inaccessible due to the presence of interstellar hydrogen, the far UV regime between 900 and 2000\,\AA\ contains the most viable diagnostics to photometrically constrain $T_\mathrm{eff}$ or infer the wind properties and ionizing fluxes if spectroscopy is available.

The ionizing fluxes of hot stars are commonly quantified by integrating the photon flux beyond a frequency $\nu_\mathrm{edge}$, namely
\begin{equation}
  \label{eq:qdef}
   Q_\mathrm{edge} := \int\limits_{\nu_\mathrm{edge}}^\infty \frac{F_\nu}{h \nu} \mathrm{d}\nu
\end{equation}
with typically $\log Q_\mathrm{edge}$ given in the literature, taking the logarithm of the value of photons per second. Commonly used are $\log Q_0 \equiv \log Q_\mathrm{H\,\textsc{i}}$ for the hydrogen-ionizing flux beyond the Lyman edge, and $\log Q_1 \equiv \log Q_\mathrm{He\,\textsc{i}}$ as well as $\log Q_2 \equiv \log Q_\mathrm{He\,\textsc{ii}}$ for the ionizing fluxes beyond the edges of He\,\textsc{i} and He\,\textsc{ii} ionization.

\section{Measuring Ionizing Fluxes}

Given the EUV inaccessibility, ionizing fluxes are mostly obtained by either measuring nebula emission lines or integrating the spectral energy distribution (SED) from an atmosphere model assigned to the source star(s) \citep[e.g.,][]{Ramachandran+2018}.

\begin{figure}[htb]
\begin{center}
   \includegraphics[angle=0,width=0.49\textwidth]{blackbody-smc-o3v-proc}
	 \hfill
   \includegraphics[angle=0,width=0.49\textwidth]{ostar-ionflux-comparison}
     \caption{Left: Comparison of the spectral energy distribution (SED) from an atmosphere model (red) for an O3 SMC dwarf to a blackbody (black) of the same $T_\mathrm{eff}$. Right: Comparison of $Q_\mathrm{H\,\textsc{i}}$ and $Q_\mathrm{He\,\textsc{i}}$ for OB star models from \citet{Smith+2002} and \citet{Martins+2005}.}
   \label{fig:ionflux-intro}
\end{center}
\end{figure}

The use of atmosphere models is crucial as estimates from blackbodies can be severely off (cf.\ Fig.\,\ref{fig:ionflux-intro}). For $Q_\mathrm{H\,\textsc{i}}$ and $Q_\mathrm{He\,\textsc{i}}$, spectral subtype calibrations vary a bit among the literature (cf.\ right panel Fig.\,\ref{fig:ionflux-intro}) with many uncertainties existing for lower metallicities.

\section{Interplay between stellar winds and He\,II ionizing flux}

In contrast to $Q_\mathrm{H\,\textsc{i}}$ and $Q_\mathrm{He\,\textsc{i}}$, the He\,\textsc{ii} ionizing flux $Q_\mathrm{He\,\textsc{ii}}$ depends strongly on the wind strength of the stars and thus also on the metallicity $Z$. In strong, dense winds, He\,\textsc{iii} recombines, making the atmosphere opaque for $Q_\mathrm{He\,\textsc{ii}}$ photons \citep[e.g.][]{Schmutz+1992}. The absorbed flux also drives the wind and is re-emitted only at longer wavelengths.

\begin{figure}[htb]
\begin{center}
   \includegraphics[angle=0,width=0.49\textwidth]{sed-ionflux-NGC346-SSN9}
	 \hfill
   \includegraphics[angle=0,width=0.49\textwidth]{sed-ionflux-R136a1}
     \caption{SED highlighting the flux beyond the ionization edges of H\,\textsc{i}, He\,\textsc{i}, and He\,\textsc{ii} for the O2 giant SSN9 in the SMC \citep[left panel, model from][]{Rickard+2022} and the WN5h star R136a1 \citep[right panel, model based on][]{Hainich+2014}.}
   \label{fig:ionflux-SSN9-R136a1}
\end{center}
\end{figure}

Apart from a high intrinsic temperature, stellar sources for $Q_\mathrm{He\,\textsc{ii}}$ thus require winds that are mainly optically thin. Consequently, most Wolf-Rayet (WR) stars are not very efficient sources of $Q_\mathrm{He\,\textsc{ii}}$, regardless of whether they are classical, helium-burning WRs or very massive stars (VMS) that might still be hydrogen-burning. This is illustrated in Fig.\,\ref{fig:ionflux-SSN9-R136a1}, where an O2 giant in the SMC provides orders of magnitude more He\,\textsc{ii} ionizing flux than the hotter and five times more luminous WN5h-star R136a1 in the LMC.

\begin{figure}[htb]
\begin{center}
   \includegraphics[angle=0,width=0.49\textwidth]{sed-ionflux-SMC-AB1}
	 \hfill
   \includegraphics[angle=0,width=0.49\textwidth]{sed-ionflux-SMC-AB8}
     \caption{Same as Fig.\,\ref{fig:ionflux-SSN9-R136a1}, but for SMC AB\,1 \citep[WN3ha, left panel, model based on][]{Hainich+2015} and the WO4-star in SMC AB\,8 \citep[right panel, model based on][]{Shenar+2016}.}
   \label{fig:ionflux-AB1-AB8}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
   \includegraphics[angle=0,width=0.49\textwidth]{sed-ionflux-azv476-prim-obs}
	 \hfill
   \includegraphics[angle=0,width=0.49\textwidth]{sed-ionflux-azv476-prim-mesa}
     \caption{SEDs for the observed, partially stripped O4 primary in the SMC AzV 476 system (left panel) and the corresponding prediction based on a binary evolution model (right panel). The underlying models for the SED plots are taken from \citet{Pauli+2022}.}
   \label{fig:ionflux-AzV476}
\end{center}
\end{figure}

Nonetheless, some WR stars actually contribute huge amounts of $Q_\mathrm{He\,\textsc{ii}}$, such as weak-winded, early-type WNs and WO stars (see Fig.\,\ref{fig:ionflux-AB1-AB8} for examples). In all cases, the wind density is key here, making predictions difficult, not least due to the uncertainties in $\dot{M}$. An example involving binary evolution is shown in Fig.\,\ref{fig:ionflux-AzV476}, where the empirically obtained wind by \citet{Pauli+2022} is stronger than predicted by the evolutionary scenario for the system, leading to a reduction of $Q_\mathrm{He\,\textsc{ii}}$ by three orders of magnitude.

In summary, stellar contributors of He\,\textsc{ii} ionizing flux require high temperatures and relatively thin winds. Known strong contributors are the most massive main sequence stars at low $Z$ as well as classical WR stars with weak winds (e.g. WN3ha, WO). A population of hot, hydrogen-depleted stars below the WR regime \citep[``stripped stars'', see e.g.][]{Goetberg+2020} could be important as well, but yet lacks observational confirmation.

\def\apj{{ApJ}}
\def\nat{{Nature}}
\def\jgr{{JGR}}
\def\apjl{{ApJ Letters}}
\def\aap{{A\&A}}
\def\aaps{{A\&A Supplement}}
\def\mnras{{MNRAS}}
\def\aj{{AJ}}
\def\pasa{{PASA}}
\def\pasp{{PASP}}
\def\ssr{{Space Science Reviews}}
\let\mnrasl=\mnras


\section{Introduction}

Over the past decade, the MAX phase materials which are mainly composed of layered ternary carbides and nitrides are getting huge attention due to their unique properties that arise from the presence of both covalent and metallic bonds in the structure \cite{BARSOUM2000201, JEITSCHKO1964133}. The MAX phases (M\textsubscript{n+1}AX\textsubscript{n}, where M is early transition metals, A is a representative elements, X is C or N, and n = 1,2,3 $\cdots$) are ternary carbides and nitrides with a family of more than 80 members \cite{HADI2020109275}. MAX consists of early transition metals and group A elements, which are highly abundant and not hazardous to the environment. The wide availability and nontoxic nature of building elements of the MAX phase make the MAX phase cost-effective and environmentally friendly. These materials have a wide range of possible applications, including structural materials for high temperature applications, protective coatings and bond-coats for gas turbines, fuel cladding, solar receivers, electrical contacts, catalysts, and many more.
Depending on the different property between M-A (metallic) and M-X (covalent) bonding, MAX phases exhibit both ceramic and metallic property, including high electrical conductivity, oxidation resistance, catalysis, electronics, gas sensing, and energy storage \cite{doi:10.1021/acs.accounts.7b00481, doi:10.1021/acsenergylett.6b00247, doi:10.1021/acsnano.9b06394}.
Some MAX phases, such as Ti\textsubscript{3}AlC\textsubscript{2}, Ti\textsubscript{2}AlC, Cr\textsubscript{2}AlC, also show exceedingly good oxidation resistance and oxidation-induced self-healing ability \cite{FARLE20171969}, and can be used as high-temperature structural materials. Moreover, MAX phase materials are expected to be applied to cladding and structural materials of the rapid reactor for molten lead cooling due to excellent resistance to neutron and He, Ar, Kr, I, Xe, Au plasma irradiation damage ability and excellent resistance to molten lead corrosion \cite{LIU201823254}. Also, some MAX phase films, for example, Ti\textsubscript{3}SiC\textsubscript{2} films, are able to be applied to the surface of metal bipolar plates to improve the corrosion resistance and conductivity of metal bipolar plates as well a considerable application prospect in commercial fuel cells\cite{HANSON2019195}. MAX phases contain two alternately stacked fundamental structural units: the non-stoichiometric transition metal carbide or nitride slabs in NaCl-type crystal structure and the close-packed A-group atomic plane \cite{BENTZEL20154107}.

MAX phase has been mostly studied with metal hydride to enhance the hydrogen storage capacity or production. 7 wt\% Ti\textsubscript{3}AlC\textsubscript{2} MAX-phase composition with Mg significantly reduce dehydrogenation temperature of the composite \cite{LAKHNIK20227274}. Liu et al. show Ti\textsubscript{3}AlC\textsubscript{2}/Pd catalyst with 3 wt\% Pd loading had a much higher capacity for hydrogen production than conventional Pd nanoparticles \cite{nano12050843}. MAX phase behaves as a support structure for palladium nanoparticles for the hydrogen generation from alkaline formaldehyde solution at room temperature. The layered structures of the MAX phase can also be useful in the hydrogen storage field. Very few studies are available on the hydrogen storage performance of the MAX phase or MXene. Hydrofluoric acid incompletely etched multilayered MXene phase shows $\sim4$ wt\% hydrogen storage capacity at ambient condition and Ti\textsubscript{2}CT\textsubscript{x} (T is a functional group) can store 8.8 wt\% at room temperature and 60 bar H\textsubscript{2} \cite{Liu2021}. Single layer MXene can store up to 8.6 wt\% hydrogen by chemisorption of the H atom (1.7 wt \%), physisorption of the H\textsubscript{2} molecule (3.4 wt \%), and Kubas-type binding of the H\textsubscript{2} molecule (3.4 wt \%) \cite{doi:10.1021/jp409585v}. Ding et al. studied the hydrogen incorporation mechanism into the MAX phase Ti\textsubscript{3}AlC\textsubscript{2} by first-principles simulation and find out Ti-Al layers are suitable for hydrogen incorporation \cite{DING20166387}.

Among all phases, Ti\textsubscript{2}AlC is one of the most light-weight and oxidation-resistant layered ternary carbides and it has been extensively studied by many researchers \cite{WANG2021109788, https://doi.org/10.1111/jace.17582, BADIE2021117025}. Light-weight material can enhance hydrogen storage capacity in any structure. In this work, we investigate the hydrogen storage capacity of Ti\textsubscript{2}AlC in MAX phase and also possible doped 211 MAX structure with higher hydrogen storage capacity. We have employed first-principles simulation along with cluster expansion approach to find out possible doped MAX structures. We find that Hydrogen can be only incorporated in Ti-Al layers. Hydrogen in Ti-C layers is unstable and easily migrate to Ti-Al layers. 83.33\% doped TiAl\textsubscript{x}Cu\textsubscript{1-x}C structure has higher hydrogen storage capacity than Ti\textsubscript{2}AlC and Ti\textsubscript{2}CuC structure. In addition, the doped TiAl\textsubscript{x}Cu\textsubscript{1-x}C structure is very stable in nature and can be synthesized experimentally.

\section{Computational details}

\subsection{\textit{ab-initio} calculations}

All first-principles calculations have been performed using the density functional theory (DFT) as implemented in the Vienna Ab initio Software Package (VASP) \cite{Kresse1996CMS, Kresse1996PRB}. For the core and valence electrons, projected augmented wave (PAW) potentials and plane-wave basis sets were used with the Perdew–Burke–Ernzerhof (PBE) functional to get the electronic energy \cite{Kresse1999, Bl1994}. Corrections were made to account for long-range interactions using the semi-empirical Grimme D2 dispersion method and for non-spherical contributions to the PAW potentials that were built into the code \cite{Grimme2006}. All energies converged within a cutoff of 500 eV. The conjugate gradient algorithm is used for structural optimization \cite{Perdew1996}. The convergence criteria for energy and force are 10\textsuperscript{-7} eV and 0.01 eV\AA\textsuperscript{-1}, respectively. In all calculations, spin polarisation was enabled, with the only exception of the isolated H\textsubscript{2} closed shell molecule.

The adsorption energy (E\textsubscript{ads}) of hydrogen on the different MAX structure were calculated by:

\begin{equation}
	E_{ads} =\frac{1}{n}\times\left(E_{MAX+nH}-\left( E_{MAX}+\frac{n}{2}E_{H_2}\right)\right)
	\label{eq:1}
\end{equation}

where E\textsubscript{MAX+nH} is the energy of the hydrogenated MAX structure, E\textsubscript{MAX} and E\textsubscript{H\textsubscript{2}} are the energies of the pristine structure and isolated H\textsubscript{2}, respectively, and n is the number of hydrogen atoms involved in the adsorption.

The atomic charges for all structures were calculated by Bader charge analysis, as suggested by Henkelman and co-workers \cite{Henkelman2006, Tang2009}. The formation energy of Ti\textsubscript{2}Al\textsubscript{x}Cu\textsubscript{1-x}C structure at different Cu concentration ere calculated by:

\begin{equation}
	E_{f} = E_{Ti\textsubscript{2}Al\textsubscript{x}Cu\textsubscript{1-x}C} - x E_{Ti\textsubscript{2}AlC} - (1-x) E_{Ti\textsubscript{2}CuC}
	\label{eq:2}
\end{equation}

where E is the energy of the pure or doped structure and x is the atomic fraction in the doped structure.

In order to estimate the release temperature, the chemical potential of the hydrogen in the solid phase and in the gas phase needs to be evaluated. The chemical potential energy of hydrogen as a function of temperature (T) and pressure (p) was obtained using standard statistical thermodynamics, via the formula \cite{Khan2016}:

\begin{equation}
	\mu(T,p)=\frac{1}{2}\left(E_{el}+E_{ZPE}+H^0(T)-H^0(0)-TS^0(T)+k_BTln\left(\frac{p}{p^0}\right)\right)
	\label{eq:u}
\end{equation}

where E\textsubscript{el} and E\textsubscript{ZPE} are the electronic and zero point energy of the hydrogen molecule as derived from DFT, H\textsuperscript{0} and S\textsuperscript{0} is  enthalpy and entropy at standard pressure p\textsuperscript{0} = 1 bar. The values of enthalpy and entropy at standard pressure can be found in JANAF thermochemical tables \cite{Chase1982}. The critical hydrogen chemical potential, $\mu^C_H$, at which the different MAX structure release hydrogen from the structure is given by $\mu^C_H=1/n(E_{MAX+nH}-E_{MAX})$ \cite{Williamson2004}.

\subsection{Calculation of formation energies of the solid solutions using the Cluster expansion approach}
A DFT-based technique cannot efficiently scan the energies of the immensely large configuration space that includes doped materials. A simple but efficient way is to construct an effective Hamiltonian using a relatively smaller database of DFT energies for these systems, which can then be used to scan a large configuration space. The cluster expansion method is an widely used method and can be used to calculate the energy of a large number of doped structures \cite{PhysRevB.49.8627, BADIE2021117025}.  The total energy of an alloy with configuration $\bm{\sigma}$ is calculated by constructing the so called cluster expansion effective Hamiltonian which is given by~\cite{CE1,CE2,CE3}

\begin{equation}
	E^{CE}(\bm{\sigma}) = J_{0} + \sum_{\alpha}J_{\alpha} \Biggl \langle \prod_{i \in \alpha} \sigma_{i} \Biggr \rangle=J_0+\sum_\alpha J_\alpha\Pi_\alpha
	\label{eq:3}
\end{equation}
Here J\textsubscript{$\alpha$} is known as the effective cluster interactions (ECI). J\textsubscript{0} represents ECI of an empty cluster. $\alpha$ can represent a single (point), pair or triplet cluster of sites. Each lattice site (i) is assigned an occupation variable $\sigma_{i}$ depending on whether it is occupied by an Al or Cu atom. $\Pi_\alpha$ are the cluster correlation functions.
To create a DFT database of Cu-doped Ti\textsubscript{2}AlC structures, the DFT energy of Cu doped 155 random structures of Ti\textsubscript{2}Al\textsubscript{x}Cu\textsubscript{1-x}C was calculated. Of these, 140 configurations were used to train the cluster expansion model and the remaining 15 were used to test the model. We performed 10-fold cross-validations to avoid over-fitting the data.
\par
To obtain the ECI's and the model we use the least absolute shrinkage and selection operator (LASSO) approach where the following minimization is performed to obtain the optimal $J_\alpha$,
\begin{equation}
{\bf J}_{opt}=\underset{{\bf J}}{\operatorname{arg min}} ||{\bf E_{DFT}}-\bf \Pi \bf J||^2_2+\lambda||J||_1
\label{lasso}
\end{equation}
Here $\lambda||{\bf J}||_1$ is the $l_1$ regularization term and $\lambda$ is an hyper-parameter. To obtain new stable Cu-doped materials we follow the following strategy: first we fit the DFT energies in a satisfactory way to the cluster expansion Hamiltonian given by Eq.\ref{eq:3}, next we generate a large pool of doped configurations (1000 ) and calculate the total energies of  such (Ti\textsubscript{2}Al\textsubscript{x}Cu\textsubscript{1-x}C) configurations using the fitted cluster expansion Hamiltonian. These energies  for different compositions are then further used to calculate the formation energies. The formation energy of all compositions was calculated using equation \ref{eq:2}.

\section{Results and Discussion}

In Ti\textsubscript{2}AlC has several interstitial sites and all of them fall into two parts: the interstitial sites in the Ti-Al and Ti-C layers. To find out possible hydrogen adsorption sites, all possible interstitial sites are considered. Figure \ref{fig:ti2alc} shows the possible H occupation sites in both Ti-Al and Ti-C layers. There are three kinds of interstitial sites in Ti-Al layers: the tetrahedral interstitial site I\textsubscript{tetr}-1 consists of three Al atoms and one Ti atom, the tetrahedral interstitial site I\textsubscript{tetr}-2 consists of three Ti atoms and one Al atom, and the octahedral interstitial site I\textsubscript{oct}-3 consists of three Al atoms and three Ti atoms. The Ti-C layers have two tetrahedral interstitial sites I\textsubscript{tetr}-4 and I\textsubscript{tetr}-5 consisting of three Ti atoms and one Ti atom from different Ti-C layers. The interstitial sites are designated as 1, 2, 3, 4, and 5, respectively in Fig. \ref{fig:ti2alc}.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=1.0\linewidth]{fig/ti2alc.png}
	\caption{Crystal structure of 211 MAX phases and interstitial sites. Blue, purple, black, and white balls correspond to Ti, Al, C, and H, respectively.}
	\label{fig:ti2alc}
\end{figure}

The Hydrogen adsorption energy at different sites is calculated by Equation \ref{eq:1} and the results are shown in Table \ref{table:1}. Hydrogen atoms are exothermically adsorbed at all interstitial sites in Ti-Al layers. Among all interstitial sites in Ti-Al layers, I\textsubscript{tetr}-2 is more stable than I\textsubscript{oct}-3 and I\textsubscript{tetr}-1. In a relaxed structure, the H atom at I\textsubscript{tetr}-1 moves to the Al layers and locates in the center of the triangular interstitial site consisting of three atoms. In one Ti\textsubscript{2}AlC formula unit, a total of 6 interstitial sites are present in Ti-Al layers. Among 6 interstitial sites, two of each interstitial type (I\textsubscript{tetr}-1, I\textsubscript{tetr}-2, and I\textsubscript{oct}-3) are present. H atom in Ti-C layers is unstable and H atom move to I\textsubscript{tetr}-2 interstitial site after relaxation from I\textsubscript{tetr}-4 and I\textsubscript{tetr}-5.

\begin{table}[htbp]
	\caption{The adsorption energy of hydrogen at different interstitial sites of Ti\textsubscript{2}AlC}
	\centering
	\begin{tabular}{c c c c c c c}
		\hline
		\textbf{Sites of hydrogen} & \multicolumn{3}{c}{\textbf{Ti-AL layer}} & & \multicolumn{2}{c}{\textbf{Ti-C layer}} \\
		\cmidrule{2-4}
		\cmidrule{6-7}
		& I\textsubscript{tetr}-1 & I\textsubscript{tetr}-2 & I\textsubscript{oct}-3 & & I\textsubscript{tetr}-4 & I\textsubscript{tetr}-5 \\
		\hline
		E\textsubscript{ads}(eV) & -0.06 & -0.40 & -0.23 & & unstable & unstable \\
		\hline
	\end{tabular}
	\label{table:1}
\end{table}

In the current study, we examined how the performance of hydrogen storage is affected when aluminium (Al) in the MAX phase is replaced with 3d-transition metals such as Ti, Cr, Mn, Fe, Co, Ni, Cu, and Zn atoms. Hydrogen adsorption energy was calculated by placing one hydrogen atom at I\textsubscript{tetr}-2 site in different Al substituted MAX phase. Figure \ref{fig:alloy} shows adsorption energy of different 211 MAX phase. Among all MAX phase Ti\textsubscript{2}CuC has the highest hydrogen adsorption energy of -0.76 eV/H. For hydrogen storage calculation, we compared hydrogen capacity between Ti\textsubscript{2}AlC and Ti\textsubscript{2}CuC phase.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=1.0\linewidth]{fig/alloy.png}
	\caption{Hydrogen adsorption energy at I\textsubscript{tetr}-2 in different 211 MAX phases.}
	\label{fig:alloy}
\end{figure}

Next, we calculated electronic structure of Ti\textsubscript{2}AlC and Ti\textsubscript{2}CuC phase. As Ti has lower electronegativity and high oxidation states electron may transfer from Ti\textsubscript{2}C layer to Al or Cu atom. Figure \ref{fig:dos} shows total electronic density of states (DOS) and projected DOS of atomic orbitals of Ti\textsubscript{2}AlC and Ti\textsubscript{2}CuC MAX phases. In Ti\textsubscript{2}AlC phase, the Fermi level is contributed by not only Ti 3d but also Al 3s and Al 3p orbitals (fig. \ref{fig:dos} (a)). With increasing atomic number in the MAX phase, the contribution of Cu atom to the Fermi level decreases. The band center of s, p, and orbitals of Cu is located much deeper than Al (fig. \ref{fig:dos} (b)). The Ti 3d orbital of Ti\textsubscript{2}CuC phase contributed more to the Fermi level.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=1.0\linewidth]{fig/dos.png}
	\caption{Electronic DOS of Ti\textsubscript{2}AlC (a) and Ti\textsubscript{2}CuC (b) structure.}
	\label{fig:dos}
\end{figure}

To find out the hydrogen adsorption mechanism, electron transfer between the hydrogen atom and MAX phase was calculated. The Bader charge analysis of Ti\textsubscript{2}AlC, Ti\textsubscript{2}CuC, and their hydrides a listed in Table \ref{table:2}. In the MAX phase, hydrogen is negatively charged because it accepts electron from neighboring atoms. The electron transfer to H is mostly contributed by Al and Cu atoms in Ti\textsubscript{2}AlC and Ti\textsubscript{2}CuC phase, respectively. In Ti\textsubscript{3}A tetrahedral, hydrogen present as H\textsuperscript{-} anion. Hydrogen in I\textsubscript{tetr}-2 site shows a greater stabilization effect through the Coulomb interaction. Hydrogen in Ti\textsubscript{2}CuC accepts lesser electrons than hydrogen in the Ti\textsubscript{2}AlC phase, which indicates the hydrogen bonding characteristic in the MAX phase changes from ionic to metallic. The bonding of hydrogen in the Ti\textsubscript{2}AlC phase is more ionic in nature, whereas in the Ti\textsubscript{2}CuC phase the bonging is metallic in nature.

\begin{table}[htbp]
	\caption{The Bader charge analysis of Ti\textsubscript{2}AlC, Ti\textsubscript{2}CuC, and their hydrides.}
	\centering
	\begin{tabular}{c c c c c c}
		\hline
		 & & \textbf{Ti} & \textbf{Al} & \textbf{C} & \textbf{H}\\
		\hline
		\multirow{2}{4em}{Ti\textsubscript{2}AlC} & Ti\textsubscript{4}Al\textsubscript{2}C\textsubscript{2} & 2.03 & -1.21 & -2.84 & \\
		 & Ti\textsubscript{4}Al\textsubscript{2}C\textsubscript{2}H\textsubscript{1/2} & 2.03 & -1.06 & -2.82 & -1.34\\
		\hline
		 & & \textbf{Ti} & \textbf{Cu} & \textbf{C} & \textbf{H}\\
		\hline
		\multirow{2}{4em}{Ti\textsubscript{2}CuC} & Ti\textsubscript{4}Cu\textsubscript{2}C\textsubscript{2} & 2.03 & -1.31 & -2.75 & \\
		 & Ti\textsubscript{4}Cu\textsubscript{2}C\textsubscript{2}H\textsubscript{1/2} & 2.04 & -1.21 & -2.77 & -0.79\\
		\hline
	\end{tabular}
	\label{table:2}
\end{table}

To incorporate hydrogen in the material, the diffusion path of hydrogen in the material is necessary as the adsorption energy of the hydrogen in the materials. Ding et al. \cite{DING20166387} shows hydrogen diffusion in Ti-A layers of Ti\textsubscript{3}AlC\textsubscript{2} phase is following. Hydrogen firstly diffused from I\textsubscript{tetr}-1 to a nearest I\textsubscript{tetr}-2 (Path I), then diffused to a nearest I\textsubscript{oct}-3 (Path II), and finally to nearest I\textsubscript{tetr}-1 site (Path III). In this study same path is followed to find out the suitable MAX phase for easier diffusion of hydrogen among Ti\textsubscript{2}AlC and Ti\textsubscript{2}CuC. Figure \ref{fig:neb} (a) shows the H energy diffusion barrier for Ti\textsubscript{2}AlC. The energy diffusion barrier between I\textsubscript{tetr}-1 and I\textsubscript{tetr}-2 site is 0.53 eV. However, the opposite direction has a higher diffusion energy barrier of 0.87 eV. The diffusion energy barrier for I\textsubscript{tetr}-2 to I\textsubscript{oct}-3 is 0.32 eV, whereas the opposite path is easier with energy barrier of 0.15 eV. The path between I\textsubscript{oct}-3 and I\textsubscript{tetr}-1 is also easier than the path I. The energy barrier between I\textsubscript{oct}-3 and I\textsubscript{tetr}-1 is 0.44 eV and the opposite path is much easier with barrier energy 0.27 eV. Based on the energy barrier value, I\textsubscript{oct}-3 behaves as an intermediate position for hydrogen adsorption and desorption in Ti\textsubscript{2}AlC. The H atom of both I\textsubscript{tetr}-1 and I\textsubscript{tetr}-2 first diffuse to nearby I\textsubscript{oct}-3, then move to other I\textsubscript{tetr}-1 and I\textsubscript{tetr}-2 site. The activation energy of hydrogen diffusion in Ti-Al layers of Ti\textsubscript{2}AlC is much lower than $\alpha$-Ti and TiC \cite{papazoglou1968diffusion, hatano2006diffusion}, which could be useful for the hydrogen insertion mechanism in Ti\textsubscript{2}AlC.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=1.0\linewidth]{fig/neb.png}
	\caption{Diffusion of H in Ti-Al layers in Ti\textsubscript{2}AlC (a) and Ti\textsubscript{2}CuC (b) structure.}
	\label{fig:neb}
\end{figure}

Figure \ref{fig:neb} (b) shows the H energy diffusion barrier for Ti\textsubscript{2}CuC. The path between I\textsubscript{tetr}-1 and I\textsubscript{tetr}-2 site has very negligible energy barrier. However, the opposite direction has a higher diffusion energy barrier of 0.62 eV. The diffusion energy barrier for I\textsubscript{tetr}-2 to I\textsubscript{oct}-3 is 0.37 eV, whereas the opposite path is easier with energy barrier of 0.10 eV. The path between I\textsubscript{oct}-3 and I\textsubscript{tetr}-1 is also easier than the path I. The energy barrier between I\textsubscript{oct}-3 and I\textsubscript{tetr}-1 is 0.28 eV and the opposite path is much easier with barrier energy 0.003 eV. Similar to Ti\textsubscript{2}AlC structure in Ti\textsubscript{2}CuC I\textsubscript{oct}-3 behave as an intermediate position. But in the Ti\textsubscript{2}CuC structure, the H atom at I\textsubscript{tetr}-1 directly move to I\textsubscript{tetr}-2 (instead of I\textsubscript{oct}-3) since the energy barrier is negligible.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=1.0\linewidth]{fig/ce.png}
	\caption{The agreement between the DFT recorded data and the cluster expansion (CE) model (a). The data points are clustered about the y = x line. Formation energy of Ti\textsubscript{2}Al\textsubscript{x}Cu\textsubscript{1-x}C at different Cu concentration (b). Both DFT (red) and CE (blue) formation energy are shown.}
	\label{fig:ce}
\end{figure}

The formation energy of Ti\textsubscript{2}CuC is higher than Ti\textsubscript{2}AlC. To find out the stable structures with copper substituted in Ti\textsubscript{2}AlC, cluster expansion (CE) Hamiltonian described in Eq.\ref{eq:3} was utilized to fit the DFT database of energies of Ti\textsubscript{2}Al\textsubscript{x}Cu\textsubscript{1-x}C alloys.  The accuracy of the CE model was verified by comparing it with DFT calculated energy of Cu doped Ti\textsubscript{2}AlC structure. The resulting CE energy is reasonably good at predicting DFT energy as the data points are clustered around the y=x line as shown in Fig. \ref{fig:ce} (a). The root-mean-square error (RMSE) of the fit was 0.0065 eV/atom. Using the fitted CE model, the formation energy of Cu doped Ti\textsubscript{2}AlC structures were calculated by increasing doped concentration by 5.56\%. Figure \ref{fig:ce} (b) shows the changes in formation energy in different Cu concentration obtained from both DFT calculation and CE. At 83.33\% doping concentration, the Ti\textsubscript{2}Al\textsubscript{x}Cu\textsubscript{1-x}C structure has the lowest formation energy, which is - 66.75 meV/unit cell (Fig. \ref{fig:h-alloy} (a)). The 83.33\% doped structure was relaxed by first-principles calculation before calculating hydrogen storage capacity and the dynamic stability of the structure was checked by phonon calculation. The 83.33\% Cu doped Ti\textsubscript{2}Al\textsubscript{x}Cu\textsubscript{1-x}C alloy structure is also dynamically stable in nature.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=1.0\linewidth]{fig/h-pure.png}
	\caption{Adsorption energy of hydrogen in Ti\textsubscript{2}AlC (a) and Ti\textsubscript{2}CuC (b) at different H concentration.}
	\label{fig:h-pure}
\end{figure}

To find out the hydrogen storage capacity of Ti\textsubscript{2}AlC, Ti\textsubscript{2}CuC, and TiAl\textsubscript{x}Cu\textsubscript{1-x}C, all interstitial sites in the Ti-A layer are filled increasingly with the hydrogen atom and calculate the adsorption energy. Figure \ref{fig:h-pure} shows the changes in adsorption energy with increasing H concentration in Ti-A layer in Ti\textsubscript{2}AlC and Ti\textsubscript{2}CuC structure. In Fig. \ref{fig:h-pure}, all tetrahedral and octahedral sites in Ti-A layers filled with hydrogen atoms is designated as a 1.0 fraction of hydrogen H in Ti-A layers and a 0.0 fraction of hydrogen H in Ti-A layers implies no hydrogen atom present in the structure. The hydrogen adsorption energy was calculated by increasing the hydrogen number continuously in all structures and at which highest fraction of hydrogen the adsorption energy negative was considered as the maximum hydrogen capacity of the structure. The hydrogen weight percentage was calculated at that fraction of hydrogen for each structure. On the basis of hydrogen adsorption energy, Ti\textsubscript{2}AlC has a maximum hydrogen storage capacity of 2.62 wt\%, which is higher than previously reported chemisorbed H capacity in MAXene \cite{doi:10.1021/jp409585v}. In the case of Ti\textsubscript{2}CuC structure, hydrogen storage capacity in Ti-Cu increases by 3.38 wt\%. Figure \ref{fig:h-alloy} (b) shows the changes in adsorption energy with increasing H concentration in Ti-Al-Cu layer in 83.33\% Cu doped Ti\textsubscript{2}Al\textsubscript{0.17}Cu\textsubscript{0.83}C structure. The Ti\textsubscript{2}Al\textsubscript{0.17}Cu\textsubscript{0.83}C has highest hydrogen storage capacity of 3.66 wt\% and Ti\textsubscript{2}Al\textsubscript{0.17}Cu\textsubscript{0.83}CH\textsubscript{5.72} phase. In the Ti\textsubscript{2}Al\textsubscript{0.17}Cu\textsubscript{0.83}CH\textsubscript{5.72} phase, all I\textsubscript{tetr}-2 sites are filled by hydrogen atoms and 94.44\% and 91.67\% of I\textsubscript{oct}-3 and I\textsubscript{tetr}-1 sites are occupied by hydrogen atoms, respectively.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=1.0\linewidth]{fig/h-alloy.png}
	\caption{Crystal structure 83.33\% Cu doped optimized Ti\textsubscript{2}Al\textsubscript{0.17}Cu\textsubscript{0.83}C structure (a). Blue, purple, black, and brown balls correspond to Ti, Al, C, and Cu, respectively. Adsorption energy of hydrogen in Ti\textsubscript{2}Al\textsubscript{0.17}Cu\textsubscript{0.83}C at different H concentration (b).}
	\label{fig:h-alloy}
\end{figure}

The thick black curve in Fig \ref{fig:rt} shows the variation of the chemical potential of hydrogen atoms with the temperature at 1 bar pressure. The value of $\mu^c_H$ for different MAX phases at 50\% storage capacity is plotted by horizontal lines. The intersection between the these lines and the hydrogen chemical potential curve indicates the corresponding critical temperature of hydrogen release from the particular surface at 1 bar pressure \cite{Williamson2004}. At 50\% storage capacity, Ti\textsubscript{2}AlC phase has a very low critical hydrogen chemical potential (-3.41 eV). Whereas, the hydrogen release temperature of Ti\textsubscript{2}CuC phase was 677 K, which is quite higher than the operating temperature of the hydrogen fuel cell \cite{Eberle2012}. The Ti\textsubscript{2}Al\textsubscript{0.17}Cu\textsubscript{0.83}C alloy structure has a 326 K hydrogen release temperature, which is ideal for hydrogen storage applications. The hydrogen release temperature of  the lowest energy structure, Ti\textsubscript{2}Al\textsubscript{0.28}Cu\textsubscript{0.72}C obtained  using CE method was also calculated. The Ti\textsubscript{2}Al\textsubscript{0.28}Cu\textsubscript{0.72}C structure has a hydrogen release temperature of 311K, which is lower than Ti\textsubscript{2}Al\textsubscript{0.17}Cu\textsubscript{0.83}C structure.

\begin{figure}[htbp]
	\centering
	\includegraphics[width=1.0\linewidth]{fig/rt.png}
	\caption{The temperature dependence of the chemical potential of hydrogen (black line). Horizontal lines marks the value of $\mu^C_H$ different MAX phase.}
	\label{fig:rt}
\end{figure}

\section{Conclusions}
In summary, we performed first-principles based cluster expansion calculation to study the hydrogen storage capacity in the MAX phase. The hydrogen atom is thermodynamically  stable in Ti-A layers of Ti\textsubscript{2}AC phase. Hydrogen atoms located in the Ti-C layer moved to the Ti-A layers. Among three interstitial sites in Ti-A layers, I\textsubscript{tetr}-2 site is most favorable for hydrogen storage. DFT calculation shows among different Ti\textsubscript{2}AC phases. Ti\textsubscript{2}CuC has the highest hydrogen adsorption energy due to the deeper p-orbital position of the Cu atom. Also, this phase is more preferable for the hydrogen release in comparison to Ti\textsubscript{2}AlC due to lower Coulomb interaction. The diffusion of hydrogen in Ti\textsubscript{2}CuC phase is more easier than Ti\textsubscript{2}AlC due to the lower energy barrier. The most stable Cu doped Ti\textsubscript{2}AlC alloy has been found by the cluster expansion model. The 83.33\% Cu doped Ti\textsubscript{2}Al\textsubscript{x}Cu\textsubscript{1-x}C alloy structure is both energetically and dynamically stable and can store 3.66 wt\% hydrogen at ambient atmospheric conditions, which is higher than both Ti\textsubscript{2}AlC and Ti\textsubscript{2}CuC phase. The 83.33\% Cu doped Ti\textsubscript{2}Al\textsubscript{x}Cu\textsubscript{1-x}C alloy structure is a promising 211 MAX phase for hydrogen storage application.

\section{Acknowledgement}
 This work was supported by the Korea Institute of Science and Technology (Grant number 2E31851), GKP (Global Knowledge Platform, Grant number 2V6760) project of the Ministry of Science, ICT and Future Planning.