# Learning Optimization Proxies for Large-Scale Security-Constrained Economic
Dispatch

Wenbo Chen, Seonho Park, Mathieu Tanneau, Pascal Van Hentenryck  Georgia
Institute of Technology
{wenbo.chen, spark719<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The Security-Constrained Economic Dispatch (SCED) is a fundamental
optimization model for Transmission System Operators (TSO) to clear real-time
energy markets while ensuring reliable operations of power grids. In a context
of growing operational uncertainty, due to increased penetration of renewable
generators and distributed energy resources, operators must continuously
monitor risk in real-time, i.e., they must quickly assess the system’s
behavior under various changes in load and renewable production.
Unfortunately, systematically solving an optimization problem for each such
scenario is not practical given the tight constraints of real-time operations.
To overcome this limitation, this paper proposes to learn an optimization
proxy for SCED, i.e., a Machine Learning (ML) model that can predict an
optimal solution for SCED in milliseconds. Motivated by a principled analysis
of the market-clearing optimizations of MISO, the paper proposes a novel ML
pipeline that addresses the main challenges of learning SCED solutions, i.e.,
the variability in load, renewable output and production costs, as well as the
combinatorial structure of commitment decisions. A novel Classification-Then-
Regression architecture is also proposed, to further capture the behavior of
SCED solutions. Numerical experiments are reported on the French transmission
system, and demonstrate the approach’s ability to produce, within a time frame
that is compatible with real-time operations, accurate optimization proxies
that produce relative errors below $0.6\%$.

##  I Introduction

The Security-Constrained Economic Dispatch (SCED) is a fundamental
optimization model for Transmission System Operators (TSOs) to clear real-time
energy markets while ensuring reliable operations of power grids [1]. In the
US, TSOs like MISO and PJM execute a SCED every five minutes, which means that
the optimization problem must be solved in an even tighter time frame, i.e.,
well under a minute [2]. Security constraints, which enforce robustness
against the loss of any individual component, render SCED models particularly
challenging for large systems [3, 4, 5] unless only a subset of contingencies
is considered. With more distributed resources and increased operational
uncertainty, such computational bottlenecks will only become more critical
[2].

This paper is motivated by the growing share of renewable generation,
especially wind, in the MISO system, which calls for risk-aware market-
clearing algorithms. One particular challenge is the desire to perform risk
analysis in real time, by solving a large number of scenarios for load and
renewable production [6]. However, systematically solving many SCED instances
is not practical given the tight constraints of real-time operations. To
overcome this computational challenge, this paper proposes to learn an
optimization proxy for SCED, i.e., a Machine Learning (ML) model that can
predict an optimal solution for SCED, within acceptable numerical tolerances
and in milliseconds.

It is important to emphasize that the present goal is not to replace
optimization-based market-clearing tools. Instead, the proposed optimization
proxy provides operators with an additional tool for risk assessment. This
allows to quickly evaluate how the system would behave under various
scenarios, without the need to run costly optimizations. In particular,
because these predictions are combined into aggregated risk metrics, small
prediction errors on individual instances are acceptable.

###  I-A Related Literature and Challenges

The combination of ML and optimization for power systems has attracted
increased attention in recent years. A first thread [7, 8, 9, 10, 11, 12, 13,
14] uses ML models to predict an optimal solution to the Optimal Power Flow
(OPF) problem. Pan et al. [7] train a Deep Neural Network (DNN) model to
predict solutions to the security-constrained DC-OPFs, and report results on
systems with no more than 300 buses. One common limitation of ML models, also
noted in [7], is that predictions are not guaranteed to satisfy the original
problem’s constraints. To address this limitation, recent efforts [8, 9]
integrate Lagrangian duality into the training of DNNs in order to capture the
physical and operational constraints of AC-OPFs. A similar approach is
followed in Velloso et al. [12] in the context of preventive security-
constrained DC-OPF. Namely, a DNN is trained using Lagrangian duality
techniques, then used as a proxy for the time-consuming master problem in a
column-and-constraint generation algorithm. More recently, Chatzos et al. [11]
embed the Lagrangian duality framework in a two-stage learning that exploits a
regional decomposition of the power network, enabling a more efficient
distributed training.

Another research thread is the integration of ML models within optimization
algorithms, in order to improve runtime performance. For instance, a number of
papers (e.g., [15, 16, 17, 18, 19]) try to identify active constraints in
order to reduce the problem complexity. In [17], the authors investigate
learning feasible solutions to the unit commitment problem. Venzka el al. [20]
use a DNN model to learn the feasible region of a dynamic SCOPF and transform
the resulting DNN into a mixed-inter linear program.

Existing research papers in ML for power systems typically suffer from two
limitations. On the one hand, most papers report numerical results on small
academic test systems, which are one to two orders of magnitude smaller than
real transmission systems. This is especially concerning, as higher-
dimensional data has an adverse impact on convergence and accuracy of machine-
learning algorithms. On the other hand, almost all papers rely on
artificially-generated data whose distribution does not capture the
variability found in actual operations. For instance, only changes in load are
considered, typically without capturing spatio-temporal correlations. Other
sources of uncertainty, such as renewable production, are not considered, nor
is the variability of economic bids. Finally, changes in commitment decisions
throughout the day are rarely addressed, although they introduce non-trivial,
combinatorial, distribution shifts.

###  I-B Contributions and Outline

To address the above limitations, the paper proposes a novel ML pipeline that
is grounded in the structure of real-world market operations. The approach
leverages the TSO’s forward knowledge of 1) commitment decisions and 2) day-
ahead forecasts for load and renewable productions. Indeed, both are available
by the time the day-ahead market has been executed. Furthermore, the paper
presents an in-depth analysis of the real-time market behavior, which informs
a novel ML architecture to better capture the nature of actual operations.

Figure 1: The proposed machine learning pipeline.

The proposed learning pipeline is depicted in Figure 1. First, commitment
decisions and day-ahead forecasts for load and renewable generation are
gathered, following the clearing of the day-ahead market. Second, this
information is used to generate a training dataset of SCED instances, using
classical data augmentation techniques. Third, specialized ML models are
trained, ideally one for each hour of the operating day, thereby alleviating
the combinatorial explosion of commitment decisions: each ML model only needs
to consider one set of commitments and focus on load/renewable scenarios
around the forecasts for that hour. Fourth, throughout the operating day, the
trained models are used in real time to evaluate a large number of scenarios.
The entire data-generation and training procedure is completed in a few hours.

The rest of the paper is organized as follows. Section II describes the
interplay between day-ahead and real-time markets in the MISO system, and
gives an overview of the real-time SCED for MISO. Section III analyzes the
behavior of the real-time market solutions, proposes a combined
Classification-Then-Regression architecture, and presents the overall ML
pipeline. Numerical experiments on a real-life system are reported in Section
IV

##  II Overview of MISO’s Market-Clearing Pipeline

This section describes the interplay between MISO’s day-ahead and real-time
markets; the reader is referred to [21] for a detailed overview of MISO’s
energy and reserve markets.

###  II-A MISO Optimization Pipeline

The day-ahead market consists of two phases, and is executed every day at
10am. First, a day-ahead security-constrained unit commitment (DA-SCUC) is
executed: it outputs the commitment and regulation status of each generator
for every hour of the following day. Then, a day-ahead SCED is executed to
compute day-ahead prices and settle the market. The results of the day-ahead
market clearing are then posted online at approximately 1pm, i.e., there is a
delay of several hours before the commitment decisions take effect. While MISO
may commit additional units during the operating day, through out-of-market
reliability studies, in practice, $99\%$ of commitment decisions are decided
in the day-ahead market [2]. Accordingly, for simplicity and without loss of
generality, this paper assumes that commitment decisions from the DA-SCUC are
not modified afterwards.

Then, throughout the operating day, the real-time market is executed every 5
minutes. This real-time SCED (RT-SCED) adjusts the dispatch of every generator
in response to variations in load and renewable production, and maximizes
economic benefit. Despite its short-time horizon, the RT-SCED must still
account for uncertainty in load and renewable production. In the MISO system,
this uncertainty mainly stems from the intermittency of wind farms, and from
increasing variability in load. The latter is caused, in part, by the growing
number of behind-the-meter distributed energy resources (DERs) such as
residential storage and rooftop solar. Therefore, operators continuously
monitor the state of the power grid, and may take preventive and/or corrective
actions to ensure safe and reliable operations.

###  II-B RT-SCED Formulation

The RT-SCED used by MISO is a DC-based linear programming formulation that co-
optimizes energy and reserves [22]. The reader is referred to [23] for the
full mathematical formulation; only its core elements are described here. The
computation of market-clearing prices is beyond the scope of this paper, and
is therefore not discussed here.

The RT-SCED model comprises, for each generator, one variable for energy
dispatch and up to four categories of reserves: regulating, spinning, online
supplemental, and offline supplemental. In addition, each generator is subject
to ramping constraints and individual limits on energy and reserve dispatch.
Energy production costs are modeled as piece-wise linear convex functions.
Each market participant submits its own production costs and reserve prices
via MISO’s market portal, and may submit different offers for each hour of the
day. For intermittent generators such as solar and wind, the latest forecast
is used in lieu of binding offers.

At the system level, line losses are estimated in real-time from MISO’s state
estimator, and incorporated in the formulation using a loss factor approach as
in [24]. Transmission constraints are modeled using PTDF matrices, where flow
sensitivities with respect to power injections and withdrawals are provided by
an external tool. Reserves are dispatched on individual generators, in order
to meet zonal and market-wide minimum requirements. Additional constraints
ensure that, in a contingency event, reserves may be deployed without tripping
transmission lines. Power balance, reserve requirements, and transmission
limits are soft, i.e., they may be violated, albeit at a reasonably high cost.

In summary, the RT-SCED receives the following inputs: the commitment
decisions from DA-SCUC, the most recent forecast for load and renewable
production, the economic limits and production costs of the generators, the
current state estimation, and the transmission constraints and reserve
requirements. The RT-SCED produces as outputs the active power and reserve
dispatch for each generator.

##  III The Learning Methodology

This section reviews the learning methodology for the RT-SCED. First, Section
III-A presents an analysis of the behavior of optimal SCED solutions in MISO’s
optimization pipeline. These patterns motivate a novel ML architecture, which
is described in Section III-B, followed by an overview of the proposed ML
pipeline in Section III-C. Further details on the data are given in Section
IV, as wells as in [25].

###  III-A Pattern Analysis of Optimal SCED solutions

Figure 2: Distribution of unique hourly commitments per day, for each month of
the year 2018. Higher values indicate higher variability in commitment
decisions.

The system at hand is the French transmission system, whose network topology
is provided by the French TSO, RTE. It contains $6{,}708$ buses, $8{,}965$
transmission lines and transformers, and $1{,}890$ individual generators,
$1{,}177$ of which are wind and solar generators and $713$ are conventional
units. The MISO optimization pipeline, described in Section II, is replicated
on this system for the entire year 2018, yielding $365$ DA-SCUC instances and
about $100{,}000$ RT-SCED instances. Relevant statistics are reported next.

First, unsurprisingly, commitment decisions display a high variability; for
example, in 2018, a total of $5{,}380$ different hourly commitments were
recorded across the total $8{,}760$ hours of the year, where each “hourly
commitment” is a binary vector of size 713 that contains the (hourly)
commitment status of conventional generators. The intra-day variability of
commitment decisions follows a seasonal pattern, which is illustrated in
Figure 2. Namely, for each month of the year, Figure 2 displays the
distribution, across every day of the month, of the number of unique hourly
commitments over a day. The higher values correspond to the higher variability
in commitment decisions, while the lower values indicate that the commitment
decisions are stable throughout the day. Typically, the variability of
commitment decisions is lower in summer and higher in winter; this behavior is
expected since more generators are online in winter, which naturally tends to
yield more diverse commitments. Indeed, in June 2018, all the days have at
least 5 and at most 19 different hourly commitments, with $50\%$ of days
having less than 7 and $50\%$ having more than 8. In contrast, in January
2018, except for two outliers, every day has at least 19 different
commitments, and 16 days had a different commitment decision every hour.
Overall this combinatorial explosion of commitment decisions has an adverse
effect on ML models, since it creates distribution shifts on unseen
commitments that is detrimental to the performance.

Figure 3: Proportion of generators at their maximum and minimum limits in
January (left) and August (right).

Second, an analysis of RT-SCED solutions reveals, also unsurprisingly, that a
majority of generators are dispatched to either their minimum or maximum
limit; such limits include economic offer data and ramping constraints.
Detailed statistics are reported in Figure 3 for January and August 2018: each
plot quantifies, across all RT-SCED instances for that month, the proportion
of generators dispatched at their minimum (at_min) or maximum (at_max) limit,
or neither (non_tight); these statistics exclude the renewable generators and
the generators for which the minimum and maximum limit are identical. In
January, the median proportion of generators being dispatched at their minimum
(resp., maximum) limit is close to $40\%$ (resp., $20\%$), while in August,
these values are around $45\%$ and $20\%$. The variability is also visibly
higher in winter, echoing the previous observations for commitment decisions.

###  III-B First Classify, then Perform Regression

Figure 4: The Proposed CTR model architecture. Left: Overall structure of the
Classification-Then-Regression (CTR) model where the outputs of the classifier
are used to inform the subsequent regressor. Thus the regressor only needs to
predict the dispatches of non-tight generators; Right: The Latent Surgical
Intervention (LSI) block for each classifier and regressor in the CTR model in
which a gate operator takes the binary indicator to filter out some
representation from fully connected layers. The result is then added into the
representations element-wisely.

The previous analysis indicates that, given the knowledge of which generators
are dispatched to their minimum or maximum limit, only a small number of
generators need to be predicted. This suggests a Classification-Then-
Regression (CTR) approach, first screening the active generators to identify
those at their bounds, and then applying a regression to predict the dispatch
of the remaining generators. The overall architecture is depicted in Figure 4.
The input of the learning task is a dataset
$\mathcal{D}=\\{(\mathbf{x}_{i},\mathbf{y}_{i},\mathbf{p}_{i})\\}_{i=1}^{N}$,
where $\mathbf{x}_{i}$, $\mathbf{y}_{i}$, $\mathbf{p}_{i}$ represent the
$i^{th}$ observation of the system state, the status indicators of the
generators (i.e., whether each generator is at its maximum/minimum limits),
and the optimal dispatches, respectively.

TABLE I: Input features of DNN model Feature | Size | Source
---|---|---
Loads | $L$ | Load forecasts
Cost of generators | $G$ | Bids
Cost of reserves | $2G$ | Bids
Previous solution | $G$ | SCED
Commitment decisions | $G$ | SCUC
Reserve Commitments | $G$ | SCUC
Generator min/max limits | $2G$ | Renewable forecasts
Line losses factor | $2B+1$ | System

The CTR model is composed of two modules; classifier and regressor. First, the
classifier component aims at classifying whether each generator is at its
minimum or maximum limit. It considers the status of the power network with
$B$ buses, $G$ generators, and $L$ loads as its inputs. Specifically, the
classifier, parameterized by $\mathbf{w}_{1}$, is a mapping
$f_{\mathbf{w}_{1}}:\mathbb{R}^{d}\xrightarrow{}\\{0,1\\}^{2G}$, where $d$ is
the dimension of the input features. The overall input features describe the
state of the power system, detailed in Table I. Also, it outputs the binary
vector of $2G$ meaning that, for each generator, there are two classification
choices: one for determining whether it is dispatched at its maximum limit and
one for determining whether it is dispatched at its minimum limit. In the
experiments presented in the subsequent section, the dimension $d$ of the
input space can be as large as $24{,}403$ in the RTE system. The optimal
trainable parameters $\mathbf{w^{*}_{1}}$ in the classifier are obtained
through minimizing the loss function as follows:

$\displaystyle\mathbf{w^{*}_{1}}=\operatorname*{arg\,min}_{\mathbf{w_{1}}}\frac{1}{N}\sum_{i=1}^{N}\mathcal{L}_{c}(\mathbf{y}_{i},f_{\mathbf{w}_{1}}(\mathbf{x}_{i})),$
(1)

where $\mathcal{L}_{c}$ denotes the cross entropy loss, i.e.,

$\displaystyle\mathcal{L}_{c}(\mathbf{y}_{i},\hat{\mathbf{y}}_{i})=-\sum_{j=1}^{2G}\mathbf{y}_{i,j}\log(\hat{\mathbf{y}}_{i,j})+(1-\mathbf{y}_{i,j})\log(1-\hat{\mathbf{y}}_{i,j}).$
(2)

The second architectural component, the regressor that is parameterized by
$\mathbf{w}_{2}$, is a mapping
$f_{\mathbf{w}_{2}}:\mathbb{R}^{d+2G}\xrightarrow{}\mathbb{R}^{G}$. The
additional $2G$ features in the input of the regressor come from the outputs
of the classifier. Given the trained classifier $f_{\mathbf{w}^{*}_{1}}$, the
optimal trainable parameters of the regressor $\mathbf{w^{*}_{2}}$ are
obtained by minimizing the loss function $\mathcal{L}_{r}$ over all training
instances as

$\displaystyle\mathbf{w^{*}_{2}}=\operatorname*{arg\,min}_{\mathbf{w_{2}}}\frac{1}{N}\sum_{i=1}^{N}\mathcal{L}_{r}\left(\mathbf{p}_{i},f_{\mathbf{w}_{2}}\left(\mathbf{x}_{i},f_{\mathbf{w}^{*}_{1}}(\mathbf{x}_{i})\right)\right),$
(3)

where $\mathcal{L}_{r}$ is the mean absolute error (MAE) loss, i.e.,
$\mathcal{L}_{r}(\mathbf{p},\hat{\mathbf{p}})=\|\mathbf{p}-\hat{\mathbf{p}}\|_{1}$.

The CTR architecture features a deep neural network (DNN). Specifically, it
uses a Latent Surgical Intervention (LSI) network [26] as its building block.
LSI is a variant of residual neural networks [27] that is augmented by binary
interventions. As illustrated in Figure 4, the LSI block exploits, via gate
operators, the binary information coming from commitment decisions and the
classifier. As mentioned earlier, because the variability in SCED solutions
primarily comes from the combinatorial nature of the commitment decisions, it
is crucial to design the DNN architecture to use commitment decisions as the
input of the model. Thus, the LSI-based CTR model makes it possible to learn
the generator dispatch from various commitment decisions, thereby allowing the
proposed approach to generalize to the cases where the models are trained over
multiple commitment decisions.

###  III-C Machine Learning Pipeline

The high variability of SCED instances in real operations, in terms of
commitment decisions and forecasts of loads and renewable energy sources,
makes it extremely challenging to learn the SCED optimization. To mitigate
this significant variability, this paper proposes a learning pipeline that
closely follows MISO’s optimization pipeline. The machine-learning pipeline is
depicted in Figure 1 and consists of three phases: data preparation, training,
and prediction.

#### Data Preparation

At 12pm of day $D-1$, the TSO outputs the commitment decisions of generators
for the next 24 hours and Monte-Carlo scenarios for day $D$. Each scenario
consists of 15 minute-level forecasts for load demands and renewable
generators. Since the SCED is solved every 5 minutes throughout a day, the ML
pipeline first linearly interpolates the forecasts at 5 minute level. Further
data augmentation may be performed, if necessary, by perturbing load and
renewable generations following the strategy described in [9]. The data is
then used as input to SCED optimization models, which generate the various
optimal dispatches. The entire dataset is then divided into subsets, each of
which spans one or a few hours, based on computational considerations. Indeed,
the goal is to strike the proper balance between the accuracy of the models
and the training time, since the models have to be available before midnight.
This data preparation process yields the input instances to the learning task
described previously.

#### Training

First, each subset is split into the traditional training/validation/test
instances. The CTR models are trained on training instances in sequence while
the validation dataset is used for hyperparameter tuning and the test dataset
is used for reporting the performance of the machine learning models. This
training step takes place in parallel for each subset.

#### Prediction

Starting from midnight on day $D$, at each time step, the corresponding ML
model takes the latest system state as input and predicts the optimal dispatch
of SCED models in real-time.

##  IV Experimental Results

###  IV-A Test Cases

The experiments replicate MISO’s operations on the French transmission system,
for four representative days in 2018, namely, February 12th, April 5th, August
26th, and October 23rd. These four days were selected at random to represent
annual seasonality. For each operating day, 2000 Monte-Carlo scenarios for
total load and renewable power production are sampled in a day-ahead fashion,
i.e., scenarios are produced at noon of the previous day. These scenarios are
illustrated in Figure 5, which displays 200 scenarios for solar production,
wind production, total load, and total net load, respectively. The total net
load in this figure represents the amount of power production expected to be
generated by the conventional generators, which is identical to the total load
minus the renewable generation.

The forecasting models for load consumption and renewable power generation in
this paper used Long Short-Term Memory (LSTM) neural networks [28] and the
scenarios are generated by a MC-dropout approach [29]. Other forecasting
methods and scenario generation approaches could be used instead, without
changing the methodology.

Figure 5: Day-ahead scenarios for solar output (top-left), wind output (top-
right), total load (bottom-left) and net load (bottom-right). 200 Monte-Carlo
scenarios are depicted.

Day-ahead commitment decisions are obtained by solving a DA-SCUC problem and
recording its solution. The SCUC formulation used in the present experiment
follows MISO’s DA-SCUC described in [21, 30]. Again, the proposed ML pipeline
is agnostic to the SCUC formulation itself and how it is solved, as it only
requires the resulting commitments.

Finally, for each scenario and each day, 288 RT-SCED instances are solved (one
every 5 minutes), yielding a total of $576{,}000$ instances. This initial
dataset is then divided into 24 hourly datasets, each containing $24{,}000$
SCED instances, since commitment decisions are hourly. To avoid information
leakage, the data is split between training/validation/testing instances as
follows: $85\%$ of scenarios are used for training, $7.5\%$ for validation,
and $7.5\%$ for testing. All the reported results are in terms of the testing
instances.

The optimization problems (SCUC and SCED) are formulated in the JuMP modeling
language [31] and solved with Gurobi 9.1 [32], using Linux machines with dual
Intel Xeon<EMAIL_ADDRESS>CPUs on the PACE Phoenix cluster [33]; the entire data
generation phase is performed in less than 2 hours. The proposed CTR model is
implemented using PyTorch [34] and trained using the Adam optimizer [35] with
a learning rate of 5e-4. A grid search is performed to choose the
hyperparameters, i.e., the dimension of the hidden layers (taken in the set
{64, 128, 256}), and the number of layers (from the set {2, 3, 4, 5, 6}) for
the fully connected layers in the LSI block. A leaky ReLU with leakiness of
$\alpha=0.01$ is used as the activation function in the LSI block. An early
stopping criterion with 20 epochs is used for training the classifier and
regressor models in order to prevent overfitting: when the loss values (Eq. 1
and Eq. 3) on the validation dataset do not decrease for 20 consecutive
epochs, the training process is terminated. Training is performed using Tesla
V100-PCIE GPUs with 16GBs HBM2 RAM, on machines with Intel CPU cores at
2.1GHz.

###  IV-B Baselines

The proposed CTR models are evaluated against the following baselines: Naive
CTR, Naive Regression (Naive Reg), and Regression (Reg). The naive baselines
replicate the behavior of the previous SCED solution, i.e., they use the
dispatch solution obtained 5 minutes earlier. The naive baselines are
motivated by the fact that, if the system does not change much between two
consecutive intervals, then the SCED solution should not be changed much
either. The Naive CTR uses the same approach as the CTR models: it first
predicts whether a given generator is at its minimum (resp., maximum) limit if
it was at its minimum (resp. maximum) limit in the previous dispatch; then,
for the remaining generators, it predicts the active dispatch using the
regressor. The naive baselines are expected to perform worse when large
fluctuations in load and renewable production are observed, which typically
occurs in the morning and evening. Note that these times of the day also
display the largest ramping needs, and are among the most critical for
reliability. The effectiveness of the CTR architecture is also demonstrated by
comparing it to Reg for the optimal active dispatch, i.e., by omitting the
classification step from the CTR.

###  IV-C Optimal Dispatch Prediction Errors

TABLE II: Average Classification Accuracy (%) of the CTR Classifier and Naive CTR on 4 Representative Days in 2018. | Dates |
---|---|---
Methods | Feb. 12 | Apr. 5 | Aug. 26 | Oct. 23 | Avg.
Naive classifier | 98.26 | 97.79 | 98.64 | 98.31 | 98.25
CTR | 99.56 | 99.18 | 99.40 | 99.24 | 99.35

Table II reports the overall classification accuracy of the CTR classifier and
the Naive CTR across four representative days in 2018. Surprisingly, the naive
classifier is a strong baseline, with an accuracy that ranges from $97.79\%$
in the Spring to $98.64\%$ in the Summer. The CTR classifier always improves
on the baseline, by around one percentage point on average and by up to $1.40$
percentage point in the Spring. More detailed results about the classifiers
are given in Appendix -A.

Figure 6: Mean Absolute Error (MAE) Over Time on Feb. 12 (top-left), Apr. 5
(top-right), Aug. 26 (bottom-left), and Oct. 23 (bottom-right). TABLE III:
Mean Absolute Error (MW) by Generator Size†.

Date | Method | Small | Medium | Large | All
---|---|---|---|---|---
Feb. 12 | Naive Reg | 0.122 | 0.465 | 1.602 | 0.374
Naive CTR | 0.084 | 0.333 | 1.128 | 0.262
Reg | 0.057 | 0.188 | 0.654 | 0.153
CTR | 0.043 | 0.141 | 0.535 | 0.117
Apr. 05 | Naive Reg | 0.242 | 0.345 | 7.480 | 0.772
Naive CTR | 0.197 | 0.220 | 4.374 | 0.463
Reg | 0.149 | 0.152 | 2.553 | 0.282
CTR | 0.105 | 0.110 | 2.291 | 0.241
Aug. 26 | Naive Reg | 0.097 | 0.256 | 7.447 | 0.637
Naive CTR | 0.080 | 0.149 | 4.045 | 0.352
Reg | 0.054 | 0.124 | 2.454 | 0.218
CTR | 0.034 | 0.064 | 2.220 | 0.190
Oct. 23 | Naive Reg | 0.176 | 0.535 | 8.425 | 0.778
Naive CTR | 0.140 | 0.341 | 4.596 | 0.434
Reg | 0.106 | 0.192 | 2.724 | 0.263
CTR | 0.076 | 0.145 | 2.525 | 0.235

†Small: 0-10MW; Medium: 10-100MW; Large: $>$100MW

Figure 6 reports the mean absolute error (MAE) for active power dispatch of
the different models for each of the considered days. Given a ground-truth
dispatch $p_{i}^{g}$ and the predicted dispatch $\hat{p}^{g}_{i}$, the MAE is
defined as

$\displaystyle
MAE=\frac{1}{N}\frac{1}{G}\sum_{i=1}^{N}\sum_{g=1}^{G}\|p^{g}_{i}-\hat{p}^{g}_{i}\|,$

where $N$ is the number of instances in the test dataset and $G$ is the number
of generators for each instance. As shown in Figure 6, the MAEs of CTR are
always lower than those of Reg, demonstrating the benefits of the classifier.
The MAEs of Naive CTR are always lower than those of Naive Reg, showing that
even a naive classifier has benefits. Moreover, the methods using DNNs for
classification and regression, i.e., CTR and Reg, are always better than their
naive counterparts, demonstrating the value of deep learning. Note that the
performance of the naive methods fluctuates during the day, mainly due to the
variability of the commitments. The naive methods are not robust with respect
to changes in commitments, contrary to the proposed CTR approach.

To investigate how the machine-learning models perform for different generator
types, Table III reports their behavior for different generator sizes. The
generators are clustered into three groups based on their actual active
dispatches, and the MAE is reported for each group separately. Small-size
generators have a capacity between $0$ and $10$MW, medium-size generators have
a capacity between $10$ and $100$MW, and large-size generators have a capacity
above $100$MW. The last column report the MAEs across all generators. Table
III further confirms that the CTR models consistently outperform the
corresponding regressors. These improvements are most notable for small and
medium generators which are expected to have higher variability: the MAEs are
decreased by up to $37.04\%$ across small generators and $48.39\%$ across
medium generators (in Aug. 26). Moreover, across four days, the CTR model
achieves a Mean Average Percentage Error (MAPE) of $0.59\%$ and $0.34\%$ for
medium and large generators.

###  IV-D Solving Time vs Inference Time

Solving the SCED using traditional optimization tools takes, on average,
$15.93$ seconds. Actual computing times fluctuate throughout the day, with
increased load and congestion leading to higher computing times. In contrast,
evaluating the optimization proxy for a batch of 288 instances takes an
average $1.5$ milliseconds. In other words, roughly $200{,}000$ scenarios may
be evaluated in less than one second. This represents an improvement of 4
orders of magnitude, even under the assumption that several hundred SCED
instances can be solved in parallel.

##  V Sensitivity Analysis

###  V-A Motivations and Experiment Settings

How many CTR models need to be trained is an important question to ponder in
practice. Preparing a single CTR model for 24 hours on Day $D$ would be most
convenient. However, as described in Section III-A, the variability in
commitments is notoriously harmful to the prediction accuracy and training
time of the CTR model. However, successive hours during a day may only have a
few differences in commitments and hence a single CTR model may be sufficient
for predicting the associated SCED optimizations.

To answer that question, the original instance data is used to produce a
variety of datasets, containing data for consecutive 2, 3, 4, 6, 8, 12, and 24
hours. For instance, the 8 hours dataset contains three sets of instances,
each grouping SCED data for 8 successive hours. Three models are then trained
using instances covering their 8 hours of data. All models use the termination
criterion presented in Section IV. The various CTR models for 1, 2, …, 24
hours are then compared for the 24 hours of the day, using the models
appropriate for each hour. It is expected that the quality of the prediction
will deteriorate with a coarser granularity, since there will be a larger
variability in commitments and net loads. However, the LSI architecture of the
CTR model and the availability of more instances during training may
compensate for this increase in variability. In the following, CTR models
trained on $i$ hours are denoted as CTRi.

###  V-B Results

Figure 7: MAEs of CTR models for representative days: the CTR model scales
well when trained over multiple hours. The performance degradation becomes
significant only when aggregating 6 hours or more.

Figure 7 illustrates the MAE values of the various CTR models. The performance
of the CTR model remains strong even when aggregating up to 4 successive
hours. As more hours are aggregated, the performance starts to degrade. To get
more insight on the performance of the various CTR models, it is useful to
consider the energy distance [36] between the empirical distributions of the
testing and training instances. Recall that, given two empirical distributions
$\\{x_{i}\\}_{i=1}^{N}$ and $\\{y_{i}\\}_{i=1}^{M}$, their energy distance is
given by

$\displaystyle\mathcal{E}(\mathcal{X},\mathcal{Y})$
$\displaystyle=\frac{2}{NM}\sum_{i=1}^{N}\sum_{j=1}^{M}\|x_{i}-y_{j}\|$
$\displaystyle-\frac{1}{N^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\|x_{i}-x_{j}\|$
$\displaystyle-\frac{1}{M^{2}}\sum_{i=1}^{M}\sum_{j=1}^{M}\|y_{i}-y_{j}\|.$

Figure 8: Scaled MAE value of the CTR models vs. energy distance between the
test datasets and the training datasets for given hour’s period on Feb. 12
(top-left), Apr. 5 (top-right), Aug. 26 (bottom-left), and Oct. 23 (bottom-
right).

It is thus interesting to study the relationship between the MAE of a CTR
model and the energy distance between its testing and training datasets.
Figure 8 reports the relationship between the scaled MAE and the energy
distance, where the scaled MAE is the MAE value divided by the best MAE value
across all experiments. Each dot captures the scaled MAE value (y-axis) for a
CTR model and the associated energy distance (x-axis) between its training and
testing datasets.

Observe first that the energy distance increases as more hourly data are
combined into the training set. In Figure 8, the dots from CTRi tend to have
smaller energy distances than those of CTRj for $j>i$. The scaled MAE also
increases as the energy distance grows. More interestingly, Figure 8 shows
that the MAE value only increases slightly as the energy distance becomes
larger. Specifically, an increase of $10^{2}$ in energy distance produces an
increase of only 10% in MAE (corresponding to a scaled MAE of 1.1). These
results confirm that the CTR models scale well when facing a reasonable amount
of variability in commitment decisions and net loads. However, when the energy
distance goes over a certain threshold, then the performance of CTR models
starts to degrade significantly even in presence of more data and more
training time. Note that the energy distance between the training and testing
datasets can be computed before the training process. Hence, the proper
granularity for the CTR model can be chosen to obtain the desired accuracy.
This quantification shows the potential of CTR to perform well with more
complex components involved in the MISO pipeline such as the Look Ahead
Commitment (LAC).

Figure 9: Training Times of CTRi Models: Regression (Left) and Classification
(Right)

Figure 9 reports the training times of the CTR models. The average training
time increases as more hourly data are included. The regression takes more
training time than the classification task. Section III-C indicated that the
training process must complete within about 12 hours after accumulating the
data instances at 12pm on day $D-1$. This is clearly possible for all the CTRi
models with $i\leq 4$. Consider CTR4 for instance. Figure 9 shows that it can
be trained for a 4-hour block within 10 hours in the worst case (4 hours for
classification and 6 hours for regression). As a result, all 6 CRT4 models can
easily be trained in parallel within the required time frame.

##  VI Conclusion

The paper presented several challenges that may hamper the use of ML models in
power systems operations, including the variability of electricity demand and
renewable production, the variations in production costs, and the
combinatorial structure of commitment decisions. To address these challenges,
the paper proposed a novel ML pipeline that leverages day-ahead forecasts and
the TSO’s knowledge of commitment decisions several hours before they take
place. The proposed pipeline consists of two main phases: 1) day-ahead data
preparation and training; and 2) real-time predictions, which fits naturally
within the operations of a TSO. Moreover, informed by the behavior of real-
time markets, the paper proposed a novel Classification-Then-Regression (CTR)
approach that leverages deep neural networks based on Latent Surgical
Interventions to capture generator commitments and generators operating at
their limits. Computational experiments that replicated MISO’s operational
pipeline on a real, large-scale transmission system demonstrated the
feasibility of the approach. In particular, the results show that optimization
proxies based on ML models have the potential to provide operators with new
tools for real-time risk monitoring.

Several extensions are possible, which will be the topic of future work. The
integration, during training, of Lagrangian-based penalties has the potential
to further improve the performance of the neural network models. Such methods
have been shown to improve the feasibility of the predictions with respect to
the original constraints. Moreover, once trained, the proposed classifier may
be used to also accelerate the SCED resolution, by providing hints as to which
variables should be fixed at their minimum or maximum limit. Finally, the
proposed CTR model is not limited to SCED models, and may be applied to any
optimization problems.

## Acknowledgments

This research is partly funded by NSF Awards 1912244 and 2112533, and ARPA-E
Perform Award AR0001136.

## References

  * [1] A. J. Conejo and L. Baringo, _Power system operations_. Springer, 2018.
  * [2] Y. Chen, F. Wang, J. Wan, and F. Pan, “Developing next generation electricity market clearing optimization software,” in _2018 IEEE Power Energy Society General Meeting (PESGM)_ , 2018, pp. 1–5.
  * [3] N. Chiang and A. Grothey, “Solving security constrained optimal power flow problems by a structure exploiting interior point method,” _Optimization and Engineering_ , vol. 16, no. 1, pp. 49–71, 2015.
  * [4] Q. Wang, J. D. McCalley, T. Zheng, and E. Litvinov, “Solving corrective risk-based security-constrained optimal power flow with Lagrangian relaxation and Benders decomposition,” _International Journal of Electrical Power & Energy Systems_, vol. 75, pp. 255–264, 2016.
  * [5] A. Velloso, P. Van Hentenryck, and E. S. Johnson, “An exact and scalable problem decomposition for security-constrained optimal power flow,” _Electric Power Systems Research_ , vol. 195, p. 106677, 2021.
  * [6] T. Werho, J. Zhang, V. Vittal, Y. Chen, A. Thatte, and L. Zhao, “Scenario generation of wind farm power for real-time system operation,” 2021.
  * [7] X. Pan, T. Zhao, and M. Chen, “Deepopf: Deep neural network for dc optimal power flow,” in _2019 IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm)_. IEEE, 2019, pp. 1–6.
  * [8] F. Fioretto, T. W. Mak, and P. Van Hentenryck, “Predicting ac optimal power flows: Combining deep learning and lagrangian dual methods,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, no. 01, 2020, pp. 630–637.
  * [9] M. Chatzos, F. Fioretto, T. W. Mak, and P. Van Hentenryck, “High-fidelity machine learning approximations of large-scale optimal power flow,” _arXiv preprint arXiv:2006.16356_ , 2020.
  * [10] X. Lei, Z. Yang, J. Yu, J. Zhao, Q. Gao, and H. Yu, “Data-driven optimal power flow: A physics-informed machine learning approach,” _IEEE Transactions on Power Systems_ , vol. 36, no. 1, pp. 346–354, 2020.
  * [11] M. Chatzos, T. W. Mak, and P. Van Hentenryck, “Spatial network decomposition for fast and scalable ac-opf learning,” _arXiv preprint arXiv:2101.06768_ , 2021.
  * [12] A. Velloso and P. Van Hentenryck, “Combining deep learning and optimization for security-constrained optimal power flow,” _arXiv preprint arXiv:2007.07002_ , 2020.
  * [13] A. S. Zamzam and K. Baker, “Learning optimal solutions for extremely fast ac optimal power flow,” in _2020 IEEE International Conference on Communications, Control, and Computing Technologies for Smart Grids (SmartGridComm)_. IEEE, 2020, pp. 1–6.
  * [14] D. Owerko, F. Gama, and A. Ribeiro, “Optimal power flow using graph neural networks,” in _ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2020, pp. 5930–5934.
  * [15] D. Deka and S. Misra, “Learning for dc-opf: Classifying active sets using neural nets,” in _2019 IEEE Milan PowerTech_. IEEE, 2019, pp. 1–6.
  * [16] S. Misra, L. Roald, and Y. Ng, “Learning for constrained optimization: Identifying optimal active constraint sets,” _INFORMS Journal on Computing_ , 2021.
  * [17] Á. S. Xavier, F. Qiu, and S. Ahmed, “Learning to solve large-scale security-constrained unit commitment problems,” _INFORMS Journal on Computing_ , vol. 33, no. 2, pp. 739–756, 2021.
  * [18] Y. Yang, Z. Yang, J. Yu, K. Xie, and L. Jin, “Fast economic dispatch in smart grids using deep learning: An active constraint screening approach,” _IEEE Internet of Things Journal_ , vol. 7, no. 11, pp. 11 030–11 040, 2020\.
  * [19] N. Guha, Z. Wang, M. Wytock, and A. Majumdar, “Machine learning for ac optimal power flow,” _arXiv preprint arXiv:1910.08842_ , 2019.
  * [20] A. Venzke, D. T. Viola, J. Mermet-Guyennet, G. S. Misyris, and S. Chatzivasileiadis, “Neural networks for encoding dynamic security-constrained optimal power flow to mixed-integer linear programs,” _arXiv preprint arXiv:2003.07939_ , 2020.
  * [21] MISO, “Energy and operating reserve markets,” Business Practices Manual Energy and Operating Reserve Markets.
  * [22] X. Ma, Y. Chen, and J. Wan, “Midwest iso co-optimization based real-time dispatch and pricing of energy and ancillary services,” in _2009 IEEE Power Energy Society General Meeting_ , 2009, pp. 1–6.
  * [23] MISO, “Real-time energy and operating reserve market software formulations and business logic,” Business Practices Manual Energy and Operating Reserve Markets Attachment D.
  * [24] B. Eldridge, R. P. O’Neill, and A. Castillo, “Marginal loss calculations for the dcopf,” 2017.
  * [25] M. Chatzos, M. Tanneau, and P. Van Hentenryck, “Data-driven Time Series Reconstruction for Modern Power Systems Research,” submitted.
  * [26] B. Donnot, I. Guyon, Z. Liu, M. Schoenauer, A. Marot, and P. Panciatici, “Latent surgical interventions in residual neural networks,” 2018.
  * [27] K. He, X. Zhang, S. Ren, and J. Sun, “Identity mappings in deep residual networks,” in _European conference on computer vision_. Springer, 2016, pp. 630–645.
  * [28] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [29] Y. Gal and Z. Ghahramani, “Dropout as a bayesian approximation: Representing model uncertainty in deep learning,” in _international conference on machine learning_. PMLR, 2016, pp. 1050–1059.
  * [30] MISO, “Day-ahead energy and operating reserve market software formulations and business logic,” Business Practices Manual Energy and Operating Reserve Markets Attachment B.
  * [31] I. Dunning, J. Huchette, and M. Lubin, “Jump: A modeling language for mathematical optimization,” _SIAM Review_ , vol. 59, no. 2, pp. 295–320, 2017.
  * [32] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2021. [Online]. Available: https://www.gurobi.com
  * [33] PACE, _Partnership for an Advanced Computing Environment (PACE)_ , 2017. [Online]. Available: http://www.pace.gatech.edu
  * [34] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “Pytorch: An imperative style, high-performance deep learning library,” in _Advances in Neural Information Processing Systems 32_ , H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, Eds. Curran Associates, Inc., 2019, pp. 8024–8035.
  * [35] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [36] M. L. Rizzo and G. J. Székely, “Energy distance,” _wiley interdisciplinary reviews: Computational statistics_ , vol. 8, no. 1, pp. 27–38, 2016.

###  -A Classification accuracy across four days

The detailed results of Naive and CTR model trained under different hourly
data in four days are shown in Table IV V VI VII, respectively. First, across
all days, CTR models always outperform the baseline model. Moreover, the
performance decay is up to $3.74\%$ for Naive model (in Apr. 5) and up to
$1.71\%$ for CTR model (CTR12 in Aug. 26) along the day. CTR models have much
smaller performance variations with the fluctuation of loads and renewables
along the four days. Second, the CTR model scales well with more hours data
included. The stark performance decay can only be observed with more than 12
hours data included.

TABLE IV: Classification accuracy of Naive and CTR models (Feb. 12)

Hour | Naive | CTR1 | CTR2 | CTR3 | CTR4 | CTR6 | CTR8 | CTR12 | CTR24
---|---|---|---|---|---|---|---|---|---
03:00 - 04:00 | 98.37 | 99.24 | 99.24 | 99.29 | 99.24 | 99.31 | 99.26 | 99.17 | 98.94
04:00 - 05:00 | 98.79 | 99.36 | 99.38 | 99.40 | 99.32 | 99.42 | 99.37 | 99.30 | 99.11
05:00 - 06:00 | 98.74 | 99.68 | 99.70 | 99.71 | 99.69 | 99.71 | 99.70 | 99.65 | 99.55
06:00 - 07:00 | 98.60 | 99.85 | 99.87 | 99.86 | 99.87 | 99.86 | 99.86 | 99.84 | 99.77
07:00 - 08:00 | 98.14 | 99.88 | 99.88 | 99.87 | 99.85 | 99.87 | 99.85 | 99.85 | 99.80
08:00 - 09:00 | 98.66 | 99.97 | 99.97 | 99.97 | 99.97 | 99.97 | 99.96 | 99.96 | 99.95
09:00 - 10:00 | 99.27 | 99.96 | 99.96 | 99.96 | 99.97 | 99.96 | 99.96 | 99.95 | 99.94
10:00 - 11:00 | 99.04 | 99.80 | 99.80 | 99.80 | 99.79 | 99.77 | 99.77 | 99.74 | 99.70
11:00 - 12:00 | 97.64 | 99.24 | 99.38 | 99.38 | 99.34 | 99.29 | 99.31 | 99.21 | 99.08
12:00 - 13:00 | 97.96 | 99.37 | 99.36 | 99.34 | 99.31 | 99.30 | 99.31 | 99.19 | 99.04
13:00 - 14:00 | 97.97 | 99.22 | 99.26 | 99.25 | 99.21 | 99.21 | 99.18 | 99.05 | 98.91
14:00 - 15:00 | 97.25 | 99.29 | 99.32 | 99.30 | 99.29 | 99.29 | 99.27 | 99.13 | 98.95
15:00 - 16:00 | 98.37 | 99.25 | 99.23 | 99.25 | 99.18 | 99.19 | 99.16 | 99.06 | 98.94
16:00 - 17:00 | 98.71 | 99.52 | 99.53 | 99.54 | 99.53 | 99.51 | 99.45 | 99.43 | 99.34
17:00 - 18:00 | 98.12 | 99.69 | 99.69 | 99.68 | 99.69 | 99.65 | 99.62 | 99.59 | 99.53
18:00 - 19:00 | 98.93 | 99.88 | 99.88 | 99.88 | 99.88 | 99.82 | 99.83 | 99.82 | 99.78
19:00 - 20:00 | 99.17 | 99.87 | 99.88 | 99.88 | 99.88 | 99.84 | 99.85 | 99.83 | 99.81
20:00 - 21:00 | 99.03 | 99.68 | 99.64 | 99.69 | 99.59 | 99.58 | 99.59 | 99.58 | 99.54
21:00 - 22:00 | 97.83 | 99.47 | 99.47 | 99.39 | 99.38 | 99.31 | 99.38 | 99.33 | 99.22
22:00 - 23:00 | 96.11 | 99.45 | 99.47 | 99.42 | 99.37 | 99.32 | 99.36 | 99.30 | 99.23
23:00 - 24:00 | 96.82 | 99.02 | 98.96 | 98.89 | 98.82 | 98.76 | 98.81 | 98.71 | 98.60

TABLE V: Classification accuracy of Naive and CTR models (Apr. 05)

Hour | Naive | CTR1 | CTR2 | CTR3 | CTR4 | CTR6 | CTR8 | CTR12 | CTR24
---|---|---|---|---|---|---|---|---|---
03:00 - 04:00 | 97.95 | 99.24 | 99.24 | 99.24 | 99.24 | 99.24 | 99.21 | 99.11 | 98.95
04:00 - 05:00 | 98.29 | 99.34 | 99.33 | 99.33 | 99.31 | 99.33 | 99.31 | 99.23 | 99.08
05:00 - 06:00 | 98.43 | 99.44 | 99.45 | 99.46 | 99.44 | 99.46 | 99.44 | 99.36 | 99.19
06:00 - 07:00 | 99.02 | 99.66 | 99.69 | 99.69 | 99.69 | 99.60 | 99.69 | 99.65 | 99.55
07:00 - 08:00 | 98.20 | 99.56 | 99.54 | 99.55 | 99.53 | 99.46 | 99.53 | 99.49 | 99.35
08:00 - 09:00 | 98.22 | 99.31 | 99.34 | 99.30 | 99.32 | 99.24 | 99.23 | 99.24 | 99.11
09:00 - 10:00 | 98.05 | 99.27 | 99.27 | 99.27 | 99.28 | 99.15 | 99.16 | 99.17 | 98.97
10:00 - 11:00 | 98.40 | 99.26 | 99.29 | 99.30 | 99.30 | 99.21 | 99.23 | 99.23 | 99.08
11:00 - 12:00 | 97.68 | 98.94 | 98.97 | 98.98 | 99.00 | 98.86 | 98.89 | 98.88 | 98.64
12:00 - 13:00 | 97.11 | 98.46 | 98.51 | 98.45 | 98.44 | 98.22 | 98.42 | 98.10 | 98.05
13:00 - 14:00 | 97.36 | 98.83 | 98.80 | 98.77 | 98.78 | 98.56 | 98.65 | 98.45 | 98.32
14:00 - 15:00 | 96.64 | 98.94 | 98.97 | 98.98 | 99.00 | 98.82 | 98.89 | 98.71 | 98.64
15:00 - 16:00 | 96.55 | 99.10 | 99.21 | 99.14 | 99.21 | 99.01 | 99.07 | 98.93 | 98.81
16:00 - 17:00 | 97.36 | 98.93 | 98.71 | 98.81 | 98.56 | 98.60 | 98.54 | 98.52 | 98.38
17:00 - 18:00 | 97.03 | 98.72 | 98.59 | 98.67 | 98.51 | 98.46 | 98.47 | 98.42 | 98.31
18:00 - 19:00 | 98.01 | 99.15 | 99.19 | 99.19 | 99.02 | 99.11 | 99.02 | 98.95 | 98.88
19:00 - 20:00 | 98.72 | 99.33 | 99.39 | 99.40 | 99.21 | 99.33 | 99.27 | 99.20 | 99.15
20:00 - 21:00 | 98.87 | 99.37 | 99.41 | 99.42 | 99.38 | 99.38 | 99.31 | 99.24 | 99.21
21:00 - 22:00 | 98.04 | 99.41 | 99.43 | 99.34 | 99.39 | 99.35 | 99.30 | 99.26 | 99.22
22:00 - 23:00 | 95.28 | 99.22 | 99.23 | 99.23 | 99.25 | 99.17 | 99.07 | 99.00 | 98.87
23:00 - 24:00 | 98.37 | 99.39 | 99.36 | 99.35 | 99.37 | 99.33 | 99.25 | 99.22 | 99.16

TABLE VI: Classification accuracy of Naive and CTR models (Aug. 26)

Hour | Naive | CTR1 | CTR2 | CTR3 | CTR4 | CTR6 | CTR8 | CTR12 | CTR24
---|---|---|---|---|---|---|---|---|---
03:00 - 04:00 | 99.73 | 99.95 | 99.94 | 99.92 | 99.95 | 99.92 | 99.91 | 99.90 | 99.88
04:00 - 05:00 | 99.73 | 99.92 | 99.91 | 99.89 | 99.88 | 99.89 | 99.88 | 99.88 | 99.85
05:00 - 06:00 | 99.66 | 99.88 | 99.90 | 99.88 | 99.87 | 99.88 | 99.87 | 99.86 | 99.84
06:00 - 07:00 | 97.93 | 99.65 | 99.67 | 99.66 | 99.67 | 99.62 | 99.65 | 99.61 | 99.50
07:00 - 08:00 | 98.94 | 99.63 | 99.63 | 99.63 | 99.63 | 99.60 | 99.60 | 99.59 | 99.51
08:00 - 09:00 | 98.94 | 99.62 | 99.62 | 99.62 | 99.59 | 99.58 | 99.46 | 99.58 | 99.49
09:00 - 10:00 | 98.87 | 99.52 | 99.55 | 99.49 | 99.52 | 99.51 | 99.39 | 99.52 | 99.41
10:00 - 11:00 | 98.91 | 99.59 | 99.58 | 99.57 | 99.58 | 99.57 | 99.50 | 99.58 | 99.50
11:00 - 12:00 | 97.98 | 99.15 | 99.18 | 99.19 | 99.19 | 99.17 | 99.15 | 99.15 | 99.10
12:00 - 13:00 | 98.08 | 99.03 | 99.07 | 99.06 | 99.07 | 98.95 | 98.83 | 98.44 | 98.82
13:00 - 14:00 | 98.31 | 98.71 | 98.87 | 98.97 | 98.98 | 98.81 | 98.64 | 98.35 | 98.65
14:00 - 15:00 | 98.05 | 98.68 | 98.84 | 98.93 | 98.97 | 98.90 | 98.56 | 98.20 | 98.60
15:00 - 16:00 | 98.15 | 98.77 | 98.90 | 98.89 | 98.98 | 98.90 | 98.67 | 98.29 | 98.68
16:00 - 17:00 | 97.85 | 98.75 | 98.75 | 98.88 | 98.62 | 98.94 | 98.32 | 98.19 | 98.66
17:00 - 18:00 | 98.96 | 99.32 | 99.36 | 99.38 | 99.38 | 99.31 | 99.37 | 99.11 | 99.31
18:00 - 19:00 | 98.77 | 99.42 | 99.46 | 99.44 | 99.43 | 99.36 | 99.39 | 99.10 | 99.28
19:00 - 20:00 | 98.68 | 99.64 | 99.64 | 99.63 | 99.61 | 99.59 | 99.58 | 99.42 | 99.53
20:00 - 21:00 | 98.42 | 99.55 | 99.50 | 99.50 | 99.53 | 99.46 | 99.46 | 99.25 | 99.38
21:00 - 22:00 | 97.96 | 99.46 | 99.44 | 99.45 | 99.47 | 99.35 | 99.36 | 99.06 | 99.25
22:00 - 23:00 | 98.63 | 99.59 | 99.61 | 99.59 | 99.62 | 99.52 | 99.53 | 99.31 | 99.46
23:00 - 24:00 | 98.89 | 99.63 | 99.64 | 99.61 | 99.62 | 99.52 | 99.54 | 99.29 | 99.46

TABLE VII: Classification accuracy of Naive and CTR models (Oct. 23)

Hour | Naive | CTR1 | CTR2 | CTR3 | CTR4 | CTR6 | CTR8 | CTR12 | CTR24
---|---|---|---|---|---|---|---|---|---
03:00 - 04:00 | 99.46 | 99.75 | 99.75 | 99.71 | 99.75 | 99.70 | 99.71 | 99.64 | 99.60
04:00 - 05:00 | 98.62 | 99.55 | 99.56 | 99.50 | 99.54 | 99.50 | 99.54 | 99.39 | 99.22
05:00 - 06:00 | 98.95 | 99.57 | 99.60 | 99.55 | 99.58 | 99.55 | 99.60 | 99.50 | 99.32
06:00 - 07:00 | 99.16 | 99.55 | 99.56 | 99.53 | 99.55 | 99.42 | 99.55 | 99.42 | 99.29
07:00 - 08:00 | 98.43 | 99.31 | 99.34 | 99.30 | 99.25 | 99.09 | 99.27 | 99.01 | 98.74
08:00 - 09:00 | 97.02 | 98.99 | 98.99 | 99.01 | 98.90 | 98.80 | 98.70 | 98.67 | 98.41
09:00 - 10:00 | 97.14 | 99.03 | 99.05 | 98.90 | 98.91 | 98.83 | 98.71 | 98.65 | 98.33
10:00 - 11:00 | 97.66 | 98.85 | 98.67 | 98.75 | 98.77 | 98.67 | 98.55 | 98.50 | 98.15
11:00 - 12:00 | 97.74 | 98.75 | 98.62 | 98.71 | 98.68 | 98.65 | 98.60 | 98.48 | 98.24
12:00 - 13:00 | 97.34 | 98.61 | 98.70 | 98.65 | 98.52 | 98.41 | 98.51 | 98.26 | 98.05
13:00 - 14:00 | 98.43 | 98.81 | 98.84 | 98.86 | 98.71 | 98.68 | 98.66 | 98.55 | 98.40
14:00 - 15:00 | 98.33 | 98.76 | 98.86 | 98.82 | 98.70 | 98.62 | 98.63 | 98.52 | 98.35
15:00 - 16:00 | 97.00 | 98.64 | 98.83 | 98.76 | 98.63 | 98.57 | 98.52 | 98.33 | 98.07
16:00 - 17:00 | 98.47 | 99.31 | 99.33 | 99.36 | 99.31 | 99.15 | 99.28 | 99.08 | 98.99
17:00 - 18:00 | 99.13 | 99.47 | 99.47 | 99.47 | 99.47 | 99.33 | 99.46 | 99.32 | 99.25
18:00 - 19:00 | 99.13 | 99.57 | 99.58 | 99.57 | 99.54 | 99.55 | 99.54 | 99.43 | 99.36
19:00 - 20:00 | 98.20 | 99.35 | 99.42 | 99.42 | 99.34 | 99.38 | 99.36 | 99.21 | 99.06
20:00 - 21:00 | 98.78 | 99.48 | 99.51 | 99.52 | 99.47 | 99.49 | 99.47 | 99.34 | 99.20
21:00 - 22:00 | 97.76 | 99.54 | 99.53 | 99.52 | 99.52 | 99.52 | 99.50 | 99.40 | 99.18
22:00 - 23:00 | 98.69 | 99.59 | 99.59 | 99.58 | 99.58 | 99.57 | 99.56 | 99.46 | 99.40
23:00 - 24:00 | 99.01 | 99.57 | 99.58 | 99.55 | 99.55 | 99.53 | 99.53 | 99.43 | 99.34

###  -B Hourly data distribution

The distribution shifts between different hours are measured using the hamming
distance of commitment decisions and the energy distance of the input
features, shown in Figure 10 and 11, respectively. Each entry indicates the
distance between the corresponding hours. Therefore, the diagonal entries are
always 0 since the distances between the same dataset are $0$. Interestingly,
the block-wise patterns can be observed for energy distances, which is more
obvious for hamming distances. Within the block which consists of several
consecutive hours, the distribution shift tends to be small. Cross the block,
the distribution shift largely. This, to some degree, explains why the CTR
models scale well when the span of the combined dataset is less than certain
hours, from the data distribution perspective.

Figure 10: Heat map of commitment decisions across hours

Figure 11: Heat map of the energy distances of hourly distribution across
hours


# Cosmic-ray ionization rate in protoplanetary disks with sheared magnetic
fields

Yuri I. Fujii Graduate School of Human and Environmental Studies,
Kyoto University, Yoshida-Nihonmatsu, Sakyo, Kyoto 606-8501, Japan  Department
of Physics, Nagoya University, Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8602
Shigeo S. Kimura Frontier Research Institute for Interdisciplinary Sciences,
Tohoku University, Sendai 980-8578, Japan  Astronomical Institute, Graduate
School of Science,
Tohoku University, Sendai 980-8578, Japan

###### Abstract

We investigate the effects of magnetic field configurations on the ionization
rate by cosmic rays in protoplanetary disks. First, we consider cosmic-ray
propagation from the interstellar medium (ISM) to the protoplanetary disks and
showed that the cosmic-ray density around the disk should be 2 times lower
than the ISM value. Then, we compute the attenuation of cosmic rays in
protoplanetary disks. The magnetic fields in the disk are stretched to the
azimuthal directions, and cosmic rays need to detour while propagating to the
midplane. Our results show that the detouring effectively enhances the column
density by about two orders of magnitudes. We employ a typical ionization rate
by cosmic rays in diffuse ISM, which is considered too high to be consistent
with observations of protoplanetary disks, and find that the cosmic rays are
significantly shielded at the midplane. In the case of the disk around IM lup,
the midplane ionization rate is very low for $r\lesssim\,100$ au, while the
value is as large as a diffuse ISM in the outer radii. Our results are
consistent with the recent ALMA observation that indicates the radial gradient
in the cosmic-ray ionization rate. The high ionization rate in the outer radii
of disks may activate the magnetorotational instability that was thought to be
suppressed due to ambipolar diffusion. These results will have a strong
influence on the dynamical and chemical evolutions of protoplanetary disks.

Cosmic rays (329), Protoplanetary disks (1300), Magnetic fields (994),
Ionization (2068)

##  1 Introduction

Protoplanetary disks are the birthplaces of planets. They consist of weakly
ionized gas, and their dynamics are mainly governed by the magnetohydrodynamic
(MHD) processes. Hence, the ionization degree inside the disks determines the
fate of the protoplanetary disks (e.g., Sano et al., 2000; Wardle, 2007;
Gressel et al., 2012, 2015; Béthune et al., 2017; Bai, 2017; Mori et al.,
2019). Ionization also affects the chemical reactions in the disk, which is
crucial to interpreting the observational data of ALMA (e.g., Walsh et al.,
2012; Eistrup et al., 2016; Aikawa et al., 2021; Law et al., 2021; Öberg et
al., 2021; Notsu et al., 2021)

Figure 1: Schematic view of our scenario. An hourglass-shaped magnetic field
configuration focuses cosmic rays from a molecular cloud core to a
protoplanetary disk. The purple oval represents the area from which the disk
can receive cosmic ray particles. The magnetic fields are stretched to
azimuthal direction by shearing motion.

Cosmic rays are an essential ionization source in protoplanetary disks.
Protostellar radiations, including ultraviolet and X-rays, are the most
powerful ionization source in the upper layer of the disk. However, they are
easily blocked by the disk gas, and their ionization rates at the midplane are
strongly suppressed (Igea & Glassgold, 1999; Ercolano & Glassgold, 2013).
Cosmic rays have a stronger penetration power, and thus, they are expected to
be the dominant source of ionization at the midplane of the inner disk (e.g.,
Gammie, 1996; Sano et al., 2000; Wardle, 2007; Okuzumi, 2009; Fujii et al.,
2011), although many previous studies on ionization degree of protoplanetary
disks only consider X-ray ionization (e.g., Glassgold et al., 1997; Ilgner &
Nelson, 2006; Fromang et al., 2002; Turner & Sano, 2008). In the deepest part
of the disk, where even cosmic rays cannot penetrate, the ionization floor is
determined by radionuclides. The ionization rate by cosmic rays in
protoplanetary disks is still controversial due to the uncertainty of cosmic-
ray intensity there. Cleeves et al. (2013, 2014, 2015) discussed the
modulation of cosmic rays by stellar winds from T-Tauri stars, i.e.,
T-Tauriosphere. Cosmic rays can be produced by protostellar activities, which
are known as stellar energetic particles. These may also be an important
ionization source, especially in the surface layer of the disk (Turner &
Drake, 2009; Padovani et al., 2015, 2016; Rab et al., 2017; Fraschetti et al.,
2018).

Most of the previous studies assumed that the energy density of cosmic rays at
the surface of a protoplanetary disk is the same as that in an interstellar
medium (e.g., Gammie, 1996; Sano et al., 2000; Furuya et al., 2013). They also
consider that cosmic rays travel to the disk midplane in the shortest path,
i.e., cosmic rays enter the protoplanetary disk straightly from the vertical
direction. In reality, the shear motion in protoplanetary disks stretches the
magnetic field, and cosmic rays travel only along the magnetic field lines due
to their small gyration radii (see Figure 1). Thus, the cosmic rays should
detour to reach into the protoplanetary disks, which changes both the cosmic-
ray energy density and the penetration depth. Padovani et al. (2018) estimated
the cosmic-ray ionization rate in the inner radii of a protoplanetary disk
with taking the effect of the disk magnetic field into account. In this
Letter, we consider the magnetic fields that are connected to the molecular
cloud and discuss the effects of sheared field lines on the cosmic-ray
ionization rate in global protoplanetary disks.

##  2 Focusing and reflecting cosmic ray particles

We consider the magnetic field configuration shown in Figure 1, where magnetic
fields threading the protoplanetary disk are connected to its parent molecular
cloud. This configuration is naturally expected if we consider disk formation
by the gravitational collapse of the cloud core. Observations of some
molecular cloud cores show an hourglass-shaped magnetic field configuration,
which may support this picture (e.g., Maury et al., 2018). This configuration
is similar to those used by Padovani & Galli (2011); Padovani et al. (2013).
In contrast, the magnetic fields in the solar wind extend horizontally with a
spiral shape and are disconnected from those in the ISM (Parker, 1958;
Peterson et al., 2019). If the magnetic field configuration around a
protoplanetary disk is similar to that of the solar wind, the cosmic-ray
density can be modulated due to T-Tauriosphere as discussed in Cleeves et al.
(2013, 2014, 2015). However, in the case that magnetic fields are connected to
the ISM, cosmic rays can travel along the field lines (see Section 5).

Gyro radius of a cosmic ray proton is $r_{\rm G}=E/(eB)\sim
0.2~{}(E/10{\rm~{}GeV})(B/10{\rm~{}\mu G})^{-1}$ au, where $E$ is the energy
of cosmic ray protons, $e$ is the elementary charge, $B$ is the magnetic field
strength, and we use $B=10\rm~{}\mu G$ in a molecular cloud (e.g. Crutcher et
al., 2010). This is much smaller than the scale of the cloud core ($\sim 0.1$
pc $\sim 2\times 10^{4}$ au), which can be regarded as the turbulence
scale.111We consider that the kinetic-scale turbulence can be damped due to
the low ionization rate in the cloud core and the disk. Also, the turbulence
should be injected in the largest scale and the power spectrum of MHD
turbulence is expected to be steep (Goldreich & Sridhar, 1995). Thus, cosmic
rays of $E\lesssim 10$ TeV travel along magnetic field lines. We assume that
cosmic-ray propagation is in the free-streaming regime and ignore the
diffusion (see Section 5). The disk bundles up the magnetic field lines from a
large area222During the star and disk formation, a part of the magnetic fields
slips out from the disk. We use the magnetic field strength after the disk
formation in our estimate.. If we assume that the cosmic-ray flux is
proportional to the magnetic field strength, which is true for an isotropic
cosmic-ray distribution and inefficient CR diffusion perpendicular to the
magnetic field lines, the cosmic-ray flux at the disk surface is enhanced by a
factor of

$\frac{S_{\rm cloud}}{S_{\rm disk}}=\frac{B_{\rm disk}}{B_{\rm cloud}}\sim
10^{3},$ (1)

where $S_{\rm cloud}$ is the surface area in the molecular cloud from which
the protoplanetary disk can collect cosmic rays (the purple area in Figure 1),
$S_{\rm disk}$ is the total area of the disk, and, $B_{\rm disk}\sim 10$ mG
and $B_{\rm cloud}\sim 10\mu$G are the typical strengths of magnetic fields in
a protoplanetary disk (at $r\sim 10$ au, $z\sim 4H$, where $H$ is the scale
height of the disk; cf. Suzuki et al., 2010)) and molecular cloud core,
respectively.

In accordance with the conservation of kinetic energy and magnetic moment of a
cosmic ray particle, pitch angles of the particles in a molecular cloud core,
$\alpha_{\rm cloud}$, and that in a protoplanetary disk, $\alpha_{\rm disk}$,
have the following relation (e.g. Padovani & Galli, 2011),

$\frac{{\rm sin}^{2}\alpha_{\rm disk}}{{\rm sin}^{2}\alpha_{\rm
cloud}}=\frac{B_{\rm disk}}{B_{\rm cloud}}.$ (2)

If we adopt the value in Equation (1), the pitch angle in the cloud core
should satisfy sin$\alpha_{\rm cloud}\leq 0.03$, which means that only the
cosmic ray particles that travel almost perfectly parallel to the field line
can enter the protoplanetary disk. All the other particles are reflected by
the magnetic mirror effect, and the cosmic-ray flux is reduced by a factor of
$\alpha^{2}/2$. Because of the balance between the focusing and reflecting
effects, the net cosmic-ray flux at the disk surface ends up with a half of
that in the parental cloud core (see Silsbee et al., 2018, for more detailed
treatment). This factor 2 is generic as long as $B_{\rm cloud}/B_{\rm disk}\ll
1$. We define this factor for diminishing CR flux as $\Gamma_{\rm dim}$.

##  3 Disk Structure

Below the disk surface, the magnetic fields are stretched to the azimuthal
direction due to the shaer motion of the disk (see Figure 1). In order to
obtain the sheared field structure and the stratified density distribution of
a protoplanetary disk, we perform magnetohydrodynamic (MHD) simulations. The
circumstance under which magnetorotational instability (MRI) generates
turbulence is still controversial (Bai & Stone, 2013; Bai, 2013; Gressel et
al., 2015; Simon et al., 2015; Mori et al., 2019). The activation of MRI is
strongly affected by the ionization-rate distribution, which will be obtained
in this work. In this Letter, we consider a non-turbulent protoplanetary disk
because it leads to the most efficient ionization by cosmic rays (see Section
5 for the effect of turbulence). We conduct low-resolution shearing-box
simulations with ideal MHD approximation (Goldreich & Lynden-Bell, 1965;
Hawley et al., 1995) to mimic the non-turbulent disk structure.

Figure 2:  The top panel shows the ratio of the azimuthal and vertical
components of the magnetic field strength at each height from the midplane.
The purple jagged line represents $\langle B_{y}/B_{z}\rangle$. We also plot
$\langle B_{y}\rangle/\langle B_{z}\rangle$ in green for comparison. The
bottom panel presents the effective column densities that cosmic ray particles
have to go through from the closer surface to reach the corresponding heights.
The dashed curves are the vertical column densities that are used in the
conventional estimates, and the dotted horizontal line represents the
attenuation length of cosmic rays.

The simulations were performed with Athena MHD code (Stone et al., 2008; Stone
& Gardiner, 2010). The benefit of the ideal MHD approximation is that
simulations can be done with normalized physical values (Keplerian frequency,
scale height of the disk, and the midplane density are set to be unity). We
can freely scale the unit to arbitrary radii in the disk. The resolution was
$30\times 60\times 240$ for a local box with the size of $\pm 0.5H$, $\pm 1H$,
and $\pm 4H$ in the $x$, $y$, and $z$ directions, respectively. The
$z$-component of the initial plasma beta in the midplane of the disk
(corresponding to the net vertical field) is parameterized as $\beta_{\rm
z}=10^{4}$ and $10^{6}$.

The top panel of Figure 2 shows the ratio of the $y$ and $z$ components of the
magnetic fields as a function of the distance from the midplane for
$\beta_{\rm z}=10^{6}$. One can see that the shear elongates the field lines
to the azimuthal direction by a factor of $100-1000$. This means that cosmic
ray particles traveling along the magnetic fields have to go through
$\Gamma_{\rm detour}=\int\langle B_{y}/B_{z}\rangle dz\sim 100-1000$ times
longer distance to arrive at a certain height than the vertically straight
path. Although the absolute values of the field strength are different between
the cases with $\beta_{\rm z}=10^{4}$ and $10^{6}$, the values of $\Gamma_{\rm
detour}$ are similar. Thus, our results can be applied to the disks that have
typical magnetic fields ($\beta_{\rm z}=10^{4}-10^{6}$). We cautiously note
that the averaging procedure will affect the detouring factor; $\langle
B_{y}/B_{z}\rangle$ is different from $\langle B_{y}\rangle/\langle
B_{z}\rangle$. When we compute the values in the brackets, we first
horizontally average the data and take a time average for 50 orbital periods
from the point of 350 orbits after the beginning of the simulation.

Taking the effect of detouring into consideration, we estimate the effective
column density that cosmic rays have to go through to enter into the disk. The
results are presented in the bottom panel of Figure 2 for the disk around IM
Lup (Zhang et al., 2021) as an example of the observed protoplanetary disks.
With conventional treatment, cosmic rays of attenuation length of $\sim 100$ g
cm-2 (Umebayashi & Nakano, 1981, see Padovani et al. (2018) for more updated
values for the high-density region) can arrive at the midplane for $r\gtrsim$
10 au unless the shielding effect by T-Tauriosphere is effective. On the other
hand, our treatment indicates that cosmic rays cannot arrive at the disk
midplane even at $\sim$ 100 au for IM Lup.

In our simulations, the toroidal component of magnetic fields is generated at
all heights due to the ideal MHD approximation. In the very dense region with
poor coupling with magnetic fields, i.e., near the midplane around a few au
from the star, the field lines would not be stretched by the velocity shear.
Then, cosmic rays do not have to detour much. In such a region, however,
cosmic rays are easily shielded, and our conclusions are unchanged.

##  4 Effects on the ionization rates

Figure 3: Top panel: the radial distribution of ionization rate at $2H$ from
the midplane based on the new estimate of the CR ionization rate (solid green)
with $\zeta_{\rm ISM}=3\times 10^{-16}$ s-1 using the effective column density
along the elongated path, together with the contribution of X-ray (orange) and
short-lived radionuclides (SLR; pink) is displayed with the thick blue curve.
The conventional estimate of $\zeta_{\rm CR}$ based on $\zeta_{\rm
ISM}=10^{-17}$ s-1 with $\Gamma_{\rm detour}=1$ (the dotted green line) makes
up the thin blue curve of the conventional ionization rate. Middle panel: the
same with the top panel, but for the midplane of the disk. The dashed pink
line is the ionization rate of the long-lived radionuclides (LLR). Bottom
panel: the total ionization rate at $z=3H,\ 2H,\ 1H$, and midplane using the
CR ionization rate of this work. The case with the absence of SLR and the
estimate using the conventional CR ionization rate are plotted in light blue
and gray, respectively. The x-axis in each plot starts from 10au.

Figure 4: The distribution of ionization rate for our new model considering
the propagation of cosmic rays along magnetic fields and that for the
conventional model considering the vertical attenuation of cosmic rays in the
disk around IM Lup. The white lines represent the scale height of the disk
(which is $\sqrt{2}$ times larger than the definition of Zhang et al. (2021))
and the green dots mark the regions where the cosmic rays are the dominant
ionization source.

We show a significant enhancement of the effective column density in Figure 2.
Due to the detouring of the cosmic rays along the sheared field lines, cosmic
ray particles stay in the surface layer for a longer time before they enter
into the disk. The cosmic-ray ionization rate at the point of interest within
the disk can be calculated as

$\displaystyle\zeta_{\rm CR}=\frac{1}{2}\Gamma_{\rm dim}\zeta_{\rm
ISM}\left\\{\exp({-\Gamma_{\rm detour}\chi_{\rm eff,u}/\chi_{\rm CR}})\right.$
$\displaystyle\left.+\exp({-\Gamma_{\rm detour}\chi_{\rm eff,l}/\chi_{\rm
CR}})\right\\},$ (3)

where $\zeta_{\rm ISM}$ is the ionization rate of cosmic rays in the
interstellar medium and $\chi_{\rm CR}=96$ g cm-2 is the attenuation length of
cosmic rays (Umebayashi & Nakano, 1981). We define the vertical column density
from the point to the upper/lower surface of the disk as $\chi_{\rm eff,u}$
and $\chi_{\rm eff,l}$, respectively, and $\Gamma_{\rm detour}$ is the ratio
of the effective and vertical column density from the point to the closer
surface, i.e. $\Gamma_{\rm detour}=1$ in conventional studies and $\Gamma_{\rm
detour}\chi_{\rm eff}=\int\langle\rho B_{y}/B_{z}\rangle dz$ in our work. Our
results suggest $\Gamma_{\rm detour}\sim 100$ when the cosmic rays travel
along the sheared magnetic fields. Intuitively, one would think that the CR
density can be enhanced by the sheared magnetic field by a factor of
$\Gamma_{\rm detour}$. However, the sheared magnetic field makes the mirroring
effect effective, and thus, the CR density does not change (see Section 2).

As an example, we estimate the ionization rate assuming the surface density
and temperature distributions of the disk around IM lup (Zhang et al., 2021).
We consider 4 ionization sources: cosmic rays, X-rays, and short-lived (SLR)
and long-lived radionuclides (LLR). The cosmic-ray ionization rate is given by
Equation (3). The typical value used in conventional studies is $\zeta_{\rm
ISM}=10^{-17}$ s-1 (Spitzer & Tomasko, 1968; Cummings et al., 2016), and we
adopt this for the conventional model. Observations of diffuse clouds in
Galactic disk suggest higher values (Padovani et al., 2020), however, those
values are not consistent with lower ionization rate in the midplane of
protoplanetary disks (Aikawa et al., 2021; Seifert et al., 2021). Here, we
employ $\zeta_{\rm ISM}=3\times 10^{-16}$ s-1(Indriolo et al., 2015) for our
new estimates. The value of $\Gamma_{\rm detour}$ is insensitive to the
averaging method; the calculations of $\langle\rho\ B_{y}/B_{z}\rangle$ and
$\langle\rho\rangle\langle B_{y}/B_{z}\rangle$ lead to similar results. We use
the X-ray ionization rate by Bai (2011) that considers both the direct and
scattered emissions based on Igea & Glassgold (1999). We use X-ray luminosity
$L_{X}=10^{30}$ erg/s (Cleeves et al., 2017) and plasma temperature $T_{X}=3$
keV. The ionization rate by SLR and LLR are $\zeta_{\rm SLR}=7.6\times
10^{-19}\rm~{}s^{-1}$ and $\zeta_{\rm LLR}=1.4\times 10^{-22}\rm~{}s^{-1}$,
respectively (Umebayashi & Nakano, 2009).

The top panel of Figure 3 shows the ionization rates as a function of $r$ at
$z=2H$. We can see that cosmic rays are attenuated above $z=2H$, and the X-ray
ionization is dominant for $r\lesssim 50$ au. On the other hand, the
conventional treatment leads to the opposite trend: the X-ray ionization is
dominant at outer radii, while cosmic-ray ionization is important in the inner
part of the disk.

The modification of the ionization rate in the midplane is summarized in the
middle panel of Figure 3. In the conventional model, cosmic rays are the
dominant ionization source for all the radii. On the other hand, the cosmic-
ray ionization rate of our model at $r\lesssim 100$ au is much lower than that
of the conventional model because of the significant reduction of the cosmic-
ray flux. The ionization is dominated by SLR there. The X-ray ionization is
sub-dominant in all radii, because of the significant attenuation.

$\zeta_{\rm ISM}$ higher than the conventional model does not affect the total
ionization rate in the upper layers of $z\gtrsim 3H$, because the X-ray
ionization rate is higher at the upper layer of the disk. When the X-ray flux
is much smaller, the ionization rate may be dominated by the adopted
$\zeta_{\rm ISM}$ even in the higher atmosphere of the inner disk.

The half-life of the most effective species among SLR, 26Al, is 7 $\times
10^{5}$ years, and they might not necessarily still be present in the
protoplanetary disks of interest (See Section 5 for the discussion on the
production of 26Al). Thus, we calculate the ionization rate of the disk in the
absence of the SLR (the bottom panel of Figure 3). The ionization is actually
led by LLR within $\sim$ 20 au in this case. Then, X-rays play an important
role for 20 au $\lesssim r\lesssim$ 70 au, and cosmic rays are dominant for
$r\gtrsim 70$ au. Umebayashi & Nakano (2009) pointed out that the ionization
by radionuclides suddenly becomes inefficient when the mean size of the dust
grains grows to about $>1$ cm because the radiation is absorbed by the grains.
In such a case, ionization would be governed by X-rays while cosmic rays are
highly attenuated.

The 2D map of the ionization rate for our new model and that for the
conventional estimate are shown in Figure 4. The ionization rate for our model
is lower than that for the conventional model by about an order of magnitude
for $z\lesssim 1H$ and $r\lesssim 100$ au due to the detouring of cosmic rays.
The ionization rate in this region is dominated by SLR in our model, while the
cosmic-ray ionization is dominant in the conventional model (see the
distribution of green dots in the figure). On the other hand, the ionization
rate for our model is higher than that for the conventional model for $z\sim
2H$ or the midplane for $r\gtrsim 120$ au, owing to the larger value of
$\zeta_{\rm ISM}$.

##  5 Discussion

The suppression of the energy density of CRs with energies $<10$ GeV is
observed in the solar system (Potgieter, 2013). The situation in the
protostellar system is more speculative. Some claim that winds from young
stars should suppress the CRs of similar or even higher energies (Cohen et
al., 2012; Rodgers-Lee et al., 2020b). However, we consider that the T-Tauri
winds are unlikely to affect the CR density in the protoplanetary disks. The
mass outflow rate from T-Tauri stars is much lower than the winds from
protoplanetary disks, and then, the T-Tauri winds cannot push out the disk
winds. Thus, the disk winds fill the region above the disk and may modulate
the CR density. We speculate that the magnetic field configuration is crucial
for the suppression of CRs. If the disk winds are driven by the magnetic
centrifugal force (Blandford & Payne, 1982; Bai et al., 2016), the magnetic
field lines anchored in the disk can be connected to the ISM, and the magnetic
field lines are not completely dominated by the toroidal component. Then, the
cosmic rays are likely able to enter the disk by the free-streaming or
diffusion parallel to the magneitc field lines. On the other hand, if the disk
winds are driven by other mechanisms, such as turbulence heating (Suzuki et
al., 2010) or photoevaporation (Owen et al., 2010; Kimura et al., 2016a;
Kunitomo et al., 2020, 2021), the magnetic field lines are unlikely connected
to the ISM and likely dominated by the toroidal component. Then, the cosmic
rays need to cross the magnetic field lines many times to reach the disk,
which may suppress low-energy cosmic rays as in the solar winds.

Stellar energetic particles are regarded as a possible ionization source in
protoplanetary disks. Depending on the settings of calculations, such as the
particle energy and the treatment of propagation, they may be the
dominant/important ionization source in the surface layer of the disk (Turner
& Drake, 2009; Rab et al., 2017; Fraschetti et al., 2018) or even at the
midplane of the disk (Rodgers-Lee et al., 2020a). However, even for high
energy particles of $E\sim 10$ GeV, the Larmor radius is much smaller than the
scale height, and they need to detour along the sheared magnetic fields as
long as the diffusion perpendicular to the magnetic field lines are
inefficient. Based on our scenario, stellar energetic particles may
significantly contribute to the ionization rate at the disk surface, but they
are unlikely to reach the midplane of the disk. Recent studies suggest 26Al
can be produced by stellar energetic particles (Gaches et al., 2020; Adams,
2021). However, since the stellar energetic particles cannot reach the disk
midplane, the newly synthesized 26Al is unlikely to affect the ionization
rates in the range of our interest.

An ALMA large program, ‘Molecules with ALMA at Planet-forming Scales’ (MAPS)
discussed the ionization rates at the midplane of protoplanetary disks (Aikawa
et al., 2021). Their results suggest $\zeta\sim 10^{-18}\rm~{}s^{-1}$ and also
show the variation among disks of a similar mass range. This might be caused
by the different configurations of magnetic-field lines or the turbulent state
of the disks. Another ALMA observation by Seifert et al. (2021) reported a
lower ionization rate in the inner disk and a higher value in the outer radii
with the transition at $\sim$80–100 au in the disk around IM lup. Our model
predicts a similar tendency as seen in Figure 4.

Some protoplanetary disks show lower values of the CO abundance compared to
the canonical ISM value of $10^{-4}$ in the layer between the CO
photodissociation layer in the surface and the CO freeze-out layer near the
midplane (e.g., Favre et al., 2013; Miotello et al., 2017). This may suggest a
high ionization rate of $\zeta\geq 10^{-17}$s-1 and low temperature ($T\sim
20-30$ K) in the disks because CO molecules are converted into less volatile
molecules in such an environment (Furuya & Aikawa, 2014; Bosman et al., 2018;
Schwarz et al., 2018). Our model with the large value of $\zeta_{\rm ISM}$
predicts higher ionization rates at the middle layer ($\sim 2H$) in the outer
disks ($\sim$100 au) compared to the canonical cosmic-ray ionization rate,
while keeping the midplane ionization rate as low as values suggested by ALMA
observations.

Figure 5: Radial distributions of Ambipolar Elsasser number at the midplane of
the disk around IM Lup.

Bai & Stone (2013); Bai (2013); Simon et al. (2015) showed that ambipolar
diffusion weakens the MRI turbulence at the outer radii of disks. We calculate
ambipolar Elsasser number Am $=v_{A}^{2}/(\eta_{A}\Omega)$, where $v_{A}$,
$\eta_{A}$, $\Omega$ are Alfven velocity, ambipolar diffusivity, and Keplerian
frequency, respectively. We use the formulae of Wardle (2007) for $\eta_{A}$
and compute the ionization rate with the simplified method of Fujii et al.
(2011). Figure 5 shows Am in the IM Lup disk. Our ionization model suggests Am
$\lesssim 1$, where we expect a purely laminar disk, within $\sim 150$ au. The
tendency that the disk is laminar at $\sim$ 100 au is consistent with the
recent ALMA observation suggesting very weak turbulence due to the small scale
height of the dust disk in the disk around Oph 163131 (e.g., Villenave et al.,
2022). In the further outer disk, where previous studies predict weak MRI
turbulence due to ambipolar diffusion (Bai & Stone, 2013; Bai, 2013; Simon et
al., 2015), our model expects Am $\gtrsim$ 10\. This may fully activate the
MRI turbulence. Depending on the location, the Hall effect might be the
dominant non-ideal MHD effect (Wardle & Salmeron, 2012). Non-ideal MHD
simulations with our ionization model are required to further investigate the
level of turbulence.

The existence of turbulence affects cosmic-ray transport. We have assumed that
the magnetic fields have no fluctuations and cosmic-ray transport is in the
free-streaming regime. If some mechanism, such as MRI, drives turbulence,
fluctuations of magnetic fields can scatter the cosmic rays. Then, the cosmic-
ray transport should be in the diffusion regime. Even in the diffusion regime,
cosmic rays move along the magnetic field lines because they have an
anisotropic diffusion coefficient with the parallel diffusion much faster than
the perpendicular diffusion, as shown in particle simulations in MRI
turbulence (Kimura et al., 2016b, 2019; Sun & Bai, 2021). Then, the particle
trajectory can be expressed by the random walk along the field line. Thus, the
path length to the disk midplane in the diffusion regime should be much longer
than in the streaming regime. In the diffusion regime, cosmic rays likely lose
their energy at a further upper region of the disk, and then, the ionization
degree around the midplane of the disk would become low. This leads to
suppression of turbulence, and cosmic rays would become able to enter the disk
again. In summary, we start the discussion from the laminar magnetic field
(implying suppression of turbulence). Then, cosmic rays will enhance the
ionization rate, and the disk will become the ideal MHD regime that drives
turbulence. The turbulence would shield the cosmic rays, which leads to the
non-ideal MHD regime in which turbulence is suppressed. To consistently
understand the evolution of protoplanetary disks, we need to solve the
magnetohydrodynamics, cosmic-ray transport, and ionization rate
simultaneously.

##  6 Summary

We have investigated the effects of the magnetic-field configurations on
cosmic-ray ionization rates in protoplanetary disks. We have first considered
the cosmic-ray propagation from ISM to the disk, where magnetic field strength
changes by several orders of magnitude (see Figure 1). The cosmic-ray density
at the disk surface should decrease by a factor of $\Gamma_{\rm dim}\sim 1/2$
from the ISM value because of the balance between the magnetic-flux focusing
and reflection by a magnetic mirror.

The magnetic fields are stretched to the azimuthal direction in the disk due
to the shearing motion, which affects the cosmic-ray distributions in the
protoplanetary disks. Cosmic rays need to detour to reach the disk midplane
due to the shearing magnetic field, and the vertical penetration depth should
decrease. We have performed a set of local MHD simulations to obtain the
detouring path length and showed that the cosmic-ray density at midplane for
$r\lesssim 100$ au can be significantly suppressed even when we employ
$\zeta_{\rm ISM}=3\times 10^{-16}$ s-1. Our results suggest a totally
different ionization rate distribution from that of the conventional model.
This may change the chemical and dynamical evolution of the protoplanetary
disks, which is crucial to understanding the star and planet formation
processes.

We acknowledge Satoshi Okuzumi, Yusuke Tsukamoto, Jiro Shimoda, Yuri Aikawa,
Kenji Furuya, Yoshihide Yamato, and Akimasa Kataoka for useful discussions. We
are grateful to Soonyoung Roh and Marco Padovani for fruitful discussions at
the beginning of this project. We also thank Kedron Silsbee and Alexei Ivlev
for useful comments. This work is partly supported by JSPS KAKENHI Grant
Number 19H05077, 22K14086 (Y.I.F), 19J00198, 21H04487, and 22K14028 (S.S.K.).
Y.I.F was supported in-part by the Program for the Development of Next-
generation Leading Scientists with Global Insight (L-INSIGHT), sponsored by
the Ministry if Education, Culture, Sports, Science and Technology (MXT),
Japan. S.S.K. acknowledge the support by the Tohoku Initiative for Fostering
Global Researchers for Interdisciplinary Sciences (TI-FRIS) of MEXT’s
Strategic Professional Development Program for Young Researchers. Numerical
computations were in part carried out on Cray XC50 at Center for Computational
Astrophysics, National Astronomical Observatory of Japan.

##  Appendix A Minimum Mass Solar Nebula Model

For the comparison with theoretical studies adopting the Minimum Mass Solar
Nebula (MMSN) model (Hayashi, 1981), we provide the effective column density
for cosmic rays to arrive at various heights from the midplane of MMSN.
Conventionally, cosmic rays are thought to easily reach the midplane at $\sim$
10 au. As one can see in Fig.6, however, they can barely penetrate into the
midplane at 100 au in our new model.

Figure 6:  The effective column densities that cosmic rays have to pass
through from the closer surface to reach each height of the MMSN at 1, 10, and
100au. The dashed black line is the vertically measured column density that is
used in conventional studies for 1au. The thin horizontal line shows the
attenuation length of cosmic rays.

## References

  * Adams (2021) Adams, F. C. 2021, ApJ, 919, 10, doi: 10.3847/1538-4357/ac1111
  * Aikawa et al. (2021) Aikawa, Y., Cataldi, G., Yamato, Y., et al. 2021, ApJS, 257, 13, doi: 10.3847/1538-4365/ac143c
  * Bai (2011) Bai, X.-N. 2011, ApJ, 739, 50, doi: 10.1088/0004-637X/739/1/50
  * Bai (2013) —. 2013, ApJ, 772, 96, doi: 10.1088/0004-637X/772/2/96
  * Bai (2017) —. 2017, ApJ, 845, 75, doi: 10.3847/1538-4357/aa7dda
  * Bai & Stone (2013) Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76, doi: 10.1088/0004-637X/769/1/76
  * Bai et al. (2016) Bai, X.-N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152, doi: 10.3847/0004-637X/818/2/152
  * Béthune et al. (2017) Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75, doi: 10.1051/0004-6361/201630056
  * Blandford & Payne (1982) Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883, doi: 10.1093/mnras/199.4.883
  * Bosman et al. (2018) Bosman, A. D., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 618, A182, doi: 10.1051/0004-6361/201833497
  * Cleeves et al. (2013) Cleeves, L. I., Adams, F. C., & Bergin, E. A. 2013, ApJ, 772, 5, doi: 10.1088/0004-637X/772/1/5
  * Cleeves et al. (2014) Cleeves, L. I., Bergin, E. A., & Adams, F. C. 2014, ApJ, 794, 123, doi: 10.1088/0004-637X/794/2/123
  * Cleeves et al. (2017) Cleeves, L. I., Bergin, E. A., Öberg, K. I., et al. 2017, ApJ, 843, L3, doi: 10.3847/2041-8213/aa76e2
  * Cleeves et al. (2015) Cleeves, L. I., Bergin, E. A., Qi, C., Adams, F. C., & Öberg, K. I. 2015, ApJ, 799, 204, doi: 10.1088/0004-637X/799/2/204
  * Cohen et al. (2012) Cohen, O., Drake, J. J., & Kóta, J. 2012, ApJ, 760, 85, doi: 10.1088/0004-637X/760/1/85
  * Crutcher et al. (2010) Crutcher, R. M., Wandelt, B., Heiles, C., Falgarone, E., & Troland, T. H. 2010, ApJ, 725, 466, doi: 10.1088/0004-637X/725/1/466
  * Cummings et al. (2016) Cummings, A. C., Stone, E. C., Heikkila, B. C., et al. 2016, ApJ, 831, 18, doi: 10.3847/0004-637X/831/1/18
  * Eistrup et al. (2016) Eistrup, C., Walsh, C., & van Dishoeck, E. F. 2016, A&A, 595, A83, doi: 10.1051/0004-6361/201628509
  * Ercolano & Glassgold (2013) Ercolano, B., & Glassgold, A. E. 2013, MNRAS, 436, 3446, doi: 10.1093/mnras/stt1826
  * Favre et al. (2013) Favre, C., Cleeves, L. I., Bergin, E. A., Qi, C., & Blake, G. A. 2013, ApJ, 776, L38, doi: 10.1088/2041-8205/776/2/L38
  * Fraschetti et al. (2018) Fraschetti, F., Drake, J. J., Cohen, O., & Garraffo, C. 2018, ApJ, 853, 112, doi: 10.3847/1538-4357/aaa48b
  * Fromang et al. (2002) Fromang, S., Terquem, C., & Balbus, S. A. 2002, MNRAS, 329, 18, doi: 10.1046/j.1365-8711.2002.04940.x
  * Fujii et al. (2011) Fujii, Y. I., Okuzumi, S., & Inutsuka, S. 2011, ApJ, 743, 53, doi: 10.1088/0004-637X/743/1/53
  * Furuya & Aikawa (2014) Furuya, K., & Aikawa, Y. 2014, ApJ, 790, 97, doi: 10.1088/0004-637X/790/2/97
  * Furuya et al. (2013) Furuya, K., Aikawa, Y., Nomura, H., Hersant, F., & Wakelam, V. 2013, ApJ, 779, 11, doi: 10.1088/0004-637X/779/1/11
  * Gaches et al. (2020) Gaches, B. A. L., Walch, S., Offner, S. S. R., & Münker, C. 2020, ApJ, 898, 79, doi: 10.3847/1538-4357/ab9a38
  * Gammie (1996) Gammie, C. F. 1996, ApJ, 457, 355, doi: 10.1086/176735
  * Glassgold et al. (1997) Glassgold, A. E., Najita, J., & Igea, J. 1997, ApJ, 480, 344, doi: 10.1086/303952
  * Goldreich & Lynden-Bell (1965) Goldreich, P., & Lynden-Bell, D. 1965, MNRAS, 130, 125
  * Goldreich & Sridhar (1995) Goldreich, P., & Sridhar, S. 1995, ApJ, 438, 763, doi: 10.1086/175121
  * Gressel et al. (2012) Gressel, O., Nelson, R. P., & Turner, N. J. 2012, MNRAS, 422, 1140, doi: 10.1111/j.1365-2966.2012.20701.x
  * Gressel et al. (2015) Gressel, O., Turner, N. J., Nelson, R. P., & McNally, C. P. 2015, ApJ, 801, 84, doi: 10.1088/0004-637X/801/2/84
  * Hawley et al. (1995) Hawley, J. F., Gammie, C. F., & Balbus, S. A. 1995, ApJ, 440, 742
  * Hayashi (1981) Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35, doi: 10.1143/PTPS.70.35
  * Igea & Glassgold (1999) Igea, J., & Glassgold, A. E. 1999, ApJ, 518, 848, doi: 10.1086/307302
  * Ilgner & Nelson (2006) Ilgner, M., & Nelson, R. P. 2006, A&A, 445, 205, doi: 10.1051/0004-6361:20053678
  * Indriolo et al. (2015) Indriolo, N., Neufeld, D. A., Gerin, M., et al. 2015, ApJ, 800, 40, doi: 10.1088/0004-637X/800/1/40
  * Kimura et al. (2016a) Kimura, S. S., Kunitomo, M., & Takahashi, S. Z. 2016a, MNRAS, 461, 2257, doi: 10.1093/mnras/stw1531
  * Kimura et al. (2016b) Kimura, S. S., Toma, K., Suzuki, T. K., & Inutsuka, S.-i. 2016b, ApJ, 822, 88, doi: 10.3847/0004-637X/822/2/88
  * Kimura et al. (2019) Kimura, S. S., Tomida, K., & Murase, K. 2019, MNRAS, 485, 163, doi: 10.1093/mnras/stz329
  * Kunitomo et al. (2021) Kunitomo, M., Ida, S., Takeuchi, T., et al. 2021, ApJ, 909, 109, doi: 10.3847/1538-4357/abdb2a
  * Kunitomo et al. (2020) Kunitomo, M., Suzuki, T. K., & Inutsuka, S.-i. 2020, MNRAS, 492, 3849, doi: 10.1093/mnras/staa087
  * Law et al. (2021) Law, C. J., Loomis, R. A., Teague, R., et al. 2021, ApJS, 257, 3, doi: 10.3847/1538-4365/ac1434
  * Maury et al. (2018) Maury, A. J., Girart, J. M., Zhang, Q., et al. 2018, MNRAS, 477, 2760, doi: 10.1093/mnras/sty574
  * Miotello et al. (2017) Miotello, A., van Dishoeck, E. F., Williams, J. P., et al. 2017, A&A, 599, A113, doi: 10.1051/0004-6361/201629556
  * Mori et al. (2019) Mori, S., Bai, X.-N., & Okuzumi, S. 2019, ApJ, 872, 98, doi: 10.3847/1538-4357/ab0022
  * Notsu et al. (2021) Notsu, S., van Dishoeck, E. F., Walsh, C., Bosman, A. D., & Nomura, H. 2021, A&A, 650, A180, doi: 10.1051/0004-6361/202140667
  * Öberg et al. (2021) Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ApJS, 257, 1, doi: 10.3847/1538-4365/ac1432
  * Okuzumi (2009) Okuzumi, S. 2009, ApJ, 698, 1122, doi: 10.1088/0004-637X/698/2/1122
  * Owen et al. (2010) Owen, J. E., Ercolano, B., Clarke, C. J., & Alexander, R. D. 2010, MNRAS, 401, 1415, doi: 10.1111/j.1365-2966.2009.15771.x
  * Padovani & Galli (2011) Padovani, M., & Galli, D. 2011, A&A, 530, A109
  * Padovani et al. (2013) Padovani, M., Hennebelle, P., & Galli, D. 2013, A&A, 560, A114
  * Padovani et al. (2015) Padovani, M., Hennebelle, P., Marcowith, A., & Ferrière, K. 2015, A&A, 582, L13
  * Padovani et al. (2018) Padovani, M., Ivlev, A. V., Galli, D., & Caselli, P. 2018, A&A, 614, A111, doi: 10.1051/0004-6361/201732202
  * Padovani et al. (2016) Padovani, M., Marcowith, A., Hennebelle, P., & Ferrière, K. 2016, A&A, 590, A8
  * Padovani et al. (2020) Padovani, M., Ivlev, A. V., Galli, D., et al. 2020, Space Sci. Rev., 216, 29, doi: 10.1007/s11214-020-00654-1
  * Parker (1958) Parker, E. N. 1958, ApJ, 128, 664, doi: 10.1086/146579
  * Peterson et al. (2019) Peterson, E. E., Endrizzi, D. A., Beidler, M., et al. 2019, Nature Physics, 15, 1095, doi: 10.1038/s41567-019-0592-7
  * Potgieter (2013) Potgieter, M. S. 2013, Living Reviews in Solar Physics, 10, 3, doi: 10.12942/lrsp-2013-3
  * Rab et al. (2017) Rab, C., Güdel, M., Padovani, M., et al. 2017, A&A, 603, A96, doi: 10.1051/0004-6361/201630241
  * Rodgers-Lee et al. (2020a) Rodgers-Lee, D., Taylor, A. M., Downes, T. P., & Ray, T. P. 2020a, MNRAS, 491, 4742, doi: 10.1093/mnras/stz3266
  * Rodgers-Lee et al. (2020b) Rodgers-Lee, D., Vidotto, A. A., Taylor, A. M., Rimmer, P. B., & Downes, T. P. 2020b, MNRAS, 499, 2124, doi: 10.1093/mnras/staa2737
  * Sano et al. (2000) Sano, T., Miyama, S. M., Umebayashi, T., & Nakano, T. 2000, ApJ, 543, 486, doi: 10.1086/317075
  * Schwarz et al. (2018) Schwarz, K. R., Bergin, E. A., Cleeves, L. I., et al. 2018, ApJ, 856, 85, doi: 10.3847/1538-4357/aaae08
  * Seifert et al. (2021) Seifert, R. A., Cleeves, L. I., Adams, F. C., & Li, Z.-Y. 2021, ApJ, 912, 136, doi: 10.3847/1538-4357/abf09a
  * Silsbee et al. (2018) Silsbee, K., Ivlev, A. V., Padovani, M., & Caselli, P. 2018, ApJ, 863, 188, doi: 10.3847/1538-4357/aad3cf
  * Simon et al. (2015) Simon, J. B., Lesur, G., Kunz, M. W., & Armitage, P. J. 2015, MNRAS, 454, 1117, doi: 10.1093/mnras/stv2070
  * Spitzer & Tomasko (1968) Spitzer, Lyman, J., & Tomasko, M. G. 1968, ApJ, 152, 971, doi: 10.1086/149610
  * Stone & Gardiner (2010) Stone, J. M., & Gardiner, T. A. 2010, ApJS, 189, 142, doi: 10.1088/0067-0049/189/1/142
  * Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
  * Sun & Bai (2021) Sun, X., & Bai, X.-N. 2021, MNRAS, 506, 1128, doi: 10.1093/mnras/stab1643
  * Suzuki et al. (2010) Suzuki, T. K., Muto, T., & Inutsuka, S.-i. 2010, ApJ, 718, 1289, doi: 10.1088/0004-637X/718/2/1289
  * Turner & Drake (2009) Turner, N. J., & Drake, J. F. 2009, ApJ, 703, 2152, doi: 10.1088/0004-637X/703/2/2152
  * Turner & Sano (2008) Turner, N. J., & Sano, T. 2008, ApJ, 679, L131, doi: 10.1086/589540
  * Umebayashi & Nakano (1981) Umebayashi, T., & Nakano, T. 1981, PASJ, 33, 617
  * Umebayashi & Nakano (2009) —. 2009, ApJ, 690, 69, doi: 10.1088/0004-637X/690/1/69
  * Villenave et al. (2022) Villenave, M., Stapelfeldt, K. R., Duchêne, G., et al. 2022, ApJ, 930, 11, doi: 10.3847/1538-4357/ac5fae
  * Walsh et al. (2012) Walsh, C., Nomura, H., Millar, T. J., & Aikawa, Y. 2012, ApJ, 747, 114, doi: 10.1088/0004-637X/747/2/114
  * Wardle (2007) Wardle, M. 2007, Ap&SS, 311, 35, doi: 10.1007/s10509-007-9575-8
  * Wardle & Salmeron (2012) Wardle, M., & Salmeron, R. 2012, MNRAS, 422, 2737, doi: 10.1111/j.1365-2966.2011.20022.x
  * Zhang et al. (2021) Zhang, K., Booth, A. S., Law, C. J., et al. 2021, ApJS, 257, 5, doi: 10.3847/1538-4365/ac1580