# Comprehending Linguistic and Emotional Cues in Demographically Diverse
Spatial Social Media Data

Salim Sazzed
Department of Computer Science,Old Dominion University, Norfolk, VA USA
Department of Computer Science,University of Memphis, Memphis, TN, USA
<EMAIL_ADDRESS>
# Comprehending Lexical and Affective Ontologies in the Demographically
Diverse Spatial Social Media Discourse

Salim Sazzed
Department of Computer Science,Old Dominion University, Norfolk, VA USA
Department of Computer Science,University of Memphis, Memphis, TN, USA
<EMAIL_ADDRESS>
###### Abstract

This study aims to comprehend linguistic and socio-demographic features,
encompassing English language styles, conveyed sentiments, and lexical
diversity within spatial online social media review data. To this end, we
undertake a case study that scrutinizes reviews composed by two distinct and
demographically diverse groups. Our analysis entails the extraction and
examination of various statistical, grammatical, and sentimental features from
these two groups. Subsequently, we leverage these features with machine
learning (ML) classifiers to discern their potential in effectively
differentiating between the groups. Our investigation unveils substantial
disparities in certain linguistic attributes between the two groups. When
integrated into ML classifiers, these attributes exhibit a marked efficacy in
distinguishing the groups, yielding a macro F1 score of approximately 0.85.
Furthermore, we conduct a comparative evaluation of these linguistic features
with word n-gram-based lexical features in discerning demographically diverse
review data. As expected, the n-gram lexical features, coupled with fine-tuned
transformer-based models, show superior performance, attaining accuracies
surpassing 95% and macro F1 scores exceeding 0.96. Our meticulous analysis and
comprehensive evaluations substantiate the efficacy of linguistic and
sentimental features in effectively discerning demographically diverse review
data. The findings of this study provide valuable guidelines for future
research endeavors concerning the analysis of demographic patterns in textual
content across various social media platforms.

##  I Introduction

Demographic data concerning user traits enable an understanding of user
behaviors, ultimately facilitating improved decision-making for various social
and business challenges. For instance, the analysis of demographic data
assists businesses in making informed decisions regarding marketing, product
development, customer experiences, and competitive positioning. In a social
context, diverse demographic review data empowers policymakers to gain
insights into the experiences and perspectives of different social groups,
which, in turn, facilitates decision-making that promotes equity, inclusive
representation, targeted interventions, evidence-based decision-making, and
accountability. Furthermore, demographically tagged social media plays a
significant role in computational social science by highlighting differences
in beliefs and behaviors among demographic groups [1].

Recent studies have shown that demographic traits can be inferred from the
linguistic characteristics of written content [2, 3, 4, 5]. The determination
of various demographic attributes of users, such as age or gender, from
written comments, has been investigated by [2, 6]. For example, Rosenthal and
McKeown [7] attempted to predict the users’ ages from blog content by
incorporating various features specific to the blog and the behavior and
interest of users. Schler et al. [8] observed significant differences in
content and style levels between male and female bloggers by analyzing a large
corpus of blogs of around 300 million words.

In addition, some studies tried to identify the English language nativeness of
the writers from demographically diverse data. Although the perspective of
their study was the second language acquisition (SLA) research, such as
contrastive analysis, syntactic or grammatical errors made by non-native
speakers [9, 10] based on corpus compiled from the sample essay of ESL
(English as a Second Language) learners such as TOEFL (Test of English as a
foreign language) [11], the international corpus of learner English [12].
Research pertaining to demographically diverse informal reviews, such as those
found in social media and prominently affected by English language nativeness
and fluency level, remains mostly unexplored [13] (except a few [14, 15, 16]).

Therefore, in this study, we aim to understand how linguistic and semantic
attributes of text vary across demographically different groups in the context
of social media, taking into account factors such as English language
nativeness, geography, and socio-culture. In particular, we aim to provide
insight into the following research questions-

  * •

RQ1: Do the variations of the linguistic features (e.g., synthetic, lexical)
render sufficient signals to distinguish diverse demographic groups when
incorporated in classical ML classifiers?

  * •

RQ2: Whether the linguistic or n-gram lexical features perform better for the
demography prediction task when incorporated into ML classifiers.

The two distinct demographic groups considered here represent individuals of
two different socio-economic cultures. More importantly, the English reviews
written by these two groups differ significantly in terms of English language
nativeness and fluency levels. The first review group represents restaurant
review data collected from Bangladesh, a country with mostly low proficient
non-native English speakers 111https://www.ef.com/wwen/epi/. In Bangladesh,
almost 98% of people are Bengali native speakers
222https://www.worldatlas.com/articles/what-languages-are-spoken-in-
bangladesh.html. The other group contains restaurant reviews written by users
located in the USA (mostly English native speakers). We extract various
statistical and syntactic features such as review length, frequency of opinion
words, and usage of POS (part-of-speech) from the reviews of both groups. We
find that linguistic features exhibit sufficient distinguishing signals to
differentiate between the two groups of reviews when used as input for
classical machine learning (CML) classifiers. Additionally, we explore the
performance of lexical word n-grams-based features for classifying review
groups by incorporating them into classical ML classifiers and transformer-
based language models. As expected, we observe that transformer-based models,
when utilizing lexical features, outperform both types of feature-based
classical ML classifiers.

###  I-A Contributions

The main contributions of this study can be summarized as follows-

  1. 1.

We identify differences in stylistic, syntactic, and statistical features
among reviews from diverse demographic groups.

  2. 2.

We show that it is possible to distinguish demographically diverse informal
social media reviews by incorporating various linguistic features into the
machine learning classifiers.

  3. 3.

Finally, we compare the differences in the performance of linguistic and
lexical n-gram-based features for distinguishing diverse review groups.

TABLE I: Sample reviews from Demography-1 and Demography-2
---
Figure 1: Overall framework of the study

##  II Demography Prediction Task

###  II-A Dataset

As mentioned earlier, in this study, we investigate restaurant reviews written
in English from two different demographic groups. The first group, which we
refer to as Demography-1, consists of reviews collected from restaurants
located in Bangladesh 333https://www.kaggle.com/tuxboy/restaurant-reviews-in-
dhaka-bangladesh. In Bangladesh, the majority of people (98%) speak standard
Bengali or one of its many dialects as their first language, while the
remaining population speaks regional or minority languages. English
proficiency in Bangladesh is generally categorized as low among secondary
speakers 444https://www.ef.com/wwen/epi/. The second review group, referred to
as Demography-2, comprises a subset of the Yelp restaurant review dataset
written by primarily English native speakers residing in the USA.

To mitigate potential domain bias in the evaluation process, the textual
content used in both groups was sourced exclusively from the same domain,
namely, restaurant reviews. Additionally, an equal number of reviews were
selected for each group to prevent the influence of class imbalance on the
classification outcomes. The final dataset comprises a total of 9974 reviews,
where each group includes 4987 reviews.

###  II-B Classical Machine Learning-based Prediction

To distinguish the two demographic groups, we employ several classical ML
classifiers: Logistic Regression (LR), Support Vector Machine (SVM), Random
Forest (RF), and Extra Trees (ET). For all classical ML classifiers, the
default parameter settings of the scikit-learn library [17] are used.

As an input of classical ML classifiers, we utilize two distinct groups of
features: i) linguistic features and ii) word unigram and bigram-based lexical
features, separately. The linguistic features further encompass several sub-
types, including grammatical and sentiment features.

####  II-B1 Linguistic Features

We extract three different types of linguistic features from the reviews,
which are utilized as input for the classical ML classifiers. Text statistics
such as review length in terms of words and sentences have been employed in
related work such as language variety identification task [18]. Grammatical
features such as the usage of PoS tags have been studied for language
nativeness identification tasks in the earlier works as these stylistic
features reflect user communication behavior and interaction style [19]. In
addition, lexicon coverage, a feature that refers to the usage of opinion or
sentiment words, is considered.

We employ the Mann-Whitney U test to determine which linguistic features show
significant differences in the two groups of reviews. The Mann-Whitney U test
is a non-parametric test of the null hypothesis, which is often used to test
the differences in the distributions of two sets of values. We utilize the
Mann-Whitney U test between two groups for all the linguistic features. The
null hypothesis states that the distribution of a specific attribute in
Demography-1 is the same as the underlying distribution of the same attribute
in Demography-2, while the alternative hypothesis suggests the opposite.

TABLE II: The statistics of various linguistic features in the reviews belong
to two groups, Demography-1 and Demography-2

Type | Feature | Demo-1 | Demo-2
---|---|---|---
| Total words (in corpus) | 147401 | 656672
Text | Total sentences (in corpus) | 16272 | 47922
Statistics | Mean review length (#words) | 29.56 | 131.70
| Mean review length (#sentences) | 3.32 | 9.611
| Mean sentence length (#words) | 8.89 | 13.70
Lexical | Total unique words (in corpus) | 10473 | 28297
Diversity |  |  |
| Total negation words (in corpus) | 1585 | 5287
Negation | (%) of negation words ( in corpus ) | 1.0% | 0.80%
| Negation words per review (mean) | 0.32 | 1.06
| Total articles (in corpus) | 6553 | 45918
| (%) of article in corpus (word) | 4.44% | 6.99%
| Articles/review (words) | 1.31 | 9.20
Grammatical | Total adjectives (in corpus) | 15632 | 57561
Feature | (%) of the adjective (in corpus) | 10.60% | 8.76%
| Number of adjectives/review | 3.13 | 11.54
| Total verbs (in corpus) | 14754 | 74881
| (%) of the verb in corpus | 10.01% | 11.40%
| Number of verbs/review | 2.95 | 15.01
| Total prepositions (in corpus) | 8466 | 52637
| (%) of the preposition (in corpus) | 5.79% | 8.00%
| Number of prepositions/review | 1.69 | 10.55
| Total SC(in corpus) | 2960 | 19402
| (%)of the SC(in corpus) | 2.0% | 2.95%
| SC per reviews | 0.59 | 3.89
| Total opinion words (Hu &Liu) | 8799 | 33691
Sentiment | Lexicon coverage (Hu & Liu) | 6.22% | 5.14%
Lexicon | Total opinion words (VADER) | 9516 | 36390
| Lexicon coverage (VADER) | 6.51% | 5.56%

#####  Text Length Features

These features represent the length of the reviews based on word and sentence
level: i) Number of words per review, ii) Number of sentences per review, iii)
Number of words per sentence. As earlier studies suggested that English texts
written by non-native English speakers are usually simpler than those of
natives [20], the sentence length could be a distinguishing factor for the
demography prediction task.

##### Grammatical and Negation Features

We consider a set of grammatical attributes that may provide signals to
discern the diverse demographic groups.

  * •

Articles: The number of articles (i.e., a, an, the) present in a review is
computed.

  * •

Adjectives: The number of adjectives present in a review is computed. The
spaCy [21] library is used to identify adjectives in a text.

  * •

Verbs: The usage of verbs in both corpora is provided. Similar to adjective
identification, spaCy [21] library is used for verb identification.

  * •

Prepositions: We calculate the number of prepositions present in the reviews
of two groups. A list of commonly used prepositions is considered (details can
be found here 555https://github.com/sazzadcsedu/LinguisticAnalysis).

  * •

Subordinating conjunctions (SC): Additionally, we take into account the
presence of subordinating conjunctions that indicate complex sentences. A
complex sentence typically consists of one or more dependent (subordinate)
clauses and one or more independent clauses. Subordinating conjunctions are
words or phrases that connect dependent clauses to independent clauses.
Examples of subordinating conjunctions include ”although,” ”as,” ”because,”
”before,” ”how,” ”if,” ”once,” ”since,” and so on. We examined the occurrence
of 50 commonly used subordinating conjunctions in each review
666https://github.com/sazzadcsedu/LinguisticAnalysis.

  * •

Negative words: The VADER [22] negative word list is used as a reference to
find the number of negative words in each group
777https://github.com/cjhutto/vaderSentiment .

##### Sentiment Lexicon Coverage

The coverage of two popular English sentiment lexicons, Opinion Lexicon [23]
and VADER [22], is computed for each of the reviews from both groups.

####  II-B2 Lexical n-gram Features

In addition, as lexical features, we extract word n-grams from the reviews of
both groups. An n-gram denotes a consecutive sequence of n items within a
given textual context. Particularly, we extract both unigrams (individual
words) and bigrams (two-word combinations) from the review texts. Following
this extraction, we calculate their respective term frequency-inverse document
frequency (tf-idf) scores and incorporate these computed scores as input
features into the CML classifiers.

Figure 2: Performances of classical ML classifiers using various types of
linguistic features, (A) text statistics features, (B) grammatical features,
(C) sentiment features, and (D) using features A, B, and C altogether Figure
3: Performances of classical ML classifiers using lexical features (i.e., word
unigrams and bigrams)

###  II-C Deep Learning-based Classification

The transformer-based pre-trained language models, such as BERT [24] and
RoBERTa [25], have shown state-of-the-art results in various text
classification tasks with limited labeled data. We fine-tune both transformer-
based pre-trained models for categorizing reviews into two demographically
diverse groups. Since this is a binary-level classification task, we utilize
the classification module of these pre-trained models. The Hugging face
library [26] is used for fine-tuning all the pre-trained models. As the
initial layers of pre-trained models primarily learn general features, they
are left unchanged during fine-tuning process. Only the last layer of the pre-
trained model is trained using new data specifically for the binary-level
classification task.

We tokenize the input data for fine-tuning the language model. As pre-trained
models typically support texts with a maximum of 512 tokens, we divide reviews
longer than 512 words into 512-word chunks. During training, all the 512-token
chunks are assigned the same class as the original review. During testing, the
final class of the review is determined by majority voting. In the event of a
tie, we consider the word length of the chunks to decide the final class
label. For training, we use a mini-batch size of 8 and a learning rate of
4*10-5. During the training process, 20% of the samples are dedicated to
validation. We optimize the pre-trained models using the Adam optimizer, with
the loss parameter set to categorical cross-entropy. The training procedure is
carried out for four epochs with an early stopping criterion set. Note that
all hyperparameters of deep learning models are determined based on empirical
evaluation.

Figure 4: The comparisons between the best results achieved by classical ML
classifiers (i.e., RF with linguistic features and SVM with lexical features),
and two transformer-based language models, for the demographic group
prediction task

##  III Results and Discussion

###  III-A Evaluation Settings

To compare different classification methods for the demography prediction
task, we utilize 5-fold cross-validation. Several key metrics, namely
precision, recall, macro F1 score, and accuracy, are employed to assess the
performance of classifiers. These metrics provide a comprehensive assessment
of the performance of the different classifiers.

TABLE III: Top adjectives(#occurrences) and verbs (#occurrences)

Top adjectives | Top verbs
---|---
Demography-1 | Demography-2 | Demography-1 | Demography-2
good | 1609 | good | 2854 | had | 1609 | had | 2987
bad | 407 | great | 1925 | have | 363 | have | 2489
worst | 331 | other | 1073 | go | 361 | get | 1818
best | 327 | little | 944 | ordered | 272 | go | 1609
great | 274 | nice | 884 | serve | 242 | got | 888
worst | 264 | more | 822 | went | 218 | ordered | 848
nice | 250 | best | 802 | love | 198 | know | 829
poor | 232 | few | 654 | served | 196 | make | 767
awesome | 214 | friendly | 645 | visit | 160 | going | 743

###  III-B Performance of classical ML Classifiers using Linguistic and
Lexical Features

The linguistic analysis reveals that several attributes of reviews, such as
review length, range of vocabulary used, coverage of opinion lexicon, and
usage of part-of-speech (POS) vary across two demographic groups that can be
attributed to the reviewer’s English language nativeness or proficiency level
and socio-culture (shown in Table II). It is observed that the usage of common
opinion terms is more apparent (i.e., high coverage of sentiment lexicon) in
the writing of non-native speakers (Demography-1). Figure 2 provides the
precision, recall, F1 score, and accuracy of various classical ML-based
methods utilizing linguistic features. We can see that when all the linguistic
features (i.e., length, grammatical, and sentiment) are utilized, all of the
four classical ML classifiers LR, SVM, RF, and ET perform similarly; They
achieve F1 scores of around 0.833 and accuracy around 83%.

We apply the Mann-Whitney test to find whether any of the stylistic features
are significantly different in the reviews of Demography-1 and Demography-2. A
$p$-value of 0.05 is used for the significance test. The Mann-Whitney test
indicates many of the linguistic features are significantly different in the
two groups; however, we observe they do not provide very high distinguishing
power when incorporated into the classical ML classifiers. For example, the
text statistic features, such as review length and sentence length, are
significantly different ($p$-value less than 0.05) in both groups; however,
when incorporated into the classical ML, they obtain a much lower F1 score of
0.813 than the lexical features. The high standard deviations (std.) of the
review length, with respect to both word and sentence, indicate that many
reviews in each group spread far away from its group mean value, which may
induce the classifier to yield wrong predictions.

Classical ML classifiers yield F1 scores between 0.89 and 0.94 when word
unigram and bigram lexical features are utilized (see Fig. 3). The better
performance of ML classifiers with lexical features indicates the presence of
some distinguishing socio-culture-specific words and named entities in the
reviews of both groups, which generate effective signals to discern the two
groups. Nevertheless, we find that reviews of both groups share some common
adjectives, such as good, great, nice, and best (Table III), which are among
the most frequently occurring adjectives in both groups. A few other
adjectives are also common to both groups.

This observation suggests that both demographic groups, irrespective of
English language nativeness, tend to use simple and commonly used adjectives
to express opinions and feelings. One dissimilarity we observe is that in
reviews of Demography-2 (written mostly by native speakers), adjectives of
quantity such as more, little, and few are more frequent than in the
Demography-1 review group. None of these adjectives of quantity appear among
the top 10 adjectives in Demography-1, and none of them occur more than 214
times. When examining the verbs used in the two groups of reviews, we observe
that the most frequent verbs are also very similar (see Table III). Although
there is a high presence of overlapping words in both groups, there are
certainly distinguishing words, such as named entities, that help the
classifiers discern between the two groups.

###  III-C Performance of Deep Learning based Classifiers

The fine-tuned transformer-based language models yield impressive results
utilizing unigram and bigram-based lexical features; they attain almost
perfect accuracy by correctly classifying around 97% instances (Fig 4). The
pre-trained language models are generated based on an enormous amount of
textual content, which helps to capture the implicit pattern of the reviews
and can effectively identify the language nativeness of reviewers. Also, the
high efficacy of lexical feature-based classical ML classifiers and
transformer-based models can be attributed to the English language proficiency
level of the people of Bangladesh, who are not known as very fluent English
speakers. Additionally, the presence of named entities with specific meanings
and socio-cultural terms in the reviews of both groups supports the lexical
n-gram-based approach, leading to better performance.

###  III-D Implications of the Study

This study provides insights into various perspectives by analyzing data from
two demographic groups, including the followings-

##### Language Landscapes and Identity

The study contributes to understanding how English language variation in
social media usage is linked to different demographic groups. It can shed
light on how individuals from diverse demographic backgrounds use language to
express their cultural and social identities, providing insights into the
complex relationship between language, ethnicity, and nationality.

##### Lexical Diversity in Online Communication

Examining the lexical diversity in online social media review data can shed
light on the richness and variety of language used by individuals. This
analysis can help researchers understand the level of vocabulary
sophistication, linguistic creativity, or the influence of cultural factors in
online communication.

##### Sociolinguistic Research and Language Policy

The study’s findings can contribute to sociolinguistic research and language
policy discussions. Understanding the relationship between language
nativeness, demographic diversity, and social media usage can inform
discussions on language rights, language maintenance, and linguistic identity
in digital spaces.

##  IV Summary, Limitations and Future Work

This study aims to distinguish reviews of two different socio-demographic
groups leveraging various linguistic and lexical features, language models,
and ML classifiers. From two demographically distinct review groups
(Demography-1 and Demography-2), various linguistic features are extracted to
train ML classifiers for the demography prediction task. In addition, two
state-of-the-art pre-trained language models, BERT and RoBERTa, are fine-tuned
with the n-gram-based features for the prediction task. We observe that
linguistic features are capable of distinguishing demographically diverse
reviews; when they are fed into classical ML classifiers, an F1 score (best
result) close to 0.85 is obtained. The pre-trained models exhibit very high
efficacy for distinguishing reviews using lexical features, which can be
attributed to the presence of name-entity and sociocultural-specific features
in the two review groups. Our analysis reveals the contrast and similarity of
implicit characteristics of reviews written by two demographically diverse
review groups. We present multiple approaches for the demography prediction
task that can help a diverse set of downstream decision-making tasks.

One of the limitations of this work is that it only considers two
demographically distinct groups. Since English proficiency levels (a
predominant factor affecting linguistic characteristics) among non-native
speakers may vary across demographics, such as geography, cultures, and
language families, it is worthwhile to analyze non-native English review data
from multiple demographics. Furthermore, it is worth exploring whether
linguistic features exhibit similar distinguishing signals when demographics
representing people with similar English language proficiency levels (e.g.,
native and highly proficient non-native speakers) are considered. In our
future work, we will focus on collecting, annotating, and analyzing review
data from diverse demographics, including variations in English fluency,
geographical regions, native language families, and socio-cultures.

## References

  * [1] Z. Wood-Doughty, M. Smith, D. Broniatowski, and M. Dredze, “How does twitter user behavior vary across demographic groups?,” in Proceedings of the Second Workshop on NLP and Computational Social Science, pp. 83–89, 2017.
  * [2] R. G. Guimaraes, R. L. Rosa, D. De Gaetano, D. Z. Rodriguez, and G. Bressan, “Age groups classification in social network using deep learning,” IEEE Access, vol. 5, pp. 10805–10816, 2017.
  * [3] X. Chen, Y. Wang, E. Agichtein, and F. Wang, “A comparative study of demographic attribute inference in twitter,” in Proceedings of the international AAAI conference on web and social media, vol. 9, pp. 590–593, 2015.
  * [4] F. Hsieh, R. Dias, and I. Paraboni, “Author profiling from facebook corpora,” in Proceedings of the Eleventh International Conference on Language Resources and Evaluation (LREC 2018), 2018.
  * [5] R. Dias and I. Paraboni, “Cross-domain author gender classification in brazilian portuguese,” in Proceedings of the 12th Language Resources and Evaluation Conference, pp. 1227–1234, 2020.
  * [6] D. Nguyen, D. Trieschnigg, A. S. Doğruöz, R. Gravel, M. Theune, T. Meder, and F. De Jong, “Why gender and age prediction from tweets is hard: Lessons from a crowdsourcing experiment,” in 25th International Conference on Computational Linguistics (COLING 2014), pp. 1950–1961, Dublin City University and Association for Computational Linguistics, 2014.
  * [7] S. Rosenthal and K. McKeown, “Age prediction in blogs: A study of style, content, and online behavior in pre-and post-social media generations,” in Proceedings of the 49th Annual Meeting of the Association for Computational Linguistics: Human Language Technologies, pp. 763–772, 2011.
  * [8] J. Schler, M. Koppel, S. Argamon, and J. W. Pennebaker, “Effects of age and gender on blogging.,” in AAAI spring symposium: Computational approaches to analyzing weblogs, vol. 6, pp. 199–205, 2006.
  * [9] S.-M. J. Wong and M. Dras, “Contrastive analysis and native language identification,” in Proceedings of the Australasian Language Technology Association Workshop 2009, pp. 53–61, 2009.
  * [10] M. Koppel, J. Schler, and K. Zigdon, “Automatically determining an anonymous author’s native language,” in International Conference on Intelligence and Security Informatics, pp. 209–217, Springer, 2005.
  * [11] D. Blanchard, J. Tetreault, D. Higgins, A. Cahill, and M. Chodorow, “Toefl11: A corpus of non-native english,” ETS Research Report Series, vol. 2013, no. 2, pp. i–15, 2013.
  * [12] S. Granger, “The international corpus of learner english: a new resource for foreign language learning and teaching and second language acquisition research,” Tesol Quarterly, vol. 37, no. 3, pp. 538–546, 2003.
  * [13] R. Sarkar, S. Mahinder, and A. KhudaBukhsh, “The non-native speaker aspect: Indian english in social media,” in Proceedings of the Sixth Workshop on Noisy User-generated Text (W-NUT 2020), pp. 61–70, 2020.
  * [14] S. Sazzed, “Impact of demography on linguistic aspects and readability of reviews and performances of sentiment classifiers,” International Journal of Information Management Data Insights, vol. 2, no. 2, p. 100135, 2022.
  * [15] S. Sazzed, “A hybrid approach of opinion mining and comparative linguistic analysis of restaurant reviews,” in Proceedings of the International Conference on Recent Advances in Natural Language Processing (RANLP 2021), pp. 1281–1288, 2021.
  * [16] S. Sazzed, “Stylometric and semantic analysis of demographically diverse non-native english review data,” in 2022 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining (ASONAM), pp. 470–476, IEEE, 2022.
  * [17] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine learning in Python,” Journal of Machine Learning Research, vol. 12, pp. 2825–2830, 2011.
  * [18] C. van der Lee and A. van den Bosch, “Exploring lexical and syntactic features for language variety identification,” in Proceedings of the fourth workshop on NLP for similar languages, varieties and dialects (VarDial), pp. 190–199, 2017.
  * [19] S. Volkova and E. Bell, “Identifying effective signals to predict deleted and suspended accounts on twitter across languages,” in Proceedings of the International AAAI Conference on Web and Social Media, vol. 11, pp. 290–298, 2017.
  * [20] G. Goldin, E. Rabinovich, and S. Wintner, “Native language identification with user generated content,” in Proceedings of the 2018 conference on empirical methods in natural language processing, pp. 3591–3601, 2018.
  * [21] M. Honnibal and I. Montani, “spacy 2: Natural language understanding with bloom embeddings, convolutional neural networks and incremental parsing,” To appear, 2017.
  * [22] C. Hutto and E. Gilbert, “Vader: A parsimonious rule-based model for sentiment analysis of social media text,” in Proceedings of the international AAAI conference on web and social media, vol. 8, pp. 216–225, 2014.
  * [23] M. Hu and B. Liu, “Mining and summarizing customer reviews,” in Proceedings of the Tenth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’04, (New York, NY, USA), pp. 168–177, ACM, 2004.
  * [24] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep bidirectional transformers for language understanding,” arXiv preprint arXiv:1810.04805, 2018.
  * [25] Y. Liu, M. Ott, N. Goyal, J. Du, M. Joshi, D. Chen, O. Levy, M. Lewis, L. Zettlemoyer, and V. Stoyanov, “Roberta: A robustly optimized bert pretraining approach,” arXiv preprint arXiv:1907.11692, 2019.
  * [26] T. Wolf, L. Debut, V. Sanh, J. Chaumond, C. Delangue, A. Moi, P. Cistac, T. Rault, R. Louf, M. Funtowicz, J. Davison, S. Shleifer, P. von Platen, C. Ma, Y. Jernite, J. Plu, C. Xu, T. L. Scao, S. Gugger, M. Drame, Q. Lhoest, and A. M. Rush, “Transformers: State-of-the-art natural language processing,” in Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations, (Online), pp. 38–45, Association for Computational Linguistics, Oct. 2020.


# A Generalized Semi-Analytic Model for Magnetar-Driven Supernovae

Conor M. B. Omand,1 Nikhil Sarin,2,3
1The Oskar Klein Centre, Department of Astronomy, Stockholm University,
AlbaNova, SE-106 91 Stockholm, Sweden
2Nordita, Stockholm University and KTH Royal Institute of Technology Hannes
Alfvéns väg 12, SE-106 91 Stockholm, Sweden
3The Oskar Klein Centre, Department of Physics, Stockholm University,
AlbaNova, SE-106 91 Stockholm, Sweden
E-mail<EMAIL_ADDRESS>
###### Abstract

Several types of energetic supernovae, such as superluminous supernovae
(SLSNe) and broad-line Ic supernovae (Ic-BL SNe), could be powered by the
spin-down of a rapidly rotating magnetar. Currently, most models used to infer
the parameters for potential magnetar-driven supernovae make several
unsuitable assumptions that likely bias the estimated parameters. In this
work, we present a new model for magnetar-driven supernovae that relaxes
several of these assumptions and an inference workflow that enables accurate
estimation of parameters from lightcurves of magnetar-driven supernovae. In
particular, in this model, we include the dynamical evolution of the ejecta,
coupling it to the energy injected by the magnetar itself while also allowing
for non-dipole spin down. We show that the model can reproduce SLSN and Ic-BL
SN light curves consistent with the parameter space from computationally
expensive numerical models. We also show the results of parameter inference on
four well-known example supernovae, demonstrating the model’s effectiveness at
capturing the considerable diversity in magnetar-driven supernova lightcurves.
The model fits each light curve well and recovers parameters broadly
consistent with previous works. This model will allow us to explore the full
diversity of magnetar-driven supernovae under one theoretical framework, more
accurately characterize these supernovae from only photometric data, and make
more accurate predictions of future multiwavelength emission to test the
magnetar-driven scenario better.

###### keywords:

supernovae: general – stars: magnetars – supernovae: individual: SN 2015bn –
supernovae: individual: SN 2007ru – supernovae: individual: ZTF20acigmel –
supernovae: individual: iPTF14gqr

††pubyear: 2023††pagerange: A Generalized Semi-Analytic Model for Magnetar-
Driven Supernovae–A

##  1 Introduction

Recent wide-field high-cadence surveys, such as the Zwicky Transient Facility
(ZTF, Bellm et al., 2019) and the Asteroid Terrestrial-impact Last Alert
System (ATLAS, Tonry, 2011; Tonry et al., 2018), can discover more than 1000
supernovae per year (Cappellaro, 2022), and next generation facilities such as
Rubin Observatory (Ivezić et al., 2019) will be able to discover more than
1000 supernovae per night (LSST Science Collaboration et al., 2009;
Stritzinger & Moriya, 2018). Core-collapse supernovae (CCSNe), caused by the
deaths of massive stars, tend to have explosion energies of $\sim$ 1051 erg
and radiated energies $\sim$ 1049 erg, which can be well explained by models
(Arnett, 1980, 1982). However, several classes of CCSNe have energies higher
than these standard models. Superluminous supernovae (SLSNe) radiate $\sim$
100 times more energy than a standard CCSN (Gal-Yam, 2012; Nicholl, 2021) and
broad-line Type Ic supernovae (SNe Ic-BL) have inferred kinetic energies
$\sim$ 10 times higher than a typical CCSN (e.g. Taddia et al., 2019). These
energetic supernovae have both been associated with long- or ultra-long gamma-
ray bursts (GRBs) (Gendre et al., 2013; Nakauchi et al., 2013; Levan et al.,
2014; Cano et al., 2017), and SLSNe and SNe Ic-BL also have similar host
galaxies (Lunnan et al., 2014; Leloudas et al., 2015a; Angus et al., 2016;
Schulze et al., 2018; Ørum et al., 2020), and have similar spectral features
at both early (Pastorello et al., 2010; Inserra et al., 2013; Nicholl et al.,
2013; Blanchard et al., 2019) and late (Milisavljevic et al., 2013; Jerkstrand
et al., 2017; Nicholl et al., 2017b) times. Other types of unusually energetic
supernovae, such as fast blue optical transients (FBOTs) and bright ultra-
stripped supernovae (USSNe), have been suggested to have similar power sources
as SLSNe and SNe Ic-BL (Liu et al., 2022; Sawada et al., 2022).

Several different models can be used to explain the high energies of one or
both of SLSNe or SNe Ic-BL. Stars with masses M∗ $\gtrsim$ 130 M⊙ can explode
as pair instability supernovae (PISNe) (Barkat et al., 1967; Heger & Woosley,
2002; Gal-Yam et al., 2009), which can generate tens of solar masses of 56Ni
and present as an extremely long-lived, luminous supernovae. Slightly less
massive stars, with 130 M⊙ $\gtrsim$ M∗ $\gtrsim$ 100 M⊙, can eject shells of
material through the pair-instability without being completely destabilized,
and collisions between the ejected shells can be as luminous as an SLSN (Heger
et al., 2003; Woosley et al., 2007; Chatzopoulos & Wheeler, 2012; Yoshida et
al., 2016; Woosley, 2017) \- these are known as pulsational pair-instability
supernovae (PPISNe). The collision of supernova ejecta and circumstellar
material surrounding the progenitor, ejected through either a steady wind,
eruptive mass loss, or binary interaction (Smith, 2014), can convert much of
the supernova kinetic energy into radiation, leading to a highly luminous
supernova (Chatzopoulos & Wheeler, 2012; Chatzopoulos et al., 2013; Villar et
al., 2017; Jiang et al., 2020). Finally, the compact remnant can inject some
energy into the ejecta; this energy may come from fallback accretion onto a
central black hole or neutron star (Dexter & Kasen, 2013; Moriya et al.,
2018), a collapsar or jet (MacFadyen & Woosley, 1999), or the spin-down energy
of a rapidly-rotating magnetar (Kasen & Bildsten, 2010; Woosley, 2010).

In the magnetar model, the spin-down energy is emitted as a highly magnetized
particle wind, which expands relativistically until it collides with the inner
edge of the supernova ejecta, create forward and reverse shocks. The expanding
wind becomes shocked by the reverse shock, accelerating the particles up to
ultrarelativistic energies, which then emit via synchrotron radiation and
inverse Compton scattering (Gaensler & Slane, 2006); this shocked wind is
known as a pulsar wind nebula (PWN). The PWN applies a pressure to the
supernova ejecta, causing it to accelerate, and the radiation from the nebula
is thermalized in the ejecta, causing the ejecta temperature and supernova
luminosity to both increase (Kasen & Bildsten, 2010). The PWN-ejecta
interaction can also cause Rayleigh-Taylor instabilities, which can shred the
inner ejecta and cause non-spherical structure to emerge in the ejecta (Suzuki
& Maeda, 2017, 2021).

The magnetar model predicts many multiwavelength signals that could be used to
identify and characterize the newborn neutron star. Once the ejecta becomes
optically thin, the non-thermal emission from the PWN can be detected
directly, either at high energy in hard X-rays or gamma rays (Kotera et al.,
2013; Murase et al., 2015; Kashiyama et al., 2016), or at low energy in radio
(Omand et al., 2018; Eftekhari et al., 2021; Murase et al., 2021); two
energetic SNe have radio detections at late times that are consistent with PWN
emission, PTF10hgi (Eftekhari et al., 2019; Law et al., 2019; Mondal et al.,
2020; Hatsukade et al., 2021) and SN2012au (Stroh et al., 2021). Dust formed
in the supernova can absorb PWN emission and re-emit that energy in infrared
(Omand et al., 2019), causing bright continuum emission recently seen in four
SLSNe (Chen et al., 2021; Sun et al., 2022). High-energy photons and Rayleigh-
Taylor induced mixing can change the chemical and ionization structure of the
ejecta, leading to unique signatures in the supernova nebular spectrum
(Chevalier & Fransson, 1992; Omand & Jerkstrand, 2023). Aspherical ejecta
caused by either hydro instabilities or an aspherical PWN can produce an
optical polarization signal (Tanaka et al., 2017), which has been detected in
some energetic SNe (Inserra et al., 2016; Saito et al., 2020; Pursiainen et
al., 2022; Poidevin et al., 2023; Pursiainen et al., 2023), but not others
(Leloudas et al., 2015b; Lee, 2019, 2020; Poidevin et al., 2022, 2023;
Pursiainen et al., 2023).

Accurate parameter estimation from the light curve around optical peak is
essential for a number of reasons. Firstly, new surveys such as the LSST will
detect SLSNe out to high redshift, where they can not all be classified with
spectroscopy. Therefore, these supernovae will have to be characterized from
their light curve data alone. Also, Predicting late-time multiwavelength
signals requires accurate parameter estimation, as different parameters that
produce similar optical light curves can produce vastly different
multiwavelength signals (e.g. Omand et al., 2018, 2019). The model used in
codes currently widely used for inference of magnetar-driven SNe (e.g Nicholl
et al., 2017b) make several assumptions to reduce computational complexity,
which are unjustified outside a small region of the parameter space. In
particular, they assume a constant ejecta velocity, which is independent of
the magnetar luminosity, although numerical simulations show that ejecta
acceleration due to PWN pressure plays a vital role in the dynamics of the
ejecta (Chen et al., 2016; Suzuki & Maeda, 2017, 2019; Chen et al., 2020;
Suzuki & Maeda, 2021). They also assume the magnetar spin down through pure
vacuum dipole emission, even though studies of Galactic pulsars (Lyne et al.,
2015; Parthasarathy et al., 2020) and putative magnetars born in GRBs (Lasky
et al., 2017; Sarin et al., 2020a, b) show that most neutron stars are
inconsistent with a pure vacuum dipole.

In this work, we present a model where these assumptions are relaxed, which
can fully explore the diversity of magnetar-driven supernovae and unite
phenomenologically different supernovae, such as SNe Ic-BL and SLSNe, under
one theoretical framework. In Section 2, we introduce our model for magnetar-
driven supernovae. In Section 3, we show the diversity of supernovae and
supernova observables resulting from differences in initial parameters. In
Section 4, we perform Bayesian inference on a few varying supernovae to show
how our model can consistently reproduce and explain them. Finally, in Section
5, we discuss the model’s implications and conclude. Throughout the paper, we
use the notation $Q_{x}=Q/10^{x}$ in cgs units unless otherwise noted.

##  2 Model

The physics of our model is based on previous magnetar-driven kilonovae models
(e.g. Yu et al., 2013; Metzger, 2019; Sarin et al., 2022), but modified to
describe supernovae. We present here a non-relativistic model description,
although the model implementation is fully relativistic. For the fully
relativistic description of the kinematics (Equations 5, 6, 14, 16, 17, 18,
and 19), see Sarin et al. (2022). We note that relativistic corrections are
largely unimportant for supernovae apart from transients with exceptionally
low ejecta masses and powerful magnetar engines.

###  2.1 Model Physics

The central magnetar spins down by releasing its rotational energy

$E_{\rm rot}=\frac{1}{2}I\Omega^{2},$ (1)

where $I$ is the moment of inertia of the magnetar and $\Omega$ is the
rotational angular frequency of the magnetar. The time derivative of this
relation gives the spin down luminosity,

$L_{\rm SD}=I\Omega\dot{\Omega},$ (2)

which, given $\dot{\Omega}\propto-\Omega^{n}$ for braking index $n$, can be
modelled generally as (Lasky et al., 2017)

$L_{\rm SD}(t)=L_{0}\left(1+\frac{t}{t_{\rm SD}}\right)^{\frac{1+n}{1-n}},$
(3)

where $L_{0}$ is the initial magnetar spin-down luminosity and $t_{\rm SD}$ is
the magnetar spin-down time. A braking index of $n=3$ corresponds to pure
vacuum dipole spin down (Ostriker & Gunn, 1969; Goldreich & Julian, 1969),
which gives a late-time spin-down down luminosity of $L\propto t^{-2}$ (Zhang
& Mészáros, 2001), while a braking index of $n=5$ corresponds to gravitational
wave spin down via magnetic deformation (Cutler & Jones, 2000), and has a
late-time spin-down luminosity of $L\propto t^{-3/2}$. We note that in
general, $n$ is expected to be variable during the early life of the magnetar
(Lander & Jones, 2018, 2020), but we keep it constant for simplicity.
Integrating Equation 3 gives the total rotational energy as a function of
braking index:

$E_{\rm rot}=\frac{n-1}{2}L_{0}t_{\rm SD}.$ (4)

This recovers $E_{\rm rot}=L_{0}t_{\rm SD}$ for vacuum dipole spin-down.

The rotational energy from the magnetar is converted into a pulsar wind. This
highly magnetized, ultrarelativistic wind collides with and pushes a shock
into the supernova ejecta, increasing its kinetic and internal energy. The
internal energy is also increased by the absorption of PWN photons by the
ejecta. The total energy of the system is the combination of the kinetic and
internal energy

$E_{\rm ej}=\frac{1}{2}M_{\rm ej}v_{\rm ej}^{2}+E_{\rm int},$ (5)

where $M_{\rm ej}$ is the ejecta mass. The evolution of this system is
governed by the energy sources, radioactive heating and magnetar spin-down
luminosity, and energy losses from radiated luminosity and adiabatic cooling
from expansion. The evolution of the internal energy of the ejecta is written
as (Kasen et al., 2016)

$\frac{dE_{\rm int}}{dt}=\xi L_{\rm SD}+L_{\rm ra}-L_{\rm
bol}-\mathcal{P}\frac{dV}{dt},$ (6)

where $L_{\rm ra}$ and $L_{\rm bol}$ are the radioactive power and emitted
bolometric luminosity, respectively, $\xi$ is the fraction of spin-down
luminosity injected into the ejecta, and $\mathcal{P}$ and $V$ are the
pressure and volume of the ejecta.

Here, we adopt the Wang et al. (2015) prescription for gamma-ray leakage used
in other models (e.g. Nicholl et al., 2017b; Sarin et al., 2022) with

$\xi=1-e^{-At^{-2}},$ (7)

where

$A=\frac{3\kappa_{\gamma}M_{\rm ej}}{4\pi v^{2}_{\rm ej}}$ (8)

is the leakage parameter and $\kappa_{\gamma}$ is the gamma-ray opacity of the
ejecta.

The radioactive power from the decay of 56Ni is given by

$L_{\rm ra}=f_{\rm Ni}M_{\rm
ej}(L_{\rm{}^{56}Ni}e^{-t/t_{\rm{}^{56}Ni}}+L_{\rm{}^{56}Co}e^{-t/t_{\rm{}^{56}Co}}),$
(9)

where $f_{\rm Ni}$ is the nickel fraction of the ejecta,
$L_{\rm{}^{56}Ni}=6.45\times 10^{43}$ erg s-1 $M_{\odot}^{-1}$ and
$L_{\rm{}^{56}Co}=1.45\times 10^{43}$ erg s-1 $M_{\odot}^{-1}$ are the decay
luminosities of 56Ni and 56Co, and $t_{\rm{}^{56}Ni}=8.8$ days and
$t_{\rm{}^{56}Co}=111.3$ days are the decay timescales for 56Ni and 56Co
(Nadyozhin, 1994).

The dynamical evolution of the ejecta is given by

$\frac{dv_{\rm ej}}{dt}=\frac{c^{2}\mathcal{P}(d\mathcal{V}/dt)}{M_{\rm
ej}v_{\rm ej}^{3}},$ (10)

where

$\displaystyle\mathcal{V}$ $\displaystyle=\frac{4}{3}\pi R_{\rm ej}^{3},$ (11)
$\displaystyle\frac{d\mathcal{V}}{dt}$ $\displaystyle=4\pi R_{\rm
ej}^{2}v_{\rm ej},$ (12) $\displaystyle\mathcal{P}$
$\displaystyle=\frac{E_{\rm int}}{3\mathcal{V}}.$ (13)

Substituting these into Equation 10 gives

$\frac{dv_{\rm ej}}{dt}=\frac{c^{2}E_{\rm int}}{M_{\rm ej}R_{\rm ej}v_{\rm
ej}^{2}}.$ (14)

The initial ejecta velocity is set by

$v_{\rm ej,0}=\sqrt{\frac{2E_{\rm SN}}{M_{\rm ej}}},$ (15)

where $E_{\rm SN}$ is the supernova explosion energy.

The bolometric radiated luminosity is (Kasen & Bildsten, 2010; Kotera et al.,
2013)

$\displaystyle L_{\rm bol}=$ $\displaystyle\frac{E_{\rm int}c}{\tau R_{\rm
ej}}=\frac{E_{\rm int}t}{t_{\rm dif}^{2}}$ $\displaystyle(t\leq t_{\tau}),$
(16) $\displaystyle=$ $\displaystyle\frac{E_{\rm int}c}{R_{\rm ej}},$
$\displaystyle(t>t_{\tau}),$ (17)

where

$\tau=\frac{\kappa M_{\rm ej}R_{\rm ej}}{\mathcal{V}}$ (18)

is the optical depth of the ejecta, $\kappa$ is the ejecta opacity,

$t_{\rm dif}=\left(\frac{\tau R_{\rm ej}t}{c}\right)^{1/2}$ (19)

is the effective diffusion time, and $t_{\tau}>t_{\rm dif}$ is the time when
$\tau=1$.

Calculating the bolometric luminosity of the magnetar-driven transient
(Equations 16 and 17) involves solving the evolution of the internal energy
and dynamics (Equations 6 and 14 respectively) using the input power sources
(Equations 3 and 9). The photospheric temperature is determined from the
bolometric luminosity and ejecta radius until the temperature reaches the
photospheric plateau temperature, as in Nicholl et al. (2017b). This can be
expressed as

$T_{\rm phot}(t)=\begin{cases}\left(\frac{L_{\rm bol}(t)}{4\pi\sigma R_{\rm
ej}^{2}}\right)^{1/4}&\text{ for }\left(\frac{L_{\rm bol}(t)}{4\pi\sigma
R_{\rm ej}^{2}}\right)^{1/4}>T_{\rm min},\\\ T_{\rm min}&\text{ for
}\left(\frac{L_{\rm bol}(t)}{4\pi\sigma R_{\rm ej}^{2}}\right)^{1/4}\leq
T_{\rm min}\end{cases}$ (20)

The spectral energy distribution (SED) of the transient is then calculated
using the cutoff blackbody used in Nicholl et al. (2017b), with
$F_{\lambda<\lambda_{\rm cut}}=F_{\lambda}(\lambda/\lambda_{\rm cut})$ and
$\lambda_{\rm cut}=3000\AA$ (Chomiuk et al., 2011; Nicholl et al., 2017a), but
this can also be switched to a simple blackbody.

###  2.2 Parameters, Priors, and Implementation

The model presented above is implemented into the open-source electromagnetic
transient fitting software package, Redback (Sarin et al., 2023). The input
parameters for the model, their default priors, and the values used in Section
3 are listed in Table 1.

Parameter | Definition | Units | Default Prior | Section 3.1 Value
---|---|---|---|---
$L_{0}$ | Initial Magnetar Spin-Down Luminosity | erg s-1 | L[$10^{40}$, $10^{50}$] | Varying
$t_{\rm SD}$ | Spin-Down Time | s | L[$10^{2}$, $10^{8}$] | Varying
$n$ | Magnetar Braking Index |  | U[1.5, 10] | 3
$f_{\rm Ni}$ | Ejecta Nickel Mass Fraction |  | L[10-3, 1] | 0
$M_{\rm ej}$ | Ejecta Mass | $M_{\odot}$ | L[10-2, 1 $\times$ 102] | Varying
$E_{\rm SN}$ | Supernova Explosion Energy | erg | L[$5\times 10^{50}$, $2\times 10^{51}$] | $10^{51}$
$\kappa$ | Ejecta Opacity | cm2 g-1 | L[0.05, 0.2] | 0.1
$\kappa_{\gamma}$ | Ejecta Gamma-Ray Opacity | cm2 g-1 | L[$10^{-4}$, $10^{4}$] | 0.1
$T_{\rm min}$ | Photospheric Plateau Temperature | K | L[$10^{3}$, $3\times 10^{4}$] | 5000
Table 1: The parameters and default priors for the generalized magnetar model,
as well as the values used for the parameter exploration in Section 3. Priors
are either uniform (U) or log-uniform (L).

Due to the acceleration of the ejecta from the PWN, $v_{\rm ej}$ is not
constant and can not be used as a free parameter, since it is coupled to both
$E_{\rm SN}$ and $L_{\rm SD}$. Velocity information from spectroscopy can be
used to weight the results, although the velocity measured from absorption
widths is not the same as the photospheric velocity which is not the same as
the ejecta velocity (Arnett, 1982), so we caution against this unless the
velocities are well calibrated (e.g. Dessart et al., 2016).

Our model differs from others (e.g. Nicholl et al., 2017b) in the choice of
input parameters used to determine the magnetar luminosity. Previous models
use the initial magnetar spin period $P_{0}$, dipole component of the pulsar
magnetic field $B$, and neutron star mass $M_{\rm NS}$; while we use $L_{0}$
and $t_{\rm SD}$ (and $n$, which is implicitly fixed to 3 in other models). We
note that as the magnetar luminosity for vacuum dipole spin-down can
determined from only $L_{0}$ and $t_{\rm SD}$ (see Equation 3), using three
parameters is unnecessary for parameter inference. Using these parameters also
avoids assumptions such as the moment of inertia of the neutron star, which
depends on the equation of state (EoS) (Lattimer & Schutz, 2005) and can vary
depending on the mass and spin period of the neutron star (Worley et al.,
2008). To recover parameters such as the magnetar spin period and magnetic
field, one can use the scalings used in other models such as

$\displaystyle L_{0}=$ $\displaystyle 2.0\times
10^{47}P_{0,-3}^{-4}B_{14}^{2},$ (21) $\displaystyle t_{\rm SD}=$
$\displaystyle 1.3\times 10^{5}P_{0,-3}^{2}B_{14}^{-2}\left(\frac{M_{\rm
NS}}{1.4M_{\odot}}\right)$ (22)

for transients consistent with $n=3$, which assumes a neutron star with moment
of inertia of $\sim$ 1.3 $\times$ 1045 g cm2 for a 1.4 $M_{\odot}$ neutron
star, which is consistent with the APR EoS (Akmal et al., 1998) or MDI EoS
with density independent nuclear symmetry energy (Das et al., 2003; Shetty et
al., 2007), which both give neutron star radii $\sim$ $11.5-12$ km (Worley et
al., 2008). These scalings become more complicated for $n\neq 3$ (e.g.,
Shapiro & Teukolsky, 1983), with other dependencies such as the ellipticity of
the neutron star or bulk viscosity (depending on the spin-down processes
involved), leaving it difficult to definitively recover a spin period or
magnetic field with such simplified models.

##  3 Results

We now explore the diversity of magnetar driven-supernovae for different
initial conditions using the model derived in Section 2. We assume the
supernova explosion energy is 1051 erg, typical for a neutrino-driven
explosion, and the supernova is entirely powered by magnetar energy after the
initial explosion, with no contribution from 56Ni. We fix both the ejecta
opacity and gamma-ray opacity to 0.1 cm2 g-1; the former is typical of
stripped-envelope supernovae (Inserra et al., 2013; Kleiser & Kasen, 2014) and
the latter will not affect the light curve peak properties unless it is
extremely low, so we use a wavelength-independent opacity here for simplicity.
We also fix the photospheric plateau temperature to 5000 K, which is typical
for SLSNe (Nicholl et al., 2017b).

First, we show the diversity of supernovae that can be produced using this
model in Section 3.1. To compare with previously derived results, we assume
vacuum dipole spin-down and use Equations 21 and 22 with a 1.4 $M_{\odot}$
neutron star to express our initial conditions in terms of $P_{0}$ and $B$;
the mapping between ($L_{0}$, $t_{\rm SD}$) and ($P_{0}$, $B$) is shown in
Figure 1. Then we show the effect of changing braking index in Section 3.2.

Figure 1: Initial pulsar luminosities $L_{0}$ (top) and spin-down times
$t_{\rm SD}$ (bottom) for different initial pulsar rotation periods and
magnetic fields.

###  3.1 Energetics, Timescales, and Observables

We first see if our model can reproduce typical observables for the observed
populations of SLSNe and SNe Ic-BL. SLSNe at optical light curve peak
typically have g-band absolute magnitudes $-23<M_{g}<-20$, with most brighter
than $-$21; spectroscopically determined photospheric velocities of $\sim$ 10
000 $-$ 15 000 km s-1; rise times of 10 $-$ 70 days; and $g-r$ of $-$0.3 $-$
0.3 (Chen et al., 2023a). SNe Ic-BL at optical light curve peak typically have
r-band absolute magnitudes around $-20<M_{r}<-18$; spectroscopically
determined photospheric velocities of $\sim$ 15 000 $-$ 30 000 km s-1; rise
times of 5 $-$ 20 days; and $g-r$ of 0 $-$ 0.5 (Taddia et al., 2019). The
inferred ejecta masses of typical SLSNe and SNe Ic-BL are both around 5
$M_{\odot}$ (Taddia et al., 2019; Chen et al., 2023b), and SLSNe show a
negative correlation between initial spin period and ejecta mass (Blanchard et
al., 2020).

Several inefficiencies can prevent all the magnetar spin-down luminosity from
eventually being emitted as supernova luminosity. Some fraction ($\xi$ from
Equation 6) will escape without interacting with the ejecta at all, and
possibly be detected as non-thermal x-rays or gamma rays. Some fraction of the
energy will also accelerate the ejecta instead of thermalizing and being re-
emitted, which will affect both the supernova luminosity and peak timescale;
this fraction can be determined by the ratio of the spin-down time and the
diffusion time (Suzuki & Maeda, 2021; Sarin et al., 2022).

Figure 2 shows the ratio of the final supernova kinetic and radiated energies
for various magnetic fields and ejecta masses for spin periods of 1 ms (close
to the mass shedding limit for neutron stars (Watts et al., 2016)) and 3 ms
(where the spin-down luminosity and explosion energy can become comparable),
and for various magnetic fields and spin periods for an ejecta mass of 10
$M_{\odot}$. The energy ratio $E_{\rm kin}$/$E_{\rm rad}$ does correlate
strongly with $\zeta=t_{\rm SD}/t_{\rm dif}$ up to large ejecta mass for low
spin periods, although at spin period increases, the correlation gets weaker
as the behaviour of the energy ratio changes; this is due to the total amount
of energy injected by the magnetar prior to the diffusion time decreasing
below the explosion energy, meaning that the ejecta dynamics are no longer
primarily determined by the magnetar spin-down energy. There is a small region
of the parameter space at low spin period, magnetic field, and ejecta mass
where the radiated energy can surpass the kinetic energy. However, for most of
the parameter space this ratio is between 1 and 100. Typical supernovae have
$E_{\rm kin}$/$E_{\rm rad}$ = 1051 erg / 1049 erg = 100, but without a
contribution from 56Ni this ratio can go much higher.

$P_{\rm 0}$ = 1 ms $P_{\rm 0}$ = 3 ms $M_{\rm ej}$ = 10 $M_{\odot}$

Figure 2: Ratio of kinetic to radiated energy for supernovae with varying
ejecta mass and $P_{0}=1$ ms (left) and $P_{0}=3$ ms (middle) and with varying
spin period and $M_{\rm ej}=10M_{\odot}$ (right). The black lines indicate
contours of constant $E_{\rm kin}/E_{\rm rad}$, while the magenta lines
indicate contours of constant $\zeta=t_{\rm SD}/t_{\rm dif}$.

Figure 3 shows the peak timescale of the bolometric luminosity and the final
ejecta velocity over the same parameter grid as Figure 2. The trends match our
theoretical expectations, as $P_{0}$ increases, the total energy of the system
decreases, causing the ejecta velocity to decrease and the peak timescale to
increase. Conversely, as $B$ increases, the spin-down time decreases, causing
the ejecta velocity to increase and the peak timescale to decrease. Finally,
as $M_{\rm ej}$ increases, the diffusion time increases, causing the ejecta
velocity to decrease and the peak timescale to increase. The final ejecta
velocity is largely insensitive to the magnetic field strength above a certain
field threshold despite a change in peak timescale. This is a product of the
coupling of dynamical evolution of the ejecta to the magnetar’s rotational
energy, as the final velocity is reached earlier in the supernova evolution
for higher magnetic field (due to their shorter spin-down timescale). For a
magnetar with rotation period of 1 ms, a magnetic field of $>$ 1015 is needed
to reduce the peak timescale below 10 days, even at an ejecta mass of 1
$M_{\odot}$, meaning that the fastest SNe Ic-BL likely require both an ejecta
mass below 1 $M_{\odot}$ and a magnetar spinning at close to breakup speeds.
Timescales of $>$ 20 days require ejecta masses below 5 $M_{\odot}$, meaning
that our model will likely estimate a lower ejecta mass for SNe Ic-BL than
Taddia et al. (2019). Timescales and ejecta velocities typical of SLSNe can be
reproduced over a large portion of the parameter space. However, a higher
ejecta mass is required for faster spinning magnetars to keep the velocities
below that of SNe Ic-BL. In contrast, a low ejecta mass is required for
faster-spinning magnetars to keep the timescales below $\sim$ 100 days and the
velocities higher than $\sim$ 10 000 km s-1, providing some phenomenological
justification for the mass-spin correlation found by Blanchard et al. (2019).

$P_{\rm 0}$ = 1 ms $P_{\rm 0}$ = 3 ms $M_{\rm ej}$ = 10 $M_{\odot}$ $t_{\rm
peak}$ $v_{\rm ej}$

Figure 3: Bolometric peak timescale $t_{\rm peak}$ (top) and final ejecta
velocity $v_{\rm ej}$ (bottom) for supernovae with varying ejecta mass and
$P_{0}=1$ ms (left) and $P_{0}=3$ ms (middle) and with varying spin period and
$M_{\rm ej}=10M_{\odot}$ (right). The black lines indicate notable values of
$t_{\rm peak}$ and $v_{\rm ej}$.

Figure 4 shows the peak $g$-band absolute magnitude and peak $g-r$ colour over
the same parameter grid. We find a portion of parameter space for low spin
period, ejecta mass, and magnetic field which produces a transient more
luminous than any previously observed superluminous supernova. However, it is
unlikely such a combination (particularly low spin period and magnetic fields)
can be conceived as magnetic-field amplification mechanisms such as the
magneto-rotational instability or Kelvin-Helmholtz instability likely amplify
most typical progenitor fields to larger poloidal fields than seen in this
parameter space (e.g., Reboul-Salze et al., 2021), meanwhile, the stability of
magnetic-field configurations in this part of the parameter space is also
questionable (e.g., Braithwaite, 2009) and so magnetars in this parameter
space may never materialise. We note that as we do not track gravitational-
wave losses, the newly born magnetar could potentially spin-down rapidly
through gravitational-wave radiation (Sarin & Lasky, 2021) in this parameter
space, depleting the energy reservoir to power such a luminous transient. To
compound this all further, it is unknown whether stellar explosions with such
small ejecta masses could harbour magnetars that are rapidly rotating but have
quite weak poloidal fields in the first place. The parameter space between
$M_{g}=-21$ and $M_{g}=-23$, where SLSN are, shifts to lower masses for higher
spin periods, where the parameter space around $M_{g}=-19$, where most SNe Ic-
BL are, require either a large ejecta mass ($\gtrsim 5M_{\odot}$), higher spin
period, or extremely high magnetic field. The parameter space where $g-r<0$
mostly overlaps with the $M_{g}<-21$ region, showing that most SLSNe should
have $-0.5<g-r<0$ at peak, which is broadly consistent with observations (Chen
et al., 2023a).

$P_{\rm 0}$ = 1 ms $P_{\rm 0}$ = 3 ms $M_{\rm ej}$ = 10 $M_{\odot}$ $M_{\rm
g}$ $g-r$

Figure 4: $g$-band absolute magnitude $M_{g}$ (top) and $g-r$ colour (bottom)
at peak for supernovae with varying ejecta mass and $P_{0}=1$ ms (left) and
$P_{0}=3$ ms (middle) and with varying spin period and $M_{\rm
ej}=10M_{\odot}$ (right). The black lines indicate notable values of $M_{g}$
and $g-r$.

###  3.2 Effect of Varying Braking Index

Light curve luminosity and morphology can vary significantly with variations
in magnetar braking index. Figure 5 shows the bolometric luminosity, absolute
$g$-band magnitude, and absolute $r$-band magnitude for several supernovae
where only the braking index $n$ is varied. The ejecta mass, spin-down time,
and total rotational energy are fixed to 10 $M_{\odot}$, 106 s, and 1052 erg,
with initial magnetar luminosity calculated from Equation 4; all other
parameters are the same as in Section 3.1. The timing of the light curve peak
can vary by a factor of $\sim$ 3 and late-time luminosities by orders of
magnitude, with large $n$ peaking later and having higher luminosities at
later times, although the variation in late-time luminosity asymptotes as $n$
increases due to the exponent in Equation 3 asymptoting to $-1$. The peak
luminosity can vary by $\sim$ 1 mag and is highest for $n\approx$ 3 in this
case, although this will likely vary depending on the energetics and diffusion
time of the supernova.

Figure 5: Bolometric luminosity (left), absolute $g$-band magnitude (middle),
and absolute $r$-band magnitude (right) for several supernovae where only the
braking index $n$ is varied.

##  4 Case Studies for Inference

As a proof of concept, we perform inference on several different classes of
supernovae to see if the model is flexible enough to recover sensible
parameters for a variety of objects. First, we validate the model using a
simulated SLSN. Then, we perform inference on SN 2015bn, an SLSN; SN 2007ru, a
SN Ic-BL; ZTF20acigmel (better known as “the Camel”), an FBOT; and iPTF14gqr,
a USSN.

Inference is performed on multiband photometry using the open-source software
package Redback (Sarin et al., 2023) with the dynesty sampler (Speagle, 2020)
implemented in Bilby (Ashton et al., 2019; Romero-Shaw et al., 2020). We
sample with a Gaussian likelihood and an additional white noise term, and
sample in flux density rather than magnitude. We use the default priors in all
cases (shown in Table 1) except for the explosion energy of iPTF14gqr, where
the lower limit is reduced to 5 $\times$ 1048 erg to capture the lower
expected explosion energies of USSNe (Suwa et al., 2015). We also sample the
unknown explosion time with a uniform prior of up to 100 days before the first
observation and an extinction term $A_{v}$, and use a constraint that the
total rotational energy of the magnetar $E_{\rm rot}\lesssim 10^{53}$ erg.

The fitted light curve and corner plot for the simulated SN are shown in
Figure 6, the fitted light curves for the other SNe are shown in Figure 7, the
input parameters for the simulated SN and recovered parameters for each SN are
shown in Table 2, and the corner plots for each of the real SNe are shown in
Appendix A.

Figure 6: Fitted light curve (left) and posteriors of key parameters (right)
for the simulated SLSN. The solid lines in the light curve plot indicate the
light curve from the model with the highest likelihood, while the shaded area
indicates the 90$\%$ credible interval. The orange dots and lines in the
posterior indicate the injected parameters.

Figure 7: Fitted light curves for the four supernovae from our case studies. The solid lines indicate the light curve from the model with the highest likelihood, while the shaded area indicates the 90$\%$ credible interval. Object | Supernova type | log($L_{0}$) [erg s-1] | log($t_{\rm SD}$) [s] | $n$ |  $M_{\rm ej}$ [$M_{\odot}$]
---|---|---|---|---|---
Simulated (Injected) |  | 46.0 | 6.0 | 3.0 | 10.0
Simulated (Recovered) |  | 46.47${}^{+0.55}_{-0.58}$ | 5.56${}^{+0.59}_{-0.60}$ | 4.70${}^{+2.51}_{-1.65}$ | 16.94${}^{+6.53}_{-5.46}$
SN2015bn | SLSN | 45.79${}^{+0.53}_{-0.07}$ | 6.98${}^{+0.15}_{-0.41}$ | 3.21${}^{+1.06}_{-0.74}$ | 9.57${}^{+2.25}_{-1.76}$
SN2007ru | Ic-BL | 46.16${}^{+1.14}_{-1.22}$ | 3.32${}^{+0.90}_{-0.83}$ | 5.47${}^{+3.32}_{-0.83}$ | 1.49${}^{+0.35}_{-0.29}$
ZTF20acigmel | FBOT | 46.91${}^{+1.11}_{-0.78}$ | 4.08${}^{+0.70}_{-1.07}$ | 3.57${}^{+2.09}_{-0.84}$ | 0.23${}^{+0.09}_{-0.07}$
iPTF14gqr | USSN | 43.09${}^{+0.23}_{-0.16}$ | 6.23${}^{+0.28}_{-0.37}$ | 3.77${}^{+2.32}_{-1.53}$ | 0.18${}^{+0.09}_{-0.06}$
Table 2: Injected parameters for the simulated supernovae and median inferred
parameter values and 1$\sigma$ uncertainties for the simulated supernova and
the four supernovae from the case studies.

###  4.1 Validation on a Simulated Supernova

To test if the model could recover parameters correctly, we simulated an SLSN
with $L_{0}=10^{46}$ erg s-1, $t_{\rm SD}=10^{6}$ s, $n=3$, and $M_{\rm
ej}=10$ $M_{\odot}$, with the other parameters the same as in Section 3. The
supernova was placed at redshift $z=0.5$ and data was generated using the
Redback simulation workflow as observed by ZTF in $g$, $r$ and $i$ band for
the first 200 days post explosion.

The light curve (Figure 6 (left)) is fit well by the model throughout its
evolution, and three of the four injected parameters are recovered to within
$1\sigma$ in the one-dimensional posteriors, with only $M_{\rm ej}$ being
slightly outside that error region, and every parameter recovered to within
$1\sigma$ in the two-dimensional posteriors due to correlations in each
parameter. The correlation between $L_{0}$ and $t_{\rm SD}$ also suggests that
the rotational energy of the magnetar is well constrained to $\sim$ 1052 erg,
i.e., the injected value.

###  4.2 SN 2015bn

SN 2015bn is an SLSN-I at $z=0.1136$ that was first discovered by the Catalina
Sky Survey on 2014 December 23. It peaked at 79 rest-frame days post
discovery, which made it one of the slowest evolving SLSNe at the time, and
had peak magnitudes of $M_{g}=-22.0\pm 0.08$ mag (AB) and $M_{U}=-23.07\pm
0.09$ mag (Vega), making it one of the most luminous as well (Nicholl et al.,
2016a). Since then, it has been followed-up extensively in optical/UV/NIR,
with photometry and spectroscopy (Nicholl et al., 2016a, b, 2018), and
polarimetry (Inserra et al., 2016; Leloudas et al., 2017), and well as in
radio (Nicholl et al., 2018; Eftekhari et al., 2021; Murase et al., 2021) and
X-rays (Inserra et al., 2017; Bhirombhakdi et al., 2018). SN 2015bn shows
strong undulations in the light curve on a timescale of 30-50 days (Nicholl et
al., 2016a; Inserra et al., 2017), and was detectable in optical/UV for more
than 1000 days (Nicholl et al., 2018), although the supernova has yet to be
detected in either radio or x-rays.

We import the observational data (Nicholl et al., 2016a, b) from the Open
Supernova Catalog (Guillochon et al., 2017). Our model is able to fit the
supernova peak very well in all bands and can reproduce $r$ and $i$ band data
for more than 400 days, although although it finds much bluer emission than
observed in the post-peak photospheric phase. The inferred rotational energy
of the magnetar is $\sim 6\times 10^{52}$ erg, which is close to the maximum
rotational energy the can be extracted from a newborn magnetar. The magnetar
also shows $n\approx 3$, meaning the vacuum dipole is a good approximation for
this object. Using the scaling relations from Equations 21 and 22 for a 1.4
$M_{\odot}$ neutron star with the same equation of state as Nicholl et al.
(2017b), we get a spin period $P_{0}\approx 0.7$ ms, extremely close to the
mass-shedding limit, and a magnetic field $B\approx 8\times 10^{12}$ G.

Comparing to the results of Nicholl et al. (2017b) shows some parameters in
agreement, but strong discrepancy between others. Using Equations 4, 21, and
22 with parameters from Nicholl et al. (2017b) show that they find a spin-down
timescale of $\sim$ 105 days, in good agreement with our $\sim$ 110 days, but
a magnetar rotational energy of only $\sim$ 1052 erg, a factor of $\sim$ 5
lower than recovered by our model. This is because the magnetar needs to
supply both the radiated and kinetic energy of the supernova in our model,
while in the Nicholl et al. (2017b) model the ejecta velocity is input
separately. We also recover a similar ejecta mass of $\sim$ 10 $M_{\odot}$.

###  4.3 SN 2007ru

SN 2007ru is an SN Ic-BL at $z=0.01546$ that was first discovered by the
Himalayan Chandra Telescope on 2007 December 2 (Sahu et al., 2009). The peak
magnitude of the supernova was $M_{V}\approx-19.06$ mag, one of the brighter
SNe Ic-BL, and estimated rise time was $8\pm 3$ days (Sahu et al., 2009). The
photospheric velocity was estimated to be around 20 000 km s-1. A nickel-
powered model estimated the $M_{\rm ej}\approx$ 1.3 $M_{\odot}$ and $M_{\rm
Ni}\approx$ 0.4 $M_{\odot}$ (Sahu et al., 2009), while a previous magnetar-
powered model estimated $P_{0}\approx 2.30$ ms, $B\approx 6.2\times 10^{15}$
G, and $M_{\rm ej}\approx 4.43M_{\odot}$ (Wang et al., 2016).

We import the observational data (Sahu et al., 2009) from the Open Supernova
Catalog (Guillochon et al., 2017). The model fits most of the data well in the
optical, although it does slightly underpredict the NIR data around peak. The
magnetar energy is $\lesssim 10^{50}$ erg for this supernova, lower than the
explosion energy, which dominates the dynamics here. This energy is also much
lower than that estimated by Wang et al. (2016), and the spin-down time is a
factor of $\sim$ 10 larger. The ejecta mass we estimate is similar to the
previous nickel-powered model, but lower than the previous magnetar-powered
model. Finally, the braking index $n$ is much higher than 3, and we can reject
vacuum dipole spin-down at at $\gtrsim 95\%$ confidence. Our posterior on the
braking index is also consistent with the magnetar spin-down being dominated
by gravitational-wave radiation (Sarin et al., 2018, 2020a). However we
caution against strong conclusions based on the measurement of $n$ due to our
simplified treatment of the neutron star spin evolution.

###  4.4 ZTF20acigmel

AT2020xnd or ZTF20acigmel, the ‘Camel’, is an FBOT at z = 0.2433 that was
first discovered by ZTF (Bellm et al., 2019) on 2020 October 12 (Perley et
al., 2021). The peak magnitude was $M_{\rm 5000\AA}\approx-20.6$ mag or
$M_{\rm 3900\AA}\approx-20.9$ mag, and estimated rise time was $\sim 2$ days
(Perley et al., 2021). The photospheric radius is already receding at 7 days,
and the photospheric temperature from 7-13 days is estimated to be 20 000
$\pm$ 2000 K. ZTF20acigmel was also found to be luminous in both radio (Ho et
al., 2022) and x-rays (Bright et al., 2022).

We used the publicly available photometric data from Perley et al. (2021) for
our fit. The model fits the data in all filters throughout the evolution of
the object. The total rotational energy of the magnetar is around 1051 erg,
comparable to the explosion energy. The braking index is consistent with
vacuum dipole, and using the scaling relations and the same 1.4 $M_{\odot}$
neutron star as for SN 2015bn, we can derive an initial spin period of $\sim$
5 ms and magnetic field of $\sim 2\times 10^{15}$ G; these values, as well as
the ejecta mass, are consistent with values found for the FBOT distribution
(Liu et al., 2022).

###  4.5 iPTF14gqr

iPTF14gqr is a USSN at $z=0.063$ that was first discovered by the intermediate
Palomar Transient Factory (iPTF) (Law et al., 2009) on 2014 October 14. The
supernova featured a bright first peak that faded within a day, followed by a
more extended light curve the rises after about $\sim$ 4 days. The first peak
can be explained by shock cooling (De et al., 2018) and the second peak by
either radioactivity (De et al., 2018) or a new born magnetar (Sawada et al.,
2022), although both peaks can also be explained by interaction alone (Khatami
& Kasen, 2023).

We used the publicly available photometric data from De et al. (2018),
although we exclude all data points within one day post-explosion, since we
only claim that the second peak could be magnetar powered. The model fits the
data in all filters throughout the evolution of the object. The explosion
energy, which has a wider prior to account for the low explosion energy of
USSNe, is constrained to $\sim 6\times 10^{49}$ erg, while the total magnetar
rotational energy is $\sim 2\times 10^{49}$ erg. The braking index is not very
well constrained; although it is consistent with vacuum dipole to within
error, the most likely value is $\sim$ 2\. The median spin-down time of
$\approx$ 20 days matches that found by Sawada et al. (2022), although our
magnetar energy is a factor of $\sim$ 5 higher; this is likely because our
explosion energy is smaller. The ejecta mass we derive is also consistent with
estimates by both De et al. (2018) and Sawada et al. (2022).

##  5 Discussion and Summary

As shown by both the exploration of the parameter space (Section 3) and the
case studies (Section 4), the model presented here is incredibly versatile.
Three of the four case studies had previous been fit by a magnetar model
(Nicholl et al., 2017b; Wang et al., 2016; Sawada et al., 2022), but each of
these studies used a different model with different assumptions. A versatile
model allows comparisons of different populations to be done self-consistently
and determine what model variations manifest in vastly different types of
supernovae, as well as probe whether a continuum between these sources could
exist or whether there are multiple distinct classes.

Much of the flexibility of our model comes from self-consistent dynamical
evolution of the ejecta, while the addition of magnetar braking index as a
parameter allows for some possible insight into the spin-down mechanism of the
newborn millisecond magnetar. Figure 8 (top) shows the posterior probability
distribution of braking index for the four case study supernova. As mentioned
above, SN 2015bn is very well approximated by $n=3$, while $n=3$ is rejected
in SN 2007ru at $>95\%$ confidence; both of these supernovae have posteriors
that are roughly Gaussian. ZTF20acigmel and iPTF14gqr both show non-
Gaussianity in their posteriors; the Camel has a tail at high $n$ but peaks
very close to $n=3$, while iPTF14gqr is not well constrained but peaks closer
to $n=2$ and has significant posterior support at $n=1.5$, the lowest value of
our prior. While making definitive statements about spin-down mechanisms will
require a much larger sample, this small sample already shows diversity in
their inferred braking index, highlighting a potential interesting question
for future studies.

Figure 8: Posterior distributions of the magnetar braking index (top) and
rotational energy (bottom) for the four supernovae from our case studies. The
median values are indicated by the blue vertical lines within the
distribution. The black dashed line indicates vacuum dipole spin-down ($n=3$).

All the objects studied in Section 4, including the simulated supernova, show
a strong negative correlation between $L_{0}$ and $t_{\rm SD}$. This shows
that the magnetar rotational energy $E_{\rm rot}$, and thus the total energy
budget of the supernova, can be constrained for these objects to within an
order of magnitude (see Figure 8 (bottom)). The ejecta masses were all found
to be similar to previous studies and the spin-down times were found to be
very close to previous models for SN 2015bn and iPTF14gqr (Nicholl et al.,
2017b; Wang et al., 2016; Sawada et al., 2022), although the total magnetar
energy budget was different in each case due to the way we treated our
dynamics. The magnitude of these discrepancies seems to vary depending on the
supernova, and a large-scale sample study on supernovae that have been
previously characterized by a magnetar model (e.g. Chen et al., 2023b) is
necessary to characterize the systematic difference between our model and
previous models.

This model has a few caveats which may prevent it from properly describing
certain transients. The first is that it is a one-zone, one-dimensional model.
This makes the treatment of the photospheric radius very simplified compared
to real supernovae. Engine-driven supernovae also show hydrodynamic
instabilities in multidimensional simulations (e.g. Chen et al., 2016, 2020;
Suzuki & Maeda, 2017, 2019, 2021) which can shred the inner ejecta, causing a
decrease in the effective optical depth of the ejecta, as well as affecting
the timescale for non-thermal leakage. If the spin-down timescale of the
magnetar is smaller than the Kelvin-Helmholtz timescale of the magnetar for
neutrino emission $t_{\rm KH,\nu}\lesssim 100$ s, then baryon loading on the
magnetized wind via the neutrino-driven wind can be relevant (Thompson et al.,
2004) and the magnetized wind can be collimated by anisotropic and hoop stress
(Bucciantini et al., 2007, 2008) and form a jet (Kashiyama et al., 2016). This
model also does not self-consistently track gravitational-wave emission and
how it depletes the overall rotational energy reservoir. The model also has no
way to explain the bumps and undulations that have recently found in a large
number of SLSNe (Hosseinzadeh et al., 2022; Chen et al., 2023b), which have
been explained by both circumstellar material (e.g. West et al., 2023; Chugai
& Utrobin, 2023) and magnetars (Chugai & Utrobin, 2022; Moriya et al., 2022;
Dong et al., 2023). The model we present can only explain a smooth light curve
with minimal fine structure; however, this structure is likely connected to
small mass ejections or binary interactions in the final few years of the life
of the progenitor and may not be strongly connected to the power source of the
supernova. Furthermore, The SED used in our model has a simplified treatment
of line blanketing that is calibrated to SLSNe (Nicholl et al., 2017b), and
may not be accurate for other transients. Although it is possible to switch
the SED to a blackbody if line blanketing is not expected to be strong, the
SED is also not a good approximation deep into the nebular phase on the
supernova, since the emission will start to be dominated by line emission
instead of photospheric emission (e.g., Schulze et al., 2023).

Although the magnetar-driven supernova model is versatile enough to fit light
curves of many different supernovae, the best way to determine whether a
magnetar is really the power source is to compare the nebular spectra,
polarization, and non-thermal emission from the supernovae with different
models. Within the magnetar model, this emission is form the PWN or its
interaction with the supernova ejecta. If the ejecta is optically thin, the
synchrotron and inverse Compton emission from the PWN can leak through and be
detected (Kotera et al., 2013; Metzger et al., 2014; Murase et al., 2015;
Omand et al., 2018), while if the ejecta is optically thick, these photons
will be absorbed and change temperature and electronic state of the ejecta,
both giving detectable signals (Chevalier & Fransson, 1992; Omand et al.,
2019; Omand & Jerkstrand, 2023). However, interaction with circumstellar
material can also produce non-thermal emission, polarization, and spectra with
high ionization lines, therefore detailed modeling is necessary to make any
strong conclusions.

In this work, we present a more flexible, inference-capable, publicly
available model for the light curves of magnetar-driven supernovae. The main
changes from previous models are the coupling of the magnetar energy injection
to the kinetic energy of the supernova and the addition of the magnetar
braking index as a free parameter, allowing the exploration of non-vacuum-
dipole spin down. We show that the model can reproduce the basic properties on
several phenomenologically different supernovae, and also fit four different
types of supernovae, retrieving parameters consistent with works using
separate models. This model will allow us to explore the full diversity of
these supernovae, better characterize these supernovae from just their light
curves, and make better predictions future multiwavelength emission to better
test the magnetar-driven scenario.

## Acknowledgements

We thank Claes Fransson and Steve Schulze for helpful comments and
discussions. N. Sarin is supported by a Nordita Fellowship, Nordita is funded
in part by NordForsk.

## Data Availability

Light curve data for SN2015bn and SN2007ru was obtained from the Open
Supernova Catalog (Guillochon et al., 2017) using Redback. Light curve data
for ZTF20acigmel and iPTF14gqr were downloaded from repositories released with
the corresponding papers (Perley et al., 2021; De et al., 2018). The model is
available for public use within Redback (Sarin et al., 2023).

## References

  * Akmal et al. (1998) Akmal A., Pandharipande V. R., Ravenhall D. G., 1998, Phys. Rev. C, 58, 1804
  * Angus et al. (2016) Angus C. R., Levan A. J., Perley D. A., et al. 2016, MNRAS, 458, 84
  * Arnett (1980) Arnett W. D., 1980, ApJ, 237, 541
  * Arnett (1982) Arnett W. D., 1982, ApJ, 253, 785
  * Ashton et al. (2019) Ashton G., Hübner M., Lasky P. D., et al. 2019, ApJS, 241, 27
  * Barkat et al. (1967) Barkat Z., Rakavy G., Sack N., 1967, Phys. Rev. Lett., 18, 379
  * Bellm et al. (2019) Bellm E. C., Kulkarni S. R., Graham M. J., Dekany R., et al. 2019, PASP, 131, 018002
  * Bhirombhakdi et al. (2018) Bhirombhakdi K., Chornock R., Margutti R., et al. 2018, ApJ, 868, L32
  * Blanchard et al. (2019) Blanchard P. K., Nicholl M., Berger E., et al. 2019, ApJ, 872, 90
  * Blanchard et al. (2020) Blanchard P. K., Berger E., Nicholl M., et al. 2020, ApJ, 897, 114
  * Braithwaite (2009) Braithwaite J., 2009, MNRAS, 397, 763
  * Bright et al. (2022) Bright J. S., Margutti R., Matthews D., et al. 2022, ApJ, 926, 112
  * Bucciantini et al. (2007) Bucciantini N., Quataert E., Arons J., et al. 2007, MNRAS, 380, 1541
  * Bucciantini et al. (2008) Bucciantini N., Quataert E., Arons J., et al. 2008, MNRAS, 383, L25
  * Cano et al. (2017) Cano Z., Wang S.-Q., Dai Z.-G., et al. 2017, Advances in Astronomy, 2017, 8929054
  * Cappellaro (2022) Cappellaro E., 2022, Nuovo Cimento Rivista Serie, 45, 549
  * Chatzopoulos & Wheeler (2012) Chatzopoulos E., Wheeler J. C., 2012, ApJ, 760, 154
  * Chatzopoulos et al. (2013) Chatzopoulos E., Wheeler J. C., Vinko J., et al. 2013, ApJ, 773, 76
  * Chen et al. (2016) Chen K.-J., Woosley S. E., Sukhbold T., 2016, ApJ, 832, 73
  * Chen et al. (2020) Chen K.-J., Woosley S. E., Whalen D. J., 2020, ApJ, 893, 99
  * Chen et al. (2021) Chen T. W., Brennan S. J., Wesson R., et al. 2021, arXiv e-prints, p. arXiv:2109.07942
  * Chen et al. (2023a) Chen Z. H., Yan L., Kangas T., et al. 2023a, ApJ, 943, 41
  * Chen et al. (2023b) Chen Z. H., Yan L., Kangas T., et al. 2023b, ApJ, 943, 42
  * Chevalier & Fransson (1992) Chevalier R. A., Fransson C., 1992, ApJ, 395, 540
  * Chomiuk et al. (2011) Chomiuk L., Chornock R., Soderberg A. M., et al. 2011, ApJ, 743, 114
  * Chugai & Utrobin (2022) Chugai N. N., Utrobin V. P., 2022, MNRAS, 512, L71
  * Chugai & Utrobin (2023) Chugai N. N., Utrobin V. P., 2023, Astronomy Letters, 49, 49
  * Cutler & Jones (2000) Cutler C., Jones D. I., 2000, Phys. Rev. D, 63, 024002
  * Das et al. (2003) Das C. B., Das Gupta S., Gale C., et al. 2003, Phys. Rev. C, 67, 034611
  * De et al. (2018) De K., Kasliwal M. M., Ofek E. O., et al. 2018, Science, 362, 201
  * Dessart et al. (2016) Dessart L., Hillier D. J., Woosley S., et al. 2016, MNRAS, 458, 1618
  * Dexter & Kasen (2013) Dexter J., Kasen D., 2013, ApJ, 772, 30
  * Dong et al. (2023) Dong X.-F., Liu L.-D., Gao H., et al. 2023, ApJ, 951, 61
  * Eftekhari et al. (2019) Eftekhari T., Berger E., Margalit B., et al. 2019, ApJ, 876, L10
  * Eftekhari et al. (2021) Eftekhari T., Margalit B., Omand C. M. B., et al. 2021, ApJ, 912, 21
  * Gaensler & Slane (2006) Gaensler B. M., Slane P. O., 2006, ARA&A, 44, 17
  * Gal-Yam (2012) Gal-Yam A., 2012, Science, 337, 927
  * Gal-Yam et al. (2009) Gal-Yam A., Mazzali P., Ofek E. O., et al. 2009, Nature, 462, 624
  * Gendre et al. (2013) Gendre B., Stratta G., Atteia J. L., et al. 2013, ApJ, 766, 30
  * Goldreich & Julian (1969) Goldreich P., Julian W. H., 1969, ApJ, 157, 869
  * Guillochon et al. (2017) Guillochon J., Parrent J., Kelley L. Z., et al. 2017, ApJ, 835, 64
  * Hatsukade et al. (2021) Hatsukade B., Tominaga N., Morokuma T., et al. 2021, ApJ, 911, L1
  * Heger & Woosley (2002) Heger A., Woosley S. E., 2002, ApJ, 567, 532
  * Heger et al. (2003) Heger A., Fryer C. L., Woosley S. E., Langer N., et al. 2003, ApJ, 591, 288
  * Ho et al. (2022) Ho A. Y. Q., Margalit B., Bremer M., et al. 2022, ApJ, 932, 116
  * Hosseinzadeh et al. (2022) Hosseinzadeh G., Berger E., Metzger B. D., et al. 2022, ApJ, 933, 14
  * Inserra et al. (2013) Inserra C., Smartt S. J., Jerkstrand A., et al. 2013, ApJ, 770, 128
  * Inserra et al. (2016) Inserra C., Bulla M., Sim S. A., et al. 2016, ApJ, 831, 79
  * Inserra et al. (2017) Inserra C., Nicholl M., Chen T. W., et al. 2017, MNRAS, 468, 4642
  * Ivezić et al. (2019) Ivezić Ž., Kahn S. M., Tyson J. A., Abel B., et al. 2019, ApJ, 873, 111
  * Jerkstrand et al. (2017) Jerkstrand A., Smartt S. J., Inserra C., et al. 2017, ApJ, 835, 13
  * Jiang et al. (2020) Jiang B., Jiang S., Ashley Villar V., 2020, Research Notes of the American Astronomical Society, 4, 16
  * Kasen & Bildsten (2010) Kasen D., Bildsten L., 2010, ApJ, 717, 245
  * Kasen et al. (2016) Kasen D., Metzger B. D., Bildsten L., 2016, ApJ, 821, 36
  * Kashiyama et al. (2016) Kashiyama K., Murase K., Bartos I., et al. 2016, ApJ, 818, 94
  * Khatami & Kasen (2023) Khatami D., Kasen D., 2023, arXiv e-prints, p. arXiv:2304.03360
  * Kleiser & Kasen (2014) Kleiser I. K. W., Kasen D., 2014, MNRAS, 438, 318
  * Kotera et al. (2013) Kotera K., Phinney E. S., Olinto A. V., 2013, MNRAS, 432, 3228
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration Abell P. A., Allison J., et al. 2009, arXiv e-prints, p. arXiv:0912.0201
  * Lander & Jones (2018) Lander S. K., Jones D. I., 2018, MNRAS, 481, 4169
  * Lander & Jones (2020) Lander S. K., Jones D. I., 2020, MNRAS, 494, 4838
  * Lasky et al. (2017) Lasky P. D., Leris C., Rowlinson A., et al. 2017, ApJ, 843, L1
  * Lattimer & Schutz (2005) Lattimer J. M., Schutz B. F., 2005, ApJ, 629, 979
  * Law et al. (2009) Law N. M., Kulkarni S. R., Dekany R. G., et al. 2009, PASP, 121, 1395
  * Law et al. (2019) Law C. J., Omand C. M. B., Kashiyama K., et al. 2019, ApJ, 886, 24
  * Lee (2019) Lee C.-H., 2019, ApJ, 875, 121
  * Lee (2020) Lee C.-H., 2020, Astronomische Nachrichten, 341, 651
  * Leloudas et al. (2015a) Leloudas G., Schulze S., Krühler T., et al. 2015a, MNRAS, 449, 917
  * Leloudas et al. (2015b) Leloudas G., Patat F., Maund J. R., et al. 2015b, ApJ, 815, L10
  * Leloudas et al. (2017) Leloudas G., Maund J. R., Gal-Yam A., et al. 2017, ApJ, 837, L14
  * Levan et al. (2014) Levan A. J., Tanvir N. R., Starling R. L. C., et al. 2014, ApJ, 781, 13
  * Liu et al. (2022) Liu J.-F., Zhu J.-P., Liu L.-D., et al. 2022, ApJ, 935, L34
  * Lunnan et al. (2014) Lunnan R., Chornock R., Berger E., et al. 2014, ApJ, 787, 138
  * Lyne et al. (2015) Lyne A. G., Jordan C. A., Graham-Smith F., et al. 2015, MNRAS, 446, 857
  * MacFadyen & Woosley (1999) MacFadyen A. I., Woosley S. E., 1999, ApJ, 524, 262
  * Metzger (2019) Metzger B. D., 2019, Living Reviews in Relativity, 23, 1
  * Metzger et al. (2014) Metzger B. D., Vurm I., Hascoët R., et al. 2014, MNRAS, 437, 703
  * Milisavljevic et al. (2013) Milisavljevic D., Soderberg A. M., Margutti R., et al. 2013, ApJ, 770, L38
  * Mondal et al. (2020) Mondal S., Bera A., Chandra P., et al. 2020, MNRAS, 498, 3863
  * Moriya et al. (2018) Moriya T. J., Nicholl M., Guillochon J., 2018, ApJ, 867, 113
  * Moriya et al. (2022) Moriya T. J., Murase K., Kashiyama K., et al. 2022, MNRAS, 513, 6210
  * Murase et al. (2015) Murase K., Kashiyama K., Kiuchi K., et al. 2015, ApJ, 805, 82
  * Murase et al. (2021) Murase K., Omand C. M. B., Coppejans D. L., et al. 2021, MNRAS, 508, 44
  * Nadyozhin (1994) Nadyozhin D. K., 1994, ApJS, 92, 527
  * Nakauchi et al. (2013) Nakauchi D., Kashiyama K., Suwa Y., et al. 2013, ApJ, 778, 67
  * Nicholl (2021) Nicholl M., 2021, Astronomy and Geophysics, 62, 5.34
  * Nicholl et al. (2013) Nicholl M., Smartt S. J., Jerkstrand A., et al. 2013, Nature, 502, 346
  * Nicholl et al. (2016a) Nicholl M., Berger E., Smartt S. J., et al. 2016a, ApJ, 826, 39
  * Nicholl et al. (2016b) Nicholl M., Berger E., Margutti R., et al. 2016b, ApJ, 828, L18
  * Nicholl et al. (2017a) Nicholl M., Berger E., Margutti R., et al. 2017a, ApJ, 835, L8
  * Nicholl et al. (2017b) Nicholl M., Guillochon J., Berger E., 2017b, ApJ, 850, 55
  * Nicholl et al. (2018) Nicholl M., Blanchard P. K., Berger E., et al. 2018, ApJ, 866, L24
  * Omand & Jerkstrand (2023) Omand C. M. B., Jerkstrand A., 2023, A&A, 673, A107
  * Omand et al. (2018) Omand C. M. B., Kashiyama K., Murase K., 2018, MNRAS, 474, 573
  * Omand et al. (2019) Omand C. M. B., Kashiyama K., Murase K., 2019, MNRAS, 484, 5468
  * Ørum et al. (2020) Ørum S. V., Ivens D. L., Strandberg P., et al. 2020, A&A, 643, A47
  * Ostriker & Gunn (1969) Ostriker J. P., Gunn J. E., 1969, ApJ, 157, 1395
  * Parthasarathy et al. (2020) Parthasarathy A., Johnston S., Shannon R. M., et al. 2020, MNRAS, 494, 2012
  * Pastorello et al. (2010) Pastorello A., Smartt S. J., Botticella M. T., et al. 2010, ApJ, 724, L16
  * Perley et al. (2021) Perley D. A., Ho A. Y. Q., Yao Y., et al. 2021, MNRAS, 508, 5138
  * Poidevin et al. (2022) Poidevin F., Omand C. M. B., Pérez-Fournon I., et al. 2022, MNRAS, 511, 5948
  * Poidevin et al. (2023) Poidevin F., Omand C. M. B., Könyves-Tóth R., et al. 2023, MNRAS, 521, 5418
  * Pursiainen et al. (2022) Pursiainen M., Leloudas G., Paraskeva E., et al. 2022, arXiv e-prints, p. arXiv:2202.01635
  * Pursiainen et al. (2023) Pursiainen M., Leloudas G., Cikota A., et al. 2023, arXiv e-prints, p. arXiv:2301.08111
  * Reboul-Salze et al. (2021) Reboul-Salze A., Guilet J., Raynaud R., Bugli M., 2021, A&A, 645, A109
  * Romero-Shaw et al. (2020) Romero-Shaw I. M., Talbot C., Biscoveanu S., et al. 2020, MNRAS, 499, 3295
  * Sahu et al. (2009) Sahu D. K., Tanaka M., Anupama G. C., et al. 2009, ApJ, 697, 676
  * Saito et al. (2020) Saito S., Tanaka M., Moriya T. J., et al. 2020, ApJ, 894, 154
  * Sarin & Lasky (2021) Sarin N., Lasky P. D., 2021, General Relativity and Gravitation, 53, 59
  * Sarin et al. (2018) Sarin N., Lasky P. D., Sammut L., et al. 2018, Phys. Rev. D, 98, 043011
  * Sarin et al. (2020a) Sarin N., Lasky P. D., Ashton G., 2020a, Phys. Rev. D, 101, 063021
  * Sarin et al. (2020b) Sarin N., Lasky P. D., Ashton G., 2020b, MNRAS, 499, 5986
  * Sarin et al. (2022) Sarin N., Omand C. M. B., Margalit B., et al. 2022, MNRAS, 516, 4949
  * Sarin et al. (2023) Sarin N., et al., 2023, Redback, https://github.com/nikhil-sarin/redback, doi:https://redback.readthedocs.io/en/latest/
  * Sawada et al. (2022) Sawada R., Kashiyama K., Suwa Y., 2022, ApJ, 927, 223
  * Schulze et al. (2018) Schulze S., Krühler T., Leloudas G., et al. 2018, MNRAS, 473, 1258
  * Schulze et al. (2023) Schulze S., Fransson C., Kozyreva A., et al. 2023, arXiv e-prints, p. arXiv:2305.05796
  * Shapiro & Teukolsky (1983) Shapiro S. L., Teukolsky S. A., 1983, Black holes, white dwarfs and neutron stars. The physics of compact objects, doi:10.1002/9783527617661.
  * Shetty et al. (2007) Shetty D. V., Yennello S. J., Souliotis G. A., 2007, Phys. Rev. C, 75, 034602
  * Smith (2014) Smith N., 2014, ARA&A, 52, 487
  * Speagle (2020) Speagle J. S., 2020, MNRAS, 493, 3132
  * Stritzinger & Moriya (2018) Stritzinger M., Moriya T., 2018, arXiv e-prints, p. arXiv:1801.06643
  * Stroh et al. (2021) Stroh M. C., Terreran G., Coppejans D. L., et al. 2021, ApJ, 923, L24
  * Sun et al. (2022) Sun L., Xiao L., Li G., 2022, MNRAS, 513, 4057
  * Suwa et al. (2015) Suwa Y., Yoshida T., Shibata M., et al. 2015, MNRAS, 454, 3073
  * Suzuki & Maeda (2017) Suzuki A., Maeda K., 2017, MNRAS, 466, 2633
  * Suzuki & Maeda (2019) Suzuki A., Maeda K., 2019, ApJ, 880, 150
  * Suzuki & Maeda (2021) Suzuki A., Maeda K., 2021, ApJ, 908, 217
  * Taddia et al. (2019) Taddia F., Sollerman J., Fremling C., et al. 2019, A&A, 621, A71
  * Tanaka et al. (2017) Tanaka M., Maeda K., Mazzali P. A., et al. 2017, ApJ, 837, 105
  * Thompson et al. (2004) Thompson T. A., Chang P., Quataert E., 2004, ApJ, 611, 380
  * Tonry (2011) Tonry J. L., 2011, PASP, 123, 58
  * Tonry et al. (2018) Tonry J. L., Denneau L., Heinze A. N., et al. 2018, PASP, 130, 064505
  * Villar et al. (2017) Villar V. A., Berger E., Metzger B. D., et al. 2017, ApJ, 849, 70
  * Wang et al. (2015) Wang S. Q., et al., 2015, ApJ, 799, 107
  * Wang et al. (2016) Wang L.-J., Han Y.-H., Xu D., et al. 2016, ApJ, 831, 41
  * Watts et al. (2016) Watts A. L., Andersson N., Chakrabarty D., et al. 2016, Reviews of Modern Physics, 88, 021001
  * West et al. (2023) West S. L., Lunnan R., Omand C. M. B., et al. 2023, A&A, 670, A7
  * Woosley (2010) Woosley S. E., 2010, ApJ, 719, L204
  * Woosley (2017) Woosley S. E., 2017, ApJ, 836, 244
  * Woosley et al. (2007) Woosley S. E., Blinnikov S., Heger A., 2007, Nature, 450, 390
  * Worley et al. (2008) Worley A., Krastev P. G., Li B.-A., 2008, ApJ, 685, 390
  * Yoshida et al. (2016) Yoshida T., Umeda H., Maeda K., et al. 2016, MNRAS, 457, 351
  * Yu et al. (2013) Yu Y.-W., Zhang B., Gao H., 2013, ApJ, 776, L40
  * Zhang & Mészáros (2001) Zhang B., Mészáros P., 2001, ApJ, 552, L35

##  Appendix A Posteriors from Case Studies

Here we show the full posteriors for the real supernovae examined as case
studies, SN 2015bn (Figure 9), SN 2007ru (Figure 10), ZTF20acigmel (Figure
11), and iPTF14gqr (Figure 12).

Figure 9: Corner plot for SN 2015bn, an SLSN. Figure 10: Corner plot for SN
2007ru, an SN Ic-BL. Figure 11: Corner plot for ZTF20acigmel (the Camel), an
FBOT. Figure 12: Corner plot for iPTF14gqr, a USSN.