# Copacabana: A Probabilistic Membership Assignment Method for Galaxy Clusters

J. H. Esteves1, M. E. S. Pereira,2 M. Soares-Santos,1 J. Annis,3 A. Farahi,4
F. Andrade-Oliveira,1 P. Barchi,5 A. Palmese,6 H. Lin,3 B. Welch,7,8 H.-Y.
Wu,9 M. Aguena,10 O. Alves,1 D. Bacon,11 S. Bocquet,12 D. Brooks,13 A. Carnero
Rosell,14,10,15 J. Carretero,16 M. Costanzi,17,18,19 L. N. da Costa,10 J. De
Vicente,20 P. Doel,13 S. Everett,21 B. Flaugher,3 J. Frieman,3,22 J. García-
Bellido,23 D. Gruen,12 R. A. Gruendl,24,25 G. Gutierrez,3 S. R. Hinton,26 D.
L. Hollowood,27 K. Honscheid,28,29 D. J. James,30 K. Kuehn,31,32 C.
Lidman,33,34 M. Lima,35,10 J. L. Marshall,36 J. Mena-Fernández,37 R.
Miquel,38,16 J. Myles,39 R. L. C. Ogando,40 A. Pieres,10,40 A. A. Plazas
Malagón,41,42 A. K. Romer,43 E. Sanchez,20 D. Sanchez Cid,20 B. Santiago,44,10
M. Schubnell,1 I. Sevilla-Noarbe,20 M. Smith,45 E. Suchyta,46 M. E. C.
Swanson,24 N. Weaverdyck,1,47 P. Wiseman,45 and M. Yamamoto48 (DES
Collaboration)
_(Affiliations are listed at the end of paper)_ Contact e-mail:
<EMAIL_ADDRESS>
###### Abstract

Cosmological analyses using galaxy clusters in optical/NIR photometric surveys
require robust characterization of their galaxy content. Precisely determining
which galaxies belong to a cluster is crucial. In this paper, we present the
COlor Probabilistic Assignment of Clusters And BAyesiaN Analysis (Copacabana)
algorithm. Copacabana computes membership probabilities for all galaxies
within an aperture centred on the cluster using photometric redshifts,
colours, and projected radial probability density functions. We use
simulations to validate Copacabana and we show that it achieves up to 89%
membership accuracy with a mild dependency on photometric redshift
uncertainties and choice of aperture size. We find that the precision of the
photometric redshifts has the largest impact on the determination of the
membership probabilities followed by the choice of the cluster aperture size.
We also quantify how much these uncertainties in the membership probabilities
affect the stellar mass–cluster mass scaling relation, a relation that
directly impacts cosmology. Using the sum of the stellar masses weighted by
membership probabilities ($\rm\mu_{\star}$) as the observable, we find that
Copacabana can reach an accuracy of 0.06 dex in the measurement of the scaling
relation. These results indicate the potential of Copacabana and
$\rm\mu_{\star}$ to be used in cosmological analyses of optically selected
clusters in the future.

###### keywords:

galaxies: clusters: general, methods: data analysis

††pubyear: 2024††pagerange: Copacabana: A Probabilistic Membership Assignment
Method for Galaxy Clusters–A

##  1 Introduction

Galaxy clusters have long been considered a promising astrophysical probes of
dark energy (Albrecht et al., 2006; Allen et al., 2011; Dodelson et al., 2016)
as their abundance as a function of redshift and mass is sensitive to the
growth rate of structures in the Universe. Clusters are complementary to
geometry based probes such as type-Ia supernovae and can be used to test
different dark energy models (Huterer, 2023). The challenge in realizing this
promise is to obtain a well-understood sample of clusters with unbiased mass
measurements of a few per cent precision across a wide range of redshifts
($z\sim 0-1$) and masses ($M\sim 10^{13}-10^{15}\mathrm{M_{\odot}}$). In
galaxy clusters, more than 80% of the mass is in the form of dark matter while
5%–15% is diffuse hot gas, and only 1% – 5% is in galaxies (Gonzalez et al.,
2007, 2013; Laganá et al., 2013; Song et al., 2017; Pratt et al., 2019;
Umetsu, 2020). Apart from a few extremely massive clusters for which direct
total mass measurements via gravitational lensing are possible, we rely on
indirect scaling relations to infer cluster masses. For instance, the hot
intracluster gas has two main observational signatures: a thermal
bremsstrahlung (X-ray) emission and the Sunyaev-Zeldovich (SZ) effect.
Although these signals correlate strongly with cluster mass, they are reliably
detectable only for clusters at the high mass end where the gas temperature
and density are highest (e.g. Sarazin, 1988; Bleem et al., 2020; Klein et al.,
2022).

Large optical imaging surveys such as the Sloan Digital Sky Survey (SDSS; York
et al., 2000) and the Dark Energy Survey (DES; The Dark Energy Survey
Collaboration, 2005) have produced samples of tens of thousands of clusters
with masses $\lesssim 10^{14}\mathrm{M_{\odot}}$ (Rykoff et al., 2014; Rykoff
et al., 2016). For these low-mass galaxy clusters, the galaxy content is
crucial to unlocking their potential for cosmology (Wu et al., 2021).
Establishing observable quantities that correlate with cluster masses is a
challenge. One such quantity is richness ($\lambda$), defined as the
probability-weighted sum of red-sequence galaxies identified and selected by
cluster finding algorithms such as redMaPPer (Rykoff et al., 2016). Richness
is an empirical mass proxy optimized to find clusters and is correlated with
the cluster mass (Rozo et al., 2009; Rykoff et al., 2012). This quantity
relies on a linear colour-magnitude relation known as the red sequence.

The formation and evolution of the red-sequence are not well understood
(Butcher & Oemler, 1984; Andreon et al., 2006; De Lucia et al., 2007; Cooper
et al., 2007; Puddu et al., 2021). This lack of understanding poses a
challenge to optical and NIR wavelength cluster cosmology programs that use
red-galaxy counts as the mass-proxy as systematic uncertainties associated
with the mass-proxy scaling relations dominate the error budget (Wu et al.,
2022). One promising avenue to address this challenge is the development of a
mass-proxy that includes all of the galaxy content of the clusters. Such a
mass-proxy has a stronger theoretical foundation and, thus, can be simulated
and studied more easily than the red sequence (e.g.: Anbajagane et al., 2020).

Mass measurements can be biased at the low richness end if selection effects
are not treated properly. Significant effort has been made to understand the
systematics uncertainties associated with richness. Projection effects
(Costanzi et al., 2019; Myles et al., 2021) and optical selection bias are the
leading terms contributing to the cluster mass systematic uncertainties
(Sunayama et al., 2020; Wu et al., 2022).

Recently, alternative mass proxies have been studied, in particular, those
derived from the stellar mass (e.g. Andreon, 2012; Pereira et al., 2018;
Bradshaw et al., 2020) and intracluster light (e.g. Huang et al., 2022;
Golden-Marx et al., 2023). In Pereira et al. (2018); Palmese et al. (2020), we
designed a novel mass proxy around the stellar mass content of the galaxies in
a cluster. Defined as the weighted sum of the stellar masses, $\rm\mu_{\star}$
uses a more complete representation of the full population of the cluster
galaxies than red-galaxy count methods. In Palmese et al. (2020), we used an
X-ray sample of clusters to compare $\rm\mu_{\star}$–mass observable and the
$\lambda$ method, finding that both present a similar scatter in total mass
(see their Figure 6). In Pereira et al. (2020), we performed a detailed weak-
lensing mass calibration of $\rm\mu_{\star}$ and we found that the precision
of the mass–$\rm\mu_{\star}$–$z$ scaling relation was comparable to the one
obtained by McClintock et al. (2019) for a mass–richness–$z$ relation. Those
results indicate that $\rm\mu_{\star}$ has the potential to become a
competitive mass proxy for cluster cosmology.

In this study, we introduce a new methodology called COlor Probabilistic
Assignment of Clusters And BAyesiaN Analysis (Copacabana). This method assigns
probabilities for all galaxies in the cluster region, independent of the
cluster finder selection. Copacabana continues to improve the methodology of
previous $\rm\mu_{\star}$-based papers (Pereira et al., 2018; Palmese et al.,
2020).

Copacabana’s membership assignment enables value-added information for cluster
finders to be produced, even for those that do not rely on galaxy catalogues,
such as SZ and X-ray. Moreover, Copacabana can be used to study the evolution
and properties of galaxy clusters, such as their mass content and galaxy
population. This algorithm is particularly useful for X-ray and SZ-selected
samples, as their stellar-mass function and the baryon content of the
Universe, can then be analyzed (Gonzalez et al., 2013; Leauthaud et al., 2012;
Kravtsov et al., 2018).

In this work, we apply Copacabana to improve and validate the estimates of
$\rm\mu_{\star}$. We use the Buzzard DES Year 3 (Y3) simulations (DeRose et
al., 2022) to validate our algorithm by quantifying the probability’s impact
on the scaling relation. In addition, we study the impact of the uncertainties
in photometric redshifts of three large photometric surveys: SDSS, DES, and
the Legacy Survey of Space and Time (LSST; LSST Science Collaboration et al.,
2009).

This paper is organized as follows: Section 2 details the membership
assignment methodology based on the photometric redshift, colour and projected
radial probability distributions; Section 3 gives an overview of the simulated
dataset employed in our study; the validation of our algorithm is shown in
Section 4, and conclusions are presented in Section 5. Throughout this
manuscript, we use logarithm in base ten ($\log$) and we adopt the
cosmological parameter values: $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, and
$h=0.7$.

##  2 Formalism

In this section, we outline the method used in the Copacabana
algorithm111https://github.com/estevesjh/ccopa for assigning membership
probabilities to galaxies for a given cluster field. The main motivation is to
produce stellar mass estimations for cluster galaxies that only have
photometric information. The method presented here improves and extends the
algorithm used in Palmese et al. (2016); Pereira et al. (2018); Pereira et al.
(2020). This work has been inspired by previous papers by George et al.
(2011); Rykoff et al. (2014); Castignani & Benoist (2016). Copacabana has two
main differences relative to those papers: colour distribution and
optimisation of the cluster aperture. The red and blue galaxy populations are
modeled simultaneously. For the cluster aperture, we estimate
$R_{\mathrm{200c}}$ which is defined as the radius containing 200 times the
critical the density of the universe (at the cluster redshift).
$R_{\mathrm{200c}}$ is inferred from the galaxy distribution around each
cluster.

###  2.1 Membership Probabilities

For a given cluster photometric field, the galaxies present belong to only two
classes: gravitationally bound systems and field galaxies, i.e. in the
background and foreground galaxies. For the bound systems, we define member
galaxies as the population inside the $R_{\mathrm{200c}}$ defined by host halo
mass. We make this distinction given that there are projected correlated
structures, filaments, infalling groups, and galaxies. In simulations we know
which are the correlated galaxies, in data we do not.

We adopt a Bayesian inference approach to estimate membership probabilities.
Within this framework, the probability of a galaxy being a member of the
cluster is in general:

${\rm P(member|data)=\frac{P(data|member)P(member)}{P(data)}}\;,$ (1)

where “data” represents the galaxy input variables: cluster-centric distance,
photo-z and colour. The term $\rm P(data|member)$ is our likelihood
distribution, which is modeled as the product of the distributions of our
input variables (described in detail in Sec. 2.1.1). The prior $\rm P(member)$
is defined as the ratio of the number of member galaxies and the total number
of galaxies, i.e., $\rm n_{C}/(n_{C}+n_{F})$, where $\rm C$ and $\rm F$
indicate cluster and field, respectively.

The denominator, $\rm P(data)$, is the probability of the union $\rm(C\cup F)$
of the two groups:

$\rm P(data)=P(data|member)P(member)+P(data|field)P(field),$ (2)

where $\rm p(data|field)$ is the observed field distribution and $\rm
P(field)$ is $\rm 1-P(member)$.

####  2.1.1 Cluster Likelihood

The likelihood depends on the joint radial, photometric redshift and colour
distributions of the cluster galaxies:

$\rm P(data|member)=P(R,z_{p},c|member)\;,$ (3)

where the term on the right side is given by the product of each variable
(radius R, photometric redshift zp and colour c):

$\rm P(R,z_{p},c|member)=P(R|member)P(z_{p}|member)P(c|z_{p},member)$

We assume that these variables are independent. This assumption is an
approximation. A potential improvement, not explored in this paper, would be
to model the joint color and radial probabilities since galaxies are redder at
the cluster centre. Combining Equation 1 and Equation 2 we find:

$\rm P(member|R,z_{p},c)=\frac{P(R,z_{p},c|member)\times P(member)}{Q},$ (4)

where Q comes from the law of total probability:

$\mathrm{Q=P(member)P(R,z_{p},c|member)+P(field)P(R,z_{p},c|field)}.$

In principle, the full membership probability with the three variables has
more potential constraining power than probabilities using fewer variables.
For some cases, though, it might be useful to separate the impact of the
colour and photo-z variables. For example, we could study the impact of
photo-z outliers in cluster galaxies or for blue BCGs. For this reason, we
also compute the probability for each model variable:

$\rm P(member|R)=\frac{P(R|member)\times
P(member)}{P(member)P(R|member)+P(R|field)P(field)}\;,$ (5) $\rm
P(member|z_{p})=\frac{P(z_{p}|member)\times
P(member)}{P(member)P(z_{p}|member)+P(field)P(z_{p}|field)}\;,$ (6) $\rm
P(member|c)=\frac{P(c|member)\times
P(member)}{P(member)P(c|member)+P(field)P(c|field)}\;.$ (7)

In this formalism, the probabilities allow flexibility for the user to drop
the selection in a given variable if needed. As an example, in the case of
group galaxies, the assumption of a Navarro-Frenk-White (NFW) profile might
not be applicable, and we might, therefore, remove the probability related to
the radius R, so removing the radial filter. In the next section, we present
the definition of each probability.

####  2.1.2 Radial Filter

We assume that the cluster galaxy radial distribution is a projected NFW
profile (Wright & Brainerd, 2000), and a constant radial distribution for the
background. In our case, the NFW profile density has the form:

$\Sigma(R)=\begin{cases}\frac{2\rho_{s}R_{s}}{r^{2}-1}\bigg{[}1-\frac{2}{\sqrt{r^{2}-1}}\arctan{\sqrt{\frac{r-1}{r+1}}}\bigg{]}&r>1\\\
\frac{2\rho_{s}R_{s}}{3}&r=1\\\
\frac{2\rho_{s}R_{s}}{r^{2}-1}\bigg{[}1-\frac{2}{\sqrt{1-r^{2}}}\operatorname{arctanh}{\sqrt{\frac{1-r}{r+1}}}\bigg{]}&r>1\;,\\\
\end{cases}$ (8)

where $r=R/R_{s}$ is the dimensionless radial distance, $\rho_{s}$ is the
density scale parameter, and $R_{s}=R_{200c}/c_{200}$, where $c_{200}$ is the
concentration parameter, is a characteristic radius. To convert the surface
mass density profile to a radial probability density function (PDF), we
compute the normalization factor

$\rm Norm={\int_{0}^{R_{200}}{2\pi
R^{\prime}\Sigma(R^{\prime})}dR^{\prime}}\;$ (9)

such that

$P(R|\rm{member})=\Sigma(R,R_{200},c_{200})/\rm{Norm}\;.$ (10)

The NFW has two free parameters, the radius $R_{200}$ and the concentration
$c_{200}$. We infer $R_{200}$ using a halo occupation distribution (HOD) model
(see Sec. 2.4) and we set $c_{200}=3.59$ as this was shown to be a good fit
for halos in this mass range selected in the DES Science Verification dataset
(Hennig et al., 2017). A common assumption is a surface constant background
density $\Sigma_{\rm field}(R)=n_{bkg}$ (Rykoff et al., 2016). As a result,
the field radial density probability, $\rm P(R|field)$, is a constant value
determined by the normalization $\rm\int_{0}^{R_{max}}2\pi
R^{\prime}\Sigma_{\rm field}dR^{\prime}=1$.

####  2.1.3 Photometric Redshift Distribution

If the galaxy is assumed to be in the cluster we can set its probability
density function to be at the cluster redshift. In the scenario where all
member’s redshifts are known, the cluster galaxy redshift distribution can be
described by a normal distribution, with mean $z_{\rm cls}$ and standard
deviation $\sigma_{\rm{cls}}$. In the context of photometric redshifts for
clusters, the cluster redshift uncertainty $\delta z$ is much larger than
$\sigma_{\rm{cls}}$. For instance, the typical redMaPPer cluster error is
$\delta z=0.01(1+z_{\rm cls})$ (Rykoff et al., 2016), as opposed to a
spectroscopically derived $\sigma_{\rm{cls}}\approx 0.001$. In this regime,
the galaxy photo-z distribution can be described analogously by:

$P\left(z_{p}|\text{member}\right)=\frac{1}{\sqrt{2\pi\delta
z^{2}}}\mathrm{e}^{-\left(\mathrm{z}_{\mathrm{cls}}-\mathrm{z}_{\mathrm{p}}\right)^{2}/2\delta
z^{2}}$ (11)

where $z_{p}$ is the galaxy photo-z. Despite the simplicity of the model, it
is a robust estimator, as shown by Castignani & Benoist (2016). For a given
photo-z sample, the galaxy population photo-z can be biased relative to the
cluster redshift. In such a case an offset to the cluster redshift can be
applied (Aguena et al., 2021).

If the galaxy is assumed to be in the field, we need the field photometric
redshift distribution. We use Gaussian Kernel Density Estimation (KDE) to
estimate the photometric redshift density distribution in a ring around the
cluster centre:

$P(z_{p}|\mathrm{field})={\mathrm{KDE}}(z_{p}|R,h).$ (12)

The KDE has a bandwidth $h$ (i.e. the width of the Gaussian kernel) as a free
parameter. The choice of a proper bandwidth depends on the shape of the
underlying distribution and the number of objects available to construct the
estimator. Assuming a normal distribution, we take the optimal bandwidth as:

$h=\left(\frac{4\hat{\sigma}}{3n}\right)^{1/5}$ (13)

where $\hat{\sigma}$ is the sample standard deviation and $n$ denotes the
number of objects. This “Scott’s rule” bandwidth is in common usage in
statistics (Scott, 1992).

####  2.1.4 Color Distribution

Generalizing from Pereira et al. (2020), we add color probabilities by using a
color distribution subtraction method. The cluster color distribution is
computed by subtracting the background color distribution from the total color
distribution using a KDE:

$\displaystyle N_{\rm cls}P(c|\rm{member})={}$ $\displaystyle N_{\rm
total}{\rm KDE}(c|h_{\rm eff},\rm{total})$ (14) $\displaystyle-N_{\rm
bkg}^{\prime}{\rm KDE}(c|h,\rm{field}).$

The KDEs are weighted by the photo-z probability weights and $h_{\rm eff}$ is
a Scott rule bandwidth divided by 10, a factor that is arbitrarily chosen to
avoid over-smoothing. We only take the excess in the subtraction, so there are
no negative values. For $h_{\rm eff}$ the number of objects is computed in
terms of an effective number

$n_{\rm eff}=\frac{(\sum{\rm weights})^{2}}{\sum{\rm weights^{2}}}.$

After the subtraction, we normalize $p(c|z,R)$ to unity. We choose one color
filter at a time to be our color model:

$c(z)=\begin{cases}(g-i)\quad\text{for }z\leq 0.35\\\ (r-z)\quad\text{for
}z>0.35.\end{cases}$ (15)

The color filter changes at $z=0.35$ due to the $4000$ Å break exiting the red
edge of the $g$ band. In comparison with Pereira et al. (2020), the addition
of color probabilities in general improved the performance.

###  2.2 Number Densities

Computing probability weighted number density is the first step in our
algorithm, and is used in all the steps were galaxies are counted, i.e. the
background subtraction, the estimation of the radius $R_{200}^{\rm HOD}$ and
the color model. The first step is to count galaxies around the cluster
center, within a cylinder of radius $R_{200}$ and height $2\sigma_{z}(1+z_{\rm
cls})$, where $z_{\rm cls}$ is the redshift of the cluster. The galaxies are
counted using a probabilistic weight. The probability of a galaxy being at the
cluster redshift, $P_{z_{0}}$, is the integral of the galaxy’s photometric
redshift distribution ($\Pi(z)$) around the cluster redshift. The limits of
integration are chosen in a window of $2\sigma_{z}(1+z_{\rm cls})$ around the
mean value:

$P_{z_{0}}=\int_{z_{-}}^{z_{+}}\Pi(z)\;dz,\qquad z_{\pm}=z_{\rm cls}\pm
2\times\sigma_{z,0}(1+z_{\rm cls})\;,$ (16)

The limits of the integration depend on the photo-z precision for a given
redshift, $\sigma_{z,0}(z)$. This counting method selects galaxies near the
cluster redshift while avoiding a sharp, arbitrary cutoff in redshift space.

This probability weight is the first step in our algorithm and is used in all
the steps where galaxies are counted, i.e. the background subtraction, the
estimation of the radius $R_{200}^{\rm HOD}$ and the color model.

###  2.3 Background Subtraction

Background subtraction is an essential step for computing membership
probabilities. There are two methods traditionally used: global and local
subtraction. The global background subtraction method assumes that the
background density only depends on the redshift (e.g. Rykoff et al., 2014).
However, this assumption must be invalid as clusters are nodes of the cosmic
web. Consequently, galaxies in the line of sight are more likely to be
assigned as members.

In our work, we choose to use a local background subtraction method that
probes the surroundings of each cluster region. Specifically, we compute the
galaxy density in an annulus centered on the cluster with inner and outer
radii of $4$ Mpc and $6$ Mpc, respectively. The inner radii are always larger
than $R_{\mathrm{200c}}$, even for the most massive clusters for which
$R_{\mathrm{200c}}$ is approximately $3$ Mpc. Although scaling the radii with
the cluster radius would be an optimal choice, we prefer to use fixed values
since our $R_{\mathrm{200c}}$ estimation depends on the background density.

###  2.4 $R_{\mathrm{200c}}$ estimator: HOD Model

Clusters don’t have obvious edges, and various investigations define cluster
apertures differently. For instance, redMaPPer (Rykoff et al., 2014) assumes
an aperture that scales with richness $\lambda$, and AMICO (Bellagamba et al.,
2018) assumes a fixed aperture corresponding to a cluster with $M_{\rm
200,c}=10^{13.5}\;M_{\odot}$. In this work, we introduce a new cluster
aperture estimator based on a Halo Occupation Distribution (HOD) model, which
is independent of our mass proxy.

Our aperture estimator uses the galaxy number density profile of a HOD model.
A given HOD model provides the number of halo galaxies as a function of mass
which allows us to convert the number density profile to a mass density
profile. Assuming spherical symmetry, we can calculate the mass density,
$\rho$. We can make a rough estimation of $R_{200}$ by interpolating the mass
density profile as a function of radii. By definition, where the mass density
profile is 200 times the critical density, we have our aperture estimation,
$R_{\rm 200c}^{\rm HOD}$. For this work, we adopt the HOD model of Tinker et
al. (2012). The model consists of a relation between the number of central
($N_{\rm cen}$) and satellite galaxies ($N_{\rm sat}$) inside a halo of given
mass ($M_{\rm 200,c}$) and below a given luminosity threshold. The occupation
function for central galaxies takes the form:

$\langle N_{\rm
cen}\rangle_{M}=\frac{1}{2}\Bigg{[}1+\operatorname{erf}{\Bigg{(}\frac{\log
M-\log M_{\rm min}}{\sigma_{\log M}}\Bigg{)}}\Bigg{]},$ (17)

where $M$ is the halo mass, $M_{\rm min}$ represents the halo mass at which
the probability of containing a central galaxy is $50$ per cent, and
$\sigma_{\log M}$ accounts for the scatter in halo mass at a fixed luminosity
of the galaxy population. The occupation function for satellites is given by a
power law:

$\langle N_{\rm sat}\rangle_{M}=\langle N_{\rm
cen}\rangle_{M}\times\Bigg{(}\frac{M}{M_{\rm sat}}\Bigg{)}^{\alpha_{\rm
sat}}\exp{\Bigg{(}-\frac{M_{\rm cut}}{M}\Bigg{)}},$ (18)

where $\alpha_{\rm sat}$ is the slope at high halo masses, with an exponential
cutoff at halo masses below $M_{\rm cut}$, and $M_{\rm sat}$ is the
characteristic halo mass for satellites. Combining the central and satellite
occupation functions produces a total occupation function of the form

$\langle N_{\rm tot}\rangle_{M}=\langle N_{\rm
cen}\rangle_{M}\times\Bigg{[}1+\Bigg{(}\frac{M}{M_{\rm
sat}}\Bigg{)}^{\alpha_{\rm sat}}\exp{\Bigg{(}-\frac{M_{\rm
cut}}{M}\Bigg{)}}\Bigg{]}.$ (19)

There are five free parameters this model. The best-fit values for these
parameters (given in Table 4 of Tinker et al. (2012)) were derived from the
SDSS dataset, using the maxBCG cluster sample (Koester et al., 2007), with an
absolute magnitude cut $M_{r}^{0.1}\leq-19.5$. We use these values as a
reference since they are close to the $0.2L_{\star}$ cut applied here (Rykoff
et al., 2012).

We found that this approach results in $R_{\mathrm{200c}}$ estimates that are
biased low. To account for this bias, we introduce a calibration factor which
can be computed as the mean ratio of our predictions and the actual
$R_{\mathrm{200c}}$ values in the simulations:

$\rm\eta_{HOD}=\frac{R_{HOD,200c}}{R_{200c,true}}\qquad\text{where},\qquad
R^{3}_{HOD,200c}=200\frac{M(N_{200c})}{4\pi\rho_{c}/3}\;.$ (20)

This calibration factor is independent of redshift, as we wish to use only one
factor for the entire population.

###  2.5 Stellar Mass Estimation: BMA

The code we call BMA222https://github.com/apalmese/BMAStellarMasses is a
Bayesian model averaging code (see e.g. Taylor et al., 2011) applied to the
output of a stellar population synthesis code, and developed into a pipeline
(Palmese et al., 2020). We use the stellar population synthesis code FSPS
(Conroy et al., 2009; Conroy & Gunn, 2010) to evaluate a 5-dimensional space
of quantities, resulting in 24 models. We choose the models evaluated at the
cluster redshift for a given galaxy with apparent magnitude, colours, and
photo-z. Then the likelihood of each model given the galaxy magnitudes,
colours, and errors is computed. The properties of interest from the models,
e.g. stellar mass, are then computed as the likelihood-weighted sum over all
models, a Bayesian model average. The code was validated on the Millennium
simulations (Springel et al., 2005) and on the COSMOS dataset (Laigle et al.,
2016). See Palmese et al. (2020) for a full description of the BMA
methodology.

###  2.6 $\mu$⋆ Estimator

A key motivation for the development of Copacabana is to improve measurements
of the mass proxy $\rm\mu_{\star}$. To accomplish that, we focus on the
membership probabilities. $\rm\mu_{\star}$ is defined in a probabilistic
manner:

$\mu_{\star}=\sum_{i}P_{{\rm mem},i}M_{\star,i}\;\;\text{for}\;\;R\leq R_{\rm
aper}\;,$ (21)

where $M_{\star}$ is the galaxy’s stellar mass, $P_{\rm mem}$ is the
membership probability (Equation 4). $R_{\rm aper}$ is the cluster radius
aperture which is $R_{\mathrm{200c}}$ if not defined otherwise. A galaxy’s
stellar mass is estimated from photometric data assuming it is at the cluster
redshift via BMA.

##  3 Validation Setup

###  3.1 Data

To validate Copacabana, we use the Buzzard v2.0 simulations (DeRose et al.,
2019), meant to correspond with the DES Y3 area. The dataset consists of
synthetic dark matter simulations with galaxy information added by the AddGals
algorithm (Wechsler et al., 2022). This procedure places galaxies onto the
dark-matter-only simulation, weighted by local dark matter density, matching
the observed luminosity function and luminosity-dependent two-point
correlation function.

For context, we briefly describe how the galaxies are pasted onto the dark
matter particles. First, AddGals creates a catalogue of galaxies based on the
luminosity function, $\phi(M_{r})$, performing subhalo abundance matching
between a small high-resolution N-body simulation and the observed SDSS
luminosity function in the r-band. The algorithm calibrates a relation (not a
HOD) for the central and the non-central galaxies, which is then evolved using
a functional form. The model of the central galaxies is a log-normal
distribution at a fixed halo mass and redshift, $P(M_{r,{\rm cen}}|M_{\rm
vir},z)$, and the model of the non-centrals is based on the local density.
These relations are then used to assign galaxies to resolved haloes or dark
matter particles in a large light-cone simulation with a lower resolution.
This modelling scheme was chosen such that it predicts the clustering in SDSS
to high precision (e.g. Conroy et al., 2006; Reddick et al., 2013; Lehmann et
al., 2017). Second, colours are assigned using a spectral energy distribution
(SED), chosen such that the simulation matches the SED distribution (at fixed
luminosity and galaxy density) measured in the SDSS data.

Figure 1: Number of galaxies inside $R_{\rm 200,c}$ ($N_{200}$) as a function
of the halo mass $M_{\rm 200,c}$, for three absolute magnitude cuts in the
r-band. The Buzzard galaxy distribution follows the Tinker et al. (2012) model
closely (red dotted line) when both apply the same magnitude selection of
$M_{r}^{0.1}\leq-19.5\,\text{ mag}$ (dark blue points).

The galaxies in Buzzard, unlike many cosmological simulations, are not placed
using a HOD prescription. Nonetheless, the halo occupation distribution on
Buzzard follows the Tinker et al. (2012) model closely for an absolute
magnitude selection of $M_{r}^{0.1}\leq-19.5$, as we can see in Figure 1. In
addition, Zacharegkas et al. (2022), using a Buzzard redMaGic galaxy
selection, found acceptable HOD model fits (Equation 17, Equation 18) to the
galaxy distribution.

####  3.1.1 Sample Selection

In order to accurately assess the performance of our code across the halo
redshift and mass ranges, we select with bin-dependent uniform probability
halos in bins of redshift and logarithmic mass, i.e. [$\log M_{\rm 200,c},z$].
The aim is to ensure we have the same number of halos for each bin, which we
do not quite achieve in high mass bins due to a lack of clusters. This
approach prevents our assessment from being biased by the low end of the halo
distribution.

The data chosen by our selection is presented in Figure 2, where 2,200 halos
are plotted in $(z,\log M_{\rm 200,c})$ space with histograms on the $x$– and
$y$–axis. The limits of the sample are $z\in\left[0.1,0.65\right]$ and $(\log
M_{\rm 200,c})>13.5\;M_{\odot}/h$. The choice of the redshift range follows
the DES cluster cosmology analysis (Abbott et al., 2020). The halo mass
threshold is similar to that adopted by other cluster finder algorithms (e.g.:
Castignani & Benoist, 2016; Bellagamba et al., 2019). Overall, there is a
uniform selection, except for the highest mass bins, where there are not
enough systems.

Figure 2: Uniform selection of 2,200 Buzzard v2.0 halos on a halo
mass–redshift grid. The upper and right panels are the redshift and halo mass
distributions, respectively.

####  3.1.2 Simulated Photo-z

To validate the Copacabana algorithm with respect to photo-z, we add offsets
to the simulated galaxy redshifts. In detail, we draw a random offset
following:

$z_{p}=\mathcal{N}\left(z,\sigma_{z,0}(1+z)\right)\;,$ (22)

where $\sigma_{z,0}$ is the photo-z precision that correspondents to a typical
photo-z error.

We simulate three different levels of uncertainty:
$\sigma_{z,0}=0.01,0.03,\text{ and }0.05$ as an ideal, a realistic, and a
pessimistic case, respectively. These choices mimic three different surveys:
LSST (LSST Science Collaboration et al., 2009), DES (Gschwend et al., 2018;
Aguena et al., 2021) and SDSS (Carliles et al., 2010), respectively.

In the context of clusters, the main differences between simulated Gaussian
photo-z’s with real data photo-z’s are the bias and the presence of outliers.
For instance, Aguena et al. (2021) using the WaZP cluster catalogue studied
the differences of redMaPPer cluster redshift with the ones derived from the
DNF photo-z algorithm in the DES Y1 (Gschwend et al., 2018). They quantified a
redshift bias that is less than $0.003\times(1+z)$. For future applications of
Copacabana on data, a description of the bias between the photo-z sample
employed and the cluster redshift must be taken into account as a bias on the
cluster photo-z distribution (Equation 11).

###  3.2 Validation Metrics

####  3.2.1 Assessing $\mu$⋆ precision

To validate our $\rm\mu_{\star}$ estimation, we compare it with the simulation
cluster member stellar masses. For this purpose, we define
$\rm\mu_{\star,true}$,

$\mu_{\star,\rm true}=\sum_{i\in{\rm
members}}M_{\star,i}\;\;\text{with}\;\;R\leq R_{\rm aper}\;,$ (23)

as the sum of the cluster members’ stellar masses, where the “true” members
are defined as the galaxies inside the three-dimensional $R_{\mathrm{200c}}$
distance from the cluster centre. In other words, we don’t consider the line-
of-sight infall galaxies, nor the gravitational status of the galaxy.

Next we can define the ratio $x\equiv\mu_{\star}/\mu_{\star,\rm true}$. Since
our data has a non-Gaussian tail at very low richnesses, $N_{\rm gal}<10$, a
robust metric is adopted, the scaled median absolute deviation (MAD):

$\sigma_{\rm
MAD}(\log(x))=1.48\times\text{Median}\left(\left|\log(x)-\text{Median}\left({\log}(x)\right)\right|\right)\;.$
(24)

Note that if $\log(x)$ follows a normal distribution $\sigma_{\rm MAD}$ is
equal to the standard deviation.

It is important to stress that our assessment is primarily on our estimator
due to membership probabilities and we do not evaluate uncertainty due to
stellar mass estimates.

####  3.2.2 Assessing the accuracy of $R_{200c}$

To validate our estimates of $R_{\mathrm{200c}}$, we use the current value
from the simulation. $R_{200c,\rm true}$ was retrieved from the Buzzard truth
table. Analogously to $\rm\mu_{\star}$, we evaluate Equation 24 with x set to
the ratio of true versus measured $R_{\mathrm{200c}}$.

####  3.2.3 Completeness and Purity

The membership probabilities play a role in thresholds that distinguish the
classes in the framework of classifying members and non-members of a given
galaxy cluster. To assess the performance, we use metrics commonly used in
statistical classification problems, purity ($P$) and completeness ($C$).These
metrics rely on true positive (TP), false positive (FP), and false negative
(FN) predictions:

$P=\frac{TP}{TP+FP}\quad\text{and}\quad C=\frac{TP}{TP+FN};,$ (25)

TP represents correct positive predictions, while FP and FN refer to incorrect
positive and negative predictions, respectively. Purity indicates the
proportion of positives that are cluster members. Completeness measures the
fraction of true members that were successfully identified among all the
selected galaxies.

The overall accuracy of a classifier can be evaluated by:

$\text{accuracy}=\frac{TP+TN}{TP+FP+TN+FN}\,,$ (26)

where TN is the true negatives, i.e. correctly identified field galaxies.

####  3.2.4 $\mu$⋆-Cluster Mass Scaling Relation

For photometric surveys, one of the main requirements for tight constraints on
cosmological parameters is a mass proxy that predicts the cluster mass with
significant accuracy and is robust against systematic effects. Here, we assess
the possible impact of the membership probabilities on deriving cosmological
results using the relation between the weighted stellar mass and the total
cluster mass, which we will refer to as the $\rm\mu_{\star}$ \- cluster mass
scaling relation.

In simulations, the $\rm\mu_{\star}$ \- cluster mass scaling relation is
accessible since the halo mass is known. The probability of a given halo of
mass $M_{\rm 200,c}$ to have a $\mu_{\star}$ value is generally modeled by a
Log-Normal relation with mean:

$\left<\log(\mu_{\star})|M_{\rm 200,c}\right>=\alpha+\beta\log(M_{\rm
200,c}/M_{p})\;,$ (27)

with an associated intrinsic error $\sigma$. Here $\alpha$ is the intercept,
$\beta$ is the slope, $M_{p}=10^{15.5}M_{\odot}$ is the pivot mass. The
inference of the model parameters is made by employing a hierarchical Bayesian
algorithm (linmix; Kelly, 2007). The linmix algorithm allows us to include the
error on the $y$–dependent variable, in our case, $\rm\mu_{\star}$.

In general the scaling relation evolves with redshift. We model redshift
evolution by fitting the observable–mass relation in different redshift bins.
We follow the standard parametrization of the evolution as proportional to the
power of the scale factor $a\equiv 1/(1+z)$ or the dimensionless Hubble
parameter $E(z)\equiv H(z)/H_{0}$ (Evrard et al., 2014).

##  4 Results

In this section, we examine the performance of Copacabana using the Buzzard
simulation.

###  4.1 Uncertainty in $\mu$⋆ Estimations

We run Copacabana on the Buzzard v2.0 simulation using the three values of
$\sigma_{z,0}$ presented in Sec. 3.1.2. We employ the photometric stellar
masses computed by BMA at the cluster redshift (Sec. 2.5). The stellar mass
$\rm\mu_{\star}$ is computed within an aperture $R_{\mathrm{200c}}$ estimated
using the HOD model presented in Sec. 2.4. The membership probabilities are
expected to depend mainly on the photo-z uncertainty.

Figure 3 shows the estimated $\rm\mu_{\star}$ versus $\rm\mu_{\star,true}$
within $R_{\mathrm{200c}}$ for three values of $\sigma_{z,0}$. The
$\rm\mu_{\star}$ values for all halo mass regimes follow the
$\rm\mu_{\star,true}$ values closely. The $\rm\mu_{\star}$ errors are roughly
within $0.20\text{ dex }\sqrt{\mu_{\star,p}/\mu_{\star}}$, where
$\mu_{\rm\star,p}=10^{12.22}M_{\odot}$ is the sample mean. Uncertainty in the
position of the galaxies along the line-of-sight adds significant noise,
especially for the low-mass halos as they have fewer galaxies. As a result,
the quality of the photometric redshifts substantially impacts the
$\rm\mu_{\star}$ measurements and particularly does at the low-mass end.
Nevertheless, a sample with accurate and precise photometric redshifts can
accurately predict $\rm\mu_{\star}$ down to the lowest bin, $5\times
10^{11}M_{\odot}$.

Figure 3: Copacabana predicted $\rm\mu_{\star}$ as a function of
$\rm\mu_{\star,true}$ (the sum of the stellar masses of the true cluster
members within $R_{\mathrm{200c}}$) for different values of $\sigma_{z,0}$,
the uncertainty in the galaxy photometric redshifts. The binned points were
slightly shifted for better visualization. The estimated values closely follow
the one-to-one relation (grey dotted line). For comparison, we also show the
range (in light grey) corresponding to a Poissonian error of $0.20\text{ dex
}\sqrt{\mu_{\star,p}/\mu_{\star}}$ where $\mu_{\rm\star,p}$ is the sample mean
of $10^{12.22}M_{\odot}$. Note that our result is robust well below the
threshold of $10^{12}M_{\odot}$ used in our previous work (Pereira et al.,
2020).

####  4.1.1 Photometric Redshift Uncertainties And Cluster Apertures

In this section, we quantify how the uncertainty in $\rm\mu_{\star}$ depends
on $\sigma_{z,0}$ and the size of the apertures used to estimate
$\rm\mu_{\star}$. We show that the photo-z quality is the main systematic on
the $\rm\mu_{\star}$ error budget.

In Figure 4, we show the $\sigma_{\rm MAD}$ (defined in Equation 24) as a
function of the cluster aperture and photo-z uncertainty. The accuracy of
$\rm\mu_{\star}$ correlates linearly with the uncertainty in the photometric
redshifts. For instance, an improvement of a factor 5 (1.7) in $\sigma_{z,0}$
reduces the $\rm\mu_{\star}$ error to $0.15\text{ dex}$ ($0.09\text{ dex}$)
when compared to $\sigma_{z,0}=0.05$. This improvement in the photometric
redshifts has a big impact, especially for low-mass halos, see Figure 3. This
result follows from the fact that the uncertainty along the line of sight is
the major source of galaxy membership contamination. For example, a redshift
error of $0.01\times(1+z)$ translates into a physical length of $\sim 40$ Mpc
which is $\sim 10-50\times R_{\rm 200,c}$.

Figure 4: Copacabana $\sigma_{\rm MAD}$ error of the $\log$
$\rm\mu_{\star}$$/$$\rm\mu_{\star,true}$ as a function of the cluster aperture
radius for three different photo-z precisions. The $\sigma_{\rm MAD}$ error is
always below $0.16$ dex, and it varies significantly with the photo-z quality
(different markers). In particular, for an LSST-like photo-z (pink markers),
$\sigma_{\rm MAD}$ is $0.09$ dex. The precision can be improved considerably
by defining a smaller cluster aperture since the density is higher in the
cluster centre. The error bars were estimated by bootstrapping the cluster
sample a thousand times.

A second source of error is the cluster aperture, which has a non-negligible
effect on the uncertainty in $\rm\mu_{\star}$. As shown in Figure 4, smaller
cluster apertures decrease the uncertainties in $\rm\mu_{\star}$. This effect
is driven by the higher galaxy density in the core. Lopes & Ribeiro (2020)
showed that the local density is a very good indicator of membership galaxies
in their well-characterized cluster sample.

The choice of the cluster aperture is important for studying cluster
cosmology. Ideally, we would expect that $\rm\mu_{\star}$ has the highest
correlation with the halo mass when computed at $R_{\mathrm{200c}}$ (Rykoff et
al., 2012; Bradshaw et al., 2020; Huang et al., 2020). In Sec. 4.4 we discuss
further the choice of the cluster aperture in terms of optimizing the scatter
of the $\rm\mu_{\star}$–$M_{\rm 200,c}$ relation.

####  4.1.2 Stellar Mass Estimation

Thus far, the $\rm\mu_{\star}$ uncertainty was computed without taking into
account the errors on stellar mass. The BMA stellar-mass error is around 0.2
dex (Palmese et al., 2020), validated using the stellar masses computed with
multi-band data in 16 filters from UV to infrared of the COSMOS deep field
(Laigle et al., 2016). The BMA stellar-mass errors are comparable to the
uncertainties induced by the photo-z errors, see Figure 4. Therefore, they
have a significant impact. To quantify this impact, we add random noise
normally distributed to the estimated stellar masses. We assume the typical
BMA error of 0.20 dex. This assumption should set an upper bound on the
uncertainty in $\rm\mu_{\star}$. Because $\rm\mu_{\star}$ is dominated by high
stellar mass (bright) galaxies that have lower mass uncertainties.

The additional scatter on $\rm\mu_{\star}$ is $0.07\pm 0.01\,\rm{dex}$ for the
three photo-z samples. This result is equivalent to adding the stellar mass
error in quadrature, $\sigma_{\rm MAD}^{2}+\sigma_{\rm BMA}^{2}$. The
implication of this additional error for the DES-like photo-z case is that the
$\rm\mu_{\star}$ error is at the same level as an SDSS-like photo-z. In future
work, it would be important to reduce the uncertainty on stellar masses to
reduce the uncertainty on $\rm\mu_{\star}$.

###  4.2 Precision of $R_{\mathrm{200c}}$ Estimations

The new $R_{\mathrm{200c}}$ estimator is based on a HOD model, which is
sensitive to the relation between the number of galaxies and the halo mass. To
perform our measurements, we applied a calibration factor $\eta_{HOD}$
(Equation 20). The derived correction factor is $0.63\pm 0.11$ for the same
HOD luminosity cut. We use all Buzzard halos with $M_{\rm 200,c}>5\times
10^{13}M_{\odot}$ for this computation. Note that only this correction factor
is set to calibrate the HOD relation.

Our estimator predicts the $R_{\mathrm{200c}}$ for halo masses probed in this
study as shown in Figure 5 with a scatter around $30\%$ after calibration.
Although the scatter is large, our estimated values correlate with the $R_{\rm
200c,true}$ for low and high-mass clusters, which makes them a good probe of
the cluster size. Interestingly, the photo-z does not have a significant
impact on the predicted values. Likewise, the accuracy of the estimated radius
is independent of cluster mass, unlike what we found for $\rm\mu_{\star}$.
These results indicate that the scatter observed in Figure 5 comes mainly from
another source, likely the intrinsic error of the HOD relation.

Figure 5: The Copacabana predicted $R_{\mathrm{200c}}$ based on a HOD model
versus the true radius ($R_{\rm 200c,true}$) for three values $\sigma_{z,0}$.

The significant uncertainties in our estimates of $R_{\mathrm{200c}}$ relate
to the uncertainty of the HOD model itself. For instance, the HOD model halo
mass scatter is $0.204\;\rm dex$. Also, optical mass proxies, in general, have
similar or higher intrinsic errors, e.g. the brightest cluster galaxy (BCG)
stellar-mass proxy intrinsic scatter is $0.20$ (Behroozi et al., 2010; Tinker
et al., 2012). In the future, if there are precise stellar mass measurements,
the stellar-to-halo mass can be incorporated into our methodology (Behroozi et
al., 2010).

###  4.3 Completeness and Purity of the Members List

In Figure 6, we plot the purity versus completeness for different values of
$\sigma_{z,0}$, the photometric redshift uncertainty. The figure is
constructed by varying the cluster membership probability threshold. The
optimal choice (coloured points) is the one that maximizes the product of both
quantities.

Figure 6: Purity (P) versus completeness (C) for different photo-z samples.
The purity of the membership selection is dependent on the accuracy of the
photometric redshifts. For instance, with an optimistic accuracy of
$\sigma_{z,0}=0.01$, the purity and completeness reach an optimal value of
$79\%$ and $94\%$ respectively. The DNF photo-z algorithm has an accuracy that
is similar to our realistic $\sigma_{z,0}=0.03$ scenario. However, it has
lower completeness due to outliers, which are not present in the artificial
Gaussian photo-z sample.

Overall, Copacabana performs well for different photo-z samples when compared
with other classifiers (George et al., 2011; Castignani & Benoist, 2016; Lopes
& Ribeiro, 2020). The product of purity (P) and completeness (C) is maximised
a $P=64\%$ and $C=93\%$ in an optimistic scenario ($\sigma_{z,0}=0.01$), and
values of $P=45\%$ and $C=80\%$ in the worst scenario ($\sigma_{z,0}=0.05$).
It is important to note that the completeness is not higher than $96\%$. The
$2\sigma$ photo-z threshold we use translates into $\sim 5\%$ of galaxies
having $P_{\rm mem}=0$. This threshold avoids outliers that might be present
in the color distributions.

The membership accuracy is higher for the smaller $0.5\times
R_{\mathrm{200c}}$ aperture (see the dashed lines in Figure 6) and an accuracy
of $89\%$ is achieved in the best scenario. For a given science case, for
instance, for studies of the red-sequence, Copacabana can provide an excellent
membership selection without relying on previous knowledge of the red-
sequence.

###  4.4 $\mu_{\star}$$-M_{200,c}$ Scaling Relation

An example of the $\mu_{\star}-M_{\rm 200,c}$ scaling relations is shown in
Figure 7, using $\rm\mu_{\star}$ and $\rm\mu_{\star,true}$ as the mass proxy
for clusters in the the lowest redshift bin, $0.47<z<0.56$, and using DES-like
photo-z accuracies, $\sigma_{z,0}=0.03$. The actual relation,
$\rm\mu_{\star,true}$$-M_{\rm 200,c}$ (grey line) is consistent (within
$2\sigma$) of with the purple line that Copacabana estimated. We note some
small differences between the two curves. For example, the intrinsic scatter
is larger, and the slope is shallower for the purple line. These two
differences hint at how membership probabilities bias our results of the
actual scaling relation. In the following section, we present and discuss the
impact of the quality of the photometric redshifts and size cluster aperture
on the fitted parameters.

Figure 7: Scaling relation, $\mu_{\star}$ versus $M_{200,c}$. In purple, the
DES-like photo-z, $\sigma_{z,0}=0.03$, and in grey, the true relation,
$\rm\mu_{\star,true}$. The fitted linear relation using a Bayesian regression
linear method (linmix) is shown in solid lines, and the $68\%$ confidence
level is shown in shaded bands.

####  4.4.1 Fitted Parameters

For cosmological parameter estimation, the scatter at a fixed $\rm\mu_{\star}$
is the important quantity describing the halo mass function. At first order,
the scatter at fixed mass-proxy ($\sigma_{\log M|\log\mu_{\star}}$) can be
written as (e.g. Evrard et al., 2014):

$\sigma_{\log M|\log\mu_{\star}}=\sigma_{\log{M_{\rm 200,c}}}/\beta\;.$ (28)

A steeper slope results in a lower $\sigma_{\log M|\log\mu_{\star}}$ just as
much as a smaller intrinsic scatter $\sigma$. For this reason, we focus on the
slope $\beta$ and the scatter $\sigma$ since they are the important parameters
for cluster cosmological analysis.

Figure 8: Scatter versus redshift for two stellar mass estimations
$\rm\mu_{\star,true}$ (True, grey points) and $\rm\mu_{\star}$(Copacabana,
coloured points). Each row displays the results for a given photo-z sample and
each column for a given cluster aperture. Down in columns, the scatter
increases as the uncertainty in the photometric redshifts increases, while the
scatter for the true relation is fixed. The coloured dashed lines show the
mean values of the scatter in each sub-plot for visual comparison with the
mean scatter of the true relation (grey dashed lines). The grey stars are the
scatter of the true relation, and therefore do not change between rows. Across
a row, the scatter in the true relation decreases.

The fitted scatter and slope values as a function of redshift are shown in
Figure 8 and Figure 9, respectively. The panel displays three different
photo-z precisions across rows and three different cluster apertures across
columns. There is an overall shift of the Copacabana from the true values,
indicated by the mean values (dashed lines). The gap, i.e. the additional
shift, increases with a poorer photo-z precision and a larger cluster
aperture.

Figure 9: As in Figure 8 but for the slope of the scaling relation. A larger
cluster aperture increases the contribution from background galaxies and
produces a shallower slope. As a result, the $0.5\times R_{\mathrm{200c}}$
aperture has the steepest slope, and it is the only case where the slope
derived by Copacabana agrees with the true one irrespective of the photo-z
sample used.

In particular, the fitted scatter does not change significantly with the
cluster aperture. In Figure 8, the mean values are about the same across the
columns for any photo-z. While the difference between dashed lines is larger,
the actual intrinsic scatter (grey stars) decreases with the cluster aperture,
which counterbalances the noise added by a larger aperture. The outcome,
stressed by the mean value, is that the intrinsic scatter does not depend
significantly on the cluster aperture. In contrast, intrinsic scatter did
depend on aperture when considering the $\rm\mu_{\star}$ precision in Sec.
4.1.1.

Regarding the slope, Figure 9 demonstrates that a decrease in the accuracy of
membership probability tends to result in a shallower slope. As discussed
earlier in Sec. 4.1.1, larger photo-z precision and a larger aperture size can
reduce the accuracy of the membership probability. This effect is most
noticeable in the bottom left panel of Figure 9, where the most significant
discrepancy with the actual value is seen due to the combination of low
photo-z precision and a large aperture. When we consider a specific survey
scenario (represented by a given row), the mean slope (colorful dashed line)
tends to be shallower when the cluster aperture is larger. This trend is
reduced with better photo-z accuracy.

The simultaneous change on the slope and the scatter imply that the scatter at
fixed $\rm\mu_{\star}$, $\sigma_{\rm logM|log\mu_{\star}}$, is affected by the
uncertainty on the membership probabilities. Using Equation 28, we can infer
$\sigma_{\rm logM|log\mu_{\star}}$ and quantify the impact of our methodology
on the scaling relation parameters. The Figure 10 shows these results
displayed similarly to that of Figure 8.

Figure 10: Scatter at fixed observable relation versus the redshift for two
independent variables $\rm\mu_{\star,true}$(True) and
$\rm\mu_{\star}$(Copacabana). Each row displays a result for a given photo-z
sample, and each column for a given cluster aperture. The photo-z precision
impacts the scatter of the recovered $\rm\mu_{\star}$$-M_{\rm 200,c}$ scaling
relation for any given cluster aperture. The scatter at fixed $\rm\mu_{\star}$
has a trend with the cluster aperture. The gap between the colourful and grey
dashed lines increases across the rows and columns. Note: the grey stars are
fixed between rows.

A smaller cluster aperture improves the $\rm\mu_{\star}$ measurement by
introducing less noise to the intrinsic scatter of the observable-mass
relation. However, for an optimistic photo-z sample such as LSST with
$\sigma_{z,0}=0.01$, using scatter computed at the $R_{\mathrm{200c}}$ is a
feasible option, as suggested by the observed trend in the mean scatter value.
Conversely, when dealing with an SDSS-like photo-z sample, where
$\sigma_{z,0}=0.05$, opting for a smaller cluster aperture may be the optimal
approach to reduce the noise introduced by uncertainty in the redshift.

The remarkable precision achieved by using all the cluster population, blue
and red galaxies, in photometric data demonstrates the power of the Copacabana
probabilities. In this work, we have obtained a scatter of $0.19$ dex for our
mass proxy $\rm\mu_{\star}$ in a DES-like photo-z sample. For comparison, the
typical scatter for $\lambda$ is $0.20$ dex, and for X-ray, luminosity, and
temperature are $0.30$ dex and $0.21$ dex, respectively. Additionally,
$\rm\mu_{\star}$ is an astrophysically motivated mass-proxy since it contains
the stellar information of the entire cluster population. Previous works show
that the stellar content has a low intrinsic scatter value according to
simulations (Anbajagane et al., 2020; Ho et al., 2023). Moreover, there are
hints that a mass proxy that only uses the BCG stellar mass has a comparable
scatter value to X-ray luminosity and temperature (Mulroy et al., 2019; Ho et
al., 2023).

Regarding the redshift evolution of the fitted parameters, we do not see a
smooth evolution with redshift as we would expect. Instead, there is a jump
around redshift $0.3$. This feature is also in the true underlying relation.
DeRose et al. (2019) showed that the switch on the simulation resolution at
$z=0.32$ impacts the matter power spectrum of the Buzzard data. This change in
resolution might explain the redshift discontinuity seen in this work.

####  4.4.2 $\rm\mu_{\star}$ as a low scatter mass-proxy

The total galaxy stellar content of halos ($\rm\mu_{\star,true}$) is a
powerful predictor of the cluster halo mass (Behroozi et al., 2010; Kravtsov
et al., 2018; Bradshaw et al., 2020). Let us consider the case where all the
true members inside $R_{\mathrm{200c}}$ are known and focus on the actual
$\rm\mu_{\star,true}$– $M_{200,c}$ relation, i.e., without photo-z errors.

Here, we can see the ideal stellar-mass-based mass proxy. The grey points in
Figure 8 show the scatter, and the $\rm\mu_{\star,true}$– $M_{200,c}$ relation
exhibits a low scatter. At the lowest redshift bins, $\sigma$ is about
$0.05\;\text{dex}$, a comparable value to the intra-cluster medium, e.g.
$Y_{\rm SZ}$, $M_{\rm gas}$ (Kravtsov et al., 2006; Mulroy et al., 2019; Bleem
et al., 2020; Ho et al., 2023), and redshift-based mass-proxies, e.g. velocity
dispersions (Ruel et al., 2014). At higher redshifts $z>0.32$, the scatter
increases, though we believe this is a simulation effect as the switch of box
resolution of the Buzzard simulation occurs at $z=0.33$.

On the other hand, the cluster aperture significantly affects the true
scatter, as in Figure 8, the mean scatter values (grey dashed lines) decrease
over the columns. The best aperture is the $1.00\times R_{\mathrm{200c}}$ as
expected by our physical knowledge of halos – the region that encloses the
overall virialized galaxy population should best predict the cluster mass. In
contrast, the smaller the aperture weakens the correlation of the visible
member stellar masses with the halo mass.

The stellar mass–halo mass scaling relation from Buzzard is in overall
agreement with four hydrodynamical simulations (Anbajagane et al., 2020). In
particular, the intrinsic scatter from the Bahamas simulations (Anbajagane et
al., 2020) at $z=0$ for high mass halos is close to $0.05$ dex, comparable
with the lowest redshift result value in Figure 10. On the other hand, the
slope has a slightly lower mean value, see Figure 9. Furthermore, a lower
slope value was also reported for richness in simulations compared to the
redMaPPer DES Y1 measured relations (DeRose et al., 2019).

##  5 Conclusions

This work presents the Copacabana algorithm, which assigns membership
probabilities to galaxies in a given cluster field. We validated the algorithm
using the Buzzard simulation. In particular, the stellar mass of cluster
galaxies weighted by the membership probabilities, $\rm\mu_{\star}$, was found
to have an accuracy up to $0.06$ dex for photometric redshifts that have
photo-z uncertainties comparable to that expected in LSST data. In addition,
we show that our methodology could precisely recover the scaling relation
between $\rm\mu_{\star}$ and cluster mass, indicating that $\rm\mu_{\star}$
can be a competitive mass-proxy for optically selected clusters in future
cosmological analysis.

  * •

Performance: we show the uncertainty in $\rm\mu_{\star}$ is mainly affected by
the quality of the photometric redshifts, followed by the cluster aperture. In
the best case, photometric redshifts with LSST-like accuracies
($\sigma_{z,0}=0.01$), we reported a $\rm\mu_{\star}$ uncertainty of $0.09$
dex. A smaller aperture, for instance, $0.5\times R_{\mathrm{200c}}$, leads to
a smaller uncertainty $\sim 0.05$ dex.

  * •

Cluster Size: We present a new method to measure cluster size,
$R_{\mathrm{200c}}$, with only photometric data. The procedure is based on the
HOD relation. We inferred a precision of $30\%$ in the context of the Buzzard
simulations. The accuracy of our estimate does not depend on the quality of
the photometric redshifts or the halo mass.

  * •

Impact on halo mass estimations: We quantify how the $\rm\mu_{\star}$
uncertainty propagates to estimates of the halo mass. Specifically, we study
the parameters of the scaling relation $\mu_{\star}-M_{200c}$ with a focus on
the slope and the intrinsic scatter. The photo-z uncertainty is the main
parameter in affecting the scaling relation parameters. For instance, in the
LSST-like scenario, we recover the parameters with no significant difference
compared to the true relation. In contrast, in the scenario of large photo-z
uncertainties and large apertures, there was a significant impact on
recovering the true parameters.

  * •

How the cluster aperture impacts the halo mass: While the accuracy of the
photometric redshifts is survey-dependent, the aperture size can be adjusted
to suit the scientific objectives. For example, a smaller cluster aperture can
significantly improve the membership probabilities. However, that improvement
does not necessarily translate into gains in predicting the halo mass. For
instance, we do not find a substantial improvement in the scatter of the
scaling relation, and the precision gained by a smaller aperture is
counterbalanced by the larger intrinsic scatter in $\mu_{\star}-M_{200c}$,
which is minimal at $R_{\mathrm{200c}}$. On the other hand, we find that the
aperture size significantly affects the slope. In particular, at an aperture
of $0.5\times R_{\mathrm{200c}}$, the recovered slope is very close to the
true slope and is insensitive to the range of photometric redshift
uncertainties that one encounters in modern photometric surveys. While at
$R_{\mathrm{200c}}$, the recovered slope can be up to $20\%$ shallower.

  * •

Membership probability performance: We report our galaxy member selection in
terms of Purity (P) and Completeness (C). In our best scenario, the accuracy
was $81\%$ with P and C of $64\%$ and $93\%$; when we consider a smaller
aperture ($0.5\times R_{\mathrm{200c}}$), these values were $89\%$, $79\%$,
and $94\%$, respectively. The membership probabilities substantially improved
with a smaller cluster aperture.

In sum, Copacabana is a powerful tool to predict the total stellar mass
content of galaxies in clusters. In future work, we will apply the stellar
mass-proxy $\rm\mu_{\star}$ in a cosmological analysis using optical data.

## Acknowledgements

JHE and MSS are funded by the U.S. Department of Energy (DOE) under grant No.
DE-SC0019193. MESP is funded by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) under Germany’s Excellence Strategy – EXC 2121
“Quantum Universe” – 390833306.

This research uses resources of the National Energy Research Scientific
Computing Center (NERSC), a U.S. Department of Energy Office of Science User
Facility located at Lawrence Berkeley National Laboratory.

Funding for the DES Projects has been provided by the U.S. Department of
Energy, the U.S. National Science Foundation, the Ministry of Science and
Education of Spain, the Science and Technology Facilities Council of the
United Kingdom, the Higher Education Funding Council for England, the National
Center for Supercomputing Applications at the University of Illinois at
Urbana-Champaign, the Kavli Institute of Cosmological Physics at the
University of Chicago, the Center for Cosmology and Astro-Particle Physics at
the Ohio State University, the Mitchell Institute for Fundamental Physics and
Astronomy at Texas A&M University, Financiadora de Estudos e Projetos,
Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro,
Conselho Nacional de Desenvolvimento Científico e Tecnológico and the
Ministério da Ciência, Tecnologia e Inovação, the Deutsche
Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy
Survey.

The Collaborating Institutions are Argonne National Laboratory, the University
of California at Santa Cruz, the University of Cambridge, Centro de
Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the
University of Chicago, University College London, the DES-Brazil Consortium,
the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH)
Zürich, Fermi National Accelerator Laboratory, the University of Illinois at
Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the
Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory,
the Ludwig-Maximilians Universität München and the associated Excellence
Cluster Universe, the University of Michigan, NSF’s NOIRLab, the University of
Nottingham, The Ohio State University, the University of Pennsylvania, the
University of Portsmouth, SLAC National Accelerator Laboratory, Stanford
University, the University of Sussex, Texas A&M University, and the OzDES
Membership Consortium.

Based in part on observations at Cerro Tololo Inter-American Observatory at
NSF’s NOIRLab (NOIRLab Prop. ID 2012B-0001; PI: J. Frieman), which is managed
by the Association of Universities for Research in Astronomy (AURA) under a
cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation
under Grant Numbers AST-1138766 and AST-1536171. The DES participants from
Spanish institutions are partially supported by MICINN under grants
ESP2017-89838, PGC2018-094773, PGC2018-102021, SEV-2016-0588, SEV-2016-0597,
and MDM-2015-0509, some of which include ERDF funds from the European Union.
IFAE is partially funded by the CERCA program of the Generalitat de Catalunya.
Research leading to these results has received funding from the European
Research Council under the European Union’s Seventh Framework Program
(FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478. We
acknowledge support from the Brazilian Instituto Nacional de Ciência e
Tecnologia (INCT) do e-Universo (CNPq grant 465376/2014-2).

This manuscript has been authored by Fermi Research Alliance, LLC under
Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of
Science, Office of High Energy Physics.

## Data Availability

The data used in this work is available upon request.

## References

  * Abbott et al. (2020) Abbott T. M. C., et al., 2020, Phys. Rev. D, 102, 023509
  * Aguena et al. (2021) Aguena M., et al., 2021, MNRAS, 502, 4435
  * Albrecht et al. (2006) Albrecht A., et al., 2006, arXiv e-prints, pp astro–ph/0609591
  * Allen et al. (2011) Allen S. W., Evrard A. E., Mantz A. B., 2011, ARA&A, 49, 409
  * Anbajagane et al. (2020) Anbajagane D., Evrard A. E., Farahi A., Barnes D. J., Dolag K., McCarthy I. G., Nelson D., Pillepich A., 2020, MNRAS, 495, 686
  * Andreon (2012) Andreon S., 2012, A&A, 548, A83
  * Andreon et al. (2006) Andreon S., Quintana H., Tajer M., Galaz G., Surdej J., 2006, MNRAS, 365, 915
  * Behroozi et al. (2010) Behroozi P. S., Conroy C., Wechsler R. H., 2010, ApJ, 717, 379
  * Bellagamba et al. (2018) Bellagamba F., Roncarelli M., Maturi M., Moscardini L., 2018, MNRAS, 473, 5221
  * Bellagamba et al. (2019) Bellagamba F., Sereno M., et al., 2019, MNRAS, 484, 1598
  * Bleem et al. (2020) Bleem L. E., Bocquet S., Stalder B., Gladders M. D., Ade P. A. R., Allen S. W., 2020, ApJS, 247, 25
  * Bradshaw et al. (2020) Bradshaw C., Leauthaud A., Hearin A., Huang S., Behroozi P., 2020, MNRAS, 493, 337
  * Butcher & Oemler (1984) Butcher H., Oemler A. J., 1984, ApJ, 285, 426
  * Carliles et al. (2010) Carliles S., Budavári T., Heinis S., Priebe C., Szalay A. S., 2010, ApJ, 712, 511
  * Castignani & Benoist (2016) Castignani G., Benoist C., 2016, A&A, 595, A111
  * Conroy & Gunn (2010) Conroy C., Gunn J. E., 2010, ApJ, 712, 833
  * Conroy et al. (2006) Conroy C., Wechsler R. H., Kravtsov A. V., 2006, ApJ, 647, 201
  * Conroy et al. (2009) Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486
  * Cooper et al. (2007) Cooper M. C., et al., 2007, MNRAS, 376, 1445
  * Costanzi et al. (2019) Costanzi M., et al., 2019, MNRAS, 482, 490
  * De Lucia et al. (2007) De Lucia G., et al., 2007, MNRAS, 374, 809
  * DeRose et al. (2019) DeRose J., Wechsler R. H., et al., 2019, arXiv e-prints, p. arXiv:1901.02401
  * DeRose et al. (2022) DeRose J., et al., 2022, Phys. Rev. D, 105, 123520
  * Dodelson et al. (2016) Dodelson S., Heitmann K., Hirata C., Honscheid K., Roodman A., Seljak U., Slosar A., Trodden M., 2016, arXiv e-prints, p. arXiv:1604.07626
  * Evrard et al. (2014) Evrard A. E., Arnault P., Huterer D., Farahi A., 2014, MNRAS, 441, 3562
  * George et al. (2011) George M. R., Leauthaud A., et al., 2011, ApJ, 742, 125
  * Golden-Marx et al. (2023) Golden-Marx J. B., et al., 2023, MNRAS, 521, 478
  * Gonzalez et al. (2007) Gonzalez A. H., Zaritsky D., Zabludoff A. I., 2007, ApJ, 666, 147
  * Gonzalez et al. (2013) Gonzalez A. H., Sivanandam S., Zabludoff A. I., Zaritsky D., 2013, ApJ, 778, 14
  * Gschwend et al. (2018) Gschwend J., et al., 2018, Astronomy and Computing, 25, 58
  * Hennig et al. (2017) Hennig C., Mohr J. J., Zenteno A., et al., 2017, MNRAS, 467, 4015
  * Ho et al. (2023) Ho M., Soltis J., Farahi A., Nagai D., Evrard A., Ntampaka M., 2023, arXiv e-prints, p. arXiv:2303.00005
  * Huang et al. (2020) Huang S., et al., 2020, MNRAS, 492, 3685
  * Huang et al. (2022) Huang S., et al., 2022, MNRAS, 515, 4722
  * Huterer (2023) Huterer D., 2023, A&ARv, 31, 2
  * Kelly (2007) Kelly B. C., 2007, ApJ, 665, 1489
  * Klein et al. (2022) Klein M., Oguri M., Mohr J. J., Grandis S., Ghirardini V., Liu T., Liu A., Bulbul E., 2022, A&A, 661, A4
  * Koester et al. (2007) Koester B. P., et al., 2007, ApJ, 660, 239
  * Kravtsov et al. (2006) Kravtsov A. V., Vikhlinin A., Nagai D., 2006, ApJ, 650, 128
  * Kravtsov et al. (2018) Kravtsov A. V., Vikhlinin A. A., Meshcheryakov A. V., 2018, Astronomy Letters, 44, 8
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration et al., 2009, arXiv e-prints, p. arXiv:0912.0201
  * Laganá et al. (2013) Laganá T. F., Martinet N., Durret F., Lima Neto G. B., Maughan B., Zhang Y. Y., 2013, A&A, 555, A66
  * Laigle et al. (2016) Laigle C., et al., 2016, ApJS, 224, 24
  * Leauthaud et al. (2012) Leauthaud A., et al., 2012, ApJ, 744, 159
  * Lehmann et al. (2017) Lehmann B. V., Mao Y.-Y., Becker M. R., Skillman S. W., Wechsler R. H., 2017, ApJ, 834, 37
  * Lopes & Ribeiro (2020) Lopes P. A. A., Ribeiro A. L. B., 2020, MNRAS, 493, 3429
  * McClintock et al. (2019) McClintock T., et al., 2019, MNRAS, 482, 1352
  * Mulroy et al. (2019) Mulroy S. L., et al., 2019, MNRAS, 484, 60
  * Myles et al. (2021) Myles J., et al., 2021, MNRAS, 505, 33
  * Palmese et al. (2016) Palmese A., Lahav O., et al., 2016, MNRAS, 463, 1486
  * Palmese et al. (2020) Palmese A., et al., 2020, MNRAS, 493, 4591
  * Pereira et al. (2018) Pereira M. E. S., Soares-Santos M., Makler M., Annis J., Lin H., Palmese A., Vitorelli A. Z., ohters 2018, MNRAS, 474
  * Pereira et al. (2020) Pereira M. E. S., Palmese A., Soares-Santos M., DES Collaboration 2020, MNRAS, 498, 5450
  * Pratt et al. (2019) Pratt G. W., Arnaud M., Biviano A., Eckert D., Ettori S., Nagai D., Okabe N., Reiprich T. H., 2019, Space Sci. Rev., 215, 25
  * Puddu et al. (2021) Puddu E., et al., 2021, A&A, 645, A9
  * Reddick et al. (2013) Reddick R. M., Wechsler R. H., Tinker J. L., Behroozi P. S., 2013, ApJ, 771, 30
  * Rozo et al. (2009) Rozo E., Rykoff E. S., et al., 2009, ApJ, 699, 768
  * Ruel et al. (2014) Ruel J., et al., 2014, ApJ, 792, 45
  * Rykoff et al. (2012) Rykoff E. S., et al., 2012, ApJ, 746, 178
  * Rykoff et al. (2014) Rykoff E. S., et al., 2014, ApJ, 785, 104
  * Rykoff et al. (2016) Rykoff E. S., Rozo E., DES Collaboration 2016, ApJS, 224, 1
  * Sarazin (1988) Sarazin C. L., 1988,
  * Scott (1992) Scott D. W., 1992, John Wiley & Sons
  * Song et al. (2017) Song H., Hwang H. S., Park C., Tamura T., 2017, ApJ, 842, 88
  * Springel et al. (2005) Springel V., et al., 2005, Nature, 435, 629
  * Sunayama et al. (2020) Sunayama T., et al., 2020, MNRAS, 496, 4468
  * Taylor et al. (2011) Taylor E. N., et al., 2011, MNRAS, 418, 1587
  * The Dark Energy Survey Collaboration (2005) The Dark Energy Survey Collaboration 2005, arXiv e-prints, pp astro–ph/0510346
  * Tinker et al. (2012) Tinker J. L., Sheldon E. S., Wechsler R. H., et al., 2012, ApJ, 745, 16
  * Umetsu (2020) Umetsu K., 2020, A&ARv, 28, 7
  * Wechsler et al. (2022) Wechsler R. H., DeRose J., Busha M. T., Becker M. R., Rykoff E., Evrard A., 2022, ApJ, 931, 145
  * Wright & Brainerd (2000) Wright C. O., Brainerd T. G., 2000, ApJ, 534, 34
  * Wu et al. (2021) Wu H.-Y., Weinberg D. H., Salcedo A. N., Wibking B. D., 2021, ApJ, 910, 28
  * Wu et al. (2022) Wu H.-Y., et al., 2022, MNRAS, 515, 4471
  * York et al. (2000) York D. G., et al., 2000, AJ, 120, 1579
  * Zacharegkas et al. (2022) Zacharegkas G., et al., 2022, MNRAS, 509, 3119

##  Appendix A Affiliations

1 Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA
2 Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, 21029 Hamburg,
Germany
3 Fermi National Accelerator Laboratory, P. O. Box 500, Batavia, IL 60510, USA
4 Departments of Statistics and Data Science, University of Texas at Austin,
Austin, TX 78757, USA
5
6 Department of Physics, Carnegie Mellon University, Pittsburgh, Pennsylvania
15312, USA
7 Observational Cosmology Lab, NASA Goddard Space Flight Center, Greenbelt, MD
20771, USA
8 Department of Astronomy, University of Maryland, College Park, 20742, USA
9 Department of Physics, Boise State University, Boise, ID 83725, USA
10 Laboratório Interinstitucional de e-Astronomia - LIneA, Rua Gal. José
Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil
11 Institute of Cosmology and Gravitation, University of Portsmouth,
Portsmouth, PO1 3FX, UK
12 University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität,
Scheinerstr. 1, 81679 Munich, Germany
13 Department of Physics & Astronomy, University College London, Gower Street,
London, WC1E 6BT, UK
14 Instituto de Astrofisica de Canarias, E-38205 La Laguna, Tenerife, Spain
15 Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife,
Spain
16 Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of
Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona) Spain
17 Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo
11, I-34131 Trieste, Italy
18 INAF-Osservatorio Astronomico di Trieste, via G. B. Tiepolo 11, I-34143
Trieste, Italy
19 Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014
Trieste, Italy
20 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas
(CIEMAT), Madrid, Spain
21 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak
Grove Dr., Pasadena, CA 91109, USA
22 Kavli Institute for Cosmological Physics, University of Chicago, Chicago,
IL 60637, USA
23 Instituto de Fisica Teorica UAM/CSIC, Universidad Autonoma de Madrid, 28049
Madrid, Spain
24 Center for Astrophysical Surveys, National Center for Supercomputing
Applications, 1205 West Clark St., Urbana, IL 61801, USA
25 Department of Astronomy, University of Illinois at Urbana-Champaign, 1002
W. Green Street, Urbana, IL 61801, USA
26 School of Mathematics and Physics, University of Queensland, Brisbane, QLD
4072, Australia
27 Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA
28 Center for Cosmology and Astro-Particle Physics, The Ohio State University,
Columbus, OH 43210, USA
29 Department of Physics, The Ohio State University, Columbus, OH 43210, USA
30 Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA
31 Australian Astronomical Optics, Macquarie University, North Ryde, NSW 2113,
Australia
32 Lowell Observatory, 1400 Mars Hill Rd, Flagstaff, AZ 86001, USA
33 Centre for Gravitational Astrophysics, College of Science, The Australian
National University, ACT 2601, Australia
34 The Research School of Astronomy and Astrophysics, Australian National
University, ACT 2601, Australia
35 Departamento de Física Matemática, Instituto de Física, Universidade de São
Paulo, CP 66318, São Paulo, SP, 05314-970, Brazil
36 George P. and Cynthia Woods Mitchell Institute for Fundamental Physics and
Astronomy, and Department of Physics and Astronomy, Texas A&M University,
College Station, TX 77843, USA
37 LPSC Grenoble - 53, Avenue des Martyrs 38026 Grenoble, France
38 Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain
39 Department of Astrophysical Sciences, Princeton University, Peyton Hall,
Princeton, NJ 08544, USA
40 Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ -
20921-400, Brazil
41 Kavli Institute for Particle Astrophysics & Cosmology, P. O. Box 2450,
Stanford University, Stanford, CA 94305, USA
42 SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA
43 Department of Physics and Astronomy, Pevensey Building, University of
Sussex, Brighton, BN1 9QH, UK
44 Instituto de Física, UFRGS, Caixa Postal 15051, Porto Alegre, RS -
91501-970, Brazil
45 School of Physics and Astronomy, University of Southampton, Southampton,
SO17 1BJ, UK
46 Computer Science and Mathematics Division, Oak Ridge National Laboratory,
Oak Ridge, TN 37831
47 Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA
94720, USA
48 Department of Physics, Duke University Durham, NC 27708, USA


Exactly solvable complex $\mathcal{PT}$ symmetry potential
$A[\operatorname{sech}(\lambda x)+i\tanh(\lambda x)]$

Wei Yang

College of Science, Guilin University of Technology, Guilin, Guangxi 541004,
China

ABSTRACT

We obtained the bound and scattering solutions of the $\mathcal{PT}$ symmetric
potential $V(x)=A[\operatorname{sech}(\lambda x)+i\tanh(\lambda x)]$, and its
bound-state energy spectrum does not turn into complex conjugate pairs, i.e.
the spontaneous breakdown of $\mathcal{PT}$ symmetry does not occur. We found
that the reflection coefficients exhibit a handedness effect. As the
asymptotically non-vanishing imaginary potential component, that the
transmission coefficient increases a complex phase factor when the direction
of the incoming wave is reversed.

Emails<EMAIL_ADDRESS>
##  1 Introduction

Since Bender and Boettcher [1] introduced $\mathcal{PT}$ symmetric quantum
mechanics in 1998, the study of non-Hermitian quantum mechanical systems has
developed rapidly. They showed that quantum systems with a non-Hermitian
Hamiltonian can have a set of eigenstates with real eigenvalues. In other
words, they found that the Hermiticity of the Hamiltonian is not a necessary
condition for the realness of its eigenvalues, and new quantum mechanics can
be constructed based on such Hamiltonians [2, 3].

In one-dimensional Schrödinger equation with $\mathcal{PT}$ symmetry, that is
the invariance of the Hamiltonian under the combined action of the space and
time inversion, due to the relation for the potential given $V(x)=V(-x)^{*}$.
Later $\mathcal{PT}$ symmetric quantum mechanics was identified as a special
case of pseudo-hermiticity [4, 5, 6]. A class of typical $\mathcal{PT}$
symmetric potentials $V(x)=x^{2}(ix)^{\epsilon}$ were studied numerically,
show these potential have the real and positive energy spectrum for
$\epsilon\geq 0$, when $\epsilon<0$, the energy eigenvalues gradually merge
pairwise and reemerge as complex conjugate pairs. This mechanism was
interpreted as the spontaneous breakdown of $\mathcal{PT}$ symmetry, because
the eigenfunctions is no longer the eigenfunctions of $\mathcal{PT}$ operator.

In addition to numerical research, there is still a lot of analytical work [7,
8, 9, 10] to solve the $\mathcal{PT}$ symmetric potentials. In fact, more
generally, we can construct $\mathcal{PT}$ symmetric potentials from Natanzon
class [11, 12]. Here we study the $\mathcal{PT}$ symmetric potential
$V(x)=A(\operatorname{sech}(\lambda x)+i\tanh(\lambda x))$, defined on the
real $x$-axis, and investigate both its bound and scattering states. This
$\mathcal{PT}$ symmetric potential has not been considered in previous
articles, and its imaginary component asymptotically non-vanishing. This
similar potential is also considered in the article[13, 14].

The work is organized as follows. In Sec. 2, we will give the exact solution
for this potential. In Sec. 3, use the obtained exact solution to building-up
scattering information. In Sec. 4, we obtained the bound states and energy
eigenvalues, further more we found the bound states are the poles of the
transmission and reflection coefficients. In Sec. 5, we will give some
conclusions for this paper.

##  2 Exact solution

In this paper, we will consider the $\mathcal{PT}$ symmetric potential
$V(x)=A[\operatorname{sech}(\lambda x)+i\tanh(\lambda x)]$, then the one-
dimension Schrödinger equation write as

$\displaystyle[\frac{d^{2}}{dx^{2}}+E-V(x)]\psi(x)=0\,.$ (1)

Define $z=ie^{\lambda x}$, which transform the variable $x\in(-\infty,\infty)$
to a new variable $z\in(0,\infty)$, and introduce a function $u(z)$, define as
$\psi(z)=(-z)^{p}u(z)$, by means of this substitution, we can rewrite the
Schrödinger equation as

$\displaystyle
z(1-z)u^{\prime\prime}(z)+[1+2p-(1+2p)z]u^{\prime}(z)+\frac{E-iA+p^{2}\lambda^{2}}{\lambda^{2}}u(z)=0\,.$
(2)

If and only if $E+iA+p^{2}\lambda^{2}=0$, Namely

$\displaystyle p=\pm\sqrt{-iA-E}/\lambda\,.$ (3)

We can write the solution by comparing with hypergeometric equation

$\displaystyle
z(1-z)u^{\prime\prime}(z)+[\gamma-(\alpha+\beta+1)z]u^{\prime}(z)-\alpha\beta
u(z)=0\,.$ (4)

After simple calculation, we have

$\displaystyle\alpha=p\pm\sqrt{iA-E}/\lambda\equiv p\pm
q\,,\qquad\beta=2p-\alpha\,,\qquad\gamma=2p+1\,.$ (5)

Because of the commutative symmetry of $\alpha$ and $\beta$, without losing
generality, we can choose $\alpha=p-q$. From the hypergeometric functions
theory, the wave function is given by

$\displaystyle\psi(z)$
$\displaystyle=(-z)^{p}[c_{1}F(\alpha,\beta,\gamma;z)+c_{2}(-z)^{1-\gamma}F(\alpha-\gamma+1,\beta-\gamma+1,2-\gamma;z)]\,.$
(6)

Where $c_{1},c_{2}$ are two arbitrary combination coefficients. Due to the
positive and negative sign of $p$ cause the exchange of left and right
travelling waves of wave function. In this we can choose
$p=\sqrt{-iA-E}/\lambda$.

##  3 Scattering states

Considering the parametrization

$\displaystyle\sqrt{E+iA}=k_{R}+ik_{I}=k\,.\,\qquad\sqrt{E-iA}=k^{\prime}_{R}+ik^{\prime}_{I}=k^{\prime}$
(7)

and assuming that $E,A$ are real numbers, the relations

$\displaystyle
E=k_{R}^{2}-k_{I}^{2}={k^{\prime}}_{R}^{2}-{k^{\prime}}_{I}^{2}\,\qquad
A=2k_{R}k_{I}=-2{k^{\prime}}_{R}{k^{\prime}}_{I}\,.$ (8)

So we have $k^{2}=[(k^{\prime})^{2}]^{*}$, which implies that $k^{\prime}=\pm
k^{*}$.

For $A>0$, choice $k^{\prime}=k^{*}$ this mean $k_{R}=k^{\prime}_{R}>0$ and
$k_{I}=-k^{\prime}_{I}>0$ (For $A<0$ the roles of $k$ and $k^{\prime}$ are
exchanged). Let us analyse the asymptotic behaviour of solutions of (6) for
$x\rightarrow-\infty$ (or $z\rightarrow 0$), we have

$\displaystyle\psi(x)_{-\infty}\sim
c_{1}(-z)^{p}+c_{2}(-z)^{-p}=(-i)^{p}c_{1}e^{ikx}+(-i)^{-p}c_{2}e^{-ikx}\,.$
(9)

So we have

$\displaystyle a_{1-}=(-i)^{p}\,\qquad a_{2-}=0\,\qquad b_{1-}=0\,\qquad
b_{2-}=(-i)^{-p},\,.$ (10)

where the notation can be found in [15] or see appendix A.

And the asymptotic behaviour of the wave function $\psi(x)$ for
$x\rightarrow\infty$ (or $z\rightarrow\infty$) is given

$\displaystyle\psi(x)_{\infty}$
$\displaystyle=[c_{1}\frac{\Gamma(1+2p)\Gamma(-2q)}{\Gamma(p-q)\Gamma(1+p-q)}+c_{2}\frac{\Gamma(1-2p)\Gamma(-2q)}{\Gamma(-p-q)\Gamma(1-p-q)}](-z)^{-q}$
$\displaystyle+[c_{1}\frac{\Gamma(1+2p)\Gamma(2q)}{\Gamma(p+q)\Gamma(1+p+q)}+c_{2}\frac{\Gamma(1-2p)\Gamma(2q)}{\Gamma(-p+q)\Gamma(1-p+q)}](-z)^{q}\,.$
(11)

where we used the formula

$\displaystyle
F(\alpha,\beta,\gamma;z)=\frac{\Gamma(\gamma)\Gamma(\beta-\alpha)}{\Gamma(\beta)\Gamma(\gamma-\alpha)}(-z)^{-\alpha}F(\alpha,\alpha-\gamma+1,\alpha-\beta+1;1/z)$
$\displaystyle+\frac{\Gamma(\gamma)\Gamma(\alpha-\beta)}{\Gamma(\alpha)\Gamma(\gamma-\beta)}(-z)^{-\beta}F(\beta,\beta-\gamma+1,\beta-\alpha+1;1/z)\,.$
(12)

We have similar notation, which are expressed

$\displaystyle b_{1+}$
$\displaystyle=(-i)^{-q}\frac{\Gamma(1+2p)\Gamma(-2q)}{\Gamma(p-q)\Gamma(1+p-q)}\,,\qquad
b_{2+}=(-i)^{-q}\frac{\Gamma(1-2p)\Gamma(-2q)}{\Gamma(-p-q)\Gamma(1-p-q)}$
$\displaystyle a_{1+}$
$\displaystyle=(-i)^{-q}\frac{\Gamma(1+2p)\Gamma(2q)}{\Gamma(p+q)\Gamma(1+p+q)}\,,\qquad
a_{2+}=(-i)^{-q}\frac{\Gamma(1-2p)\Gamma(2q)}{\Gamma(-p+q)\Gamma(1-p+q)}\,.$
(13)

From Eq. (10) and Eq. (13), we have the transmission and reflection
coefficients for a wave coming from the right are

$\displaystyle R_{R\rightarrow
L}=(-i)^{2q}\frac{\Gamma(-p-q)\Gamma(1-p-q)\Gamma(2q)}{\Gamma(-2q)\Gamma(-p+q)\Gamma(1-p+q)}\,,$
$\displaystyle T_{R\rightarrow
L}=(-i)^{q-p}\frac{\Gamma(-p-q)\Gamma(1-p-q)}{\Gamma(1-2p)\Gamma(-2q)}\,.$
(14)

The transmission and reflection coefficients for a wave coming from the left
are also given

$\displaystyle R_{L\rightarrow
R}=-e^{ip\pi}\frac{\Gamma(1+2p)\Gamma(-p-q)\Gamma(1-p-q)}{\Gamma(p-q)\Gamma(1+p-q)\Gamma(1-2p)}\,,$
$\displaystyle T_{L\rightarrow R}={\frac{p}{q}}T_{R\rightarrow L}\,.$ (15)

we see that $R_{R\rightarrow L}\neq R_{L\rightarrow R}$, this reflects the
handedness of the problem. As the imaginary part of potential asymptotically
non-vanishing, that have $T_{L\rightarrow R}={\frac{p}{q}}T_{R\rightarrow L}$.

Here we point out for $A>0$, if let $k^{\prime}=-k^{*}$ this mean
$k_{R}=-k^{\prime}_{R}>0$ and $k_{I}=k^{\prime}_{I}>0$ (For $A<0$ the roles of
$k$ and $k^{\prime}$ are exchanged). The corresponding transmission and
reflection coefficients is expressed as

$\displaystyle R_{R\rightarrow L}=e^{i\pi
q}\frac{\Gamma(-p+q)\Gamma(1-p+q)\Gamma(-2q)}{\Gamma(2q)\Gamma(-p-q)\Gamma(1-p-q)}\,,$
$\displaystyle T_{R\rightarrow
L}=e^{i\pi(p+q)/2}\frac{\Gamma(-p+q)\Gamma(1-p+q)}{\Gamma(1-2p)\Gamma(2q)}\,,$
$\displaystyle R_{L\rightarrow
R}=-e^{ip\pi}\frac{\Gamma(1+2p)\Gamma(-p+q)\Gamma(1-p+q)}{\Gamma(p+q)\Gamma(1+p+q)\Gamma(1-2p)}\,,$
$\displaystyle T_{L\rightarrow R}=-{\frac{p}{q}}T_{R\rightarrow L}\,.$ (16)

We can see that we just need to change $q$ to $-q$ in Eq. (16), the expression
will return to the case of $k^{\prime}=-k^{*}$.

##  4 Bound states

Let us find the solutions of the eigenvalue problem associated to this
potential, One way to obtain those solutions is to analyze the behavior of the
eigenfunctions in the limits $x\rightarrow-\infty$. Without loss of
generality, we choose $\text{Re}p<0$, this corresponds to the case of
$k_{I}=-k^{\prime}_{I}>0$, it is required that eigenfunctions square
integrable, so we have $\psi(x)\sim
z^{-p}F(\alpha-\gamma+1,\beta-\gamma+1,2-\gamma;z)$. Similarly, when
$x\rightarrow\infty$, there is bound states only when $-p\pm q=-n$

$\displaystyle\psi_{n}(x)\sim
z^{-p}F(-n,\beta-\gamma+1,2-\gamma;z)={\frac{z^{-p}n!}{(1-2p)_{n}}}P_{n}^{-2p,-1}(1-2z)\,.$
(17)

where $P_{n}^{a,b}(z)$ is Jacobi polynomials. And energy eigenvalues

$\displaystyle
E_{n}={\frac{A^{2}}{n^{2}\lambda^{2}}}-{\frac{n^{2}\lambda^{2}}{4}}\,.$ (18)

we can see the conditions for bound states $-p\pm q=-n$ are the poles of the
transmission and reflection coefficients. By the way, we choose
$\text{Re}p>0$, this corresponds to the case of $k_{I}=-k^{\prime}_{I}<0$, we
will have conditions for bound states $p\pm q=-n$, the wave functions give

$\displaystyle\psi_{n}(x)\sim
z^{p}F(-n,\beta,\gamma;z)={\frac{z^{p}n!}{(2p+1)_{n}}}P_{n}^{2p,-1}(1-2z)\,.$
(19)

And we also have the same energy eigenvalues
$E_{n}={\frac{A^{2}}{n^{2}\lambda^{2}}}-{\frac{n^{2}\lambda^{2}}{4}}$.

##  5 Conclusion

In this work, we consider solving Schrödinger equations with the potential
$V(x)=A[\operatorname{sech}(\lambda x)+i\tanh(\lambda x)]$, and obtained for
the transmission and reflection coefficients of this $\mathcal{PT}$ symmetric
potential. Due to the asymmetric structure of the potential, the asymptotic
wave numbers differ in both cases for $x\rightarrow\infty$ and
$x\rightarrow-\infty$, and the wave numbers are complex for the $\mathcal{PT}$
symmetric potential. Changing the direction of the incoming wave will changes
the transmission and reflection coefficients. The transmission coefficient
changes by a complex phase factor, while the reflection coefficient changes
more drastically, corresponding to the handedness effect observed for local
$\mathcal{PT}$ symmetric scattering potentials in general. By the way, we also
can consider the $\mathcal{PT}$ symmetric potential
$V(x)=A[\operatorname{sech}(\lambda x)-i\tanh(\lambda x)]$, that only need
redefine $z=-ie^{\lambda x}$, and change the $iA$ to $-iA$ in the solutions of
(6). Obviously, the energy eigenvalues corresponding to this potential is also
$E_{n}={\frac{A^{2}}{n^{2}\lambda^{2}}}-{\frac{n^{2}\lambda^{2}}{4}}$.

## Acknowledgement

This work was supported by the Guangxi Scientific Programm Foundation under
grant No. 2020AC20014, the Scientific Research Foundation of Guilin University
of Technology under grant No. GUTQDJJ2019206.

##  Appendix A Transmission and reflection coefficients

Consider the time independent Schrödinger equation satisfied by

$\displaystyle H\psi(x)=[-\frac{d^{2}}{dx^{2}}+V(x)]\psi(x)=k^{2}\psi(x)\,.$
(20)

With the asymptotic states (finite-range local potential) can be expressed at
$x\rightarrow\pm\infty$ as

$\displaystyle\begin{split}\psi(x)=\left\\{\begin{array}[]{ll}A_{-}e^{ikx}+B_{-}e^{-ikx},&x\rightarrow-\infty\\\
A_{+}e^{ikx}+B_{+}e^{-ikx},&x\rightarrow-\infty\end{array}\right.\end{split}$

Schrödinger equation admits a general solution written as a linear combination
of two independent solutions, $F_{1}(x)$ and $F_{2}(x)$, with non-zero
Wronskian, whose asymptotic expressions are both of the form:

$\displaystyle\lim_{x\rightarrow\pm\infty}F_{1}=a_{1\pm}e^{ikx}+b_{1\pm}e^{-ikx}\,,$
$\displaystyle\lim_{x\rightarrow\pm\infty}F_{2}=a_{2\pm}e^{ikx}+b_{2\pm}e^{-ikx}\,.$
(21)

which related to the asymptotic amplitudes $A_{\pm}$ and $B_{\pm}$

$\displaystyle A_{\pm}=\alpha a_{1\pm}+\beta a_{2\pm}\,,$ $\displaystyle
B_{\pm}=\alpha b_{1\pm}+\beta b_{2\pm}\,.$ (22)

The transmission and reflection coefficients of the wave moving from left to
right are written as

$\displaystyle T_{L\rightarrow
R}={\frac{a_{2+}b_{1+}-a_{1+}b_{2+}}{a_{2-}b_{1+}-a_{1-}b_{2+}}}\,,$
$\displaystyle R_{L\rightarrow
R}={\frac{b_{1+}b_{2-}-b_{1-}b_{2+}}{a_{2-}b_{1+}-a_{1-}b_{2+}}}\,.$ (23)

The transmission and reflection coefficients of the wave moving from right to
left are

$\displaystyle T_{R\rightarrow
L}={\frac{a_{2-}b_{1-}-a_{1-}b_{2-}}{a_{2-}b_{1+}-a_{1-}b_{2+}}}\,,$
$\displaystyle R_{L\rightarrow
R}={\frac{a_{1+}a_{2-}-a_{1-}a_{2+}}{a_{2-}b_{1+}-a_{1-}b_{2+}}}\,.$ (24)

## References

  * [1] C. M. Bender and S. Boettcher, Phys. Rev. Lett. 80, 5243-5246 (1998)
  * [2] C. M. Bender, D. C. Brody and H. F. Jones, Phys. Rev. Lett. 89, 270401 (2002) [erratum: Phys. Rev. Lett. 92, 119902 (2004)]
  * [3] C. M. Bender, Rept. Prog. Phys. 70, 947 (2007)
  * [4] A. Mostafazadeh, J. Math. Phys. 43, 205-214 (2002)
  * [5] A. Mostafazadeh, J. Math. Phys. 43, 2814-2816 (2002)
  * [6] A. Mostafazadeh, J. Math. Phys. 43, 3944-3951 (2002)
  * [7] M. Znojil, Phys. Lett. A 259, 220-223 (1999)
  * [8] M. Znojil, J. Phys. A 33, L61-L62 (2000)
  * [9] M. Znojil and G. Levai, J. Phys. A 35, 8793 (2002)
  * [10] Z. Ahmed, Phys. Lett. A 282, 343-348 (2001)
  * [11] A.K. Bose, Phys. Lett. 7, 245 (1963).
  * [12] Natanzon, G. A.Theoretical and Mathematical Physics 38 (1979): 146-153.
  * [13] Géza, Lévai and E. Magyari, J. Phys. A, Math. Theor. 42, 195302 (2009)
  * [14] Géza, Lévai, Int J Theor Phys 50, 997-1004 (2011)
  * [15] F. Cannata, J. P. Dedonder and A. Ventura, Annals Phys. 322, 397-433 (2007)