DOI: 10.5281/zenodo.18603404
title: Pulsational Instability of Quasi-Stars: Interpreting the Variability of Little Red Dots

description: # Pulsational Instability of Quasi-Stars


Companion data and code repository for:


> **Pulsational Instability of Quasi-Stars: Interpreting the Variability of Little Red Dots**

>

> Matteo Cantiello, Jake B. Hassan, Rosalba Perna, Philip J. Armitage, Mitchell C. Begelman, Yan-Fei Jiang, Taeho Ryu, Richard H. D. Townsend

>

> Submitted to ApJL ([arXiv:2512.17997](https://arxiv.org/abs/2512.17997))


## Overview


This repository contains the MESA stellar evolution models, GYRE pulsation analysis, and Jupyter notebooks needed to reproduce the results of the paper. We model quasi-stars -- supermassive stellar envelopes powered by accretion onto a central black hole -- across a mass range of 10^4 to 5x10^5 solar masses. We perform linear non-adiabatic pulsation calculations to identify a "Quasi-Star Instability Strip" where the kappa-mechanism drives radial pulsations with periods of ~20-180 years.


## Repository Structure


```

quasistars-zenodo/

|

|-- quasistars-mesa-template/ MESA work directory template for quasi-star models

|-- pulsations/ GYRE pulsation analysis and Jupyter notebooks

|-- figures/ Publication figures (PDF)

|

|-- 1e4Msun/ Quasi-star model: M = 10^4 Msun

|-- 2e4MsunIC/ Quasi-star model: M = 2x10^4 Msun

|-- 5e4MsunIC/ Quasi-star model: M = 5x10^4 Msun

|-- 1e5Msun/ Quasi-star model: M = 10^5 Msun

|-- 2e5Msun/ Quasi-star model: M = 2x10^5 Msun

|-- 5e5MsunIC/ Quasi-star model: M = 5x10^5 Msun

|

|-- 2e5Msun_283000_567/ MESA hydro run: unstable model (Teff ~ 4682 K)

|-- 2e5Msun_33000_67/ MESA hydro run: stable model (Teff ~ 5445 K)

```


## Software Requirements


- **MESA** (Modules for Experiments in Stellar Astrophysics), version `r24.08.1`

- https://docs.mesastar.org

- The MESA SDK is required to compile and run MESA.

- **GYRE** version 8.1 (stellar oscillation code, included in MESA)

- https://gyre.readthedocs.io

- **Python 3** with the following packages:

- `numpy`, `scipy`, `pandas`, `matplotlib`, `seaborn`, `h5py`

- `mesa_reader` -- for reading MESA output files

- `pygyre` -- for reading GYRE HDF5 output

- `astropy` -- for physical constants and units


## Quasi-Star MESA Template


The `quasistars-mesa-template/` directory is a complete MESA work directory configured to evolve a quasi-star model. Key features of the setup:


- **Central black hole**: modeled as a point mass (`M_center`) that grows by accreting from the stellar envelope. The black hole luminosity `L_center` is set by the accretion rate and a radiative efficiency parameter (`x_ctrl(1)`).

- **TOV correction**: General-relativistic corrections to gravity are enabled via `use_other_cgrav = .true.` (see `src/other_cgrav.inc`).

- **Custom physics**: The Fortran source in `src/run_star_extras.f90` implements the quasi-star energy source, black hole growth, and optional eruptive mass loss.

- **GYRE output**: MESA is configured to write GYRE-compatible pulse data with each saved profile (`write_pulse_data_with_profile = .true.`).

- **Zero metallicity**: Models use primordial composition (Z = 0, Y = 0.3).


### Running a Model


The simulation proceeds in two stages:


1. **Initial model** (`inlist_initial`): Creates a pre-main-sequence model of the desired mass and saves it as `start.mod`.

2. **Evolution** (`inlist_evolve`): Loads `start.mod`, relaxes a central point mass (the seed black hole), and evolves the quasi-star.


To change the quasi-star mass, edit `initial_mass` in `inlist_common`. To change the initial black hole mass, edit `new_core_mass` in `inlist_evolve`. For example:


```fortran

! In inlist_common:

initial_mass = 1d5 ! Total quasi-star mass in Msun


! In inlist_evolve:

new_core_mass = 1d2 ! Initial seed BH mass in Msun

```


To compile and run:


```bash

cd quasistars-mesa-template

./mk # Compile

./rn # Run (uses inlist_initial then inlist_evolve)

```


### Key Control Parameters (in `inlist_common`)


| Parameter | Description |

|-----------|-------------|

| `x_ctrl(1)` | Radiative efficiency of BH accretion (default: 0.1) |

| `x_ctrl(2)` | Convective efficiency parameter (default: 0.1) |

| `x_ctrl(4)` | Efficiency factor for eruptive mass loss (default: 0.1) |

| `x_logical_ctrl(1)` | Enable accretion onto the black hole |

| `x_logical_ctrl(2)` | Enable mass loss from rest-energy conversion |


## Pre-Computed MESA Models


Each model directory (e.g., `1e5Msun/`) contains:


- `LOGS/history.data` -- time-resolved stellar evolution output (luminosity, Teff, BH mass, etc.)

- `LOGS/profiles.index` -- maps model numbers to profile numbers

- `LOGS/profileN.data` -- detailed stellar structure snapshots

- `LOGS/profileN.data.GYRE` -- GYRE-compatible pulse data files

- `instability_strip_filtered.csv` -- GYRE pulsation results (growth rates, periods) for selected profiles


The two MESA hydrodynamic runs (`2e5Msun_283000_567/` and `2e5Msun_33000_67/`) are restarts from specific snapshots of the `2e5Msun` model, evolved with MESA's implicit hydrodynamics to verify the pulsational instability predicted by linear analysis.


## Pulsation Analysis


The `pulsations/` directory contains the analysis pipeline and four Jupyter notebooks that reproduce the paper figures.


### Notebooks


| Notebook | Description | Paper Figures |

|----------|-------------|---------------|

| `quasistar_HRD.ipynb` | HR diagram with evolutionary tracks colored by M_BH | Fig. 1 |

| `pulsation_analysis.ipynb` | GYRE mode analysis: eigenfunctions, work integrals, opacity derivatives | Figs. 2, 3 |

| `instability_strip.ipynb` | Constructs the quasi-star instability strip with RBF interpolation; growth rates and periods vs. Teff | Figs. 4, 5 |

| `MESA_hydro.ipynb` | Analyzes MESA hydrodynamic runs: HRD loops, radius pulsations, surface Mach number | Fig. 6 |


### GYRE Configuration


The file `gyre.in` is the GYRE input template used for both adiabatic and non-adiabatic radial mode calculations (l=0). Key settings:


- Frequency scan: 0.01 to 20 cycles/day (real axis), -10 to 2 (imaginary axis)

- Solver: `MAGNUS_GL2` differencing scheme

- Non-adiabatic search seeded from adiabatic solutions (`nad_search = 'AD'`)

- Output includes work integrals (`dW_dx`), opacity derivatives (`kap_T`, `kap_rho`), and eigenfunctions


### Batch GYRE Processing


The script `run_gyre_batch.py` automates running GYRE across many MESA profiles:


```python

import run_gyre_batch as gb


# List of GYRE-compatible profile files

file_list = ['path/to/profile10.data.GYRE', 'path/to/profile20.data.GYRE', ...]


# Run GYRE on each profile and save results

gb.process_gyre_profiles(file_list, output_filename='instability_strip_filtered.csv')

```


**Note**: You must edit the `GYRE_CMD` path in `run_gyre_batch.py` to point to your local GYRE installation.


### Supporting Python Files


- `astro_setup.py` -- imports, physical constants, and plot styling (loaded by all notebooks via `from astro_setup import *`)

- `Constants.py` -- CGS physical constants class

- `strip.py` -- instability strip analysis utilities


### Pre-Computed Data Files


- `summary_ad.h5`, `summary_nad.h5` -- GYRE adiabatic and non-adiabatic mode summaries for the 2x10^5 Msun reference model

- `mode_nad.l0.n+*.h5` -- detailed non-adiabatic eigenfunctions for individual modes

- `instability_strip.csv` -- compiled pulsation data across all models

- `quasistars_data.pkl` -- cached processed data from all MESA models (used by `instability_strip.ipynb` to avoid re-reading large ASCII files; set `LOAD_FROM_PICKLE = False` to regenerate)

- `hrd_data_cache.pkl` -- cached HRD track data (used by `quasistar_HRD.ipynb`)


## Reproducing the Paper


1. **Install prerequisites**: MESA r24.08.1, Python packages listed above.

2. **Inspect pre-computed models**: The MESA model outputs are included in this repository. To re-run a model from scratch, copy `quasistars-mesa-template/` to a new directory, adjust the mass and BH parameters, compile with `./mk`, and run with `./rn_nomodfiles inlist_evolve_header`

3. **Run the notebooks**: Open the Jupyter notebooks in `pulsations/` and execute them in the following order:

- `quasistar_HRD.ipynb` (Fig. 1)

- `pulsation_analysis.ipynb` (Figs. 2-3)

- `instability_strip.ipynb` (Figs. 4-5)

- `MESA_hydro.ipynb` (Fig. 6)

4. **Re-run GYRE analysis** (optional): To recompute the pulsation mode data from the MESA profiles, use `run_gyre_batch.py` as described above, or run GYRE manually with the provided `gyre.in` template.