Case $M=VW(\lambda y.P)$.

$\begin{array}[]{rclcrcl}([\lambda x.N/k](VW(\lambda
y.P)))^{\times}&=&VW(\lambda y.[\lambda x.N/k]P)^{\times}\\\
&=&V^{\times\\!\\!\\!\times}(W^{\times\\!\\!\\!\times},y.([\lambda
x.N/k]P)^{\times})\\\
&=&V^{\times\\!\\!\\!\times}(W^{\times\\!\\!\\!\times},y.[P^{\times}\backslash
x]N^{\times})&\text{(by IH)}\\\
&=&[V^{\times\\!\\!\\!\times}(W^{\times\\!\\!\\!\times},y.P^{\times})\backslash
x]N^{\times}\\\ &=&[(VW(\lambda y.P))^{\times}\backslash
x]N^{\times}\end{array}$

∎

Theorem 5. ${\scriptstyle{CNF}}\cong{\scriptstyle{CPS}}$.

  1. 1.

For all $M,V\in{\scriptstyle{CNF}}$, ${M^{-}}^{+}=M$ and
${M^{\wr}}^{\times}=M$ and ${V^{\sim}}^{\times\\!\\!\\!\times}=V$.

  2. 2.

For all $P,M,V\in{\scriptstyle{CPS}}$, ${P^{+}}^{-}=P$ and
${M^{\times}}^{\wr}=M$ and ${V^{\times\\!\\!\\!\times}}^{\sim}=V$.

  3. 3.

If $M_{1}\to M_{2}$ in ${\scriptstyle{CNF}}$ then $M_{1}^{\wr}\to M_{2}^{\wr}$
in ${\scriptstyle{CPS}}$ (hence $M_{1}^{-}\to M_{2}^{-}$ in
${\scriptstyle{CPS}}$).

  4. 4.

If $M_{1}\to M_{2}$ in ${\scriptstyle{CPS}}$ then $M_{1}^{\times}\to
M_{2}^{\times}$ in ${\scriptstyle{CNF}}$. Hence If $P_{1}\to P_{2}$ in
${\scriptstyle{CPS}}$ then $P_{1}^{+}\to P_{2}^{+}$ in ${\scriptstyle{CNF}}$.

###### Proof.

Item 1 is straightforward by simultaneous induction on $M$ and $V$. Similarly,
item 2 is straightforward by simultaneous induction on $P,M$ and $V$. Item 3
follows by induction on the relation $M_{1}\to M_{2}$ in
${\scriptstyle{CNF}}$. We just show the base case.

Case $(\beta_{v})$.

$\begin{array}[]{rcll}((\lambda y.M)(V,x.P))^{\wr}&=&(\lambda y.\lambda
k.M^{\wr})V^{\sim}(\lambda x.P^{\wr})\\\ &\to_{\beta_{v}}&[\lambda
x.P^{\wr}/k][V^{\sim}/y]M^{\wr}\\\ &=&[\lambda
x.P^{\wr}/k]([V/y]M)^{\wr}&\text{(by Lemma \ref{lem:subs-neg})}\\\
&=&([[V/y]M\backslash x]P)^{\wr}&\text{(by Lemma \ref{lem:col-
neg})}\end{array}$

Analogously, item 4 follows by induction on the relation $M_{1}\to M_{2}$ in
${\scriptstyle{CPS}}$. We just show the base case.

Case $(\beta_{v})$.

$\begin{array}[]{rcll}((\lambda y.\lambda x.M)W(\lambda
x.N))^{\times}&=&(\lambda y.M^{\times})(W^{\times\\!\\!\\!\times},x.N^{+})\\\
&\to_{\beta_{v}}&[[W^{\times\\!\\!\\!\times}/y]M^{\times}\backslash x]N^{+}\\\
&=&[([W/y]M)^{\times}\backslash x]N^{+}&\text{(by Lemma \ref{lem:subs-
cps})}\\\ &=&([\lambda x.N/k][W/y]M)^{\times}&\text{(by Lemma \ref{lem:col-
cps})}\end{array}$

∎


11institutetext: Natural Language Processing Lab
The University of British Columbia
<EMAIL_ADDRESS><EMAIL_ADDRESS>
# BERT-Based Arabic Social Media
Author Profiling

Chiyu Zhang  Muhammad Abdul-Mageed

###### Abstract

We report our models for detecting age, language variety, and gender from
social media data in the context of the Arabic author profiling and deception
detection shared task (APDA) [32]. We build simple models based on pre-trained
bidirectional encoders from transformers (BERT). We first fine-tune the pre-
trained BERT model on each of the three datasets with shared task released
data. Then we augment shared task data with in-house data for gender and
dialect, showing the utility of augmenting training data. Our best models on
the shared task test data are acquired with a majority voting of various BERT
models trained under different data conditions. We acquire 54.72% accuracy for
age, 93.75% for dialect, 81.67% for gender, and 40.97% joint accuracy across
the three tasks.

###### Keywords:

author profiling identification, BERT, Arabic, social media

##  1 Introduction

The proliferation of social media has made it possible to collect user data in
unprecedented ways. These data can come in the form of usage and behavior
(e.g., who likes what on Facebook), network (e.g., who follows a given user on
Instagram), and content (e.g., what people post to Twitter). Availability of
such data have made it possible to make discoveries about individuals and
communities, mobilizing social and psychological research and employing
natural language processing methods. In this work, we focus on predicting
social media user age, dialect, and gender based on posted language. More
specifically, we use the total of 100 tweets from each manually-labeled user
to predict each of these attributes. Our dataset comes from the Arabic author
profiling and deception detection shared task (APDA) [32]. We focus on
building simple models using pre-trained bidirectional encoders from
transformers (BERT) [12] under various data conditions. Our results show (1)
the utility of augmenting training data, and (2) the benefit of using majority
votes from our simple classifiers.

In the rest of the paper, we introduce the dataset, followed by our
experimental conditions and results. We then provide a literature review and
conclude.

##  2 Data

For the purpose of our experiments, we use data released by the APDA shared
task organizers. The dataset is divided into train and test by organizers. The
TRAIN set is distributed with labels for the three tasks of age, dialect, and
gender. Following the standard shared tasks set up, the test set is
distributed without labels and participants were expected to submit their
predictions on test. The shared task predictions are expected by organizers at
the level of users. The distribution has 100 tweets for each user, and so each
tweet is distributed with a corresponding user id. As such, in total, the
distributed training data has 2,250 users, contributing a total of 225,000
tweets. The official task test set contains 720,00 tweets posted by 720 users.
For our experiments, we split the training data released by organizers into
90% TRAIN set (202,500 tweets from 2,025 users) and 10% DEV set (22,500 tweets
from 225 users). The age task labels come from the tagset {under-25,
between-25 and 34, above-35}. For dialects, the data are labeled with 15
classes, from the set {Algeria, Egypt, Iraq, Kuwait, Lebanon-Syria, Lybia,
Morocco, Oman, Palestine-Jordan, Qatar, Saudi Arabia, Sudan, Tunisia, UAE,
Yemen}. The gender task involves binary labels from the set {male, female}.

##  3 Experiments

As explained earlier, the shared task is set up at the user level where the
age, dialect, and gender of each user are the required predictions. In our
experiments, we first model the task at the tweet level and then port these
predictions at the user level. For our core modelling, we fine-tune BERT on
the shared task data. We also introduce an additional in-house dataset labeled
with dialect and gender tags to the task as we will explain below. As a
baseline, we use a small gated recurrent units (GRU) model. We now introduce
our tweet-level models.

###  3.1 Tweet-Level Models

####  3.1.1 Baseline GRU.

Our baseline is a GRU network for each of the three tasks. We use the same
network architecture across the 3 tasks. For each network, the network
contains a layer unidirectional GRU, with 500 units and an output linear
layer. The network is trained end-to-end. Our input embedding layer is
initialized with a standard normal distribution, with $\mu=0$, and $\sigma=1$,
i.e., $W\sim N(0,1)$. We use a maximum sequence length of 50 tokens, and
choose an arbitrary vocabulary size of 100,000 types, where we use the 100,000
most frequent words in TRAIN. To avoid over-fitting, we use dropout [43] with
a rate of 0.5 on the hidden layer. For the training, we use the Adam [22]
optimizer with a fixed learning rate of $1e-3$. We employ batch training with
a batch size of 32 for this model. We train the network for 15 epochs and save
the model at the end of each epoch, choosing the model that performs highest
accuracy on DEV as our best model. We present our best result on DEV in Table
1. We report all our results using accuracy. Our best model obtains 42.48% for
age, 37.50% for dialect, and 57.81% for gender. All models obtain best results
with 2 epochs.

####  3.1.2 BERT.

For each task, we fine-tune on the BERT-Base Muultilingual Cased model relesed
by the authors [12] 111https://github.com/google-
research/bert/blob/master/multilingual.md. The model was pre-trained on
Wikipedia of 104 languages (including Arabic) with 12 layer, 768 hidden units
each, 12 attention heads, and has 110M parameters in entire model. The
vocabulary of the model is 119,547 shared WordPices. We fine-tune the model
with maximum sequence length of 50 tokens and a batch size of 32. We set the
learning rate to $2e-5$ and train for 15 epochs. We use the same network
architecture and parameters across the 3 tasks. As Table 1 shows, comparing
with GRU, BERT is 3.16% better for age, 4.85% better for dialect, and 2.45%
higher for gender.

####  3.1.3 Data Augmentation.

To further improve the performance of our models, we introduce in-house
labeled data that we use to fine-tune BERT. For the gender classification
task, we manually label an in-house dataset of 1,100 users with gender tags,
including 550 female users, 550 male users. We obtain 162,829 tweets by
crawling the 1,100 users’ timelines. We combine this new gender dataset with
the gender TRAIN data (from shared task) to obtain an extended dataset, to
which we refer as EXTENDED_Gender. For the dialect identification task, we
randomly sample 20,000 tweets for each class from an in-house dataset gold
labeled with the same 15 classes as the shared task. In this way, we obtain
298,929 tweets (Sudan only has 18,929 tweets). We combine this new dialect
data with the shared task dialect TRAIN data to form EXTENDED_Dialect. For
both the dialect and gender tasks, we fine-tune BERT on EXTENDED_Dialect and
EXTENDED_Gender independently and report performance on DEV. We refer to this
iteration of experiments as BERT_EXT. As Table 1 shows, BERT_EXT is 2.18%
better than BERT for dialect and 0.75% better than BERT for gender. 222We note
that it was not possible for us to use external age-labeled data and hence we
do not report on the age task with this data augmentation setting.

Table 1: Tweet level results on DEV | Age | Dialect | Gender
---|---|---|---
GRU | 42.48 | 37.50 | 57.81
BERT | 45.64 | 42.35 | 60.26
BERT_EXT | - | 44.53 | 61.01

###  3.2 User-Level Models

Our afore-mentioned models identify user’s profiling on the tweet-level,
rather than directly detecting the labels of a user. Hence, we follow the work
of Zhang & Abdul-Mageed [47] to identify user-level labels. For each of the
three tasks, we use tweet-level predicted labels (and associated softmax
values) as a proxy for user-level labels. For each predicted label, we use the
softmax value as a threshold for including only highest confidently predicted
tweets. Since in some cases softmax values can be low, we try all values
between 0.00 and 0.99 to take a softmax-based majority class as the user-level
predicted label, fine-tuning on our DEV set. Using this method, we acquire the
following results at the user level: BERT models obtain an accuracy of 55.56%
for age, 96.00% for dialect, and 80.00% for gender. BERT_EXT models achieve
95.56% accuracy for dialect and 84.00% accuracy for gender.

###  3.3 APDA@FIRE2019 submission

First submission. For the shared task submission, we use the predictions of
BERT_EXT as out first submission for gender and dialect, but only BERT for age
(since we have no BERT_EXT models for age, as explained earlier). In each
case, we acquire results at tweet-level first, then port the labels at the
user-level as explained in the previous section. For our second and third
submitted models, we also follow this method of going from tweet to user
level. Second submission. We combine our DEV data with our EXTENDED_Dialect
and EXTENDED_Gender data, for dialect and gender respectively, and train our
second submssions for the two tasks. For age second submsision, we concatenate
DEV data to TRAIN and fine-tune the BERT model. We refer to the settings for
our second submission models collectively as BERT_EXT+DEV.

Third submission. Finally, for our third submission, we use a majority vote of
(1) first submission, (2) second submission, and (3) predictions from our
user-level BERT model. These majority class models (i.e., our third
submission) achieve best results on the official test data. We acquire 54.72%
accuracy for age, 81.67% accuracy for gender, 93.75% accuracy for dialect, and
40.97% joint accuracy.

Table 2: Results of our submissions on official test data (user level) | Exp. Condition | Age | Dialect | Gender | Joint
---|---|---|---|---|---
Submission 1 | BERT_EXT | 54.72 | 93.33 | 77.08 | 38.75
Submission 2 | BERT_EXT+DEV | 54.72 | 92.64 | 81.67 | 40.97
Submission 3 | MAJ_CLASS | 54.72 | 93.75 | 81.67 | 40.97

##  4 Related Works

Arabic. Arabic is a term that refers to a collection of languages, varieties,
and dialects. The standard variety, Modern Standard Arabic (MSA), is the one
usually used in formal communication and educational settings. Arabic also has
a wide range of under-studied varieties and dialects that classically used to
be categorized in a coarse-grained fashion (e.g., Levantine, North African)
[17, 46, 1, 13, 2]. More recent treatments focus on fine-grained
categorizations such as country and city levels [25, 39, 3, 40, 31, 47].
Differences between varieties of Arabic happen at various linguistic levels,
including including phonological, morphological, lexical, and syntactic [19,
6, 27, 1].

Social Media Author Profiling. Author profiling is the term usually used to
refer to detecting a host of attributes of (often social media) users. This
include identifying attributes such as age, gender, educational level,
economic class, stance or ideology [10, 30], personality [42, 23, 7, 18],
moral traitstraits [21, 28], and other socilogical and psychological
constructs. Author profiling based on text [15, 5] is rooted in computational
stylometry [16, 44] and has traces in the early work of Holmes [20]. The task
of author profiling has also been approached from network perspective where
cues based on friending, following, mentioning, and commenting have been
leveraged for identifying author attributes [24] for author profiling. In
addition, the PAN author profiling shared task [34, 33, 37, 36, 35] was
established to advance related work. More information about PAN can be found
in [29].

Age and Gender. A number of studies have been conducted on English-based age
and gender detection, including [38, 14, 8, 45, 11]. Many of these works use
feature engineering such as text n-gram and topic models [42, 41]. In these
works, age is either cast as a multi-class classification task with, e.g.,
labels from the set {10-19, 20-29, 30-39} or as a regression task [26]. Other
works model age with both classification and regression combined [9]. With
rare exceptions [35, 4], we do not know of work on Arabic targeting age and
gender.

##  5 Conclusion

In this work, we described our submitted models to the Arabic author profiling
and deception detection shared task (APDA) [32]. We focused on detecting age,
dialect, and gender using BERT models under various data conditions, showing
the utility of additional, in-house data on the task. We also showed that a
majority vote of our models trained under different conditions outperforms
single models on the official evaluation. In the future, we will investigate
automatically extending training data for these tasks as well as better
representation learning methods.

##  6 Acknowledgement

We acknowledge the support of the Natural Sciences and Engineering Research
Council of Canada (NSERC), the Social Sciences Research Council of Canada
(SSHRC), and Compute Canada (www.computecanada.ca).

## References

  * [1] Abdul-Mageed, M.: Subjectivity and sentiment analysis of Arabic as a morophologically-rich language. Ph.D. thesis, Indiana University (2015)
  * [2] Abdul-Mageed, M.: Modeling arabic subjectivity and sentiment in lexical space. Information Processing & Management (2017)
  * [3] Abdul-Mageed, M., Alhuzali, H., Elaraby, M.: You tweet what you speak: A city-level dataset of arabic dialects. In: Proceedings of the Eleventh International Conference on Language Resources and Evaluation (LREC-2018) (2018)
  * [4] Alrifai, K., Rebdawi, G., Ghneim, N.: Arabic tweeps gender and dialect prediction. In: CLEF (Working Notes) (2017)
  * [5] Argamon, S.E.: Register in computational language research. Register Studies 1(1), 100–135 (2019)
  * [6] Bassiouney, R.: Arabic sociolinguistics. Edinburgh University Press (2009)
  * [7] Bleidorn, W., Hopwood, C.J.: Using machine learning to advance personality assessment and theory. Personality and Social Psychology Review p. 1088868318772990 (2018)
  * [8] Burger, J.D., Henderson, J., Kim, G., Zarrella, G.: Discriminating gender on twitter. In: Proceedings of the conference on empirical methods in natural language processing. pp. 1301–1309. Association for Computational Linguistics (2011)
  * [9] Chen, J., Cheng, L., Yang, X., Liang, J., Quan, B., Li, S.: Joint learning with both classification and regression models for age prediction. In: Journal of Physics: Conference Series. vol. 1168, p. 032016. IOP Publishing (2019)
  * [10] Colleoni, E., Rozza, A., Arvidsson, A.: Echo chamber or public sphere? predicting political orientation and measuring political homophily in twitter using big data. Journal of communication 64(2), 317–332 (2014)
  * [11] Daneshvar, S., Inkpen, D.: Gender identification in twitter using n-grams and lsa. In: Proceedings of the Ninth International Conference of the CLEF Association (CLEF 2018) (2018)
  * [12] Devlin, J., Chang, M.W., Lee, K., Toutanova, K.: Bert: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint arXiv:1810.04805 (2018)
  * [13] Elaraby, M., Abdul-Mageed, M.: Deep models for arabic dialect identification on benchmarked data. In: Proceedings of the Fifth Workshop on NLP for Similar Languages, Varieties and Dialects (VarDial 2018). pp. 263–274 (2018)
  * [14] Flekova, L., Preoţiuc-Pietro, D., Ungar, L.: Exploring stylistic variation with age and income on twitter. In: Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 2: Short Papers). vol. 2, pp. 313–319 (2016)
  * [15] Gamon, M.: Linguistic correlates of style: authorship classification with deep linguistic analysis features. In: Proceedings of the 20th international conference on Computational Linguistics. p. 611. Association for Computational Linguistics (2004)
  * [16] Goswami, S., Sarkar, S., Rustagi, M.: Stylometric analysis of bloggers’ age and gender. In: Third international AAAI conference on weblogs and social media (2009)
  * [17] Habash, N.Y.: Introduction to arabic natural language processing. Synthesis Lectures on Human Language Technologies 3(1), 1–187 (2010)
  * [18] Hinds, J., Joinson, A.: Human and computer personality prediction from digital footprints. Current Directions in Psychological Science p. 0963721419827849 (2019)
  * [19] Holes, C.: Modern Arabic: Structures, functions, and varieties. Georgetown University Press (2004)
  * [20] Holmes, D.I.: The evolution of stylometry in humanities scholarship. Literary and linguistic computing 13(3), 111–117 (1998)
  * [21] Johnson, K., Goldwasser, D.: Classification of moral foundations in microblog political discourse. In: Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers). pp. 720–730 (2018)
  * [22] Kingma, D., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * [23] Matz, S.C., Kosinski, M., Nave, G., Stillwell, D.J.: Psychological targeting as an effective approach to digital mass persuasion. Proceedings of the national academy of sciences 114(48), 12714–12719 (2017)
  * [24] Mitrou, L., Kandias, M., Stavrou, V., Gritzalis, D.: Social media profiling: A panopticon or omniopticon tool? In: Proc. of the 6th Conference of the Surveillance Studies Network. Barcelona, Spain (2014)
  * [25] Mubarak, H., Darwish, K.: Using twitter to collect a multi-dialectal corpus of arabic. In: Proceedings of the EMNLP 2014 Workshop on Arabic Natural Language Processing (ANLP). pp. 1–7 (2014)
  * [26] Nguyen, D., Smith, N.A., Rosé, C.P.: Author age prediction from text using linear regression. In: Proceedings of the 5th ACL-HLT Workshop on Language Technology for Cultural Heritage, Social Sciences, and Humanities. pp. 115–123. Association for Computational Linguistics (2011)
  * [27] Palva, H.: Dialects: classification. Encyclopedia of Arabic Language and Linguistics 1, 604–613 (2006)
  * [28] Pang, D., Eichstaedt, J.C., Buffone, A., Slaff, B., Ruch, W., Ungar, L.H.: The language of character strengths: Predicting morally valued traits on social media. Journal of personality (2019)
  * [29] Potthast, M., Rosso, P., Stamatatos, E., Stein, B.: A decade of shared tasks in digital text forensics at pan. In: European Conference on Information Retrieval. pp. 291–300. Springer (2019)
  * [30] Preoţiuc-Pietro, D., Liu, Y., Hopkins, D., Ungar, L.: Beyond binary labels: political ideology prediction of twitter users. In: Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers). pp. 729–740 (2017)
  * [31] Qwaider, C., Saad, M., Chatzikyriakidis, S., Dobnik, S.: Shami: A corpus of levantine arabic dialects. In: Proceedings of the Eleventh International Conference on Language Resources and Evaluation (LREC-2018) (2018)
  * [32] Rangel, F., Rosso, P., Charfi, A., Zaghouani, W., Ghanem, B., Sánchez-Junquera, J.: Overview of the track on author profiling and deception detection in arabic. In: Mehta P., Rosso P., Majumder P., Mitra M. (Eds.) Working Notes of the Forum for Information Retrieval Evaluation (FIRE 2019). CEUR Workshop Proceedings. In: CEUR-WS.org, Kolkata, India, December 12-15 (2019)
  * [33] Rangel, F., Rosso, P., Chugur, I., Potthast, M., Trenkmann, M., Stein, B., Verhoeven, B., Daelemans, W.: Overview of the 2nd author profiling task at pan 2014. In: CLEF 2014 Evaluation Labs and Workshop Working Notes Papers, Sheffield, UK, 2014. pp. 1–30 (2014)
  * [34] Rangel, F., Rosso, P., Koppel, M., Stamatatos, E., Inches, G.: Overview of the author profiling task at pan 2013. In: CLEF Conference on Multilingual and Multimodal Information Access Evaluation. pp. 352–365. CELCT (2013)
  * [35] Rangel, F., Rosso, P., Potthast, M., Stein, B.: Overview of the 5th author profiling task at pan 2017: Gender and language variety identification in twitter. Working Notes Papers of the CLEF (2017)
  * [36] Rangel, F., Rosso, P., Verhoeven, B., Daelemans, W., Potthast, M., Stein, B.: Overview of the 4th author profiling task at pan 2016: cross-genre evaluations. In: Working Notes Papers of the CLEF 2016 Evaluation Labs. CEUR Workshop Proceedings/Balog, Krisztian [edit.]; et al. pp. 750–784 (2016)
  * [37] Rangel Pardo, F.M., Celli, F., Rosso, P., Potthast, M., Stein, B., Daelemans, W.: Overview of the 3rd author profiling task at pan 2015. In: CLEF 2015 Evaluation Labs and Workshop Working Notes Papers. pp. 1–8 (2015)
  * [38] Rao, D., Yarowsky, D., Shreevats, A., Gupta, M.: Classifying latent user attributes in twitter. In: Proceedings of the 2nd international workshop on Search and mining user-generated contents. pp. 37–44. ACM (2010)
  * [39] Sadat, F., Kazemi, F., Farzindar, A.: Automatic identification of arabic language varieties and dialects in social media. Proceedings of SocialNLP p. 22 (2014)
  * [40] Salameh, M., Bouamor, H.: Fine-grained arabic dialect identification. In: Proceedings of the 27th International Conference on Computational Linguistics. pp. 1332–1344 (2018)
  * [41] Sap, M., Park, G., Eichstaedt, J., Kern, M., Stillwell, D., Kosinski, M., Ungar, L., Schwartz, H.A.: Developing age and gender predictive lexica over social media. In: Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP). pp. 1146–1151 (2014)
  * [42] Schwartz, H.A., Eichstaedt, J.C., Kern, M.L., Dziurzynski, L., Ramones, S.M., Agrawal, M., Shah, A., Kosinski, M., Stillwell, D., Seligman, M.E., et al.: Personality, gender, and age in the language of social media: The open-vocabulary approach. PloS one 8(9), e73791 (2013)
  * [43] Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., Salakhutdinov, R.: Dropout: a simple way to prevent neural networks from overfitting. The Journal of Machine Learning Research 15(1), 1929–1958 (2014)
  * [44] Verhoeven, B.: Two authors walk into a bar: studies in author profiling. Ph.D. thesis, University of Antwerp (2018)
  * [45] Verhoeven, B., Daelemans, W., Plank, B.: Twisty: a multilingual twitter stylometry corpus for gender and personality profiling. In: Proceedings of the 10th Annual Conference on Language Resources and Evaluation (LREC 2016)/Calzolari, Nicoletta [edit.]; et al. pp. 1–6 (2016)
  * [46] Versteegh, K.: The arabic language. Edinburgh University Press (2014)
  * [47] Zhang, C., Abdul-Mageed, M.: No army, no navy: Bert semi-supervised learning of arabic dialects. In: Proceedings of the Fourth Arabic Natural Language Processing Workshop. pp. 279–284 (2019)


# Universal Upper End of the Stellar Initial Mass Function in the Young and
Compact LEGUS clusters

Dooseok Escher Jung Department of Astronomy, University of Massachusetts
Amherst, 710 North Pleasant Street, Amherst, MA 01003, USA  Daniela Calzetti
Department of Astronomy, University of Massachusetts Amherst, 710 North
Pleasant Street, Amherst, MA 01003, USA  Matteo Messa Observatoire de Genève,
Université de Genève, 51 Ch. des Maillettes, CH-1290 Versoix, Switzerland
Mark Heyer Department of Astronomy, University of Massachusetts Amherst, 710
North Pleasant Street, Amherst, MA 01003, USA  Mattia Sirressi Department of
Astronomy, Oscar Klein Centre, Stockholm University, AlbaNova, SE-106 91
Stockholm, Sweden  Sean T. Linden Department of Astronomy, University of
Massachusetts Amherst, 710 North Pleasant Street, Amherst, MA 01003, USA
Angela Adamo Department of Astronomy, Oscar Klein Centre, Stockholm
University, AlbaNova, SE-106 91 Stockholm, Sweden  Rupali Chandar Ritter
Astrophysical Research Center, University of Toledo, Toledo, OH 43606, USA
Michele Cignoni Department of Physics, University of Pisa, Largo Pontecorvo,
3, I-56127 Pisa, Italy  David O. Cook Caltech/IPAC, 1200 E. California
Boulevard, Pasadena, CA 91125, USA  Clare L. Dobbs School of Physics and
Astronomy, University of Exeter, Stocker Road, Exeter EX4 4QL, UK  Bruce G.
Elmegreen IBM T. J. Watson Research Center, 1101 Kitchawan Road, Yorktown
Heights, NY 10598, USA  Aaron S. Evans Department of Astronomy, University of
Virginia, 530 McCormick Road, Charlottesville, VA 22904, USA  Michele
Fumagalli Dipartimento di Fisica G. Occhialini, Università degli Studi di
Milano-Bicocca, Piazza della Scienza 3, 20126 Milano, Italy  INAF -
Osservatorio Astronomico di Trieste, via G. B. Tiepolo 11, 34143 Trieste,
Italy  John S. Gallagher III Department of Astronomy, University of Wisconsin,
475 N. Charter Street, Madison, WI, 53706, USA  Deidre A. Hunter Lowell
Observatory, 1400 West Mars Hill Road, Flagstaff, AZ 86001, USA  Kelsey E.
Johnson Department of Astronomy, University of Virginia, 530 McCormick Road,
Charlottesville, VA 22904, USA  Robert C. Kennicutt, Jr. Department of
Astronomy and Steward Observatory, University of Arizona, Tucson, AZ 85721,
USA  Department of Physics and Astronomy, Texas A&M University, College
Station, TX 77843, USA  Institute of Astronomy, University of Cambridge,
Madingley Road, Cambridge CB3 0HA, UK  Mark R. Krumholz Research School of
Astronomy and Astrophysics, Australian National University, Canberra, ACT
2601, Australia  Daniel Schaerer Observatoire de Genève, Université de Genève,
Chemin Pegasi 51, 1290 Versoix, Switzerland  CNRS, IRAP, 14 Avenue E. Belin,
31400 Toulouse, France  Elena Sabbi Space Telescope Science Institute, 3700
San Martin Drive, Baltimore, MD, 21218, USA  Linda J. Smith Space Telescope
Science Institute, 3700 San Martin Drive, Baltimore, MD, 21218, USA  Monica
Tosi INAF – OAS Osservatorio di Astrofisica e Scienza dello Spazio, I-40129
Bologna, Italy  Aida Wofford Instituto de Astronomia, Universidad Nacional
Autonoma de Mexico, Unidad Académica en Ensenada, 22860 Ensenada, Mexico

(Received April 10, 2023; Revised July 25, 2023; Accepted July 26, 2023)

###### Abstract

We investigate the variation in the upper end of stellar initial mass function
(uIMF) in 375 young and compact star clusters in five nearby galaxies within
$\sim 5$ Mpc. All the young stellar clusters (YSCs) in the sample have ages
$\lesssim 4$ Myr and masses above 500 $M_{\odot}$, according to standard
stellar models. The YSC catalogs were produced from Hubble Space Telescope
images obtained as part of the Legacy ExtraGalactic UV Survey (LEGUS) Hubble
treasury program. They are used here to test whether the uIMF is universal or
changes as a function of the cluster’s stellar mass. We perform this test by
measuring the H$\alpha$ luminosity of the star clusters as a proxy for their
ionizing photon rate, and charting its trend as a function of cluster mass.
Large cluster numbers allow us to mitigate the stochastic sampling of the
uIMF. The advantage of our approach relative to previous similar attempts is
the use of cluster catalogs that have been selected independently of the
presence of H$\alpha$ emission, thus removing a potential sample bias. We find
that the uIMF, as traced by the H$\alpha$ emission, shows no dependence on
cluster mass, suggesting that the maximum stellar mass that can be produced in
star clusters is universal, in agreement with previous findings.

galaxies: nearby galaxies – galaxies: star clusters: general – galaxies: star
formation – stars: initial mass function – stars: massive

††journal: ApJ††facilities: HST(ACS and WFC3)††software: SExtractor (Bertin &
Arnouts, 1996), Astropy (The Astropy Collaboration et al., 2022)

##  1 Introduction

The stellar Initial Mass Function (IMF) is a valuable tool that enables the
derivation of physical quantities such as star formation rates (SFRs)
(Kennicutt, Jr., 1998) and star formation histories (Tolstoy et al., 2009)
from calibrations of light emission. Similarly, the assumed IMF affects
theoretical predictions of chemical evolution models on galactic scales
(Romano et al., 2005) and the stellar feedback into the
interstellar/circumgalactic media (Krumholz et al., 2019). To effectively use
this tool, it is imperative to establish whether the functional form of the
stellar IMF is universal or varies with local or global environments.

Salpeter (1955) was the first to introduce the concept of the stellar IMF, the
mass distribution of stars at birth, as a single power-law with a slope
$\Gamma\sim 1.35$ in the mass range $\sim 3\mbox{--}15\>M_{\odot}$, in the
form:

$\Phi(\mathrm{log}\,m)=\frac{dN}{d\,\mathrm{log}\,m}\propto m^{-\Gamma}$ (1)

where $\mathrm{log}\,m$ is the logarithmic mass of a star and $N$ is the
number of stars between $\mathrm{log}\,m$ and
$\mathrm{log}\,m+d\>\mathrm{log}\,m$.

Miller & Scalo (1979) found that the slope is shallower at masses lower than
those examined by Salpeter (1955), and fitted the IMF to both a segmented
power law and a log-normal function. Kroupa (2001) revised the multi-segmented
power law to have $\Gamma=-0.7$ for substellar objects ($M<0.08\>M_{\odot}$),
$\Gamma=0.3$ in the lower-stellar mass regime
($0.08\>\mbox{--}\>0.5\>M_{\odot}$) and $\Gamma=1.3$ for stars with
$M>0.5\>M_{\odot}$. Chabrier (2003) also fitted the IMF to a log-normal
function:

$\phi(m)\sim\mathrm{exp}(-\,\frac{\mathrm{log}\,m-\mathrm{log}\,m_{c}}{2\sigma^{2}})$
(2)

where $m_{c}$ is a characteristic mass and $\sigma$ is a dispersion of the
log-normal distribution. The Chabrier-type IMF is almost indistinguishable
from a Kroupa-type IMF (Dabringhausen et al., 2008). In the context of this
paper, therefore, the IMF is termed ‘universal’ if its shape, slopes in the
case of the functional form in Kroupa (2001), and upper/lower stellar limits
are independent of the local or global environment.

The debate on whether or not the IMF depends on other parameters, such as
metallicity, environmental properties, etc., has been ongoing for several
decades. Some authors argue that the stellar IMF is universal (Massey, 2003;
Bastian et al., 2010; Oey, 2011; Weisz et al., 2012), while other studies find
non-universal IMFs, with variations that are dependent on the galactic
environment (Fardal et al., 2007; Wilkins et al., 2008; Van Dokkum & Conroy,
2011; Kroupa et al., 2012; Cappellari et al., 2012; Geha et al., 2013). For
example, using star counts, Schneider et al. (2018) found that the slope of
the high-mass IMF in 30 Doradus, a giant H II region in the Large Magellanic
Cloud (LMC), is shallower than the Salpeter/Kroupa value. Counting stars is
generally considered the most accurate method to derive the IMF.

However, counting stars individually is challenging, especially when studying
potential IMF variations in the massive star regime due to the evolution of
stars and star clusters. The most massive stars are short lived ($\lesssim 3$
Myr) and their numbers are significantly lower than those of low-mass stars.
Moreover, they are mainly found in the dense environments of star clusters,
where crowding becomes an issue (Kalari et al., 2022). Mass segregation that
occurs in gravitationally bound systems makes more massive stars sink toward
the center of the stellar clusters over time, and causes low-mass stars to be
dynamically ejected. As a result, in the crowded areas of the cluster centers,
massive stars are difficult to distinguish individually, and low-mass stars
can remain undetected due to their faintness or the ejection from the cluster
itself (Maíz-Apellniz, 2008; Ascenso et al., 2009). Thus, the problem of
crowding complicates attempts at resolving high-mass stars even with the
angular resolution of the Hubble Space Telescope (HST). The proximity of the
Magellanic Clouds (distances $\sim 50\mbox{--}60$ kpc) has helped overcome
some of these complications and stars heavier than $\sim 150\>M_{\odot}$ have
been found in the core of R136, the central cluster in 30 Doradus (Crowther et
al., 2010; Brands et al., 2022; Kalari et al., 2022). By using semi-resolved
star clusters of intermediate masses
($10^{3}\>M_{\odot}\mbox{--}\>10^{4}M_{\odot}$) in M31 ($\sim 1$ Mpc), Weisz
et al. (2015) found that the slope of the high-mass IMF ($\gtrsim
2\mbox{--}3\>M_{\odot}$) in the clusters is close to the Kroupa value
($\Gamma=1.45^{+0.03}_{-0.06}$) for the same stellar mass range.

An alternative method to individual star counting consists of measuring the
integrated or summed IMF of a young stellar population. This method avoids the
issues of the resolved-star approach; however, it comes with challenges of its
own. Ideally, the population needs to be very young ($\lesssim 3$ Myr, to
contain massive stars), coeval (to provide one representation of the IMF), and
homogeneous (to sample the same IMF across the population). Young star
clusters (YSCs) provide, at least in principle, such a population. The
challenges consist in (1) deriving accurate masses and ages for the star
clusters, and (2) deriving the IMF in an unambiguous way from the integrated
light of the young cluster population. The latter is challenging even for
massive stars, which tend to dominate the light output of YSCs. Furthermore,
an accurate census of massive stars can be hindered by feedback, which causes
high infant mortality rates in stellar clusters; less than 20% of cluster
systems survive longer than 10 Myr (Lada & Lada, 2003). A dispersed cluster
becomes difficult to recognize as a single-age entity for the purpose of
measuring the IMF.

Dispersed clusters contribute their stars to the field and produce an
integrated mass function, sometimes called the present-day mass function
(PDMF) (Miller & Scalo, 1979). Because high-mass stars disappear before the
low-mass stars have scattered away from their birthplace in the galaxy, the
PDMF in any limited region has a complicated relationship to the IMF, and they
are not equivalent. Also, if the SFR increases with time, then the summed mass
function for all the remaining stars will be shallower than the original
cluster IMFs (Elmegreen & Scalo, 2006). Even for a steady SFR, and accounting
all the stars that ever formed in, e.g., a young, pre-supernova region, the
summed IMF from many clusters is not necessarily the same as the IMFs of the
individual clusters. They are the same only if stellar mass is randomly
sampled from a universal IMF with a power law shape to infinitely high mass,
or if the cluster mass function (CMF) has a slope equal to $-2$ or shallower
(Elmegreen, 2006).

If the maximum mass of the star(s) depends physically on the cluster mass,
increasing with cluster or molecular cloud mass for example (Larson, 1982;
Kroupa & Weidner, 2003) and not just randomly by the size of sample effect
(Elmegreen, 1983), then the summed IMF can be steeper than the individual IMF.
This is because a lot of clusters sum together to yield low-mass stars but
only the highest-mass clusters sum together to yield the high-mass stars. To
get purely stochastic sampling (the randomly sampled population of massive
stars that can occur in clusters with mass $<10^{5}\>M_{\odot}$, i.e. with too
low a mass to fully sample the IMF) with a summed IMF slope equal to the
individual cluster’s slope, some low-mass clusters have to produce high-mass
stars. Because cluster-forming cloud complexes are usually much more massive
than individual stars, there is typically enough gas to do this. Available
observations of the Milky Way bulge, large star fields and whole galaxies
confirm that the summed IMF from composite star forming regions is about equal
to the IMF from individual star forming regions (Elmegreen, 2006).

Pflamm-Altenburg et al. (2009) considered the implications of a physical,
rather than stochastic, relation between the mass of the the most massive star
($m_{\mathrm{max}}$) and the mass of the embedded cluster ($M_{\mathrm{ecl}}$)
that hosts the star(s) with $m_{\mathrm{max}}$, generally expressed as a
$m_{\mathrm{max}}\mbox{--}M_{\mathrm{ecl}}$ relation. Embedded clusters are
those clusters that are still retaining their original gas, and are therefore
about twice as massive as their exposed counterparts, which we term YSCs
(Kroupa & Boily, 2002). Pflamm-Altenburg et al. (2009) developed models for an
integrated galactic IMF (IGIMF) to test the
$m_{\mathrm{max}}\mbox{--}M_{\mathrm{ecl}}$ relation on galactic scales. They
found that this relation can explain observations such as the systematic
dearth of H$\alpha$ emission in low-mass galaxies, if low-mass, low-density
galaxies mostly form low-mass clusters which, in turn, form a smaller numbers
of massive stars than expected from a fully sampled IMF. Weidner et al. (2010)
summarized several theoretical studies that predict the
$m_{\mathrm{max}}\mbox{--}M_{\mathrm{ecl}}$ relation. In addition, they
reported several observational results that show that $m_{\mathrm{max}}$ and
$M_{\mathrm{ecl}}$ are correlated in resolved stars and clusters, suggesting
that the $m_{\mathrm{max}}\mbox{--}M_{\mathrm{ecl}}$ relation has general
validity. Oey (2011) further extended this concept to explain the steepness of
the Salpeter IMF in star clusters as the result of an initial clump mass
function with a slope of $-2$ combined with an inability of low-mass clumps to
form massive stars. Similar considerations can be used to explain the steeper
H II region luminosity function in the interarm regions relative to the spiral
arms of galaxies (Oey & Clarke, 1998). Yan et al. (2017) found that the IGIMF
model explains well the correlation between SFR and the mass of the most
massive star cluster ($M_{\mathrm{ecl,max}}$) in the host galaxy, and produce
a prediction for the expected supernova rate in very-low SFR galaxies.
Jeřábková et al. (2018) applied the IGIMF framework to calculate grids of
models that depend on both metallicity and SFR. Their IGIMF model can
potentially explain several observational results, including the evolution of
the time-varying IMF in massive elliptical galaxies.

Conversely, Fumagalli et al. (2011), suggested that the IGIMF cannot fully
describe the observed large scatter in the H$\alpha$/FUV ratio for low-
luminosity galaxies regardless of their SFR. After adding stochasticity and
time evolution to a Kroupa (2001) IMF and convolving the IMF with the CMF,
Fumagalli et al. (2011) reproduced the observed
$L_{\mathrm{H}\alpha}/L_{\mathrm{FUV}}$ ratios in dwarf galaxies without the
need to implement the IGIMF framework. Similar results were obtained by
Hermanowicz et al. (2013) for the $L_{\mathrm{H}\alpha}/L_{\mathrm{FUV}}$
ratios of a large sample of H II regions extracted from eight nearby star-
forming galaxies.

In an attempt to discriminate among these different results, Calzetti et al.
(2010) used the integrated light from the YSC population to test the
universality of the upper end of the stellar IMF (uIMF), bypassing the
traditional method of individual star counts. Based on Corbelli et al. (2009),
individual YSCs were divided in bins of stellar mass and the properties of the
clusters were combined together within each mass bin to emulate a high-mass
cluster. One of the combined properties is the luminosity of the hydrogen
recombination-line intensity, H$\alpha$, as a tracer of the ionizing photon
rate, which is a proxy for the number of massive and ionizing stars present in
the star clusters located within H II regions (Oey & Clarke, 1998). This
method allows one to investigate the universality of the uIMF in star clusters
across a wide range of masses, including those below $\sim 3000\>M_{\odot}$,
which are subject to significant stochastic sampling (Cervino et al., 2002).
The application of this method to the cluster populations in the nearby
galaxies NGC 5194, NGC 4214 and M 83 (Calzetti et al., 2010; Andrews et al.,
2013, 2014) shows that the clusters with even the lowest mass down to
$500\>M_{\odot}$ produce high-mass stars up to $120\>M_{\odot}$. Thus,
according to these authors, the maximal star mass is not dependent on the mass
of its parent cluster, implying a disagreement with a physically-based
$m_{\mathrm{max}}\mbox{--}M_{\mathrm{ecl}}$ relation (Weidner et al., 2010).

Intriguingly, Weidner et al. (2014) argue that, in the case of NGC 4214, the
addition of stochastic sampling helps reconcile the trends observed in the
cluster population of this galaxy with the
$m_{\mathrm{max}}\mbox{--}M_{\mathrm{ecl}}$ relation, opposite to the
conclusions of Andrews et al. (2013). However, NGC 4214 alone is not a good
test sample for the $m_{\mathrm{max}}\mbox{--}M_{\mathrm{ecl}}$ relation due
to its relatively low SFR ($0.2\>M_{\odot}\,\mathrm{yr}^{-1}$), which implies
low numbers of observed clusters, and thus, low-number statistics.

Progress on this controversy can be aided by the data products from the Legacy
ExtraGalactic UV Survey (LEGUS), a Hubble treasury program that has obtained
imaging data for $\sim 50$ nearby star-forming galaxies (distances $\sim
3.5\mbox{--}18$ Mpc) in five bands from NUV to I (Calzetti et al., 2015a). As
part of the effort, star cluster catalogs for 31 LEGUS galaxies have been
produced, which include visual identification, location and photometry for
each star cluster; for clusters detected in more than four bands, physical
parameters: age, mass, and color excess, $E(B-V)$, derived from spectral
energy distribution fitting (Adamo et al., 2017), are included as well.

The LEGUS collaboration has already investigated the general properties of the
YSC populations in the observed galaxies, including the CMF and Cluster
Luminosity Functions (CLF). Adamo et al. (2017), Messa et al. (2018a) and
(Cook et al., 2019) generalized the $-2$ slope power law behavior of the CMFs
and CLFs to a much larger sample of galaxies than previously analyzed
including dwarf galaxies. Adamo et al. (2017) and Messa et al. (2018a, b)
determined the value of the upper mass truncation to the CMF, and investigated
its dependency on galactic environment. Adamo et al. (2017) also showed that
the different evolutionary timescales between YSCs and compact associations
track the hierarchical process of star formation (Grasha et al., 2015, 2017a,
2017b). The investigation of the spatial correlation between the LEGUS YSCs
and molecular clouds in M51 and NGC 7793 showed that the younger the clusters,
the closer they tend to be located to the molecular clouds, spiral arms and
the galactic center (Grasha et al., 2018, 2019). Orozco-Duarte et al. (2022)
developed a photometric library of synthetic star clusters with stochastically
sampled IMFs; by comparing their models with the young star clusters in NGC
7793 from the LEGUS catalogs, those authors determine that physical quantities
(ages, masses, extinctions) are highly model-dependent and most clusters have
multi-peaked age probability distributions.

The present study builds on the results of earlier studies, with the goal of
testing whether or not the uIMF varies as a function of the stellar mass of
YSCs, using the large cluster catalogs of the LEGUS collaboration (Calzetti et
al., 2015a; Adamo et al., 2017). We combine the star clusters from the
catalogs of selected LEGUS galaxies with archival H$\alpha$ imaging data to
repeat the analysis of Andrews et al. (2013, 2014). Our approach differs from
those papers in that the clusters are selected independently of whether they
have detected line emission or not. In this way, we do not bias our analysis,
since we do not exclude low-mass clusters that may be missing massive,
ionizing stars due to stochastic sampling. Section 2 explains the data
collection and the image processing. Section 3 describes how we use our data
to get to a description of the uIMF. In Section 4, we discuss our results and
their implications.

##  2 Data

Figure 1: Left: the locations of the selected LEGUS clusters by class (class 1
as red, class 2 as green and class 3 as blue circles) with ages within 1–4 Myr
in NGC 4449 are shown on the F555W image. Right: a zoomed-in image of the
selected area (the black solid square) in the left panel. Candidates
identified as class 4 are excluded here and the morphology-based cluster
classification is explained in Section 2.1.

The Legacy ExtraGalatic UV Survey (LEGUS) is a Hubble Space Telescope (HST)
Treasury Program (Cycle 21) which collected five-band images (UV, U, B, V and
I) of 50 star-forming galaxies within 12 Mpc with the Wide-Field Camera 3
(WFC3) using the filters F275W, F336W, F438W, F555W and F814W (Calzetti et
al., 2015a). Archival images obtained with the Advanced Camera for Surveys
(ACS) wide field in the filters F435W, F606W and F814W are also used when they
are available to complement the LEGUS observations. The pixel scale of the
ACS/WFC is slightly larger, 0′′.05 pixel-1, than that of WFC3 (0′′.0396
pixel-1), but the LEGUS collaboration aligned and resampled the images to the
WFC3 ones.

YSC catalogs were produced using a combination of SExtractor (Bertin &
Arnouts, 1996) for the automatic search of extended sources and human
identification. Spectral energy distribution (SED) fits for all sources
detected with signal-to-noise (S/N) of 3 or better in at least four bands were
performed, and the sources assigned ages, masses, and extinctions, with
uncertainties. The YSC catalogs are complete out to $\sim 200$ Myr in age and
down to masses between 500 and $5000\>M_{\odot}$, depending on the age of the
source. Details on the catalog production and characteristics can be found in
Adamo et al. (2017). We supplement the LEGUS observations with archival images
in the WFC3/F657N or ACS/F658N and the adjacent WFC3/F547M or ACS/F550M (as
available) filters to measure the hydrogen recombination line (H$\alpha$).

###  2.1 Galaxy and Cluster Selection

As our goal is to measure the ionizing photon rate of each YSC, we need to
separate the ionization regions of neighboring YSCs from each other. This sets
an upper limit of about 5 Mpc to the maximum distance of the galaxies whose
YSC-produced ionization regions we can resolve with HST. Within this distance
limit, individual clusters can be detected down to $\sim 500\>M_{\odot}$
(Andrews et al., 2014) which allows our study to investigate whether massive
stars are able to be formed in low-mass clusters. We also require that the
YSCs exhibit low dust color excess, $E(B-V)<0.1$ mag to avoid heavy impact
from extinction corrections. This choice is likely to select against the
youngest (most embedded) YSCs, but it also enables accurate accounting of the
ionizing photon rate. The clusters in the LEGUS catalogs show on average low
extinction values, typically $E(B-V)\lesssim 0.2$ mag (Adamo et al., 2017;
Messa et al., 2018a). After applying these selection criteria, our sample
includes five nearby star-forming galaxies with distance $\lesssim$ 5 Mpc. The
galaxies span an order of magnitude in SFRs
($0.1\mbox{--}1.2\>M_{\odot}\>\mathrm{yr}^{-1}$) and stellar masses
($10^{8}\mbox{--}10^{9}\>M_{\odot}$); these and other global properties are
listed in Table LABEL:tab:object.

Table 1: Properties of the LEGUS galaxies in this study.

Namea | RAa | Deca | Morph.a | Dist.b (Mpc) |  $M_{\bigstar}^{b}$ [$M_{\odot}$] | $E(B-V)^{a}$ | $Z^{c}$ | SFR (UV)b [$M_{\odot}\>\mathrm{yr}^{-1}$]
---|---|---|---|---|---|---|---|---
NGC 1313 | 03:18:30.05 | -66:28:45.0 | SBd | 4.39 | 2.6$\times 10^{9}$ | 0.0965 | 0.008 | 1.15
NGC 1705 | 04:54:13.7 | -53:21:40.9 | SA0/BCG | 5.1 | 1.3$\times 10^{8}$ | 0.0071 | 0.008 | 0.11
NGC 4395 | 12:25:56.60 | 33:31:10.1 | SAm | 4.3 | 6.0$\times 10^{8}$ | 0.0152 | 0.004 | 0.34
NGC 4449 | 12:28:12.06 | 44:05:42.3 | IBm | 4.31 | 1.1$\times 10^{9}$ | 0.0171 | 0.004 | 0.94
NGC 7793 | 23:57:40.86 | -32:35:20.6 | SAd | 3.44 | 3.2$\times 10^{9}$ | 0.0171 | 0.008 | 0.52

a Galaxy Name, coordinates (J2000), morphological type and the foreground
color excess ($A_{V}/R_{V}$) due to the dust in our Milky Way, as listed in
the NASA Extragalactic Database (NED).

b Distance, stellar mass and SFR of each galaxy from Calzetti et al. (2015a).

c Metallicity of each galaxy, from Sabbi et al. (2018).

The cluster candidates in the LEGUS catalogs have been classified according to
four morphological classes as defined in Adamo et al. (2017). Class 1
candidates are symmetric and concentrated. Class 2 candidates look less
symmetric and rather elongated relative to class 1. Class 3 candidates are
asymmetric, less compact and have multiple peaks joined by diffuse emission in
between the peaks. Candidates in both class 1 and class 2 are considered YSCs
while class 3 candidates are “compact associations”. Class 4 candidates are
“non clusters”: artifacts, background galaxies, single stars or asterisms. In
our analysis, we only use the candidates of classes 1, 2, and 3 as real star
clusters; class 4 sources will not be considered any further.

Since we are trying to characterize the uIMF, we require our YSCs to have
undergone minimal evolution. We choose clusters with ages $\leq$ 4 Myr and
treat them as coeval objects that have not yet undergone supernova events or
strong gas expulsion (Weidner et al., 2010; Relaño et al., 2012). Leakage of
ionizing photons from H II regions can still occur, and can affect our
results, as we discuss below. We require cluster masses above
$500\>M_{\odot}$. This ensures that our YSCs are not dominated by mis-
classified single massive stars (Andrews et al., 2014), but are still probing
stochastic sampling effects. After all selections are applied, we have 375
clusters in total for this study (Table LABEL:tab:class). Figure 1 shows the
locations of the clusters identified by class 1, 2, and 3 in NGC 4449.

Table 2: Statistics of the LEGUS star clusters by class and by galaxy for age within 1–4 Myr, $E(B-V)\leq 0.2$, and mass over $500\>M_{\odot}$ used in this study. Name | Class 1 | Class 2 | Class 3 | Total
---|---|---|---|---
|  |  |  | (each galaxy)
NGC 1313 | 27 | 51 | 88 | 166
NGC 1705 | 0 | 1 | 1 | 2
NGC 4395 | 20 | 42 | 53 | 115
NGC 4449 | 6 | 5 | 2 | 13
NGC 7793 | 15 | 34 | 30 | 79
Total | 68 | 133 | 174 | 375

###  2.2 Photometry

In order to measure the H$\alpha$ luminosity of each star cluster, we retrieve
from the Hubble Legacy Archive the images of the five galaxies in the
F547M/F550M (V-short) and in the F657N/F658N (H$\alpha$+[N II]) filters
obtained with WFC3/ACS. We give preference for the medium band F550M or F547M
over the F555W filter to avoid the strong [O III] emission line which affects
the broader filter. We align the images to the F814W (I-band) images of the
LEGUS data for aperture photometry and additional corrections. Aperture
photometry is performed on all 375 clusters in all three bands using the
PHOTUTILS package in ASTROPY with a 5 pixel aperture radius, which corresponds
to 0.35–0.5 R${}_{\rm Str\ddot{o}mgren}$ (Calzetti et al., 2010) and using an
annulus between 7 to 9 pixel radius for local background. The positions of the
selected clusters are taken from the LEGUS catalogs.

The flux density in each band is obtained by converting counts to
erg$\cdot$cm${}^{-2}\cdot$s-1 after local background subtraction. Flux
densities in the H$\alpha$ in each cluster are measured from the difference
between the observed flux density in the H$\alpha$+[N II] filters
(F657N/F658N) and that of the underlying stellar continuum in the narrow band
obtained from the interpolation between the F547M/F550M and the F814W filters.
The flux conversion keyword PHOTFLAM from the images’ header is then used to
convert the counts to physical units. Finally, we convert the flux densities
to fluxes in H$\alpha$+[N II] by multiplying our measurements by the FWHMs of
the narrow-band filters. In the case of undetected H$\alpha$, the detection
limits of each cluster in the F657N/F658N filters are reported as upper
limits.

Aperture corrections are applied to all photometry measurements. The continuum
subtracted H$\alpha$+[N II] (F657M/F658M) images were used to select isolated
sources with strong emission to construct curves of growth. The final aperture
correction is derived from the average calculated for the isolated sources.
The typical correction is a factor $\sim 2.3$. As our sources are all
extremely young, the nebular emission, when present, is generally concentrated
and coincident with the location of the star cluster’s stellar emission. Thus,
we do not find large variations in the correction factor. Removal of the
contamination by [N II] in the F657M/F658M filters is performed by applying
the average ratios of [N II] to H$\alpha$ of each galaxy to the flux densities
in H$\alpha$+[N II] (Kennicutt, Jr. et al., 2008).

Extinction correction for $L_{\mathrm{H}\alpha}$ is performed considering two
contributions; the foreground Milky Way color excess,
$E(B-V)_{\mathrm{MilkyWay}}$, and the color excess intrinsic to the cluster,
$E(B-V)_{\mathrm{cluster}}$. The latter is determined from the SED fitting of
the cluster’s photometry (Adamo et al., 2017). For the intrinsic attenuation
of the nebular gas, we adopt the recipe of Calzetti et al. (1994) and Calzetti
et al. (2000), which establish:
$E(B-V)_{\mathrm{gas}}=E(B-V)_{\mathrm{cluster}}/0.44$. The attenuation
affecting the star cluster can then be expressed as:

$F_{\mathrm{int}}(\lambda)=F_{\mathrm{obs}}(\lambda)\>10^{\,0.4\,E(B-V)\,k(\lambda)}$
(3)

where $k(\lambda)=2.54$ for H$\alpha$ and
$E(B-V)=E(B-V)_{\mathrm{MilkyWay}}+E(B-V)_{\mathrm{cluster}}/0.44$. Figure 2
shows an example of the spatial distribution of the selected clusters for this
study in NGC 4449, indicating whether the H$\alpha$ emission is present or
absent, and assigning color-coding to their ages and masses.

Figure 2: The spatial distribution of the star clusters in NGC 4449 selected
for this study, and showing the presence/absence of H$\alpha$ (left), their
masses (middle) and their ages (right). The values of the masses
($10^{2}\mbox{--}10^{5}\>M_{\odot}$) and the ages (1–4 Myr) are shown by the
color bars in the middle and right panels, respectively. The gray dots are all
star clusters in NGC 4449 from the LEGUS catalogs.

##  3 Analysis and results

The number of massive stars in a cluster can be inferred from the hydrogen
ionizing photon rate, $Q(\mathrm{H}^{0})$, due to the strong ionizing emission
from young and hot OB stars. The H$\alpha$ luminosity, $L_{\mathrm{H}\alpha}$,
of a cluster is a tracer of $Q(\mathrm{H}^{0})$, such that we can test whether
or not the uIMF depends on cluster mass by normalizing $L_{\mathrm{H}\alpha}$
by the cluster mass ($M_{\mathrm{cl}}$) (Calzetti et al., 2010).

The distributions of the ratio $L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}$ for our
cluster samples are shown as histograms of the relative frequencies for
individual galaxy and for all samples combined in Figure 3; we plot separately
the H$\alpha$ detected clusters and the non-detected ones as upper limits. The
$\log(L_{\mathrm{H}\alpha}/M_{\mathrm{cl}})$ values cover the range 31–35 when
all 375 YSCs are considered. The ratio $L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}$
of the selected clusters in the five galaxies are shown individually in Figure
4. We can hardly see any clear trend between
$L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}$ and $M_{\mathrm{cl}}$, suggesting a
weak correlation between $m_{\mathrm{max}}$ and $M_{\mathrm{cl}}$.

Figure 3: The distributions of $L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}$ values
for each galaxy and for all samples combined. The blue bars are
H$\alpha$-emission detected sources and the red bars are undetected ones, for
which the indicated values are upper limits. For the sources detected in
H$\alpha$ the typical 2-sided 1$\sigma$ uncertainties are about 1/2 bin in
size. The SFR of each galaxy is given under the name of the galaxy.

The clusters with undetected H$\alpha$ emission represent 54% of the total
number of sources and are located $\sim 1.8$ orders of magnitude on average
below the sources with H$\alpha$ emission. Given that the majority of our star
clusters have masses below $\sim 3000\>M_{\odot}$ (Figure 4), the large number
of sources without H$\alpha$ emission is easily explainable as an effect of
the stochastic sampling of the IMF in this regime. Thus, Figure 3 supports
that the LEGUS sample selection is unbiased relative to presence or absence of
nebular emission. The gap between the two distributions is likely an effect of
our measurement technique coupled with evolutionary effects of the clusters.
For instance, very young clusters hosting massive stars can develop energetic
stellar winds over short timescales that remove the nebular gas surrounding
the cluster (MacKenty et al., 2000; Krumholz et al., 2019). Thus, we cannot
exclude that, at least for some of the sources without detected H$\alpha$
emission, the lack of detection is due to purely morphological effects because
the H$\alpha$ emission may actually be present but located in a shell around
the cluster (MacKenty et al., 2000; Calzetti et al., 2010; Relaño et al.,
2012; Andrews et al., 2013; Calzetti et al., 2015b; Belfiore et al., 2022). An
extreme example is the case of young, UV-bright starbursts with high SFRs and
many massive O stars but which lack of ionized gas due to ionizing photon
leakage (Marques-Chaves et al., 2022). In this case, our photometric
measurements provide a lower limit to the true nebular emission. We cannot
correct for this effect as most of our sources are located in crowded regions,
where their H$\alpha$ emission overlaps with that of neighboring sources and
shells of nebular emission cannot be easily identified and measured.

Conversely, we expect many of the low-mass sources undetected in H$\alpha$ to
be subject to stochastic sampling effects, meaning that their low mass
prevents them from forming massive stars in significant numbers. In these
cases, the sources cannot ionize the surrounding gas and our reported
H$\alpha$ values are upper limits. As we cannot discriminate, for the
H$\alpha$-undetected sources, between the upper limits from stochastic
sampling and the lower limits from expelled gas (see above), we will consider
both cases in the analysis that follows.

Before binning by mass to mitigate stochasticity, we check whether the
physical parameter we use, the cluster mass, may be influenced by the way the
SED fitting performed on the photometry of the clusters in the LEGUS catalogs
(Adamo et al., 2017). In other words, we want to make sure that the values of
the cluster masses are not biased by potential covariances in the SED fitting
results due to photometric uncertainties. In order to perform this assessment,
we carry out the same procedure as in Figure 4 but using purely observational
quantities, i.e., the I-band (F814W) flux density instead of
$M_{\mathrm{cl}}$. The choice of this band is driven by its low sensitivity to
extinction. Figure 5 shows the distribution of the clusters’
$L_{\mathrm{H}\alpha}/L_{\mathrm{I}}$ as a function of $L_{\mathrm{I}}$.
Figure 5 has, in the low-luminosity regime ($L_{\mathrm{I}}<38$), a comparable
scatter to Figure 4 in the low-mass regime; there is a gap between the
averages of $L_{\mathrm{H}\alpha}/L_{\mathrm{I}}$ for sources detected and
undetected in H$\alpha$, of comparable magnitude to the gap in Figure 4. The
errors in $L_{\mathrm{H}\alpha}/L_{\mathrm{I}}$ do not show a trend as a
function of $L_{\mathrm{I}}$, again similar to the data in Figure 4. This
suggests that the clusters’ $L_{\mathrm{H}\alpha}$ and their masses are not
strongly correlated with each other, and there is no bias in using the cluster
mass as our reference parameter.

Figure 4: Left: the ratio $L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}$ of all the
selected clusters as a function of stellar mass, as listed in the LEGUS
catalogs. Star markers are the clusters with H$\alpha$ emission, and dots are
those without H$\alpha$ emission. For these, the indicated values of
$L_{\mathrm{H}\alpha}$ are upper limits. Dashed horizontal lines are the
average of $L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}$ for the sources detected in
H$\alpha$ ($\langle L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle_{\bigstar}$)
and those undetected in the line ($\langle
L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle_{\bullet}$), respectively. Right:
the same plot as the left panel, but including the error bars of each cluster.
Here the upper limits are clearly marked.

We test the $m_{\mathrm{max}}\mbox{--}M_{\mathrm{ecl}}$ relation with the
approach suggested by Calzetti et al. (2010). The IMF and the IGIMF are
equivalent under an assumption of a universal IMF. The only effect at play is
stochastic sampling, which is predominant at low cluster masses (Hermanowicz
et al., 2013). A low-mass cluster will be missing stars at randomly sampled
masses; since, for a typical IMF functional shape, there are far more low-mass
stars than high-mass stars, random sampling will affect the latter more than
the former. In order to be truly stochastic, the incomplete sampling needs to
be erased when accumulating multiple low-mass clusters together. That is, the
shape of the IMF resulting from the sum of 100 individual clusters with each
mass of $1000\>M_{\odot}$ will be equivalent to the IMF of a single
$10^{5}\>M_{\odot}$ cluster. (Elmegreen, 2000, 2006; Bastian et al., 2010;
Calzetti et al., 2010). So, we sum all of the values of
$L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}$ as

$\langle\frac{L_{\mathrm{H}\alpha}}{M_{\mathrm{cl}}}\rangle=\frac{\Sigma_{i}\,L_{\mathrm{H}\alpha,i}}{\Sigma_{i}\,M_{\mathrm{cl},i}},$
(4)

which we divide into several mass bins of equal total mass each (Table
LABEL:tab:bin) to test any trend of the uIMF with mass (Andrews et al., 2013,
2014). We include in the sum above also clusters that are undetected in
H$\alpha$ to avoid biasing our results Weidner et al. (2014).

Figure 5: The ratio $L_{\mathrm{H}\alpha}/L_{\mathrm{I}}$ of the cluster
samples in this study as a function of $L_{\mathrm{I}}$. This figure has the
same structure with Figure 4 but $L_{\mathrm{I}}$ instead of
$M_{\mathrm{cl}}$. Figure 6: Left: the location of $\langle
L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$
($=\Sigma\,L_{\mathrm{H}\alpha}\,/\,\Sigma\,M_{\mathrm{cl}}$) with their
standard deviations (Table LABEL:tab:bin) as a function of the mean of the
samples’ mass in each bin. Different colors and shapes identify different
galaxies. The dark-gray squares are the $\langle
L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ ratios of all galaxies combined,
as a function of binned mass. The horizontal bar associated with each symbol
shows the size of each mass bin. The horizontal dotted lines show the
evolution of model $\langle L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ as a
function of age (1–6 Myr) with fixed $m_{\mathrm{max}}=120\>M_{\odot}$ from
SB99. Right: the same data points as the left panel shown with a range of
model upper mass values $m_{\mathrm{max}}$ averaged within the 1–4 Myr age
produced by SB99 from $20\>M_{\odot}$ to $120\>M_{\odot}$ in $10\>M_{\odot}$
intervals (horizontal dotted lines).

Figure 6 shows the $\langle L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$
values in bins of cluster mass, for each galaxy and for the entire sample
combined. Data are divided into several mass bins to check the
$m_{\mathrm{max}}\mbox{--}M_{\mathrm{ecl}}$ relation except for NGC 1705,
which has only one data point located at $M_{\mathrm{cl}}\sim
10^{3.4}\>M_{\odot}$ due to the small cluster catalog (only two YSCs, Table
LABEL:tab:class). The value of the maximum mass bin in NGC 7793 is $\sim 10$
times lower than those of other galaxies, so $\langle
L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ of NGC 7793 does not cover the
higher mass range ($M_{\mathrm{cl}}\gtrsim 10^{3.5}\>M_{\odot}$) of the YSCs
covered by the other four galaxies.

Table 3: Statistics on the binned data by mass in each individual galaxy and
in all galaxies combined.

Name & Bin No. |  $\langle M\rangle$ [$M_{\odot}$] | Mass range [$M_{\odot}$] | No. of clusters |  $M_{total}$ [$M_{\odot}$] | $\langle L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ |  $\sigma_{\langle L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle}$a
---|---|---|---|---|---|---
NGC 1313 |  |  |  |  |  |
Bin 1 | $1.235\times 10^{3}$ | $5.211\times 10^{2}-2.348\times 10^{3}$ | 127 | $1.569\times 10^{5}$ | 33.819 | 0.0256
Bin 2 | $4.685\times 10^{3}$ | $2.358\times 10^{3}-1.260\times 10^{4}$ | 33 | $1.546\times 10^{5}$ | 33.876 | 0.0687
Bin 3 | $2.658\times 10^{4}$ | $1.452\times 10^{4}-4.957\times 10^{4}$ | 6 | $1.595\times 10^{5}$ | 33.224 | 0.237
NGC 1705 |  |  |  |  |  |
Bin 1 | $2.342\times 10^{3}$ | $1.091\times 10^{3}-3.592\times 10^{3}$ | 2 | $4.683\times 10^{3}$ | 33.475 | 0.131
NGC 4395 |  |  |  |  |  |
Bin 1 | $8.448\times 10^{2}$ | $5.366\times 10^{2}-1.422\times 10^{3}$ | 11 | $9.292\times 10^{3}$ | 34.500 | 0.0486
Bin 2 | $2.544\times 10^{4}$ | $2.902\times 10^{3}-4.797\times 10^{4}$ | 2 | $5.087\times 10^{4}$ | 33.532 | 0.488
NGC 4449 |  |  |  |  |  |
Bin 1 | $1.193\times 10^{3}$ | $5.205\times 10^{2}-3.102\times 10^{3}$ | 64 | $7.637\times 10^{4}$ | 34.292 | 0.0388
Bin 2 | $5.911\times 10^{3}$ | $3.259\times 10^{3}-1.881\times 10^{4}$ | 13 | $7.685\times 10^{4}$ | 34.604 | 0.480
Bin 3 | $5.260\times 10^{4}$ | $1.952\times 10^{4}-8.567\times 10^{4}$ | 2 | $1.052\times 10^{5}$ | 33.839 | 0.103
NGC 7793 |  |  |  |  |  |
Bin 1 | $7.684\times 10^{2}$ | $5.172\times 10^{2}-1.181\times 10^{3}$ | 69 | $5.302\times 10^{4}$ | 33.834 | 0.0489
Bin 2 | $1.643\times 10^{3}$ | $1.193\times 10^{3}-2.445\times 10^{3}$ | 33 | $5.422\times 10^{4}$ | 34.123 | 0.0397
Bin 3 | $4.067\times 10^{3}$ | $2.471\times 10^{3}-1.197\times 10^{4}$ | 13 | $5.287\times 10^{4}$ | 34.396 | 0.0745
All galaxies |  |  |  |  |  |
Bin 1 | $1.087\times 10^{3}$ | $5.172\times 10^{2}-2.317\times 10^{3}$ | 293 | $3.186\times 10^{5}$ | 34.038 | 0.0160
Bin 2 | $4.368\times 10^{3}$ | $2.348\times 10^{3}-1.452\times 10^{4}$ | 73 | $3.189\times 10^{5}$ | 34.176 | 0.119
Bin 3 | $3.522\times 10^{4}$ | $1.881\times 10^{4}-8.567\times 10^{4}$ | 9 | $3.170\times 10^{5}$ | 33.900 | 0.0786

aPropagated in log10 scale.

In order to compare the observational uIMFs to theoretical ones, we generate
models using STARBURST99 (SB99; Leitherer et al. 1999) with a standard Kroupa
IMF up to mmax=120 M⊙, Padova stellar evolution tracks with AGB stars and
metallicity $Z=0.004$. The choice of metallicity for the models is consistent
with the metallicity of the sample galaxies (Table LABEL:tab:object). The data
points in Figure 6 are compared with cluster models with both fixed
$m_{\mathrm{max}}$ but different ages (left panel) and different
$m_{\mathrm{max}}$ averaged over the age range 1–4 Myr (right panel). After
distributing the clusters into 3 mass bins, the total equivalent masses in
each bin are $3.186\times 10^{5}\>M_{\odot}$, $3.189\times 10^{5}\>M_{\odot}$
and $3.170\times 10^{5}\>M_{\odot}$ for the mass bins centered at $1.087\times
10^{3}\>M_{\odot}$, $4.368\times 10^{3}\>M_{\odot}$ and $3.522\times
10^{4}\>M_{\odot}$, respectively. The mass range of each bin is $5.172\times
10^{2}\>M_{\odot}\,\mbox{--}\,2.317\times 10^{3}\>M_{\odot}$, $2.348\times
10^{3}\>M_{\odot}\,\mbox{--}\,1.452\times 10^{4}\>M_{\odot}$ and $1.881\times
10^{4}\>M_{\odot}\,\mbox{--}\,8.567\times 10^{4}\>M_{\odot}$, respectively.
More details on the binned data for each galaxy and for the combined sample
are given in Table LABEL:tab:bin. To test the impact of upper limits on our
binned values, we replace all upper limits with zero luminosity in H$\alpha$
($L_{\mathrm{H}\alpha}=0$). We find that that impact is negligible due to the
small contribution ($\lesssim 1.6\%$) of the upper limits to the mean value of
$L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}$ in Eq.(4) (Figure 4).

The distribution of $\langle L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ for
the combined sample shows a plateau as a function of mass rather than
following the increasing trend for increasing mass described by the
$m_{\mathrm{max}}\mbox{--}M_{\mathrm{ecl}}$ relation (Weidner et al., 2010).
Taking $\langle L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle\sim
10^{34}\,$erg$\cdot$s-1$\cdot$$M_{\odot}^{-1}$ as our fiducial value at high
cluster masses, the expected value at $M_{\mathrm{cl}}=10^{3}\>M_{\odot}$
($M_{\mathrm{ecl}}\sim 2\times 10^{3}\>M_{\odot}$) should be about 0.5 dex
lower according to the IGIMF models. We don’t observe such a decrease, which
is several standard deviations below our fiducial value at that mass (Figure
6). Furthermore, the $\langle L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$
values in the highest mass bins, $M_{\mathrm{cl}}\sim 10^{4.5}\>M_{\odot}$,
are lower than those at lower mass, similarly to what observed in NGC 4214
(Andrews et al., 2013).

The most likely reason for this decrease is that massive star clusters
statistically host massive stars more frequently than lower-mass clusters,
implying that feedback from those massive stars can effectively expel gas from
the cluster location, an effect magnified in clusters residing in low H I
density regions (Pellegrini et al., 2012). Gas expulsion implies that there is
less gas nearby to produce H$\alpha$. This is observed in the most massive
cluster in NGC 4214, where the presence of P-Cygni profiles for UV lines like
Si IV ($\lambda\sim$1400 Å) and C IV ($\lambda\sim$1550 Å) indicates young
ages, $\sim$4 Myr (Leitherer et al., 1996, 2002), but the cluster occupies a
cavity in ionized gas (MacKenty et al., 2000).

The horizontal dotted lines in the left panel (Figure 6) show the theoretical
$L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}$ as a function of cluster age. We expand
the range of ages to 6 Myr which is broader than in our sample selection,
because we include the uncertainties in the age determinations of the
clusters, estimated to be 1–2 Myr (Adamo et al., 2017). There is almost no
difference in the theoretical values of the ratios among 1–3 Myr located at
$L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\sim 34.8$. The
$L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}$ becomes lower when the age of the
cluster increases to 6 Myr ($L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\sim 33.7$).
Nearly all data points are above $L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}$ at 6
Myr except those of the highest mass bins in NGC 1313 and NGC 4395 where most
of the sources undetected in H$\alpha$ are located in shell-like or diffused H
II regions, possibly due to massive stellar feedback.

The right panel in Figure 6 shows the comparison of the observed $\langle
L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ and the theoretical values of
$L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}$ as a function of $m_{\mathrm{max}}$.
The theoretical uIMFs are calculated using SB99 with different
$m_{\mathrm{max}}$ values from $20\>M_{\odot}$ to $120\>M_{\odot}$ in
$10\>M_{\odot}$ steps. The H$\alpha$ luminosities are averaged within the 1–4
Myr age range for each $m_{\mathrm{max}}$ bin. The uIMFs of all data points
are distributed below $m_{\mathrm{max}}\sim\>70\>M_{\odot}$ at
$m_{\mathrm{max}}\sim 30\mbox{--}40\>M_{\odot}$. As already remarked, this may
point to the ages of our sample clusters being somewhat older than our assumed
$\sim$4 Myr, and consistent with ages as old as $\sim$6 Myr.

Although the distribution of $\langle
L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ is relatively flat as a function
of cluster mass in the mass range $\sim 500\mbox{--}50,000\>M_{\odot}$ (Figure
6), we find that the ionizing photon flux is significantly lower than expected
for young ($\leq 4$ Myr) star clusters populated with a universal Kroupa IMF
(Kroupa, 2001) with upper mass cut off heavier than $100\>M_{\odot}$. In other
words, the combined sample data (gray squares) are located in correspondence
of the $\langle L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ values expected
for populations with ages between 5 Myr and 6 Myr. In addition to the
uncertainties in the ages, as already discussed above, we consider here
ionizing photon leakage out of the star clusters. Although we attempt to
recover as much of the $L_{\mathrm{H}\alpha}$ as possible, we cannot measure
ionizing photons leakage directly. The escape fraction of ionizing photons
from H II regions is estimated to be around 40–60% by previous studies of
diffuse ionized gas (DIG) in nearby galaxies (Oey et al., 2007; Pellegrini et
al., 2012; Belfiore et al., 2022). We have some indirect evidence of this
process from the difference in the average
$L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}$ values for NGC 1313 and NGC 4449. The
average values for the clusters in NGC 4449 are a factor of $\sim$ 5
systematically higher than those of NGC 1313, for all mass bins. NGC 4449 is a
higher density galaxy than NGC 1313; comparing their SFR surface densities as
a proxy for gas density, previous studies show $\Sigma_{\mathrm{SFR}}\sim
0.04\>M_{\odot}\cdot\mathrm{yr}^{-1}\cdot\mathrm{kpc}^{-2}$ for NGC 4449
(Annibali et al., 2011) and $\Sigma_{\mathrm{SFR}}\sim
0.01\>M_{\odot}\cdot\mathrm{yr}^{-1}\cdot\mathrm{kpc}^{-2}$ for NGC 1313
(Silva-Villa & Larsen, 2011; Messa et al., 2021). Thus, to the extent that
lower density regions are also more porous (Pellegrini et al., 2012), our data
support significant ionizing photon leakage out of the H II regions
surrounding the star clusters in NGC 1313. This finding is in line with other
results, indicating that high-SFR regions, like, e.g., the central starburst
in the spiral M 83, may be less porous than low-SFR regions and contain
clusters that do not tend to leak ionizing photons (Della Bruna et al., 2022).

In the binned mass plot (Figures 6), NGC 7793 is the only galaxy with an
obvious trend of increasing $\langle
L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ for increasing mass, which may
support the $m_{\mathrm{max}}\mbox{--}M_{\mathrm{ecl}}$ model. However, the
other four galaxies do not display a comparable trend, and NGC 7793 does not
contain clusters with mass $>1.2\times 10^{4}\>M_{\odot}$, which leaves any
inference about trends inconclusive. Furthermore, in the lowest mass bin in
NGC 7793, 75% of the clusters have $\lesssim 10^{3}M_{\mathrm{cl}}$ and
$L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}<10^{33.4}$. This characteristic produces
a mass bin with a value of $\langle
L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ that is much lower for NGC 7793
than for other galaxies. As a counterpoint, only 23% of the clusters that make
the NGC 4449’s lowest-mass bin have $\lesssim 10^{3}M_{\mathrm{cl}}$, and this
flattens considerably the trend of $\langle
L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ with cluster mass for this
galaxy.

Finally, irrespective of the observed trend, there are numerous low-mass
clusters, $\sim 10^{3}\>M_{\odot}$ and lower, in all the sample galaxies but
NGC 1705, that have significant emission in H$\alpha$ (Figure 4), even
stronger emission than the sources at masses $\gtrsim 5000\>M_{\odot}$. The
presence of H$\alpha$ emitting low-mass clusters indicates that stochasticity
is a major contributor to the ionized gas emission of star clusters, further
supporting our results.

##  4 Summary and Conclusion

We have investigated the possibility of a non-universal upper-end for the
stellar IMF (uIMF) using catalogs of young star clusters. Star clusters
represent prime targets for this test due to the quasi-single-age nature of
their stellar populations.

We choose five nearby galaxies within $\sim 5$ Mpc with a broad range of SFRs
as calculated from the dust-corrected UV
($0.1\mbox{--}1.2\>M_{\odot}\>\mathrm{yr}^{-1}$) from the LEGUS data (Calzetti
et al., 2015a). The cluster catalogs of these galaxies (Adamo et al., 2017)
are used to test variations of uIMF and we limit our analysis to young stellar
clusters (YSCs) and compact associations with ages within 1–4 Myr and masses
over $500\>M_{\odot}$ (as determined from SED fits) to mitigate uncertainties
due to the evolution of massive stars out of the main sequence in older
clusters. We include sources without H$\alpha$ emission to prevent our results
from being affected by biases, due to the stochastic nature of the IMF
sampling, which disfavors presence of massive stars in low-mass clusters.

Our analysis includes photometry of the H$\alpha$ emission in correspondence
of all young star clusters identified in the five galaxies. The H$\alpha$
emission is obtained from imaging in the HST H$\alpha$+[N II] (F657N or F658N)
filters, after stellar continuum subtraction using the interpolated stellar
emission between V-short and I, as well as removal of the [N II] contamination
in the F657N/F658N filters. We apply aperture corrections to account for the
ionized gas outside of the photometric apertures by using growth curves
derived from isolated H$\alpha$-emitting sources.

Our results, shown in Figure 6, are consistent with those of Andrews et al.
(2013, 2014), in that the mass-normalized ionizing photon flux is basically
constant with cluster mass. This implies that the IMF from which the cluster
populations are drawn is universal. The $\langle
L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ values in Figure 6 do not
strongly follow the trend expected in the
$m_{\mathrm{max}}\mbox{--}M_{\mathrm{ecl}}$ relation of Weidner et al. (2010).
The three mass bins of the combined clusters are centered at $1\times
10^{3}\>M_{\odot}$, $4.4\times 10^{3}\>M_{\odot}$ and $3.5\times
10^{4}\>M_{\odot}$, shown as gray squares in Figure 6, and their $\langle
L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ values mark a relatively constant
trend, around $\sim 30\mbox{--}40\>M_{\odot}$ for the most massive stars. This
value can be explained if our sample, which is selected to include clusters up
to $\sim$ 4 Myr, actually includes clusters potentially as old as $\sim$ 6
Myr, due to uncertainties in the age determinations, as well as significant
ionizing photon leakage out of the H II regions. Leakage can lower the
measured H$\alpha$ luminosity by about a factor of 2, on average (Oey et al.,
2007; Pellegrini et al., 2012; Belfiore et al., 2022).

Looking at individual galaxies, only NGC 7793 shows a trend consistent with
the $m_{\mathrm{max}}\mbox{--}M_{\mathrm{ecl}}$ relation (Figures 6). However,
this galaxy lacks clusters more massive than $1.2\times 10^{4}\>M_{\odot}$,
which leaves the issue unsettled. In addition, several low-mass clusters in
NGC 7793, down to $\sim 10^{3}\>M_{\odot}$, show significant H$\alpha$
emission (Figure 4), implying the low-mass clusters host massive stars.

The $\langle L_{\mathrm{H}\alpha}/M_{\mathrm{cl}}\rangle$ values of the
heaviest mass bin at $\sim 10^{4.5}\>M_{\odot}$ have lower values than those
of other mass bins both for individual galaxies and for the combined sample.
This finding had been reported in Andrews et al. (2013, 2014) and is likely
caused by feedback effects from the faster dynamical evolution in higher mass
clusters.

We have confirmed that the uIMF is consistent with a universal trend, and
deviations are due to stochasticity, in good agreement with previous studies
(Calzetti et al., 2010; Andrews et al., 2013, 2014). More cluster samples,
especially from galaxies spanning a larger range of metallicities and SFRs
than those explored here, would be desirable to provide a definite test of the
dependence of the maximum mass stars on the cluster masses where they are
embedded.

Furthermore, infrared observations would allow us to obtain a census of the
clusters that are still deeply embedded in their natal clouds, and are missing
from our accounting. These tend to be the youngest among our clusters,
typically with ages $\leq 3$ Myr (Messa et al., 2021). These embedded clusters
are the most promising for sampling the youngest stars and the unexpanded gas
regions they ionize. Future observations with the James Webb Space Telescope
(JWST) are likely to provide such samples of embedded clusters even in more
distant galaxies then 5 Mpc, thanks to Webb’s better resolution in the
infrared than the Hubble Space Telescope.

We thank the LEGUS Survey team for providing us with the high-quality
observational data and catalogs to support this research.

## References

  * Adamo et al. (2017) Adamo, A., Ryon, J. E., Messa, M., et al. 2017, Astrophys. J., 841, 131, doi: 10.3847/1538-4357/aa7132
  * Andrews et al. (2013) Andrews, J. E., Calzetti, D., Chandar, R., et al. 2013, Astrophys. J., 767, 51, doi: 10.1088/0004-637X/767/1/51
  * Andrews et al. (2014) —. 2014, Astrophys. J., 793, 4, doi: 10.1088/0004-637X/793/1/4
  * Annibali et al. (2011) Annibali, F., Tosi, M., Aloisi, A., & Van Der Marel, R. P. 2011, Astron. J., 142, doi: 10.1088/0004-6256/142/4/129
  * Ascenso et al. (2009) Ascenso, J., Alves, J., & Lago, M. T. 2009, Astron. Astrophys., 495, 147, doi: 10.1051/0004-6361:200809886
  * Bastian et al. (2010) Bastian, N., Covey, K. R., & Meyer, M. R. 2010, Annu. Rev. Astron. Astrophys., 48, 339, doi: 10.1146/annurev-astro-082708-101642
  * Belfiore et al. (2022) Belfiore, F., Santoro, F., Groves, B., et al. 2022, Astron. Astrophys., 659, doi: 10.1051/0004-6361/202141859
  * Bertin & Arnouts (1996) Bertin, E., & Arnouts, S. 1996, Astron. Astrophys., 117, 393.  https://arxiv.org/abs/arXiv:1011.1669v3
  * Brands et al. (2022) Brands, S. A., De Koter, A., Bestenlehner, J. M., et al. 2022, Astron. Astrophys., 663, doi: 10.1051/0004-6361/202142742
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, Astrophys. J., 533, 682, doi: 10.1086/308692
  * Calzetti et al. (2010) Calzetti, D., Chandar, R., Lee, J. C., et al. 2010, Astrophys. J. Lett., 719, 158, doi: 10.1088/2041-8205/719/2/L158
  * Calzetti et al. (1994) Calzetti, D., Kinney, L. A., & Storchi-Bergmann, T. 1994, Astrophys. J., 429, 582
  * Calzetti et al. (2015a) Calzetti, D., Lee, J. C., Sabbi, E., et al. 2015a, Astron. J., 149, 51, doi: 10.1088/0004-6256/149/2/51
  * Calzetti et al. (2015b) Calzetti, D., Johnson, K. E., Adamo, A., et al. 2015b, Astrophys. J., 811, 75, doi: 10.1088/0004-637X/811/2/75
  * Cappellari et al. (2012) Cappellari, M., McDermid, R. M., Alatalo, K., et al. 2012, Nature, 484, 485, doi: 10.1038/nature10972
  * Cervino et al. (2002) Cervino, M., Valls-Gabaud, D., Luridiana, V., & Mas-Hesse, J. M. 2002, Astron. Astrophys., 381, 51, doi: 10.1051/0004-6361
  * Chabrier (2003) Chabrier, G. 2003, Publ. Astron. Soc. Pacific, 115, 763, doi: 10.1086/376392
  * Cook et al. (2019) Cook, D. O., Lee, J. C., Adamo, A., et al. 2019, Mon. Not. R. Astron. Soc., 484, 4897, doi: 10.1093/mnras/stz331
  * Corbelli et al. (2009) Corbelli, E., Verley, S., Elmegreen, B. G., & Giovanardi, C. 2009, Astron. Astrophys., 495, 479, doi: 10.1051/0004-6361:200811086
  * Crowther et al. (2010) Crowther, P. A., Schnurr, O., Hirschi, R., et al. 2010, Mon. Not. R. Astron. Soc., 408, 731, doi: 10.1111/j.1365-2966.2010.17167.x
  * Dabringhausen et al. (2008) Dabringhausen, J., Hilker, M., & Kroupa, P. 2008, Mon. Not. R. Astron. Soc., 386, 864, doi: 10.1111/j.1365-2966.2008.13065.x
  * Della Bruna et al. (2022) Della Bruna, L., Adamo, A., Mcleod, A. F., et al. 2022, Astron. Astrophys., 666, A29, doi: 10.1051/0004-6361/202243395
  * Elmegreen (1983) Elmegreen, B. G. 1983, Mon. Not. R. Astron. Soc., 203, 1011
  * Elmegreen (2000) —. 2000, ESA SP, 445, 265
  * Elmegreen (2006) —. 2006, Astrophys. J., 648, 572, doi: 10.1086/505785
  * Elmegreen & Scalo (2006) Elmegreen, B. G., & Scalo, J. 2006, Astrophys. J., 636, 149, doi: 10.1086/497889
  * Fardal et al. (2007) Fardal, M. A., Katz, N., Weinberg, D. H., & Davé, R. 2007, Mon. Not. R. Astron. Soc., 379, 985, doi: 10.1111/j.1365-2966.2007.11522.x
  * Fumagalli et al. (2011) Fumagalli, M., Da Silva, R. L., & Krumholz, M. R. 2011, Astrophys. J. Lett., 741, 1, doi: 10.1088/2041-8205/741/2/L26
  * Geha et al. (2013) Geha, M., Brown, T. M., Tumlinson, J., et al. 2013, Astrophys. J., 771, 29, doi: 10.1088/0004-637X/771/1/29
  * Grasha et al. (2015) Grasha, K., Calzetti, D., Adamo, A., et al. 2015, Astrophys. J., 815, 93, doi: 10.1088/0004-637X/815/2/93
  * Grasha et al. (2017a) —. 2017a, Astrophys. J., 840, 113, doi: 10.3847/1538-4357/aa6f15
  * Grasha et al. (2017b) Grasha, K., Elmegreen, B. G., Calzetti, D., et al. 2017b, Astrophys. J., 842, 25, doi: 10.3847/1538-4357/aa740b
  * Grasha et al. (2018) Grasha, K., Calzetti, D., Bittle, L., et al. 2018, Mon. Not. R. Astron. Soc., 481, 1016, doi: 10.1093/MNRAS/STY2154
  * Grasha et al. (2019) Grasha, K., Calzetti, D., Adamo, A., et al. 2019, Mon. Not. R. Astron. Soc., 483, 4707, doi: 10.1093/mnras/sty3424
  * Hermanowicz et al. (2013) Hermanowicz, M. T., Kennicutt, R. C., & Eldridge, J. J. 2013, Mon. Not. R. Astron. Soc., 432, 3097, doi: 10.1093/mnras/stt665
  * Jeřábková et al. (2018) Jeřábková, T., Hasani Zonoozi, A., Kroupa, P., et al. 2018, Astron. Astrophys., 620, 1, doi: 10.1051/0004-6361/201833055
  * Kalari et al. (2022) Kalari, V. M., Horch, E. P., Salinas, R., et al. 2022, Astrophys. J., 935, 162, doi: 10.3847/1538-4357/ac8424
  * Kennicutt, Jr. (1998) Kennicutt, Jr., R. C. 1998, Astrophys. J., 498, 541, doi: 10.1086/305588
  * Kennicutt, Jr. et al. (2008) Kennicutt, Jr., R. C., Lee, J. C., Funes, S. J., J. G., Sakai, S., & Akiyama, S. 2008, Astrophys. J. Suppl. Ser., 178, 247, doi: 10.1086/590058
  * Kroupa (2001) Kroupa, P. 2001, Mon. Not. R. Astron. Soc., 322, 231, doi: 10.1046/j.1365-8711.2001.04022.x
  * Kroupa & Boily (2002) Kroupa, P., & Boily, C. M. 2002, Mon. Not. R. Astron. Soc., 336, 1188, doi: 10.1046/j.1365-8711.2002.05848.x
  * Kroupa & Weidner (2003) Kroupa, P., & Weidner, C. 2003, Astrophys. J., 598, 1076, doi: 10.1086/379105
  * Kroupa et al. (2012) Kroupa, P., Weidner, C., Pflamm-Altenburg, J., et al. 2012, Planets, Stars Stellar Syst. Vol. 5 Galact. Struct. Stellar Popul., 5, 115, doi: 10.1007/978-94-007-5612-0_4
  * Krumholz et al. (2019) Krumholz, M. R., McKee, C. F., & Bland-Hawthorn, J. 2019, Annu. Rev. Astron. Astrophys., 57, 227, doi: 10.1146/annurev-astro-091918-104430
  * Lada & Lada (2003) Lada, C. J., & Lada, E. A. 2003, Annu. Rev. Astron. Astrophys., 41, 57, doi: 10.1146/annurev.astro.41.011802.094844
  * Larson (1982) Larson, R. B. 1982, Mon. Not. R. Astron. Soc., 200, 159.  http://eprints.uanl.mx/5481/1/1020149995.PDF
  * Leitherer et al. (2002) Leitherer, C., Calzetti, D., & Martins, L. P. 2002, Astrophys. J., 574, 114, doi: 10.1086/340902
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 1, 3
  * Leitherer et al. (1996) Leitherer, C., Vacca, W. D., Conti, P. S., et al. 1996, Astrophys. J., 465, 717\.  https://arxiv.org/abs/تم
  * MacKenty et al. (2000) MacKenty, J. W., Maiz-Apellaniz, J. S., Pickens, C. E., Norman, C. A., & Walborn, N. R. 2000, Astron. J., 2, 894
  * Maíz-Apellniz (2008) Maíz-Apellniz, J. 2008, Astrophys. J., 677, 1278, doi: 10.1088/0004-637X/699/2/1938
  * Marques-Chaves et al. (2022) Marques-Chaves, R., Schaerer, D., Alvarez-Marquez, J., et al. 2022, Mon. Not. R. Astron. Soc., 517, 2972.  https://arxiv.org/abs/2210.02392
  * Massey (2003) Massey, P. 2003, Annu. Rev. Astron. Astrophys., 41, 15, doi: 10.1146/annurev.astro.41.071601.170033
  * Messa et al. (2018a) Messa, M., Adamo, A., Ostlin, G., et al. 2018a, Mon. Not. R. Astron. Soc., 473, 996, doi: 10.1093/mnras/stx2403
  * Messa et al. (2018b) Messa, M., Adamo, A., Calzetti, D., et al. 2018b, Mon. Not. R. Astron. Soc., 477, 1683, doi: 10.1093/mnras/sty577
  * Messa et al. (2021) Messa, M., Calzetti, D., Adamo, A., et al. 2021, Astrophys. J., 909, 121, doi: 10.3847/1538-4357/abe0b5
  * Miller & Scalo (1979) Miller, G. E., & Scalo, J. M. 1979, Astrophys. J. Suppl. Ser., 41, 513, doi: 10.1086/190629
  * Oey (2011) Oey, M. S. 2011, Astrophys. J. Lett., 739, L46, doi: 10.1088/2041-8205/739/2/L46
  * Oey & Clarke (1998) Oey, M. S., & Clarke, C. J. 1998, Astron. J., 115, 1543.  http://iopscience.iop.org/1538-3881/115/4/1543
  * Oey et al. (2007) Oey, M. S., Meurer, G. R., Yelda, S., et al. 2007, Astrophys. J., 661, 801, doi: 10.1086/517867
  * Orozco-Duarte et al. (2022) Orozco-Duarte, R., Wofford, A., Vidal-García, A., et al. 2022, Mon. Not. R. Astron. Soc., 509, 522, doi: 10.1093/mnras/stab2988
  * Pellegrini et al. (2012) Pellegrini, E. W., Oey, M. S., Winkler, P. F., et al. 2012, Astrophys. J., 755, 40, doi: 10.1088/0004-637X/755/1/40
  * Pflamm-Altenburg et al. (2009) Pflamm-Altenburg, J., Weidner, C., & Kroupa, P. 2009, Mon. Not. R. Astron. Soc., 395, 394, doi: 10.1111/j.1365-2966.2009.14522.x
  * Relaño et al. (2012) Relaño, M., Kennicutt, R. C., Eldridge, J. J., Lee, J. C., & Verley, S. 2012, Mon. Not. R. Astron. Soc., 423, 2933, doi: 10.1111/j.1365-2966.2012.21107.x
  * Romano et al. (2005) Romano, D., Karakas, A. I., Tosi, M., & Matteucci, F. 2005, Astron. Astrophys., 430, 491, doi: 10.1051/0004-6361/201014483
  * Sabbi et al. (2018) Sabbi, E., Calzetti, D., Ubeda, L., et al. 2018, Astrophys. J. Suppl. Ser., 235, 23, doi: 10.3847/1538-4365/aaa8e5
  * Salpeter (1955) Salpeter, E. E. 1955, Astrophys. J., 121, 161.  https://ui.adsabs.harvard.edu/abs/1955ApJ...121..161S/abstract
  * Schneider et al. (2018) Schneider, F. R. N., Sana, H., Evans, C. J., et al. 2018, Science (80-. )., 361, 69, doi: 10.1126/science.aat6506
  * Silva-Villa & Larsen (2011) Silva-Villa, E., & Larsen, S. S. 2011, Astron. Astrophys., 529, doi: 10.1051/0004-6361/201016206
  * The Astropy Collaboration et al. (2022) The Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, Astrophys. J., 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Tolstoy et al. (2009) Tolstoy, E., Hill, V., & Tosi, M. 2009, Annu. Rev. Astron. Astrophys., 47, 371, doi: 10.1146/annurev-astro-082708-101650
  * Van Dokkum & Conroy (2011) Van Dokkum, P. G., & Conroy, C. 2011, Astrophys. J. Lett., 735, L13, doi: 10.1088/2041-8205/735/1/L13
  * Weidner et al. (2010) Weidner, C., Kroupa, P., & Bonnell, I. A. 2010, Mon. Not. R. Astron. Soc., 401, 275, doi: 10.1111/j.1365-2966.2009.15633.x
  * Weidner et al. (2014) Weidner, C., Kroupa, P., & Pflamm-Altenburg, J. 2014, Mon. Not. R. Astron. Soc., 441, 3348, doi: 10.1093/mnras/stu640
  * Weisz et al. (2012) Weisz, D. R., Johnson, B. D., Johnson, L. C., et al. 2012, Astrophys. J., 744, 44, doi: 10.1088/0004-637X/744/1/44
  * Weisz et al. (2015) Weisz, D. R., Johnson, L. C., Foreman-Mackey, D., et al. 2015, Astrophys. J., 806, 198, doi: 10.1088/0004-637X/806/2/198
  * Wilkins et al. (2008) Wilkins, S. M., Hopkins, A. M., Trentham, N., & Tojeiro, R. 2008, Mon. Not. R. Astron. Soc., 391, 363, doi: 10.1111/j.1365-2966.2008.13890.x
  * Yan et al. (2017) Yan, Z., Jerabkova, T., & Kroupa, P. 2017, Astron. Astrophys., 607, 1, doi: 10.1051/0004-6361/201730987