# Asteroseismology of evolved stars to constrain the internal transport of angular momentum

J. den Hartogh, P. Eggenberger, S. Deheuvels

<figure><img src="image_1.png" /></figure>

## To cite this version:

J. den Hartogh, P. Eggenberger, S. Deheuvels. Asteroseismology of evolved stars to constrain the internal transport of angular momentum. Astronomy & Astrophysics - A&A, 2020, 634, pp.L16. 10.1051/0004-6361/202037568. [hal-04938492](https://hal.science/hal-04938492v1)

HAL Id: hal-04938492

https://hal.science/hal-04938492v1

Submitted on 9 Aug 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
LETTER TO THE EDITOR

# Asteroseismology of evolved stars to constrain the internal transport of angular momentum

## III. Using the rotation rates of intermediate-mass stars to test the Fuller-formalism

J. W. den Hartogh<fnref n="1" />,<fnref n="2" />,*, P. Eggenberger<fnref n="3" />, and S. Deheuvels<fnref n="4" />

Received 24 January 2020 / Accepted 7 February 2020

## ABSTRACT

*Context.* The internal characteristics of stars, such as their core rotation rates, are obtained via asteroseismic observations. A comparison of core rotation rates found in this way with core rotation rates as predicted by stellar evolution models demonstrate a large discrepancy. This means that there must be a process of angular momentum transport missing in the current theory of stellar evolution. A new formalism was recently proposed to fill in for this missing process, which has the Tayler instability as its starting point (by Fuller et al. 2019, MNRAS, 485, 3661, hereafter referred to as “Fuller-formalism”).

*Aims.* We investigate the effect of the Fuller-formalism on the internal rotation of stellar models with an initial mass of 2.5 $M_\odot$.

*Methods.* Stellar evolution models, including the Fuller-formalism, of intermediate-mass stars were calculated to make a comparison between asteroseismically obtained core rotation rates in the core He burning phase and in the white dwarf phase.

*Results.* Our main results show that models including the Fuller-formalism can match the core rotation rates obtained for the core He burning phases. However, these models are unable to match the rotation rates obtained for white dwarfs. When we exclude the Fuller-formalism at the end of the core He burning phase, the white dwarf rotation rates of the models match the observed rates.

*Conclusions.* We conclude that in the present form, the Fuller-formalism cannot be the sole solution for the missing process of angular momentum transport in intermediate-mass stars.

**Key words.** stars: evolution – stars: rotation – stars: interiors – stars: oscillations

## 1. Introduction

Recent developments in asteroseismology have led to a wealth of information on the internal properties of mainly low- and intermediate-mass stars (recently reviewed in Aerts et al. 2019). In particular, the analysis of mixed modes have been used to infer core rotation rates in post-main sequence stars (Beck et al. 2012; Deheuvels et al. 2012, 2014, 2015, 2017; Mosser et al. 2012; Di Mauro et al. 2016, 2018; Gehan et al. 2018; Tayar et al. 2019). A large discrepancy exists between the predicted values of stellar evolution models and observed values (Eggenberger et al. 2012; Ceillier et al. 2013; Marques et al. 2013; Cantiello et al. 2014), which confirms the findings of earlier works (Pinsonneault et al. 1989; Chaboyer et al. 1995; Eggenberger et al. 2005; Suijs et al. 2008; Denissenkov et al. 2010).

This discrepancy between predicted and observed core rotation rates shows that an additional, dominant process that transports angular momentum (AM) from the core to the outer layers is active in stars, but missing from stellar evolution codes. There is still no consensus on the physics behind this discrepancy.

Several physical processes are currently being investigated. One is presented in the work by Kissin & Thompson (2015), who argues that magnetic torques leads to solid-body rotation in radiative regions, and differential rotation in convective regions. Another physical process being investigated is transport of AM by internal gravity waves; see Rogers et al. (2013), and an extensive overview in Aerts et al. (2019). Magnetised winds have been tested by Tayar & Pinsonneault (2018) and Tayar et al. (2019). Recently, a re-investigation of the Tayler instability was published by Fuller et al. (2019) (hereafter referred to as the “Fuller-formalism”). There is no consensus on which of these options, or whether any of them at all, are actually present in stars and acts as a dominant process in terms of AM transport.

In this study, we focus on testing the Fuller-formalism on a small grid of intermediate-mass stars. Previous tests of this formalism highlighted that it indeed transports AM and is dominant over the meridional and shear instabilities Eggenberger et al. (2019a,b). These tests, which have only concerned low-mass stars, show that the Fuller-formalism in its current form is not the full solution to the missing process of AM transport. In particular, Eggenberger et al. (2019b) shows that the Fuller-formalism has trouble reproducing the core rotation rates of subgiant stars

<footnote n="1">¹ Konkoly Observatory, CSFK, Konkoly Thege Miklós út 15-17, 1121 Budapest, Hungary
e-mail: jacqueline.den.hartogh@csfk.mta.hu</footnote>

<footnote n="2">² Astrophysics Group, Lennard-Jones Laboratories, Keele University, ST5 5BG Newcastle, UK</footnote>

<footnote n="3">³ Observatoire de Genève, Université de Genève, 51 Ch. des Maillettes, 1290 Sauverny, Switzerland</footnote>

<footnote n="4">⁴ IRAP, Université de Toulouse, CNRS, CNES, UPS, Toulouse, France</footnote>

<footnote>* NuGrid Collaboration, http://www.NuGridstars.org</footnote>
and argues that different values of free parameter $\alpha$ are needed to match consecutive stellar evolutionary phases.

In this study, we focus on the comparison of the Fuller-formalism with the asteroseismic observations of intermediate-mass stars (Deheuvels et al. 2015). The seven core He burning stars in that study all show weak radial differential rotation and slow core rotation rates. Recently, these findings were confirmed to be a general characteristic of core He burning stars with a mass between 2 and 3 $M_\odot$ and not merely an observational bias, by Tayar et al. (2019). This particular initial mass range is important as it is known to be the site for the main component of the slow neutron capture process and, thus, responsible for about 50% of all elements heavier than iron (Herwig 2005; Käppeler et al. 2011; Karakas & Lattanzio 2014). Finally, many of the known white dwarfs whose rotational periods have been measured (Kawaler 2015; Hermes et al. 2017) are believed to have progenitors in the same mass range of 2–3 $M_\odot$. In the following section we introduce the physics of our stellar evolution models and our results are shown in Sect. 3. Our conclusions are summarised in Sect. 4.

## 2. Stellar evolution models

The stellar evolution models are generated with MESA revision 8845 Paxton et al. (2011, 2013, 2015). The initial settings are the same as in den Hartogh et al. (2019a,b), as is the method we used to get through the convergence issues that can occur at the end of the asymptotic giant branch (AGB) phase. The implementation of the Fuller-formalism is the same as in Fuller et al. (2019), so that only the formalism’s effect on the transport of angular momentum is included in the stellar evolution models and not its effect on the mixing of chemical elements. No other processes linked to rotation are included in the calculations in our approach, which is also in agreement with Fuller et al. (2019).

In Spruit (2002), the magnetic torque directly follows from the physics of the Tayler instability itself as it interacts with differential rotation. In Fuller et al. (2019), the torque follows from a turbulence formalism that focusses on energy dissipation. The resulting energy dissipation in the Fuller-formalism is then lower than in the TS-dynamo, leading to stronger magnetic fields and Maxwell stresses. As a result, the Fuller-formalism is more efficient in transporting AM than the TS-dynamo. The equations for the transport of AM in the Fuller-formalism are:

$$v_F = \alpha^3 r^2 \Omega \left(\frac{\Omega}{N_{\text{eff}}}\right)^2, \qquad (1)$$

$$N_{\text{eff}}^2 = \frac{\eta}{K} N_T^2 + N_\mu^2, \qquad (2)$$

$$q_{\text{min}} = \alpha^{-3} \left(\frac{N_{\text{eff}}}{\Omega}\right)^{5/2} \left(\frac{\eta}{r^2 \Omega}\right)^{3/4}, \qquad (3)$$

with $\alpha$ being a dimensionless factor (discussed below), $\eta$ and $K$ the magnetic and thermal diffusivities, $N_{\text{eff}}$ the effective Brunt-Väisälä frequency, and $q_{\text{min}}$ the new minimum threshold for shear ($q = -\frac{\partial \ln \Omega}{\partial \ln r}$). The Fuller-formalism is only active and transporting AM when the shear is larger than $q_{\text{min}}$.

$$\left(\frac{\partial \Omega}{\partial t}\right)_m = \frac{1}{j} \left(\frac{\partial}{\partial m}\right)_t \left[ (4\pi r^2 \rho)^2 j D_{\text{am}} \left(\frac{\partial \Omega}{\partial m}\right) - \frac{2\Omega}{r} \left(\frac{\partial r}{\partial t}\right)_m \left(\frac{1}{2} \frac{d\ln j}{d\ln r}\right) \right], \qquad (4)$$

where $\Omega$ is the angular velocity, $j$ the specific angular momentum, and the total diffusion coefficient for AM transport is $D_{\text{am}}$.

Following Fuller et al. (2019), (1) the only process included in the AM transport is the Fuller-formalism and, thus, $D_{\text{am}} = v_F$; and (2) there is no mixing of chemical elements due to rotation included in the models presented here.

The small set of models calculated for this work have an initial mass of 2.5 $M_\odot$, an initial metallicity of $Z = 0.01$, and an initial rotation rate of 50 km s$^{-1}$; this follows a direct continuation of den Hartogh et al. (2019a). The Fuller-formalism has only one free parameter<fnref n="1" />: $\alpha$, which we vary in order to match the core rotation rates of the seven stars analysed in Deheuvels et al. (2015).

These seven core He burning stars were chosen because both the core and surface rotation rate were known, which allows for an improved constraining of the initial parameters. The small size of this sample made it difficult to determine whether it was sufficiently representative of intermediate-mass core He burning stars, especially as their core rotation rates are used to place them among the slowest rotating core He burning stars, see Fig. 4 in Aerts et al. (2019). Recently, 99 new core rotation rates for core He burning stars were published by Tayar et al. (2019), in which most stars have a mass of 2–3 solar masses. The core rotation rates of the seven stars of Deheuvels et al. (2015) fit well in the new data set and they are, thus, representative values. We include all the known white dwarfs rotation rates, as presented in Kawaler (2015) and Hermes et al. (2017).

## 3. Results

The evolution of the core rotation rates of the 2.5 $M_\odot$ models, including the Fuller-formalism, is shown in the top panel of Fig. 1. The core rotation rates correspond to the mean value determined in the g-mode cavity (following Goupil et al. 2013) and are plotted against the surface gravity. Black dots are added to indicate the core He burning phase, where interval between dots accounts for 10% of the total duration of the core He burning phase. The first dot (located around $\log_{10}(g/cm \text{s}^{-2}) \approx 1.8$) indicates the start of the core He burning phase, while the final black dot (located around $\log_{10}(g/cm \text{s}^{-2}) \approx 2.3$) indicates the end of the core He burning phase. These markers show that the stellar models spend the majority of the phase in the part of the parameter space of $\Omega_c$, where the observed rotation rates are located. A slightly lower initial mass is needed in the stellar evolution models to reach all the corresponding $\log_{10}g$ values.

In order to reach the values by Deheuvels et al. (2015), we had to include a value for $\alpha$ of two to four. This increase enhances the transport of angular momentum, leading to lower core rotation rates. The grey data points in the top panel are the core rotation rates of stars with a mass between 2 and 3 $M_\odot$ of Mosser et al. (2012) and Tayar et al. (2019), which cover a similar space in the figure. This indicates that these values for $\alpha$ are needed for all stars in this mass range, possibly even higher values for the data points just below $\Omega_c = 10^2$ nHz.

The increase of $\alpha$ also results in a decrease of the differential rotation, which improves the match of the models to the ratio of core to surface rotation rate, as shown in the bottom panel of Fig. 1. The need for these higher values of $\alpha$ is in contrast to Fuller et al. (2019), who suggests that $\alpha$ should be around 1.

These results are different in several aspects from the ones obtained by den Hartogh et al. (2019a), using a constant additional, artificial viscosity to match the seven stars from Deheuvels et al. (2015). The models with the Fuller-formalism experience solid body rotation until a $\log_{10}(g/cm \text{s}^{-2}) \approx 2.9$ is

<footnote n="1">¹ There are a few more free parameters, but those only concern the smoothing of $v_F$. Their values are the same as in Fuller et al. (2019).</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Comparison of models with different alphas to all observed core rotation rates (top panel) within 2–3 $M_{\odot}$ range. Highlighted in colour are the core He burning stars of Deheuvels et al. (2015), showing the core (red diamonds) and surface (orange stars) values in the top panel and the ratio of the two values (purple squares) in the bottom panel. The core rotation rates of Mosser et al. (2012) and Tayar et al. (2019) are shown in grey. Values of $\alpha > 1$ are needed to reach the observed range of core rotation rates and degree of differential rotation. One surface rotation rate is plotted (dotted line) in the left panel. The surface rotation rates of the other models are identical to the one plotted and, therefore, they are not included in this figure.</figcaption></figure>

reached, while the models in den Hartogh et al. (2019a) diverge from solid body rotation at the end of the main sequence, at $\log_{10}(g/cm\ s^{-2}) \approx 4.0$ (see their Figs. A.2 and A.3). The elon- gation of the solid body phase is also visible in Fig. 3 of Fuller et al. (2019). Another difference is, of course, that the vis- cosity due to the Fuller-formalism is not a constant. When com- paring our results to Fuller et al. (2019), and specifically their Fig. 4, some differences can be identified. In particular, in Fig. 4 of Fuller et al. (2019), the range of clump rotational periods is matched with $\alpha$ of 0.5–2, while we need a larger range for $\alpha$ to match all observed values. This can be explained by two points: the models presented here are calculated with a higher initial mass than in Fuller et al. (2019): namely, here it is $2.5 M_{\odot}$, there it is $1.5 M_{\odot}$. And, as shown in Fig. 6 of Fuller et al. (2019), their models (all using $\alpha = 1$) with initial mass of 2.2 and $3.0 M_{\odot}$ range show lower periods during their clump phase when com- pared to their models with an initial mass of 1.2 and $2.0 M_{\odot}$ (respectively, 40 days and 100–200 days). These lower periods would not match the observed range as included in their Fig. 4, which is in agreement with our results.

The stellar evolution models are continued into the white dwarf phase by using the same mass loss parameters as in

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Comparison of observed white dwarf rotational periods and white dwarf periods from models presented here (star symbols). Only the period of the $\alpha = 1$ model falls within the observed range of peri- ods, but this model does not reach the observed core rotation rates in the core He burning phase. All other stellar evolution models reach the DAV temperature range with a rotational period higher than those observed.</figcaption></figure>

den Hartogh et al. (2019a,b). The comparison between all known white dwarf rotational periods (Kawaler 2015; Hermes et al. 2017) and the rotational periods of the models when reaching the DAV temperature range are shown in Fig. 2. This histogram shows that most observed white dwarfs have rotational periods of two days or less, with the longest observed rotational period being 109.1 h. We note that the observation method of Hermes et al. (2017) is sensitive to larger rotation periods up to 15 days, but the longest period they found was about 4.5 days.

The white dwarf rotational periods of our models start at higher values. Specifically, the stellar models that match the rota- tion rates of Deheuvels et al. (2015), with $\alpha = 2-4$, lead to white dwarf rotational periods that are all higher than the observed periods.

Similarly to the case with our results of the core He burn- ing phase, we compare this result to the results presented in Fuller et al. (2019). In their Fig. 7, they compare their white dwarf rotational periods (all with $\alpha = 1$) to the observed values. Our results agree that these models match the observations, however, as we explained before: the models with $\alpha = 1$ do not match the low rotation rates of the clump stars analysed by Deheuvels et al. (2015), and other clump stars in the same initial mass range of $2-3 M_{\odot}$ by Mosser et al. (2012) and Tayar et al. (2019).

We further investigate the post-core He burning evolution of the core rotation rates by analysing the $\Omega$ and specific angular momentum $j$ profiles throughout the stellar evolution. These pro- files are shown in Fig. 3, from the start of the main sequence “start MS” until the final model number in the calculation “Final”, which is after the calculation has gone through the DAV temperature range. The transport of AM due to the Fuller- formalism leads to the decrease of a factor of 15 in the core rotation rate (left panel of Fig. 3) between the end of the main sequence and the start of core He burning. This decrease is also visible the profiles of $j$. In contrast to the findings of Cantiello et al. (2014) and den Hartogh et al. (2019a) however, there is also AM transport happening after the core He burn- ing phase. This decrease of $j$ in the core of about one order of magnitude reduces the core rotation rate and, thus, increases the rotational periods we calculate for the white dwarfs. We there- fore calculate stellar evolutionary models that include the Fuller- formalism until the end of the core He burning phase and exclude it going forward. The resulting white dwarf rotational periods
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Profiles of $\Omega$ and $j$ for stellar evolutionary phases taken from the model with $\alpha = 4$. The eight phases are the start and end of the main sequence (respectively “start MS” and “end MS”), the start and end of the core He burning phases (“start He burn”, and “end He burn”), the first thermal pulse in the AGB phase (“AGB TP1”), the end of the AGB phase (“Envelope gone”), the start of the white dwarf phase (“WD”), and the final profile in the calculation (“Final”), which is located lower on the white dwarf cooling track than the DAV phase. The figures show that there is AM transport after the core He burning phases.</figcaption></figure>

<table><caption>**Table 1.** Rotational periods of our white dwarf models, calculated in the DAV temperature range.</caption><thead><tr><th>α</th><th>P<sub>DAV</sub> (h)</th></tr></thead><tbody><tr><td>1</td><td>81.0</td></tr><tr><td>2</td><td>143</td></tr><tr><td>3</td><td>212</td></tr><tr><td>4</td><td>280</td></tr><tr><td>1*</td><td>28.4</td></tr><tr><td>4*</td><td>70.0</td></tr></tbody></table>

**Notes.** The first four models are also shown in Fig. 2, the final two models include the Fuller-formalism from the main sequence to the end of core He burning, and exclude the formalism from then onward. The latter two periods are both within the observed range.

for stellar evolution models with $\alpha = 1$ and $\alpha = 4$ are shown in Table 1 (as “1*” and “4*”, respectively). These new periods are 35% and 25%, respectively, of the original period calculated with the inclusion of the Fuller-formalism throughout the evolution. This decrease is large enough for the new periods to fall within the observed range of white dwarf rotational periods.

## 4. Conclusions

In this Letter, we test the Fuller-formalism in stellar evolutionary models with an initial mass of $2.5 M_{\odot}$. We investigate whether this formalism can lead to core rotation rates that match the observational range of rotation rates during the core He burning phase and the DAV temperature range. Our findings are:

- When including the Fuller-formalism in our stellar models, we need an $\alpha$ value of at least 2 to reach the core rotation rates of the intermediate stars analysed by Deheuvels et al. (2015), while lower values are needed to match core rotation rates of low mass star (Eggenberger et al. 2019b);

- After continuing with these stellar models with $\alpha = 2-4$ to the white dwarf phase, we found that the resulting white dwarf rotational periods are in the range of 81–280 h;

- The models that match the range of core He burning stars of Deheuvels et al. (2015) do not match the observed white dwarf rotation rates, and vice versa;

- The range of observed white dwarf rotation rates is matched when the Fuller-formalism is only included from the start of the main-sequence to the end of core He burning.

Therefore, it is important to use caution when adopting the Fuller-formalism in stellar evolution models. The formalism, based on Eggenberger et al. (2019a,b) and this work, is not likely to be the sole missing piece in the process of AM transport, but could still could be part of the solution. It is possible that a physical reason behind the exclusion of the Fuller-formalism at the end of core He burning may be found, which would lead to a more favourable judgement of the formalism when included in stellar evolution models of intermediate-mass stars.

**Acknowledgements.** We thank Henk Spruit for useful discussions. This work has been supported by the European Research Council (ERC-2016-CO Grant 724560), and the EU COST Action CA16117 (ChETEC). PE acknowledges support by the European Research Council (ERC grant agreement No. 833925). SD acknowledges support from the project BEAMING ANR-18-CE31-0001 of the French National Research Agency (ANR).

## References

Aerts, C., Mathis, S., & Rogers, T. M. 2019, ARA&A, 57, 35

Beck, P. G., Montalban, J., Kallinger, T., et al. 2012, Nature, 481, 55

Cantiello, M., Mankovich, C., Bildsten, L., Christensen-Dalsgaard, J., & Paxton, B. 2014, ApJ, 788, 93

Ceillier, T., Eggenberger, P., García, R. A., & Mathis, S. 2013, A&A, 555, A54

Chaboyer, B., Demarque, P., & Pinsonneault, M. H. 1995, ApJ, 441, 865

Deheuvels, S., García, R. A., Chaplin, W. J., et al. 2012, ApJ, 756, 19

Deheuvels, S., Doğan, G., Goupil, M. J., et al. 2014, A&A, 564, A27

Deheuvels, S., Ballot, J., Beck, P., & Mosser, B. 2015, A&A, 580, A96

Deheuvels, S., Ouazzani, R. M., & Basu, S. 2017, A&A, 605, A75

den Hartogh, J. W., Eggenberger, P., & Hirschi, R. 2019a, A&A, 622, A187

den Hartogh, J. W., Hirschi, R., Lugaro, M., et al. 2019b, A&A, 629, A123

Denissenkov, P. A., Pinsonneault, M., Terndrup, D. M., & Newsham, G. 2010, ApJ, 716, 1269

Di Mauro, M. P., Ventura, R., Cardini, D., et al. 2016, ApJ, 817, 65

Di Mauro, M. P., Ventura, R., Corsaro, E., & Lustosa De Moura, B. 2018, ApJ, 862, 9

Eggenberger, P., Maeder, A., & Meynet, G. 2005, A&A, 440, L9

Eggenberger, P., Montalbán, J., & Miglio, A. 2012, A&A, 544, L4

Eggenberger, P., Buldgen, G., & Salmon, S. J. A. J. 2019a, A&A, 626, L1

Eggenberger, P., den Hartogh, J. W., Buldgen, G., et al. 2019b, A&A, 631, L6

Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS, 485, 3661

Gehan, C., Mosser, B., Michel, E., Samadi, R., & Kallinger, T. 2018, A&A, 616, A24

Goupil, M. J., Mosser, B., Marques, J. P., et al. 2013, A&A, 549, A75

Hermes, J. J., Gänsicke, B. T., Kawaler, S. D., et al. 2017, ApJS, 232, 23

Herwig, F. 2005, ARA&A, 43, 435

Käppeler, F., Gallino, R., Bisterzo, S., & Aoki, W. 2011, Rev. Mod. Phys., 83, 157

Karakas, A. I., & Lattanzio, J. C. 2014, PASA, 31, e030

Kawaler, S. D. 2015, in 19th European Workshop on White Dwarfs, eds. P. Dufour, P. Bergeron, & G. Fontaine, ASP Conf. Ser., 493, 65

Kissin, Y., & Thompson, C. 2015, ApJ, 808, 35

Marques, J. P., Goupil, M. J., Lebreton, Y., et al. 2013, A&A, 549, A74

Mosser, B., Goupil, M. J., Belkacem, K., et al. 2012, A&A, 548, A10

Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3

Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4

Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15

Pinsonneault, M. H., Kawaler, S. D., Sofia, S., & Demarque, P. 1989, ApJ, 338, 424

Rogers, T. M., Lin, D. N. C., McElwaine, J. N., & Lau, H. H. B. 2013, ApJ, 772, 21

Spruit, H. C. 2002, A&A, 381, 923

Suijs, M. P. L., Langer, N., Poelarends, A.-J., et al. 2008, A&A, 481, L87

Tayar, J., & Pinsonneault, M. H. 2018, ApJ, 868, 150

Tayar, J., Beck, P. G., Pinsonneault, M. H., García, R. A., & Mathur, S. 2019, ApJ, 887, 203