††thanks: Equal contributions††thanks: Equal contributions

# Cu doping effects on the electronic structure of Fe1-xCuxSe

S. Huh  Department of Physics and Astronomy, Seoul National University (SNU),
Seoul 08826, Republic of Korea  Center for Correlated Electron Systems,
Institute for Basic Science (IBS), Seoul 08826, Republic of Korea  Z. Lu
Beijing National Laboratory for Condensed Matter Physics and Institute of
Physics, Chinese Academy of Sciences, Beijing 100190, China  School of
Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049,
China  Y. S. Kim  Department of Physics and Astronomy, Seoul National
University (SNU), Seoul 08826, Republic of Korea  Center for Correlated
Electron Systems, Institute for Basic Science (IBS), Seoul 08826, Republic of
Korea  D. Kim  Department of Physics and Astronomy, Seoul National University
(SNU), Seoul 08826, Republic of Korea  Center for Correlated Electron Systems,
Institute for Basic Science (IBS), Seoul 08826, Republic of Korea  S. B. Liu
Beijing National Laboratory for Condensed Matter Physics and Institute of
Physics, Chinese Academy of Sciences, Beijing 100190, China  School of
Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049,
China  M. W. Ma  Beijing National Laboratory for Condensed Matter Physics and
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
Songshan Lake Materials Laboratory, Dongguan, Guangdong 523808, China  L. Yu
<EMAIL_ADDRESS>Beijing National Laboratory for Condensed Matter Physics and
Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 100049, China  Songshan Lake Materials Laboratory, Dongguan, Guangdong
523808, China  F. Zhou  Beijing National Laboratory for Condensed Matter
Physics and Institute of Physics, Chinese Academy of Sciences, Beijing 100190,
China  School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 100049, China  Songshan Lake Materials Laboratory, Dongguan, Guangdong
523808, China  X. L. Dong  Beijing National Laboratory for Condensed Matter
Physics and Institute of Physics, Chinese Academy of Sciences, Beijing 100190,
China  School of Physical Sciences, University of Chinese Academy of Sciences,
Beijing 100049, China  Songshan Lake Materials Laboratory, Dongguan, Guangdong
523808, China  C. Kim<EMAIL_ADDRESS>Department of Physics and
Astronomy, Seoul National University (SNU), Seoul 08826, Republic of Korea
Center for Correlated Electron Systems, Institute for Basic Science (IBS),
Seoul 08826, Republic of Korea  Z. X. Zhao  Beijing National Laboratory for
Condensed Matter Physics and Institute of Physics, Chinese Academy of
Sciences, Beijing 100190, China  School of Physical Sciences, University of
Chinese Academy of Sciences, Beijing 100049, China  Songshan Lake Materials
Laboratory, Dongguan, Guangdong 523808, China

###### Abstract

Using angle-resolved photoemission spectroscopy (ARPES), we studied the
evolution of the electronic structure of Fe1-xCuxSe from x = 0 to 0.10. We
found that the Cu dopant introduces extra electron carriers. The hole bands
near the $\Gamma$ point are observed to steadily shift downward with
increasing doping and completely sink down below the Fermi level ($E_{F}$) for
x $>$ 0.05. Meanwhile, the electron pocket near the M point becomes larger but
loses the spectral weight near $E_{F}$. We also observed that effective mass
of the electron band near the M point increases with doping. Our result
explains why superconductivity disappears and metal insulator transition (MIT)
like behavior occurs upon Cu doping in terms of electronic structure, and
provide insight into emergent magnetic fluctuation in Fe1-xCuxSe.

##  I I. Introduction

Chemical substitution or doping studies have been widely conducted in
correlated materials to control the electronic structure as well as physical
properties. For example, it can change the chemical pressure, electron
correlation, Fermi energy, carrier density. As such, investigating the
chemical substitution effects can provide critical clues to understanding the
underlying mechanism for various phenomena such as emergent magnetism, metal-
insulator transition (MIT), spin/charge density wave, and superconductivity
[1, 2, 3, 4, 5, 6, 7].

Numerous chemical substitution studies have been conducted for iron based
superconductors (IBSs) to understand the superconductivity mechanism and the
interplay between the superconductivity and emergent phenomena [8, 9]. Among
IBSs, chemical substitution studies have been intensively done on FeSe. The
resultant changes of superconductivity and relevant electronic and magnetic
properties through doping have provided rich insights of the underlining
nature. Moreover, exotic phenomena such as Bardeen-Cooper-Schrieffer (BCS)
Bose-Einstein condensation (BEC) crossover and topological superconductivity
have been observed in substituted FeSe1-xSx and FeSe1-xTex systems [10, 11,
12, 13]. While these studies expanded the boundary of the physics in FeSe,
most of the studies have been limited to substitution of Se with isovalent
atoms. There are relatively less doping studies on the Fe site due to the
challenge in synthesizing high quality single crystals.

Investigating the Fe-site substitution may provide new opportunities for a
better understanding of superconductivity and exploration of emergent
phenomena. Cu doped FeSe is a suitable candidate to perform such studies,
because of the recent success of high quality single crystal synthesis [14].
Previous studies on Cu doped FeSe reported that the structural transition and
superconductivity are suppressed by a small amount of Cu doping and that a
metal-insulator transition (MIT) like behavior appears with further doping
[14, 15]. This unexpected MIT in FeSe call for systematic spectroscopic
studies on the Cu doping eﬀect [16]. In addition, emergence of magnetic
fluctuation [17] and restoration of TC near 30 K under hydrostatic pressure
[18] make this Fe1-xCuxSe even more interesting and require deeper
understanding of its properties in microscopic level. In elucidating the Cu
doping eﬀect and mechanism of new emergent phenomena, electronic structure
information should be vital, which have yet been lacking.

In this paper, we report the angle-resolved photoemission spectroscopy (ARPES)
measurements of the electronic structure of Fe1-xCuxSe (x = 0, 0.02, 0.03,
0.05, and 0.10). We ﬁnd that the Fermi surface is strongly Cu doping
dependent. The size of the hole (electron) pocket near the $\Gamma$ (M) point
decreases (increases) with doping. Upon the doping increased to x $>$ 0.05,
the hole bands near $\Gamma$ point sink down below the Fermi level ($E_{F}$),
which changes the overall Fermi surface topology. In addition, we observe the
spectral weight is reduced near $E_{F}$ as well as an increased effective mass
of the electron band near the M point. Our observation reveals for the first
time the electronic structure evolution of Fe1-xCuxSe system in detail. It
attributes the MIT behavior to electronic structure instead of merely impurity
scattering. Our results may provide microscopic experimental foundation for
understanding the emergent magnetic ﬂuctuation in this doping range.

##  II II. Experiment

High quality single crystals of Fe1-xCuxSe were synthesized by the chemical
vapor transport method using mixture of AlCl3 and KCl as transport agent [14].
In-plane electrical transport measurements were performed with a Quantum
Design PPMS-9 system. ARPES measurements were performed using a home lab based
system at Seoul National University. All spectra were acquired using a VG-
Scienta electron analyzer and a discharge lamp from Fermi instruments. 21.218
eV photon energy was used in the experiment. Instrumental energy resolution
was better than 20 meV. Samples were cleaved in an ultrahigh vacuum better
than $8\times 10^{-11}$ Torr. To minimize the aging effect, all data were
taken within 10 hours after the cleaving.

##  III III. Results and discussion

Figure 1: (Color online) (a) Optical image of Fe1-xCuxSe single crystals. (b)
Temperature dependent resistivity between 2 and 300 K. The inset shows the
derivative of the resistivity (x = 0 and 0.02) as a function of temperature.
The black arrow indicates structural transition (TS).

Figure 1(a) shows an optical image of Fe1-xCuxSe single crystals. Sizable
single crystals with various doping concentration used in our ARPES
experiments were obtained. The chemical composition of crystals was determined
by inductively coupled plasma atomic emission spectroscopy. Temperature
dependent in-plane resistivity for various doping is shown in Figure 1(b).
Consistent with previous results [14, 15], TC and structural transition
temperature (TS) are completely suppressed when Cu is only slightly doped (x
$<$ 0.02). For a better comparison, we plot the first derivative of the
resistivity in the inset. There is an anomaly near 86 K for x = 0 sample which
indicates TS, while no signature of such anomaly is observable at x = 0.02.
The heat capacity measurement gives a consistent TS value [14]. In addition to
the suppression of the phase transition, MIT like behavior is observed from
the doping dependence of resistivity data. An upturn in the resistivity curve
was observed for x = 0.05 and 0.10 while metallic behavior is observed over
the entire temperature range for lower doping concentrations. Similar MIT
behavior is reported for Cu doped iron pnictide systems [19, 20], which
implies that the Cu doping effect may be universal in IBSs.

We performed ARPES experiments to investigate the Cu doping effect on the
electronic structure of Fe1-xCuxSe. Figures 2(a)-(e) show Fermi surface maps
of Fe1-xCuxSe for different Cu concentrations. FeSe data were taken at 100 K
to avoid the change in the band structure from the nematic phase transition
which occurs near 90 K [21], while other data were taken at 20 K since no
structural transition is observed for these concentrations as mentioned above.
This allows us to systematically compare the band structures of Fe1-xCuxSe
with various x values. With increasing Cu doping, the hole pocket near
$\Gamma$ gradually decreases and disappear eventually when x $>$ 0.05. On the
other hand, the size of the electron pocket near $M_{1}$ and $M_{2}$ point
increases. The overall change of the Fermi surface topology suggests that Cu
substitution seems to systematically dope extra electron carrier into the FeSe
system.

A number of studies reported that Cu doping to IBSs may result in hole doping,
but it is certainly not the case for FeSe [22, 23]. In addition, the observed
electron doping behavior is similar to that of Ba(Fe1-xCux)2As2 and
NaFe1-xCuxAs [26, 24, 25]. It is noteworthy that the shape of the electron
pocket becomes more circular in comparison to the two perpendicularly crossing
elliptical shape Fermi surface pocket of pristine FeSe. Such big circular
electron pockets in absence of hole pocket at the center suggest the
interesting similarity of the FS topology for heavily electron doped FeSe
systems such as intercalated FeSe [27] and 1 monolayer FeSe grown on SrTiO3
[28].

Figure 2: (Color online). (a)-(e) Fermi surface map of Fe1-xCuxSe for
different Cu concentrations. x = 0 data was taken at 100 K in order to avoid
change in the band structure from the nematic phase transition, while others
were taken at 20 K. All the data were taken by the same ARPES spectrometer
with identical experimental setting parameters. Figure 3: (Color online). (a)
High symmetry cut and (b) schematic electronic structure of FeSe near the
$\Gamma$ point. (c)-(f) High symmetry cut and (g)-(j) its second derivative
data near the $\Gamma$ point for various x. Energy distribution curve (EDC) of
(k) raw and (l) second derivative data. Red, green and blue arrows indicate
the positions of the $\alpha$, $\beta$ and $\gamma$ bands, respectively. (m)
Band positions as function of doping extracted from the EDCs plotted in panels
(k) and (l).

In order to illustrate the evolution of electronic structure in details, we
compare band dispersions along high symmetry directions. Figures 3(a) and (b)
show a high symmetry cut of FeSe near the $\Gamma$ point and its schematic
electronic structure, respectively. Similar to previous FeSe studies [29, 30,
31], so-called $\alpha$ and $\beta$ hole bands with $d_{xz/yz}$ orbital
characters and $\gamma$ band with $d_{xy}$ orbital are well defined in our
measurement. The effect of Cu doping can be resolved in this high symmetry
spectra (Figures 3(c)-3(f)) and their corresponding second derivative spectra
(Figures 3(g)-(j)). They clearly show that hole bands shift downward with
doping. Especially, $\alpha$ hole band locates completely below $E_{F}$ from x
$>$ 0.05, leaving strong reduction of spectral weight near the Fermi-level
(Figure 3(k)), as a consequence of band gap opening. Such dramatic change of
Fermi surface topology and spectral weight near $E_{F}$ can partially account
for the insulating behavior appearing at doping x $>$ 0.05.Note that the weak
residual hole pocket feature near the $\Gamma$ point in the Fermi surface maps
of x = 0.05 and 0.10 in Figure 2 could come from broadening of dispersion tail
possibly due to impurity scattering. In Figures 3(k) and (l), we plot energy
distribution curves (EDC) of the raw data and their second derivative at the
$k=0.0{\text{\AA}}^{-1}$, respectively. Red, green and blue arrows indicate
the topmost energy positions of the $\alpha$, $\beta$ and $\gamma$ bands,
respectively. A summary on the band positions is shown in Figure 3(m), which
clearly shows downward shift of hole bands through Cu doping. It needs to be
emphasized that the doping effect seems more severe for $\alpha$ and $\beta$
bands compared to the $\gamma$ band. In other words, $d_{xz/yz}$ orbitals are
more affected by doping than $d_{xy}$ orbital, implying an orbital selective
doping dependence near the $\Gamma$ region.

Figure 4: (Color online). (a)-(e) High symmetry spectra and (g)-(j) the
corresponding second derivatives counterpart near the M2 point. (k) Doping
dependent angle integrated spectra near the M2 point. Momentum integration
window is ($-0.4{\text{\AA}}^{-1}$, $0.4{\text{\AA}}^{-1}$). (l) Effective
mass of the electron band as a function of doping. $m^{*}$ is the mass of free
electron.

Furthermore, the doping dependence of the band structures around the M point
had also been examined. Figures 4(a)-(e) show high symmetry spectra data near
the M2 point measured at 20 K. As indicated in Figure 4(a), the splitting of
the hole band (two red dash lines), as well known spectroscopic signature for
the nematic phase, is not observed in Figures 4(b)-(e). This observation
supports that nematic phase is suppressed upon Cu doping as mentioned above.
For a better visualization, we also plotted the second derivative spectra in
Figures 4(f)-(j). It is seen in the data that the electron pocket size
increases with Cu doping. Here, only one electron band is resolved, because
the low visibility of the $d_{xy}$ electron due to the matrix element effect
[29, 32, 33]. Figure 4(k) shows doping dependent angle-integrated spectra over
a momentum window of ($-0.4{\text{\AA}}^{-1}$, $0.4{\text{\AA}}^{-1}$). The
spectral weight near $E_{F}$ gradually decreases as Cu is doped. In addition,
eﬀective mass of electron band is determined by a parabolic ﬁt of the
dispersion as shown in Figures 4(b)-(e). The effective mass almost doubly
enhanced by Cu doping as shown in Figure 4(l). As mentioned above, the
electron Fermi surface pocket near the M point looks similar to that of
heavily electron doped FeSe systems which also show an enhanced effective mass
compared to the pristine FeSe. Especially, the eﬀective mass for x = 0.10 is
similar to that of 1 monolayer FeSe grown on SrTiO3 [34]. Although Cu doped
FeSe (x = 0.10) has a similar Fermi surface topology and correlation to 1
monolayer FeSe grown on SrTiO3, the reason for its lack of superconductivity
can be deduced from the impurity effect from Cu doping.

Nevertheless, our experimental observations suggest that the spectra weigh
loss near $E_{F}$ around both $\Gamma$ and M points can account for the MIT in
Cu doped FeSe, instead of simple impurity effect. Moreover, the enhanced
correlation with Cu doping (larger effective mass) should also be considered
as observed for the electron band at the M point. Emergence of a Mott
insulating state in heavily Cu doped NaFeAs [19, 35] may also point to the
importance of electron correlation in Cu doped FeSe system. Assuming that each
Cu substitution results in doping of an electron, MIT appears when the total
electron number of Fe 3d orbital is about 6.15 [36]. Therefore, the standard
Mott transition theory cannot account for insulating state of this system
because the electron filling number is not an integer [37]. Additional
electronic/magnetic order or other factors from electron correlation should
trigger the transition to an insulating phase.

Moreover, previous studies reported that magnetic fluctuation exists from x =
0.04 [17]. On the other hand, our data show that the hole band near $\Gamma$
point is located below $E_{F}$ for a similar doping level. Thus, the origin of
this magnetic fluctuation may not be solely explained within a Fermi surface
nesting picture which is mostly believed to be responsible for the magnetism
in IBS [38]. Instead, changed Fermi surface topology and enhanced correlation
should be certainly included to form theoretical model to explain the
magnetism in Fe1-xCuxSe. Concrete understanding requires further
theoretical/experimental studies.

##  IV IV. CONCLUSION

We report results of a detailed electronic structure study on Fe1-xCuxSe. We
observe that the size of the hole (electron) pocket near the $\Gamma$ (M)
point decreases (increases) with doping. With the increase of doping, we find
that the hole band at the $\Gamma$ point is located below $E_{F}$, which
changes the overall Fermi surface topology. In addition, it is found that the
effective mass of the electron band at the M point increases upon Cu doping.
Our observation explains, in terms of electronic structure, why the Fe1-xCuxSe
shows a MIT behavior and may provide clues to the mechanism of the magnetic
fluctuation.

##  V Acknowledgment

Authors would like to thank Prof. J. P. Hu, S. Ma, M. S. Kim, J. Y. Kwon and
W. S. Kyung for helpful discussions. This work is supported by IBS-R009-G2
through the IBS Center for Correlated Electron Systems, and also supported by
the National Natural Science Foundation of China (Grant Nos. 11834016,
11888101, 12061131005), the National Key Research and Development Projects of
China (Grant No. 2017YFA0303003), the Key Research Program and Strategic
Priority Research Program of Frontier Sciences of the Chinese Academy of
Sciences (Grant Nos. QYZDY-SSW-SLH001, XDB33010200 and XDB25000000).

## References

  * [1] J. Paglione and R. L. Greene, Nat. Phys. 6, 645 (2010).
  * [2] M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
  * [3] B. G. Jang, G. Han, I. Park, D. Kim, Y. Y. Koh, Y. Kim, W. Kyung, H. -D. Kim, C. -M. Cheng, K. -D. Tsuei, K. D. Lee, N. Hur, J. H. Shim, C. Kim, and G. Kotliar, Nat. Commun. 12, 1208 (2021).
  * [4] K. E. Wagner, E. Morosan, Y. S. Hor, J. Tao, Y. Zhu, T. Sanders, T. M. McQueen, H. W. Zandbergen, A. J. Williams, D. V. West, and R. J. Cava, Phys. Rev. B. 78, 104520 (2008).
  * [5] H. Miao, G. Fabbris, R. J. Koch, D. G. Mazzone, C. S. Nelson, R. A. -Esteves, G. D. Gu, Y. Li, T. Yilimaz, K. Kaznatcheev, E. Vescovo, M. Oda, T. Kurosawa, N. Momono, T. Assefa, I. K. Robinson, E. S. Bozin, J. M. Tranquada, P. D. Johnson, and M. P. M. Dean, npj Quantum Materials 6, 1 (2021).
  * [6] X. Chen, J. L. Schmehr, Z. Islam, Z. Porter, E. Zoghlin, K. Finkelstein, J. P. C. Ruff, and S. D. Wilson, Nat. Commun. 9, 103 (2018).
  * [7] J. A. Sobota, Y. He, and Z.-X. Shen, Rev. Mod. Phys. 93, 025006 (2021).
  * [8] Q. Si, R. Yu, and E. Abrahams, Nature Review Materials 1, 16017 (2016).
  * [9] X. Chen, P. Dai, D. Feng, T. Xiang, and F. -C. Zhang, National Science Review 1, 371 (2014).
  * [10] S. Rinott, K. B. Chashka, A. Ribak, E. D. L. Rienks, A. T. -Ibrahimi, P. L. Fevre, F. Bertran, M. Randeria, A. Kanigel, Sci. Adv. 3, e1602372 (2017).
  * [11] T. Hashimoto, Y. Ota, A. Tsuzuki, T. Nagashima, A. Fukushima, S. Kasahara, Y. Matsuda, K. Matsuura, Y. Mizukami, T. Shibauchi, S. Shin, K. Okazaki, Sci. Adv. 6, eabb9052 (2020).
  * [12] P. Zhang, K. Yaji, T. Hashimoto, Y. Ota, T. Kondo, K. Okazaki, Z. Wang, J. Wen, G. D. Gu, H. Ding, S. Shin, Science 360, 182 (2018).
  * [13] Z. Wang, P. Zhang, G. Xu, L. K. Zeng, H. Miao, X. Xu, T. Qian, H. Weng, P. Richard, A. V. Fedorov, H. Ding, X. Dai, and Z. Fang, Phys. Rev. B. 92, 115119 (2015).
  * [14] H. Li, M. W. Ma, S. B. Liu, F. Zhou, and X. L. Dong, Chinese Physics B 29, 127404 (2020).
  * [15] C. Gong, S. Sun, S. Wang, and H. Lei, Phys. Rev. B. 103, 174510 (2021).
  * [16] Y. Mizuguchi, F. Tomioka, S. Tsuda, T. Yamaguchi, and Y. Takano, J. Phys. Soc. Jpn. 78, 074712 (2009).
  * [17] S. I. Shylin, V. Ksenofontov, P. G. Naumov, S. A. Medvedev, and C. Felser, Journal of Superconductivity and Novel Magnetism 31, 763-769 (2018).
  * [18] L. M. Schoop, S. A. Medvedev, V. Ksenofontov, A. Williams, T. Palasyuk, I. A. Troyan, J. Schmitt, F. Casper, C. Wang, M. Eremets, R. J. Cava, and C. Felser, Phys. Rev. B. 84, 174505 (2011).
  * [19] Y. Song, Z. Yamani, C. Cao, Y. Li, C. Zhang, J. S. Chen, Q. Huang, H. Wu, J. Tao, Y. Zhu, W. Tian, S. Chi, H. Cao, Y. -B. Huang, M. Dantz, T. Schmitt, R. Yu, A. H. Nevidomskyy, E. Morosan, Q. Si, and P. Dai, Nat. Commun. 7, 13879 (2016).
  * [20] C. Ye, W. Ruan, P. Cai, X. Li, A. Wang, X. Chen, and Y. Wang, Phys. Rev. X. 5, 021013 (2015).
  * [21] S. Margadonna, Y. Takabayashi, M. T. McDonald, K. Kasperkiewicz, Y. Mizuguchi, Y. Takano, A. N. Fitch, E. Suard, and K. Prassides, Chem. Commun. 2008, 5607 (2008).
  * [22] Y. J. Yan, P. Cheng, J. J. Ying, X. G. Luo, F. Chen, H. Y. Zou, A. F. Wang, G. J. Ye, Z. J. Xiang, J. Q. Ma, and X. H. Chen, Phys. Rev. B. 87, 075105 (2013).
  * [23] V. K. Anand, P. Kanchana Perera, Abhishek Pandey, R. J. Goetsch, A. Kreyssig, and D. C. Johnston, Phys. Rev. B. 85, 214523 (2012).
  * [24] S. T. Cui, S. Kong, S. L. Ju, P. Wu, A. F. Wang, X. G. Luo, X. H. Chen, G. B. Zhang, and Z. Sun, Phys. Rev. B. 88, 245112 (2013).
  * [25] S. T. Cui, S. Y. Zhu, A. F. Wang, S. Kong, S. L. Ju, X. G. Luo, X. H. Chen, G. B. Zhang, and Z. Sun, Phys. Rev. B. 86, 155143 (2012).
  * [26] S. Ideta, T. Yoshida, I. Nishi, A. Fujimori, Y. Kotani, K. Ono, Y. Nakashima, S. Yamaichi, T. Sasagawa, M. Nakajima, K. Kihou, Y. Tomioka, C. H. Lee, A. Iyo, H. Eisaki, T. Ito, S. Uchida, and R. Arita Phys. Rev. Lett. 110, 107007 (2013).
  * [27] Y. Zhang, L. X. Yang, M. Xu, Z. R. Ye, F. Chen, C. He, H. C. Xu, J. Jiang, B. P. Xie, J. J. Ying, X. F. Wang, X. H. Chen, J. P. Hu, M. Matsunami, S. Kimura, and D. L. Feng, Nat. Materials 10, 273-277 (2011).
  * [28] S. Tan, Y. Zhang, M. Xia, Z. Ye, F. Chen, X. Xie, R. Peng, D. Xu, Q. Fan, H. Xu, J. Jiang, T. Zhang, X. Lai, T. Xiang, J. Hu, B. Xie, and D. Feng, Nat. Materials 12, 634-640 (2013).
  * [29] S. S. Huh, J. J. Seo, B. S. Kim, S. H. Cho, J. K. Jung, S. Kim, C. I. Kwon, J. S. Kim, Y. Y. Koh, W. S. Kyung, J. D. Denlinger, Y. H. Kim, B. N. Chae, N. D. Kim, Y. K. Kim, and C. Kim, Commun. Phys. 3, 52 (2020).
  * [30] D. Liu, C. Li, J. Huang, B. Lei, L. Wang, X. Wu, B. Shen, Q. Gao, Y. Zhang, X. Liu, Y. Hu, Y. Xu, A. Liang, J. Liu, P. Ai, L. Zhao, S. He, L. Yu, G. Liu, Y. Mao, X. Dong, X. Jia, F. Zhang, S. Zhang, F. Yang, Z. Wang, Q. Peng, Y. Shi, J. Hu, T. Xiang, X. Chen, Z. Xu, C. Chen, and X. J. Zhou, Phys. Rev. X. 8, 031033 (2018).
  * [31] T. Shimojima, Y. Suzuki, A. Nakamura, N. Mitsuishi, S. Kasahara, T. Shibauchi, Y. Matsuda, Y. Ishida, S. Shin, and K. Ishizaka, Nat. Commun. 10, 1946 (2019).
  * [32] M. D. Watson, T. K. Kim, L. C. Rhodes, M. Eschrig, M. Hoesch, A. A. Haghighirad, and A. I. Coldea, Phys. Rev. B. 94, 201107 (2016).
  * [33] C. Cai, T. T. Han, Z. G. Wang, L. Chen, Y. D. Wang, Z. M. Xin, M. W. Ma, Y. Li, and Y. Zhang, Phys. Rev. B. 101, 180501(R) (2020).
  * [34] D. Liu, W. Zhang, D. Mou, J. He, Y. -B. Ou, Q. -Y. Wang, Z. Li, L. Wang, L. Zhao, S. He, Y. Peng, X. Liu, C. Chen, L. Yu, G. Liu, X. Dong, J. Zhang, C. Chen, Z. Xu, J. Hu, X. Chen, X. Ma, Q. Xue, and X.J. Zhou, Nat. Commun. 3, 931 (2012).
  * [35] C. E. Matt, N. Xu, B. Lv, J. Ma, F. Bisti, J. Park, T. Shang, C. Cao, Y. Song, A. H. Nevidomskyy, P. Dai, L. Patthey, N. C. Plumb, M. Radovic, J. Mesot, and M. Shi, Phys. Rev. Lett. 117, 097001 (2016).
  * [36] R. Yu, H. Hu, E. M. Nica, J. -X. Zhu, and Q. Si, Frontiers in Physics 9, 578347 (2021).
  * [37] J. Kwon, M. Kim, D. Song, Y. Yoshida, J. Denlinger, W. Kyung, and C. Kim, Phys. Rev. Lett. 123, 106401 (2019).
  * [38] I. I. Mazin, D. J. Singh, M. D. Johannes, and M. H. Du, Phys. Rev. Lett. 101, 057003 (2008).


# The effect of selection – a tale of cluster mass measurement bias induced by
correlation and projection

Yuanyuan Zhang1,2, James Annis2
1 Mitchell Institute for Fundamental Physics and Astronomy and Department of
Physics and Astronomy, Texas A&M University, College Station, TX 77843-4242
2 Fermilab, Cosmic Physics Center, Kirk& Pine Road, Batavia, IL, 60510
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Cosmology analyses using galaxy clusters by the Dark Energy Survey have
recently uncovered an issue of previously unknown selection effect affecting
weak lensing mass estimates. In this letter, we use the Illustris-TNG
simulation to demonstrate that selecting on galaxy counts induces a selection
effect because of projection and correlation between different observables. We
compute the weak-lensing-like projected mass estimations of dark matter halos,
and examine their projected subhalo counts. In the 2-D projected space, halos
that are measured as more massive than truth have higher subhalo counts. Thus,
projection along the line of sight creates cluster observables that are
correlated with cluster mass measurement deviations, which in turn creates a
mass measurement bias when the clusters are selected by this correlated
observable. We demonstrate that the bias is predicted in a forward model using
the observable-mass measurement correlation.

###### keywords:

(cosmology:) large-scale structure of Universe – galaxies: clusters: general

††pubyear: 2015††pagerange: The effect of selection – a tale of cluster mass
measurement bias induced by correlation and projection–The effect of selection
– a tale of cluster mass measurement bias induced by correlation and
projection

##  1 Introduction

Galaxy clusters have long been utilized as sensitive probes of cosmology in
optical wide-field surveys such as the Sloan Digital Sky Survey (Rozo et al.,
2010; Zu et al., 2014; Costanzi et al., 2019b), the Dark Energy Survey (Abbott
et al., 2020) and also planned for the upcoming Legacy Survey of Space and
time at the Rubin Observatory (The LSST Dark Energy Science Collaboration et
al., 2018). These optical survey programs study galaxy clusters by their
optical observables (e.g., Rykoff et al., 2014; Palmese et al., 2020) and weak
lensing mass signals (Simet et al., 2017; McClintock et al., 2019; Pereira et
al., 2020), often to a lower mass threshold than X-ray and CMB experiments
(e.g., Vikhlinin et al., 2009; Mantz et al., 2015; Planck Collaboration et
al., 2016; Bocquet et al., 2019). However recent discoveries from the Dark
Energy Survey in Abbott et al. (2020) point to additional systematic effects
that plague the accuracy of cluster cosmological constraints, which stem from
previously undetected selection effects biasing cluster weak lensing mass
measurements. Ongoing and future cosmic surveys will need to address this
challenge in order to achieve the full potential of galaxy cluster cosmology
analyses.

Specifically, the analysis in Abbott et al. (2020) shows that the weak-lensing
measured masses of galaxy clusters appear to have deviated from their model
values significantly, which, in the end, affect cosmological parameters
derived from modeling the cluster mass distribution and abundance. Abbott et
al. (2020) and further Sunayama et al. (2020) find that the bias of cluster
weak-lensing measurements likely originate from cluster selection. Their
analyses based on simulations have demonstrated that galaxy clusters selected
by their richness observable, defined as a weighted number count of red
sequence cluster galaxies, is a biased population when compared to a unbiased
population selected by their unbiased truth quantities. Sunayama et al. (2020)
has further determined that the projection of cosmic structures along the line
of sight, as well as cluster orientation and shapes (Dietrich et al., 2014;
Osato et al., 2018, and Z. Zhang in prep.) may be causing this selection bias.

These effects may cause additional correlation between the cluster’s
observable and its weak-lensing mass measurement, and thus we can consider
alleviating their effects through modeling the correlations ( as in, e.g.,
Zhang et al., 2019; Grandis et al., 2020, 2021b). In this paper, we
demonstrate the potential effectiveness of this correlation approach by
examining the projected observables of cluster-sized dark matter halos and
their correlations with weak-lensing like mass measurements in simulation.
Through doing so, we validate previous conclusions about the "projection"
origin of the previously over-looked selection effect.

##  2 Simulation Data

Figure 1: Projected mass surface density profiles of the massive dark matter
halos studied in this paper. The color coding indicate their weak-lensing-like
mass measurement deviations shown in Figure 2. The solid and dashed black
lines show the mean and the 1 $\sigma$ uncertainty range of the distributions.

In this work, we use products from the IllustrisTNG simulation suite, in
particular, from the IllustrisTNG 300-1 simulation (Pillepich et al., 2018;
Nelson et al., 2018; Springel et al., 2018; Naiman et al., 2018; Marinacci et
al., 2018; Nelson et al., 2019), the largest hydro-dynamic simulation box,
with high-resolution baryonic processes incorporated. We analyze dark matter
halos in the $79th$ snapshot, corresponding to redshift 0.27, which is the
intermediate redshift of the lowest redshift bin used in DES cluster cosmology
analysis. We have also analyzed the simulation snapshots that correspond to
redshift 0.42 and 0.58, and find that similar effects reported in this letter
are also present.

We analyze halos with $M_{\mathrm{200m}}$ above $10^{13.75}M_{\odot}/h$ 111Our
analysis of lower mass halos, especially those below $10^{13.5}M_{\odot}/h$,
shows that the adopted halo models in this paper do not work well, which may
be worthy of a separate study itself.. For each of the halo, we derive their
density profiles using the dark matter, gaseous and stellar particles in the
simulation snapshot. These particles are projected onto the $X-Y$ plane of the
simulation with a projection depth along the $Z$ direction of 120
$\mathrm{cMpc}/h$ 222The letter $c$ indicates comoving distance as a notation
adopted in IllustrisTNG simulation.. Given that the the simulation box size is
205 $\mathrm{cMpc}/h$ on each side, to maximize the range of projection depth,
we only use the particles that are located on one side of the halo’s $Z$-axis
up to 120 $\mathrm{cMpc}/h$ from the halo center. The projected density
profile is derived according to the dark matter, stellar and gaseous
particles’ projected distance to the halo center on the $X-Y$ plane, and then
multiplied by a factor of 2 assuming the halo to be symmetrical on the
$Z$-axis to recover its total density profile. A projected background density
is estimated as the averaged density of the dark matter, stellar and gaseous
particles projected onto the $X-Y$ plane, with a thickness of 240 cMpc/h along
the $Z$-axis. To avoid running into the simulation boundaries, we exclude
halos that are within 20 cMpc/h of the simulation box, or 20 $\mathrm{cMpc/h}$
within the the $Z$-axis middle plane ($z=102.50\mathrm{~{}cMpc/h}$). In the
end, we analyze a total of 254 dark matter halos.

We have investigated the derived halo matter profiles are sensitivities to (1)
baryonic effects (2) simulation resolution and (3) projection depth. In the
IllustrisTNG dark-matter only simulation, TNG-Dark-300-1, the density profiles
of halos above $10^{14}M_{\odot}/h$ show a relative difference of up to $\sim
10\%$ in the central 200 $\mathrm{cKpc}/h$ region when compared to halo
profiles from the TNG 300-1 hydrodynamic simulation. Similarly, when we
compare the halo density profiles of our fiducial measurements from the
TNG300-1 hydrodynamic simulation to profiles from the lower resolution TNG
300-3, TNG 300-2 simulations, we find differences up to $\sim 5\%$ within the
central 100 $\mathrm{cKpc}/h$ regions. Given the accuracies of previous
cluster mass calibration studies (Becker & Kravtsov, 2011; Grandis et al.,
2021a), we would like profiles insensitive to baryonic effects or simulation
resolution at the 10% or ideally 5% level. Thus we exclude the central 200
$\mathrm{cKpc}/h$ regions of the halo profiles during the further analysis.

Our investigation of projection depth sets the outer radius. We compute the
projected halo density profiles and corresponding projected background
densities for projection length, 75, 100 and 135 $\mathrm{cMpc}/h$ from the
TNG300-1 and TNG300-3 simulations. In the central 2 $\mathrm{cMpc}/h$ regions,
the halos’ profiles agree within $\sim 10\%$ for projection depth
$>100\mathrm{cMpc}/h$. Outside the 2 $\mathrm{cMpc}/h$ radial range, using
different projection depths causes noticeable fluctuations, and the amplitude
of fluctuations decreases as halo masses increase. It is unclear if the halo
density profiles agree to within $10\%$ relative difference. Therefore, in
this work, we only analyze the 0.2 to 2 $\mathrm{cMpc}/h$ radial range of the
halo density profiles.

The derived halo density profiles are shown in Fig. 1. We also derive the
covariance matrix as the halo-to-halo variation after normalizing the halo
profiles by their radii, the $R_{200m}$ value of each halo (see studies in Wu
et al., 2019). This covariance matrix is later used during the derivation of
halo masses from their projected density profiles.

##  3 Halo Mass Observables

Figure 2: The recovered halo masses, and their 1 $\sigma$ uncertainties
(standard deviation from the MCMC posterior sampling) by fitting the projected
halo mass density profiles with analytical models. Note that there is a
fraction of dark matter halos below $M_{200m}$ of
$10^{14}\mathrm{M_{\odot}}/h$ that we are unable to reliably recover their
masses as the fitting values of concentration and mass has failed to converge,
indicated by the data points with large error bars (light-colored data
points). For the rest of the halos (dark-colored data points), the analytical
models recover the true masses with an average of 1.4% bias.

We first derive halo masses in the projected $X-Y$plane. We fit the projected
halo profiles to a model corresponding to both a halo’s gravitationally-bound
matter distribution and the contribution from large scale structures, using
the following form:

$\begin{split}\Sigma(r|M_{200m},c,z)&=\mathrm{max}(\Sigma_{\mathrm{NFW}}(r|M_{200m},c),\Sigma_{\mathrm{2-halo}}(r|M_{200m},z)),\\\
\Sigma_{\mathrm{2-halo}}(r|M_{200m},z)&=b(M_{200m},z)\Sigma_{\mathrm{nl}}(r|M_{200m},z)\\\
\Sigma_{\mathrm{nl}}(r|M_{200m},z)&=\int^{+\infty}_{-\infty}\mathrm{d}\chi\rho_{m}\xi_{\mathrm{nl}}(\sqrt{r^{2}+\chi^{2}}).\end{split}$
(1)

Here we combine the projected NFW halo profiles, the so-called "one-halo"
term, and the nonlinear matter correlation function weighted by halo bias
$b(M_{200m},z)$, the "two-halo" term. In this equation, $\xi_{\mathrm{nl}}$ is
the non-linear matter correlation function, which is the 3D Fourier
transformation of the non-linear power spectrum. These fitting models are
inspired by those adopted in Simet et al. (2017); McClintock et al. (2019) and
tested to give similar mass estimations within $\sim 1\%$ but with improved
speed. During the fitting procedure, the model masses and concentrations are
treated as varying parameters sampled through Markov Chain Monte Carlo (MCMC),
which is implemented using EMCEE (Foreman-Mackey et al., 2013) with a
likelihood constructed from the $\chi^{2}$ value between the model and each
halos’s measurements. The constrained values of the halo’s masses are shown in
Fig. 2 and compared to each halo’s truth values.

There are some dark matter halos below $\mathrm{log}M_{200m}$ of
$10^{14}\mathrm{M_{\odot}}/h$ for which we are unable to recover their masses
as the fittings of concentration and mass fail to converge to reasonable
values. In Fig. 2, those halos have large mass uncertainties, exceeding 1.0
dex. We were, however, able to recover the masses of those halos using their
3D density profiles, and those halos also tend to have a projected density
profile significantly below the average, as shown in Figure 1. This indicates
that they may have under-dense environments and thus their profiles are not
represented by the analytical models in the 3D to 2D projection process. For
the rest of the halos, the analytical models recover the true masses with an
average bias of 1.4% dex.

Figure 3: Relations between subhalo counts, $N_{\mathrm{3D}}$(left),
$N_{\mathrm{2D}}$ (right) and halo’s masses, color coded by the halo’s mass
measurement deviations (red, green and blue). $N_{\mathrm{2d}}$ displays signs
of correlation with the halo’s mass deviations – halos that have positive mass
deviations (blue squares) tend to have higher subhalo counts than the average,
and vice versa for the halos that have negative mass deviations (red rounds).
$N_{\mathrm{3D}}$ displays no sign of the correlation. The dashed lines in the
upper panels indicate the fitted model predictions for the halos with mass
deviations in different value ranges, and the lower panels show the residuals
between the measured subhalo counts and their mean model values, with the
dotted lines showing the modelled scatters.

In this section, we compute the number of subhalos in dark matter halos as a
halo mass proxy. In optical studies, the cluster selection process often
relies on a galaxy over-density above a given threshold. In this exercise, we
aim to reproduce this selection by counting the number of massive subhalos
with masses above $5\times 10^{9}\mathrm{M_{\odot}}/h$, within a radial
aperture around the halo centers. We consider two different kinds of apertures
based on the 1 $\mathrm{Mpc}$ physical distance radius (close to the average
$R_{\mathrm{200m}}$ analyzed of the halo samples, which is 0.87 cMpc/h):

  * •

A 3D-radial aperture of 1 $\mathrm{cMpc}/h$. With this aperture, any massive
subhalos that have a 3-D distance less than 1 $\mathrm{cMpc}/h$ to the halo
center are counted. The number of subhalos satisfying this criteria is
designated $N_{\mathrm{3D}}$.

  * •

A 2D-radial aperture of 1 $\mathrm{cMpc}/h$. With this aperture, massive
subhalos that have a 2-D distance (in the $X-Y$ plane of the simulation) less
than 1 $\mathrm{cMpc}/h$ to the halo center are counted. Only subhalos within
a given projection length selected in the same way of the mass profile
projection in Sect. 2 are considered. The number of subhalos is designated
$N_{\mathrm{2D}}$.

In Fig. 3, we demonstrate that the halo’s mass measurement deviations are
correlated with subhalo number counts, when controlling the halo’s truth mass.
In this figure, the subhalo number counts of $N_{\mathrm{3D}}$ and
$N_{\mathrm{2D}}$ are shown against the halo’s truth mass, showing a clear
relation between the two: more massive halos host more subhalos. More
interestingly, the relation between the subhalo counts and halo mass is
further correlated with the halo’s mass measurement deviations (derived in the
previous sub-section in the case of $N_{\mathrm{2D}}$). Halos that are
measured as more massive than their truth masses have a higher
$N_{\mathrm{2D}}$ values (blue data points), while halos that are measured to
be less massive than their truth masses have lower $N_{\mathrm{2D}}$ values
(red data points).

To further quantify this correlation, we paramterize the relation between
subhalo counts and the halo’s masses using the following statistical model for
counts as a function of mass and mass deviation:

$\begin{split}\mathrm{log}N&=N(\mu(M_{\mathrm{true}},M_{\mathrm{fit}}),\sigma^{2}(M_{\mathrm{true}},M_{\mathrm{fit}}))\\\
\Delta&=\mathrm{log}M_{\mathrm{fit}}-\mathrm{log}M_{\mathrm{true}}\\\
\mu(M_{\mathrm{true}},M_{\mathrm{fit}})&=a\times(\mathrm{log}M_{\mathrm{true}}-14.0)+b+(\alpha+\beta\times\Delta)\times\Delta\\\
\sigma(M_{\mathrm{true}},M_{\mathrm{fit}})&=\sigma_{0}+q\times(\mathrm{log}M_{\mathrm{true}}-15.5)\end{split}$
(2)

This model has a linear relation between the subhalo counts and halo masses at
the log scale as adopted in some cluster cosmology studies (e.g., Rozo et al.,
2010; Zu et al., 2014), although newer studies suggest a more complicated
functional forms to be more accurate (Costanzi et al., 2019a, b). The subhalo
counts deviations from the mean are described by a Gaussian scatter, and the
Gaussian scatter depends on the halo mass in a linear form as adopted in
Murata et al. (2018). However in order to quantify the additional correlation
between $N_{\mathrm{2D}}$ and the halo mass measurement deviation, we further
incorporate first-order and second-order dependencies on the mass measure
deviations
$\Delta=\mathrm{log}M_{\mathrm{fit}}-\mathrm{log}M_{\mathrm{true}}$, in the
mean linear relation, controlled by the values of parameters $\alpha$ and
$\beta$. We have attempted incorporating these dependencies in the
$N_{\mathrm{2D}}$ Gaussian scatters, but do not find significant dependencies
and therefore do not adopt them here. The free parameters in the relations,
$a$, $b$, $\alpha$, $\beta$, $\sigma_{0}$ and $q$ are constrained using Markov
Chain Monte Carlo sampling with flat Bayesian priors.

Table 1: Constraints on the parameters of Eq. 2, the halo’s subhalo counts observable to mass relation, for different subhalo count definitions. Model | $a$ | $b$ | $\alpha$ | $\beta$ | $\sigma_{0}$ | q |  $\chi^{2}$ ($n=254$)
---|---|---|---|---|---|---|---
$N_{\mathrm{3D}}$ | $0.65\pm 0.04$ | $1.195\pm 0.011$ | $0.009\pm 0.033$ | $-0.003\pm 0.006$ | $0.096\pm 0.038$ | $-0.052\pm 0.026$ | 245.38
$N_{\mathrm{2D}}$ | $0.439\pm 0.037$ | $1.491\pm 0.011$ | $0.259\pm 0.031$ | $0.036\pm 0.006$ | $0.056\pm 0.032$ | $-0.069\pm 0.022$ | 245.35
$N_{\mathrm{r3D}}$ | $0.916\pm 0.035$ | $1.23\pm 0.011$ | $0.011\pm 0.032$ | $-0.002\pm 0.006$ | $0.043\pm 0.033$ | $-0.086\pm 0.022$ | 245.71
$N_{\mathrm{r2D}}$ | $0.750\pm 0.035$ | $1.56\pm 0.011$ | $0.268\pm 0.031$ | $0.038\pm 0.006$ | $0.041\pm 0.030$ | $-0.070\pm 0.020$ | 245.42
$N_{\mathrm{iterR}}$ | $0.735\pm 0.056$ | $1.53\pm 0.018$ | $0.414\pm 0.053$ | $0.062\pm 0.010$ | $0.056\pm 0.044$ | $-0.146\pm 0.030$ | 246.43

The posterior values of the constrained parameters are listed in Table 1.
Interestingly, $N_{\mathrm{2D}}$ have significant correlation with the mass
measurement deviation, given the positive values of $\alpha$ and $\beta$. On
the other hand, $N_{\mathrm{3D}}$ does not have a significant correlation with
the mass measurement deviation and the values of $\alpha$ and $\beta$ are
consistent with 0. This result is in agreement with the qualitative
observations in Fig. 3.

We further investigate this effect. In observational studies, the radial
apertures used to derive a cluster galaxy over-density observable is often
iteratively adjusted according to the cluster’s “sizes”/“masses”. Thus we
further test the correlations with the following subhalo counting definitions
that use apertures based on halo masses:

  * •

A 3D-radial aperture of $R_{\mathrm{200m}}$, derived from the halo’s
$M_{\mathrm{200m}}$. Any massive subhalos that have a 3-D distance less than
$R_{\mathrm{200m}}$ to the halo centers are counted. The number of subhalos
satisfying this criteria is designated $N_{\mathrm{r3D}}$.

  * •

A 2D-radial aperture of $R_{\mathrm{200m}}$, derived from the halo’s
$M_{\mathrm{200m}}$. Any massive subhalos that have a 2-D distance (in the X-Y
plane of the simulation) less than $R_{\mathrm{200m}}$ to the halo centers are
counted. The number of subhalos satisfying this criteria is designated
$N_{\mathrm{r2D}}$.

  * •

An iterative 2D-radial aperture $R_{\mathrm{iter}}$, which is adjusted
according to the derived subhalo counts $N_{\mathrm{iterR}}$ until they
satisfy the relation
$N_{\mathrm{iterR}}=26.514\times(R_{\mathrm{iter}}/\mathrm{[cMpc/h]})^{2.414}$.
This relation is the mean relation between $R_{\mathrm{200m}}$ and
$N_{\mathrm{r2D}}$ defined above.

The last is designed to model the usual case in optical observations that a
cluster’s true $R_{\mathrm{200m}}$ is not known, but the galaxy counting
apertures can be iteratively improved from the counts.

The parameterization of the relations between those subhalo counts, and halo
masses are also listed in Table 1. Again, the 3D-radial aperture based
observable, $N_{\mathrm{r3D}}$, shows no correlation with the halo mass
measurement deviation. On the other hand, $N_{\mathrm{iterR}}$ has the
strongest correlation with the mass measurement deviation. We expect being in
a high-density environment, or other similar factors, affect both halo mass
measurements and the apertures used to select galaxy over-densities. This
creates the greater correlation between the two quantities as seen with
$N_{\mathrm{r3D}}$.

##  4 Mass Bias of the Selected Halos

When there exist correlations between the halos’s observed mass and their
selection observable, the selected halo sample exhibits a mass measurement
bias. Assuming that the halo’s observed mass is a fractional deviation from
the halo’s truth mass, the mathematical expression of the bias can be derived
by examining the average masses of the halos selected by observable,

$\begin{split}M_{M}(N)&=\frac{1}{P(N)}\int_{-\infty}^{\infty}\mathrm{d}\Delta\int_{10^{13.75}}^{\infty}\mathrm{d}M(M10^{\Delta})P(N|M,\Delta)P(\Delta|M)P(M)\\\
M_{T}(N)&=\frac{1}{P(N)}\int_{-\infty}^{\infty}\mathrm{d}\Delta\int_{10^{13.75}}^{\infty}\mathrm{d}M(M)P(N|M,\Delta)P(\Delta|M)P(M)\\\
\text{Bias}&=\mathrm{log}\frac{M_{M}(N)}{M_{T}(N)}\\\ \end{split}$ (3)

In those equations, $M_{M}$ represents the average measured mass of the halos
selected by the observable $N$, while $M_{T}$ represents their truth average
mass. In computing those biases, $P(N|M,\Delta)$ is the halo’s mass-observable
relation, incorporating $N$’s correlation to halo masses and measurement
deviations studied in Section 3, while $P(M)$ is the halo’s true mass
distribution, derived from the theoretical halo mass function. $P(\Delta|M)$
models the distribution of individual halo mass measurement deviations without
selection, which has been studied previously in literature in the context of
"mass calibration" for cluster weak lensing studies (Becker & Kravtsov, 2011).
Here, we model it as a Log-Normal distribution, dependent on the halo’s truth
mass as described in Sect. 3. We impose a requirement that
$\int_{-\infty}^{\infty}\mathrm{d}\Delta\mathrm{10}^{\Delta}P(\Delta|M)=1$, by
adjusting the mean value of $\Delta$ based on the measured scatter of $\Delta$
so that when there is no correlation between $N$ and $\Delta$, $P(N|M,\Delta)$
reduces to $P(N|M)$, and $M_{M}(N)=M_{N}(N)$. In the situation that $N$
correlates to $\Delta$ in $P(N|M,\Delta)$, then $M_{M}(N)$ and $M_{N}(N)$ are
no longer necessarily equivalent.

We quantify the mass measurement bias of halos selected by $N$ as
$\mathrm{log}\frac{M_{M}(N)}{M_{T}(N)}$. Using quantities derived in previous
sections and a Tinker mass function (Tinker et al., 2008) as $P(M)$, we show
the bias prediction $\mathrm{log}\frac{M_{M}}{M_{T}}$ in Fig. 4. When the
selection observable $N$ has no significant correlation with the halo mass
measurement deviation, as is the case with $N_{\mathrm{3D}}$ and
$N_{\mathrm{r3D}}$, the derived biases are consistent with 0 (shaded red and
yellow bands respectively). When the selection observable are correlated, as
are the cases with $N_{\mathrm{2D}}$, $N_{\mathrm{r2D}}$, and
$N_{\mathrm{iterR}}$, the theoretically derived biases (shaded green, blue and
magenta bands) are no longer consistent with 0. Note that because of the
artificially-imposed mass selection cut at $10^{13.75}\mathrm{M_{\odot}}/h$,
the mass measurement bias turns negative at the lower mass end in both
measurement and theoretical calculation because of the lack of up-scattered
halos from lower mass ranges. This trend would disappear if we removed the
halo mass selection threshold (dashed lines). In general, the more correlated
observable yield more prominent mass biases when selected upon, and the biases
also vary with halo masses.

Figure 4: Bias in the measurement of the clusters’ average masses while the
clusters are selected by different mass observables. Clusters selected by the
projected subhalo counts ($N_{\mathrm{2D}}$, dash-dotted green line), which
has strong observable-mass deviation correlation, show a high level of
measurement bias, while the clusters selected by 3D subhalo counts
($N_{\mathrm{3D}}$, dashed red line) or true masses
($\mathrm{log}M_{\mathrm{200m}}$, solid black line) show no significant
biases.

Do those expectation match what we would observe in observations? We further
measure the mass measurement biases of the halos selected according to the
different subhalo counts. We group the halos according to their subhalo counts
and measure their average masses, a technique often employed by cluster
lensing analyses. The average masses are measured by fitting Equations 1 to
their averaged radial profiles with a single set of mass and concentration
values. This method does assume that the averaged halo profiles are well
described by the model of a single halo, but this has been tested to be
accurate at the 5% level (Melchior et al., 2017; McClintock et al., 2019) for
recovered halo masses.

Fig. 4 shows the results of the recovered average halo mass selected by the
different observables, in terms of their measurement biases vs their truth
average mass. The halos selected by $N_{\mathrm{2D}}$, $N_{\mathrm{r2D}}$, and
$N_{\mathrm{iterR}}$ show biases in their averaged mass measurements. This is
in contrast to the halos that are grouped according to their truth masses
$\mathrm{log}M_{\mathrm{200m}}$ and $N_{3D}$, $N_{\mathrm{r3D}}$, which show a
negligible level of biases. Those biases qualitatively match the quantitative
expectations from the parameterized models, and in the high mass end, also
match the effect of selection biases discovered in Abbott et al. (2020);
Sunayama et al. (2020) at the $\sim 10\%$ level.

In observational applications, the exact forms of those correlations will not
need to be precisely known. In cluster cosmological analysis, the cluster
mass-observable relations are modeled as nuisance parameters together with
cosmological parameters, which can be expanded to include correlations and
simultaneously constrained with the cluster cosmological parameters. Our
analysis indicate that including the correlation parameters has the potential
to minimize the selection bias plaguing the precision of cluster cosmology
analysis as discussed in Abbott et al. (2020). An alternative for tackling the
selection bias is to rely on cluster observables that are either not or are
less affected by projection effect to select cluster samples, as demonstrated
by the mass bias results based on $N_{3D}$, $N_{\mathrm{r3D}}$.

##  5 Summary

In this letter, we use the TNG simulation to explore the effects of cluster
selection on weak-lensing-like mass measurements. We demonstrate that when
dark matter halo are selected by their projected galaxy-like observables,
these observables become correlated with the halo weak-lensing-like mass
measurement deviations because weak lensing measurements are made in projected
space. As a result, the dark matter halos selected and ranked upon those
observables display a collective mass measurement biases. This bias is also
predicted in models that account for the observable-mass correlation and the
mass measurement scatter.

Our study differs from previous ones, for example Sunayama et al. (2020), as
we prioritize simplicity and a direct modeling solution over fidelity to the
cluster finding process. In our assumption that cluster galaxies are modeled
by dark matter subhalos in the massive halo, we preserve the shape and
orientation of each halo. Our procedure illuminates the effect of projection
along the line of sight on both the subhalo counts and measured halo mass.

This analysis provides insight into the selection effect that plagues recent
cluster cosmology analyses. Although the efficacy of the forward-modeling
approach in our work still need to be tested with observational data sets and
a set of more realistic simulations with light-cone realization and high-
resolution mass measurement observable, our analysis demonstrates that the
selection effect can be predicted with assumptions about the the correlation
between the selection observable and the cluster weak lensing mass measurement
deviations. It would be prudent to further understand the physical origins of
those correlations, and develop a well-motivated functional form for the
correlation. We recommend developing cluster selection observables that are
less affected by projected observations on the plane of the sky, thus reducing
the level of correlations between cluster selection and weak lensing mass
measurement.

##  6 Acknowledgements

We thank the anonymous referee for the very helpful comments and suggestions.
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. The IllustrisTNG simulations were
undertaken with compute time awarded by the Gauss Centre for Supercomputing
(GCS) under GCS 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), as well as on the machines of the Max Planck Computing and Data
Facility (MPCDF) in Garching, Germany.

##  7 Data Availability

The data underlying this article were accessed from the Illustris-TNG
database. The derived data generated in this research will be shared on
reasonable request to the corresponding author.

## References

  * Abbott et al. (2020) Abbott T. M. C., et al., 2020, Phys. Rev. D, 102, 023509
  * Becker & Kravtsov (2011) Becker M. R., Kravtsov A. V., 2011, ApJ, 740, 25
  * Bocquet et al. (2019) Bocquet S., et al., 2019, ApJ, 878, 55
  * Costanzi et al. (2019a) Costanzi M., et al., 2019a, MNRAS, 482, 490
  * Costanzi et al. (2019b) Costanzi M., et al., 2019b, MNRAS, 488, 4779
  * Dietrich et al. (2014) Dietrich J. P., et al., 2014, MNRAS, 443, 1713
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Grandis et al. (2020) Grandis S., et al., 2020, MNRAS, 498, 771
  * Grandis et al. (2021a) Grandis S., Bocquet S., Mohr J. J., Klein M., Dolag K., 2021a, arXiv e-prints, p. arXiv:2103.16212
  * Grandis et al. (2021b) Grandis S., et al., 2021b, MNRAS, 504, 1253
  * Mantz et al. (2015) Mantz A. B., et al., 2015, MNRAS, 446, 2205
  * Marinacci et al. (2018) Marinacci F., et al., 2018, MNRAS, 480, 5113
  * McClintock et al. (2019) McClintock T., et al., 2019, MNRAS, 482, 1352
  * Melchior et al. (2017) Melchior P., et al., 2017, MNRAS, 469, 4899
  * Murata et al. (2018) Murata R., Nishimichi T., Takada M., Miyatake H., Shirasaki M., More S., Takahashi R., Osato K., 2018, ApJ, 854, 120
  * Naiman et al. (2018) Naiman J. P., et al., 2018, MNRAS, 477, 1206
  * Nelson et al. (2018) Nelson D., et al., 2018, MNRAS, 475, 624
  * Nelson et al. (2019) Nelson D., et al., 2019, Computational Astrophysics and Cosmology, 6, 2
  * Osato et al. (2018) Osato K., Nishimichi T., Oguri M., Takada M., Okumura T., 2018, MNRAS, 477, 2141
  * Palmese et al. (2020) Palmese A., et al., 2020, MNRAS, 493, 4591
  * Pereira et al. (2020) Pereira M. E. S., et al., 2020, MNRAS, 498, 5450
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 475, 648
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A24
  * Rozo et al. (2010) Rozo E., et al., 2010, ApJ, 708, 645
  * Rykoff et al. (2014) Rykoff E. S., et al., 2014, ApJ, 785, 104
  * Simet et al. (2017) Simet M., McClintock T., Mandelbaum R., Rozo E., Rykoff E., Sheldon E., Wechsler R. H., 2017, MNRAS, 466, 3103
  * Springel et al. (2018) Springel V., et al., 2018, MNRAS, 475, 676
  * Sunayama et al. (2020) Sunayama T., et al., 2020, MNRAS, 496, 4468
  * The LSST Dark Energy Science Collaboration et al. (2018) The LSST Dark Energy Science Collaboration et al., 2018, arXiv e-prints, p. arXiv:1809.01669
  * Tinker et al. (2008) Tinker J., Kravtsov A. V., Klypin A., Abazajian K., Warren M., Yepes G., Gottlöber S., Holz D. E., 2008, ApJ, 688, 709
  * Vikhlinin et al. (2009) Vikhlinin A., et al., 2009, ApJ, 692, 1060
  * Wu et al. (2019) Wu H.-Y., Weinberg D. H., Salcedo A. N., Wibking B. D., Zu Y., 2019, MNRAS, 490, 2606
  * Zhang et al. (2019) Zhang Y., et al., 2019, MNRAS, 487, 2578
  * Zu et al. (2014) Zu Y., Weinberg D. H., Rozo E., Sheldon E. S., Tinker J. L., Becker M. R., 2014, MNRAS, 439, 1628


PIP-DualP | 40 | 90.70 | 87.15 | 87.17 | 87.43 | 86.98 | 86.75 | 85.17 | 84.05 | 83.26 | 82.11 | 86.08
Fed-L2P | 20 | 61.00 | 59.45 | 63.23 | 64.85 | 65.10 | 65.03 | 62.39 | 60.35 | 59.19 | 57.15 | 61.77
Fed-DualP | 20 | 73.90 | 75.55 | 75.10 | 74.28 | 72.44 | 72.75 | 71.51 | 69.45 | 68.78 | 70.71 | 72.45
PIP-L2P | 20 | 80.20 | 78.15 | 77.43 | 77.83 | 76.12 | 76.85 | 74.14 | 73.64 | 71.77 | 67.18 | 75.33
PIP-DualP | 20 | 89.20 | 87.90 | 87.63 | 88.48 | 87.52 | 86.68 | 84.86 | 83.00 | 83.31 | 82.14 | 86.07
Table A17. Complete numerical results on TinyImageNet (T=10, local clients =
3) with smaller rounds in one-seeded run i.e. 2021