
\section{Introduction}
In recent years, the Sunyaev-Zeldovich (SZ) effect \citep{Zeldovich1969,Sunyaev1970} has become a powerful tool for studying the warm-hot gas in the universe \citep[][for a recent review]{Planck2014XX, Pratt2019}. Advances in high-resolution SZ observations promise to improve our understanding of thermodynamic and kinetic properties of the warm hot universe \citep[see][for a recent review]{Mroczkowski2019}. The thermal and kinematic SZ effects have been widely used to measure thermal pressure \citep[e.g.,][]{Mroczkowski2009,Planck2013V} and the peculiar velocity of clusters \citep[e.g.,][]{Sayers2019}, with applications for both astrophysical and cosmological studies.

A promising next step lies in the study of the relativistic SZ (rSZ) corrections \citep{Sazonov1998, Nozawa2006, Chluba2012} which depend on the SZ temperature. Despite many attempts in single cluster observations \citep[e.g.,][]{Hansen2002, Hansen2004, Zemcov2012, Prokhorov2012, Chluba2013, Butler2021} and stacking analyses \citep[e.g.,][]{Hurier2016,Erler2018}, the detection significance of the rSZ effect remains low, and to date direct constraints on the SZ temperatures have been limited.
However, the continued advances in both angular resolution and sensitivities, as well as upcoming high-resolution, multi-frequency SZ observations, such as CCAT-Prime \citep{Stacey2018}, NICA2 \citep{NIKA2}, TolTEC \citep{TolTEC}, and The Simons Observatory \citep{TheSimonsObservatoryCollaboration2018}, promise to measure ICM temperature using the rSZ effect for individual galaxy clusters and stacked samples of groups, especially in high-redshift ($z\gtrsim 1$) objects \citep{Morandi2013} without any aid of X-ray observations.
Moreover, even when the SZ temperatures cannot be measured with sufficient statistical certainty, the omission of rSZ effects leads to an overestimation of the Compton-$y$ parameter, especially for the most massive clusters. This can impact ongoing and future cosmological analysis, showing that it is imperative to identify ways of including rSZ in the signal modeling \citep{Remazeilles2018, Remazeilles2019}.

Modern hydrodynamical cosmological simulations have become an indispensable tool for understanding the ICM structures and evolution and their impact on SZ observables \citep[e.g.,][]{Nagai2006,Nagai2007b,Battaglia2012,Kay2012,Pike2014,Yu2015,Planelles2017,LeBrun2017,Henden2018,Henden2019,Pop2022a,Pop2022b}. However, these simulations are also known to exhibit significant variations among different hydrodynamic codes \citep{Rasia2014, Sembolini2016}, which gives rise to differences in a variety of cluster observables, such as the hydrostatic mass derived from X-ray mocks \citep{Rasia2006,Nagai2007a}.  
It is hence important to develop complementary approaches for measuring the ICM temperature based on SZ effect observations. The rSZ effect with upcoming SZ observations promises to provide insights into the thermodynamic structure and evolution \citep{Lee2020}.

Due to the inherent variation of temperature within clusters and groups, both being neither isothermal nor spherically distributed \citep[e.g.,][]{Nagai2007b, Vikhlinin2009cosmoparams}, the temperature used in cluster analysis depend son the context in which they are applied. This has been discussed in detail for SZ measurements alongside X-ray observations \citep[e.g.,][]{Pointecouteau1998,Hansen2004,Kay2008,Lee2020}.
While simulations encompass many variations, it is exciting to explore what can be learnt from common behaviour between simulations. In particular, we find the SZ temperatures demonstrate a consistency between simulations, allowing for predictions of the rSZ corrections that can be detected observationally, in addition to determining what observations may tell us about the underlying physics.

In this work, we study a sample of clusters and groups extracted from 4 different hydrodynamical cosmological simulations in a follow up to earlier work aiming at establishing detailed rSZ temperature-mass relations \citep{Lee2020}. In particular, we use the hydrodynamical simulations $\textsc{Bahamas}$ \citep{McCarthy2017, McCarthy2018};  $\textsc{Macsis}$ \citep{Barnes2017}; {\sc The Three Hundred Project} \citep{Cui2018}; {\sc Magneticum Pathfinder} \citep{Hirschmann2014, Bocquet2016}; and {\sc IllustrisTNG} \citep{Springel2018,Pillepich2018FirstResults,Nelson2018,Marinacci2018, Naiman2018}. These simulations provide large samples of groups and clusters 
over 5 redshifts between $z=0$ to $z=1.5$. In each simulation, we study three different temperature measures: the spectroscopic-like temperature \citep{Mazzotta2004}; the mass-weighted temperature, which is closely related to the compton-$y$ parameter; and finally the y-weighted temperature, which is a close approximation for the averaged rSZ temperature, needed to account for the rSZ corrections \citep{Hansen2004, Remazeilles2019}.

Our work is an extension to a growing body of literature that uses ensembles of the latest hydrodynamical simulations to estimate theoretical uncertainties in cluster scaling relations. These uncertainties arise from our lack of knowledge of the true astrophysical mechanisms at play in our Universe. Either some, or all, of the simulations we use in this work have been used previously to marginalize over astrophysics models and estimate scaling relation uncertainties for different cluster properties, such as the thermal gas pressure \citep{Lim2021}, the central and satellite galaxy stellar properties \citep{Anbajagane2020}, and the dark matter and satellite galaxy velocity dispersions \citep{Anbajagane2022}. Here we extend the discussion to the modeling of the rSZ effect.

This paper is structured as follows. In Section~\ref{sec:Temperature_Weightings}, we discuss the mathematical background behind each temperature measure, alongside some of the formalism we use. We discuss the simulation and our halo samples in Section~\ref{sec:Simulations}. In Section~\ref{sec:TS}, we examine the behaviour of the cluster-averaged temperatures in each simulation and how they vary with mass, redshift, radius, Compton-$y$ parameter and other temperature measures. Section~\ref{sec:cross_sim_fits} contains an exploration of the cross-simulation averaged properties and fits. We discuss applications of our results in Section~\ref{sec:applications}. We summarize our findings in Section~\ref{sec:Conclusion}.

\vspace{-4mm}
\section{Theory}
\label{sec:Temperature_Weightings}
In this paper, we follow much of the same formalism introduced in \cite{Lee2020}; however, for completeness, we briefly review it here. We start by examining each temperature measure. We then briefly discuss observed temperature scaling relations, and finally mention a number of specifics relevant to calculating temperatures consistently within simulations.

In general, when we consider groups and clusters within simulations, we are in fact calculating the signal from a sphere co-located with the dark matter halo, of a certain radius, $R_\Delta$. These radii are defined as the radii containing a certain averaged density, and thus mass $M_\Delta$, as will be explored further in Sect.~\ref{sec:MvR_theory}.

These haloes evidently do not have constant densities nor temperatures throughout the observed sphere, so any observed temperature will instead be an averaged quantity. These averages are, however, always defined by the weighting procedure, which depends on the observable at hand. We can therefore define
\begin{equation}
\langle T\rangle \equiv \frac{\int wT\,{\rm d}V}{\int w\,{\rm d}V},
\end{equation}
where $w$ represents the weighting. As discussed later in this section, for the mass-weighted, $y$-weighted and spectroscopic-like temperatures one respectively has $w=n$, $w=nT$ and $w=n^2T^{-\alpha}$ with $\alpha\simeq0.75$ and $n$ and $T$ are the electron density and temperature, respectively.

\vspace{-4mm}
\subsection{SZ Temperatures}
\topic{Compton-$y$ parameter}: The classical thermal SZ (tSZ) signal has an amplitude proportional to the Compton-$y$ parameter. This Compton-$y$ parameter is proportional to the integrated electron pressure, $P_{\rm e}$:
\begin{equation}
\label{eq:yparameter}
    y \equiv \int \frac{{k_{\rm B}} {T_{\rm e}}}{{m_{\rm e}} c^2}\,{\rm d}\tau = \frac{{\sigma_{\rm T}}}{{m_{\rm e}} c^2}\int P_{\rm e}\,{\rm d}l = 
    \frac{{\sigma_{\rm T}} {k_{\rm B}}}{{m_{\rm e}} c^2}\int nT\, {\rm d}l.
\end{equation}
Here, $\tau$ is the Thomson scattering optical depth; ${\rm d}l$ parametrizes the integral along the line of sight. All the other constants have their usual meaning.

The non-relativistic tSZ effect is then written as \citep{Zeldovich1969}
\begin{equation}
\label{eq:non-rel-tSZ}
\Delta I_\nu \equiv I_0\, y\, g(x) = I_0\, y \frac{x^4 \expf{x}}{(\expf{x}-1)^2}\left(x\frac{\expf{x}+1}{\expf{x}-1}-4\right),
\end{equation}
where $g(x)$, the tSZ spectral function, is defined implicitly. We also use the dimensionless frequency $x=h\nu/{k_{\rm B}} T_{\rm CMB}$, with the present-day CMB temperature $T_{\rm CMB}=2.7255\,{\rm K}$. The intensity normalisation constant is $I_0 = 270.33 (T_{\rm CMB}/2.7255{\rm K})^3$~MJy/sr.

\vspace{2mm}
\topic{Mass-weighted temperature}:
Eq.~\eqref{eq:yparameter} and \ref{eq:non-rel-tSZ} directly motivate the use of a mass-weighted or $\tau$-weighted temperature:
\begin{equation}
    T_{\rm m} \equiv \frac{\int T\,{\rm d}m}{\int {\rm d}m} = \frac{\int nT\,{\rm d}V}{\int n\,{\rm d}V},
\end{equation}
with $m$ the mass of the electron gas. Hence the volume integrated Compton-$y$ parameter, $Y$, is
\begin{equation}
    Y = \frac{{\sigma_{\rm T}} {k_{\rm B}}}{{m_{\rm e}} c^2}\int nT\, {\rm d}V \propto M_{\rm gas}T_{\rm m}
\end{equation}
where $M_{\rm gas}$ is the total gas mass within the halo. From observations, this temperature measure can be estimated by combing tSZ and X-ray measurements, where the latter is used to obtain a mass/$\tau$ estimate. However, since the X-ray temperature does not have the same weighting (see below) and because non-thermal pressure contributions can affect the inference, this leads to a mass bias \citep[e.g.,][]{Arnaud2005, Nagai2007a, Battaglia2012SZphysics, Nelson2012, Nelson2014b, Shi2016}.

\vspace{2mm}
\topic{$y$-weighted temperature}: The relativistic corrections to the SZ effect lead to a temperature-dependent modification to the spectral function which can be accurately modelled using {\tt SZpack} \citep{Chluba2012,Chluba2013}. This implies that we should write for the relativistically-corrected SZ signal, $\Delta I_\nu = I_0 y f(x,{T_{\rm e}})$, in an isothermal cluster. However, since the temperature varies within each cluster, we can write an expansion about an arbitrary temperature pivot $\bar{T}$ \citep{Chluba2013,Remazeilles2018}. To second order in $\Delta T={T_{\rm e}}-\bar{T}$, we have
\begin{equation}
    \frac{\Delta I}{I_0}\simeq yf(x,\bar{T})+y^{(1)}(\bar{T})f^{(1)}(x,\bar{T})+\frac{1}{2}y^{(2)}(\bar{T})f^{(2)}(x,\bar{T}),
\end{equation}
where $f^{(k)}=\partial_T^{k}f(x,{T_{\rm e}})$ and $y^{(k)}=\langle(T-\bar{T})^k y\rangle$. Here, $\langle X\rangle$ is indicating the cluster averaged value. This motivates a relativistic temperature, which removes the first-order correction (i.e., $y^{(1)}(\bar{T})=0$) and we call the $y$-weighted temperature
\begin{equation}
    T_{\rm y} \equiv \frac{\int y T\,{\rm d}V}{\int y\,{\rm d}V} =\frac{\int nT^2\,{\rm d}V}{\int nT\,{\rm d}V}. 
\end{equation}
We note that there will still be contributions from higher order temperature terms, however in \citet{Lee2020}, these are determined to lead to $\simeq0.1\%$ corrections to the SZ signal peak, and as such we neglect further study of them here. We find that this $y$-weighted temperature is systematically higher than the mass-weighted and X-ray temperatures, indicating that, especially for the largest clusters in the Universe, the relativistic corrections will be relevant and may bias cosmological inferences if not modelled.

\vspace{-4mm}
\subsection{X-ray Temperatures}
X-ray emission, from hot clusters (${k_{\rm B}} {T_{\rm e}}\gtrsim 3$~keV), is dominated by bremstrahllung radiation within the ICM, while cooler clusters have increasing contributions from emission lines.
However, in this hot regime, temperature weightings have been determined by fitting the X-ray spectra with a thermal emission model \citep{Mazzotta2004, Vikhlinin2006}. This has led to the \topic{spectroscopic-like temperature},
\begin{equation}
    T_{\rm sl} \equiv \frac{\int n^2T^{1-\alpha}\,{\rm d}V}{\int n^2T^{-\alpha}\,{\rm d}V}
\end{equation}
with $\alpha=0.75$ for high temperatures (${k_{\rm B}} {T_{\rm e}}>3.5$~keV), which matches well with observations from both {\it Chandra} and {\it XMM-Newton}. X-ray temperatures have also been calibrated within simulations, to determine the differences between different X-ray temperatures and to confirm the weighting for $T_{\rm sl}$ \citep[e.g.,][]{Rasia2014}.

\vspace{-4mm}
\subsection{Halo definition and redshift dependence}
\label{sec:MvR_theory}
\topic{Radius definitions}: As previously mentioned, to define the extent of groups and clusters within our simulations, we consider spheres of a given radius, co-located with the halo, i.e., centred on the minimum of potential. Then a halo of radius $R_\Delta$, is defined so that it contains the mass $M_\Delta$.

In this work, we use five different radii, in particular, $R_{\rm 500c}$, $R_{\rm 200c}$, $R_{\rm 500m}$, $R_{\rm 200m}$ and $R_{\rm vir}$. For $R_{\rm \Delta c}$, these are defined in terms of an overdensity of $\Delta$ times the critical density of the universe, $\rho_{\rm crit}$. i.e.,
\begin{equation}
    M_{\rm \Delta c} = \frac{4\pi}{3}R_{\rm \Delta c}^3 \rho_{\rm crit}(z) \Delta,
\end{equation}
while the $R_{\rm \Delta m}$ are defined in terms of overdensities of $ \rho_{\rm mean}\Delta$, the mean density of the universe. The critical density and mean density are defined as,
\begin{equation} \begin{split}
    \rho_{\rm crit} &\equiv \frac{3H_0^2}{8\pi G}\left[\Omega_m(1+z)^3+\Omega_\Lambda\right] = \frac{3H_0^2}{8\pi G}E^2(z)\\
    \rho_{\rm mean} &\equiv \frac{3H_0^2}{8\pi G}\left[\Omega_m(1+z)^3\right]
\end{split} \end{equation}
for a flat $\Lambda$CDM universe. We have also here implicitly defined $E(z)$, which describes the redshift evolution of the Hubble parameter.
The virial radius $R_{\rm vir}$ is defined by $\rho_{\rm crit}$ and $\Delta_{\rm vir}$ calculated using the approximation in \cite{Bryan1998}, that is, $\Delta_{\rm vir}\approx102$ at $z=0$.

In this work, we predominantly define our clusters as being within the radius $R_{\rm 200c}$, a common proxy for the observational radii used in many SZ calculations.
We note that here we are always using the {\it true} (simulation) mass, rather than any proxy for the observed mass (i.e., the hydrostatic mass), which could introduce observational biases \citep[e.g.,][]{Nagai2007a,Lau2009,Lau2013,Nelson2014a,Shi2016,Biffi2016,Barnes2017}.

\topic{Redshift dependence}: As the critical density has a redshift dependence, we can accordingly expect an evolution of the cluster properties. A simple geometric consideration and assumption of isothermality within the virialized sphere would lead to a temperature dependence akin to
\begin{equation}
    T_\Delta \propto M_\Delta/R_\Delta.
\end{equation}
Then, by using the $\rho_{\rm crit}$ redshift dependence, we would expect a self-similar evolution of
\begin{equation} \begin{split}
    M_\Delta &\propto E^2(z) R_\Delta^3,\\
    T_\Delta &\propto E^{2/3}(z) M_\Delta^{2/3}.
\end{split} \end{equation}
As such we can consider any deviation from this evolution as being inherent evolution of the temperatures, due to non-gravitational processes in the clusters themselves.

\vspace{-4mm}
\section{Simulations}
\label{sec:Simulations}
\begin{table}
\caption{Cosmological parameters used in each of the simulations.}
\begin{tabular}{l |@{}c |@{}c |@{}c |@{}c |@{}c |@{}c |@{}c }
    \hline
    Simulation & $\Omega_\Lambda$ & $\Omega_\mathrm{m}$ & $\Omega_\mathrm{b}$ & $\sigma_8$ & $n_\mathrm{s}$ &$h^\dagger$\\
    \hline
    $\textsc{Bahamas}$ & 0.6825 & 0.3175 & 0.0490 & 0.8340 & 0.9624 & 0.6711 \\
    $\textsc{Macsis}$ & 0.6930 & 0.3070 & 0.0482 & 0.8288 & 0.9611 & 0.6777 \\
    $\textsc{The300}$ & 0.6930 & 0.3070 & 0.0480 & 0.8230 & 0.9600 & 0.6780 \\
    $\textsc{Magneticum}$ & 0.7280 & 0.2720 & 0.0457 & 0.8090 & 0.9630 & 0.7040 \\
    $\textsc{TNG}$ & 0.6911 & 0.3089 & 0.0486 & 0.8159 & 0.9667 & 0.6774 \\
    \hline \noalign{\smallskip}
    \multicolumn{7}{l}{$^\dagger$ where $h \equiv H_0/(100 \,\mathrm{km\,s}^{-1}\,\mathrm{Mpc}^{-1})$ }
\end{tabular}
\label{tab:SimulationCosmology}
\end{table}
\begin{table*}
	\centering
	\begin{tabular}{l|c|c|c|c|c|l}
		Simulation & $L$ [Mpc] & $N$ & $\epsilon_{\rm \,DM}^{z=0} \,\rm [kpc]$ & $m_{\rm DM} \, [\msol]$ & $m_{\rm gas} \, [\msol]$ & Calibration \\
        \hline
            $\textsc{Bahamas}$ & $596$ & $2 \times 1024^3$ & $5.96$ & $6.7 \times 10^9$ & $1.2 \times 10^9$ & GSMF, CL $f_{\rm gas}$ \\
            $\textsc{Macsis}$ & ---   & ---    & $5.90$ & $6.5 \times 10^9$ & $1.2 \times 10^9$ & GSMF, CL $f_{\rm gas}$ \\
            $\textsc{The300}$ & ---   & ---    & $9.6$  & $1.9 \times 10^9$ & $3.5 \times 10^8$ & See \citet{Cui2018} \\
            $\textsc{Magneticum}$ & $500$ & $2 \times 1584^3$& $5.33$ & $9.8 \times 10^8$ & $2.0 \times 10^8$ & SMBH, Metals, CL $f_{\rm gas}$ \\
            $\textsc{TNG}$ & $303$ & $2 \times 2500^3$ & $1.48$ & $5.9 \times 10^7$ & $1.1 \times 10^7$ & See \citet{Pillepich2018SimulatingGalaxyFormation} \\
            \hline
    \end{tabular}
    \caption{Simulation characteristics for each sample. Table adapted from \citet[Table 1]{Anbajagane2022}. From left to right, we show: (i) simulation, (ii) comoving box size, (iii) the total particle count, (iv) force softening scale, (v) mass of DM particles, (vi) mass of gas particles and (vii) empirical sources used for tuning sub-grid parameters of each simulation, which consist of the Galaxy Stellar Mass Function (GSMF), Supermassive Black Hole scaling (SMBH), Metallicity scaling (Metals), and cluster hot gas mass fraction $<R_{500c}$ (CL $f_{\rm gas}$). We do not quote the box sizes or particle numbers for $\textsc{Macsis}$ and $\textsc{The300}$ as they are zoom-in simulations.}
	\label{tab:SimulationParams}
\end{table*}

We use four samples of haloes from the following simulations: (i) a superset of {\sc Bahamas} and {\sc Macsis}, (ii) {\sc The Three Hundred Project}, (iii) {\sc Magneticum Pathfinder}, and (iv) {\sc IllustrisTNG}. A summary of the cosmological properties of each simulation can be found in Table~\ref{tab:SimulationCosmology} and the simulation parameters in Table~\ref{tab:SimulationParams}.

We take haloes within each simulation with $M_{\rm 500c} >10^{13}$~$\msol$ and consider the redshifts $z \simeq \{0.0, 0.25, 0.5, 1.0, 1.5\}$. The population counts at each redshift can be found in Table~\ref{tab:SimulationPopulations}.
In each simulation, we calculate our properties using a temperature cut -- that is, only considering particles with temperatures, $T \geq 10^{5.2}$~K. Moreover, when calculating the {\it spectroscopic-like} temperature, we use a core excision procedure as has been shown to make X-rays a better proxy for mass \citep[e.g.,][]{Pratt2009}, and exclude all particles with $r<0.15\,R_{\rm 500c}$. We also note in Appendix~\ref{app:relaxation}, that using a relaxed subsample does not seem to meaningfully change our results at $z=0$, so this is not further examined here.

\begin{table*}
    \begin{tabular}{l|cc|cc|cc|cc|cc}
        Nominal redshift & \multicolumn{2}{c|}{$z=0.00$} & \multicolumn{2}{c|}{$z=0.25$} & \multicolumn{2}{c|}{$z=0.50$} & \multicolumn{2}{c|}{$z=1.00$} & \multicolumn{2}{c}{$z=1.50$} \\
         & $z_{\rm sim}$ & $N_{\rm haloes}$ & $z_{\rm sim}$ & $N_{\rm haloes}$ & $z_{\rm sim}$ & $N_{\rm haloes}$ & $z_{\rm sim}$ & $N_{\rm haloes}$ & $z_{\rm sim}$ & $N_{\rm haloes}$ \\ \hline
        $\textsc{Bahamas+Macsis}$ & 0.00 & 21361 & 0.25/0.24 & 19610 & 0.50/0.46 & 17016 & 1.00 & 10928 & 1.50/1.56 & 5750 \\
        $\textsc{The300}$ & 0.00 & 8465 & 0.25 & 8433 & 0.49 & 8300 & 0.99 & 7349 & 1.48 & 5627 \\
        $\textsc{Magneticum}$ & 0.066 & 10429 & 0.25 & 9812 & 0.52 & 8265 & 0.96 & 5659 & 1.47 & 3080 \\
        $\textsc{TNG}$ & 0.00 & 2548 & 0.24 & 2333 & 0.50 & 2015 & 1.00 & 1290 & 1.50 & 700
    \end{tabular}
    \caption{The number of haloes with $M_{\rm 500c}>10^{13}$~$\msol$ within each simulation, at each considered redshift. The nominal redshift is the how we will refer to each sample, while $z_{\rm sim}$ indicates the `true' simulation redshift used. For the $\textsc{Bahamas+Macsis}$ sample, where the redshifts are recorded as `x/y', the first number refers to $z_{\rm sim}$ in the $\textsc{Bahamas}$ simulation, while the second that in the $\textsc{Macsis}$ simulation.}
    \label{tab:SimulationPopulations}
\end{table*}

\vspace{-4mm}
\subsection{Bahamas and Macsis}
We use a supersample of $\textsc{Bahamas}$ \citep{McCarthy2017, McCarthy2018} and its zoom-in counterpart $\textsc{Macsis}$ \citep{Barnes2017}, which we will here refer to as $\textsc{Bahamas+Macsis}$, or in figures simply $\textsc{Bahamas}$+.
$\textsc{Bahamas}$ is a calibrated version of the model used in the cosmo-OWLS simulations \citep{LeBrun2014}. Both $\textsc{Bahamas}$ and $\textsc{Macsis}$ were run using a version of the smoothed particle hydrodynamics (SPH) code {\tt GADGET-3} last publicly discussed in \cite{Springel2005}, but modified as part of the OWLS Project \citep{Schaye2010}. The $\textsc{Macsis}$ simulation was developed to extend the $\textsc{Bahamas}$ sample to higher mass haloes. It is a sample of 390 clusters, generated through a zoomed simulation from a large Dark Matter only (DMO) simulation -- a periodic cube with a side length of 3.2~Gpc. Individual separate volumes including high Friends-of-Friends (FoF) mass (clusters with $M_{\rm FoF}>10^{15}$~$\msol$) were then re-simulated with a full hydrodynamical prescription aligned with the $\textsc{Bahamas}$ simulation. Haloes within the two simulations are then identified by a FoF algorithm, and subhaloes with the {\sc Subfind} algorithm \citep{Springel2001, Dolag2009}.

\vspace{-4mm}
\subsection{The Three Hundred Project}
{\sc The Three Hundred Project}, here $\textsc{The300}$, \citep{Cui2018} comprises of massive haloes formed within 324 spherical regions each of 22 Mpc radius and each centred on the most massive clusters at redshift zero as identified in the {\sc MultiDark Planck 2} N-body, DMO simulation \citep{Klypin2016}, which has a cube of side length $1.476$~Gpc with $3840^3$ dark matter particles and used the {\tt L-GADGET-2} solver. 
The haloes in {\sc MultiDark Planck 2} were identified using the {\sc Rockstar} halo finder \citep{Behroozi2013}.
These 324 spherical regions of radius 22 (comoving) Mpc were then resimulated using a full hydrodynamics prescription \citep{Rasia2015} with the {\tt GADGET-X} SPH solver \citep{Beck2016}. Haloes and subhaloes were identified with Amiga's Halo Finder \citep{Knollmann2009}, which has a binding energy criterion similar to {\sc Subfind}, but uses and adaptive mesh refinement grid to represent the density field/contours. 

While $\textsc{The300}$ is only mass-complete above $M_{\rm 200c} \gtrsim 10^{15}$~$\msol$ at $z=0$, it resolves many haloes below this mass. Since the scaling relations derived from these lower mass clusters are in agreement with our other simulations, we note that the selection effects do not in general appear to have an impact on our temperature scaling results, with one exception. That is, within large radii (i..e, $R_{\rm vir}$ and $R_{\rm 200m}$) the scatter of temperatures is amplified, however, the population mean behaviors remain unaffected.

In Section~\ref{sec:Res_Feedback}, we also use two different simulations from $\textsc{The300}$, which we will refer to as {\sc Music} and {\sc Gizmo}. {\sc Music} uses the solver {\tt GADGET-MUSIC} \citep{Sembolini2013} in place of {\tt GADGET-X}. A summary of the differences between these simulations can be found in \cite{Cui2018}, but in brief, the {\tt GADGET-MUSIC} and {\tt GADGET-X} solvers use similar but subtly different gas treatments and stellar formation and stellar feedback mechanisms. Most notably, however, {\tt GADGET-MUSIC} does not include any Black Hole or AGN feedback, while {\tt GADGET-X} does.

{\sc Gizmo} \citep{Cui2022} on the other hand uses {\tt GIZMO-SIMBA} \citep{Dave2019} in its Meshless Finite Mass solver mode instead of {\tt GADGET-X} which uses SPH. A detailed discussion of the differences can be found in \citet[their Table 2]{Cui2022}. In general, {\sc Gizmo} uses a different feedback model, with, among other things, significantly stronger kinetic feedback calibrated for high mass haloes which causes large differences in the observed gas fractions between the {\tt GIZMO-SIMBA} and {\tt GADGET-X} runs.

\vspace{-4mm}
\subsection{Magneticum Pathfinder}
{\sc Magneticum Pathfinder} \citep{Hirschmann2014} is a suite of magneto-hydrodynamics simulations using a version of the SPH solver {\tt GADGET-3} developed independently to that used in the $\textsc{Bahamas+Macsis}$ versions. We used the {\tt Box2 hr} run, and will henceforth refer to it as the $\textsc{Magneticum}$ sample. Haloes are once again identified using a FoF algorithm, and subhaloes using {\sc Subfind}.
We note that since the $\textsc{Magneticum}$ simulations are based on WMAP7 cosmology \citep{Komatsu2011} they have the most distinct cosmology to the other simulations which all use Planck cosmologies \citep{Planck2014XVI, Planck2016XIII}.
Finally, we observe that while the other simulations show little variation between the core-excised and non-core excised values for $T_{\rm m}$ and $T_{\rm y}$ \citep[as discussed for $\textsc{Bahamas+Macsis}$ in][]{Lee2020}, this is not true in $\textsc{Magneticum}$. As such, we use the core-excised values for all three temperature measures obtained from $\textsc{Magneticum}$. This is explored in more detail in Appendix~\ref{app:core-excision}.

\vspace{-4mm}
\subsection{IllustrisTNG}
The {\sc IllustrisTNG} project, here $\textsc{TNG}$, \citep{Springel2018, Pillepich2018FirstResults, Nelson2018, Marinacci2018, Naiman2018} is a follow up to the {\sc Illustris} project \citep{Vogelsberger2014}. It uses the moving mesh code {\tt AREPO} \citep{Springel2010} and uses a full magneto-hydrodynamics treatment with galaxy formation models as detailed in \citet{Pillepich2018SimulatingGalaxyFormation, Weinberger2017}. We use TNG300-1, the highest resolution run from the suite, but reference the two lower resolution runs in Appendix~\ref{sec:Res_Resolution}. Haloes are identified using a FoF algorithm and subhaloes with {\sc Subfind}, as for the $\textsc{Magneticum}$ sample.

We estimate all TNG properties using the FoF particle set. The FoF linking length of $b = 0.2$ was chosen so that the FoF group on average contains all particles within $R_{\rm 200c}$ of the halo center. Consequently, properties computed within significantly larger apertures (primarily $R_{\rm 200m}$) will miss the contribution from particles in the far outskirts (as these are not included in the FoF, whose particle set is incomplete at such distances) and will thus incur a minor bias. However, our main analysis and results focus on $R_{\rm 200c}$ and $R_{\rm 500c}$ and are thus unaffected by this bias.

\begin{figure}
    \includegraphics[width=\linewidth]{Figures/03_TvM_z0_grid_nohatch.pdf}\\
    \caption{The three temperature measures in each simulation at $z=0$. Here the data is binned into logarithmically spaced mass bins. The shaded regions show the 16 and 84 percentile regions within these bins, while the lines show the medians. The dotted black line indicates the self-similar scaling $T\propto M^{2/3}$, whose amplitude is arbitrarily set for visual clarity. The horizontal line in the lower panel indicates the temperature above which $T_{\rm sl}$ is considered a good proxy for the X-ray temperature.}
    \label{fig:03_TvM_z0}
\end{figure}
\vspace{-4mm}
\section{Temperature Scaling Relations}
\label{sec:TS}
In this section we briefly discuss the temperature measures at redshift $z=0$, before examining how they vary with redshift (Section~\ref{sec:TS_Redshift}) and the choice of the radial cut-off (Section~\ref{sec:TS_Radial}). We also examine how the temperatures scale with respect to the Compton-$y$ parameter, and the $T_{\rm y}$ offset against $T_{\rm m}$ itself, in Sections~\ref{sec:TS_Y} and \ref{sec:TS_DeltaT}. In Section~\ref{sec:Res}, we discuss how the specifics of the simulations themselves have affected our results. 

\begin{figure*}
    \includegraphics[width=0.8\linewidth]{Figures/04_fix_Red_RelTvM_grid.pdf}\\
    \caption{Redshift evolution of the temperature measures within each simulation. As in Fig.~\ref{fig:03_TvM_z0}, the data is sorted into mass bins, with the shaded regions showing the 16 to 84 percentile region and the lines the medians. Here each redshift for each sample is divided by the bin median at $z=0$. We see there is reduced evolution in $T_{\rm y}$ (the right panels) than in $T_{\rm m}$ (those on the left).}
    \label{fig:04_TvM_redshift}
\end{figure*}
In general, we can see a good agreement between the the  4 samples for each temperature measure (Fig.~\ref{fig:03_TvM_z0}). As we will continue to use a similar plotting convention throughout this paper, we briefly explain the process here. We have sorted the data into logarithmically spaced mass bins, and within each bin then plotted the median, and 16\% and 84\% percentile regions for each sample. That is, the solid lines indicate the median of the data distribution, while the shaded regions show the 1$\sigma$ intrinsic scatter within each data set. Where there are fewer than 10 clusters within a mass bin, we have only plotted the median and not calculated the percentiles. Within each sample, this region accounts for an intrinsic variation of around $\pm 5.5\%$ for $T_{\rm m}$ and $T_{\rm y}$ (with this in general being marginally larger in $T_{\rm y}$ than $T_{\rm m}$), and $\pm 9\%$ for $T_{\rm sl}$.

We can see that for the mass-weighted temperature, $T_{\rm m}$, all four simulations show a close agreement at $z=0$, while $T_{\rm y}$ shows slightly more variation, but is, nonetheless, overall consistent. The spectroscopic-like temperature shows the most variation between the four samples, particularly at lower masses/temperatures. However, it is worth noting, that for $T_{\rm sl} \gtrsim 3.5$~keV, where $T_{\rm sl}$ is considered a good proxy for the X-ray temperature, the samples agree well.
In general, however, there is far more intrinsic scatter within the $T_{\rm sl}$ measures (within $R_{\rm 200c}$) than for the other temperature measures. We can see in particular that $\textsc{Magneticum}$ has a large scatter in $T_{\rm sl}$, especially at lower masses. This is driven by the warm gas ($T<10^6$~K) in the low mass haloes, and is discussed in more detail in Appendix~\ref{app:T-cut}. Due to our limited $T_{\rm sl}$ data where it is an appropriate proxy for the X-ray temperature, particularly at higher redshifts, we will not consider $T_{\rm sl}$ in detail for the rest of this section.

In Fig.~\ref{fig:03_TvM_z0}, we have also plotted an arbitrary indicative line to show self-similar scaling, i.e., $T\propto M^{2/3}$. This allows us to see that in simulations all three temperature measures appear to scale at slightly less than $2/3$, with $T_{\rm sl}$ lying closest to this. However, at higher masses and temperatures, $T_{\rm y}$ and $T_{\rm m}$ appear to tend to this self similarity, while at lower masses and temperatures, the scaling relation seems shallower. Conversely, this also indicates there is some mild curvature within the $T_{\rm y}$ and $T_{\rm m}$ scaling relations. The two and three parameter fits for each of these simulations can be found tabulated in Table~\ref{apptab:fits_median_z0_samples} for a more detailed comparison, and will be discussed further in Section~\ref{sec:cross_sim_fits}.
The agreement at high masses may come from the decreased relative effects of feedback in this regime. That is, at lower masses the gas is hotter than expected from solely gravitational heating due to feedback processes, pushing the equilibrium away from self-similarity, while at higher masses the potential well ensures more gas is retained in haloes; this is supported by the results of \citet{Farahi2018}, who use kernel-localized linear regression \citep{Farahi2022} to show that the gas mass in \textsc{Bahamas+} clusters approaches a self-similar scaling at higher halo masses. That the lower masses lie higher than expected from self similarity for $T_{\rm y}$ and $T_{\rm m}$, would indicate that generally the feedback is leading to hotter gas in the halo.

\vspace{-4mm}
\subsection{Redshift evolution}
\label{sec:TS_Redshift}
\begin{table}
\caption{The median value for $E(z)^{-2/3}T(z)/\bar{T}_{z=0}$ as shown in Fig.~\ref{fig:04_TvM_redshift}, for $z=0.5$ and 1.0. These can be understood as a broad mass-independent redshift correction. However, we have tabulated in fact the median correction to the 2-parameter fit for each simulation (described in Section~\ref{sec:Dis_effects} and tabulated themselves in Table~\ref{apptab:fits_median_z0_samples} in Appendix~\ref{app:tables}).} \centering
\begin{tabular}{lccc}
 $z=0.5$ &  $E(z)^{-2/3}T_{\rm y}$ &  $E(z)^{-2/3}T_{\rm m}$ & $E(z)^{-2/3}T_{\rm sl}$ \\ 
 \hline \smallskip
 $\textsc{Bahamas+Macsis}$ & $0.981^{+0.002}_{-0.002}$ & $0.951^{+0.001}_{-0.001}$ & $0.929^{+0.005}_{-0.011}$ \\ \smallskip
 $\textsc{The300}$  & $0.968^{+0.006}_{-0.004}$ & $0.953^{+0.003}_{-0.002}$ & $0.933^{+0.007}_{-0.005}$ \\ \smallskip
 $\textsc{Magneticum}$  & $0.969^{+0.002}_{-0.005}$ & $0.932^{+0.003}_{-0.003}$ & $0.728^{+0.023}_{-0.014}$ \\ \smallskip
 $\textsc{TNG}$  & $0.984^{+0.007}_{-0.005}$ & $0.972^{+0.004}_{-0.005}$ & $0.746^{+0.018}_{-0.038}$ \\
 \hline
 $z=1.0$ &  $E(z)^{-2/3}T_{\rm y}$ &  $E(z)^{-2/3}T_{\rm m}$ & $E(z)^{-2/3}T_{\rm sl}$ \\ 
 \hline \smallskip
 $\textsc{Bahamas+Macsis}$ & $0.955^{+0.003}_{-0.002}$ & $0.896^{+0.002}_{-0.002}$ & $0.815^{+0.006}_{-0.011}$ \\ \smallskip
 $\textsc{The300}$  & $0.938^{+0.004}_{-0.004}$ & $0.903^{+0.002}_{-0.002}$ & $0.856^{+0.004}_{-0.005}$ \\ \smallskip
 $\textsc{Magneticum}$  & $0.966^{+0.004}_{-0.003}$ & $0.889^{+0.003}_{-0.005}$ & $0.590^{+0.020}_{-0.013}$ \\ \smallskip
 $\textsc{TNG}$  & $0.997^{+0.007}_{-0.007}$ & $0.941^{+0.006}_{-0.005}$ & $0.501^{+0.018}_{-0.012}$ \\
 \hline
\end{tabular}
\label{tab:fits_redshift_offset_samples}
\end{table}
As discussed in Sect.~\ref{sec:MvR_theory}, from self similarity we expect the cluster temperatures to scale as $T(z) \propto E(z)^{2/3}$ at fixed mass. In Fig.~\ref{fig:04_TvM_redshift}, we can see that this is not quite the case within our sample.
We examine our temperature measures at 5 different redshifts, $z=0.0$, 0.25, 0.5, 1.0 and 1.5. In general we find that, the temperature measures increase (with increasing redshift) {\it slower} than self similarity would suggest. That is, in Fig.~\ref{fig:04_TvM_redshift}, were self similar evolution to occur, we would expect all 5 redshifts to align -- as is nearly true for $T_{\rm y}$. We note that here we have used the cosmological parameters and true redshift for each individual simulation to calculate this redshift behaviour. Due to the larger differences in cosmology in the $\textsc{Magneticum}$ simulation, this allows for more consistent results between simulations.

We found that $T_{\rm sl}$, while not plotted here, evolves with the least accordance to self similarity, while $T_{\rm y}$ diverges the least from it. That is, graphically, we can see the spread in the 5 redshift means is {\it smaller} for $T_{\rm y}$ than for $T_{\rm m}$. At $z=1$, we can see that $T_{\rm m}$ has a median lowered to around the 1$\sigma$ intrinsic scatter at $z=0$. $T_{\rm y}$ however shows almost no redshift dependence in $\textsc{Magneticum}$ and $\textsc{TNG}$ and only mild evolution in $\textsc{Bahamas+Macsis}$ and $\textsc{The300}$.
Physically, we can motivate this aspect, as at higher redshifts, clusters on average had shorter cooling times, leading to more cool dense gas, which down-weights $T_{\rm sl}$. On the other hand, $T_{\rm y}$ is determined by the gas pressure, which, assuming hydrostatic equilibrium, would be fixed to match the size of the potential well itself, leading to a more self-similar temperature measure.

It is also interesting to note that in general our samples do not show overt mass dependence in the redshift evolution for the two SZ temperatures, i.e., graphically, the mean lines are roughly horizontal for all of our samples. Here, it is worth noting that we are considering the masses at the redshift the cluster temperature is measured. In Table~\ref{tab:fits_redshift_offset_samples}, we quantify these offsets for all three temperature measures in each sample at $z=0.5$ and $1.0$. Here we can more quantitatively see that there is increasing divergence from self-similarity (value of unity), comparing $T_{\rm y}$ to $T_{\rm m}$ to $T_{\rm sl}$ in all cases. We can also see the variety in redshift behaviour for each sample. We also note that the larger errors in the $T_{\rm sl}$ offsets arise from there being more mass dependence to the redshift evolution (particularly for $\textsc{Magneticum}$ and $\textsc{TNG}$) than in $T_{\rm y}$ and $T_{\rm m}$.

\subsection{Radial dependence}
\label{sec:TS_Radial}
\begin{figure}
    \includegraphics[width=\linewidth]{Figures/05_Rad_TvR_bySim_grid.pdf}\\
    \caption{A depiction of the radial dependence of each simulation. Here every cluster's temperature and radius is divided by the {\it same cluster's} radius and $T_{\rm m}$ within $R_{\rm 200c}$, for the five radii considered. These are, from left to right $R_{\rm 500c}$, $R_{\rm 200c}$, $R_{\rm 500m}$, $R_{\rm vir}$ and $R_{\rm 200m}$. Hence the error bars show the intrinsic scatter within each simulation for the clusters temperature profiles, fixed at $R_{\rm 200c}$ and $T_{\rm m,200c}$.}
    \label{fig:05_radial_scatter}
\end{figure}
\begin{figure*}
    \includegraphics[width=0.8\linewidth]{Figures/06_Rad_TvM_grid.pdf}\\
    \caption{The variation in temperature scalings when averaged over three different radii, at $z=0$. This figure is arranged as in Fig.~\ref{fig:03_TvM_z0}, with $T_{\rm m}$ on the left and $T_{\rm y}$ on the right.}
    \label{fig:06_radial_scalings}
\end{figure*}

Next, we study the dependence of the temperature measures on the radius of the sphere that we average over. Fig.~\ref{fig:05_radial_scatter} shows how the averaged temperatures vary over the five radii we consider at $z=0$. Here, on a cluster-by-cluster basis, we divide each radius and temperature of the cluster, by the $R_{\rm 200c}$ and $T_{\rm m,200c}$ values for that cluster, and then we show the bulk averages of these values. These values can be found for each sample for convenient reference in Table~\ref{apptab:R_err_sample} in the Appendix.

We find that the profiles are very similar within the four samples, and are consistent out to $R_{\rm 500m}$. The variations beyond this are likely driven by the particulars of the simulations rather than reflecting any inconsistency. We discuss in Section \ref{sec:Simulations} that when calculating our $\textsc{TNG}$ values, only particles which are linked to the FoF group are included, which may bias the large radius temperatures high. $\textsc{The300}$ shows a consistent profile with the other simulations, but a larger scatter due to the sample selection for the low mass haloes. \citet{Anbajagane2022} found similar scatter amplification in the velocity dispersion of low mass clusters in \textsc{The300}, and explicitly showed this was generated by the selection effect (see their Figure B1).
We can also see that the offset between $T_{\rm y}$ and $T_{\rm m}$ is larger in $\textsc{Bahamas+Macsis}$ and $\textsc{Magneticum}$ than the other two simulations -- this will be explored more in Section~\ref{sec:TS_DeltaT}. This offset also seems to marginally increase at larger radii.

The difference in profile steepness can be appreciated in Fig.~\ref{fig:06_radial_scalings}. Here we show the temperature-mass relations, at $R_{\rm 500c}$, $R_{\rm 200c}$ and $R_{\rm 200m}$, so we can see how the variation in profiles depend on the masses of the clusters. For example, we find that the differences in the profiles between $\textsc{The300}$ and $\textsc{Bahamas+Macsis}$ are largely driven by lower mass haloes, and agree well at all radii for higher masses.

It also becomes evident that, for both $T_{\rm m}$ and $T_{\rm y}$, the four models agree best within the $R_{\rm 200c}$. At lower radii, the shallow profile of $\textsc{TNG}$ corresponds to temperatures below those in the other simulations, and at larger radii to temperatures above. Similarly, the steeper profile of $\textsc{Magneticum}$ leads to its temperatures sinking with respect to the other samples, as the radius increases. $T_{\rm sl}$, while not shown here, has a larger variation at each radius -- and it is harder to determine the agreement between simulations due to the difference in intrinsic scatter in each sample. Within $R_{\rm 500c}$, $T_{\rm sl}$ in $\textsc{Magneticum}$ agrees far better with the other simulations -- however, $T_{\rm sl}$ in all the other simulations agree slightly worse with each other than at $R_{\rm 200c}$. 

\subsection{Temperature--Y scalings}
\label{sec:TS_Y}
\begin{figure}
    \includegraphics[width=\linewidth]{Figures/07_TvY_z0_grid.pdf}\\
    \caption{The variation of each temperature measure with respect to the volume averaged Compton-$y$ within $R_{\rm 200c}$. These are the measures at $z=0$ and the figure is arranged as in Fig.~\ref{fig:03_TvM_z0}, grouped into logarithmically spaced Y bins, rather than mass bins. The dotted black line is an indicative line with scaling $T\propto Y^{2/5}$, i.e., a self-similar scaling.}
    \label{fig:07_TvY_scaling}
\end{figure}
While temperature-mass scaling relations are important and relevant for comparison with X-ray observations, a consideration of the Compton-$y$ parameter in relation to temperature could lead to a way to self-calibrate SZ observables \citep{Lee2020}. In Fig.~\ref{fig:07_TvY_scaling}, we show how the SZ temperatures scale within our samples with $Y_{\rm 200c}$.
We can immediately see that there is slightly poorer agreement within these quantities than when using $M_{\rm 200c}$, however, these do still predominantly agree between the different samples. In particular, as $Y_{\rm 200c}$ increases, the agreement between $T_{\rm m}$ in each sample improves, and as we will see later in Section~\ref{sec:Res}, most of this variation correlates with the variation of $f_{\rm gas}$ in each simulation. 

We note that the $Y_{\rm 200c}-M_{\rm 200c}$ relationship agrees very well between simulations for all masses $\gtrsim 10^{14}$~$\msol$, which corresponds to $Y_{\rm 200c} \gtrsim 10^{-6}$~Mpc$^{-2}$. It is also worth considering that since our haloes are selected with a mass cut-off, for the lowest values of $Y$ ($\lesssim 10^{-7}$~Mpc$^{-2}$) the data we have are not necessarily complete for each value of $Y$, so may be biased slightly high.
Self-similarity would suggest a scaling relation of $T \propto Y^{2/5}$, but we see shallower scaling relations $\simeq 0.33$ in all our samples. In particular, $\textsc{Bahamas+Macsis}$ has a shallower scaling relation ($\simeq 0.31$), for both $T_{\rm y}$ and $T_{\rm m}$ than the other three samples. 

\subsection{$\Delta T$--Temperarature scalings}
\label{sec:TS_DeltaT}
\begin{figure}
    \includegraphics[width=\linewidth]{Figures/09_DTvT_z0_grid.png}\\
    \caption{The relative variation in temperature measures, $\Delta T= T_{\rm y}-T_{\rm m}$ with respect to the mass-weighted temperature, $T_{\rm m}$. Here each point represents a single halo within the simulation, the solid line represents the medians within logarithmically spaced $T_{\rm m}$ bins, while the dashed lines show the 16 and 84\% percentiles for the same bins. The discrepancy in the scatter for $\textsc{The300}$ around $T\simeq 3$~keV is discussed below in Section~\ref{sec:TS_DeltaT}.}
    \label{fig:09_DTvT_scaling}
\end{figure}
As already mentioned, we find an offset between the temperature measures, which we now explore in more detail. For this, in Fig.~\ref{fig:09_DTvT_scaling}, we plot the fractional temperature difference, $\Delta T/T_{\rm m} = (T_{\rm y}-T_{\rm m})/T_{\rm m}$ against $T_{\rm m}$ itself. We already identified that the mass-weighted temperature is a good proxy for the mass itself, so it is worth noting that this temperature difference has a similar, if subtly different, relationship than when plotted against mass.

We can see immediately that the systematic offset is different in each of the four samples. However, it is interesting to note that this offset is subject to a significant skew -- which is to say, at the simplest level the offset holds on a cluster by cluster basis, and $T_{\rm y}$ is {\it always} greater than $T_{\rm m}$ within clusters.
This can be seen visually, as here we have plotted every halo within our samples in Fig.~\ref{fig:09_DTvT_scaling}, and it is evident that within each cluster, there is a minimum difference, greater than zero, between all the clusters. Taking the 1st percentile for each cluster sample we can estimate this minimum offset as being highest in $\textsc{Magneticum}$ with $\Delta T \geq 0.178\, T_{\rm m}$, and smallest in $\textsc{The300}$ where $\Delta T \geq 0.066\, T_{\rm m}$.

Here, it is worth briefly discussing why the $\textsc{The300}$ sample has such a non-uniform intrinsic scatter compared to the other samples. We note that this divergence is worst at $T_{\rm m} \simeq 3$~keV, reaching convergence again at high temperatures, ($T_{\rm m} \gtrsim 4$~keV). This is an artifact of the selection process for low mass clusters within the $\textsc{The300}$ sample. That is, since all `low mass' ($M_{\rm 200c}\lesssim 0.9\times10^{15}$~$\msol$) clusters, exist in the region of larger haloes, this biases the temperatures of the clusters. As such, this region is less relevant for a mass-complete understanding of the temperature differences here.

\subsection{Resilience of results}
\label{sec:Res}
\begin{figure}
    \includegraphics[width=\linewidth]{Figures/10_red_fgasvM_grid.pdf}\\
    \caption{A comparison of the gas fraction between each simulation. Here the gas fraction, $f_{\rm gas, 200c}=M_{\rm gas,200c}/M_{\rm 200c}$ has been normalised by the relevant $\Omega_{\rm B}$ and $\Omega_{\rm M}$ for each simulation. This figure is arranged as in Fig.~\ref{fig:03_TvM_z0}, with $z=0$ in the top panel and $z=1$ in the lower panel.}
    \label{fig:10_GasFraction_vM}
\end{figure}
All the simulations use different physical models and numerical methods to generate the halo populations we study. Moreover, they are all calibrated to different measurements, as such it is remarkable that we see the agreement we have found across the four different samples. However, in this section we will focus more specifically on the effects these differences have caused on our observed results.

Firstly, we briefly consider the effect of resolution in simulation -- this is presented in more detail in Appendix~\ref{sec:Res_Resolution}. We have found inconclusive results: We used three different resolution runs from $\textsc{TNG}$ where all the other parameters were held the same, and found minimal effects on $T_{\rm m}$ and $T_{\rm y}$. $T_{\rm sl}$ did decrease in value with increasing resolution, most likely due to its dependence on resolving small dense clumps in clusters. However, we also compared two different resolution runs $\textsc{Magneticum}$, where each was independently calibrated (although these results are not plotted in this work). Here the lower resolution run also resulted in lower values for $T_{\rm m}$ and $T_{\rm y}$, with a shift around that of the intrinsic variance within each simulation. Within the context of \citet{Schaye2015}, this indicates that the temperatures show `strong' resolution convergence, but show a diminished `weak' convergence.

Since our temperature measures are dependent on gas density, it is important to consider how the differences in gas fraction (i.e., $f_{\rm gas, 200c} = M_{\rm gas, 200c}/M_{\rm 200c}$.) vary between our simulations. In Fig.~\ref{fig:10_GasFraction_vM}, we show the gas fractions within $R_{\rm 200c}$ (scaled by the cosmic baryonic fraction in the simulations) at $z=0$ and 1. In general, we often consider $f_{\rm gas}$ to be a probe of the feedback used in clusters, with lower values for $f_{\rm gas}$ indicating more effects from feedback. Hence generally speaking, at lower masses, $f_{\rm gas}$ is lower as the potential well for haloes are shallower, allowing for more gas to be ejected from haloes from feedback.

It is evident that there are significant variations in $f_{\rm gas}$ among the simulations considered in this work. At the highest masses, $\textsc{Magneticum}$, $\textsc{The300}$ and $\textsc{Bahamas+Macsis}$ start to agree, with values $>0.9$, but $\textsc{TNG}$ falls greatly below this. Similarly, at the lowest masses, at $z=0$, $\textsc{Magneticum}$, $\textsc{The300}$ and $\textsc{TNG}$ roughly agree, while $\textsc{Bahamas+Macsis}$ is significantly smaller. 

Furthermore, the samples exhibit varying level of redshift evolution. The gas fraction rises generally, indicating that as clusters evolve, feedback results in gas being ejected from haloes, so younger clusters (i.e., in general, clusters at higher redshifts), will have higher values for $f_{\rm gas}$. However, the slopes in each of the samples decrease by differing values -- $\textsc{TNG}$ most notably, with its values reducing dramatically everywhere except the lowest mass haloes.

This diverging behaviour in $f_{\rm gas, 200c}$, lies in contrast to the strong agreement we find between the samples for both the $T_{\rm m}-$mass relations, and the $Y-$mass relations, although it is indicated by the $T-Y$ variation. Since $Y$ is a measure of the gas pressure within clusters, we could expect it to be self-calibrating -- that is, a cluster in hydrostatic equilibrium, would lead to a certain gas pressure to counter the gravitational well of its own mass. As such we would expect haloes to have a strong $Y-M$ relationship.
 
However, we can see a reflection of the $f_{\rm gas}$ variation subtly in $T_{\rm y}$. That is, we can see that the the curvature of the $T_{\rm y}-$mass relation as seen in, e.g., Fig.~\ref{fig:03_TvM_z0}, is similar (albeit inversely) to the curvature in the $f_{\rm gas}-$mass relation. In fact, we can find that this is driven by the changes in $M_{\rm gas}$ as $M_{\rm gas} T_{\rm y}$ forms a tight relation with $Y$ and a slightly more consistent relation (than $T_{\rm y}$ alone) with mass. This is perhaps unsurprising, as we already know that on a cluster by cluster level, $T_{\rm y}$ and $T_{\rm m}$ form a strong relationship, and since $M_{\rm gas} T_{\rm m} \propto Y$, we can consider $M_{\rm gas} T_{\rm y}$ as a kind of relativistic equivalent to $Y$. None the less, this means a more detailed understanding of $M_{\rm gas}$ (or equivalently $f_{\rm gas}$) in clusters would lead to more assurance in the exact $T_{\rm y}-$mass and $T_{\rm y}-Y$ relationships.

In Appendix~\ref{sec:Res_Feedback}, we compare three different feedback model runs used in the $\textsc{The300}$ project. In general, we find that $T_{\rm m}$ is barely affected by the different models, only perturbed in regions where the feedback models are extreme. In contrast, $T_{\rm y}$ matches our suppositions from studying $f_{\rm gas}$, where again $T_{\rm y}$ is largely consistent across the different feedback models, but shows variations inversely related to the variation of $f_{\rm gas}$. That is to say, stronger feedback in general leads to higher values of $T_{\rm y}$, but variations in the specific implementation of feedback lead to more complicated effects. $T_{\rm sl}$ varies significantly in amplitude and gradient between the three feedback runs, making it less reliable as a temperature proxy.

\section{Cross simulation averaged results}
\label{sec:cross_sim_fits}
Due to the broad agreement of all of the samples in each temperature, we can consider the effects of averaging across simulations to obtain `simulation independent' predictions for these values from simulations. Here we obtain these by joining all the samples to form one large population of halos, which we then sort into mass bins as before, and find sample fits to the mass-binned averages.\footnote{This approach obtains results consistent with those obtained by, for instance, taking the mass-binned values for each sample, and joining these together.} We note that this means our averages will be weighted more by those simulations with larger populations of halos. We also provide fits to individual simulations in Table \ref{apptab:fits_median_z0_samples}, and the variations in simulation predictions can be easily used to estimate the theoretical, astrophysics-driven uncertainty in the mean temperature-mass relations, for example via the method described in \citet[see Section 4.3]{Anbajagane2022}.

When fitting this data to obtain temperature-mass relations, we then consider both a two- and three-parameter fit of the form
\begin{equation} \label{eqn:fit_format}
    T = E(z)^{2/3} A \left(\frac{M}{10^{14}\;\msol}\right)^{B+C\log_{10}(M/10^{14}\msol)}\,{\rm keV},
\end{equation}
where for the two-parameter fit we set $C=0$. The cross-simulation averaged results at $z=0$ for all three temperatures can be found in Table~\ref{tab:fits_median_z0_all} -- the fits for each individual sample at $z=0$ can be found in Table~\ref{apptab:fits_median_z0_samples}. The errors are obtained through bootstrap techniques and show the error within the mean. Here we have also calculated a measure of the scatter through the root mean squared dispersion around these mean values as
\begin{equation} \label{eqn:sigma_def}
    \sigma_{\log_{10}T}=\sqrt{\frac{1}{N}\sum_i\left[\log_{10}\left(\frac{T_i}{T_{\rm fit}(M_i)}\right)\right]^2},
\end{equation}
with $i$ indexing over all the halos at a given redshift. We note that this measure is weighted more by the lower mass clusters and groups, due to the larger number of them in each sample than higher mass haloes. This runs opposite to our fits where we have minimised this bias by fitting to the averages gained from a selection of mass bins. 

\begin{table}
\caption{The two and three parameter fits [Eq.~\eqref{eqn:fit_format}] and a measure of the intrinsic scatter [Eq.~\eqref{eqn:sigma_def}] in the fits against mass within $R_{\rm 200c}$ at $z=0$ for the cross simulation averaged sample. The errors here are the errors in the fits found through bootstrapping.} \centering
\begin{tabular}{lcccc}\smallskip
 & $A$ & $B$ & $C$ &$\sigma_{\log_{10}T}$ \\ 
 \hline \smallskip
 $T_{\rm y}$  & $1.465^{+0.002}_{-0.002}$ & $0.586^{+0.003}_{-0.002}$ & 0 & $0.1025^{+0.0003}_{-0.0002}$\\ \smallskip
 $T_{\rm y}$  & $1.426^{+0.006}_{-0.007}$ & $0.566^{+0.001}_{-0.001}$ & $0.024^{+0.005}_{-0.004}$ & $0.1011^{+0.0001}_{-0.0001}$\\ 
 \hline \smallskip
 $T_{\rm m}$  & $1.210^{+0.001}_{-0.001}$ & $0.591^{+0.003}_{-0.003}$ & 0 & $0.0805^{+0.0002}_{-0.0001}$\\ \smallskip
 $T_{\rm m}$  & $1.207^{+0.005}_{-0.005}$ & $0.589^{+0.001}_{-0.001}$ & $0.003^{+0.004}_{-0.005}$ & $0.0804^{+0.0001}_{-0.0001}$\\
 \hline \smallskip
 $T_{\rm sl}$ & $1.135^{+0.003}_{-0.003}$ & $0.601^{+0.006}_{-0.011}$ & 0 & $0.2067^{+0.0016}_{-0.0008}$ \\ \smallskip
 $T_{\rm sl}$ & $1.196^{+0.020}_{-0.009}$ & $0.641^{+0.003}_{-0.003}$ & $-0.048^{+0.007}_{-0.018}$ & $0.2028^{+0.0002}_{-0.0001}$\\
 \hline
\end{tabular}
\label{tab:fits_median_z0_all}
\end{table}

The first aspect we can note is the exceptional lack of curvature in the $T_{\rm m}-$ Mass relationship, where even when we allow for curvature, $C$, is fully consistent with zero, albeit with a slight tendency towards positive curvature. This is true even if we vary the pivot mass of the fit. We also see that $T_{\rm y}$ has definitive positive curvature, as follows from our early discussion. On the other hand, $T_{\rm sl}$ shows some significant negative curvature, but this may be an artifact of the lower mass haloes' scatter, and may not be representative of the behaviour of observed X-ray temperatures.

It is also important to note that the intrinsic variance is slightly larger in the combined sample than it is in each simulation, due to small bulk offsets between each sample. This is {\it not} contained within the errors tabulated in Table~\ref{tab:fits_median_z0_all}, which show the errors in the central fit parameter. At $z=0$, the intrinsic variance can be found to be around $\pm6.0\%$ in $T_{\rm m}$, $\pm7\%$ for $T_{\rm y}$ and around $\pm10\%$ in $T_{\rm sl}$. However, a more detailed understanding of the intrinsic variance can be found by looking at the cross-simulation fits to the 16 and 84 percentile lines which can be found in Tables~\ref{apptab:fits_median_z0_all} and \ref{apptab:fits_3_median_z0_all}.

In general, it is also important to note that there is still a fixed offset between $T_{\rm y}$ and $T_{\rm m}$ to be found in the combined sample -- we find a median offset of around 22\% (see Table~\ref{tab:R_err_all}). However, as with the discussion in Section~\ref{sec:TS_DeltaT}, we can consider the minimal offset to be determined by the 1\% percentile line -- in this case giving a minimum offset of $8.4\%$. This large variance in $T_{\rm y}/T_{\rm m}$ is of course, largely driven by the slight disagreement between clusters, which could potentially be broken with future measurements of $f_{\rm gas}$.

We also use this cross-simulation sample to calculate the redshift evolution of these quantities, and the radial dependence of these results. The first order redshift and radial evolution can be seen in Tables~\ref{tab:fits_redshift_offset_all} and \ref{tab:R_err_all} akin to those discussed in Section~\ref{sec:TS}. The mean 2-parameter fits can be found in Table~\ref{apptab:fits_median_Redshift_all}, alongside the full set of 2- and 3-parameter fits in Tables~\ref{apptab:fits_R5_median_z0_all} and \ref{apptab:fits_3_R5_median_z0_all} for the three temperatures within $R_{\rm 500c}$ at $z=0$.

\begin{table}
\caption{The median correction to the amplitude, $A$, alone of the two parameter fit at each redshift for the temperature measures in the cross-simulation averaged simulation. That is, we are tabulating the averaged value $A^\star$, so that $E(z)^{-2/3} T_z = A^\star T_{z=0}$ The full two parameter fits for each redshift can be found in Table~\ref{apptab:fits_median_Redshift_all}.} \centering
\begin{tabular}{lccc}
 &  $E(z)^{-2/3}T_{\rm y}$ &  $E(z)^{-2/3}T_{\rm m}$ & $E(z)^{-2/3}T_{\rm sl}$ \\ 
 \hline \smallskip
 $z=0.00$ & 1.000 & 1.000 & 1.000 \\ \smallskip
 $z=0.25$ & $0.982^{+0.002}_{-0.002}$ & $0.979^{+0.002}_{-0.001}$ & $0.960^{+0.017}_{-0.011}$ \\ \smallskip
 $z=0.50$ & $0.976^{+0.002}_{-0.002}$ & $0.952^{+0.001}_{-0.001}$ & $0.929^{+0.004}_{-0.007}$ \\ \smallskip
 $z=1.00$ & $0.944^{+0.002}_{-0.002}$ & $0.897^{+0.002}_{-0.002}$ & $0.837^{+0.007}_{-0.011}$ \\ \smallskip
 $z=1.50$ & $0.919^{+0.002}_{-0.003}$ & $0.851^{+0.002}_{-0.002}$ & $0.749^{+0.004}_{-0.004}$ \\
 \hline
\end{tabular}
\label{tab:fits_redshift_offset_all}
\end{table}
\begin{table}
\caption{The shifts and errors in the radius, mass, $T_{\rm m}$ and $T_{\rm y}$ over the combined sample against those values within $R_{\rm 200c}$. These values are calculated on a cluster-by-cluster basis, and then the averages are found within these. The central value here is the median with the errors given by the 16 and 84 percentiles.} 
\centering
\begin{tabular}{lcccc}\smallskip
 & $R/R_{\rm 200c}$ & $M/M_{\rm 200c}$ & $T_{\rm m}/T_{\rm m,200c}$ & $T_{\rm y}/T_{\rm m,200c}$ \\ 
 \hline \smallskip
 $R_{\rm 500c}$  & $0.66^{+0.01}_{-0.02}$ & $0.71^{+0.05}_{-0.07}$ & $1.20^{+0.04}_{-0.07}$ & $1.40^{+0.16}_{-0.12}$\\  \smallskip 
 $R_{\rm 200c}$  & $1.00$ & $1.00$ & $1.00$ & $1.22^{+0.11}_{-0.07}$\\  \smallskip 
 $R_{\rm 500m}$  & $1.11^{+0.01}_{-0.01}$ & $1.08^{+0.03}_{-0.02}$ & $0.95^{+0.01}_{-0.01}$ & $1.18^{+0.10}_{-0.06}$\\  \smallskip 
 $R_{\rm vir}$ & $1.33^{+0.03}_{-0.02}$ & $1.22^{+0.08}_{-0.05}$ & $0.87^{+0.04}_{-0.02}$ & $1.11^{+0.09}_{-0.06}$\\  \smallskip 
 $R_{\rm 200m}$ & $1.64^{+0.06}_{-0.04}$ & $1.39^{+0.15}_{-0.10}$ & $0.79^{+0.07}_{-0.04}$ & $1.05^{+0.09}_{-0.06}$\\
 \hline
\end{tabular}
\label{tab:R_err_all}
\end{table}

We can see in Table~\ref{tab:fits_redshift_offset_all} that the sample averaged redshifts show the same variation we expected from considering each sample independently. That is, all of the temperatures diverge from self similarity, with $T_{\rm sl}$ showing the greatest departure, and $T_{\rm y}$ the smallest. However, when looking at high redshift haloes, some departure is none the less to be expected. Again, we reiterate that the errors here are the errors in the fit, and as such do not encapsulate any variation in the intrinsic scatter from combining samples. It is also interesting to note that, although not shown here, the $\sigma_{\log_{10}T}$ values for both these 1-parameter corrections and the 2-parameter fits that are tabulated in \ref{apptab:fits_median_Redshift_all} are broadly the same. Although, as previously noted, these will be largely weighted by the lower mass haloes, this does indicate that most of the evolution is in the amplitude of the scaling relations, and not in their power law.

This redshift evolution is also not linear with respect to redshift, changing faster at lower redshifts. However, we can find a parameterisation of this redshift evolution in the form of 
\begin{equation}
    \log_{10}(A^\star) = p \log_{10}(1+z)+ q \log^2_{10}(1+z),
\end{equation}
where $A^\star$ is the redshift correction factor, so $T_z = A^\star E(z)^{2/3} T_{z=0}$. We find for $T_{\rm y}$, $T_{\rm m}$ and $T_{\rm sl}$ that the values for $[p,q]$ are $[-0.05,-0.11]$, $[-0.08,-0.24]$ and $[-0.09, -0.57]$, respectively. However, it is worth noting that at higher redshifts, $T_{\rm sl}$ becomes an increasingly poor proxy for $T_{\rm X}$ as the number of haloes we have at temperatures above 3.5~keV diminishes. 

The radial corrections tabulated in Table~\ref{tab:R_err_all} {\it show} the effects of combining the samples, where the errors represent the intrinsic variance. We can see that the combined sample only slightly increases this intrinsic radial variance. While the effects of changing the radius are more complex than a single number can fully encapsulate, this is still useful for an indicative understanding of the effects on viewing clusters through different apertures. However, it is interesting to note that the value of $T_{\rm y}/T_{\rm m}$ at each radius, increases at higher radii. That is, as we increase the radius of interest, $T_{\rm m}$ decreases faster than $T_{\rm y}$ or equivalently, $T_{\rm y}$ has a shallower profile. The equivalent values for each sample to those in Table~\ref{tab:R_err_all}, can be found in Table~\ref{apptab:R_err_sample}.

\begin{figure}
    \includegraphics[width=\linewidth]{Figures/15_TvM_obs_grid_M500.pdf}\\
    \caption{A comparison of the cross simulation-averaged temperature scaling relations, with observational results. Note these quantities are all measured within $R_{\rm 500c}$. Here the observational lines are plotted only within the mass range of each study. Here A05, V09, M16 and B19 refer to the fits from \citet{Arnaud2005,Vikhlinin2009xrayobs,Mantz2016} and \citet{Bulbul2019}, respectively. A05 and V09 both depend on the hydrostatic mass bias, for which we here assumed $1-b=0.8$.}
    \label{fig:15_observations}
\end{figure}
\subsection{Comparison to X-ray observations}
There have been many observational studies aiming to constrain the X-ray temperature, $T_{\rm X}$, mass relationship. Here, we will consider 4 such studies to compare to our cross-simulation averaged results. In particular, we have considered \citet[A05]{Arnaud2005}, \citet[V09]{Vikhlinin2009xrayobs}, \citet[M16]{Mantz2016} and \citet[B19]{Bulbul2019}, all of which provide 2 parameter fits for $T_{\rm X}$-$M_{\rm 500c}$ relation. 

A05 uses 10 low redshift ($z<0.15$) clusters from {\it XMM-Newton}, and obtains masses through fitting an NFW-type profile assuming hydrostatic equilibrium. Here we use their fits for the subsample of 6 hot ($>3.5$~keV) clusters.
V09 uses {\it Chandra} measurements of 85 clusters (over redshifts of $z\lesssim0.9$), obtaining through fits of $\beta$-profiles, again assuming hydrostatic equilibrium. As such, both A05 and V09's masses are offset from the `true' mass according to the hydrostatic mass bias. This has been measured many times, both in simulations and using observational techniques, and in this work we will use $(1-b)=0.8$.
M16 uses 40 hot ($>5$~keV), relaxed clusters observed by {\it Chandra} (with redshifts $z\lesssim 1$). This work calibrates its masses using weak lensing measurements, so can be taken to be measures of the true mass.
Finally, B19 uses 59 clusters from {\it XMM-Newton}, over redshifts $0.2<z<1.5$. The masses are obtained through SZ measurements by the South Pole Telescope. This work also provides two different fits for $T_{\rm X}$ and here we use their core-excised fits, to compare best with our other measures.

In Fig.~\ref{fig:15_observations}, we show these 4 fits against our whole sample fits, here all within $R_{\rm 500c}$. The length of each line matches the mass range of the data sets used within each observation. It is first worth noticing, that there is no strong consensus between X-ray observations about the details of the $T_{\rm X}$-$M$ relationship. This may be exacerbated by the different techniques used in each study. For B19 and M16 an intrinsic scatter of around 13\% is given, while V09 obtains a slightly higher 15-20\%. These are calculated as $\sigma_{{\rm ln}T}$ and lie a little smaller than those values predicted from our combined simulations.

Here A05 and B19 both use data sets from {\it XMM-Newton} and agree with each other best, as well as holding comparable results to the whole-sample $T_{\rm sl}$ measurement. However, B19 has a far steeper gradient (at $B=0.80^{+0.11}_{-0.08})$ than any of the other studies, which is at a huge contrast to all of our simulation fits which indicate gradients always less than self-similarity (that is $B<2/3$). The other three studies, also all have gradients steeper than that obtained through our simulations, all of which are consistent with self-similarity. It is interesting to note, that the full sample for A05, (rather than the hot sample displayed here) does also show a shallower gradient more compatible with our simulation results.

M16 shows the most extreme result, tending towards far higher temperatures than seen in the other scalings, and higher even than the $T_{\rm y}$ values we have found within simulations. This may be a side effect of the sample selection towards hot clusters, pushing the average up, or may indicate the intrinsic spread within the $T_{\rm X}$ relationship. V09, the other measurement using {\it Chandra} observations, also lies higher than the {\it XMM-Newton} values, indicating an X-ray temperature more consistent with our simulation $T_{\rm m}$ values, than $T_{\rm sl}$.

It is still unclear whether the variation observed here between the observational methods and simulations, indicates that simulations are not capturing some facet of real galaxy clusters, or are miscalibrated, or if this comes from the difficulties in obtaining `true' masses from observations. A factor which could be exacerbated by the suggested intrinsic spread in the $T_{\rm X}$ relationship and the comparatively (when compared to the data sets from simulations) small data sets and mass ranges used within observational studies.
Either way, the $T_{\rm y}$ temperature is not directly accessible by any of these X-ray studies and we thus instead recommend using the scalings derived here in future modeling of SZ observables. 

\section{Discussion} 
\label{sec:applications} 
\subsection{$T_{\rm y}$--$Y$ self-calibration}
\label{sec:self-calibration}

The $Y$--$T_{\rm y}$ relation can be used in the SZ signal modeling once the $Y$-parameter is determined. This allows one to refine the SZ model (with the relativistic SZ corrections) even if the data for individual systems cannot provide a direct constraint on $T_{\rm y}$. The observed Y can also differ from the true Y due to angular resolution (caused by the instrument beam) but when the resolution is $\lesssim$ 1 arcmin, these smoothing effects in $Y$ are negligible \citep[][]{Yang2022}, and so the true $Y$--$T_{\rm y}$ relation presented here can be used directly on observed data without needing any further processing. Similarly, one can easily use the $Y$--$T_{\rm y}$ relation in simulations of the SZ sky, where from the simulation the cluster temperature is not directly available, e.g., in the {\tt WebSky} \citep{Stein2020}.

Here we will briefly discuss the scaling relationship we have determined, alongside the effects we may be able to determine with these results. As noted, we find some complexity in this relationship, with varying levels of curvature, especially at small values of $Y$. As such, we have tabulated the 2- and 3-parameter fits for the $T_{\rm m}$-$Y$ and $T_{\rm y}$-$Y$ relationships in Table~\ref{apptab:fits_Y_percentiles_z0_all}. The equivalents for the medians for each sample can also be found in Table~\ref{apptab:fits_Y_median_z0_sample}.

\begin{table}
\caption{The two parameter $T_{\rm y}$-$Y$ fits [Eq.~\eqref{eqn:Yfit_format}] for haloes with $Y_{\rm 200c}>10^{-6}$~Mpc, and a measure of the intrinsic scatter [Eq.~\eqref{eqn:sigma_def}] in the fits across redshifts for the median of the cross simulation averaged sample within $R_{\rm 200c}$.} \centering
\begin{tabular}{lccc}
 & $A$ & $B$ & $\sigma_{\log_{10}T}$ \\ 
 \hline \smallskip
 $z=0.00$ & $2.614^{+0.006}_{-0.006}$ & $0.368^{+0.003}_{-0.010}$ & $0.1874^{+0.0039}_{-0.0152}$ \\ \smallskip
 $z=0.25$ & $2.593^{+0.007}_{-0.007}$ & $0.372^{+0.002}_{-0.005}$ & $0.1685^{+0.0036}_{-0.0074}$ \\ \smallskip
 $z=0.50$ & $2.593^{+0.008}_{-0.008}$ & $0.381^{+0.003}_{-0.002}$ & $0.1625^{+0.0034}_{-0.0033}$ \\ \smallskip
 $z=1.00$ & $2.585^{+0.013}_{-0.012}$ & $0.382^{+0.004}_{-0.004}$ & $0.1388^{+0.0044}_{-0.0045}$ \\ \smallskip
 $z=1.50$ & $2.597^{+0.052}_{-0.031}$ & $0.373^{+0.020}_{-0.011}$ & $0.1263^{+0.0185}_{-0.0091}$ \\
 \hline
\end{tabular}
\label{tab:fits_Yhot_median_redshift_all}
\end{table}

We can, however, create more stable 2-parameter fits if we exclude the smallest halos with $Y<10^{-6}$, which are unlikely to meaningfully contribute to SZ observations, and which, even when they can be observed, will have the smallest rSZ corrections. With this restriction and 
\begin{equation} \label{eqn:Yfit_format}
    T = E(z)^{2/5} A \left(\frac{Y}{10^{-5}\;\msol}\right)^{B}\,{\rm keV},
\end{equation}
we can find fits for the combined sample at each redshift as given in Table~\ref{tab:fits_Yhot_median_redshift_all} (a full form including the 16 and 84 percentiles can be found in Table~\ref{apptab:fits_Yhot_percentiles_redshift_all}). We can immediately see that there is little redshift dependence beyond the expected self-similar evolution.

We note that even here we have more intrinsic variation in this relationship than within the $T_{\rm y}-M$ relationships; however, it can still be used as a reliable proxy to estimate the relativistic effects to haloes. Moreover, as our understanding of $f_{\rm gas}$ in clusters improves, we may be able to reduce this uncertainty between clusters, as this would allow a greater understanding of the true $M_{\rm gas}$ values we should expect for haloes, and thus a more precise idea of exactly how we expect the $T_{\rm y}-Y$ relationship to behave.

With those comments in mind, we can immediately gain a sense of the corrections we would expect for a given cluster. For instance, from our simulations we can estimate that a cluster with a true $Y_{\rm 200c} = 10^{-4}$ at $z=0$, would have a temperature of $6.1^{+0.8}_{-1.1}$~keV, where the errors here are driven almost entirely by the disagreement between the simulations. This would lead to a fractional change in the amplitude at 353~GHz of $0.907^{+0.016}_{-0.011}$ -- that is, a 10\% underestimation of $Y$ by using the non-relativistic approximation.

It is also important to remember that the relativistic corrections to the kinetic SZ (kSZ) effect, will be driven by $T_{\rm m}$ and not $T_{\rm y}$. We have not discussed the kSZ effect in much detail in this work, but nonetheless detailed studies of kSZ signals, should also consider the impact that this signal may have. These relativistic kSZ effects, are, however, broadly speaking an order of magnitude smaller that the relativistic tSZ correction \citep[e.g.,][]{Sazonov1998, Nozawa2006, Chluba2012}. Still, the relations given here can be used to model the effects.

\subsection{Applications to current and future SZ analyses}
\label{sec:Dis_effects}
In this section, we briefly discuss where the obtained scaling relations may have applications in current and future studies of the SZ effect. First and foremost, due to lack of sensitivity or spectral coverage, current SZ measurements are still not directly sensitive to the rSZ effect. Nevertheless, rSZ can already affect the inference of cosmological parameters if it is neglected in the modeling. One important example is related to the SZ power spectrum analysis, for which rSZ can lead to an underestimating of the $yy$-power spectrum and thus cause a systematic shift in the inferred value of $\sigma_8$ \citep{Remazeilles2019}. They would similarly impact cross-correlations of the SZ field with large-scale structure fields such as galaxy positions or cosmic shear, and such measurements have been recently used to infer physics like the redshift-dependent mean thermal pressure of the Universe as well as the energetics of feedback in groups and clusters \citep[e.g.,][]{Osato2018, Osato2020, Pandey2019, Pandey2021, Gatti2021}.

The size of the rSZ bias depends directly on the experimental configuration. Similarly, different experimental setups are more or less prone to mis-modeling of rSZ. For example, many of the planned CMB experiments have channels at $\nu \lesssim 220\,{\rm GHz}$. In this case, the degeneracy of the SZ signal with respect to the $y$-parameter and $T_{\rm y}$ cannot be easily broken in the observation. However, our new scaling relations will allow us to accurately include the rSZ corrections in the theoretical modeling of the data.

Similar to SZ power-spectrum analyses, the SZ cluster number count analyses is expected to be affected by the presence of rSZ. Here, two aspects are important: at a given mass, the SZ flux is diminished due to rSZ. This means that clusters are assigned to a lower signal-to-noise bin. In addition, for the extraction of the clusters, the multi-match-filtering (MMF) method \citep{Melin2006} is not optimally tuned to the correct spectral shape, leading to a mis-estimation of the noise. The scaling-relations introduced here can be directly used to inform the MMF and eliminate the impact rSZ may have. In particular the $T_{\rm y}-Y$ relation will be of significance here, since it allows to construct an iterative MMF approach \citep[e.g.,][]{Zubeldia2022} to incorporate the rSZ effect based solely on SZ observables. This should allow for a more robust comparison of theory and observation and also improve the constraining power of the obtained SZ catalogs.

The rSZ also impacts studies on the thermodynamics of cluster outskirts ($R \gtrsim R_{\rm 200m}$) --- especially non-thermal features like cosmological accretion shocks \citep{Aung2021, Baxter2021ShocksSZ} --- which contain astrophysical and cosmological information \citep[see][for a review]{Walker2019} and have only recently been observationally explored. \citet{Anbajagane2022c} performed the first large population-level analysis of tSZ profile outskirts and found signs of cosmological shocks, which manifest as decrements in the profile. Others have seen similar signs using samples of up to ten clusters \citep[e.g.,][]{Hurier2019ShocksSZPlanck, Pratt2021ShocksPlanck}. However, as was described above, the rSZ effect \textit{also} causes a tSZ signal reduction. Thus, a better understanding of the rSZ, and the magnitude of the tSZ decrement it causes, will be needed to robustly infer the non-thermal \citep{Nelson2014b,Aung2021} and plasma physics in these outskirts \citep{Rudd2009,Avestruz2015}. 

As we have seen in Fig.~\ref{fig:04_TvM_redshift}, $T_{\rm y}$ has a scaling relation, close to self-similar evolution with redshift. However, nonetheless significant redshift evolution is still expected from self-similarity, and the minor corrections we observe. This effect can have important consequences for cosmological inferences relying on the redshift-independence of the SZ effect. One example is attempts to use the SZ effect to measure the CMB temperature-redshift relation \citep{Rephaeli1980, Luzzi2009}, another relates to SZ-derived values of the Hubble constant \citep{Birkinshaw1991, Mauskopf2000, Wan2021}. Similarly, applications of SZ measurements to constrain possible time-variations of the fine structure constant \citep{Bora2021} are prone to the redshift-dependence of rSZ effect. With the derived $T_{\rm y}-Y$ relations, the effect of rSZ on these inferences can be readily marginalized over. 

As a final application of rSZ we mention predictions of the all-sky averaged SZ and rSZ effects \citep{Hill2015}. This signal is one of the targets of future CMB spectroscopy \citep{Chluba2019} and can inform us about feedback processes in cluster physics \citep[e.g.,][]{Thiele2022}. The temperature relations obtained here allow us to refine these prediction, giving a simulation-averaged view on the expected signal. In a similar manner, the relations can be used to refine the calculations of the radio \citep{Holder2022, Lee2022} and CIB SZ \citep{Sabyr2022, Acharya2022} effects. These signals might become important targets for future radio and sub-mm observations, allowing us to probe the evolution and origin of the cosmic radio and infra-red backgrounds.

\section{Conclusions}
\label{sec:Conclusion}
In this work, we present detailed comparisons of three cluster temperature measures: (a) the average rSZ temperature, (b) the mass-weighted temperature relevant for the thermal SZ (tSZ) effect, and (c) X-ray spectrosopic temperature using the $\textsc{Bahamas}$ \& $\textsc{Macsis}$, {\sc Illustris-TNG}, $\textsc{Magneticum}$, and {\sc The Three Hundred Project} simulations. We analyze gas temperature scaling relations of galaxy groups and clusters with $M_{\rm 500c}>10^{13}$~$\msol$, over 5 redshifts between $z=0$ to $z=1.5$. We provide fits to multiple scaling relations for individual simulations and for a combined cross-simulation sample, with the former ensemble of results also providing an estimate of the theoretical, astrophysics-driven uncertainty in the relations.
Our main results are summarized as follows:
\begin{itemize}
    \item There is exceptionally strong agreement for $T_{\rm m}$ between all four simulations with $M_{\rm 200c}$. $T_{\rm y}$ is consistently larger than $T_{\rm m}$ (by an average of $\simeq22\%$), which is generally a little above $T_{\rm sl}$ at $z=0$. $T_{\rm y}$ has a good agreement between the simulations, however, there is variation in the exact magnitude of the offset between $T_{\rm y}$ and $T_{\rm m}$ between simulations. $T_{\rm sl}$ also shows agreement, although it is also subject to a great deal more intrinsic scatter than the other two temperatures. All three temperature measures exhibit different mass scalings and vary differently with redshift and radius.
    \item All three temperature measures exhibit deviation from the self-similar evolution. At higher redshifts, they all fall below the expected $E(z)^{-2/3}$ scaling indicating all haloes will have lower temperatures at higher redshifts than the self-similar model prediction. However, $T_{\rm y}$ evolves very closely to self-similarity, while the other temperature measures departs further from self-similarity, so that at higher redshifts, all three temperature measures increasingly diverge. $T_{\rm y}$ has an increasing correction relative to $T_{\rm m}$ and an even larger correction to $T_{\rm sl}$.
    \item The temperature measures all agree best within $R_{\rm 200c}$ between simulations. Each simulation and temperature measure has a different radial profile leading to varied results. While the temperature measures still agree well within $R_{\rm 500c}$, the improvement at $R_{\rm 200c}$ indicates this may indeed be an optimal radius to study SZ science.
    \item The gas fraction, feedback methods and resolution all vary significantly between the simulations. In light of this, the level of agreement we see is startling, and indicates that much of the SZ gas physics is sufficiently calibrated by the microphysical constraints (i.e., stellar properties). However, when we examine these in more detail, we can find that the gas fraction is correlated with the variation in $T_{\rm y}$ between simulations. As such, if this can be determined with more accuracy in future observations, the strength of these predictions may increase. In general, when we study resolution and feedback within equivalently calibrated simulations, we find little variation to our observed SZ temperature measures, while $T_{\rm sl}$ is affected slightly more.
    \item We have created a cross-simulation sample and found the fitted values for our temperatures. In general, $T_{\rm m}$ shows limited tendency towards curvature (i.e., a mass dependent slope), while $T_{\rm y}$ has positive curvature, and $T_{\rm sl}$ negative curvature. We have provided a simple regime for calculated the redshift corrections to temperatures within our redshift range, and clarified the broad effects of varying the radial aperture which we use to define haloes. When compared to observational results, we find that the temperatures broadly agree. However, there is more variation within X-ray results than our predictions, making it difficult to draw out strong conclusions. Nonetheless, observations all suggest steeper scaling relations than we have found in the simulations.
    \item In general while these temperatures will be difficult to directly measure, they give rise to the possibility of self-calibrating SZ observations, to allow for relativistic corrections to be used within the determination of SZ measurements themselves. Our simulations suggest that for a halo with a true $Y_{\rm 200c}\simeq 10^{-4}$, we would measure a 10\% underestimation of $Y$ by neglecting relativistic effects. We provide $T_{\rm y}-Y$ relations to allow for further detailed modelling.
\end{itemize}

Future and ongoing experiments such as CCAT-Prime \citep{Stacey2018}, NICA2 \citep{NIKA2}, TolTEC \citep{TolTEC}, and The Simons Observatory \citep{TheSimonsObservatoryCollaboration2018} offer an exciting potential for measuring the ICM temperature using the rSZ effect. This will enable comprehensive analyses of the ICM structure and evolution, especially for high-redshift clusters where X-ray temperature are difficult to obtain. These will also provide an observational test for the validity of the simulation results presented here. Where rSZ temperatures cannot be directly measured, we have provided temperature scaling relations that can be used widely to estimate the impact and potential constraining power of the relativistic SZ effects in future measurements.


{\small
\section*{Acknowledgements}
The authors would like to thank Ian McCarthy for use of the $\textsc{Bahamas}$ and $\textsc{Macsis}$ simulations.
This work was supported by the ERC Consolidator Grant {\it CMBSPEC} (No.~725456).
PS gratefully acknowledges support from the YCAA Prize Postdoctoral Fellowship and helpful discussions with Antonio Ragagnin, Alex Saro and Veronica Biffi.
EL was supported by the Royal Society on grant No.~RGF/EA/180053.
DA is supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE 1746045.
JC was furthermore supported by the Royal Society as a Royal Society University Research Fellow at the University of Manchester, UK (No.~URF/R/191023).
WC is supported by the STFC AGP Grant ST/V000594/1 and the Atracci\'{o}n de Talento Contract no. 2020-T1/TIC-19882 granted by the Comunidad de Madrid in Spain. He further acknowledges the science research grants from the China Manned Space Project with NO. CMS-CSST-2021-A01 and CMS-CSST-2021-B01.
KD acknowledges support by the COMPLEX project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679 as well as support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC-2094 - 390783311. 
GY acknowledges financial support from the MICIU/FEDER (Spain) under project grant PGC2018-094975-C21. WC and GY would like to thank Ministerio de Ciencia e Innovacion (Spain) for financial support under research grant  PID2021-122603NB-C21

This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure
The \textsc{Magneticum} simulations were carried out at the Leibniz Supercomputer Center (LRZ) under the project pr83li.
The Illustris$\textsc{TNG}$ simulations were run with compute time granted by the Gauss Centre for Supercomputing (GCS) under Large-Scale Projects GCS-ILLU and GCS-DWAR on the GCS share of the supercomputer Hazel Hen at the High Performance Computing Center Stuttgart (HLRS). 
$\textsc{The300}$ project has received financial support from the European Union’s Horizon 2020 Research and Innovation programme under the Marie Sklodowskaw-Curie grant agreement number 734374, i.e. the LACEGAL project. We would like to thank The Red Espa{\~n}ola de Supercomputaci{\'o}n for granting us computing time at the MareNostrum Supercomputer of the BSC-CNS where most of the 300 cluster simulations have been performed. The MDPL2 simulation has been performed at LRZ Munich within the project pr87yi. The CosmoSim database (https://www.cosmosim.org) is a service by the Leibniz Institute for Astrophysics Potsdam (AIP). Part of the computations with \texttt{GADGET-X} have also been performed at the ‘Leibniz-Rechenzentrum’ with CPU time assigned to the Project ‘pr83li’.
}
\section*{Data Availability}
Public data releases exist for \textsc{IllustrisTNG} \citep{Nelson2019} and \textsc{Magneticum} \citep{Ragagnin2017}, and can be found at {\tt https://www.tng-project.org/} and {\tt http://magneticum.org/data.html} respectively. The thermodynamic properties presented in this work are not part of the public catalogs but can be provided on reasonable request. 
The data for {\sc The300}, {\sc Bahamas} and {\sc Macsis} are not available in a public repository, but can be provided on request.

{\small
\bibliographystyle{mnras}
