# Investigation of Phase Shift and Travel Time of Acoustic Waves in the Lower
Solar Atmosphere Using Multi-Height Velocities

Hirdesh Kumar Udaipur Solar Observatory, Physical Research Laboratory, Dewali,
Badi Road, Udaipur 313004 Rajasthan, India.  Brajesh Kumar  Udaipur Solar
Observatory, Physical Research Laboratory, Dewali, Badi Road, Udaipur 313004
Rajasthan, India.  Shibu K. Mathew  Udaipur Solar Observatory, Physical
Research Laboratory, Dewali, Badi Road, Udaipur 313004 Rajasthan, India.  A.
Raja Bayanna  Udaipur Solar Observatory, Physical Research Laboratory, Dewali,
Badi Road, Udaipur 313004 Rajasthan, India.  S. P. Rajaguru  Indian Institute
of Astrophysics, II Block, Koramangala, Bengaluru 560 034, India.

###### Abstract

We report and discuss phase-shift and phase travel time of low-frequency
($\nu$ $<$ 5.0 mHz) acoustic waves estimated within the photosphere and
photosphere-chromosphere interface regions, utilizing multi-height velocities
in the quiet Sun. The bisector method has been employed to estimate seven
height velocities in the photosphere within the Fe I 6173 Å line scan, while
nine height velocities are estimated from the chromospheric Ca II 8542 Å line
scan observations obtained from the Narrow Band Imager instrument installed
with the Multi-Application Solar Telescope operational at the Udaipur Solar
Observatory, India. Utilizing fast Fourier transform at each pixel over the
full field-of-view, phase shift and coherence have been estimated. The
frequency and height-dependent phase shift integrated over the regions having
an absolute line-of-sight magnetic field of less than 10 G indicates the non-
evanescent nature of low-frequency acoustic waves within the photosphere and
photosphere-chromosphere interface regions. Phase travel time estimated within
the photosphere shows non-zero values, aligning with previous simulations and
observations. Further, we report that the non-evanescent nature persists
beyond the photosphere, encompassing the photospheric-chromospheric height
range. We discuss possible factors contributing to the non-evanescent nature
of low-frequency acoustic waves. Additionally, our observations reveal a
downward propagation of high-frequency acoustic waves indicating the
refraction from higher layers in the solar atmosphere. This study contributes
valuable insights into the understanding of the complex dynamics of acoustic
waves within different lower solar atmospheric layers, shedding light on the
non-evanescent nature and downward propagation of the acoustic waves.

Solar oscillations – Helioseismology – Solar atmosphere – Quiet-Sun

##  1 Introduction

In 1962, Leighton et al. (1962) discovered velocity oscillations on the
photosphere of the Sun, depicting a dominant periodicity of the order of 5
minutes. These observed oscillations were theoretically explained by Ulrich
(1970), and Leibacher & Stein (1971) as standing acoustic waves trapped in the
acoustic cavities below the photosphere. Ulrich (1970) further suggested that
these standing waves may exist only along discrete lines in the diagnostic
diagram of horizontal wavenumber ($k_{h}$) versus frequency ($\nu$). Deubner
(1975) observationally confirmed the theoretical prediction of the ridges in
$k_{h}$ vs $\nu$ diagnostic diagram, and that finally formed the basis for the
development of helioseismology, which revealed the dynamics and internal
structures of the Sun utilizing the observed oscillations on the photosphere
(Christensen-Dalsgaard, 2002). We can infer the vast properties and dynamics
inside the Sun through global helioseismology. In contrast, local
helioseismology allows us to detect inhomogeneities in the sub-surface layers
by measuring the travel times of acoustic waves propagating inside the Sun (di
Mauro, 2008).

Meridional circulation refers to the north-south flow of solar material on the
Sun, which is a key component in providing an explanation of the solar cycle
(Babcock, 1961; Hathaway, 2015). Various methods have been utilized to measure
meridional flows (Komm et al., 1993; Hathaway, 1996; Roudier et al., 2018);
nevertheless, there is significant uncertainty in estimating the deep
circulation leading to unclear pictures (Schad et al., 2013; Zhao et al.,
2013; Jackiewicz et al., 2015; Rajaguru & Antia, 2015; Böning, 2017; Chen &
Zhao, 2017; Gizon et al., 2020). Understanding of the deep meridional
circulations in the Sun is based on helioseismic measurements. With the help
of helioseismic inversion techniques, information about the velocity
distribution at different depths below the solar surface can be derived
(Christensen-Dalsgaard (2002); Gizon & Birch (2005) and references therein).
In addition to the challenges posed by the intricate nature of the inversion
techniques, particularly when applied to active regions (Kosovichev (2012) and
references therein), the reliability of helioseismic measurements itself is
impeded by inherent uncertainties and systematic effects. Duvall (2003)
highlighted a substantial decrease in travel times with increasing latitude
for waves that had covered the same travel distances within the Sun. The
notable finding of Duvall & Hanasoge (2009) was the recognition that a
systematic effect was impacting the measurements, specifically along the
equator, during the examination of asymmetries between wave travel times in
opposite travelling directions. Utilizing the time-distance helioseismic
analysis technique, Zhao et al. (2012) found systematic centre-to-limb
variation in helioseismic travel times, and suggested that it should be taken
into consideration while determining the solar internal meridional flows,
whereas the underlying cause of this systematic effect was not well
understood. The analysis performed by Zhao et al. (2016) to examine the role
of foreshortening in the context of their role in affecting the measurements
of acoustic travel times through the time distance approach reveals that
foreshortening is not responsible for the systematic centre-to-limb effect in
the measured acoustic travel-time differences. Further, they suggested that
the formation height of the spectral line, different viewing angles for the
convective cells at different disk locations, and maybe other factors are
probably among the main causes of the centre-to-limb effect in many
helioseismic analyses (Zhao et al., 2016). The study by Chen & Zhao (2018)
revealed a remarkable variation in the centre-to-limb effect with frequency.
The effect changes its sign at a frequency near 5.4 mHz, close to the disk
centre, and reaches its maximum at around 4.0 mHz before the sign reversal.
Zhao & Chen (2020) indicated that asymmetries found in the computation of
travel time in the case of centre-to-limb variation and in the magnetized
regions are possibly due to different atmospheric heights, where oscillatory
signals are observed for estimating the acoustic travel times. Recently, Zhao
et al. (2022) utilizing the long duration multi-height velocities estimated
within Fe I 6173 Å line scan observations obtained from the Interferometric
BI-dimensional Spectrometer (IBIS; Cavallini (2006)) installed at the Dunn
Solar Telescope (DST) at Sacramento Peak, New Mexico, estimated the phase
shift and travel time in the quiet-Sun region, surrounding a small sunspot, of
the evanescent acoustic waves. It is to be noted that the helioseismic waves
observed in the photosphere are evanescent waves in the quiet Sun. The
acoustic cut-off frequency of the quiet-Sun photosphere is around 5.0 mHz
(Jiménez et al., 2011), implying that waves of frequency less than the cut-off
frequency do not propagate into the higher atmosphere even the strongest
oscillations near the 3.0 mHz. While the acoustic waves above 5.0 mHz
propagate into the higher atmosphere with increasing amplitude. It is
generally expected that the phase of evanescent waves does not change with
height in the atmosphere. The phase shift and travel time analysis performed
by Zhao et al. (2022) reveal that the acoustic waves ($<$ 5.0 mHz) carry
unexpected phase shift resulting in travel time of the order of around 1 s in
the quiet-Sun, suggesting that it may help to understand the centre-to-limb
variation of the travel time, which should be taken into account while
deriving the inferences of the deep meridional flows. They have also observed
a change in a sign around 3.0 mHz in the phase, having the presence of sparse
magnetic fields and moat flow around the sunspot; the influence of the same
was not assessed in the analysis. Waidele et al. (2023) investigated phase
shift and phase travel time of evanescent waves in the quiet-Sun, using 3D
radiative simulations and synthesized Fe I 6173 Å spectral line profiles. In
the simulation, they also found a non-zero phase shift of acoustic waves and
reported a substantial difference between the phase measured from Doppler
velocity and true velocity. Even though they have also found a non-zero phase
shift between simulated velocities for acoustic waves $\nu<$5 mHz, which is
quite unexpected. The trend of phase shift vs frequency obtained from true and
Doppler velocity is qualitatively similar to the phase shift derived from the
IBIS observations. However, the important point is that in the simulation
data, sign reversal occurs around 4 mHz, whereas it occurs near 3 mHz in the
observational data. Therefore, it is important to understand the phase shift
and phase travel time of evanescent acoustic waves within the photosphere of
the quiet-Sun, having negligible magnetic fields utilizing multi-height
velocities for a better understanding of the physical causes of non-evanescent
nature of the low-frequency acoustic waves. Examining the phase shift of the
evanescent waves from the photosphere to the chromosphere is imperative,
emphasizing the variation of phase shift into the higher solar atmosphere.

In this article, utilizing multi-height velocities estimated from the
photospheric Fe I 6173 Å and chromospheric Ca II 8542 Å lines over a quiet Sun
region in the disk centre, we present an analysis of the phase shift
measurement of the evanescent acoustic waves, to better understand the non-
evanescent nature of these waves within the photosphere, and photosphere-
chromosphere interface regions and to further rule out any possible role of
sparse magnetic fields and flows around the sunspot, which was a matter of
concern in the earlier analysis. This article is organized as follows: first,
we present observational data and reduction in the Section 2, followed by
analysis and results in the Section 3. The discussion and conclusions of our
results are highlighted in the Section 4.

##  2 The Observational Data

We have utilized data obtained from the Multi-Application Solar Telescope
(MAST: Mathew (2009); Venkatakrishnan et al. (2017)) operational at the Island
observatory in the middle of Lake Fatehasagar of Udaipur Solar Observatory
(USO), Udaipur, Rajasthan, India, and space-based Solar Dynamics Observatory
(SDO: Pesnell et al. (2012)) spacecraft to carry out the analysis. First, we
provide observational details of the ground-based telescope, followed by the
observations taken from the spacecraft.

###  2.1 MAST Data

The photospheric and chromospheric line-scan observations were obtained from
the MAST, which is a 50-cm off-axis Gregorian Solar telescope capable of
providing near-simultaneous observations of the photosphere and chromosphere
of the Sun. The maximum field-of-view (FOV) of the MAST is 6 arcmin, and the
theoretical diffraction limit of the telescope at 5000 Å wavelength is 0.252
arcsec. To obtain near-simultaneous observations, an imager optimized with two
or more wavelengths is integrated with this telescope. The Narrow Band Imager
(NBI: Raja Bayanna et al. (2014); Mathew et al. (2017)) has been developed
with the combination of two voltage-tunable lithium niobate Fabry-Perot
etalons along with a set of interference-blocking filters. These etalons are
used in tandem for photospheric observations in the Fe I 6173 Å (hereafter Fe
I) line. However, only one of the etalons is used for the chromospheric
observations in the Ca II 8542 Å (hereafter Ca II) line. Utilizing the
capability of the NBI instrument with the MAST, we observed a quiet-magnetic
network region located in the solar disc centre on January 02, 2023, from
05:39:54 UT to 07:32:52 UT duration. The Fe I line has been scanned at 35
wavelength positions at an equal sampling of 10 mÅ within 25 s, and there is a
delay of 21 s between two consecutive Fe I line scan observations. However,
the Ca II line has been scanned at 27 wavelength positions on both sides from
the centre of the line profile at [-750.0, -610.0, -510.0, -420.0, -340.0,
-270.0, -210.0, -160.0, -120.0, -85.0, -55.0, -30.0, -10.0, 0.0, 10.0, 30.0,
55.0, 85.0, 120.0, 160.0, 210.0, 270.0, 340.0, 420.0, 510.0, 610.0, 750.0] mÅ
positions within 40 s, and there is a delay of 6 s between two consecutive Ca
II line scan observations. Combining the time delay between two consecutive
line scans, we have the Fe I and Ca II observations together at a cadence of
46 s. Figure 1 represents the sample dark and flat corrected photospheric
intensity image in Fe I observations near the continuum and chromospheric line
core intensity image from Ca II line scan observations. The right panel of
Figure 1 clearly shows the significant presence of quiet regions as well as
magnetic network and internetwork regions. The left panel of Figure 3 shows
the Fe I line profile constructed from the average over full FOV and black
triangles denote the locations where the line profile has been scanned,
whereas the right panel of Figure 3 shows the Ca II line profile and the
triangles depict the locations where line profile has been scanned.

###  2.2 SDO Data

The observed FOV also contains magnetic network regions where the magnetic
field strength is of the order of kG (Solanki, 1993). To get the information
on magnetic fields in the network regions, we have utilized photospheric line-
of-sight magnetograms obtained from the Helioseismic Magnetic Imager (HMI:
Schou et al. (2012)) instrument onboard the SDO spacecraft, which observes the
photosphere of the Sun in Fe I 6173 Å line. The HMI instrument provides
photospheric observables: continuum intensity, line-depth, line width, line-
of-sight Dopplergrams, and line-of-sight magnetograms at a temporal cadence of
45 s with a plate scale of 0.504 arcsec per pixel.

Figure 1: Sample maps of dark and flat corrected photospheric continuum
intensity (left panel) obtained from the photospheric Fe I line scan and
chromospheric Ca II line core intensity (right panel), respectively, as
observed with NBI/MAST for a quiet Sun magnetic network region. Figure 2:
Sample map of photospheric line-of-sight magnetic fields obtained from the
HMI/SDO spacecraft at the beginning of the observations of the same region as
shown in Figure 1. It has been saturated between $\pm$ 50 G to enhance small-
scale magnetic features. Figure 3: Sample Fe I (left panel) and Ca II (right
panel) spectral line profiles reconstructed from the average intensity over
full FOV, where triangles denote the locations where line profiles have been
scanned, respectively. The horizontal lines connect the blue and red wings
with bisector points. Figure 4: Sample maps of running difference of
photospheric (left panel) and chromospheric (right panel) line-of-sight
Dopplergrams at PV2 and CV2 levels in the photospheric and chromospheric line
profiles, respectively, at the start time of the observations. Figure 5:
Velocity oscillations averaged over 20$\times$20 pixels in the middle of FOV
as shown in Figure 1 from the photospheric to the chromospheric layers for an
initial 20 min duration of the observations, showing the similarity between
these oscillations at various layers in the lower solar atmosphere.

##  3 Analysis and Results

In order to measure the phase shift of low-frequency acoustic (evanescent)
waves in the quiet Sun, we essentially require multi-height velocity
observations covering the photospheric-chromospheric layers of the Sun. The
bisector method (Gray, 2005) is considered to be a fast and simple tool to
infer line-of-sight velocity at multiple heights from a single line scan
observations. González Manrique et al. (2020) have utilized the bisector
method on photospheric synthetic Si I 10827 Å line profile in the quiet-Sun
and estimated multi-height velocities at different intensities. They compared
it with the simulated velocity derived at various optical depths and reported
good correlations with the velocity estimated at different heights in the
solar atmosphere. Beck & Choudhary (2020), and Beck et al. (2020) used the
bisector method to infer line-of-sight velocity at different line depths on
the Ca II 8542 Å and Hα 6563 Å line profiles in a sunspot atmosphere to
investigate the inverse Evershed flow. Grant et al. (2022) have used the
bisector method on Ca II 8542 Å line profile to estimate multi-height line-of-
sight velocity in the atmosphere of magnetic pores to study the propagation of
coherent waves. Kumar et al. (2023) utilized the bisector method on the Ca II
8542 Å line-scan observations to estimate line-of-sight velocity over a quiet-
Sun magnetic network region to study the propagation of low-frequency
magnetoacoustic waves in the aforementioned magnetized regions of the Sun.
Therefore, to obtain multi-height velocities from the photosphere to the
chromosphere, we have also applied the well-tested bisector method on
photospheric Fe I and chromospheric Ca II spectral line profiles over a quiet-
Sun magnetic network region. Before utilizing the bisector method, the line
scan observations have been subjected to dark and flat corrections. The
corrected line-scan intensity images are further co-aligned with their next
line-scan counterparts and are subjected to the pre-filter correction. In
order to do the pre-filter correction, we have normalized the Fe I and Ca II
line scan profiles with the available pre-filter profile. The pre-filter
profile is clipped between the measured line scan wavelength range by
extracting the central wavelength of the pre-filter. To extract the central
wavelength of the pre-filter, observed profiles and the profile obtained from
the BASS2000 (https://bass2000.obspm.fr/home.php) spectrum are compared by
shifting the pre-filter profile to get the same order of change at the wings.
The bisector method was used on the pre-filter corrected co-aligned line-scan
intensities, and thus we obtained maps of $\delta\lambda$ by taking an average
of 4 bisector points at a single intensity level (approximately to different
geometrical height). Following this, we have estimated $\delta\lambda$ cubes
($\delta\lambda(x,y,n)$) at seven (n = 7) heights in the Fe I line and nine (n
= 9) heights in the Ca II line at equal intensity intervals as shown in the
left and right panels of Figure 3, respectively. The red asterisk denotes the
centre of the horizontal lines connecting the blue and red wings of the
spectral line profiles. Further, we have also corrected the field-dependent
wavelength shift (centre to edge). For this purpose, we have used a diffuser
to take the line-scan observations, and these images are further subjected to
dark and flat corrections. Further, utilizing the bisector method, we obtain a
map of $\delta\lambda_{shift}$. This map was subtracted from the
$\delta\lambda$ cubes, and finally we got the maps of line-of-sight velocity
($v_{LOS}$) using the expression;
$v_{LOS}=\frac{\delta\lambda}{\lambda_{0}}c$, where $\lambda_{0}$ is the rest
wavelength, whereas $c$ is the speed of light. Norton et al. (2006) has
examined the formation height of the Fe I line profile using the VAL-C model
(Vernazza et al., 1981) and suggested that continuum forms around 16 km,
whereas line core forms around 302 km above $\tau_{c}$=1 region. Moreover, the
core of the Ca II line forms around 1300 km in the middle chromosphere,
whereas filtergrams at -600 mÅ sample the photosphere at a height of
approximately 200 km (Ortiz et al., 2014). Thus, we believe that after
estimating seven velocity pairs in the Fe I line starting from PV6 to PV0,
where PV6 samples lower height and PV0 samples upper height in the Fe I
spectral line profile, respectively, and nine bisector velocity pairs in Ca II
line ranging from CV8 to CV0; where CV8 samples lower height and CV0 samples
upper height in the Ca II spectral line profile, respectively, approximately
we have multi-height line-of-sight velocities covering from the photosphere to
the middle chromosphere of the Sun. Figure 4 shows the sample maps of line-of-
sight velocities at PV2 and CV2 pair after taking the running difference of
Dopplergrams. We have taken running differences of Dopplergrams to remove
slowly varying signals. Sample velocity oscillations are depicted in Figure 5
ranging from the lowest height in the Fe I line profile to the core of the Ca
II line profile, which has been averaged over a FOV of 20$\times$20 pixels2,
indicating similarities between oscillations at various heights in the lower
solar atmosphere.

###  3.1 Phase shift and phase travel time within the photospheric layers

In order to estimate the phase difference between two height velocity–velocity
pairs from these multi-height velocities, first of all, we have interpolated
all the velocity signals to a common temporal cadence of 45 s incorporating
the time lag present in the start time of the observations and used the
following expression to estimate the cross-spectra.

The phase shift and coherence between two evenly sampled time series ($A$ and
$B$) at each pixel over the full FOV have been estimated in the following way:

$X_{AB}(\nu)=I_{A}(\nu)\times I^{*}_{B}(\nu)$ (1)

Where I’s are the Fourier transforms, and $*$ denotes the complex conjugate.
The phase difference between two-time series is estimated from the phase of
complex cross product ($X_{AB}(\nu)$).

$\delta\phi(\nu)=tan^{-1}(Im(<X_{AB}(\nu)>)R(<X_{AB}(\nu)>)).$ (2)

Here, positive $\delta\phi(\nu)$ means that signal 1 leads 2, i.e., a wave is
propagating from lower height to upper height and vice-versa, and the
magnitude of $X_{AB}(\nu)$ is used to estimate the coherence (C), which ranges
between 0 to 1. Coherence is a measurement of the linear correlation between
two signals, where 0 indicates no correlation, and 1 means perfect
correlation. It can be estimated using the following expression:

$C(\nu)=sqrt(|<X_{AB}(\nu)>|^{2}/<|I_{A}|^{2}><|I_{B}|^{2}>).$ (3)

Figure 6: Average coherence over the full FOV estimated from two height
velocity–velocity pairs from lowest height difference to highest height
difference range, clearly depicting the decreasing coherence as height
difference ($\Delta H$) increases. Figure 7: Phase shifts measured from PV5,
PV3 and PV0 relative to PV6, displayed as a function of frequency (left panel)
integrated over the pixels having $|B_{LOS}|<$ 10 G and coherence greater than
0.5 over the full FOV. Standard deviations are displayed only for one of these
curves, while the errors for other curves are similar. Right panel: Same as
left panel, but relative time shifts are displayed. Figure 8: Relative phase
shifts within the photospheric height range (left panel) from PV6-PV5,
PV6-PV4, PV6-PV3, PV6-PV2, PV6-PV1, PV6-PV0, two height V-V pairs, displayed
as a function of height difference ($\Delta H$) for different frequency bands.
Standard deviations are displayed only for the 4.0 and 6.0 mHz frequency
bands, while the errors for other curves are within the plotted error bars.
Right panel: Same as left panel, but the relative time shifts are displayed.

Before estimating the phase difference $\delta\phi(\nu)$ between two height
V–V pairs, we have estimated the coherence using Equation 3 between various
two height V–V pairs. Figure 6 represents mean coherence over the full FOV
estimated between two height velocities as indicated within the plots. We
observe that as the height difference ($\Delta H$) between two height velocity
pairs is small, estimated coherence is very high (c.f., PV6–PV5, V–V pair),
while coherence decreases as the $\Delta H$ between two height V–V pairs
increases (c.f., PV6–CV0, V–V pair). The main focus of the analysis is on the
low-frequency acoustic waves i.e. $2-5$ mHz band, in which we found coherence
is greater than 0.5. The phase shift estimated between various two height V–V
pairs within the photospheric height range is shown in the left panel of
Figure 7. We have selected only those pixels where the absolute line-of-sight
magnetic field strength is less than 10 G to rule out the effect of magnetic
fields on the phase shift measurements. Magnetized pixels excluded with this
criterion occupy about 26$\%$ area within the full FOV being used in our
analysis. Further, only those pixels have been considered which depict
coherence greater than 0.5. Here, we observe that $\delta\phi$ is not only a
function of $\nu$ but also a function of $\Delta H$ between two height V–V
pairs (c.f., left panel of Figure 7), clearly indicating that as $\Delta H$
increases phase shift is also increasing. For the upward propagating waves, it
is usually expected that phase shift is small for neighbouring layers while
higher for increasing height difference layers. However, for the non-
evanescent waves, it should be zero or negligible. We also calculate phase
travel time from the phase difference as follows,
$t_{ph}=\frac{\delta\phi}{2\pi\nu}$, which gives information about the time
lag corresponding to phase shift. The right panel of Figure 7 shows the travel
time estimated within the quiet-Sun showing $t_{ph}$ around 1 s at 2 mHz,
nearly -2 s at 2.5 mHz 0 s around 3.2 mHz and then increases up to 2 s above 3
mHz and then again starts decreasing around 5 mHz. It is noteworthy that a
change in sign around 2 mHz and 3.2 mHz is clearly observed.

Furthermore, to better understand the height evolution of phase shift and
phase travel time, we have also examined phase shift and phase travel time as
a function of $\Delta H$ estimated at 2, 3, 4, 5, 6, and 7 mHz, which is
integrated over a certain frequency bin i.e. $\pm$ 0.25 mHz from the centre of
the mentioned frequencies. Figure 8 shows the phase difference and phase
travel time as a function of the $\Delta H$ between the photospheric height
range. It is generally expected that phase shift and phase travel time should
be zero or negligible for the frequencies ($\nu$) $<$ 5.0 mHz for the
evanescent waves. However, we note that the phase difference and phase travel
time at 3 mHz is close to 0, while at 2 mHz, the phase difference is rather
having a very small value, which gives non-zero phase travel time. On the
other hand, phase difference and phase travel time at 4, 5, 6, and 7 mHz
increases as a function of $\Delta H$ within the photosphere.

###  3.2 Phase shift and phase travel time within the photosphere-chromosphere
interface regions

Figure 9: Phase shifts measured from PV0, CV8, CV4, CV2, and CV0 relative to
PV6, displayed as a function of frequency (left panel) integrated over the
pixels having $|B_{LOS}|<$ 10 G and coherence greater than 0.5 over the full
FOV. Standard deviations are displayed only for one of these curves, while the
errors for other curves are similar. Right panel: Same as left panel, but the
relative time shifts are displayed. Figure 10: Phase shifts measured from CV8,
CV6, CV4, CV2 and CV0 relative to PV0, displayed as a function of frequency
(left panel) integrated over the pixels having $|B_{LOS}|<$ 10 G and coherence
greater than 0.5 over the full FOV. Standard deviations are displayed only for
one of these curves, while the errors for other curves are similar. Right
panel: Same as left panel, but the relative time shifts are displayed.

Utilizing the multi-height velocities estimated within the Ca II line along
with velocities derived from the Fe I line scan, we estimate phase difference
and phase travel time within the photosphere-chromosphere interface regions,
which are shown in Figure 9. Interestingly, we note that phase difference and
phase travel time estimated from the combination of photospheric and
chromospheric velocity pairs, as illustrated in the plots, indicate that
chromospheric oscillations significantly lag photospheric oscillations in
$2-5$ mHz band, and at around 5 mHz phase difference and phase travel time
start decreasing, showing a change in sign around 6 mHz, whereas within the
photosphere it occurs around 3.2 mHz. We also notice that when the $\Delta H$
between photospheric and chromospheric two height V-V pairs is small i.e.
PV6-CV8, as shown in Figure 9 by red colour, the plot depicts that broadening
is small and peaks around 4 mHz as compared to the higher height difference
photospheric and chromospheric velocity pairs i.e. PV6-CV0, which shows a
slight broadening compared to PV6-CV8 velocity pairs, because higher height
difference velocities pairs receive higher frequency contributions.
Nevertheless, we also observe that a small amount of change in phase
difference gives a significant amount of phase travel time as shown by the red
colour plot in the right panel of Figure 9. We have also estimated phase
difference and phase travel time by using photospheric PV0 velocity as the
lowest height velocity and all chromospheric velocities, shown in Figure 10.
We observe that the phase difference and phase travel time are almost similar
to the phase difference and phase travel time estimated from the PV6 and all
chromospheric velocity pairs as depicted in Figure 9 with some visual
differences.

##  4 Discussion and Conclusions

In this article, we present an analysis of the phase shift and phase travel
time of acoustic waves in the quiet-Sun located in the solar disc centre,
focusing within the photosphere and photosphere-chromosphere interface regions
of the Sun. Our investigation utilizes multi-height velocities estimated by
employing the bisector method within the Fe I and Ca II line scan observations
obtained from the NBI/MAST. Using multi-height velocities, our analysis shows
a non-zero phase shift of low-frequency acoustic waves in the quiet Sun.
Notably, a change in sign around 3 mHz is observed in the phase difference
measured within the photospheric velocity pairs, indicating both evanescent
and propagating characteristics of low-frequency ($\nu$ $<$ 5.0 mHz) acoustic
waves.

Recently, Zhao et al. (2022) utilizing high-resolution observations of the
photosphere in the Fe I spectral line obtained from the IBIS/DST, found a
similar non-zero phase shift and phase travel time of evanescent waves in the
quiet-Sun region surrounding a small sunspot. They also found a change in sign
around 3 mHz. Waidele et al. (2023) utilizing 3D simulations and from the
synthesis of Fe I line also found non-zero phase shift of evanescent waves.
However, the change in sign in the simulation was observed around 4 mHz.
Contrary to the suggestions of Chen & Zhao (2018) that the center-to-limb
effect is also responsible for the change in sign, which occurs around 5.5
mHz, our study and Zhao et al. (2022) observed the change in sign around 3
mHz. Nevertheless, in our analysis exclusion of pixels with absolute line-of-
sight magnetic fields greater than 10 G in the quiet-Sun observations yielded
results consistent with Zhao et al. (2022), reinforcing the notion that the
observed change in sign is not influenced by moat flows or sparse magnetic
fields present surrounding the sunspot as pointed out by Zhao et al. (2022).
Despite differences in spatial resolutions and modulation transfer functions
between our observations using NBI/MAST and Zhao et al. (2022) using IBIS/DST,
we obtained similar results. This finding suggests that the change in sign at
3 mHz is likely due to some physical phenomenon rather than observational or
instrumental factors. The most promising explanation of this non-zero
evanescent nature of low-frequency acoustic waves comes from the fact that the
Fe I line has asymmetry. Therefore, the velocity derived at various intensity
levels within the line itself may lead to a non-zero phase shift. Another
strong possibility is the non-adiabatic nature of the solar atmosphere, as
suggested by Zhao et al. (2022) and Waidele et al. (2023). For evanescent
waves, where phase changes in velocities may be subtle, the temperature
responses and observed intensity changes at different heights may not occur
simultaneously due to height-dependent time delays resulting from the
nonadiabatic nature of the atmosphere. Overall, the complicated interplay
between the nonadiabatic properties of the atmosphere and the spectral line’s
characteristics influences height- and frequency-dependent phase shifts in
Doppler velocities.

Additionally, our analysis incorporated chromospheric Ca II line scans,
revealing a decrease in phase shift around 5 mHz and a change in sign in the
5.5-5.8 mHz range. We speculate that the decrease in phase shift is possibly
due to the refraction and downward propagation of high-frequency acoustic
waves from the higher solar atmospheric layer. The inclusion of different
photospheric height velocities and chromospheric velocities further
demonstrates the height and frequency-dependent variation in the occurrence of
the change in sign of the phase difference of the acoustic waves in the solar
atmosphere. This finding aligns with the trends observed by Zhao et al.
(2022), where they have used line core and continuum intensity observations
obtained from the photospheric Fe I 6173 Å line from the IBIS/DST instrument.

In conclusion, our comprehensive analysis supports the existence of a non-zero
phase shift and phase travel time of low-frequency acoustic waves in the
quiet-Sun region, with the observed change in sign at 3 mHz within the
photosphere, attributed to physical phenomena rather than observational or
instrumental effects. The inclusion of chromospheric data and consideration of
height differences in velocity pairs contribute to a refined understanding of
the frequency-dependent characteristics of these waves. In future, multiline
inversion of simultaneous spectropolarimetric observations of the photosphere
and the chromosphere would be beneficial to further understand height-
dependent phase shift variations in the quiet Sun.

This work utilized data obtained from the 50 cm aperture Multi Application
Solar Telescope operated by the Udaipur Solar Observatory, Physical Research
Laboratory, Dept. of Space, Govt. of India. We are thankful to the MAST team
for providing the data used in this work. We acknowledge the use of data from
the HMI instrument onboard the Solar Dynamics Observatory spacecraft of NASA.
We are thankful to the SDO team for their open data policy. The work carried
out at Udaipur Solar Observatory, Physical Research Laboratory, Udaipur,
India, is supported by the Dept. of Space, Govt. of India. We are thankful to
the reviewer for the suggestions that improved the presentation of results as
well as the inclusion of relevant discussions in our manuscript.

## References

  * Babcock (1961) Babcock, H. W. 1961, ApJ, 133, 572, doi: 10.1086/147060
  * Beck & Choudhary (2020) Beck, C., & Choudhary, D. P. 2020, ApJ, 891, 119, doi: 10.3847/1538-4357/ab75bd
  * Beck et al. (2020) Beck, C., Choudhary, D. P., & Ranganathan, M. 2020, ApJ, 902, 30, doi: 10.3847/1538-4357/abb3bf
  * Böning (2017) Böning, V. G. A. 2017, PhD thesis, Kiepenheuer Institute for Solar Physics, Freiburg
  * Cavallini (2006) Cavallini, F. 2006, Sol. Phys., 236, 415, doi: 10.1007/s11207-006-0103-8
  * Chen & Zhao (2017) Chen, R., & Zhao, J. 2017, ApJ, 849, 144, doi: 10.3847/1538-4357/aa8eec
  * Chen & Zhao (2018) —. 2018, ApJ, 853, 161, doi: 10.3847/1538-4357/aaa3e3
  * Christensen-Dalsgaard (2002) Christensen-Dalsgaard, J. 2002, Reviews of Modern Physics, 74, 1073, doi: 10.1103/RevModPhys.74.1073
  * Deubner (1975) Deubner, F. L. 1975, A&A, 44, 371
  * di Mauro (2008) di Mauro, M. P. 2008, Astrophysics and Space Sciences Transactions, 4, 13, doi: 10.5194/astra-4-13-2008
  * Duvall & Hanasoge (2009) Duvall, T. L., J., & Hanasoge, S. M. 2009, in Astronomical Society of the Pacific Conference Series, Vol. 416, Solar-Stellar Dynamos as Revealed by Helio- and Asteroseismology: GONG 2008/SOHO 21, ed. M. Dikpati, T. Arentoft, I. González Hernández, C. Lindsey, & F. Hill, 103, doi: 10.48550/arXiv.0905.3132
  * Duvall (2003) Duvall, Thomas L., J. 2003, in ESA Special Publication, Vol. 517, GONG+ 2002. Local and Global Helioseismology: the Present and Future, ed. H. Sawaya-Lacoste, 259–262
  * Gizon & Birch (2005) Gizon, L., & Birch, A. C. 2005, Living Reviews in Solar Physics, 2, 6, doi: 10.12942/lrsp-2005-6
  * Gizon et al. (2020) Gizon, L., Cameron, R. H., Pourabdian, M., et al. 2020, Science, 368, 1469, doi: 10.1126/science.aaz7119
  * González Manrique et al. (2020) González Manrique, S. J., Quintero Noda, C., Kuckein, C., Ruiz Cobo, B., & Carlsson, M. 2020, A&A, 634, A19, doi: 10.1051/0004-6361/201937274
  * Grant et al. (2022) Grant, S. D. T., Jess, D. B., Stangalini, M., et al. 2022, ApJ, 938, 143, doi: 10.3847/1538-4357/ac91ca
  * Gray (2005) Gray, D. F. 2005, The Observation and analysis of Stellar Photospheres, 3rd edn. (Cambridge University Press)
  * Hathaway (1996) Hathaway, D. H. 1996, ApJ, 460, 1027, doi: 10.1086/177029
  * Hathaway (2015) —. 2015, Living Reviews in Solar Physics, 12, 4, doi: 10.1007/lrsp-2015-4
  * Jackiewicz et al. (2015) Jackiewicz, J., Serebryanskiy, A., & Kholikov, S. 2015, ApJ, 805, 133, doi: 10.1088/0004-637X/805/2/133
  * Jiménez et al. (2011) Jiménez, A., García, R. A., & Pallé, P. L. 2011, ApJ, 743, 99, doi: 10.1088/0004-637X/743/2/99
  * Komm et al. (1993) Komm, R. W., Howard, R. F., & Harvey, J. W. 1993, Sol. Phys., 147, 207, doi: 10.1007/BF00690713
  * Kosovichev (2012) Kosovichev, A. G. 2012, Sol. Phys., 279, 323, doi: 10.1007/s11207-012-9996-6
  * Kumar et al. (2023) Kumar, H., Kumar, B., Rajaguru, S. P., Mathew, S. K., & Bayanna, A. R. 2023, Journal of Atmospheric and Solar-Terrestrial Physics, 247, 106071, doi: 10.1016/j.jastp.2023.106071
  * Leibacher & Stein (1971) Leibacher, J. W., & Stein, R. F. 1971, Astrophys. Lett., 7, 191
  * Leighton et al. (1962) Leighton, R. B., Noyes, R. W., & Simon, G. W. 1962, ApJ, 135, 474, doi: 10.1086/147285
  * Mathew (2009) Mathew, S. K. 2009, in Astronomical Society of the Pacific Conference Series, Vol. 405, Solar Polarization 5: In Honor of Jan Stenflo, ed. S. V. Berdyugina, K. N. Nagendra, & R. Ramelli, 461
  * Mathew et al. (2017) Mathew, S. K., Bayanna, A. R., Tiwary, A. R., Bireddy, R., & Venkatakrishnan, P. 2017, Solar Phys., 292, 106, doi: 10.1007/s11207-017-1127-y
  * Norton et al. (2006) Norton, A. A., Graham, J. P., Ulrich, R. K., et al. 2006, Sol. Phys., 239, 69, doi: 10.1007/s11207-006-0279-y
  * Ortiz et al. (2014) Ortiz, A., Bellot Rubio, L. R., Hansteen, V. H., de la Cruz Rodríguez, J., & Rouppe van der Voort, L. 2014, ApJ, 781, 126, doi: 10.1088/0004-637X/781/2/126
  * Pesnell et al. (2012) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2012, Solar Phys., 275, 3, doi: 10.1007/s11207-011-9841-3
  * Raja Bayanna et al. (2014) Raja Bayanna, A., Mathew, S. K., Venkatakrishnan, P., & Srivastava, N. 2014, Solar Phys., 289, 4007, doi: 10.1007/s11207-014-0557-z
  * Rajaguru & Antia (2015) Rajaguru, S. P., & Antia, H. M. 2015, ApJ, 813, 114, doi: 10.1088/0004-637X/813/2/114
  * Roudier et al. (2018) Roudier, T., Švanda, M., Ballot, J., Malherbe, J. M., & Rieutord, M. 2018, A&A, 611, A92, doi: 10.1051/0004-6361/201732014
  * Schad et al. (2013) Schad, A., Timmer, J., & Roth, M. 2013, ApJ, 778, L38, doi: 10.1088/2041-8205/778/2/L38
  * Schou et al. (2012) Schou, J., Scherrer, P. H., Bush, R. I., et al. 2012, Solar Phys., 275, 229, doi: 10.1007/s11207-011-9842-2
  * Solanki (1993) Solanki, S. K. 1993, Space Sci. Rev., 63, 1, doi: 10.1007/BF00749277
  * Ulrich (1970) Ulrich, R. K. 1970, ApJ, 162, 993, doi: 10.1086/150731
  * Venkatakrishnan et al. (2017) Venkatakrishnan, P., Mathew, S. K., Srivastava, N., et al. 2017, Current Science, 113, 686, doi: 10.18520/cs/v113/i04/686-690
  * Vernazza et al. (1981) Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635, doi: 10.1086/190731
  * Waidele et al. (2023) Waidele, M., Zhao, J., & Kitiashvili, I. N. 2023, ApJ, 949, 99, doi: 10.3847/1538-4357/accae2
  * Zhao et al. (2013) Zhao, J., Bogart, R. S., Kosovichev, A. G., Duvall, T. L., J., & Hartlep, T. 2013, ApJ, 774, L29, doi: 10.1088/2041-8205/774/2/L29
  * Zhao & Chen (2020) Zhao, J., & Chen, R. 2020, in Astrophysics and Space Science Proceedings, Vol. 57, Dynamics of the Sun and Stars; Honoring the Life and Work of Michael J. Thompson, ed. M. J. P. F. G. Monteiro, R. A. García, J. Christensen-Dalsgaard, & S. W. McIntosh, 123–124, doi: 10.1007/978-3-030-55336-4_14
  * Zhao et al. (2012) Zhao, J., Nagashima, K., Bogart, R. S., Kosovichev, A. G., & Duvall, T. L., J. 2012, ApJ, 749, L5, doi: 10.1088/2041-8205/749/1/L5
  * Zhao et al. (2022) Zhao, J., Rajaguru, S. P., & Chen, R. 2022, ApJ, 933, 109, doi: 10.3847/1538-4357/ac722d
  * Zhao et al. (2016) Zhao, J., Stejko, A., & Chen, R. 2016, Sol. Phys., 291, 731, doi: 10.1007/s11207-016-0864-7


Date of publication xxxx 00, 0000, date of current version xxxx 00, 0000.
10.1109/ACCESS.2017.DOI

This work has been partially supported by the Swiss National Science
Foundation under grant 200020_182050.

Corresponding author: Laurent Pagnier (e-mail:
laurentpagnier@math.arizona.edu).

# Towards Model Reduction for Power System Transients with Physics-
Informed PDE

LAURENT PAGNIER1  JULIAN FRITZSCH2,3  PHILIPPE JACQUOD2,3  and Michael
Chertkov1  Program in Applied Mathematics, University of Arizona, Tucson,
United States.  School of Applied Sciences of Western Switzerland HES-SO,
Sion, Switzerland.  Department of Quantum Matter Physics, University of
Geneva, Geneva, Switzerland.

###### Abstract

This manuscript reports the first step towards building a robust and efficient
model reduction methodology to capture transient dynamics in a transmission
level electric power system. Such dynamics is normally modeled on seconds-to-
tens-of-seconds time scales by the so-called swing equations, which are
ordinary differential equations defined on a spatially discrete model of the
power grid. Following Seymlyen (1974) and Thorpe, Seyler, and Phadke (1999),
we suggest to map the swing equations onto a linear, inhomogeneous Partial
Differential Equation (PDE) of parabolic type in two space and one time
dimensions with time-independent coefficients and properly defined boundary
conditions. We illustrate our method on the synchronous transmission grid of
continental Europe. We show that, when properly coarse-grained, i.e., with the
PDE coefficients and source terms extracted from a spatial convolution
procedure of the respective discrete coefficients in the swing equations, the
resulting PDE reproduces faithfully and efficiently the original swing
dynamics. We finally discuss future extensions of this work, where the
presented PDE-based modeling will initialize a physics-informed machine
learning approach for real-time modeling, $n-1$ feasibility assessment and
transient stability analysis of power systems.

###### Index Terms:

Power system dynamics, Disturbance propagation, Electromechanical waves,
Inter-area oscillations, Physics-informed machine learning

=-15pt

##  I Introduction

This manuscript is focused on building a computationally efficient and
sufficiently accurate model describing the transient response of a
transmission level electric power system to a significant perturbation – for
example the disconnection and/or reconnection of a large generator. We
consider the dynamics of the transmission level of power systems on a
continental scale and focus on sub-minute transients on time scales ranging
from one second to few tens of seconds. We follow an approach that is standard
in power system studies and assume that the so-called swing equations [1, 2],
giving the time-evolution of the voltage angles at all nodes on the power
grid, provide a sufficiently accurate representation of the power system
dynamics within the considered spatio-temporal scales. Stated differently in
the language of modern machine learning, the spatio-temporal integration of
the swing equations provide a high-fidelity representation of the ground
truth. There are two competing aspects of the swing equations. On the one
hand, they are based on physically meaningful quantities and parameters such
as line capacities, machine inertia and damping. Accordingly they are expected
to correctly capture the physics of the system. On the other hand, integrating
these equations on a large, continental scale grid can be computationally very
expensive, even for a single run corresponding to a specific initial
condition. Obviously, it becomes even more expensive if the task is to screen
many possible initial conditions, and often prohibitively expensive when the
screening need to be repeated numerous times, testing many possible control
actions. Model reduction for this type of online applications [3] comes as a
way to strike a balance between accuracy and computational complexity. Central
to this optimization is that the transient dynamics of interest occurs over
time scales up to few tens of seconds while the goal is to numerically resolve
multiple scenarios of initial conditions and various controls faster than real
time.

How does model reduction work? In the current era of deep learning, many model
reduction techniques rely on neural networks and other tools of modern data
science and machine learning, see e.g. [4, 5, 6, 7]. The idea is to use the
ground truth model – the swing equations in our case – to produce dynamical
data, and then to train a pre-selected reduced model on these datasets to
fine-tune the parameters of the model. If the reduced model is of an
application agnostic type, as is customary in mainstream machine learning, the
scheme relies on very large datasets. However, recall that running our ground
truth model is computationally expensive. Then, if producing the needed
training datasets is not an option, can we still hope to build a reliable
reduced model? Our only hope in this case is to inject the relevant,
application-specific information – in our case information about the power
system physics – into the model reduction framework. Physics-Informed Machine
Learning (PIML) is the modern approach to resolve the model reduction
bottleneck – that is to compensate for the lack of data (typical of online
applications) by building models that are aware of the underlying physics [8],
as, e.g., expressed in terms of differential equations [9, 10, 11]. (See also
[12, 13] for discussion of the application of PIML to power systems.)

Why is Partial Differential Equation (PDE) modeling a sound option for power
system model reduction? In this manuscript, we propose a first step towards
developing PIML for general online applications and advancing model reduction
of PIML applications to power systems, as e.g. developed earlier by members of
our team [13, 14]. Similarly to [12], we take advantage of the PIML approach
and construct an online framework for simulating power system dynamics faster
than real time. We are however aiming to capture the transient dynamics in a
very large, continental-scale power system, a goal that has not been addressed
by any earlier related approach we are aware of. Accordingly, we choose to
build our reduced model on the continuous PDE approach to modeling power
system dynamics pioneered by Semlyen [15] and later extended by Thorpe,
Seyler, and Phadke [16] (see also [17]). These works were however restricted
to spatially-continuous one dimensional, time-dependent systems, i.e. with 1+1
dimensional PDE. Our PDE approach to power systems, to be presented below,
inherits all the relevant physics of the original swing Ordinary Differential
Equations (ODEs), accordingly it is 2+1 dimensional. Thus, it resolves power
grid dynamics over a spatially-continuous two-dimensional domain associated
with the power system’s geographical area of service. Approximating the swing
ODEs by a PDE may seem strange at first sight, as naively, this transition
dramatically increases the number of degrees of freedom. However, this naive
thinking is not quite right for several reasons. First, numerical solutions of
linear 2+1 dimensional PDE assume spatial regularization via a two-dimensional
grid, where the grid size can be chosen according to the desired spatial
resolution. Therefore, the number of grid points may eventually be comparable
to or even smaller than the number of nodes in the original grid. Second,
numerical operations, such as matrix inversion, can be performed much more
efficiently on a regular grid than on a complex meshed graph. Third, and most
importantly, the number of physical parameters in the original power grid
model (line capacities, machine inertia and damping coefficients) may be
reduced significantly within the PDE approach. Indeed, within the reduced
model, we want to faithfully capture only the long-wavelength components of
the swing dynamics. This justifies using a coarse-grained/filtered expression
for all the coefficients in the linear PDE, therefore representing the
coefficients via only a few long-wavelength harmonics.

Our Contribution: In this manuscript we make the first steps towards a novel
online methodology for multi-scenario testing and control based on modeling
the dynamics of a large, continental scale power system within a novel 2+1 PDE
modeling framework. We show how a properly coarse-grained PDE model faithfully
captures the power grid transient responses to disturbances of a high fidelity
model. Our approach is rather heuristic, but it is backed by physical
understanding of how perturbations propagate over the grid. Therefore, we
“prove by example” – illustrating our model reduction methodology on the
PanTaGruEl model of the synchronous grid of continental Europe introduced in
[18, 19]. Specifically, PanTaGruEl simulates power flow and swing equations
with high fidelity to produce the ground truth data. The latter in their turn
are used to infer a spatially continuous 2+1 dimensional PDE model. The
quality of the reconstruction is judged, first, by its ability to mimic power
system dynamics and, second, by a faithful reconstruction of spatially coarse-
grained and physically meaningful static – spatial distribution of line
impedances – and dynamic – spatial distribution of damping and inertia –
parameters. We conclude the manuscript with a suggestion for a path towards
using the PDE based reduced modeling framework for efficient online screening
of multiple failure scenarios on large transmission grids, faster than real
time.

##  II Problem Formulation

###  II-A Power Flow and Swing Equations (system of ODEs)

AC Power Flow (PF) equations describe steady distributions of electric power
flows over an AC power grid. The equations connect complex power injections
$\\{s_{i}\equiv p_{i}+\mathfrak{i}\,q_{i}\\}$ to complex voltages
$\\{{V}_{i}\equiv v_{i}\exp(\mathfrak{i}\,\theta_{i})\\}$, where $p_{i}$,
$q_{i}$, $v_{i}$ and $\theta_{i}$ denote the active and reactive power
injections, and the voltage magnitude and angle at node $i\in\mathcal{V}$
respectively:

$\displaystyle p_{i}$
$\displaystyle=\sum_{j}v_{i}v_{j}\Big{[}g_{ij}\cos\big{(}\theta_{i}-\theta_{j}\big{)}+b_{ij}\sin\big{(}\theta_{i}-\theta_{j}\big{)}\Big{]},\\!\\!\\!\\!$
(1a) $\displaystyle\ q_{i}$
$\displaystyle=\sum_{j}v_{i}v_{j}\Big{[}g_{ij}\sin\big{(}\theta_{i}-\theta_{j}\big{)}-b_{ij}\cos\big{(}\theta_{i}-\theta_{j}\big{)}\Big{]}.\\!\\!\\!\\!$
(1b)

Here, $b_{ij}$ and $g_{ij}$ are elements of the susceptance and conductance
matrices, see e.g. [20] for more details.

Suppose that a steady solution of the PF Eqs. (1) is perturbed, for example by
a fast disconnection and reconnection of a large generator or load. Such a
fault induces a transient voltage angle and amplitude dynamics. It is
customary to assume a time-scale separation between voltage amplitudes and
angles. On time scales ranging from sub-seconds to few tens of seconds,
voltage amplitudes remain constant, while voltage angles evolve according to
the swing equations [2, 21],

$m_{i}\ddot{\theta}_{i}+d_{i}\dot{\theta}_{i}=p_{i}-\sum_{j}v_{i}v_{j}b_{ij}(\theta_{i}-\theta_{j})\,.$
(2)

The voltage amplitudes $v_{i}$ and $v_{j}$ are considered constant, already
stabilized to the steady-state solution of Eqs. (1) and $m_{i}$ and $d_{i}$
denote the inertia and the damping (i.e., primary control) of the generators.

Figure 1: Evolution of the system parameters for the continuous PDE model with
increasing Gaussian kernel. (See text for details.)

Eq. (2) describes the relaxation dynamics of voltage angles towards a steady-
state solution, $\ddot{\theta}_{i}=\dot{\theta}_{i}=0$, corresponding to the
lossless, $g_{ij}=0$, linearized version of Eqs. (1a). Two comments are in
order here. First, the linearized approach used here is in practice quite
accurate to reproduce the transient dynamics following not too strong
perturbations - a problem called small signal stability [2]. Nevertheless, the
approach to be presented below can be extended to the nonlinear case, with
$(\theta_{i}-\theta_{j})\rightarrow\sin(\theta_{i}-\theta_{j})$ in Eq. (2).
Second, the swing equation approach is not restricted to the just discussed
case of a fast disconnection-reconnection fault, but also captures the voltage
angle dynamics following a fault which is not immediately cleared, such as the
removal of a generator or a load without reconnection. In such cases, the
final relaxed state is not balanced, i.e. $\sum_{i}p_{i}\neq 0$, and the power
mismatch is compensated by the second, damping term in Eq. (2), leading to an
AC frequency shift $\dot{\theta}_{i}=\omega_{pf}\,\forall i$, with
$\sum_{i}p_{i}=\omega_{pf}\sum_{i}d_{i}$, where $\omega_{pf}$ denotes the
post-fault synchronous frequency.

In the next paragraph, we construct a reduced model by mapping the discrete
system of ODEs (2) into a continuous PDE. Before we do that, we re-emphasize
why a model reduction is needed at all. The motivation was lucidly expressed
in Ref. [22] as follows: “The focus is on the construction of low-order models
which closely approximate the global behavior of the hybrid nonlinear system.
There is a growing recognition of the strong need for rapid and reliable
computation of the system dynamics.” Comprehensive discussions of model
reduction in a general context as well as for specific applications to slow
coherency and inter-area oscillations can further be found in Ref. [23].

###  II-B Mapping the power system to a two-dimensional continuum: the Swing
PDE Model

Consider a two-dimensional domain ${\bm{\Omega}}\subset\mathbb{R}^{2}$, with
coordinates ${\bm{r}}=(x,y)$, inside which the discrete, planar or quasi-
planar network is embedded. The boundary of the domain is denoted by
$\bm{\partial\Omega}$ and $\forall{\bm{r}}\in\bm{\partial\Omega}$,
$\bm{n}\equiv(n_{x},n_{y})$ denotes the normal vector to the boundary at
${\bm{r}}$. Imagine that the swing Eqs. (2) are derived by discretizing a PDE
describing the dynamics of a scalar field $\theta(t;{\bm{r}})$ over an
irregular mesh which corresponds to the original network. Then, following
[16], one naturally asks: what is the PDE corresponding to the swing Eqs. (2)?
We answer this question by writing the following, most general form of the
swing PDE on ${\bm{\Omega}}$:

$\displaystyle m(\bm{r})\frac{\partial^{2}}{\partial t^{2}}\theta(t;\bm{r})$
$\displaystyle+d(\bm{r})\frac{\partial}{\partial
t}\theta(t;\bm{r})=p(t;\bm{r})$
$\displaystyle+\sum\limits_{\alpha,\beta=1,2}\partial_{{r}_{\alpha}}b_{\alpha\beta}({\bm{r}})\partial_{{r}_{\beta}}\theta(t;\bm{r}),$
(3)

where ${r}_{1}=x,\ {r}_{2}=y$. One of our main task is to map the physical
parameters of (2) into the continuum as follows

$\displaystyle\forall i:\quad$
$\displaystyle\theta_{i}(t)\rightarrow\theta(t;\bm{r}),\ m_{i}\rightarrow
m(\bm{r}),\ d_{i}\rightarrow d(\bm{r})$ $\displaystyle p_{i}(t)\rightarrow
p(t;\bm{r}),\ b_{ij}\rightarrow b_{\alpha\beta}(\bm{r}),\
\forall\alpha,\beta=1,2.$ (4)

We discuss a procedure for initializing these continuous parameters in Section
III.

Figure 2: (a) Assessment of the local propagation speed as
$c({\bm{r}})=\sqrt{b({\bm{r}})/m({\bm{r}})}$. (b)-(d) Fronts of the
perturbation at incremental time intervals of $\Delta t=0.6$s, after a fault
in Greece (violet star), for inhomogeneous (red) and average parameters
(blue).

Next, the swing PDE (3) must be constrained with physically appropriate
boundary conditions. In our case, they are Neumann boundary conditions

$\displaystyle\forall t,\ \forall{\bm{r}}\in{\bm{\partial}\Omega}:\
\sum\limits_{\alpha,\beta=1,2}n_{\alpha}(\bm{r})b_{\alpha\beta}(\bm{r})\partial_{{r}_{\beta}}\theta(t;{\bm{r}})=0,$
(5)

corresponding to a vanishing normal derivative of the angle field on the
domain boundary $\bm{\partial\Omega}$. These boundary conditions directly
follow from the condition that post-perturbation frequencies in the continuous
model correspond to those in the original swing equations, i.e.
$\omega(t;\bm{r})\equiv\frac{\partial}{\partial t}\theta(t;\bm{r})=\omega_{\rm
pf}\,,\;\forall\bm{r}\in\Omega$. This condition translates into

$\begin{split}\omega_{\rm pf}\intop_{\bm{\Omega}}d(\bm{r}){{\rm
d}\bm{r}}&=\intop_{\bm{\Omega}}p(t;\bm{r}){{\rm d}\bm{r}}\\\
&+\sum\limits_{\alpha,\beta=1,2}\;\intop_{\bm{\partial\Omega}}n_{\alpha}({\bm{r}})b_{\alpha\beta}(\bm{r})\partial_{{r}_{\beta}}\theta(t;\bm{r}){{\rm
d}\bm{r}}\,.\end{split}$ (6)

This directly corresponds to the frequency shift $\sum d_{i}\omega_{\rm
pf}=\sum_{i}p_{i}$ in the swing Eqs. (2) if the second term in the right-hand
side of Eq. (6) vanishes identically, which is guaranteed by the Neumann
boundary conditions (5).

In the following we simplify our PDE model, assuming that the $b$-tensor is
diagonal $b_{12}=b_{21}=0$, accordingly, we use a shorter notation, $b_{11}\to
b_{x},\ b_{22}\to b_{y}$. We will shortly show that it can be even more
simplified to $b_{y}(\bm{r})\approx b_{x}(\bm{r})=:b(\bm{r})$.

Figure 3: Comparison of ground-truth and continuous PDE steady state
solutions. (a) One-to-one comparison of local voltage angle: for each bus in
the discrete model the nearest node in the continuous mesh is selected. The
red line indicates a perfect match. The outliers marked in orange, red and
green correspond to the points marked on the map in (b). The square markers
correspond the solution after adjusting the susceptances as outlined in the
text. (c) Continuous solution $\theta(\bm{r})$ after adjustment. (d) Discrete
solution $\theta^{\rm disc}$.

##  III Initialization of PDE Parameters and calibration

Once the general structure of the swing PDE (3) is established, we need to
initialize the PDE parameters $b_{x,y}({\bm{r}}),m({\bm{r}})$ and
$d({\bm{r}})$. This is achieved by applying a smoothing to the original
parameters defined on the discrete grid of Eq. (2). Obviously, there are many
different choices for this coarse-graining/filtering procedure. Therefore, it
is crucial to develop a validation criteria. We calibrate and validate via
post factum tests, described in the following section, where we compare the
dynamics of a fault in the original, discrete swing equations with that in the
spatially continuous model. In future work, this initialization procedure will
be complemented by a machine learning scheme performing tailored adjustments
to increase the accuracy of the model even more.

This simple smoothing procedure was proposed in Ref. [15, 16] focusing on a
1+1 dimensional PDE representation of a linear power network, where all
parameters in the 1+1 dimensional PDE were chosen to be spatially constant.
This homogeneous smoothing procedure was improved in [17], where non-uniform
parameters of the 1+1 PDE system were derived by means of a convolution with a
fixed Gaussian kernel.

We introduce a slight generalization of this Gaussian smoothing process. We
apply an Artificial Diffusion (AD) to spatial distributions of the physical
coefficients, $m({\bm{r}}),d({\bm{r}})$ or $b_{\alpha\beta}({\bm{r}}$). It
starts by assigning discrete physical quantities to the nearest nodes
discretizing the PDE (3), then they diffuse over the lattice. The longer the
diffusion is allowed, the broader the convolution kernel. The diffusion is
stopped when parameters satisfy some smoothness criterion. This generalization
is advantageous because it allows the optimal width of the Gaussian kernel to
be self-determined (no additional criteria are required).

Figure 4: Frequency response of generators in (a) Bulgaria, (b) Poland, (c)
France, and (d) Spain to a 900 MW loss of power in Greece. Data from the
continuous PDE and the ground-truth ODE models are displayed as dotted and
solid lines respectively. The map (e) shows the locations of the fault (marked
by the purple star) and of the generators considered in panels (a)-(d) (marked
by dots with the color corresponding to the respective frequency plots).

We illustrate the AD process on PanTaGruEl which consists of $3809$ buses,
$618$ generators and $4944$ lines. There are $3221$ nodes in the
discretization of our continuous model. Even though this number is not
significantly smaller than the number of buses in the discrete system, the
continuum model can be parametrized efficiently with a much fewer parameters
than what is required in the discrete model as we will show shortly. Fig. 1
illustrates the AD process on the susceptances, the damping and inertia
coefficients and the power injections. These numerical results suggest in
particular that the resulting spatial distribution of the diagonal part of the
susceptance tensor is isotropic, i.e. $b_{x}(\bm{r})\approx b_{y}(\bm{r})$.
This, combined with the assumption that the off-diagonal terms of the
$b$-matrix are much smaller than the diagonal ones (thus dropped), means that
the entire susceptance tensor may be approximated by a scalar function,
$b_{\alpha\beta}(\bm{r})\approx\delta_{\alpha\beta}b({\bm{r}})$, where
$\delta_{\alpha\beta}$ is the Kronecker symbol.

We conclude this section by showing how already this relatively simple
initialization procedure provides a useful and intuitive insight into the
system behavior. Indeed, once the continuous values of the susceptance tensor,
${\bm{b}}({\bf r})$ and of the inertia vector, ${\bm{m}}({\bm{r}})$, are
determined and validated, we can immediately use them to build a spatial map
of the electro-mechanical wave velocity,
$c({\bm{r}})=\sqrt{b({\bm{r}})/m({\bm{r}})}$, shown in Fig. (2). This figure
also illustrates how knowing the velocity, $c({\bm{r}})$, allows us to
reconstruct the dynamical spread of a localized perturbation. Notice, however
that these illustrations, even though encouraging, also suggest that one needs
to be careful in extending the approach to the grids with strong degree of
heterogeneity. The extension is certainly possible, however to achieve an
accurate representation of the actual grid dynamics by the PDE model will
require introducing more trainable parameters representing higher degree of
spatial inhomogeneity.

##  IV Numerical Validation

In this Section, we juxtapose our PDE model, with the parameters found through
the AD process explained in the preceding section, against the original swing
model considered as the ground truth. Therefore, our first validation task is
to compare the steady state solution (of the static PF equations) with the
solution of a Poisson problem associated with the static version of the PDE
model. We then compare responses, within our PDE model vs the ground truth
model, to a sufficiently large perturbation - a power outage. Finally, we show
that the number of parameters describing the PDE model can be reduced
dramatically without loss in accuracy.

###  IV-A Steady State Experiments

We start with a comparative analysis of the steady-state solutions of the PDE
and of the ground-truth model. To do that, we find the grid points closest to
the location of the buses within PanTaGruEl and terminate the AD process when
the voltage angles of the two steady-state solutions are as close to one
another as possible, $\theta_{i}^{\rm cont}\simeq\theta_{i}^{\rm disc}$.
Comparison of the two solutions is shown in Fig. 3 (a). They are clearly in a
good agreement overall, even though not without some discrepancies. The
outliers were highlighted in different colors in Fig 3 (b). We conjecture that
the discrepancies are largely due to misrepresentation of parameters in the
part of the grid with strong aspect ratio, e.g. at the Italian peninsula which
is long and narrow. Specifically, in this case the boundary conditions we set
in the continuous model may be too restrictive, effectively forcing parameters
in the part of the grid with large aspect ratios to become much smaller than
what we observe in the respective part of the discrete model.

To verify the hypothesis, we simply modify susceptances in the parts of the
grid corresponding to the outliers. Specifically, we increase susceptance
uniformly within the Italian peninsula (as the aspect ratio there is large)
and reduce it over the Iberian peninsula and Transylvania. As seen in Fig. 3,
this simple and admittedly ad-hoc adjustment was sufficient to improve the
agreement (the outlier effect was reduced significantly). We anticipate that a
more accurate and automatic tuning (via ML tools which are work in progress)
will produce even better results.

We conclude by mentioning that some of the discrepancies just discussed can be
attributed to transformers present in the ground truth model, but absent in
the PDE model. These and other strongly localized effects cannot be, properly
represented in the continuous model. However, these discrepancies are expected
to weaken at larger (spatial) scales – that is at in the coarse-grained
picture which is in the focus of our model reduction analysis. Moreover, we do
not expect the transformers to play a significant role in analysis of
transients, we are now switching our attention to.

Figure 5: Distribution of the grid parameters before and after a Fourier low
pass filter with a cut-off frequency of 30% (of the maximum frequency)
applied.

###  IV-B Dynamic Experiments

Here we discuss how well the PDE model reproduces the ground truth dynamics
observed in response to an abrupt removal of a 900 MW power plant in Greece.
Fig 4 shows comparison of the frequency response at four generators across
Europe. The agreement between the reduced and ground truth model is very good.
It is especially accurate far away from the fault. (We remind the reader that
reproducing fatefully coarser picture, and not the small scale details, is
exactly our goal.) We also report that the propagation time and the frequency
map of the inter-area oscillations are well reproduced too. Not surprisingly,
we also see that higher harmonics present close to the fault location (and
gradually disappearing as we move away from the fault) are over-estimated by
the PDE model. We expect that a more accurate – machine-learning trained – re-
parametrization of the PDE model will be able to correct the problem, however
on the expense of introducing higher-degree of parameter inhomogeneity in the
PDE model. In the end, this is a matter of a trade-off decision which a
designer of the reduced model should make — this is a trade off between
accuracy of prediction and degree of the model reduction.

###  IV-C Moving towards Model Reduction

As we mentioned earlier in the text, the number of harmonics in the parameters
of the PDE model was set to be slightly less than the number of nodes in the
original grid model. Our next step is to see if we can reduce the number of
harmonics in the PDE model coefficients even further without any significant
loss of accuracy in the spatio-temporal dynamics coarse-grained at the
resolution of 50 km or larger. In other words we now study how our PDE model
performs once we apply a low-pass filter to the coefficients.

Figure 6: Comparison of the steady state solution before (a) and after (b)
application of the Fourier low-pass filter. The two solutions are in good
general agreement, with some deviations in peripheral regions, e.g. in the
Balkans or Spain. The same color scale as in Fig. 3 is used in panel (a) and
(b). (c) Frequency response of the system to same fault as Fig. 4 (same
configuration), frequencies obtained with the filtered model are displayed
with dashed lines.

To investigate this matter, we perform an additional Fourier filtering to the
results of the AD procedure. Specifically, we choose the cut-off spatial
frequency equal to 30% of the largest spatial frequency set in the bare
(unfiltered) version of the PDE model. The results of this filtering
experiments are shown in Figs. (5,6). We observe that some compression
artifacts, similar to “ripples”, are present, but they seem to show little to
no effect on the system dynamics, see in particular Fig. 6 (c). Here we would
like to emphasize (again) that a much more accurate filtering can be achieved
with a ML approach, however it is encouraging to see that with a relatively
simple tuning of parameters we were able to achieve reproduction of the
principal features with such a good quality, even though out filtering
(parameter fitting) procedure was certainly not optimal. Let us also notice
that the coarse adjustments, mentioned in Section IV-A, are clearly visible.

##  V Critical Evaluations, Discussions and Future Work

Our main accomplishment in this manuscript is the construction and the
validation of the continuous PDE model of the swing equations. The
construction included the accurate resolution of the boundary conditions and
the development of an efficient and flexible parameter filtering procedure
based Artificial Diffusion (AD) and Fourier filtering. The validation
proceeded via static and dynamic comparisons of the continuous PDE model and
the original discrete model.

We also made a number of interesting observations which are clearly
preliminary. Comparing the computation times for dynamical simulations of the
PDE model on a regular grid of size comparable to the size of the original
discrete graph we found that the PDE model is faster by at least a factor of
ten. This is consistent with what we stated in the introduction. We also
observed that in many instances a significant filtering of the dynamical
parameters through a rather large coarse-graining scale (via wide Gaussian
kernels and Fourier filtering) does not impact the accuracy of the dynamics of
voltage frequency waves on sufficiently large scales (hundreds of kilometers)
and sufficiently long times (seconds).

Finally, this manuscript’s most important message is that the model reduction
presented here is only the starting point for an upcoming PIML methodology,
where the functional maps for $m({\bm{r}}),d({\bm{r}})$ and
$b_{\alpha\beta}({\bm{r}})$ will be modeled as Neural Networks. Specifically,
future work will focus on using the approach developed here as a warm start
for learning physical parameters of the PDE model (3). Indeed, we envision
modeling the functional maps for, $m({\bm{r}}),d({\bm{r}})$ and
$b_{\alpha\beta}({\bm{r}})$ as Neural Networks. Hence, the AD procedure is
still expected to be useful to initialize the future PIML schemes.

Planning for this work, we are aiming for keeping training process of the
Machine Learning schemes under control in terms of computation time –
consistently with the goal of making them capable of achieving the goal of
evaluating in parallel multiple perturbation scenarios in the time which is
comparable or faster than dynamic simulations of the ground truth (swing)
model.

## Acknowledgment

The authors thank Robert Ferrando and Christopher Koh for participating in the
discussions which led to the manuscript.

##  Appendix A Details on the discretization of the PDE and its numerical
integration

We use the same spatial increment $\Delta$ for x and y axes, subsequently
$\bm{r}=(i\Delta,j\Delta)$.

$\displaystyle b_{x}(\bm{r})$
$\displaystyle=\frac{b^{x}_{i,j-1}+b^{x}_{i,j}}{2}\,,$
$\displaystyle\partial_{x}b_{x}$
$\displaystyle=\frac{b^{x}_{i,j}-b^{x}_{i,j-1}}{\Delta}+O(\Delta^{3})\,,$
$\displaystyle\partial_{x}\theta$
$\displaystyle=\frac{\theta_{i+1,j}-\theta_{i-1,j}}{2\Delta}+O(\Delta^{3})\,,$
$\displaystyle\partial_{x}^{2}\theta$
$\displaystyle=\frac{\theta_{i-1,j}-2\theta_{i,j}+\theta_{i+1,j}}{\Delta^{2}}+O(\Delta^{4})\,.$

Similar expressions are obtained for the y-axis. Then discretization of the
last term in Eq.(3) becomes

$\displaystyle\partial_{x}b_{x}\partial_{x}\theta+b_{x}\partial_{x}^{2}\theta+\partial_{y}b_{y}\partial_{y}\theta+b_{y}\partial_{y}^{2}\theta\approx\big{(}b^{x}_{i,j-1}\theta_{i,j-1}$
$\displaystyle+b^{x}_{i,j}\theta_{i,j+1}+b^{y}_{i-1,j}\theta_{i-1,j}+b^{y}_{i,j}\theta_{i+1,j}-\beta\theta_{i,j}\big{)}\\!\Big{/}\\!\Delta^{2},$
(7)

where $\beta=b^{x}_{i,j-1}+b^{x}_{i,j}+b^{x}_{i-1,j}+b^{y}_{i,j}$. In order to
make our numerical scheme more efficient we vectorize (re-index) the grid, and
the field, $\theta(t;{\bm{r}})$ defined over the grid, according to
$\theta_{i,j}\rightarrow\tilde{\theta}_{k}$, where $k=N_{y}(i-1)+j$. It
results in the following re-indexing of the grid-neighbors: $i-1,j\rightarrow
k-1,\ i+1,j\rightarrow k+1,\ i,j-1\rightarrow k-N_{y},\ i,j+1\rightarrow
k+N_{y}$. This results in reformulation of the principal part of Eq. (7) in
terms of a matrix $\bm{\Xi}$ acting on the vector $\bm{\tilde{\theta}}$.
Furthermore, with the convention that inner nodes, i.e. nodes that aren’t on
the boundary layer, have a zero normal vector, $n_{x}=0$ and $n_{y}=0$, and
introducing $\eta_{\pm}(x)=\\{1\text{ if }\pm x\geq 0\,;\;0\text{
otherwise}\\}$, we rewrite $\bm{\Xi}$, therefore accounting for the Neumann
boundary conditions (5),

$\displaystyle\Xi_{kl}=$
$\displaystyle-\tilde{\beta}_{k}\delta_{k,l}+\eta_{+}(n_{x})\tilde{b}^{x}_{k-N_{y}}\delta_{k-N_{y},l}+\eta_{-}(n_{x})\tilde{b}^{x}_{k}\delta_{k+N_{y},l}$
$\displaystyle+\eta_{+}(n_{y})\tilde{b}^{y}_{k-1}\delta_{k-1,l}+\eta_{+}(n_{y})\tilde{b}^{y}_{k}\delta_{k+1,l}\,,$
(8)

where
$\tilde{\beta}_{k}=\eta_{-}(n_{x})\tilde{b}^{x}_{k}+\eta_{+}(n_{x})\tilde{b}^{x}_{k-N_{y}}+\eta_{-}(n_{y})\tilde{b}^{y}_{k}+\eta_{+}(n_{y})\tilde{b}^{y}_{k-1}$
and $\delta_{\cdot,\cdot}$ is the Kronecker product. It is important that the
method used for the numerical integration of the PDE is a finite volume
method. This class of methods is conservative. This means that there is zero
flux leakage at the boundary by construction which, in particular, guaranties
that the post-fault system frequency is indeed at the value it is expected to
be.

Finally, we use the Crank–Nicolson method [24] to integrate PDE (3). At each
time step we solve the following system of linear equations

$\bm{A}\,\left[\\!\\!\begin{array}[]{c}\bm{\tilde{\theta}}(t+\Delta t)\\\
\bm{\tilde{\omega}}(t+\Delta
t)\end{array}\\!\\!\right]=\bm{B}\,\left[\\!\\!\begin{array}[]{c}\bm{\tilde{\theta}}(t)\\\
\bm{\tilde{\omega}}(t)\end{array}\\!\\!\right]+\bm{C}\,,$ (9)

where

$\displaystyle\bm{A}$
$\displaystyle=\left[\\!\\!\begin{array}[]{cc}\mathbb{1}&-\frac{\Delta
t}{2}\mathbb{1}\\\ -\frac{\Delta
t}{2}\,\bm{M}^{-1}\,\bm{\Xi}&\mathbb{1}+\frac{\Delta
t}{2}\,\bm{\Gamma}\end{array}\\!\\!\right]\,,$ $\displaystyle\bm{B}$
$\displaystyle=\left[\\!\\!\begin{array}[]{cc}\mathbb{1}&\frac{\Delta
t}{2}\mathbb{1}\\\ \frac{\Delta
t}{2}\,\bm{M}^{-1}\,\bm{\Xi}&\mathbb{1}-\frac{\Delta
t}{2}\,\bm{\Gamma}\end{array}\\!\\!\right]\,,$ $\displaystyle\bm{C}$
$\displaystyle=\left[\\!\\!\begin{array}[]{c}\mathbb{0}\\\ \frac{\Delta
t}{2}\bm{M}^{-1}\Big{(}\bm{\tilde{p}}(t+\Delta
t)+\bm{\tilde{p}}(t)\Big{)}\end{array}\\!\\!\right]\,,$

with $\bm{M}={\rm diag}\big{(}\bm{\tilde{m}}\big{)}$ and $\bm{\Gamma}={\rm
diag}\big{(}\bm{\tilde{m}}^{-1}\bm{\tilde{d}}\,\big{)}$.

## References

  * [1] P. Kundur, _Power System Stability and Control_. MacGraw-Hill, Inc., 1993.
  * [2] J. Machowski, J. Bialek, J. R. Bumby, and J. Bumby, _Power system dynamics and stability_. John Wiley & Sons, 1997\.
  * [3] E. Hazan, “Introduction to online convex optimization,” _Foundations and Trends in Optimization_ , vol. 2, no. 3-4, pp. 157–325, 2016. [Online]. Available: http://dx.doi.org/10.1561/2400000013
  * [4] C. Ma, J. Wang, and W. E, “Model reduction with memory and the machine learning of dynamical systems,” _Communications in Computational Physics_ , vol. 25, no. 4, pp. 947–962, 2018. [Online]. Available: http://global-sci.org/intro/article_detail/cicp/12885.html
  * [5] R. Swischuk, L. Mainini, B. Peherstorfer, and K. Willcox, “Projection-based model reduction: Formulations for physics-based machine learning,” _Computers & Fluids_, vol. 179, pp. 704–717, 2019. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0045793018304250
  * [6] K. Bhattacharya, B. Hosseini, N. B. Kovachki, and A. M. Stuart, “Model Reduction And Neural Networks For Parametric PDEs,” pp. 121–157, 2021. [Online]. Available: https://smai-jcm.centre-mersenne.org/articles/10.5802/smai-jcm.74/
  * [7] W. Chen, Q. Wang, J. S. Hesthaven, and C. Zhang, “Physics-informed machine learning for reduced-order modeling of nonlinear problems,” _Journal of Computational Physics_ , vol. 446, p. 110666, 2021. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0021999121005611
  * [8] R. King, O. Hennigh, A. Mohan, and M. Chertkov, “From Deep to Physics-Informed Learning of Turbulence: Diagnostics,” _arXiv:1810.07785_ , 2018.
  * [9] I. Lagaris, A. Likas, and D. Fotiadis, “Artificial neural networks for solving ordinary and partial differential equations,” _IEEE Transactions on Neural Networks_ , vol. 9, no. 5, p. 987–1000, 1998. [Online]. Available: http://dx.doi.org/10.1109/72.712178
  * [10] M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics informed deep learning (part ii): Data-driven discovery of nonlinear partial differential equations,” _arXiv:1711.10566_ , 2017.
  * [11] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, and L. Yang, “Physics-informed machine learning,” _Nature Reviews Physics_ , vol. 3, pp. 422–440, 2021.
  * [12] G. S. Misyris, A. Venzke, and S. Chatzivasileiadis, “Physics-informed neural networks for power systems,” _arXiv:1911.03737_ , 2020.
  * [13] L. Pagnier and M. Chertkov, “Physics-informed graphical neural network for parameter & state estimations in power systems,” _arXiv:2102.06349_ , 2021\.
  * [14] ——, “Embedding power flow into machine learning for parameter and state estimation,” _arXiv:2103.14251_ , 2021.
  * [15] A. Semlyen, “Analysis of disturbance propagation in power systems based on a homogeneous dynamic model,” _IEEE Transactions on Power Apparatus and Systems_ , pp. 676–684, 1974.
  * [16] J. S. Thorp, C. E. Seyler, and A. G. Phadke, “Electromechanical wave propagation in large electric power systems,” _IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications_ , vol. 45, no. 6, pp. 614–622, 1998.
  * [17] M. Parashar, J. S. Thorp, and C. E. Seyler, “Continuum modeling of electromechanical dynamics in large-scale power systems,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 51, no. 9, pp. 1848–1858, 2004.
  * [18] L. Pagnier and P. Jacquod, “Inertia location and slow network modes determine disturbance propagation in large-scale power grids,” _PloS one_ , vol. 14, no. 3, p. e0213550, 2019.
  * [19] M. Tyloo, L. Pagnier, and P. Jacquod, “The key player problem in complex oscillator networks and electric power grids: Resistance centralities identify local vulnerabilities,” _Science advances_ , vol. 5, no. 11, p. eaaw8359, 2019.
  * [20] A. R. Bergen and V. Vittal, _Power systems analysis_. Prentice Hall, 2000.
  * [21] T. van Cutsem and C. Vournas, _Voltage Stability of Electric Power Systems_. Springer, 1998.
  * [22] P. Parrilo, S. Lall, F. Paganini, G. Verghese, B. Lesieutre, and J. Marsden, “Model reduction for analysis of cascading failures in power systems,” in _Proceedings of the 1999 American Control Conference (Cat. No. 99CH36251)_ , vol. 6, 1999, pp. 4208–4212 vol.6.
  * [23] J. H. Chow, _Power system coherency and model reduction_. Springer, 2013, vol. 84.
  * [24] J. Crank and P. Nicolson, “A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type,” _Mathematical Proceedings of the Cambridge Philosophical Society_ , vol. 43, no. 1, p. 50–67, 1947.

|  Laurent Pagnier (Member, IEEE) received the M.S. and Ph.D. degrees in
theoretical physics from EPFL, Lausanne, Switzerland, in 2014 and 2019,
respectively. He is a Visiting Assistant Professor at the University of
Arizona. His research interest includes developing new modeling and monitoring
methods for power systems. He is currently aiming at applying Machine Learning
techniques to power systems. He is particularly interested in reinforcing the
interpretablility and trustworthiness of ML methods which are paramount to
increase their acceptance and usage inside this power system community.
---|---
|  Julian Fritzsch (Graduate Student Member, IEEE) received the B.Sc. and
M.Sc. degrees in physics from the University of Würzburg, Germany, in 2017 and
2020, respectively. He is currently pursuing the Ph.D. degree with the
Department of Quantum Matter Physics, Univer- sity of Geneva, Switzerland, and
the Electrical Energy Efficiency Group, University of Applied Sciences of
Western Switzerland—HES-SO, Sion, Switzerland. His research interest includes
dynamics of power systems.
---|---
|  Philippe Jacquod (Member, IEEE) received the Diplom degree in theoretical
physics from the ETHZ, Zürich, Switzerland, in 1992, and the PhD degree in
natural sciences from the University of Neuchâtel, Switzerland, in 1997. He is
a professor with the Department of Quantum Matter Physics, University of
Geneva, Switzerland and with the engineering department, University of Applied
Sciences of Western Switzerland, Sion, Switzerland. From 2003 to 2005 he was
an assistant professor with the theoretical physics department, University of
Geneva, Switzerland and from 2005 to 2013 he was a professor with the physics
department, University of Arizona, Tucson, USA. His main research topics is in
power systems and how they evolve as the energy transition unfolds. He has
published about 100 papers in international journals, books and conference
proceedings.
---|---
|  Michael Chertkov  (Senior Member, IEEE) received the M.Sc. degree in
physics from Novosibirsk State University, Novosi- birsk, Russia, in 1990, and
the Ph.D. degree in physics from the Weizmann Institute of Science, Rehovot,
Is- rael, in 1996. After his Ph.D., he spent three years at Princeton
University as a R.H. Dicke Fellow in the Department of Physics. He joined Los
Alamos Na- tional Lab in 1999, initially as a J. R. Oppenheimer Fellow in the
Theoretical Division, and continued as a Technical Staff Member leading
projects in physics of algorithms, energy grid systems, physics, and engi-
neering informed data science and machine learning for turbulence. In 2019, he
moved to Tucson to lead Interdisciplinary Graduate Program in Applied Math-
ematics, University of Arizona, continuing to work for LANL part time. He has
published more than 200 papers. His areas of interest include mathematics and
statistics applied to physical, engineering, and data sciences. He is a Fellow
of the American Physical Society.
---|---


* Sanner and Slotine (1995) Robert M. Sanner and Jean-Jacques E. Slotine.  Stable adaptive control of robot manipulators using “neural” networks.  _Neural Computation_ , 7(4):753–790, 1995.  doi: 10.1162/neco.1995.7.4.753.  URL https://doi.org/10.1162/neco.1995.7.4.753.
  * Shi et al. (2019) Bin Shi, Simon S. Du, Weijie J. Su, and Michael I. Jordan.  Acceleration via symplectic discretization of high-resolution differential equations.  _arXiv:1902.03694_ , 2019.
  * Slotine and Coetsee (1986) J.-J.E. Slotine and J.A. Coetsee.  Adaptive sliding controller synthesis for non-linear systems.  _International Journal of Control_ , 43(6):1631–1651, 1986.  doi: https://doi.org/10.1080/00207178608933564.
  * Slotine and Li (1991) Jean-Jacques Slotine and Weiping Li.  _Applied Nonlinear Control_.  Prentice Hall, 1991.
  * Slotine (2003) Jean-Jacques E. Slotine.  Modular stability tools for distributed computation and control.  _International Journal of Adaptive Control and Signal Processing_ , 17(6):397–416, 2003.  doi: 10.1002/acs.754.
  * Slotine and Li (1987) Jean-Jacques E. Slotine and Weiping Li.  On the adaptive control of robot manipulators.  _The International Journal of Robotics Research_ , 6(3):49–59, 1987.  doi: https://doi.org/10.1177/027836498700600303.
  * Soudry et al. (2018) Daniel Soudry, Elad Hoffer, Mor Shpigel Nacson, Suriya Gunasekar, and Nathan Srebro.  The implicit bias of gradient descent on separable data.  _J. Mach. Learn. Res._ , 19(1):2822–2878, 2018\.  ISSN 1532-4435.
  * Su et al. (2016) Weijie Su, Stephen Boyd, and Emmanuel J. Candès.  A differential equation for modeling nesterov’s accelerated gradient method: Theory and insights.  _Journal of Machine Learning Research_ , 17(153):1–43, 2016.
  * Sussillo and Abbott (2009) D. Sussillo and L. F. Abbott.  Generating coherent patterns of activity from chaotic neural networks.  _Neuron_ , 63(4):544–557, 2009.
  * Tibshirani (1996) Robert Tibshirani.  Regression shrinkage and selection via the lasso.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 58(1):267–288, 1996.  doi: 10.1111/j.2517-6161.1996.tb02080.x.
  * Tyukin et al. (2007) I. Y. Tyukin, D. V. Prokhorov, and C. van Leeuwen.  Adaptation and parameter estimation in systems with unstable target dynamics and nonlinear parametrization.  _IEEE Transactions on Automatic Control_ , 52(9):1543–1559, 2007.  ISSN 2334-3303.  doi: 10.1109/TAC.2007.904448.
  * Tyukin (2003) I. Yu. Tyukin.  Adaptation algorithms in finite form for nonlinear dynamic objects.  _Automation and Remote Control_ , 64(6):951–974, 2003.  ISSN 1608-3032.  doi: https://doi.org/10.1023/A:1024141700331.
  * Tyukin (2011) Ivan Tyukin.  _Adaptation in Dynamical Systems_.  Cambridge University Press, 2011.
  * Wensing and Slotine (2018) P. M. Wensing and J. Slotine.  Cooperative adaptive control for cloud-based robotics.  In _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 6401–6408, 2018.  doi: 10.1109/ICRA.2018.8460856.
  * Wensing et al. (2018) P. M. Wensing, S. Kim, and J. E. Slotine.  Linear matrix inequalities for physically consistent inertial parameter identification: A statistical perspective on the mass distribution.  _IEEE Robotics and Automation Letters_ , 3(1):60–67, 2018.  ISSN 2377-3774.  doi: 10.1109/LRA.2017.2729659.
  * Wensing and Slotine (2018) Patrick M. Wensing and Jean-Jacques E. Slotine.  Beyond convexity – contraction and global convergence of gradient descent.  _arXiv:1806.06655_ , 2018.
  * Wibisono et al. (2016) Andre Wibisono, Ashia C. Wilson, and Michael I. Jordan.  A variational perspective on accelerated methods in optimization.  _Proceedings of the National Academy of Sciences_ , 113(47):E7351–E7358, 2016.  ISSN 0027-8424.  doi: 10.1073/pnas.1614734113.
  * Wilson et al. (2016) Ashia C. Wilson, Benjamin Recht, and Michael I. Jordan.  A lyapunov analysis of momentum methods in optimization.  _arXiv:1611.02635_ , 2016.
  * Zhang et al. (2016) Chiyuan Zhang, Samy Bengio, Moritz Hardt, Benjamin Recht, and Oriol Vinyals.  Understanding deep learning requires rethinking generalization.  _arXiv:1611.03530_ , 2016.
  * Zhang et al. (2014) Sixin Zhang, Anna Choromanska, and Yann LeCun.  Deep learning with elastic averaging sgd.  _arXiv:1412.6651_ , 2014.