\section{INTRODUCTION}

Tracing the evolution of abundance of galaxy clusters over cosmic time
is a promising approach to constrain cosmological parameters, such as the amplitude of mass density fluctuations, the equation of state of dark energy, and possible extensions to the standard $\Lambda$CDM model
\citep[e.g., see][for review]{2011ARA&A..49..409A}.
The current cosmological constraints derived from cluster abundance
hinge on the accuracy of mass estimation and the systematic uncertainty due to inaccurate mass estimates already dominates the statistical error in current cluster-based cosmological constraints 
\citep[e.g.][]{
2009ApJ...692.1060V, 
2013JCAP...07..008H, 
2014MNRAS.440.2077M, 
2016A&A...594A..24P, 
2016ApJ...832...95D}.

Gravitational lensing is a powerful technique for accurately measuring the cluster masses, because the lensing effect is purely induced by gravity, independent from the physical nature or the dynamical state of the gravitating mass.  Thanks to deep galaxy imaging surveys, it has become feasible to 
measure the mass density distribution around individual galaxy clusters
through gravitational lensing measurements \citep[e.g.][]{2008ApJ...685L...9B, 2008JCAP...08..006M, 2010PASJ...62..811O, 2012MNRAS.420.3213O, 2014ApJ...795..163U}. 

However, accurate cluster mass estimation requires controlling various systematic uncertainties in modelling the mass distribution of galaxy clusters. Baryonic processes, such as radiative cooling, star formation, and feedback from supernovae and active galactic nuclei, can modify the mass profiles of galaxy clusters \citep[e.g.,][]{2008ApJ...672...19R, 2010MNRAS.405.2161D,2014MNRAS.442.2641V,2015MNRAS.451.1247S, 2015MNRAS.452..343S}. The change in the mass distribution due to baryonic effects can be characterized by the change in the mass concentration of the cluster halo. This effect has been modelled using simple adiabatic contraction in response to the baryonic dissipation in the cluster center \citep[e.g.,][]{2011arXiv1108.5736G, 2012MNRAS.424.1244F}. However, these adiabatic contraction models did not take into the account the dependence of the cluster mass profiles (via halo mass concentration) on the mass assembly histories (MAH). 

The mass concentration of a dark matter (DM) halo has been shown to be correlated with its MAH \citep[e.g.,][]{2002ApJ...568...52W}. Recent work by \citet{2013MNRAS.432.1103L} showed that statistically the mass profile of a given DM halo is directly related to its MAH, through the inside-out halo growth: inner regions of the halo are formed earlier, while the outer regions are governed by late-time accretion. Disentangling the effects of MAH and baryonic effects on halo mass concentration are important for accurate weak lensing mass estimates. 

In this work, we extend the relation between mass concentration and MAH of cluster-size halo of \citet{2013MNRAS.432.1103L}, by incorporating the changes in mass concentration due to baryonic physics. We calibrate this model using the {\em Omega500} hydrodynamical cosmological simulations that include various baryonic physical processes. We then apply our model to stacked weak lensing measurements to constrain the effects of baryons on the mass concentration and weak lensing mass estimate. Our model provides a physically motivated framework for simultaneously constraining both baryonic physics and cluster masses with optical cluster surveys, such as the upcoming Large Synoptic Survey Telescope (LSST) survey. 

Our paper is organized as follows.
Section~\ref{sec:model} summarizes a theoretical model of the mass distribution and concentration of DM halo and its relation to the MAH, including the effects of baryons. A brief description of our hydrodynamical simulations is given in Section~\ref{sec:sim}.
In Section~\ref{sec:res}, we present the relation of \change{the} mass profile with MAH in the presence of baryonic processes and a detailed comparison with our model and simulations, and investigate baryonic effects on the halo concentration
for future weak lensing measurements of galaxy clusters. 
Our main results are summarized in Section~\ref{sec:conclusions}.

Throughout this paper, $\log$ represents the logarithm with base 10, while $\ln$ is the natural logarithm.


\section{Theoretical model} \label{sec:model}


\subsection{Mass profile}

To characterize the mass distribution of gravitationally bounded objects, we assume averaged mass density profile $\rho(r)$ is spherically symmetric and is only a function of the cluster-centric radius $r$. 
Numerical simulations have shown that spherically averaged density profiles of DM halos are well approximated by scaling a simple formula proposed in \citet{Navarro1996}.
The profile is referred as the NFW profile and it is given by
\begin{eqnarray}
\frac{\rho_{\rm NFW}(r)}{\rho_{\rm crit}}=\frac{\delta_{c}}{(r/r_s)(1+r/r_s)^2},
\label{eq:def_NFW}
\end{eqnarray}
where 
$\rho_{\rm crit}(z)$
is the critical density of the universe at redshift of $z$, 
$r_s$ is a scale radius, 
and
$\delta_c$ is a dimensionless characteristic density of the halo. 
These parameters $\delta_c$ and $r_s$ can be condensed into one parameter, 
the concentration $c_{\Delta}=r_{\Delta}/r_s$, through the definitions of the halo mass 
\begin{eqnarray}
M_{\Delta} &\equiv& \frac{4\pi}{3} \Delta \, \rho_{\rm crit}(z) r^3_{\Delta}, \label{eq:SOmass_def} \\
M_{\Delta} &=& \int_{0}^{r_{\Delta}}\, {\rm d}x\, 4\pi x^{2}\rho_{\rm NFW}(x), \label{eq:SOmass_NFW}
\end{eqnarray}
where $\Delta$ is the over-density parameter, and 
Eqs~(\ref{eq:SOmass_def}) and (\ref{eq:SOmass_NFW}) reduce to
\begin{eqnarray}
\delta_{c} = \frac{\Delta}{3}\frac{c^3_{\Delta}}{\left[\ln (1+c_{\Delta})-c_{\Delta}/(1+c_{\Delta})\right]}.
\end{eqnarray}
\change{In the following, we simply write the density concentration parameter as $c$ and 
use $c_{\Delta}$ as necessary.}

\change{For the NFW density profile with a given $\Delta$, we have the expression of the mass 
profile, denoted as $M(\langle \rho \rangle)$, in terms of the enclosed mean density, $\langle \rho \rangle$,}
by solving the following equations;
\begin{eqnarray}
\frac{M(x)}{M_{\Delta}} &=& \frac{Y(c x)}{Y(c)},
\,\,\,\,\,\,\,\,\,\,
\frac{\langle \rho \rangle(x)}{\Delta \rho_{\rm crit}} 
= \frac{1}{x^3}\frac{Y(cx)}{Y(c)}, \label{eq:NFW_scaled}
\end{eqnarray}
where 
$x=r/r_{\Delta}$ and $Y(u)=\ln(1+u)-u/(1+u)$.

\subsection{Relation between mass profile and accretion history}

Recently, \citet{2013MNRAS.432.1103L} 
showed that there exists a one-to-one relation 
between mass profile $M(\langle \rho \rangle)$ and the mass accretion history (MAH) of the halo.
The MAH is commonly defined as
\change{the spherical over-density mass}
of the main progenitor as a function of redshift, 
denoted as $M_{\rm prog}(z)$.
\change{Analogously to Eq.~(\ref{eq:NFW_scaled}), we assume the average MAH can be given by the following
parametric expression:}
\begin{eqnarray}
\frac{M_{\rm prog}(z)}{M_{\rm prog}(z_o)}
= \frac{Y(c_{\rm MAH}q)}{Y(c_{\rm MAH})},
\,\,\,\,\,\,\,\,\,\,
\frac{\rho_{\rm crit}(z)}{\rho_{\rm crit}(z_o)}
= \frac{1}{q^3}\frac{Y(c_{\rm MAH}q)}{Y(c_{\rm MAH})}, \label{eq:MAH_scaled}
\end{eqnarray}
where 
\change{$q$ is a parameter of this model} and $z_o$ is the redshift for galaxy clusters of interest.
The assumption of
Eq.~(\ref{eq:MAH_scaled}) works for DM halos 
with various masses and redshifts on average 
\change{in dark-matter-only simulations}
\citep{2013MNRAS.432.1103L}.
In addition, the ``concentration'' parameter $c_{\rm MAH}$ is found to be 
tightly correlated with the mass concentration of $c$.
This tight correlation between $c_{\rm MAH}$ and $c$ can be summarized as 
\begin{eqnarray}
\langle \rho \rangle (r_s) \propto \rho_{\rm crit}(z_{s}), \label{eq:CDM_prof_MAH_relation}
\end{eqnarray}
where a scaled redshift $z_{s}$ is given by
\begin{eqnarray}
M_{\rm prog}(z_s) = M(<r_s). \label{eq:scaled_redshift}
\end{eqnarray}

\subsection{Modelling baryonic effects on mass concentration}

To model the effects of baryons on the mass concentration, we generalize the relation of Eq.~(\ref{eq:CDM_prof_MAH_relation}) as 
\begin{eqnarray}
\log\left(\frac{\langle\rho\rangle (r_s)}{\rho_{o}}\right)
= \alpha_{0} + \alpha_{1} \log\left(\frac{\rho_{\rm crit} (z_s)}{\rho_{o}}\right),
\label{eq:general_mass_prof_MAH_relation}
\end{eqnarray}
where $\rho_{o}$ is the critical density of the universe at $z=z_{o}$.
Throughout this paper, we refer to
Eq.~(\ref{eq:general_mass_prof_MAH_relation}) as the mass-profile-MAH relation.

\change{
Substituting Eqs~(\ref{eq:NFW_scaled}) and (\ref{eq:MAH_scaled}) into Eq.~(\ref{eq:general_mass_prof_MAH_relation}),
one can find a unique relation 
between $c$ and $c_{\rm MAH}$, and it is independent of redshift $z_o$. 
According to \citet{2013MNRAS.432.1103L}, 
the relation can be approximately expressed as
\begin{eqnarray}
c = a_1 \left( 1+a_{2} c_{\rm MAH}\right)^{a_3},
\label{eq:cs_relation}
\end{eqnarray}
where $a_{1}, a_{2}$ and $a_3$ are constants, which can be evaluated for given $\alpha_0$ and $\alpha_1$.
}

\change{
In this paper, we assume that baryonic effects will change $\alpha_0$ and $\alpha_1$ in Eq.~(\ref{eq:general_mass_prof_MAH_relation}), but $c_{\rm MAH}$ is
determined primarily by gravitational structure formation alone and insensitive to baryonic effects.
In other words, the impact of baryonic effects on mass concentration can be controlled
by three parameters $a_{1}, a_{2}$ and $a_{3}$ in Eq.~(\ref{eq:cs_relation}),
while $c_{\rm MAH}$ can be responsible for mass-redshift-cosmology dependence.
We then follow the procedure developed in \citet{2014MNRAS.441..378L} to compute 
the mass concentration $c$ for arbitrary masses, redshifts, cosmologies given Eq.~(\ref{eq:cs_relation})\footnote{
We also summarize the procedure to compute 
the mass concentration $c$ in Appendix~\ref{app:mah_model}.}.
}

In the following section, we will examine these assumptions with high-resolution hydrodynamical cosmological
simulations with varying baryonic physics.


 




\section{SIMULATIONS}\label{sec:sim}


\subsection{The {\em Omega500} simulation}
In this work, we analyze the mass-limited sample of galaxy clusters extracted from the {\it Omega500} simulation series. 
The simulations assume 
flat $\Lambda$CDM model with the WMAP five-year results \citep{2009ApJS..180..330K}:
$\Omega_{\rm m0} = 0.27$ (matter density), $\Omega_{\rm b0} = 0.0469$ (baryon density),
\change{$H_{0} \equiv 100h = 70\, {\rm km}\, {\rm s}^{-1}{\rm Mpc}^{-1}$} (Hubble constant),
and $\sigma_{8} = 0.82$ (the mass variance within a sphere with a radius of 8 $h^{-1}\, {\rm Mpc}$).
The simulation is performed using the Adaptive Refinement Tree (ART)
$N$-body+gas-dynamics code 
\citep{1999PhDT........25K, 2002ApJ...571..563K, 2008ApJ...672...19R}, 
which is an Eulerian code that uses adaptive refinement in space and time and non-adaptive refinement in mass to achieve the dynamic range necessary to resolve the cores of halos formed 
in self-consistent cosmological simulations.  
The simulation volume has a comoving box length 
of 500~$h^{-1}\, {\rm Mpc}$, resolved using a uniform $512^3$ root grid and 8 levels of mesh refinement, implying a maximum comoving spatial resolution of 3.8~$h^{-1}\, {\rm kpc}$.

The {\it Omega500} simulation series consists of three runs under different baryonic physics;
non-radiative (NR) hydrodynamics as in \citet{2014ApJ...782..107N},
the run with additional baryonic physics, such as radiative cooling, star formation (CSF), 
and further including feedback from active galactic nuclei (AGN).
In the NR run, we have treated the ICM as a non-radiative gas and ignored additional baryonic physics.
The CSF run includes metallicity-dependent radiative cooling,
star formation, thermal supernova feedback, metal enrichment and advection, 
which are based on the same subgrid physics modules in 
\citet{2007ApJ...668....1N}, to which we refer the reader for more detail. 
The AGN run includes the CSF physics, mentioned above, 
plus a subgrid thermal AGN feedback module, similar to the one adopted in in \citet{booth_etal09}. 
First, supermassive blackholes (SMBH) are seeded as particles with an initial mass of $10^5h^{-1}M_\odot$ 
at the centres of DM halos with $M_{500c} > 2 \times 10^{11}h^{-1}M_\odot$. 
These SMBH grow via mergers and gas accretion with a rate given by a modified Bondi accretion model and return the feedback energy as a fraction of the accreted rest mass energy ($\epsilon = 0.2$) into the environment in the form of thermal energy.  
We also impose a minimum heating temperature, 
$T_{\rm min} = 10^7$~K, requiring that the SMBH store enough feedback energy 
until they accumulate enough energy to heat neighbouring gas cells, each by an amount of $T_{\rm min}$ to keep the injected thermal feedback energy from radiating away. 

Cluster-sized halos are identified in the simulation 
using a spherical overdensity halo 
finder described in \citet{2014ApJ...782..107N}. 
We define the three-dimensional mass of cluster-size halos using the spherical overdensity criterion with $\Delta=500$ in Eq.~(\ref{eq:SOmass_def}). We select DM halos with 
$M_{\rm 500c}\geq 3 \times 10^{14}\, h^{-1}M_{\odot}$ at 
$z=0$ and re-simulate the box with higher resolution DM particles 
in regions of the selected halos with the ``zoom-in'' technique \citep{2001ApJ...554..903K}, resulting in an effective mass resolution of $2048^{3}$, corresponding to a DM 
particle mass of  $1.09\times10^{9}\, h^{-1} M_{\odot}$, inside a spherical region with the cluster-centric radius of three times the virial radius for each halo. 
After the selection with $M_{\rm 500c}\geq 3 \times 10^{14}\, h^{-1}M_{\odot}$, we found 63, 80, 82 halos in the NR, CSF, and AGN runs, respectively.

\begin{figure}
\centering
\includegraphics[width=0.80\columnwidth, bb=0 0 511 550]
{./figure/mass_profile_diff.pdf}
\caption{
	The baryonic effect on mass density profile of 
	cluster-sized DM halos.
	The colored points show the average mass profile 
	for the mass-selected sample with 
	$M_{\rm 500c}\ge3\times10^{14}\, h^{-1}M_{\odot}$ 
	at $z=0$.
	Black, cyan and red points represent the results from the NR, CSF, and AGN runs.
	Residuals from the NR case are shown in the bottom panel.
	The dashed lines at the bottom highlight the level of $\pm0.1$ dex.
	}
\label{fig:profiles_sim}
\end{figure} 

Note that our CSF simulation suffers from the well-known 
``overcooling'' problem, where the simulation
over-predicts the amount of central stellar mass by a factor of $\sim 2$. As such, the results of our NR and CSF runs can be used to {\em bracket} uncertainties associated with baryonic effects.
The effects of baryonic physics, such as radiative gas cooling, 
star formation and energy feedback from supernovae and AGN are important in the cluster core regions.

Figure~\ref{fig:profiles_sim} highlights the impact of baryonic physics on the mass density profile around cluster-sized halos.
In the figure, we show the mass density profile averaged over 
the mass-selected sample in three {\it Omega500} runs.
We found the additional radiative effects can increase the mass density profile in the region of $r/r_{\rm 200c}<0.1$.
We also confirmed these contractions in the mass profile are 
consistent with the modified adiabatic contraction model
as in \citet{2011arXiv1108.5736G} with appropriate parameters
(see Appendix~\ref{app:mass_contraction_model} for details).

\subsection{Fitting Density Profiles and Mass Accretion Histories}

\subsubsection{Mass density profile as a function of radius}

In this section, we examine the relation of the enclosed mass profile $M(\langle \rho \rangle)$ and MAH \citep{2013MNRAS.432.1103L} in the presence of baryonic physics using the {\em Omega500} hydrodynamic cosmological simulations. 

We first obtain the best representation of
spherically averaged mass density profile $\rho(r)$ and then 
set a scale radius $r_s$ to define the relation as in Eq.~(\ref{eq:CDM_prof_MAH_relation}). 
Specifically, we compute the spherically averaged density profiles 
of cluster-sized halos in the {\em Omega500} simulations.
We then construct radial profile by $\log$-spaced binning 
in the range of $10 <r /(h^{-1} {\rm kpc})<10^4$ with $99$ bins.
To find the best representation of $\rho(r)$,
we fit using the non-linear 
least-squares Levenberg-Marquardt algorithm \citep{1992nrfa.book.....P}
by minimizing the $\chi^{2}$-fitting metric defined as 
\begin{eqnarray}
\chi^{2} = \sum_{i=1}^{N_{\rm bin}}
\left[\ln\rho(r_{i})-\ln\rho_{\rm NFW}(r_{i}|\, r_s, \, \delta_c)\right]^2,
\end{eqnarray}
where $N_{\rm bin}$ is the number of radial bins.
When performing $\chi^2$ fitting, we consider the radial range of 
$0.15<r/r_{200c}<1$, avoiding the core region in order to mitigate uncertainties associated 
with physical modelling of cluster core regions (e.g., see Figure~\ref{fig:profiles_sim}).

After finding the best-fitted parameter of $r_s$ for each halo,
we derive the scaled redshift $z_s$ from the information of MAH
and using Eq.~(\ref{eq:scaled_redshift}).

\subsubsection{Mass profile as a function of density}
Next, we consider the enclosed mass profile as a function of density.
For the matter density profile around cluster-sized halos at $z=0$, 
we introduce the expression of the enclosed mass profile as a function of the enclosed mass density, denoted as $M(\langle \rho \rangle)$.
Assuming the NFW profile, we fit the profile $M(\langle \rho \rangle)$
with single parameter $c_{\Delta}$ once we properly normalize
$M$ and $\langle \rho \rangle$ as shown in Eqs~(\ref{eq:NFW_scaled}).
To find the best representation of $M(\langle \rho \rangle)$,
we use a $\chi^{2}$-fitting metric defined as
\begin{eqnarray}
\chi^{2} = \sum_{i=1}^{N_{\rm bin}}
\left[\ln m(\hat{\rho}_{i})-\ln m_{\rm NFW}(\hat{\rho}_{i}|\, c_{\Delta})\right]^2,
\end{eqnarray}
where $m=M/M_{\Delta}$, $\hat{\rho}=\langle \rho \rangle/(\Delta \rho_{\rm crit}(z=0))$,
and $m_{\rm NFW}(\hat{\rho})$ 
is given by the solution of the parametric form given by 
Eq.~(\ref{eq:NFW_scaled}). 
\change{We use the logarithmic binning in the range of $1\le\hat{\rho}\le10^2$ 
with the bin size of $\Delta \log \hat{\rho} = 0.2$.}
When performing $\chi^2$ fitting on the mass profile, 
we adopt $\Delta=200$ and limit the density range of 
\change{$1 <\langle \rho \rangle/[200 \, \rho_{\rm crit}(z=0)] < 35$}, 
where the latter condition corresponds to the excision of cluster core regions with $r/r_{200c}\simlt0.15$.

\subsubsection{Mass Accretion History}
Finally, we measure MAH of individual halos as a function of redshift. Here, we define the cluster mass by $M_{\rm 500c}$. Using Eq.~(\ref{eq:MAH_scaled}), we perform a $\chi^2$ fitting to the MAH by minimizing the following metric as
\begin{eqnarray}
\chi^{2} = \sum_{i=1}^{N_{\rm bin}}
\left[\ln m_{\rm prog}(\hat{\rho}_{{\rm prog}, i})
-\ln m_{\rm prog, NFW}(\hat{\rho}_{{\rm prog}, i}|\, c_{\rm MAH})\right]^2,
\end{eqnarray}
where 
$m_{\rm prog} = M_{\rm prog}(z)/M_0$,
$\hat{\rho}_{\rm prog} = \langle \rho \rangle/\rho_{\rm crit}(z=0)$, and
$m_{\rm prog, NFW}(\hat{\rho}_{\rm prog})$ 
is given by the solution of Eq.~(\ref{eq:MAH_scaled}).
\change{We use the logarithmic binning in the range of $1\le\hat{\rho}_{\rm prog}\le10$ 
with the bin size of $\Delta \log \hat{\rho}_{\rm prog}= 0.2$.}


\section{Fitting the Mass-profile-MAH model}\label{sec:res}


\subsection{Impacts of baryonic physics on the Mass-profile-MAH relation}
\label{subsec:mass_prof_MAH}


\begin{figure}
\centering
\includegraphics[width=0.80\columnwidth, bb=0 0 519 471]{./figure/rho_MAH_scatter_plot_v2.pdf}
\caption{
	Relation between mass profiles at $z=0$ 
	and accretion histories for mass-limited haloes 
	($M_{\rm 500c}\ge3\times10^{14}\, h^{-1}M_{\odot}$)
	in our simulations with various baryonic physics.
    \change{
    The vertical axis shows the mean enclosed densities within the scale radius, 
    $r_{s}$, of the NFW profile. 
    The horizontal axis is for the critical density of the universe at the time 
	when the main progenitor's mass equals the present mass enclosed within $r_s$.
	Each point represents densities normalized 
	by the critical density in the universe at $z=0$, denoted as $\rho_{o}$.
    }

   



	Black, cyan, and red points indicate individual halos from the NR, CSF, and AGN simulations, respectively, while the colored dashed lines are the best-fitted scaling relations for each run.
}
\label{fig:rho_MAH}
\end{figure} 

First, we examine the correlation between the enclosed mass density $\langle \rho \rangle$ and the critical density $\rho_{\rm crit}$ in MAH for individual halos. Figure~\ref{fig:rho_MAH} shows the scatter plot of $\langle \rho \rangle$ at $r=r_s$ and the critical density at $z=z_s$, 
where we define $z_s$ as in Eq.~(\ref{eq:scaled_redshift}).
In this figure, the different coloured points represent the results for three different physics runs (NR, CSF, AGN) in {\em Omega500}.
Performing a linear least-square fitting to 63 clusters in
our NR sample at $z=0$, the best-fit normalization and slope for the scaling relation between
$\log \langle \rho \rangle (r_s)$ and $\log \rho_{\rm crit}(z_s)$ (Eq.~\ref{eq:general_mass_prof_MAH_relation}) is
$(\alpha_0, \alpha_1) = (3.057\pm 0.289, 0.996\pm0.528)$, where the error values
indicate the $1\sigma$ errors. The best-fit slope for NR haloes is consistent with 
the results obtained using the dark-matter-only simulations by \citet{2013MNRAS.432.1103L}.

Once including baryonic effects associated with gas cooling and star formation in the CSF run, we find $\langle \rho \rangle (r_s)$ tends to become higher for a given $\rho_{\rm crit}$.
The best-fit parameters for our CSF sample is $(\alpha_0, \alpha_1) = (3.063\pm 0.354, 1.513\pm0.528)$, which indicates marginal $1\sigma$ deviation from the DM-only linear relation given by the Eq.~(\ref{eq:CDM_prof_MAH_relation}).
For the AGN run, the best-fit parameters are $(\alpha_0, \alpha_1) = (3.036\pm 0.320, 1.212\pm0.504)$.
Although the mass-profile-MAH relation for individual haloes
has a large scatter, as shown in Figure~\ref{fig:rho_MAH},
we find that the mean mass-profile-MAH relation is sensitive to baryonic physics, as noted by the significant change in slope of the mean relation, $\alpha_{1}$ in Eq~(\ref{eq:general_mass_prof_MAH_relation}).

We also quantify the intrinsic scatter of the $\langle \rho \rangle (r_s)-\rho_{\rm crit}(z_s)$ relation as
\begin{eqnarray}
\sigma^2_{\rm int} = \frac{1}{N-1}\sum_{i=1}^{N}
\left[\log \langle \rho \rangle_{i} - \log \langle \rho \rangle_{\rm fit}
\left(\rho_{\rm crit, i}\right)\right]^2,
\end{eqnarray}
where $N$ is the number of halos and $\langle \rho \rangle_{\rm fit}$
is the best-fit relation for each run. The best-fit parameters and scatters are summarized in Table~\ref{tab:alphas}.

\begin{figure*}
\centering
\includegraphics[width=1.5\columnwidth, bb=0 0 570 416]
{./figure/mass_profile_MAH_comp_revise_rhole35.pdf}
\caption{
 	Mass profile as a function of mean enclosed density
	and accretion history as a function of the critical 
    density of the universe. 
    Top left panels show the mass profile, 
    while the MAH is presented in top right.
	We consider the mass-limited sample with 
    $M_{\rm 500c}\ge3\times10^{14}\, h^{-1}M_{\odot}$
	at $z=0$ for three different simulations.
    In each panel, black, cyan, and red point or line are the results of the NR, CSF, and AGN runs, respectively.
    In both top panels, the points show the average values 
    over the sample with different baryonic physics runs.
    Solid lines in the top left panel show 
    the NFW profiles fitted to the average mass profiles for each run,
	while the dashed lines in the right panel 
    represent MAH predicted
	by the relation between mass profiles and accretion histories as in
    Figure~\ref{fig:rho_MAH}.
    For comparison, the cyan and red dotted lines in the top right panel
    show the model predicted by 
    Eq~\ref{eq:general_mass_prof_MAH_relation} with $\alpha_1=1$.
	The bottom left panel shows the residual of the measured mass profile 
    from best-fit NFW profile.
	In the bottom right panel, we show the difference between 
    the measured and predicted MAH.
    Note that the thin red lines in both top panels show 
    the result for individual halos in the AGN run, illustrating
    the typical scatter in the mass profiles.
	}
\label{fig:mass_MAH_comp}
\end{figure*} 

\begin{table*}
\centering
\caption{
	Relation of concentration parameters 
	in mass profile and accretion history
	based on the mass-profile-MAH relations 
	$ (\alpha_{0},\alpha_1,a_1,a_2a_3)$
	and the intrinsic scatter around the best fit 
	$\sigma_{\rm int}$,
	measured in the {\it Omega500} simulations.
    Note that the parameters are valid for $\Delta=200$ in mass profile and when MAH is defined as the spherical over-density mass with $\Delta=500$.
	}
\begin{tabular}{@{}lcccccccl}
\hline
\hline
Run & $\alpha_{0}$ & $\alpha_1$ & $a_1$ & $a_2$ & $a_3$  & $\sigma_{\rm int}$ \\ \hline 
NR  & $3.057\pm 0.289$ & $0.995\pm 0.528$ & $2.840$ & $0.646$ & $1.010$  & $0.206$ \\ \hline 
CSF & $3.063\pm 0.354$ & $1.513\pm 0.528$ & $3.342$ & $0.836$ & $1.438$  & $0.217$ \\ \hline 
AGN & $3.036\pm 0.320$ & $1.212\pm 0.504$ & $2.864$ & $0.764$ & $1.190$ & $0.209$ \\ \hline 
\end{tabular}\label{tab:alphas}
\end{table*}

\change{
Since $\langle \rho \rangle/\rho_o \propto c^3 [\ln(1+c)-c/(1+c)]^{-1}$ for the NFW profile,
it is expected that there exists a scaling relation between $c$ and $\rho_{\rm crit}(z_s)$.
We confirmed this expectation with {\it Omega500} simulations.
The scaling relation can be approximated as
\begin{eqnarray}
\log c
= \beta_{0} + \beta_{1} \log\left(\frac{\rho_{\rm crit} (z_s)}{\rho_{o}}\right), \label{eq:c_MAH_relation}
\end{eqnarray}
where $(\beta_0, \beta_1) = (0.402\pm 0.147, 0.467\pm 0.279)$ for the NR run,
$(0.405\pm 0.164, 0.651\pm 0.243)$ for the CSF run, 
and $(0.384\pm 0.165, 0.572\pm 0.261)$ for the AGN run.
Here the error values indicate the $1\sigma$ errors and $c$ is defined by $r_{200c}/r_s$ at $z=0$.
The intrinsic scatter around Eq.~(\ref{eq:c_MAH_relation}) in base 10 logarithmic space
is found to be 0.106, 0.103, 0.107 for the NR, CSF, and AGN runs, respectively.
}

\subsection{Testing the mass-profile-MAH relation model with hydrodynamical simulations}
\label{sec:detailed_test_relation}

Next, we test the accuracy of our mass-profile-MAH relation model by checking whether the model can reproduce the average MAH from the mass profile. To do this, we first measure the average mass profile for our mass-limited sample of simulated clusters and fit the average mass profile $M(\langle \rho \rangle)$ with the NFW profile. We then infer the MAH by using this average mass profile and the mass-profile-MAH relation given in Eq.~(\ref{eq:general_mass_prof_MAH_relation}), and compare the result with the true MAH directly measured from the simulation. 

Figure~\ref{fig:mass_MAH_comp} shows the comparison of the average MAH (left panel) and the predicted MAH from the best-fit NFW profile (right panel) for the mass-limited sample with $M_{\rm 500c}\ge3\times10^{14}\, h^{-1}M_{\odot}$ in the {\it Omega500} runs with different baryonic physics. 
We find that the predicted MAH from mass profile is in reasonable agreement with the average MAH for all three runs. 
The accuracy of our model is at a level of 0.1 dex in the range of $1< \rho_{\rm crit}/\rho_{o} < 10$, corresponding to $0<z\simlt2.5$.

We also quantify the importance of baryonic physics in estimating the MAH from mass profile. In Figure~\ref{fig:mass_MAH_comp}, the differences between dashed and dotted lines in the top right highlight the importance of including the baryonic effects on the $\langle \rho \rangle (r_s)-\rho_{\rm crit}(z_s)$ relation.
If applying a simple relation of $\langle \rho \rangle (r_s) \propto \rho_{\rm crit}(z_s)$ to CSF/AGN runs, we find an estimate of the MAH at $\rho_{\rm crit}/\rho_{o}\simeq3$ or $z\sim1$ will be biased at the level of $\sim0.4$ and 0.2 dex for CSF and AGN runs, respectively. 

As observed in the three different physics runs in the {\em Omega500} simulation, baryonic effects can modify the concentration in mass density profile without affecting the overall MAH represented by $\rho_{\rm crit}(z_s)$. This can be understood by the relative importance of gravitational physics versus baryonic physics during different phases of the MAH of the halo. In the early phase of the MAH, where the cluster halo is rapidly accreting its mass, baryonic effects are expected to be sub-dominant to gravitational physics. This means that $\rho_{\rm crit}(z_s)$ in the mass-profile-MAH relation, which represents a typical density of the halo at the early formation phase, is relatively insensitive to the details of baryonic processes. At the later stage, where the halo is slowly accreting its mass, most of the accreted mass is deposited onto the outer regions of the halo. Baryonic physics dominates in the relatively quiescent cluster core region, rendering the average density within the scale radius, $\langle \rho \rangle (r_s)$, to become more
susceptible to the baryonic effects. 
As shown in Figure~\ref{fig:profiles_sim}, 
$\langle \rho \rangle (r_s)$ increases for runs with larger baryon dissipation. 
Therefore, baryonic effects increase $\langle \rho \rangle (r_s)$ without changing $\rho_{\rm crit}(z_s)$\footnote{
\change{Our simulations indicate that baryonic effects change $\langle \rho \rangle (r_s)$ and $r_s$ so that
Eq~(\ref{eq:scaled_redshift}) remains valid.}
}.

\if0

\subsection{From mass profile to MAH}
\label{subsec:From_prof_to_MAH}

We examine how accurately we can estimate the MAH on average from the mass-profile information.
We take the following steps to study the average mass-profile-MAH relation: (i) measure the average mass profile over mass-limited sample; (ii) fit the average mass profile $M(\langle \rho \rangle)$ with the NFW profile; (iii) estimate the correspond MAH inferred from average mass profile with the mass-profile-MAH relation (Eq.~\ref{eq:general_mass_prof_MAH_relation}); (iv) compare the estimated MAH with the average MAH in simulated halos.

In step (iii), we can have a NFW expression of 
$M_{\rm prog}$ as a function of $\rho_{\rm crit}(z)$ 
for a given $\rho_{\rm crit} (z_s)$ 
so that $M_{\rm prog}(z)$ will set to be $M_{\rm 500c}$ at $z=0$ and $M(<r_s)$ at $z=z_s$.
The above procedures clarify the existence of 
underlying mass-profile-MAH relation,
since the number of simulated clusters ranges from 60 to $\sim80$
and the average mass profile will be determined with high precision. \comment{DN: this paragraph is a bit hard to follow and could benefit from streamlining.}

\begin{figure}
\centering
\includegraphics[width=0.9\columnwidth, bb=0 0 570 416]
{./figure/mass_profile_MAH_NR.pdf}
\caption{
	Mass profile as a function of mean enclosed density
	and accretion history as a function of the critical density of the universe.
	Top left panels show the mass profile, while the MAH is presented in
	top right.
	We here consider the mass-limited sample with 
	$M_{\rm 500c}\ge3\times10^{14}\, h^{-1}M_{\odot}$
	at $z=0$ in the NR simulation.
	In both top panels, gray lines represent the results for individual halos,
	while red points show the average values over the sample.
	The red solid line in the left panel shows the NFW profiles fitted to the average mass profiles,
	while the red dashed line in the right panel represent a MAH predicted
	by the relation between mass profiles and accretion histories as in Figure~\ref{fig:rho_MAH}.
	Left bottom panel shows the residual of the measured mass profile from best-fit NFW profile.
	In the right bottom panel, we show the difference between the measured and predicted MAH.
	\comment{EL: combine this figure with Figure~\ref{fig:mass_MAH_CSF_AGN} into one figure (still with the same number of panels), and showing the median values for the 3 runs with thick lines with different colors, and showing individiual MAHs/profiles for the AGN run only with thin lines. }
	}
\label{fig:mass_MAH_NR}
\end{figure} 

\begin{figure}
\centering
\includegraphics[width=0.9\columnwidth, bb=0 0 570 416]
{./figure/mass_profile_MAH_CSF.pdf}
\includegraphics[width=0.9\columnwidth, bb=0 0 570 416]
{./figure/mass_profile_MAH_AGN.pdf}
\caption{
	Similar to Figure~\ref{fig:mass_MAH_NR}, but for different baryonic effects.
	The top panels show the case in radiative run, 
	while the bottom panels are for AGN run.
	In top right panels, green dashed line shows the predicted MAH 
	by the mass-profile-MAH
	relation inferred from NR run (i.e. Eq.~\ref{eq:mean_rho_MAH_NR}), 
	while red dashed line is for appropriate relations as in 
	Eqs~(\ref{eq:mean_rho_MAH_CSF}) 
	and Eqs~(\ref{eq:mean_rho_MAH_AGN}).
	}
\label{fig:mass_MAH_CSF_AGN}
\end{figure} 

Figure~\ref{fig:mass_MAH_NR} summarizes the above procedures in
NR run.
The left top panel shows the mass profile over mass-limited sample 
with $M_{\rm 500c}\ge3\times10^{14}\, h^{-1}M_{\odot}$.
Gray lines represent the mass profiles for individual halos, 
while red points are for the average profile.
After performing $\chi^2$ fitting with the NFW profile, we obtain
the best-fit NFW profile as shown in the red line in the left top.
The left bottom panel presents the residual of the average mass profile
with the best-fit NFW profile, showing the best-fit profile gives
a good representation of the average mass profile in simulation.
In the right panels, we summarize the analyses of MAH.
Gray lines in the right top panel show the MAH for each halo,
while the red points are for the average MAH over the mass-limited sample.
The red dashed line represents a predicted MAH inferred from the best-fit NFW profile as shown in the left and 
the best-fit Eq.~(\ref{eq:general_mass_prof_MAH_relation}).
The right bottom panels show the difference between the average MAH and prediction in our model. \comment{DN: can trim some of these descriptions of figures to minimize duplication with the figure caption.}
We find that the predicted MAH from mass profile
is in reasonable agreement with the average MAH in NR run.
The accuracy of our model is a level of 0.1 dex in the range of $1< \rho_{\rm crit}/\rho_{o} < 10$, corresponding to $0<z\simlt2.5$.

We then quantify how important complex baryonic physics are 
in estimation of the MAH from mass profile.
Figure~\ref{fig:mass_MAH_CSF_AGN}
shows the similar analysis to in Figure~\ref{fig:mass_MAH_NR},
but we consider the different {\it Omega500} runs, named CSF and AGN.
In the top portion, we show the results of the CSF run, while
the bottom is for the AGN run.
The legend in figure~\ref{fig:mass_MAH_CSF_AGN} is same as
Figure~\ref{fig:mass_MAH_NR}, but the green dashed line in the top left
panel shows the predicted MAH based on the calibrated relation for NR run. 
The differences between green dashed line and red points in the top right
highlight the importance of including the baryonic effects on 
the $\langle \rho \rangle (r_s)-\rho_{\rm crit}(z_s)$ relation.
If applying a simple relation of $\langle \rho \rangle (r_s) \propto
\rho_{\rm crit}(z_s)$ to CSF/AGN runs,
we find an estimate of the MAH 
at $\rho_{\rm crit}/\rho_{o}\simeq3$ 
or $z\sim1$ will be biased with a level of $\sim0.4$ 
and 0.2 dex for CSF and AGN runs, respectively.
Our model still reproduces the average MAH 
in presence of baryonic effects, in our CSF and AGN runs.
Therefore, the difference of the slope in mass-profile-MAH relations
for various baryonic effects cannot be explained by artificial issues.

\fi

\subsection{Predicting halo concentration}
\label{subsec:From_MAH_to_prof}

\if0
\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth, bb=0 0 507 550]
{./figure/MAH_diff.pdf}
\caption{
	The baryonic effect on MAH of 
	cluster-sized DM halos.
	The colored points show the average mass profile 
	for the mass-selected sample with 
	$M_{\rm 500c}\ge3\times10^{14}\, h^{-1}M_{\odot}$ 
	at $z=0$, whereas the error bars are for the standard deviations.
	Black, cyan, and red points represent the results of NR, CSF, and AGN runs, respectively.
	Gray line shows an NFW prediction with concentration of 0.5.
	Residuals from the NFW prediction are shown in the bottom panel.
	}
\label{fig:MAH_sims}
\end{figure} 
\fi

\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth, bb=0 0 512 475]
{./figure/cprof_Mass.pdf}
\caption{
	Mass dependence in concentration parameters of mass profile. 
    The small, open, colored points are the best-fitted concentration parameters
	from NFW fits to mass profiles and the average concentrations 
	are shown in the large filled points.
	The different curves represent our model assuming
	the mass-profile-MAH relation and the universal MAH proposed by \citet{2002MNRAS.331...98V}.
	The baryonic effect on halo concentration in simulated mass profiles
	can be explained primarily by the difference in the mass-profile-MAH relations.
	}
\label{fig:cMAH_cprof_Mass}
\end{figure} 

\begin{figure}
\centering
\includegraphics[width=0.8\columnwidth, bb=0 0 509 475]
{./figure/c_M_z0_34_comp_U16_v3.pdf}
\caption{
	Comparison of concentration of cluster-sized DM halos
	with our model and the observational constraint from 
	Hubble Space Telescope and Subaru Telescope
	\citep{2016ApJ...821..116U}.
	Different colored lines show the baryonic effect 
	on the halo concentrations, while the solid, dashed, and dotted lines
	reflect the dependence on cosmological parameters.
	The solid and dashed lines are the best-fit cosmological models 
    consistent with WMAP \citep{Hinshaw2013} and Planck \citep{2016A&A...594A..13P}, respectively.
	The dotted lines correspond to the cosmological model that
	is in agreement with 
	recent cosmic shear measurement \citep{2017MNRAS.465.1454H}. 
	The gray hatched region shows the $2\sigma$ constraints
	obtained from the combined strong and weak lensing analyses 
	of X-ray selected clusters
	\citep[see][for details]{2016ApJ...821..116U}. 
	We consider the redshift of 0.34 and $h=0.7 h_{70}$
	for the comparison. 
	}
\label{fig:cprof_Mass_comp_obs_model}
\end{figure} 

Armed with the well-tested mass-profile-MAH relation, we can predict the concentration of galaxy clusters by taking into account the effects of baryon physics on their total mass profiles. Figure~\ref{fig:cMAH_cprof_Mass} shows the comparison between $c_{200c}$ predicted by our model and those measured directly in simulations, demonstrating that our model can explain a trend of increasing halo concentration in the runs with enhanced baryon dissipation. Specifically, for the mass-limited sample with $M_{\rm 500c}\ge3\times10^{14}\, h^{-1}M_{\odot}$ at $z=0$, our model predicts the concentration parameter of 3.77, 5.52, and 4.21 for NR, CSF, and AGN runs, respectively.
These values are consistent with the halo concentrations measured directly by fitting NFW profiles to the $\it Omega500$ simulations with $1\sigma$ uncertainty as shown in Figure~\ref{fig:cMAH_cprof_Mass}.
Note also that our model assumes that the baryonic effects on the average MAH of cluster-sized DM halos are negligible, and we confirmed that this assumption is tenable in cluster-sized halos (see the right panel in Figure~\ref{fig:mass_MAH_comp}). The mass-profile-MAH relation in Section~\ref{subsec:mass_prof_MAH} plays a central role in explaining 
\change{the sensitivity of the concentrations found in hydrodynamical simulations to the baryonic physics implemented.}

Figure~\ref{fig:cprof_Mass_comp_obs_model}
summarizes the predictions of our theoretical model of 
$c_{\Delta}-M_{\Delta}$ relation at $z=0.34$.
Here, we consider three representative cosmological models
that are consistent with recent CMB measurements by WMAP
\citep{Hinshaw2013}
and Planck
\citep{2016A&A...594A..13P},
and cosmic shear measurement by Kilo-Degree Survey
\citep{2017MNRAS.465.1454H}. 
We consider a different set of $\Omega_{\rm m0}$
and $\sigma_{8}$ as indicated in the legend.
The different colored lines are the results based on the mass-profile-MAH
relations calibrated by a set of three simulations with different baryonic physics. 
We find that $\sim10\%$ difference of $\Omega_{\rm m0}$ is unimportant for determining the $c_{\Delta}-M_{\Delta}$ relation,
while $\sim4\%$ difference of $\sigma_8$ is more prominent.
Instead, the impact of baryonic physics on $c_{\Delta}$
is found to be critical under the concordance $\Lambda$CDM cosmologies.
Hence, the precise measurement of $c_{\Delta}-M_{\Delta}$ relation
will play an essential role \change{in discerning how baryons affect the properties of galaxy clusters}.
For comparison, the hatched gray region in the figure 
also shows recent observational constraint 
of the $c_{\Delta}-M_{\Delta}$ relation for X-ray selected high mass clusters
\citep{2016ApJ...821..116U}.
Although the current constraint is in agreement 
with our models that include three different baryonic physics,
Figure~\ref{fig:cprof_Mass_comp_obs_model} demonstrates
a possibility of distinguishing these models 
if the measurement of $c_{\Delta}$
at lower mass scales is performed.


\section{Constraining Baryonic Effects with Stacked Weak Lensing Measurements}
\label{sec:application}


Our theoretical framework for predicting the mass profiles of gravitationally bound objects can be applied to various cosmological analyses. In this section, we consider how future stacked lensing analysis of galaxy clusters can constrain the effects of baryons on the concentration and mass profiles of galaxy clusters. 

\begin{figure*}
\centering
\includegraphics[width=1.2\columnwidth, bb=0 0 519 501]
{./figure/fisher_const_table.pdf}
\caption{
	The expected parameter 
    constraints on stacked 
    lensing signals in ongoing and
    future imaging surveys. 
    Red circles show 
    a 68\% confidence level 
    in parameter spaces of 
    interest for ongoing 
    Subaru HSC-like survey,
    while green ones are 
    for future surveys like LSST.
    For the HSC-like survey, 
    we assume the sky coverage of 
    $1400\, {\rm deg}^{2}$ and 
    the source number density of 
    $30\, {\rm arcmin}^{-2}$.
    For the LSST-like survey,
    we assume 
    the survey area of 
    $20000\, {\rm deg}^{2}$ 
    with the source number density 
    of $50\, {\rm arcmin}^{-2}$.
    In both cases, 
    we consider a cluster sample with 
    mass of $M_{\rm 200m}=10^{14.5}\,h^{-1} M_{\odot}$ 
    at redshift of $0.3$, 
    the source redshift of 1,
    and the intrinsic scatter 
    in galaxy shapes of 0.3.
    }
\label{fig:fisher_circles}
\end{figure*}

\subsection{Modelling stacked weak lensing signal}
\label{subsec:stack_lens_exp}

Stacked lensing analysis of the large-scale structure of the universe is commonly defined by the cross correlation with the shape of background galaxies and the position of foreground objects. This statistical quantity provides a unique method 
of measuring the average mass distribution of foreground objects.
The observable in stacked lensing analysis is 
azimuthally averaged profile of tangential shapes of background galaxies
with respect to the foreground objects denoted as 
$\langle \gamma_{+} \rangle (\theta)$.
Suppose that the intrinsic shape of source galaxies is randomly oriented,
the expected value of $\langle \gamma_{+} \rangle (\theta)$
can be decomposed into three parts:
\begin{eqnarray}
\langle \gamma_{+} \rangle
= f_{\rm cen} \langle \gamma_{+} \rangle_{\rm cen} + (1-f_{\rm cen})\langle \gamma_{+} \rangle_{\rm off} + \langle \gamma_{+} \rangle_{\rm 2halo},
\label{eq:shear_profile}
\end{eqnarray}
where the first and second terms on the right hand side of Eq~(\ref{eq:shear_profile})
express the contribution from single dark-matter haloes hosting foreground objects,
while the third term arises from the clustering of neighbouring haloes.
The difference between $\langle \gamma_{+} \rangle_{\rm cen}$
and $\langle \gamma_{+} \rangle_{\rm off}$ comes from 
the imperfect knowledge of centres of host DM haloes.

In optical surveys, the position of the brightest cluster galaxy (BCG) is assumed to be
the centre of the host halo, while some observational studies have shown
that BCGs do not always sit in the centres of haloes 
\citep[e.g.,][]{2010MNRAS.405.2215O, 2012MNRAS.426.2944Z}.
In this paper, we take into account the mis-centering effect of proxy of the halo centre
by introducing \change{the distribution function of the offset between the halo centre and BCG position} as
\begin{eqnarray}
P(R_{\rm off}) = \frac{R_{\rm off}}{\sigma^2_{\rm off}}\exp\left(-\frac{R^{2}_{\rm off}}{2\sigma^2_{\rm off}}\right),
\label{eq:P_miscenter}
\end{eqnarray}
where 
$R_{\rm off}$ represents the offset in comoving length scale and $\sigma_{\rm off}$ is the scatter in the probability. 
\change{
In Eq.~(\ref{eq:shear_profile}), the factor $f_{\rm cen}$ represents the fraction of objects located at the true halo centre.}

To predict the first and second terms in Eq~(\ref{eq:shear_profile}),
we assume a smoothly truncated NFW profile as proposed in \citet{2009JCAP...01..015B},
\begin{eqnarray}
\rho_{\rm BMO}(r) = \rho_{\rm NFW}(r)
\left\{\frac{\tau^2_{\rm 200m}}{\left(r/r_{\rm 200m}\right)^2+\tau^2_{\rm 200m}}\right\}^{2}, \label{eq:BMO}
\end{eqnarray}
where $r_{\rm 200m}$ is the spherical over-density radius with respect to
$200$ times mean matter density in the Universe and $\tau_{\rm 200m}$ is the characteristic halo truncation radius normalized by $r_{\rm 200m}$, which depends on halo mass and cosmology \citep[e.g.][]{2015MNRAS.449.4264G}.
\citet{2011MNRAS.414.1851O} found that Eq~(\ref{eq:BMO}) can provide 
better agreement with simulated stacked signals compared to simple NFW profile. 
Using these ansatzes, we compute the corresponding lensing signal of $\langle \gamma_{+} \rangle_{\rm cen}$,
while the off-centering term is given by \citep[see][for the derivation]{2011PhRvD..83b3008O} 
\begin{eqnarray}
\langle \gamma_{+} \rangle_{\rm off}(\theta)
= \int \frac{\ell {\rm d}\ell}{2\pi} \, 
\tilde{\kappa}_{\rm BMO}(\ell)\exp\left(-\frac{\sigma^2_{\rm off}\ell^2}{2(1+z_{l})^2 D^2_{l}}\right) J_{2}(\ell \theta),
\end{eqnarray}
where 
$\tilde{\kappa}_{\rm BMO}$ is the Fourier component of lensing convergence 
from Eq~(\ref{eq:BMO}),
$z_{l}$ is the redshift of foreground objects,
$D_{l}$ represents the angular diameter distance to foreground objects,
and $J_{2}(x)$ is the second-order Bessel function.

For the third term in Eq~(\ref{eq:shear_profile}), 
we apply the halo-based approach as in \citet{2002PhR...372....1C}
and use the following expression \citep[see also e.g., ][]{2011MNRAS.414.1851O}
\begin{eqnarray}
\langle \gamma_{+} \rangle_{\rm 2halo}(\theta)= 
\frac{b_{h} \Omega_{\rm m0}}{\Sigma_{\rm crit} D^2_{l}}
\int \frac{\ell {\rm d}\ell}{2\pi}\, 
P_{\rm lin}\left(\frac{\ell}{(1+z_{l})D_{l}}, z_{l} \right) J_{2}(\ell \theta),
\end{eqnarray}
where $b_{h}$ is the linear halo bias,
$P_{\rm lin}(k, z)$ is the linear matter power spectrum 
for the wavenumber of $k$ at redshift $z$,
and $\Sigma_{\rm crit}$ is the critical mass density for lensing.
In this paper, we adopt the model of $b_{h}$ developed by \citet{2010ApJ...724..878T}.

Hence, the expected signal in stacked lensing analysis contains
the information about the mass density profile on scales smaller than virial radius of halos,
while the cosmological information can be extracted from large-scale signals.
Within our framework,
the relevant parameters in Eq~(\ref{eq:shear_profile}) are: (i) the halo mass ($M_{\rm 200m}$); (ii) the parameters in the mass-profile-MAH relation ($\alpha_{0}$ and $\alpha_{1}$ in Eq~(\ref{eq:general_mass_prof_MAH_relation})); (iii) the parameters for mis-centering probability ($f_{\rm cen}$ and $\sigma_{\rm off}$ in Eq~(\ref{eq:P_miscenter})); (iv) cosmological parameters, especially $\sigma_{8}$.

\subsection{Detectability of baryonic mass contraction in galaxy clusters}
\label{subsec:detect}

Using the theoretical model developed in the previous Section~\ref{subsec:stack_lens_exp}, we forecast the expected observational constraints on the mass-profile-MAH relation with Fisher matrix analysis, which provides an estimate of parameter uncertainties for a given statistical measurement. 

We first introduce the Fisher matrix of stacked lensing signals as 
\begin{eqnarray}
F_{\alpha \beta} &=& 
\frac{1}{2}{\rm Tr}\left[ A_{\alpha}A_{\beta} + C^{-1} H_{\alpha \beta}\right],
\label{eq:Fisher}
\\
A_{\alpha} &=& C^{-1}\partial C/\partial p_{\alpha},
\\
H_{\alpha \beta} &=& 2 (\partial \langle \gamma_{+} \rangle/\partial p_{\alpha})
(\partial \langle \gamma_{+} \rangle/\partial p_{\beta}),
\end{eqnarray}
where $\langle \gamma_{+} \rangle$ is the azimuthally averaged profile of tangential shapes of the background lensed galaxies defined in Eq.~(\ref{eq:shear_profile}) that takes into account the off-centering and the two-halo term contributions,
$C$ is the covariance of $\langle \gamma_{+} \rangle$, 
and $p_{\alpha}$ describes the parameter of interest. 
We consider the following parameters:
$\bd{p}=\{\log (M_{\rm 200m} \, [h^{-1}M_{\odot}]), 
\alpha_{1}, f_{\rm cen}, \sigma_{\rm off} \, [h^{-1}{\rm Mpc}], \sigma_{8}\}$.
Here $\alpha_1$ is the slope in the mass-profile-MAH relation in Eq.~(\ref{eq:general_mass_prof_MAH_relation}), $f_{\rm cen}$ and $\sigma_{\rm off}$ characterize the mis-centering of haloes: $f_{\rm cen}$ is the fraction of haloes with centres correctly located in observation (see Section~\ref{subsec:stack_lens_exp}), and $\sigma_{\rm off}$ is the scatter in the mis-centering probability in Eq.~(\ref{eq:P_miscenter}).  
We fix $\alpha_{0}=3$ in Eq~(\ref{eq:general_mass_prof_MAH_relation})
and the lens redshift to be $z_{l}=0.3$.
The marginalized error for each parameter $p_{\alpha}$ 
is given by $\sqrt{F^{-1}_{\alpha \alpha}}$, 
while the un-marginalized error is $1/\sqrt{F_{\alpha \alpha}}$.
We also assume the standard $\Lambda$CDM parameters in
\citet{2016A&A...594A..13P} except for $\sigma_{8}$.
The fiducial parameters of $\bd{p}$ are set to be
$\bd{p}_{\rm fid}
=\{14.5, 1.2, 0.8, 0.42, 0.83\}$.
For the covariance matrix, we adopt the model developed in \citet{2015MNRAS.449.4264G}.
The model includes the statistical uncertainty due to the intrinsic shape, 
the uncorrelated large-scale structure, and mass-dependent term due to intrinsic variations in the mass density profiles. 
When computing the covariance, we perform the linear binning in angular scale $\theta$
from 1 arcmin to 30 arcmin with bin width of 2 arcmin.
To compute the Fisher matrix, 
we only consider the second term in Eq~(\ref{eq:Fisher})
because the mass dependence of the covariance is expected 
to be weak \citep{2015MNRAS.449.4264G}.

\begin{table*}
\centering
\caption{
	A $1\sigma$-level statistical uncertainty
    of our model parameters in hypothetical
    HSC and LSST-like galaxy imaging surveys.
    The number with a round bracket represents
    un-marginalized uncertainty,
    while the marginalized one is for the number
    without bracket.
    The survey parameters (e.g. the sky coverage)
    are the same as 
    in Figure~\ref{fig:fisher_circles}.
    \label{tab:onesigma}
	}
\begin{tabular}{@{}lccccccl}
\hline
\hline
Survey & $\log(M_{\rm 200m}[h^{-1}M_{\odot}])$ & $\alpha_1$ & $f_{\rm cen}$ & $\sigma_{\rm off}\, [h^{-1}{\rm Mpc}]$ & $\sigma_8$ \\ \hline 
HSC-like  & 
0.207 (0.0423) & 0.503 (0.0636) & 0.482 (0.0358) & 0.611 (0.109) & 0.666 (0.0462)  \\ \hline 
LSST-like & 
0.0511 (0.0110) & 0.126 (0.0161) & 0.116 (0.00914) & 0.149 (0.0271) & 0.163 (0.0117)  \\ \hline 
\end{tabular}
\end{table*}

Figure~\ref{fig:fisher_circles} summarizes
the expected constraints of our model parameters
with stacked lensing analysis in ongoing and future imaging surveys.
In this figure, we consider two representative examples, including an ongoing imaging Subaru Hyper-Suprime Cam (HSC) and an upcoming Large Synoptic Survey Telescope (LSST)-like survey. For the HSC-like survey, we assume the sky coverage 
of $1400\, {\rm deg}^{2}$ with the source number density of $30\, {\rm arcmin}^{-2}$.
On the other hand, the LSST-like survey is assumed to cover the sky of $20000\, {\rm deg}^{^2}$
with the source number density of $50\, {\rm arcmin}^{-2}$.
In either case, we set the source redshift to be 1
and the scatter of intrinsic shape in sources of 0.3.

Figure~\ref{fig:fisher_circles} shows the underlying degeneracies among 
our model parameters. 
For instance, more massive haloes tend to show larger signals around the virial radius,
while increasing $\sigma_8$ induces a similar effect on the signal through 
the larger clustering of neighbouring haloes.
Interestingly, the $\alpha_{1}$ parameter 
is found to be almost independent of other parameters, but there exists a degeneracy between $\alpha_1$ and the off-centering effect.
According to our model, 
the $\alpha_1$ parameter can modify the lensing signals at the scale smaller than the virial radius by changing halo concentration, while the off-centering effect can smear the signal at the similar scales. 

After marginalizing over the relevant parameters, we find that the HSC-like survey cannot provide a meaningful constraints on the mass-profile-MAH relation.
Nevertheless, we expect that future wide-area and deep imaging surveys like LSST can measure the mass-profile-MAH relation with good accuracy.
Table~\ref{tab:onesigma} summarizes the expected constraints on the model parameters after marginalization, demonstrating that it is possible to distinguish two models with $\alpha_1=1$ and $\alpha_1=1.2$ at the $1.6\sigma$ level. 
It is worth noting that detectability of baryonic mass contraction in clusters will depend 
on the fiducial value of $\alpha_1$. 
When the fiducial value of $\alpha_1$ is set to be 1.4,
we confirm that LSST-like survey can distinguish the models with $\alpha_1=1$ and $\alpha_1=1.4$ at the $2.8\sigma$ level. 
The constraints as in Table~\ref{tab:onesigma} are found to be insensitive to the choice of 
fiducial $\alpha_1$.
Note that the model with $\alpha_{1}=1$ corresponds to the prediction of dark-matter-only simulations. Hence, the stacked lensing analyses based on future surveys should enable us to measure the baryonic effect on the concentration of galaxy clusters.

\subsection{Bias in parameter estimation}
\label{subsec:bias_param}

Finally, we examine how the parameters in our model will be biased by assuming different baryonic physics, which is encapsulated in the $\alpha_{1}$ parameter. The bias $\delta p_{\alpha}$ for the parameter $p_{\alpha}$ can be expressed as \citep[e.g.,][]{Huterer2006}
\begin{eqnarray}
\delta p_{\alpha} = 
\sum_{\beta}F^{-1}_{\alpha \beta}
\sum_{i,j}
C^{-1}_{ij} \Delta \langle \gamma_{+} \rangle(\theta_{i})
\frac{\partial \langle \gamma_{+} \rangle(\theta_{j})}{\partial p_{\beta}},
\end{eqnarray}
where $i$ and $j$ are the indices for binned lensing profiles,
and $\Delta \langle \gamma_{+} \rangle$
describes the difference between the assumed and true models.
In our case, we are interested in how the parameters will be biased when we vary $\alpha_{1}$ from $1.0$, corresponding to no baryonic effects, to $1.2$ the best-fit value from our {\em Omega500} simulations with AGN feedback. Hence we define
\begin{eqnarray}
\Delta \langle \gamma_{+} \rangle (\theta)
\equiv \langle \gamma_{+} \rangle (\theta| \alpha_{1}=1.2)-\langle \gamma_{+} \rangle (\theta| \alpha_{1}=1),
\end{eqnarray}
where we use the same parameters as in Section~\ref{subsec:detect} except for $\alpha_{1}$.

Assuming the source redshift to be 1, the scatter of the intrinsic shape in source galaxies is 0.3, and the source number density is $50\, {\rm arcmin}^{-2}$, we find that the bias in $\log M_{\rm 200m}$ is $1.3\times10^{-3}$, and the bias in $\sigma_{8}$ is $3.5\times10^{-3}$.
These biases can be compared to the expected statistical uncertainty given 
by the Fisher matrix.
Assuming the covariance scales with the number of foreground objects,
we find the marginalized statistical uncertainty in $\log M_{\rm 200m}$
is given by $0.85\times \sqrt{10/N_{\rm stack}}$, 
where $N_{\rm stack}$ is the number of foreground objects in the stacked analysis.
Thus, $\log M_{\rm 200m}$ will have a potential bias
on the order of $5\%$ of the statistical uncertainty with $N_{\rm stack}=10,000$.
On the other hand, the marginalized statistical uncertainty in $\sigma_8$ 
is approximated as $2.5\times \sqrt{10/N_{\rm stack}}$,
meaning that the bias in $\sigma_{8}$ is only at the level of $4\%$ 
of the statistical uncertainty for $N_{\rm stack}=10,000$.
Note that the mass of $M_{200m}=10^{14.5}\, h^{-1}M_{\odot}$
approximately corresponds to the typical mass of optically selected galaxy clusters 
and $\sim10,000$ clusters are already available \citep[e.g.,][]{2014ApJ...785..104R}.
Hence, we conclude that the systematic uncertainty due to baryonic effects on stacked lensing analysis will be {\em smaller} than the statistical uncertainty of future galaxy-imaging surveys, such as WFIRST, LSST, and EUCLID. 


\section{CONCLUSIONS}
\label{sec:conclusions}


In this paper, we present a simple model to constrain the effects of baryonic physics in the mass distributions of galaxy clusters, currently one of the major systematic uncertainties in using clusters to constrain cosmology.  
Our model is based on three key assumptions: (1) baryonic processes change the mass distribution via the mass concentration of the halo; (2) there exists a tight `mass profile-MAH' relation, between the mass concentration and the mass accretion history (MAH) for a statistical sample of clusters; and (3) the MAH is determined by mass, redshift and cosmological models, insensitive to baryonic processes. We validated and calibrated our model  using a suite of the {\em Omega500} cosmological hydrodynamical simulations with varying input baryonic physics: gas cooling, star formation and AGN feedback. Our findings are summarized as follows:

\begin{enumerate}

\vspace{2mm}
\item We showed that there exists a fairly tight statistical relation between the enclosed mass density measured within the NFW scaled radius $r_s$, and the critical density of the universe at the time when the main-progenitor mass equals to the mass enclosed within $r_s$, consistent with previous results based on dark-matter-only cosmological simulations. Furthermore, we find that this `mass profile-MAH' relation remains valid in the presence of baryonic physics (cooling, star formation and feedback). Baryonic physics is responsible for changing the slope of the relation (see Figure~\ref{fig:rho_MAH}).

\vspace{2mm}
\item
The changes induced by baryonic physics in the mass profile-MAH relation originate from the increase of the mass concentration, but not from the MAH. In fact, we showed that the MAH is insensitive to baryonic effects (see the coloured points in the right top panel in Figure ~\ref{fig:mass_MAH_comp}). The change in the mass concentration predicted by our model provides a good match to the recent observational constraints from X-ray selected clusters \citep[see Figure~\ref{fig:cprof_Mass_comp_obs_model}]{2016ApJ...821..116U}.

\vspace{2mm}
\item
We showed that cluster mass concentration depends on baryonic physics much more sensitively than on the background cosmology. Precise measurements of the mass concentration with upcoming LSST-like optical survey will allow us to constrain baryonic physics in galaxy clusters (see Figure~\ref{fig:fisher_circles}).

\vspace{2mm}
\item
We applied our model to weak lensing mass measurements and showed that the baryonic effects on mass concentration can introduce the bias in weak-lensing mass estimates at a level of $\sim5\%$ for the realistic marginalized constraint with the stacked lensing measurements for $\sim 10,000$ clusters.
Our model provides a theoretical framework for {\em simultaneously} constraining baryonic physics and cluster masses with upcoming optical surveys. 
\end{enumerate}

This paper presents a theoretical framework for constraining baryonic physics in weak lensing mass, using the mass profile-MAH relation. In particular, the parameters of the relation (i.e., slope and amplitude of the relation) can be used as priors in weak lensing cluster masses. 
\change{The upper limit of the slope can be set to be 1.5, which is expected from our simulation with radiative cooling and star formation (without AGN feedback). 
The lower limit will be 1.0, corresponding to the case without any baryonic effects.}
These two values bracket the uncertainty due to baryonic physics in weak lensing mass estimates. 

Future work should focus on developing and analyzing a larger sample of simulated clusters to verify the model in higher-redshift and lower-mass regimes. In addition, how to relate MAH and other structural properties in mass distribution, e.g., asphericity, still remain uncertain. Addressing these issues is the critical step in understanding the remaining astrophysical uncertainties and hence being able to make accurate and robust interpretations of upcoming cluster surveys.


\section*{acknowledgments} 
We thank Keiichi Umetsu for sharing observational constraints of the mass-concentration relation and useful discussions. We also thank Hironao Miyatake and anonymous referee for helpful comments on the manuscript.
DN and EL acknowledge the support from the NSF grant AST-1412768 and the facilities and staff of the Yale Center for Research Computing.
Numerical computations presented in this paper were in part carried out on the general-purpose PC farm at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan.


\bibliographystyle{mnras}
