# Self-Supervised Neuron Segmentation with Multi-Agent Reinforcement Learning

Yinda Chen1,2 Wei Huang1 Shenglong Zhou1 Qi Chen1 Zhiwei
Xiong1,2,111Corresponding Author 1University of Science and Technology of
China
2Institute of Artificial Intelligence, Hefei Comprehensive National Science
Center {cyd0806, weih527, slzhou96<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The performance of existing supervised neuron segmentation methods is highly
dependent on the number of accurate annotations, especially when applied to
large scale electron microscopy (EM) data. By extracting semantic information
from unlabeled data, self-supervised methods can improve the performance of
downstream tasks, among which the mask image model (MIM) has been widely used
due to its simplicity and effectiveness in recovering original information
from masked images. However, due to the high degree of structural locality in
EM images, as well as the existence of considerable noise, many voxels contain
little discriminative information, making MIM pretraining inefficient on the
neuron segmentation task. To overcome this challenge, we propose a decision-
based MIM that utilizes reinforcement learning (RL) to automatically search
for optimal image masking ratio and masking strategy. Due to the vast
exploration space, using single-agent RL for voxel prediction is impractical.
Therefore, we treat each input patch as an agent with a shared behavior
policy, allowing for multi-agent collaboration. Furthermore, this multi-agent
model can capture dependencies between voxels, which is beneficial for the
downstream segmentation task. Experiments conducted on representative EM
datasets demonstrate that our approach has a significant advantage over
alternative self-supervised methods on the task of neuron segmentation. Code
is available at https://github.com/ydchen0806/dbMiM.

##  1 Introduction

Figure 1: A comparison of the reconstruction effectiveness of our proposed
method with MAE. The first column shows the original EM image, the second
column shows the image after masking 85% of the voxels, the third column shows
the reconstruction using the MAE method, and the fourth column shows the
reconstruction using our method.

Neuron segmentation is a crucial task for neuroscientists that allows for the
analysis of the distribution and morphology of neurons, providing valuable
insights into the connectomics research Sheridan et al. (2022); Krasowski et
al. (2017). Electron microscopy (EM) is the mainstream method for accurately
identifying neural structures, but the dense nature of neurons and the
presence of artifacts and deformations in EM images make the labeling process
costly and decrease the credibility of existing annotation data Deng et al.
(2022); Zhou et al. (2019); Chen et al. (2023a). Therefore, fully supervised
neuron segmentation methods meet great challenges, especially when applied to
large scale EM data.

Self-supervised methods have emerged as a solution to the limitations of fully
supervised methods, which can be roughly divided into two categories:
contrastive learning-based approach and mask image model (MIM)-based approach.
The former requires a large number of positive and negative samples Chen et
al. (2020); Grill et al. (2020); You et al. (2022); Chen et al. (2023b) and
relies heavily on data augmentation Caron et al. (2021), making it a high-cost
option for 3D biomedical images. The latter aims to learn useful structural
information in images by masking and recovering certain voxels, which has been
recently applied to pretraining biomedical images, showing improvements in
downstream tasks Zhou et al. (2022); Tang et al. (2022); Huang et al. (2022a).
However, the highly localized and structured nature of EM data, as well as the
existence of considerable noise, make it inefficient to directly employ the
existing MIM in extracting useful information for neuron segmentation. It has
also been observed that the masking ratio and masking strategy of MIM are
highly sensitive and the optimal ones vary greatly across different datasets.
Adjusting these configurations to train large models requires significant
manual efforts and resources.

In this paper, targeting the neuron segmentation task, we propose a novel
decision-based MIM relying on multi-agent reinforcement learning (MARL)
Littman (1994) for automatically selecting the appropriate masking ratio and
masking strategy, which consists of a target network and a policy network. Our
approach partitions the input EM volume into patches and treats each patch as
a basic control unit. The overall multi-agent task is modeled as a search
process for patch masking strategies, where the action space for each patch is
to either keep the original voxels or mask them. The feedback of the target
network, in the form of the reconstruction loss, serves as the team reward
signal for guiding the policy network to adaptively learn masking strategies
that are beneficial for the pretraining task Foerster et al. (2016). To
improve the stability of training and achieve optimal joint decision-making
for the entire volume, all agent networks share parameters and are trained in
parallel. Furthermore, we introduce the HOG feature as an additional
reconstruction loss to enable the target network to learn more structure
information. Finally, following a UNETR decoder design Hatamizadeh et al.
(2022), we add a segmentation head to the pretrained target network in the
finetuning stage. Experimental results in Figure 1 show that our decision-
based MIM achieves clearer reconstruction results than the original MAE He et
al. (2022) in the pretraining phase.

Overall, our main contribution lies in the following aspects:

1) We propose an efficient self-supervised method, named decision-based MIM,
for EM neuron segmentation using unlabeled EM data. To the best of our
knowledge, it is the first effort that large-scale transformer pretraining is
conducted on this task.

2) We propose a MARL-based approach for searching the optimal masking ratio
and masking strategy by treating each patch as an agent with a shared policy,
effectively exploring the search space and capturing dependencies between
voxels.

3) We introduce the HOG feature as an additional reconstruction loss of our
decision-based MIM, improving the convergence speed of network training and
the performance of the downstream segmentation task.

4) We comprehensively demonstrate the effectiveness of our proposed method on
two representative EM datasets, especially against alternative self-supervised
methods on the task of neuron segmentation.

##  2 Related Work

###  2.1 Neuron Instance Segmentation

In the field of EM image processing, neuron instance segmentation is an
important task. Turaga et al. (2010) first proposed a convolutional neural
network based on affinity generation, followed by clustering voxels in the
affinity graph into instances based on post-processing methods such as
watershed and LMC Beier et al. (2017). In recent years, there have been more
advanced networks for the affinity-based approach. Funke et al.Funke et al.
(2018) introduced the MALIS lossBriggman et al. (2009) during the training
process to encourage the network to generate correct topological segmentation.
Huang et al. (2022b) introduced an embedding pyramid module to simulate
affinity at different scales. Liu et al. (2022) incorporated both embedding
and affinity information and combined it with the graph neural network to
further improve the distinguishability of adjacent objects in the feature
space. However, due to the anisotropic resolution of 3D EM images in lateral
and axial directions, the usage of transformer structures remains unexplored
in the field of neuron segmentation. In this paper, we use an affinity-based
setup and upgrade the decoder of UNETR Hatamizadeh et al. (2022) to adapt to
the anisotropic EM features.

Figure 2: Our proposed network architecture is divided into two main
components: (a) the decision-based Mask Image Model (MIM) pretraining process,
which employs our proposed decision module to select appropriate patches for
masking and then utilizes a 3D Vision Transformer (ViT) encoder to encode the
visible patches. The resulting tokens are then passed through a lightweight
decoder to reconstruct the original voxels and histograms of oriented gradient
(HOG) features. (b) The fine-tuning process for the downstream segmentation
task utilizes the encoder weights from the pretraining process and adds a
UNETR segmentation head to output the affinity map. The final segmentation
results are obtained through post-processing methods such as waterz.

###  2.2 Mask Image Model (MIM)

The MIM is an important branch of self-supervised learning. Masked
autoencoders (MAE) He et al. (2022) used an asymmetric encoding-decoding
structure, encoding only the unmasked patches and using a lightweight decoder
to recover the masked patches, which greatly reduced the resources required
for computation and quickly became a mainstream structure for MIM.
Feichtenhofer et al. (2022) and Tong et al. (2022) separately validated the
effectiveness of MAE on video datasets and proved that higher image masking
ratios can be used in 3D datasets. Zhou et al. (2022) was the first to
introduce MAE to the medical image field. Bachmann et al. (2022) introduced a
multi-modal multi-task adaptation to MAE, resulting in better performance than
the original MAE. Gao et al. (2022) proposed a multi-scale mixed convolution
to encode images and achieved improved results on fine-grained downstream
tasks such as image segmentation. Wei et al. (2022) focused on the
reconstruction target of the decoder and found that reconstructing artificial
features such as HOG and SIFT can facilitate the network to classify and
localize objects.

In summary, existing MIM-based methods mainly focus on the design of the
encoder architecture, prediction head, and prediction target. Although many
works have demonstrated the impact of masking strategies on downstream tasks,
there has been little research along this line, and the latest method also
requires a search-based masking strategy with a fixed masking ratio Bandara et
al. (2022). This paper proposes a novel MARL-based approach for adaptively
learning the optimal masking ratio and masking strategy, making the pretrained
model more robust and achieving better performance on downstream tasks.

###  2.3 Multi-Agent Reinforcement Learning (MARL)

In the field of deep RL, multiple agents working together can improve the
efficiency and robustness of the model due to the limited observation and
action space of a single agent. Given the complexity of computer vision tasks,
MARL is often used to interact with the common environment to make decisions.
Liao et al. (2020) proposed a method for iteratively refining medical image
segmentation using interactive MARL, where rough image segmentation is
initially provided and the network is iteratively refined through user
feedback in the form of rewards until the segmentation is sufficient. Lin et
al. (2021a) proposed a method for augmenting images through the use of blocks,
with each block acting as an agent and working together to produce optimal
data augmentation. Still, the application of MARL in computer vision suffers
from the large decision space and difficulty in obtaining rewards. In this
paper, to reduce the complexity of the state and the difficulty of searching
for RL policies, we first segment the image into patches through a
transformer-encoder and treat each patch as an agent. The action space for
each patch is limited to only two options, masking or keeping the original
voxels, and the rewards are obtained through the reconstruction loss of MAE.
Therefore, the masking decision can be naturally modeled as a Markov process
and optimized through MARL.

##  3 Proposed Method

Our decision-based MIM consists of two stages of training: pretraining and
finetuning. Figure 2 illustrates the overall flow of the network training. In
this section, we will first introduce some basic theories of Vision
Transformer (ViT) and MARL, and then explain in detail the encoders, decoders,
and loss functions of the two stages, as well as specific MARL modeling
methods.

###  3.1 Encoder-Decoder Design

We use ViT Dosovitskiy et al. (2021) as the backbone architecture for
decision-based MIM pretraining and downstream segmentation tasks. To represent
high-dimensional data in a ViT, we must transform it into a sequence of
patches. Given an input 3D volume $\mathbf{x}\in\mathbb{R}^{H\times W\times
D\times C}$, where $C$ is the number of channels and $(H,W,D)$ is the
resolution, we reshape it into a sequence of flattened 3D patches
$\mathbf{x}_{p}\in\mathbb{R}^{N\times(P^{3}\cdot C)}$. The patch resolution is
given by $(P/4,P,P)$, and the number of patches is calculated as
$N=\frac{4HWD}{P^{3}}$. These patches are then transformed into patch
embeddings via a trainable linear projection.

Consistent with the MAE setup, we divide the image patches into visible and
masked groups. The encoder in the MAE and ViT architectures processes only the
visible blocks. To enhance the performance of our decoder, we incorporated a
histogram-of-oriented-gradients (HOG) feature Dalal and Triggs (2005), which
has been shown to improve pretraining. This is achieved by providing the
decoder with various markers, including a patch representation from the
encoder and learnable position embeddings. By incorporating positional
embedding in all input markers, we enable the decoder to simultaneously
recover both the HOG feature and the original voxels, resulting in superior
performance compared to the MAE method.

In our pretraining process, we utilize a loss function that combines both the
mean squared error (MSE) loss for reconstructing the original voxels and the
HOG loss for recovering the HOG features. The HOG feature can be calculated
using the following equation

$\text{HOG}_{i,j}=\frac{\sum_{x\in S_{i,j}}{w(x)g(x)}}{\sum_{x\in
S_{i,j}}{w(x)}},$ (1)

where HOGi,j is the histogram of oriented gradients for the cell located at
position $(i,j)$, $S_{i,j}$ is the set of voxels in the cell $(i,j)$, $w(x)$
is a weighting function that assigns a weight to each voxel $x$ in the cell,
and $g(x)$ is the gradient orientation of the voxel $x$.

Our overall loss function can be expressed as

$\mathcal{L}_{pretrain}=\lambda_{1}\mathcal{L}_{MSE}+\lambda_{2}\mathcal{L}_{HOG},$
(2)

where $\lambda_{1}$ and $\lambda_{2}$ denote the weights assigned to the MSE
and HOG losses, respectively. We set $\lambda_{1}$ to 0.1 and $\lambda_{2}$ to
1.

###  3.2 Decision Module

In our proposed decision module, we model the masking strategy of image
patches in MAE as a multi-agent cooperative decision-making problem and adopt
a multi-agent reinforcement learning method to solve it. Here, we will
introduce in detail the modeling methods of model states, observations,
actions, the design of multi-agent team rewards, and the learning method used
to update the policy network.

As shown in Figure 3, given the original input batch of image $\mathbf{x}$,
according to the ViT setup, we divide the image into equal-sized and non-
overlapping patches. Our MARL policy aims to determine the overall joint
masking policy based on the current input state and the observation of each
agent. Our basic setup is as follows.

####  3.2.1 State.

The enhancement policy for each patch is closely related to the information of
the context, so the decision-making process of MARL requires perceiving the
semantic information of the image rather than directly inputting the patches
of the raw image. In order to ensure the consistency and convergence of
training, we use the target network, i.e. ViT, as the backbone to extract deep
semantic features of the image. The global state at time step $t$ is denoted
as $S^{t}$ and is visible to all agents.

####  3.2.2 Observation.

In addition to capturing global information $S^{t}$, each agent needs to make
a masking policy decision based on its own local observation. In general, the
observation in MARL tasks is often a part of the global state. Considering
these factors, we take the deep feature ViT$(P_{i})$ of patch $P_{i}$ as the
current observation value $O_{i}^{t}$ for the i-th agent. The feature
extractor for local features is the same as the one for global features, both
using the ViT backbone.

####  3.2.3 Action.

The action of the i-th agent aims to output whether patch $P_{i}$ needs to be
masked. We define the action of the i-th agent as a vector $A_{i}$. The joint
action space can be represented as $A=\\{A_{1},A_{2},...,A_{N}\\}$, where $N$
represents the total number of patches. Since the action space only has two
possibilities, masking or keeping the original voxels, the dimension of
$A_{i}$ is 2. Given the current state $S^{t}$ and observation $O_{i}^{t}$,
each agent $i$ will determine an action $a_{i}(S^{t},O^{t}_{i})\in A_{i}$
based on the current policy. The final output is the global joint action
$\boldsymbol{a^{t}}=\\{a^{t}_{1}\cup a^{t}_{2},...,\cup a^{t}_{N}\\}$. After
all patches have taken their corresponding actions through the decision
policy, the time step is updated to $t=t+1$ and we obtain the enhanced volume
through the decision module.

####  3.2.4 Rewards.

Rewards are introduced in our MARL decision-making process in order to guide
the agents to learn expected behaviors that will improve the main task’s
performance through more reasonable masking ratios and masking strategies,
allowing the target net to better learn semantic information in the volume.
Previous works Lin et al. (2021b); Zhang et al. (2019) attempted to increase
the training loss of the target network based on rewards in order to generate
deeper, more difficult-to-learn features. Inspired by these works, we refine
the reward design based on the MAE pretraining paradigm. By comparing the loss
difference between the data obtained from the masked data
$\mathbf{x}\boldsymbol{\cdot a^{t-1}}$ generated by the previous time step’s
decision module and the data obtained from the current decision module
$\mathbf{x}\boldsymbol{\cdot a^{t}}$, we compute the reward for the MARL
policy. This encourages higher training loss during the MARL decision-making
process, as shown by the equation

$r^{t}=\mathcal{L}_{pretrain}(\phi(\mathbf{x}\boldsymbol{\cdot
a^{t}}))-\mathcal{L}_{pretrain}(\phi(\mathbf{x}\boldsymbol{\cdot a^{t-1}})).$
(3)

In the above equation, $\mathcal{L}_{pretrain}$ represents the reconstruction
loss generated by the target network, and $\phi$ denotes the target network.
The accumulated reward of one sequence is

$R^{t}=\sum_{i=t-T+1}^{t}\gamma^{i-1}\bar{r}^{i},$ (4)

where $T$ represents the desired time step length to be calculated, the
discount factor $\gamma$ takes a value in (0, 1], and $\bar{r}^{t}$ is the
mean rewards at time $t$.

####  3.2.5 Policy Learning.

Considering that the action space for MARL decisions is discrete, we utilize
the widely used Asynchronous Advantage Actor-Critic (A2C) algorithm Mnih et
al. (2016) to perform MARL policy learning. Since the search space for actions
is not large, we use simple convolution and pooling modules to adjust the
Actor and Critic networks to the structure shown in Figure 3. The policy
network is divided into an Actor and a Critic, which are adapted to the RL
training algorithm. The Actor network learns a discrete control policy
$\pi(a^{t}_{i}|S^{t},O^{t}_{i})$, while the Critic network aims to estimate
the value of the state $V_{\pi}(S^{t})$. We model the centralized action value
function $Q$, which takes in the state information $S$ and the actions of all
agents, and outputs a $Q$ value for the team, given by

$Q^{\boldsymbol{\pi}}(S^{t},\boldsymbol{a^{t}})=E_{\boldsymbol{\pi}}\left[R_{t}\mid
S^{t},a^{t}_{1},\cdots,a^{t}_{N}\right],$ (5)

where $\boldsymbol{a^{t}}$ represents the joint action of all agents, defined
as $\boldsymbol{a}=\\{a_{i},\cdots,a_{N}\\}$, and $R_{t}$ is the long-term
discounted reward, given by equation 4. The advantage function on the policy
is then given by

$A^{\pi}(S^{t},\boldsymbol{a^{t}})=Q^{\pi}(S^{t},\boldsymbol{a^{t}})-V^{\pi}(S^{t}),$
(6)

where $A^{\pi}(S^{t},\boldsymbol{a^{t}})$ is the advantage of taking action
$\boldsymbol{a^{t}}$ given state $S^{t}$ at time step $t$, $V^{\pi}(S^{t})$ is
the current state estimate output by Critic. It indicates that the actual
accumulated reward is independent of the state and reduces the variance of the
gradient. We use $\theta_{p}$ and $\theta_{v}$ to denote the parameters of the
Actor and Critic, respectively. The squared value of the dominance function
$A^{\pi}$ is taken as the loss function to update $\theta_{v}$ as

$\mathcal{L}(\theta_{v})=A^{\pi}(S^{t},\boldsymbol{a^{t}})^{2}.$ (7)

To further achieve cooperative capability, the loss function of the updated
Actor $\theta_{p}$ is defined as

$\mathcal{L}(\theta_{p})=-\log\pi_{\theta}(\boldsymbol{a^{t}}\mid
S^{t})A^{\pi}(S^{t},\boldsymbol{a^{t}}),$ (8)

where $\pi_{\theta}(\boldsymbol{a^{t}}\mid S^{t})$ is the Actor output, that
is, the probability of taking each action $a^{t}_{i}$. The Actor and Critic
are jointly trained in an end-to-end manner.

Figure 3: The framework of our proposed decision module. The first step in the
process is to extract features from the encoder of the target net. In the
second step, the policy network is used to decide whether or not the current
patch needs to be masked. The output of the joint action constitutes the final
decision.

###  3.3 Neuron Instance Segmentation Method

The UNETR model is specifically designed for 3D segmentation tasks, such as
organ segmentation. It is built upon the pretrained ViT encoder from decision-
based MIM and a randomly-initialized convolutional decoder. The UNETR
architecture is inspired by the U-Net model Ronneberger et al. (2015), with
the inclusion of skip connections between the encoder and decoder at multiple
resolutions. The input to the UNETR decoder is a sequence of representations,
denoted as ${z_{3},z_{6},z_{9},z_{12}}$, which are reshaped to restore the
spatial dimension of $\frac{H}{P}\times\frac{W}{P}\times\frac{D}{P}\times K$,
where $K$ is the feature channel.

Starting from the deepest feature, $z_{12}$, each representation is processed
through a varying number of deconvolutional blocks to increase its resolution
by a specific factor. For example, $z_{12}$ and $z_{9}$ are upsampled by a
factor of $\times 2$, $z_{6}$ is upsampled by a factor of $\times 4$, and
$z_{3}$ is upsampled by a factor of $\times 8$. Then, representations at the
same spatial resolution, such as $z_{12}$ and $z_{9}$, are concatenated and
further upsampled to match the shape of a shallower feature. This process of
concatenation and upsampling is repeated until the full resolution of the
input is restored. Finally, the output layer combines the upsampled feature
and the original full-resolution input to predict the segmentation map.

In this paper, for the segmentation task of EM neurons, we use a combination
of affinity map-based and post-processing approaches. We first create an
affinity graph based on the voxel affinity. This graph serves as the
foundation for our post-processing techniques, which utilize both waterz Funke
et al. (2018) and LMC Beier et al. (2017) to cluster the affinity map and
produce the final neuron segmentation results. Affinity-based methods have
proven to be highly effective in accurately segmenting and analyzing these
complex structures.

##  4 Experiments

###  4.1 Training Strategy

Our strategy consists of two phases: pretraining and fine-tuning. In the
pretraining phase, to improve the training efficiency of the framework, we
first pretrain decision-based MIM for 100k iterations, synchronously updating
the parameters of MAE and the policy network in the decision module, and then
fix the parameters of policy network and only update the parameters of MAE for
another 100k iterations. In the fine-tuning phase, we load the pretrained ViT
weights into the model for the downstream task and train for 200k iterations.

We use the Adam optimizer in both the pretraining and fine-tuning phases, with
$\beta_{1}=0.9,\beta_{2}=0.999$. The only difference lies in the pretraining
process, where we set the learning rate to 0.0001 and perform batch size 16
pretraining on 8 RTX 3090s. In the fine-tuning phase, we adopt a Layer-wise
Learning Rate Decay (LLRD) training method, which adjusts the learning rate
layer by layer during training. We set the learning rate of the last layer’s
parameters to 0.001 and the learning rate of the previous layer’s parameters
to 0.95 times the learning rate of the next layer’s parameters. We conduct
batch size 8 fine-tuning on 2 RTX 3090s.

###  4.2 Datasets and Evaluation Metrics

|  |  waterz Funke et al. (2018) |  LMC Beier et al. (2017)
---|---|---|---
Dataset | Method |  VOIS $\downarrow$ |  VOIM $\downarrow$ |  VOI $\downarrow$ |  ARAND $\downarrow$ |  VOIS $\downarrow$ |  VOIM $\downarrow$ |  VOI $\downarrow$ |  ARAND $\downarrow$
|  superhuman Lee et al. (2017) | 0.443 | 0.320 | 0.763 | 0.132 | 0.578 | 0.272 | 0.850 | 0.131
|  MALA Funke et al. (2018) | 0.478 | 0.627 | 1.105 | 0.459 | 0.574 | 0.303 | 0.878 | 0.146
|  UNETR Hatamizadeh et al. (2022) | 0.495 | 0.359 | 0.854 | 0.153 | 0.630 | 0.278 | 0.908 | 0.141
|  MAE He et al. (2022)+UNETR | 0.463 | 0.324 | 0.787 | 0.144 | 0.563 | 0.268 | 0.831 | 0.139
|  Dino Caron et al. (2021)+UNETR | 0.482 | 0.351 | 0.833 | 0.150 | 0.580 | 0.287 | 0.867 | 0.147
CREMI A | ours+UNETR | 0.411 | 0.331 | 0.743 | 0.131 | 0.537 | 0.260 | 0.797 | 0.134
|  superhuman Lee et al. (2017) | 0.668 | 0.409 | 1.076 | 0.082 | 0.959 | 0.232 | 1.191 | 0.060
|  MALA Funke et al. (2018) | 0.797 | 0.512 | 1.309 | 0.147 | 1.060 | 0.264 | 1.324 | 0.084
|  UNETR Hatamizadeh et al. (2022) | 0.937 | 0.397 | 1.333 | 0.103 | 1.194 | 0.230 | 1.424 | 0.116
|  MAE He et al. (2022)+UNETR | 0.776 | 0.391 | 1.167 | 0.099 | 0.994 | 0.224 | 1.218 | 0.102
|  Dino Caron et al. (2021)+UNETR | 0.916 | 0.396 | 1.312 | 0.104 | 1.167 | 0.229 | 1.396 | 0.111
CREMI B | ours+UNETR | 0.642 | 0.381 | 1.023 | 0.092 | 0.893 | 0.220 | 1.113 | 0.097
|  superhuman Lee et al. (2017) | 0.943 | 0.385 | 1.328 | 0.134 | 1.176 | 0.260 | 1.436 | 0.125
|  MALA Funke et al. (2018) | 0.901 | 0.621 | 1.522 | 0.169 | 1.137 | 0.289 | 1.426 | 0.127
|  UNETR Hatamizadeh et al. (2022) | 0.996 | 0.423 | 1.419 | 0.158 | 1.417 | 0.236 | 1.653 | 0.148
|  MAE He et al. (2022)+UNETR | 1.001 | 0.298 | 1.299 | 0.120 | 1.272 | 0.214 | 1.486 | 0.113
|  Dino Caron et al. (2021)+UNETR | 1.011 | 0.412 | 1.423 | 0.156 | 1.364 | 0.234 | 1.598 | 0.146
CREMI C | ours+UNETR | 0.925 | 0.276 | 1.201 | 0.107 | 1.194 | 0.204 | 1.398 | 0.112
|  superhuman Lee et al. (2017) | 0.721 | 0.295 | 1.016 | 0.187 | 0.770 | 0.343 | 1.113 | 0.110
|  MALA Funke et al. (2018) | 0.734 | 0.385 | 1.119 | 0.305 | 0.832 | 0.357 | 1.189 | 0.108
|  UNETR Hatamizadeh et al. (2022) | 0.908 | 0.337 | 1.245 | 0.316 | 1.007 | 0.340 | 1.347 | 0.134
|  MAE He et al. (2022)+UNETR | 0.791 | 0.306 | 1.097 | 0.254 | 0.888 | 0.298 | 1.186 | 0.120
|  Dino Caron et al. (2021)+UNETR | 0.889 | 0.329 | 1.218 | 0.298 | 1.001 | 0.314 | 1.315 | 0.129
AC4 | ours+UNETR | 0.647 | 0.285 | 0.931 | 0.243 | 0.795 | 0.284 | 1.079 | 0.113

Table 1: Results on the CREMI dataset, VOIS represents split error, VOIM
represents merge error, and VOI is the sum of the two. The final segmentation
results are generated by using two classic post-processing methods, waterz and
LMC.

##### FAFB.

The Full Adult Fly Brain (FAFB) dataset Zheng et al. (2018) is a highly
valuable resource for neuroinformatics research, offering a comprehensive and
detailed view of the neural architecture of the Drosophila melanogaster (fruit
fly) brain. With a size of approximately 40 terabytes, this dataset features
high-resolution 3D images with a resolution of approximately 4 nanometers per
pixel, as well as manually annotated segmentation data identifying various
brain structures such as neurons, glial cells, blood vessels, synapses, and
other neuropil regions. In our work, we utilize the FAFB dataset as a key
component in our pretraining process. To optimize the efficiency of our model,
we first downsample the original dataset by a factor of 4, carefully curating
a selection of 60G images that exhibit exceptional imaging quality from the
entire dataset. This strategic selection ensures that our model is trained on
the most accurate and reliable data possible, setting the foundation for its
future performance.

##### CREMI.

The CREMI dataset is derived from the FAFB dataset, which has three manually
labeled subvolumes from drosophila brain, of which CREMI A has more regular
forms, CREMI C has higher size disparities, and CREMI B has in-between
segmentation difficulty. Each sub-volume has 125 slices of 1250$\times$1250
images, and we choose the first 60 slices for training, 15 slices for
validation, and the remaining 50 slices for testing, which are utilized to
validate our method’s performance on segmentation tasks of varying difficulty.

##### AC3/AC4.

AC3/AC4 Kasthuri et al. (2015) are mouse somatosensory cortex datasets with
256 and 100 successive EM images (1024$\times$1024), respectively. The first
80 slices of AC3 are used as the training set, the following 20 slices as the
validation set, and the first 50 slices of AC4 are used as the testing set.

##### Evaluation Metrics.

We are more interested in the model’s performance in the downstream task in
the self-supervised training. To assess the influence of segmentation on EM
neurons, we primarily use Variation of Information (VOI) Nunez-Iglesias et al.
(2013) and Adapted Rand Error (Arand) Arganda-Carreras et al. (2015) metrics.
Smaller VOI and ARAND values represent better segmentation results.

###  4.3 Experimental Results

##### Decision Making Process.

Decision-based MIM is pretrained using the FAFB dataset. During pretraining,
the policy network in the decision module is also updated, and the real-time
change in masking ratio is recorded as in Figure 4. The decision-making
process starts with random decisions, resulting in an overall masking ratio of
around 0.5. Then, using the reconstruction loss as a reward, the decision-
making principles of the agents are updated. After 50k iterations, the
decision module converges and produces better reconstruction results with less
lost information. During pretraining, it is found that starting with a lower
masking ratio is more beneficial for the EM dataset, and gradually increasing
the masking ratio helps the network’s learning process progress from easy to
difficult datasets, which is effective for both upstream reconstruction and
downstream segmentation tasks.

Figure 4: The variation of the masking ratio in the MARL decision process,
where the final convergence result shows that the optimal masking ratio
fluctuates around 0.83. Figure 5: The visualization results in the CREMI C
dataset reveal the areas of over-segmentation and under-segmentation produced
by methods such as Superhuman, MALA, and UNETR.

##### Results on CREMI.

We compare our method to the original MAE pretraining method He et al. (2022)
and the Dino pretraining method Caron et al. (2021) based on contrastive
learning, as well as two affinity-based Unet structures for instance
segmentation, superhuman Lee et al. (2017) and MALA Funke et al. (2018). The
outcomes of the trials are provided in Table 1 after applying both waterz and
LMC post-processing procedures. Figure 5 shows the visualization results.
Compared to existing self-supervised methods, our decision-based MIM approach
shows a significant improvement for the downstream task and it also performs
better than commonly used CNN-based methods in neuron segmentation, as can be
seen in Table 1.

##### Results on AC3/AC4.

In order to verify the robustness and generalizability of our pretraining
method, we further conduct experiments on mouse cortical neural cells. The
results in Table 1 indicate that our approach leads to a clear improvement in
the downstream segmentation task, even when applied across different species.

###  4.4 Ablation Study

We use UNETR to perform a comprehensive ablation study on the AC4 dataset.

##### Effectiveness of the multi-task reconstruction.

We conduct ablation experiments on the model’s reconstruction targets. As
shown in Table 2, it demonstrates that utilizing HOG and MSE losses together
for reconstruction in the MIM outperforms using MSE or HOG alone.

|  | waterz | LMC
---|---|---|---
MSE | HOG |  VOI $\downarrow$ |  ARAND $\downarrow$ |  VOI $\downarrow$ |  ARAND $\downarrow$
$\surd$ |  | 1.097 | 0.254 | 1.186 | 0.120
| $\surd$ | 0.987 | 0.251 | 1.103 | 0.127
$\surd$ | $\surd$ | 0.957 | 0.245 | 1.089 | 0.117
Table 2: Ablation results of the multi-task reconstruction.

##### Effectiveness of the decision module.

We compare our proposed method to a straightforward solution that manually
adjusts the masking ratio. The ablation results, as shown in Table 3,
demonstrate that our approach not only eliminates the need for manual
adjustment of masking ratios but also outperforms the best results achieved
through manual adjustment.

|  | waterz | LMC
---|---|---|---
Rate | Decision |  VOI $\downarrow$ |  ARAND $\downarrow$ |  VOI $\downarrow$ |  ARAND $\downarrow$
0.65 |  | 1.065 | 0.264 | 1.196 | 0.131
0.75 |  | 0.997 | 0.258 | 1.132 | 0.129
0.85 |  | 0.957 | 0.245 | 1.097 | 0.121
0.95 |  | 1.005 | 0.260 | 1.121 | 0.127
/ | $\surd$ | 0.931 | 0.243 | 1.079 | 0.113
Table 3: Ablation results of the decision module.

##  5 Conclusion

In this paper, we propose a decision-based MIM approach for neuron
segmentation. Our method eliminates the need for manual adjustment of masking
ratios and masking strategies, using multi-agent cooperation to search for the
optimal solution. Additionally, during the pretraining process, we incorporate
a multi-task reconstruction and utilize HOG features to enhance the model’s
learning ability. Our method is validated on a variety of EM neuron datasets
to demonstrate its generalizability.

## Acknowledgements

This work was supported in part by the National Natural Science Foundation of
China under Grant 62021001.

## References

  * Arganda-Carreras et al. [2015] Ignacio Arganda-Carreras, Srinivas C Turaga, Daniel R Berger, Dan Cireşan, Alessandro Giusti, Luca M Gambardella, Jürgen Schmidhuber, Dmitry Laptev, Sarvesh Dwivedi, Joachim M Buhmann, et al.  Crowdsourcing the creation of image segmentation algorithms for connectomics.  Frontiers in neuroanatomy, page 142, 2015.
  * Bachmann et al. [2022] Roman Bachmann, David Mizrahi, Andrei Atanov, and Amir Zamir.  MultiMAE: Multi-modal multi-task masked autoencoders.  In European Conference on Computer Vision, 2022.
  * Bandara et al. [2022] Wele Gedara Chaminda Bandara, Naman Patel, Ali Gholami, Mehdi Nikkhah, Motilal Agrawal, and Vishal M Patel.  Adamae: Adaptive masking for efficient spatiotemporal learning with masked autoencoders.  arXiv preprint arXiv:2211.09120, 2022.
  * Beier et al. [2017] Thorsten Beier, Constantin Pape, Nasim Rahaman, Timo Prange, Stuart Berg, Davi D Bock, Albert Cardona, Graham W Knott, Stephen M Plaza, Louis K Scheffer, et al.  Multicut brings automated neurite segmentation closer to human performance.  Nature methods, 14(2):101–102, 2017.
  * Briggman et al. [2009] Kevin Briggman, Winfried Denk, Sebastian Seung, Moritz Helmstaedter, and Srinivas C Turaga.  Maximin affinity learning of image segmentation.  Advances in Neural Information Processing Systems, 22, 2009.
  * Caron et al. [2021] Mathilde Caron, Hugo Touvron, Ishan Misra, Hervé Jégou, Julien Mairal, Piotr Bojanowski, and Armand Joulin.  Emerging properties in self-supervised vision transformers.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 9650–9660, 2021.
  * Chen et al. [2020] Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton.  A simple framework for contrastive learning of visual representations.  In International conference on machine learning, pages 1597–1607. PMLR, 2020.
  * Chen et al. [2023a] Yinda Chen, Wei Huang, Xiaoyu Liu, Qi Chen, and Zhiwei Xiong.  Learning multiscale consistency for self-supervised electron microscopy instance segmentation.  arXiv preprint arXiv:2308.09917, 2023.
  * Chen et al. [2023b] Yinda Chen, Che Liu, Wei Huang, Sibo Cheng, Rossella Arcucci, and Zhiwei Xiong.  Generative text-guided 3d vision-language pretraining for unified medical image segmentation.  arXiv preprint arXiv:2306.04811, 2023.
  * Dalal and Triggs [2005] Navneet Dalal and Bill Triggs.  Histograms of oriented gradients for human detection.  In 2005 IEEE computer society conference on computer vision and pattern recognition (CVPR’05), volume 1, pages 886–893. Ieee, 2005.
  * Deng et al. [2022] Shiyu Deng, Wei Huang, Chang Chen, Xueyang Fu, and Zhiwei Xiong.  A unified deep learning framework for sstem image restoration.  IEEE Transactions on Medical Imaging, 41(12):3734–3746, 2022.
  * Dosovitskiy et al. [2021] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  In International Conference on Learning Representations, 2021.
  * Feichtenhofer et al. [2022] Christoph Feichtenhofer, Yanghao Li, Kaiming He, et al.  Masked autoencoders as spatiotemporal learners.  Advances in neural information processing systems, 35:35946–35958, 2022.
  * Foerster et al. [2016] Jakob Foerster, Ioannis Alexandros Assael, Nando De Freitas, and Shimon Whiteson.  Learning to communicate with deep multi-agent reinforcement learning.  Advances in neural information processing systems, 29, 2016.
  * Funke et al. [2018] Jan Funke, Fabian Tschopp, William Grisaitis, Arlo Sheridan, Chandan Singh, Stephan Saalfeld, and Srinivas C Turaga.  Large scale image segmentation with structured loss based deep learning for connectome reconstruction.  IEEE transactions on pattern analysis and machine intelligence, 41(7):1669–1680, 2018.
  * Gao et al. [2022] Peng Gao, Teli Ma, Hongsheng Li, Ziyi Lin, Jifeng Dai, and Yu Qiao.  Mcmae: Masked convolution meets masked autoencoders.  In Advances in Neural Information Processing Systems, 2022.
  * Grill et al. [2020] Jean-Bastien Grill, Florian Strub, Florent Altché, Corentin Tallec, Pierre Richemond, Elena Buchatskaya, Carl Doersch, Bernardo Avila Pires, Zhaohan Guo, Mohammad Gheshlaghi Azar, et al.  Bootstrap your own latent-a new approach to self-supervised learning.  Advances in neural information processing systems, 33:21271–21284, 2020.
  * Hatamizadeh et al. [2022] Ali Hatamizadeh, Yucheng Tang, Vishwesh Nath, Dong Yang, Andriy Myronenko, Bennett Landman, Holger R Roth, and Daguang Xu.  Unetr: Transformers for 3d medical image segmentation.  In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pages 574–584, 2022.
  * He et al. [2022] Kaiming He, Xinlei Chen, Saining Xie, Yanghao Li, Piotr Dollár, and Ross Girshick.  Masked autoencoders are scalable vision learners.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 16000–16009, 2022.
  * Huang et al. [2022a] Wei Huang, Chang Chen, Zhiwei Xiong, Yueyi Zhang, Xuejin Chen, Xiaoyan Sun, and Feng Wu.  Semi-supervised neuron segmentation via reinforced consistency learning.  IEEE Transactions on Medical Imaging, 2022.
  * Huang et al. [2022b] Wei Huang, Shiyu Deng, Chang Chen, Xueyang Fu, and Zhiwei Xiong.  Learning to model pixel-embedded affinity for homogeneous instance segmentation.  In Proc. AAAI Conf. Artif. Intell, 2022.
  * Kasthuri et al. [2015] Narayanan Kasthuri, Kenneth Jeffrey Hayworth, Daniel Raimund Berger, Richard Lee Schalek, José Angel Conchello, Seymour Knowles-Barley, Dongil Lee, Amelio Vázquez-Reina, Verena Kaynig, Thouis Raymond Jones, et al.  Saturated reconstruction of a volume of neocortex.  Cell, 162(3):648–661, 2015.
  * Krasowski et al. [2017] NE Krasowski, Thorsten Beier, GW Knott, Ullrich Köthe, Fred A Hamprecht, and Anna Kreshuk.  Neuron segmentation with high-level biological priors.  IEEE transactions on medical imaging, 37(4):829–839, 2017.
  * Lee et al. [2017] Kisuk Lee, Jonathan Zung, Peter Li, Viren Jain, and H Sebastian Seung.  Superhuman accuracy on the snemi3d connectomics challenge.  arXiv preprint arXiv:1706.00120, 2017.
  * Liao et al. [2020] Xuan Liao, Wenhao Li, Qisen Xu, Xiangfeng Wang, Bo Jin, Xiaoyun Zhang, Yanfeng Wang, and Ya Zhang.  Iteratively-refined interactive 3d medical image segmentation with multi-agent reinforcement learning.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 9394–9402, 2020.
  * Lin et al. [2021a] Shiqi Lin, Tao Yu, Ruoyu Feng, Xin Li, Xin Jin, and Zhibo Chen.  Local patch autoaugment with multi-agent collaboration.  arXiv preprint arXiv:2103.11099, 2021.
  * Lin et al. [2021b] Shiqi Lin, Zhizheng Zhang, Xin Li, Wenjun Zeng, and Zhibo Chen.  Selectaugment: Hierarchical deterministic sample selection for data augmentation.  arXiv preprint arXiv:2112.02862, 2021.
  * Littman [1994] Michael L Littman.  Markov games as a framework for multi-agent reinforcement learning.  In Machine learning proceedings 1994, pages 157–163. Elsevier, 1994\.
  * Liu et al. [2022] Xiaoyu Liu, Wei Huang, Yueyi Zhang, and Zhiwei Xiong.  Biological instance segmentation with a superpixel-guided graph.  In IJCAI, 2022.
  * Mnih et al. [2016] Volodymyr Mnih, Adria Puigdomenech Badia, Mehdi Mirza, Alex Graves, Timothy Lillicrap, Tim Harley, David Silver, and Koray Kavukcuoglu.  Asynchronous methods for deep reinforcement learning.  In International conference on machine learning, pages 1928–1937. PMLR, 2016.
  * Nunez-Iglesias et al. [2013] Juan Nunez-Iglesias, Ryan Kennedy, Toufiq Parag, Jianbo Shi, and Dmitri B Chklovskii.  Machine learning of hierarchical clustering to segment 2d and 3d images.  PloS one, 8(8):e71715, 2013.
  * Ronneberger et al. [2015] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In International Conference on Medical image computing and computer-assisted intervention, pages 234–241. Springer, 2015.
  * Sheridan et al. [2022] Arlo Sheridan, Tri M Nguyen, Diptodip Deb, Wei-Chung Allen Lee, Stephan Saalfeld, Srinivas C Turaga, Uri Manor, and Jan Funke.  Local shape descriptors for neuron segmentation.  Nature Methods, pages 1–9, 2022.
  * Tang et al. [2022] Yucheng Tang, Dong Yang, Wenqi Li, Holger R Roth, Bennett Landman, Daguang Xu, Vishwesh Nath, and Ali Hatamizadeh.  Self-supervised pre-training of swin transformers for 3d medical image analysis.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 20730–20740, 2022.
  * Tong et al. [2022] Zhan Tong, Yibing Song, Jue Wang, and Limin Wang.  VideoMAE: Masked autoencoders are data-efficient learners for self-supervised video pre-training.  In Advances in Neural Information Processing Systems, 2022.
  * Turaga et al. [2010] Srinivas C Turaga, Joseph F Murray, Viren Jain, Fabian Roth, Moritz Helmstaedter, Kevin Briggman, Winfried Denk, and H Sebastian Seung.  Convolutional networks can learn to generate affinity graphs for image segmentation.  Neural computation, 22(2):511–538, 2010.
  * Wei et al. [2022] Chen Wei, Haoqi Fan, Saining Xie, Chao-Yuan Wu, Alan Yuille, and Christoph Feichtenhofer.  Masked feature prediction for self-supervised visual pre-training.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 14668–14678, 2022.
  * You et al. [2022] Chenyu You, Yuan Zhou, Ruihan Zhao, Lawrence Staib, and James S Duncan.  Simcvd: Simple contrastive voxel-wise representation distillation for semi-supervised medical image segmentation.  IEEE Transactions on Medical Imaging, 41(9):2228–2237, 2022.
  * Zhang et al. [2019] Xinyu Zhang, Qiang Wang, Jian Zhang, and Zhao Zhong.  Adversarial autoaugment.  arXiv preprint arXiv:1912.11188, 2019.
  * Zheng et al. [2018] Zhihao Zheng, J Scott Lauritzen, Eric Perlman, Camenzind G Robinson, Matthew Nichols, Daniel Milkie, Omar Torrens, John Price, Corey B Fisher, Nadiya Sharifi, et al.  A complete electron microscopy volume of the brain of adult drosophila melanogaster.  Cell, 174(3):730–743, 2018.
  * Zhou et al. [2019] Shenglong Zhou, Zhiwei Xiong, Chang Chen, Xuejin Chen, Dong Liu, Yueyi Zhang, Zheng-Jun Zha, and Feng Wu.  Fast and accurate electron microscopy image registration with 3d convolution.  In Medical Image Computing and Computer Assisted Intervention–MICCAI 2019: 22nd International Conference, Shenzhen, China, October 13–17, 2019, Proceedings, Part I 22, pages 478–486. Springer, 2019\.
  * Zhou et al. [2022] Lei Zhou, Huidong Liu, Joseph Bae, Junjun He, Dimitris Samaras, and Prateek Prasanna.  Self pre-training with masked autoencoders for medical image analysis.  arXiv preprint arXiv:2203.05573, 2022.


2023

[1]Evan H. Anders

1]Center for Interdisciplinary Exploration and Research in Astrophysics
(CIERA), Northwestern University, 1800 Sherman Ave, Evanston, 60201, IL, USA

2]Department of Engineering Sciences & Applied Mathematics, Northwestern
University, 2145 Sheridan Road, Evanston, 60208, IL, USA

3]Center for Computational Astrophysics, Flatiron Institute, 160 Fifth Avenue,
New York, 10010, New York, USA

4]Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane,
Princeton, 08544, NJ, USA

5]Department of Mathematics, Massachusetts Institute of Technology, 77
Massachusetts Avenue, Cambridge, 02139, MA, USA

6]Department of Physics & Astronomy, Northwestern University, 2145 Sheridan
Road, Evanston, 60208, IL, USA

7]Department of Astronomy, University of Wisconsin-Madison, 475 N. Charter
Street, Madison, 53706, WI, USA

8]Department of Astrophysical & Planetary Sciences, University of Colorado
Boulder, 391 UCB, Boulder, 80309, CO, USA

9]School of Mathematics, University of Edinburgh, James Clerk Maxwell
Building, Edinburgh, EH9 3FD, UK

10]Department of Physics & Astronomy, Bates College, 2 Andrews Road, Lewiston,
04240, ME, USA

#  The photometric variability of massive stars due to gravity waves excited
by core convection

<EMAIL_ADDRESS>Daniel Lecoanet  Matteo Cantiello  Keaton J.
Burns  Benjamin A. Hyatt  Emma Kaufman  Richard H. D. Townsend  Benjamin P.
Brown  Geoffrey M. Vasil  Jeffrey S. Oishi  Adam S. Jermyn  [  [  [  [  [  [
[  [  [  [

###### Abstract

Massive stars die in catastrophic explosions, which seed the interstellar
medium with heavy elements and produce neutron stars and black holes.
Predictions of the explosion’s character and the remnant mass depend on models
of the star’s evolutionary history. Models of massive star interiors can be
empirically constrained by asteroseismic observations of gravity wave
oscillations. Recent photometric observations reveal a ubiquitous red noise
signal on massive main sequence stars; a hypothesized source of this noise is
gravity waves driven by core convection. We present the first 3D simulations
of massive star convection extending from the star’s center to near its
surface, with realistic stellar luminosities. Using these simulations, we make
the first prediction of photometric variability due to convectively-driven
gravity waves at the surfaces of massive stars, and find that gravity waves
produce photometric variability of a lower amplitude and lower characteristic
frequency than the observed red noise. We infer that the photometric signal of
gravity waves excited by core convection is below the noise limit of current
observations, so the red noise must be generated by an alternative process.

###### pacs:

[

UAT Keywords]Stellar Photometry (1620), Internal waves (819), Stellar
Oscillations (1617), Hydrodynamical Simulations (767), Stellar Convection
Zones (301), Stellar Cores (1592)

##  1 Introduction

The oxygen we breathe was generated in the cores of massive stars
(rauscher_etal_2002, ) and expelled into the interstellar medium in violent
explosions (woosley_etal_2002, ) before mixing with the molecular cloud that
formed our solar system (kennicutt_evans_2012, ). In addition to producing the
elements that enable life, massive stars leave behind compact remnants, whose
subsequent mergers have provided a new window into the universe through
gravitational wave astronomy (abbott_etal_2017, ). Predictions of the
elemental yield of a massive star and the nature of its remnant are sensitive
to many processes including convection (kaiser_etal_2020, ), winds
(smith_2014, ), nuclear reaction rates (fields_etal_2018, ), and numerical
algorithms (farmer_etal_2019, ). Empirical constraints upon the interior
structures of massive stars (pedersen_etal_2021, ) could reduce uncertainties
in these processes, which would improve stellar evolution calculations, and in
turn models of star formation (mckee_ostriker_2007, ), galaxy formation
(pillepich_etal_2018, ), and the reionization of the early universe
(bromm_larson_2004, ).

It was recently shown that photometric light curves of hot, massive stars
contain a ubiquitous red noise signal (bowman_etal_2019, ;
bowman_etal_2019_nat, ; bowman_etal_2020, ; dornwallenstein_etal_2020, ;
szewczuk_etal_2021, ; bowman_dornwallenstein_2022, ). Theories for the driving
mechanism of red noise include gravity waves stochastically excited by core
convection, or turbulence from subsurface convection zones
(bowman_etal_2019_nat, ; cantiello_etal_2021, ; schultz_etal_2022, ). Of
particular interest is the hypothesis that these signals are waves driven by
core convection, which could empirically constrain interior models of massive
stars.

We present here the first simulations of core convection that include all the
relevant physics to make an accurate prediction of the stellar photometric
variability from convectively-generated waves. Our simulations build upon
pioneering studies which examined the properties of the turbulent core
convection (e.g., browning_etal_2004, ; gilet_etal_2013, ;
augustson_etal_2016, ; higl_etal_2021, ; baraffe_etal_2023, ), the
characteristics of the gravity waves in the radiative envelope
(rogers_etal_2013, ; rogers_mcelwaine_2017, ; couston_etal_2018, ;
edelmann_etal_2019, ; horst_etal_2020, ; lecoanet_etal_2021, ;
lesaux_etal_2022, ; blouin_etal_2023, ), and the possible observational
features of those waves (lecoanet_etal_2021, ; samadi_etal_2010, ;
shiode_etal_2013, ; aerts_rogers_2015, ; lecoanet_etal_2019, ;
thompson_etal_2023, ).

Previous studies focused on the shape of the frequency spectrum of gravity
waves from convection (aerts_rogers_2015, ; edelmann_etal_2019, ;
horst_etal_2020, ; thompson_etal_2023, ), but could not constrain the
amplitude of the waves because their simulations do not extend to the surface
and use boosted convective luminosities. Here we present simulations with
realistic stellar luminosities which we use to predict both the amplitude and
shape of the frequency spectrum of photometric variability due to gravity
waves excited by core convection.

Predicting this frequency spectrum is challenging because of the range of
length scales, time scales, and physical processes important in the star
(jermyn_etal_2022_atlas, ). Waves are generated on fast convective timescales
($\sim 14$ days), but as they propagate to the surface they excite resonances
with much longer lifetimes ($\gg 4\,{\rm yrs}$ from Kepler data
(kurtz_etal_2015, ); $\sim 10^{5}\,{\rm yrs}$ from theoretical estimates
(shiode_etal_2013, )). It is unfeasible to simultaneously resolve these
timescales. To determine the wave signal at the surface, we appeal to an
acoustic analogy.

The character of music depends both on the sound waves produced by musicians
and on the acoustics of the environment where it is played (e.g.,
pentcheva_abel_2017, ). Music is recorded in special studios with walls that
absorb or diffuse waves to minimize the influence of the environment on the
sound and retrieve the “pure sound” of the musicians. To experience music in a
different environment, it is not necessary to physically transport the
musicians; instead, one can apply a filter to the recording, mimicking the
effects of the new environment.

Figure 1:  Snapshots of equatorial slices through 3D spherical simulations of
our fiducial 15 $M_{\odot}$ star. (Left Panel) Temperature fluctuations around
the mean profile are shown in our “Wave Propagation” simulation; this
simulation spans 93% $R_{*}$, and the hatched region displays the $0.07R_{*}$
that we do not simulate. The gravity waves cause perturbations of $\pm 1$ K,
while the core convective temperature perturbations are of order $\pm 30$ K.
(Middle Panel) The radial velocity field in our “Wave Generation” simulation
is shown, normalized by the standard deviation of the radial velocity at each
radius. The Wave Generation simulation radially includes 55% of the 15
$M_{\odot}$ star, and the outer portion of the simulation contains a wave
damping layer (indicated by the grey shading). (Right Panel) A zoomed in view
of the radial velocity in the Wave Generation simulation at the same instant.
The waves are not visible because their velocities are much smaller than the
convective velocities. An animated version of this figure is available in the
supplemental materials.

Our strategy for determining the photometric variability from gravity waves is
analogous. We run short simulations of wave generation by convection and
“record” the waves as they leave the convection zone. We also run a fully
self-consistent simulation of convection and wave propagation in a truncated
stellar model. To test our strategy, we apply to our wave recording a filter
describing the effects of wave propagation in the truncated model. The surface
variability of the self-consistent simulation matches the filtered wave
signal. Having verified our strategy, we then apply a different filter
associated with the real star to determine the stellar photometric
variability.

##  2 Results

To predict the photometric variability caused by gravity waves in massive
stars, we used Dedalus to run the first 3D simulations of core convection that
simultaneously solve the fully compressible equations, produce realistically
low-Mach-number flows, and include the full radial extent of the core
convection zone, including $r=0$. We study two types of simulations: (1) “Wave
Generation” simulations where we record the waves (couston_etal_2018, ), and a
“Wave Propagation” simulation where we test our method of applying a filter to
the “wave recording” to mimic the effects of the environment
(lecoanet_etal_2021, ). The Wave Generation simulations span twice the radial
extent of the core convection zone (55% of the star by radius in our fiducial
model) and include a damping region in the outer $\sim$10% of the simulation
domain by radius to prevent wave reflections. The Wave Propagation simulation
spans 93% of the star by radius (99.99925% by mass) and is time-integrated for
hundreds of convective overturn timescales to allow power in standing wave
modes to saturate. While both types of simulations include wave generation and
propagation, these simulation designs allow us to isolate and study (1) the
wave luminosity generated by the convective core, and (2) the surface
amplitude of gravity waves that have propagated through the radiative
envelope. Most of our simulations are based on a 15 $M_{\odot}$ MESA stellar
model with LMC metallicity at the zero-age main sequence (ZAMS); this model is
representative of a star with no near-surface convection
(jermyn_etal_2022_window, ).

In Fig. 1, we display temperature fluctuations for the Wave Propagation
simulation (left panel) and radial velocity for a Wave Generation simulation
(middle & right panels). Wave excitation can be visibly seen as an overdensity
of gravity waves near the core convection zone, e.g., in the bottom right of
the middle panel. In the right panel, we zoom in on the core convection zone;
the RMS velocity of the core convection is $5.8\times 10^{4}$ cm/s, similar to
the MLT prediction of $6.7\times 10^{4}$ cm/s from the MESA model. The
convection in our simulation is turbulent with an average Reynolds number of
$9.5\times 10^{3}$ and Mach number of $9.5\times 10^{-4}$. These Wave
Propagation and Wave Generation simulations use spherical coordinates with
respectively 512 and 1024 resolution elements across the convective core (see
Supplemental Materials Sec. 2.4).

Figure 2:  (Left panel) We show the predicted amplitude spectrum of gravity
waves from our fiducial 15 $M_{\odot}$ simulation (purple line) and Lorentzian
fits to red noise observations of stars with similar spectroscopic luminosity
$\mathscr{L}$ and effective temperature $T_{\rm eff}$ (bowman_etal_2020, )
(pink lines; the colors of the lines correspond to the colors of the points in
the right panel). (Middle panel) A comparison of predicted amplitude spectra
of gravity waves for the three stars that we simulate in this work. Note the
different y-axis scaling between the left and middle panels. (Right panel) A
spectroscopic HR diagram of the upper main sequence. The grey lines in the
background show the main sequence evolution (from lower left to upper right)
of nonrotating solar metallicity stars of [3, 5, 10, 15, 20, 40] $M_{\odot}$
from the MIST grids (mist0, ; mist1, ). The observed stars for which red noise
fits are plotted in the left panel are shown as pink circles (darkened in
color with distance from the MIST ZAMS in the
$\log_{10}\mathscr{L}-\log_{10}T_{\rm eff}$ plane). The three simulated stars
from the middle panel are displayed as star symbols.

The purpose of the Wave Propagation simulation is to test our method for
applying a filter to the waves to mimic how the stellar environment modifies
the waves (lecoanet_etal_2019, ; lecoanet_etal_2021, ). To determine the wave
amplitude in a star, we need (1) a recording of the waves generated by core
convection and (2) a transfer function (“filter”) describing how waves
propagate through the radiative envelope. The wave luminosity spectrum is
excited by stochastic, turbulent fluctuations in the convective core. We
measure (“record”) the wave luminosity of the traveling waves in Wave
Generation simulations that have different resolutions and turbulent
intensities. We find the wave luminosity is similar to theoretical predictions
(lecoanet_quataert_2013, ) and past simulations (couston_etal_2018, ;
lecoanet_etal_2021, ; lesaux_etal_2022, ), and is independent of resolution
provided there are at least 512 resolution elements across the convection zone
(Supplemental Materials Sec. 3.2). We separately evolve the Wave Propagation
simulation for a very long time, and in doing so self-consistently evolve the
generation and propagation of gravity waves in an environment similar to the
star. We build a transfer function (“filter”) using the gravity wave
eigenvalues and eigenfunctions associated with the wave cavity of the Wave
Propagation simulation (Supplemental Materials Sec. 5.2). We synthesize a wave
signal by convolving the wave luminosity from the Wave Generation simulation
with this filter, and find good agreement between wave perturbations measured
at the surface of the Wave Propagation simulation and our synthesized signal,
verifying our method, see Supplemental Materials Sec. 5.3.

The Wave Propagation simulation is not identical to the star, because it has
different diffusivities and a different wave cavity due to only including 93%
of the star, see Supplemental Materials Sec. 2.6. Our Wave Propagation
simulation has mode lifetimes of $\lesssim 10$ years, allowing us to capture
both convection and wave propagation, which is not possible for real stars
where mode lifetimes are $\sim 10^{5}$ years. Having verified our method using
the Wave Propagation simulation, we then create a separate “full” transfer
function based on the true MESA stellar stratification and the gravity waves
associated with the star’s wave cavity, which we calculate using GYRE , see
Supplemental Materials Sec. 6. We pass the wave luminosity signal from the
Wave Generation simulation through this “full” transfer function to predict
the photometric variability due to gravity waves.

Figure 3:  Shown are the magnitude $\alpha_{0}$ (left panel) and
characteristic frequency $\nu_{\rm char}$ (right panel) obtained from fitting
Eqn. 1 to our simulated gravity wave spectra (star symbols). We also include
best-fits of red noise observations (circles) (bowman_etal_2020, ). All points
are colored by $\log_{10}T_{\rm eff}$ and plotted against
$\log_{10}\mathscr{L}/\mathscr{L}_{\odot}$.

In the left panel of Fig. 2, we plot the predicted signal of gravity waves
generated by core convection in our 15 $M_{\odot}$ star. We make predictions
for both the shape of the wave signal and its amplitude. We find a wave signal
characterized by a broad peak at a frequency near 0.1 d-1. The low frequency
side of the peak is smooth and drops sharply, while the high frequency side of
the peak decreases gradually and includes high amplitude, narrow peaks
associated with resonant gravity modes. The peak amplitude of the signal is
roughly ${0.06\,\mu\text{mag}}$. We also display the publicly-available red
noise observation fits (bowman_etal_2020, ), which are flat at low frequency,
decrease sharply above $\sim$1 d-1, and have amplitudes $\gtrsim 10$ $\mu$mag.
We note that these fits underwent a pre-whitening procedure which iteratively
removed the standing modes; the raw observations do include some peaks from
standing modes. We include these fits primarily to provide a direct comparison
between the typical _amplitude_ of red noise observations and our simulation-
based predictions.

In the middle panel of Fig. 2, we compare the simulated wave signals of three
ZAMS stars of different masses. These signals are similar, characterized by a
broad peak at low frequency that has a sharp decrease on its low-frequency
side, due to radiative damping, and a gradual decrease at high-frequency, due
to less effective convective wave excitation and cancellation effects for high
spherical harmonic degree (see Supplemental Materials Sec. 6.2). The high
frequency portions of the spectra are dominated by resonant wave peaks, and
lower mass stars have more peaks and peaks of higher quality factor (see
Supplemental Materials Figs. 8 & 9). We find that the overall signal shifts to
higher amplitude and lower frequency as the stellar mass increases. The stars
considered in this work are displayed in a spectroscopic HR diagram in the
right panel of Fig. 2. Our fiducial 15 $M_{\odot}$ LMC star is shown as a
purple star, and we also simulate solar metallicity ZAMS stars of $3M_{\odot}$
(green star) and $40M_{\odot}$ (orange star). The stars whose red noise
signals we plotted are shown as pink circles, and they have
$\log_{10}\mathscr{L}$ and $\log_{10}T_{\rm eff}$ which respectively differ by
no more than $0.2$ from our fiducial star.

We fit our simulation spectra with Lorentzian profiles (bowman_etal_2020, )
(Supplemental Materials Sec. 6.3),

$\alpha(\nu)=\frac{\alpha_{0}}{1+\left(\frac{\nu}{\nu_{\rm
char}}\right)^{\gamma}}.$ (1)

In Fig. 3 we plot the amplitude $\alpha_{0}$ and characteristic frequency
$\nu_{\rm char}$ of both the gravity wave signals from our simulations and
observed photometric variability (red noise). The amplitude $\alpha_{0}$ of
the gravity wave signal increases in magnitude with increasing stellar
luminosity and mass, similar to the amplitude of red noise. The characteristic
frequency $\nu_{\rm char}$ of the gravity wave signal decreases with
increasing stellar luminosity and mass, whereas the characteristic frequency
of red noise increases with increasing effective temperature but otherwise
remains roughly constant as luminosity changes.

##  3 Discussion

Our models predict surface signals of internal gravity waves generated by
massive star core convection that are inconsistent with red noise
observations. The photometric variability due to gravity waves is orders of
magnitude lower than the observed red noise. Our results suggest that a
mechanism other than core-driven gravity waves must be responsible for the
ubiquitous red noise signal. Our results are corroborated by recent 2D
simulations with realistic luminosities lesaux_2023 , which found a wave
luminosity spectrum with the same shape as ours, and which also found that
wave propagation through the stellar envelope can be correctly modeled using
linear theory.

While our simulations are the most realistic simulations of gravity wave
generation by core convection to date, they do not include the following
physical effects: magnetic fields, stellar models beyond the ZAMS, near
surface convection zones, and rotation.

  * •

While the effects of magnetism are unclear, there are no currently known
theoretical mechanisms by which they would substantially enhance wave
generation by core convection.

  * •

As stars age away from the ZAMS, they develop composition $\mu$-gradients
outside of the core. These $\mu$-gradients modify the wave cavity and would
affect the precise shape of the transfer function that we use to mimic the
waves in the environment of a full star. We do not expect composition
gradients to have a large effect on the generation of waves. Therefore, while
compositional gradients could affect the shape of the observed spectra, there
are no currently known theoretical mechanisms by which they would increase the
amplitude of surface fluctuations.

  * •

Our fiducial 15 $M_{\odot}$ LMC model is not expected to have near-surface
convection, so our photometric variability predictions for this star could be
directly compared to LMC data in this regime where near-surface convective
zones are expected to be absent (jermyn_etal_2022_window, ). However, the
observed galactic stars that we compare our results to are all expected to
have near-surface convection, which could pollute or modify the wave signals
we predict here. The 3 $M_{\odot}$ and 40 $M_{\odot}$ Wave Generation
simulations we conducted are based on MESA models with solar metallicities and
are expected to have near surface convection zones. Recent analyses of near-
surface convection (schultz_etal_2022, ) suggest that turbulence driven by
subsurface convection could manifest as red noise, and future work should
further examine the nature of turbulence generated by near-surface convection
zones.

  * •

Rotation affects the generation of waves, and could enhance the surface
perturbations (augustson_etal_2020, ). We performed one 15 $M_{\odot}$ Wave
Generation simulation with a moderate rotation period $P_{\rm rot}=10$ d
(Supplemental Materials Sec. 7.2). We find that this rotation boosts the
photometric signal amplitude to $\alpha_{0}=0.21$ $\mu$mag, still orders of
magnitude weaker than the observed red noise. While this work cannot rule out
the possibility that convection excites gravity waves to observable amplitudes
in more rapidly rotating stars, the red noise signal is ubiquitous in both
slow and rapid rotators, and there is no correlation between the stellar
rotation rate and the red noise amplitude (see Supplementary Materials Sec.
7). Thus, the ubiquitous red noise signal cannot be due to gravity waves
excited by core convection. While rotational wave splitting would further
modify the shape of the observed spectra, it would lower the peak wave
amplitude compared to our non-rotating predictions.

While we have shown the photometric variability from convectively-excited
gravity waves is not directly detectable, these waves may mix chemicals
(varghese_etal_2023, ) or transport angular momentum (rogers_etal_2013, ),
leading to observable signals in more traditional $g$-mode asteroseismology
(decat_aerts_2002, ; pedersen_etal_2021, ; pedersen_2022, ). Even though the
red noise signal is not the surface manifestation of gravity waves, it still
carries valuable information about both the near-surface structure of massive
stars as well as their masses and ages (bowman_etal_2020, ;
bowman_dornwallenstein_2022, ).

Methods

All of the calculations in this work are based upon stellar models that we
computed using Modules for Experiments in Stellar Astrophysics (MESA), version
r21.12.1 (mesa1, ; mesa2, ; mesa3, ; mesa4, ; mesa5, ; mesa6, ), see
Supplemental Materials Sec. 1.

3D numerical simulations are performed using Dedalus (burns_etal_2020, )
version 3; the specific version of the code run was obtained from the master
branch of the Dedalus Github repository
(https://github.com/dedalusProject/dedalus) at the commit with short-sha
29f3a59. We simulate the fully-compressible equations and assume the fluid is
composed of a calorically perfect, uniform composition ideal gas. We use a
grid with high radial resolution from $r=0$ to $r=1.1R_{\rm core}$, where
$R_{\rm core}$ is the radius of the core convection zone, and we use a lower
radial resolution grid spanning the radiative envelope (See Supplemental
Materials table 3). The angular variation of all simulation variables are
expanded using a basis of spin-weighted spherical harmonics; the radial
variation of all simulations variables are expanded using a basis of radially-
weighted Zernike polynomials in the convective core (vasil_etal_2019_tcalc, ;
lecoanet_etal_2019_tcalc, ), and Chebyshev polynomials in the radiative
envelope. See Supplemental Materials, Sec. 2 for full simulation details.

Most of our dynamical simulations are “Wave Generation” simulations, which
extend from $r=0$ to $2R_{\rm core}$. These simulations include a damping
layer near the outer boundary of the simulation. In these simulations, we
measure the enthalpy and radial velocity perturbations on spherical shells at
various radii throughout the simulation radiative zone at 30 minute intervals
and use those measurements to calculate the wave luminosity as function of
frequency and spherical harmonic degree. Wave Generation simulations are time-
evolved for roughly 50 convective overturn timescales to build up sufficient
statistics to measure the wave luminosity. We find that the wave luminosity is
insensitive to the resolution and Reynolds number of calculations. See
Supplemental Materials, Sec. 3 for details on our measurements of the wave
luminosity in Wave Generation simulations.

In order to connect the convective driving to observable magnitude
fluctuations at the stellar surface, we require a transfer function. We derive
a transfer function closely following the procedure laid out in appendices A
and B of Ref. (lecoanet_etal_2019, ) with a few modifications, see Sec. 4 of
the Supplemental Materials. We note that the transfer function includes an
$\mathcal{O}(1)$ amplitude correction factor $A_{\rm corr}$ which must be
calibrated using 3D numerical simulations before it can be used to precisely
predict the magnitude fluctuations at the surface of a star. In this work, we
use $A_{\rm corr}=0.4$.

We run one “Wave Propagation” simulation, pictured in the left panel of Fig.
1. This simulation includes 93% by radius of the fiducial 15 $M_{\odot}$
stellar model, does not include a damping layer near the outer boundary, and
is time-evolved for hundreds of convective overturn timescales in order to
allow standing mode amplitudes to saturate. We evolve this simulation for 5
years, which is longer than the mode lifetime of almost all waves (the longest
mode lifetime is $\sim$10 years, see Supplemental Materials Fig. 5). We
measure the entropy perturbations at the outer boundary of this Wave
Propagation simulation at 30 min intervals to determine the spectrum of
gravity waves in a self-consistent, nonlinear simulation which includes most
of the star. We verify that our prediction built from the transfer function
and wave luminosity reproduces the measured surface perturbations well. We
furthermore ran one more Wave Propagation simulation at a lower resolution
(256 resolution elements across the convective core) and found that the same
amplitude correction factor of $A_{\rm corr}=0.4$ was required to align the
transfer function prediction and surface perturbations, suggesting this factor
is independent of the resolution, and thus applicable to stellar parameters.
See Supplemental Materials, Sec. 5, for details on the surface spectrum of the
Wave Propagation simulation, the verification of the transfer function, and
the calibration of the amplitude correction factor.

We predict the magnitude perturbation at the surface of a star from the
nonadiabatic gravity wave eigenvalues and eigenvectors calculated using GYRE
version 7.0 (townsend_teitler_2013, ; townsend_goldstein_zweibel_2018, ). We
retrieve the GYRE g-modes with the standard VACUUM boundary conditions, which
impose zero pressure perturbation at the stellar surface; we tested other
boundary conditions formulations (DZIEM, UNNO, JCD) and found that our results
were insensitive to this choice. We compute the limb-darkened, disk-averaged
differential flux functions from Eqns. 12-14 of Ref. (townsend2003, ),
accounting for an arbitrary observing angle as in Eqn. 8 of Ref.
(townsend2002, ), using MSG, version 1.1.2 (townsend_lopez_2023, ), which
synthesizes stellar spectra and convolves those spectra with an instrumental
passband to account not just for intrinsic luminosity changes but also changes
in the observed magnitude that a chosen telescope would see. We use the
OSTAR2002 grid in our $15M_{\odot}$ and $40M_{\odot}$ calculations
(lanz_hubeny_2003, ) and an extended version of the MSG demo grid for the
$3M_{\odot}$ star (castelli_kurucz_2003, ) and we select for our instrumental
passband the TESS passband from the SVO Filter Profile Service
(rodrigo_solano_2020, ). After synthesizing the differential flux functions,
we calculate the magnitude perturbation eigenfunction associated with each
eigenvalue using the GYRE solutions and Eqn. 11 of Ref. (townsend2002, ). We
then use those magnitude eigenfunctions to generate a transfer function, which
we use alongside the wave luminosity from Wave Generation simulations to
predict the photometric variability of gravity waves. The stellar transfer
functions, synthesized spectra, and red noise fits are discussed and shown in
Supplemental Materials, Sec. 6.

We additionally run one Wave Generation simulation with a rotation period of
$P_{\rm rot}=10$ d; we assume uniform, rigid rotation and include rotational
effects by adding the Coriolis term to the momentum equation. We measure the
wave luminosity in the same way as in the non-rotating simulation, and for
simplicity we use the non-rotating transfer function. See Supplemental
Materials, Sec. 7, for details of this simulation and a discussion of the
importance of rotation.

All plots were produced using matplotlib v3.5.2 (Hunter:2007, ; mpl3.5.2, ),
and the numpy v1.22.4 (harris2020array, ), scipy v1.8.1 (virtanen_etal_2020,
), and astropy v5.1 (astropy:2013, ; astropy:2018, ; astropy:2022, ) packages
were used in our simulation pre-processing and post-processing. MESA profiles
I/O was handled using https://github.com/wmwolf/py_mesa_reader.

Supplementary information

We refer the reader to the supplementary materials, which contain a full
description of the Dedalus simulation implementation, usage details of MESA,
details of the transfer function derivation, and figures showing wave
luminosity spectra, wave power spectra at the surface of the Wave Propagation
simulation, transfer functions, and rotating simulation results.

Code Availability

All Python scripts, MESA inlists, and GYRE inlists used in this work are in
the git repository located at
https://github.com/evanhanders/gmode_variability_paper. This repository is
additionally backed up in a zenodo repository (zenodo, ). Most of the
underlying logic that implements the equations in Dedalus and generates
transfer functions from Dedalus or GYRE outputs are separately located in the
pip-installable github repository
https://github.com/evanhanders/compressible_stars, version 0.1.

Data Availability

The data used to create all of the figures in this paper, as well as the
scripts used to generate the simulations and figures, are available online in
this zenodo repository: (zenodo, ).

Corresponding Author Please address correspondence and requests for materials
to Evan H. Anders.

Acknowledgments

We thank Conny Aerts, Mathias Michielsen, Timothy Van Reeth, Cole Johnston,
Kyle Augustson, Meng Sun, Will Schultz, Lars Bildsten, Dominic Bowman, Falk
Herwig, and Paul Woodward for helpful discussions which improved the work in
this manuscript. We are grateful to the two anonymous referees for taking the
time to review our work and for their valuable feedback which improved the
clarity of this manuscript and the supplemental materials. Computations were
conducted with support by the NASA High End Computing Program through the NASA
Advanced Supercomputing (NAS) Division at Ames Research Center on Pleiades
with allocation GIDs s2276. This work improved thanks to fruitful discussions
at the “Probes of Transport in Stars” KITP program in Fall 2021.

Funding

EHA was supported by a CIERA Postdoctoral Fellowship. This work was supported
in part by NASA HTMS grant 80NSSC20K1280, NASA SSW grant 80NSSC19K0026, and
NASA OSTFL grant 80NSSC22K1738. RHDT acknowledges support from NSF grants
ACI-1663696, AST-1716436 and PHY-1748958, and NASA grant 80NSSC20K0515. The
Center for Computational Astrophysics at the Flatiron Institute is supported
by the Simons Foundation. This research was also supported in part by the
National Science Foundation (NSF) under Grant Number NSF PHY-1748958.

Individual Contributions EHA, DL, and MC conceived and designed the
experiments. EHA performed the experiments. EHA, DL, MC, EK, and BPB analyzed
the data. EHA, DL, KJB, BAH, RHDT, GMV, and ASJ contributed materials/analysis
tools. EHA, DL, MC, BAH, and BPB wrote the paper and supplemental materials.

## References

  *   * (1) Rauscher, T., Heger, A., Hoffman, R.D., Woosley, S.E.: Nucleosynthesis in Massive Stars with Improved Nuclear and Stellar Physics. ApJ 576(1), 323–348 (2002) arXiv:astro-ph/0112478 [astro-ph]. https://doi.org/10.1086/341728
  * (2) Woosley, S.E., Heger, A., Weaver, T.A.: The evolution and explosion of massive stars. Reviews of Modern Physics 74(4), 1015–1071 (2002). https://doi.org/10.1103/RevModPhys.74.1015
  * (3) Kennicutt, R.C., Evans, N.J.: Star Formation in the Milky Way and Nearby Galaxies. ARA&A 50, 531–608 (2012) arXiv:1204.3552 [astro-ph.GA]. https://doi.org/10.1146/annurev-astro-081811-125610
  * (4) Abbott, B.P., et al., (LIGO Scientific Collaboration, Virgo Collaboration): GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral. Phys. Rev. Lett. 119(16), 161101 (2017) arXiv:1710.05832 [gr-qc]. https://doi.org/10.1103/PhysRevLett.119.161101
  * (5) Kaiser, E.A., Hirschi, R., Arnett, W.D., Georgy, C., Scott, L.J.A., Cristini, A.: Relative importance of convective uncertainties in massive stars. MNRAS 496(2), 1967–1989 (2020) arXiv:2006.01877 [astro-ph.SR]. https://doi.org/10.1093/mnras/staa1595
  * (6) Smith, N.: Mass Loss: Its Effect on the Evolution and Fate of High-Mass Stars. ARA&A 52, 487–528 (2014) arXiv:1402.1237 [astro-ph.SR]. https://doi.org/10.1146/annurev-astro-081913-040025
  * (7) Fields, C.E., Timmes, F.X., Farmer, R., Petermann, I., Wolf, W.M., Couch, S.M.: The Impact of Nuclear Reaction Rate Uncertainties on the Evolution of Core-collapse Supernova Progenitors. ApJS 234(2), 19 (2018) arXiv:1712.06057 [astro-ph.SR]. https://doi.org/10.3847/1538-4365/aaa29b
  * (8) Farmer, R., Renzo, M., de Mink, S.E., Marchant, P., Justham, S.: Mind the Gap: The Location of the Lower Edge of the Pair-instability Supernova Black Hole Mass Gap. ApJ 887(1), 53 (2019) arXiv:1910.12874 [astro-ph.SR]. https://doi.org/10.3847/1538-4357/ab518b
  * (9) Pedersen, M.G., Aerts, C., Pápics, P.I., Michielsen, M., Gebruers, S., Rogers, T.M., Molenberghs, G., Burssens, S., Garcia, S., Bowman, D.M.: Internal mixing of rotating stars inferred from dipole gravity modes. Nature Astronomy 5, 715–722 (2021) arXiv:2105.04533 [astro-ph.SR]. https://doi.org/10.1038/s41550-021-01351-x
  * (10) McKee, C.F., Ostriker, E.C.: Theory of Star Formation. ARA&A 45(1), 565–687 (2007) arXiv:0707.3514 [astro-ph]. https://doi.org/10.1146/annurev.astro.45.051806.110602
  * (11) Pillepich, A., Springel, V., Nelson, D., Genel, S., Naiman, J., Pakmor, R., Hernquist, L., Torrey, P., Vogelsberger, M., Weinberger, R., Marinacci, F.: Simulating galaxy formation with the IllustrisTNG model. MNRAS 473(3), 4077–4106 (2018) arXiv:1703.02970 [astro-ph.GA]. https://doi.org/10.1093/mnras/stx2656
  * (12) Bromm, V., Larson, R.B.: The first stars. Annual Review of Astronomy and Astrophysics 42(1), 79–118 (2004) https://doi.org/10.1146/annurev.astro.42.053102.134034. https://doi.org/10.1146/annurev.astro.42.053102.134034
  * (13) Bowman, D.M., Aerts, C., Johnston, C., Pedersen, M.G., Rogers, T.M., Edelmann, P.V.F., Simón-Díaz, S., Van Reeth, T., Buysschaert, B., Tkachenko, A., Triana, S.A.: Photometric detection of internal gravity waves in upper main-sequence stars. I. Methodology and application to CoRoT targets. A&A 621, 135 (2019) arXiv:1811.08023 [astro-ph.SR]. https://doi.org/10.1051/0004-6361/201833662
  * (14) Bowman, D.M., Burssens, S., Pedersen, M.G., Johnston, C., Aerts, C., Buysschaert, B., Michielsen, M., Tkachenko, A., Rogers, T.M., Edelmann, P.V.F., Ratnasingam, R.P., Simón-Díaz, S., Castro, N., Moravveji, E., Pope, B.J.S., White, T.R., De Cat, P.: Low-frequency gravity waves in blue supergiants revealed by high-precision space photometry. Nature Astronomy 3, 760–765 (2019) arXiv:1905.02120 [astro-ph.SR]. https://doi.org/10.1038/s41550-019-0768-1
  * (15) Bowman, D.M., Burssens, S., Simón-Díaz, S., Edelmann, P.V.F., Rogers, T.M., Horst, L., Röpke, F.K., Aerts, C.: Photometric detection of internal gravity waves in upper main-sequence stars. II. Combined TESS photometry and high-resolution spectroscopy. A&A 640, 36 (2020) arXiv:2006.03012 [astro-ph.SR]. https://doi.org/%****␣article.bbl␣Line␣350␣****10.1051/0004-6361/202038224
  * (16) Dorn-Wallenstein, T.Z., Levesque, E.M., Neugent, K.F., Davenport, J.R.A., Morris, B.M., Gootkin, K.: Short-term Variability of Evolved Massive Stars with TESS. II. A New Class of Cool, Pulsating Supergiants. ApJ 902(1), 24 (2020) arXiv:2008.11723 [astro-ph.SR]. https://doi.org/10.3847/1538-4357/abb318
  * (17) Szewczuk, W., Walczak, P., Daszyńska-Daszkiewicz, J.: Variability of newly identified B-type stars observed by Kepler. MNRAS 503(4), 5894–5928 (2021) arXiv:2103.06146 [astro-ph.SR]. https://doi.org/10.1093/mnras/stab683
  * (18) Bowman, D.M., Dorn-Wallenstein, T.Z.: Photometric detection of internal gravity waves in upper main-sequence stars. III. Comparison of amplitude spectrum fitting and Gaussian process regression using CELERITE2. A&A 668, 134 (2022) arXiv:2211.08347 [astro-ph.SR]. https://doi.org/10.1051/0004-6361/202243545
  * (19) Cantiello, M., Lecoanet, D., Jermyn, A.S., Grassitelli, L.: On the Origin of Stochastic, Low-Frequency Photometric Variability in Massive Stars. ApJ 915(2), 112 (2021) arXiv:2102.05670 [astro-ph.SR]. https://doi.org/%****␣article.bbl␣Line␣425␣****10.3847/1538-4357/ac03b0
  * (20) Schultz, W.C., Bildsten, L., Jiang, Y.-F.: Stochastic Low-frequency Variability in Three-dimensional Radiation Hydrodynamical Models of Massive Star Envelopes. ApJ 924(1), 11 (2022) arXiv:2110.13944 [astro-ph.SR]. https://doi.org/10.3847/2041-8213/ac441f
  * (21) Browning, M.K., Brun, A.S., Toomre, J.: Simulations of Core Convection in Rotating A-Type Stars: Differential Rotation and Overshooting. ApJ 601(1), 512–529 (2004) arXiv:astro-ph/0310003 [astro-ph]. https://doi.org/10.1086/380198
  * (22) Gilet, C., Almgren, A.S., Bell, J.B., Nonaka, A., Woosley, S.E., Zingale, M.: Low Mach Number Modeling of Core Convection in Massive Stars. ApJ 773(2), 137 (2013). https://doi.org/10.1088/0004-637X/773/2/137
  * (23) Augustson, K.C., Brun, A.S., Toomre, J.: The Magnetic Furnace: Intense Core Dynamos in B Stars. ApJ 829(2), 92 (2016) arXiv:1603.03659 [astro-ph.SR]. https://doi.org/10.3847/0004-637X/829/2/92
  * (24) Higl, J., Müller, E., Weiss, A.: Calibrating core overshooting parameters with two-dimensional hydrodynamical simulations. A&A 646, 133 (2021) arXiv:2012.05262 [astro-ph.SR]. https://doi.org/10.1051/0004-6361/202039532
  * (25) Baraffe, I., Clarke, J., Morison, A., Vlaykov, D.G., Constantino, T., Goffrey, T., Guillet, T., Le Saux, A., Pratt, J.: A study of convective core overshooting as a function of stellar mass based on two-dimensional hydrodynamical simulations. MNRAS 519(4), 5333–5344 (2023) arXiv:2301.02604 [astro-ph.SR]. https://doi.org/10.1093/mnras/stad009
  * (26) Rogers, T.M., Lin, D.N.C., McElwaine, J.N., Lau, H.H.B.: Internal Gravity Waves in Massive Stars: Angular Momentum Transport. ApJ 772(1), 21 (2013) arXiv:1306.3262 [astro-ph.SR]. https://doi.org/10.1088/0004-637X/772/1/21
  * (27) Rogers, T.M., McElwaine, J.N.: On the Chemical Mixing Induced by Internal Gravity Waves. ApJ 848(1), 1 (2017) arXiv:1709.04920 [astro-ph.SR]. https://doi.org/10.3847/2041-8213/aa8d13
  * (28) Couston, L.-A., Lecoanet, D., Favier, B., Le Bars, M.: The energy flux spectrum of internal waves generated by turbulent convection. Journal of Fluid Mechanics 854, 3 (2018) arXiv:1808.04638 [physics.flu-dyn]. https://doi.org/10.1017/jfm.2018.669
  * (29) Edelmann, P.V.F., Ratnasingam, R.P., Pedersen, M.G., Bowman, D.M., Prat, V., Rogers, T.M.: Three-dimensional Simulations of Massive Stars. I. Wave Generation and Propagation. ApJ 876(1), 4 (2019) arXiv:1903.09392 [astro-ph.SR]. https://doi.org/10.3847/1538-4357/ab12df
  * (30) Horst, L., Edelmann, P.V.F., Andrássy, R., Röpke, F.K., Bowman, D.M., Aerts, C., Ratnasingam, R.P.: Fully compressible simulations of waves and core convection in main-sequence stars. A&A 641, 18 (2020) arXiv:2006.03011 [astro-ph.SR]. https://doi.org/10.1051/0004-6361/202037531
  * (31) Lecoanet, D., Cantiello, M., Anders, E.H., Quataert, E., Couston, L.-A., Bouffard, M., Favier, B., Le Bars, M.: Surface manifestation of stochastically excited internal gravity waves. MNRAS 508(1), 132–143 (2021) arXiv:2105.04558 [astro-ph.SR]. https://doi.org/10.1093/mnras/stab2524
  * (32) Le Saux, A., Guillet, T., Baraffe, I., Vlaykov, D.G., Constantino, T., Pratt, J., Goffrey, T., Sylvain, M., Réville, V., Brun, A.S.: Two-dimensional simulations of solar-like models with artificially enhanced luminosity. II. Impact on internal gravity waves. A&A 660, 51 (2022) arXiv:2202.00801 [astro-ph.SR]. https://doi.org/10.1051/0004-6361/202142569
  * (33) Blouin, S., Mao, H., Herwig, F., Denissenkov, P., Woodward, P.R., Thompson, W.R.: 3D hydrodynamics simulations of internal gravity waves in red giant branch stars. arXiv e-prints, 2303–07332 (2023) arXiv:2303.07332 [astro-ph.SR]
  * (34) Samadi, R., Belkacem, K., Goupil, M.J., Dupret, M.-A., Brun, A.S., Noels, A.: Stochastic excitation of gravity modes in massive main-sequence stars. Ap&SS 328(1-2), 253–258 (2010). https://doi.org/10.1007/s10509-009-0215-3
  * (35) Shiode, J.H., Quataert, E., Cantiello, M., Bildsten, L.: The observational signatures of convectively excited gravity modes in main-sequence stars. MNRAS 430(3), 1736–1745 (2013) arXiv:1210.5525 [astro-ph.SR]. https://doi.org/10.1093/mnras/sts719
  * (36) Aerts, C., Rogers, T.M.: Observational Signatures of Convectively Driven Waves in Massive Stars. ApJ 806(2), 33 (2015) arXiv:1505.06648 [astro-ph.SR]. https://doi.org/10.1088/2041-8205/806/2/L33
  * (37) Lecoanet, D., Cantiello, M., Quataert, E., Couston, L.-A., Burns, K.J., Pope, B.J.S., Jermyn, A.S., Favier, B., Le Bars, M.: Low-frequency Variability in Massive Stars: Core Generation or Surface Phenomenon? ApJ 886(1), 15 (2019) arXiv:1910.01643 [astro-ph.SR]. https://doi.org/10.3847/2041-8213/ab5446
  * (38) Thompson, W., Herwig, F., Woodward, P.R., Mao, H., Denissenkov, P., Bowman, D.M., Blouin, S.: 3D hydrodynamic simulations of massive main-sequence stars II. Convective excitation and spectra of internal gravity waves. arXiv e-prints, 2303–06125 (2023) arXiv:2303.06125 [astro-ph.SR]. https://doi.org/10.48550/arXiv.2303.06125
  * (39) Jermyn, A.S., Anders, E.H., Lecoanet, D., Cantiello, M.: An Atlas of Convection in Main-sequence Stars. ApJS 262(1), 19 (2022) arXiv:2206.00011 [astro-ph.SR]. https://doi.org/10.3847/1538-4365/ac7cee
  * (40) Kurtz, D.W., Shibahashi, H., Murphy, S.J., Bedding, T.R., Bowman, D.M.: A unifying explanation of complex frequency spectra of $\gamma$ Dor, SPB and Be stars: combination frequencies and highly non-sinusoidal light curves. MNRAS 450(3), 3015–3029 (2015) arXiv:1504.04245 [astro-ph.SR]. https://doi.org/10.1093/mnras/stv868
  * (41) Pentcheva, B.V., Abel, J.S.: Icons of sound: Auralizing the lost voice of hagia sophia. Speculum 92(S1), 336–360 (2017) https://doi.org/10.1086/693439. https://doi.org/10.1086/693439
  * (42) Jermyn, A.S., Anders, E.H., Cantiello, M.: A Transparent Window into Early-type Stellar Variability. ApJ 926(2), 221 (2022) arXiv:2201.10567 [astro-ph.SR]. https://doi.org/10.3847/1538-4357/ac4e89
  * (43) Dotter, A.: MESA Isochrones and Stellar Tracks (MIST) 0: Methods for the Construction of Stellar Isochrones. ApJS 222(1), 8 (2016) arXiv:1601.05144 [astro-ph.SR]. https://doi.org/10.3847/0067-0049/222/1/8
  * (44) Choi, J., Dotter, A., Conroy, C., Cantiello, M., Paxton, B., Johnson, B.D.: Mesa Isochrones and Stellar Tracks (MIST). I. Solar-scaled Models. ApJ 823(2), 102 (2016) arXiv:1604.08592 [astro-ph.SR]. https://doi.org/10.3847/0004-637X/823/2/102
  * (45) Lecoanet, D., Quataert, E.: Internal gravity wave excitation by turbulent convection. MNRAS 430(3), 2363–2376 (2013) arXiv:1210.4547 [astro-ph.SR]. https://doi.org/%****␣article.bbl␣Line␣925␣****10.1093/mnras/stt055
  * (46) Le Saux, A., Baraffe, I., Guillet, T., Vlaykov, D.G., Morison, A., Pratt, J., Constantino, T., Goffrey, T.: Two-dimensional simulations of internal gravity waves in a 5 M⊙ zero-age-main-sequence model. MNRAS (2023) arXiv:2304.02508 [astro-ph.SR]. https://doi.org/10.1093/mnras/stad1067
  * (47) Augustson, K.C., Mathis, S., Astoul, A.: A Model of Rotating Convection in Stellar and Planetary Interiors. II. Gravito-inertial Wave Generation. ApJ 903(2), 90 (2020) arXiv:2009.10473 [astro-ph.SR]. https://doi.org/10.3847/1538-4357/abba1c
  * (48) Varghese, A., Ratnasingam, R.P., Vanon, R., Edelmann, P.V.F., Rogers, T.M.: Chemical Mixing Induced by Internal Gravity Waves in Intermediate-mass Stars. ApJ 942(1), 53 (2023) arXiv:2211.06432 [astro-ph.SR]. https://doi.org/10.3847/1538-4357/aca092
  * (49) De Cat, P., Aerts, C.: A study of bright southern slowly pulsating B stars. II. The intrinsic frequencies. A&A 393, 965–981 (2002). https://doi.org/10.1051/0004-6361:20021068
  * (50) Pedersen, M.G.: Internal Rotation and Inclinations of Slowly Pulsating B Stars: Evidence of Interior Angular Momentum Transport. ApJ 940(1), 49 (2022) arXiv:2208.14497 [astro-ph.SR]. https://doi.org/10.3847/1538-4357/ac947f
  * (51) Paxton, B., Bildsten, L., Dotter, A., Herwig, F., Lesaffre, P., Timmes, F.: Modules for Experiments in Stellar Astrophysics (MESA). ApJS 192(1), 3 (2011) arXiv:1009.1622 [astro-ph.SR]. https://doi.org/10.1088/0067-0049/192/1/3
  * (52) Paxton, B., Cantiello, M., Arras, P., Bildsten, L., Brown, E.F., Dotter, A., Mankovich, C., Montgomery, M.H., Stello, D., Timmes, F.X., Townsend, R.: Modules for Experiments in Stellar Astrophysics (MESA): Planets, Oscillations, Rotation, and Massive Stars. ApJS 208(1), 4 (2013) arXiv:1301.0319 [astro-ph.SR]. https://doi.org/10.1088/0067-0049/208/1/4
  * (53) Paxton, B., Marchant, P., Schwab, J., Bauer, E.B., Bildsten, L., Cantiello, M., Dessart, L., Farmer, R., Hu, H., Langer, N., Townsend, R.H.D., Townsley, D.M., Timmes, F.X.: Modules for Experiments in Stellar Astrophysics (MESA): Binaries, Pulsations, and Explosions. ApJS 220(1), 15 (2015) arXiv:1506.03146 [astro-ph.SR]. https://doi.org/10.1088/0067-0049/220/1/15
  * (54) Paxton, B., Schwab, J., Bauer, E.B., Bildsten, L., Blinnikov, S., Duffell, P., Farmer, R., Goldberg, J.A., Marchant, P., Sorokina, E., Thoul, A., Townsend, R.H.D., Timmes, F.X.: Modules for Experiments in Stellar Astrophysics (MESA): Convective Boundaries, Element Diffusion, and Massive Star Explosions. ApJS 234(2), 34 (2018) arXiv:1710.08424 [astro-ph.SR]. https://doi.org/10.3847/1538-4365/aaa5a8
  * (55) Paxton, B., Smolec, R., Schwab, J., Gautschy, A., Bildsten, L., Cantiello, M., Dotter, A., Farmer, R., Goldberg, J.A., Jermyn, A.S., Kanbur, S.M., Marchant, P., Thoul, A., Townsend, R.H.D., Wolf, W.M., Zhang, M., Timmes, F.X.: Modules for Experiments in Stellar Astrophysics (MESA): Pulsating Variable Stars, Rotation, Convective Boundaries, and Energy Conservation. ApJS 243(1), 10 (2019) arXiv:1903.01426 [astro-ph.SR]. https://doi.org/10.3847/1538-4365/ab2241
  * (56) Jermyn, A.S., Bauer, E.B., Schwab, J., Farmer, R., Ball, W.H., Bellinger, E.P., Dotter, A., Joyce, M., Marchant, P., Mombarg, J.S.G., Wolf, W.M., Wong, T.L.S., Cinquegrana, G.C., Farrell, E., Smolec, R., Thoul, A., Cantiello, M., Herwig, F., Toloza, O., Bildsten, L., Townsend, R.H.D., Timmes, F.X.: Modules for Experiments in Stellar Astrophysics (MESA): Time-Dependent Convection, Energy Conservation, Automatic Differentiation, and Infrastructure. arXiv e-prints, 2208–03651 (2022) arXiv:2208.03651 [astro-ph.SR]. https://doi.org/10.48550/arXiv.2208.03651
  * (57) Burns, K.J., Vasil, G.M., Oishi, J.S., Lecoanet, D., Brown, B.P.: Dedalus: A flexible framework for numerical simulations with spectral methods. Physical Review Research 2(2), 023068 (2020) arXiv:1905.10388 [astro-ph.IM]. https://doi.org/10.1103/PhysRevResearch.2.023068
  * (58) Vasil, G.M., Lecoanet, D., Burns, K.J., Oishi, J.S., Brown, B.P.: Tensor calculus in spherical coordinates using jacobi polynomials. part-i: Mathematical analysis and derivations. Journal of Computational Physics: X 3, 100013 (2019). https://doi.org/10.1016/j.jcpx.2019.100013
  * (59) Lecoanet, D., Vasil, G.M., Burns, K.J., Brown, B.P., Oishi, J.S.: Tensor calculus in spherical coordinates using jacobi polynomials. part-ii: Implementation and examples. Journal of Computational Physics: X 3, 100012 (2019). https://doi.org/10.1016/j.jcpx.2019.100012
  * (60) Townsend, R.H.D., Teitler, S.A.: GYRE: an open-source stellar oscillation code based on a new Magnus Multiple Shooting scheme. MNRAS 435(4), 3406–3418 (2013) arXiv:1308.2965 [astro-ph.SR]. https://doi.org/10.1093/mnras/stt1533
  * (61) Townsend, R.H.D., Goldstein, J., Zweibel, E.G.: Angular momentum transport by heat-driven g-modes in slowly pulsating B stars. MNRAS 475(1), 879–893 (2018) arXiv:1712.02420 [astro-ph.SR]. https://doi.org/10.1093/mnras/stx3142
  * (62) Townsend, R.H.D.: A semi-analytical formula for the light variations due to low-frequency g modes in rotating stars. MNRAS 343(1), 125–136 (2003). https://doi.org/10.1046/j.1365-8711.2003.06640.x
  * (63) Townsend, R.H.D.: Photometric modelling of slowly pulsating B stars. MNRAS 330(4), 855–875 (2002). https://doi.org/10.1046/j.1365-8711.2002.05135.x
  * (64) Townsend, R., Lopez, A.: Msg: A software package for interpolating stellar spectra in pre-calculated grids. Journal of Open Source Software 8(81), 4602 (2023). https://doi.org/10.21105/joss.04602
  * (65) Lanz, T., Hubeny, I.: A Grid of Non-LTE Line-blanketed Model Atmospheres of O-Type Stars. ApJS 146(2), 417–441 (2003) arXiv:astro-ph/0210157 [astro-ph]. https://doi.org/10.1086/374373
  * (66) Castelli, F., Kurucz, R.L.: New Grids of ATLAS9 Model Atmospheres. In: Piskunov, N., Weiss, W.W., Gray, D.F. (eds.) Modelling of Stellar Atmospheres, vol. 210, p. 20 (2003). https://doi.org/10.48550/arXiv.astro-ph/0405087
  * (67) Rodrigo, C., Solano, E.: The SVO Filter Profile Service. In: XIV.0 Scientific Meeting (virtual) of the Spanish Astronomical Society, p. 182 (2020)
  * (68) Hunter, J.D.: Matplotlib: A 2d graphics environment. Computing in Science & Engineering 9(3), 90–95 (2007). https://doi.org/10.1109/MCSE.2007.55
  * (69) Caswell, T.A., Droettboom, M., Lee, A., Sales De Andrade, E., Hoffmann, T., Klymak, J., Hunter, J., Firing, E., Stansby, D., Varoquaux, N., Hedegaard Nielsen, J., Root, B., May, R., Elson, P., Seppänen, J.K., Dale, D., Lee, J.-J., McDougall, D., Straw, A., Hobson, P., Hannah, Gohlke, C., Vincent, A.F., Yu, T.S., Ma, E., Silvester, S., Moad, C., Kniazev, N., Ernest, E., Ivanov, P.: matplotlib/matplotlib: REL: v3.5.2. Zenodo (2022). https://doi.org/10.5281/zenodo.6513224
  * (70) Harris, C.R., Millman, K.J., van der Walt, S.J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N.J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M.H., Brett, M., Haldane, A., del Río, J.F., Wiebe, M., Peterson, P., Gérard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., Oliphant, T.E.: Array programming with NumPy. Nature 585(7825), 357–362 (2020). https://doi.org/10.1038/s41586-020-2649-2
  * (71) Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S.J., Brett, M., Wilson, J., Millman, K.J., Mayorov, N., Nelson, A.R.J., Jones, E., Kern, R., Larson, E., Carey, C.J., Polat, İ., Feng, Y., Moore, E.W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E.A., Harris, C.R., Archibald, A.M., Ribeiro, A.H., Pedregosa, F., van Mulbregt, P., SciPy 1.0 Contributors: SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods 17, 261–272 (2020). https://doi.org/10.1038/s41592-019-0686-2
  * (72) Astropy Collaboration: Astropy: A community Python package for astronomy. A&A 558, 33 (2013) arXiv:1307.6212 [astro-ph.IM]. https://doi.org/10.1051/0004-6361/201322068
  * (73) Astropy Collaboration: The Astropy Project: Building an Open-science Project and Status of the v2.0 Core Package. AJ 156(3), 123 (2018) arXiv:1801.02634 [astro-ph.IM]. https://doi.org/10.3847/1538-3881/aabc4f
  * (74) Astropy Collaboration: The Astropy Project: Sustaining and Growing a Community-oriented Open-source Project and the Latest Major Release (v5.0) of the Core Package. apj 935(2), 167 (2022) arXiv:2206.14220 [astro-ph.IM]. https://doi.org/10.3847/1538-4357/ac7c74
  * (75) Anders, E.H., Lecoanet, D., Cantiello, M., Burns, K.J., Hyatt, B.A., Kaufman, E., Townsend, R.H.D., Brown, B.P., Vasil, G.M., Oishi, J.S., Jermyn, A.S.: Code, Data, and Supplemental Materials for ”The photometric variability of massive stars due to gravity waves excited by core convection”. Zenodo (2023). https://doi.org/10.5281/zenodo.7764997. https://doi.org/10.5281/zenodo.7764997