# A study of sunspot 3 minute oscillations using ALMA and GST

Yi Chai  Center for Solar-Terrestrial Research
New Jersey Institute of Technology
Newark, NJ 07102, USA  Dale E. Gary  Center for Solar-Terrestrial Research
New Jersey Institute of Technology
Newark, NJ 07102, USA  Kevin P. Reardon  National Solar Observatory, Boulder,
CO, 80303  Department of Astrophysical and Planetary Sciences, University of
Colorado, Boulder, CO, 80303  Vasyl Yurchyshyn  Big Bear Solar Observatory,
Center for Solar-Terrestrial Research
New Jersey Institute of Technology
Newark, NJ 07102, USA

###### Abstract

Waves and oscillations are important solar phenomena, not only because they
can propagate and dissipate energy in the chromosphere, but also because they
carry information about the structure of the atmosphere in which they
propagate. The nature of the three-minute oscillations observed in the umbral
region of sunspots is considered to be an effect of propagation of
magnetohydrodynamic (MHD) waves upward from below the photosphere. We present
a study of sunspot oscillations and wave propagation in NOAA AR 12470 using an
approximately one-hour long data set acquired on 2015 December 17 by the
Atacama Large Millimeter/submillimeter Array (ALMA), the Goode Solar Telescope
(GST) operating at the Big Bear Solar Observatory (BBSO), the Atmospheric
Imaging Assembly (AIA) on board the Solar Dynamics Observatory (SDO), and the
Interface Region Imaging Spectrograph (IRIS). The ALMA data are unique in
providing a time-series of direct temperature measurements in the sunspot
chromosphere. The two-second cadence of ALMA images allows us to well resolve
the three-minute periods typical of sunspot oscillations in the chromosphere.
Fourier analysis is applied to ALMA Band 3 ($\sim$100 GHz, $\sim$3 mm) and GST
H$\alpha$ data sets to obtain power spectra as well as oscillation phase
information. We analysed properties of the wave propagation by combining
multiple wavelengths that probe physical parameters of solar atmosphere at
different heights. We find that the ALMA temperature fluctuations are
consistent with that expected for a propagating acoustic wave, with a slight
asymmetry indicating non-linear steepening.

Chromosphere, Sunspot, Oscillation

††journal: ApJ

##  1 Introduction

Sunspot oscillations are a frequently studied wave phenomenon in the solar
atmosphere. These oscillations are directly connected with the propagation of
solar magnetohydrodynamic (MHD) waves and may play a role in coronal and
chromospheric heating. Moreover, they can also serve as a probe of the
structure of the solar atmosphere. Oscillations in sunspot umbrae (e.g.,
Beckers & Tallant, 1969; Beckers & Schultz, 1972; Giovanelli, 1972) were
reported shortly after the discovery of p-mode oscillations. There are two
types of characteristic oscillations in sunspots, the photospheric 5-minute
oscillations driven by the p-modes and 3-minute oscillations representing the
resonant mode of the sunspot itself (Thomas, 1981; Bogdan & Judge, 2006;
Khomenko & Collados, 2015).

Studies of sunspot oscillations have been carried out by multiple instruments
over the years, including the Solar Optical Telescope onboard Hinode
(Nagashima et al., 2007), the Nobeyama Radioheliograph (Reznikova et al.,
2012), the Atmospheric Imaging Assembly on board SDO (Reznikova et al., 2012;
Sych & Nakariakov, 2014), the Interface Region Imaging Spectrograph (Tian et
al., 2014; Yurchyshyn et al., 2015), the Goode Solar Telescope (GST) operating
at Big Bear Solar Observatory (BBSO) (Maurya et al., 2013; Yurchyshyn et al.,
2015; Su et al., 2016) and other ground-based telescopes (e.g., Centeno et
al., 2006; Löhner-Böttcher et al., 2016; Anan et al., 2019; Felipe et al.,
2020).

The sunspot oscillations are a specific manifestation of the significant
amount of acoustic energy generated through oscillatory motion in near-surface
layers. In the upper photosphere, most of the upward propagating waves at
lower frequencies (below 5 mHz, or periods longer than 3-4 minutes) are
reflected back downward while higher frequency waves continue to propagate
upward into the chromosphere. A variety of observations have revealed how
energy is deposited in the upper atmosphere due to such wave motions. For
example, a study by Reardon et al. (2008) using Ca II (854.2 nm) line
observations from the Interferometric Bidimensional Spectrometer (IBIS)
revealed the presence of a power-law distribution of significant oscillatory
power up to 25 mHz, suggestive of the widespread presence of turbulence from
shock dissipation in the chromosphere.

Waves also carry information about the structure of the atmosphere in which
they propagate. Since the sound speed is a function of temperature in the
chromosphere, pressure disturbances with frequencies above the acoustic cut-
off frequency can provide direct information about the atmosphere at different
heights and the accompanying temperature and velocity perturbations.

Over many years, numerous efforts have being made in studying the particular
oscillatory signals seen in sunspot umbrae. Using H$\alpha$ filters on the 12
inch solar telescope in Culgoora, Giovanelli (1972) was the first to measure
sunspot oscillations in H$\alpha$ velocity. Later on, more detail has been
revealed using H$\alpha$ line wing data. Phillis (1975) showed the intensity
variation in red and blue wings of H$\alpha$ at $\pm 0.3$Å that demonstrate
the line-of-sight velocity field as an oscillation on a upward flow. Uexkuell
et al. (1983) compared the oscillations measured from several lines (Ca II, Na
D1 and D2, Ni I, H$\alpha$) using power spectra analysis and studied the phase
relation between them. With the high cadence filtergrams from Universal
Birefringent Filter on the Dunn Solar Telescope (DST), Christopoulou et al.
(2000, 2001) studied the relation between running penumbra waves and umbral
oscillations in different layers of the solar atmosphere marked by different
lines from H$\alpha$ line center and wings. Tziotziou et al. (2007) studied
multiple umbral flashes in one sunspot using Ca II and H$\alpha$ intensity
images and revealed the coexistence of more than one oscillating mode,
suggesting different physical conditions existing in the umbra.

However, studies of wave propagation in the chromosphere are not
straightforward. For example, Mein & Mein (1976) and Fleck & Deubner (1989)
found no phase difference between the Doppler shifts of the Ca II infrared
triplet lines at 854.2 nm and 849.8 nm, which are formed at different heights,
and used this as evidence of a non-propagating component in the chromospheric
wave field. However, this conclusion was questioned by Skartlien et al.
(1994), whose simulations of propagating disturbances predicted little or no
phase difference between these lines, similar to the above observations, even
though propagating shocks were present. This shows that in interpreting
observations it is necessary to consider the radiative transfer effects and
the atmospheric structure over the entire range of heights contributing to the
observed emissions.

Previous observations have shown the capability of detecting oscillations in
sunspots using radio instruments such as the Nobeyama Radio Observatory
(Shibasaki, 2001) and the Very Large Array (Nindos et al., 2002). Loukitcheva
et al. (2004) has demonstrated the feasibility of measuring chromospheric
oscillations in the mm range. Such observations have the distinct advantage of
directly probing the plasma temperature at different heights in the solar
chromosphere depending on the observing frequency.

Early research Leibacher & Stein (1981); Zhugzhda & Dzhalilov (1982); Zhugzhda
(1984) on the oscillations in the chromosphere suggested that the oscillatory
modes are produced by a resonant cavity. However, in more recent studies
Felipe et al. (2010); Chae & Goode (2015), a different scenario has been
proposed that the 3 minutes oscillation is a direct indication of vertically
propagating waves that travel through the gravitationally stratified medium.
Despite the growing knowledge and complexity of theoretical models, there
remain several open questions related with the fundamental physical mechanisms
in the sunspot region.

The need for better understanding of the fine structure of the 3 min
oscillation and its time evolution in sunspots has intensified with the
development of better observating tools. Among modern observatories, the
Atacama Large Millimeter/submillimeter Array (ALMA) opens up a new era of
solar radio observation due to its high spatial and temporal resolution and
image quality (e.g. Nindos et al., 2018; Molnar et al., 2019; Nindos et al.,
2021). When combined with other cutting-edge instruments, such as BBSO/GST,
SDO/AIA and IRIS, ALMA can provide unique electron temperature diagnostics
that clarify the behavior of the solar chromosphere to propagating waves (see
Patsourakos et al., 2020; Nindos et al., 2021, for a discussion of $p$-mode
oscillations seen with ALMA).

The aim of this paper is to present an analysis of 3-minute oscillations in a
sunspot umbra using data gathered by multiple instruments and demonstrate the
presence of such oscillatory motion in ALMA Band 3 data. With the help of
high-cadence and high-resolution ALMA data as well as the linear dependence
between ALMA image intensity and brightness temperature of the plasma
(Wedemeyer et al., 2016), we report some of the first observations of
spatially resolved temperature oscillations in a sunspot umbra.

##  2 Observations

We analysed the oscillation signal measured in the western part of active
region NOAA 12470 during one-hour observing period (18:42 UT – 19:48 UT on
2015 December 17), which was part of an ALMA solar commissioning campaign
whose data are in the public domain. The ALMA dataset we used was provided by
the ALMA observatory as part of its Science Verification data release.
Previous publications based on this ALMA data release are Bastian et al.
(2017), Iwai et al. (2017), Loukitcheva et al. (2017), and Shimojo et al.
(2017). Only Shimojo et al. (2017) used data from the same day (2015 December
17) as considered here, and that study was focused on a small brightening
rather than the sunspot itself. Our observations were conducted using the band
3 receiver (centered at 100 GHz) in a single-pointing (snap-shot) mode. Due to
the necessity of phase calibration, blocks of on-target observations were run
for approximately 615 s at a cadence of 2 s, followed by a 220-s off-target
phase calibrator scan, which resulted in total of five solar scans to be
analyzed. In this paper, we will use designations $t1$, $t2$, $t3$, $t4$, and
$t5$ to represent each solar scan arranged in time increase order. The ALMA
antenna configuration consisted of total 31 antennas with twenty two 12-m
antennas and nine 7-m antennas. The longest baseline of the antenna
configuration was $\sim$300 m, which resulted in a synthesized beam of 6.3″
$\times$ 2.3″ (Shimojo et al., 2017).

We followed the general calibration method of radio interferometry using the
Common Astronomy Software Applications (CASA) package (Petry & CASA
Development Team, 2012) to create images at 2 s cadence, and found that the
result is heavily influenced by phase errors produced by the temporal and
spatially variable water content of the Earth’s atmosphere, which can be
severe enough to render sunspot oscillations undetectable. In order to counter
this effect, we carried out a self-calibration process, which uses the solar
signal itself to correct for antenna-based phase and amplitude errors, similar
to a method by Shimojo et al. (2017).

The scheme of self-calibration can be described as follows: a) split the data
that have been corrected using standard calibration into 10 min sub-datasets
$t1$-$t5$ according to the scan index. b) For each sub-dataset, make an image
from the entire 10-min period using the CLEAN task in CASA, to use as a model
for self-calibration. The model is stored in the dateset for later use.
Atmospheric fluctuations will largely average out in such a long integration.
c) Determine phase corrections for each 2-s time period relative to the model
given by step b using the gaincal task in CASA, and generate a self-
calibration table. d) Apply the phase corrections to each image in the sub-
dataset using the applycal task in CASA and generate an intermediately
corrected sub-dataset. e) Iteratively repeat steps b-d, substituting in each
intermediate sub-dataset, until no further correction is needed, as determined
by checking the calibration table for reduction of the residuals.

In our case, we applied this phase self-calibration cycle four times,
resulting in final residual phase fluctuations of around 10 degrees standard
deviation. We applied the scheme to each of the measured 10-min sub-datasets
to complete the self-calibration procedure. The result shows a clear
improvement in removing the fluctuations caused by the Earth’s atmosphere and
allows the regular variations due to sunspot oscillations to be measured.

The smallest spatial scale of the antenna configuration is roughly 70″, which
is not sufficient to measure the solar background brightness temperature over
the band 3 field of view (FOV; $\sim$80″ in Figure 1). To overcome this
drawback, the single-dish observations of the full Sun were taken using three
other antennas to measure the overall emission close in time to when ALMA was
observing the target (White et al., 2017). The information from both
interferometer and single-dish antennas then were combined using the
“feathering” process in CASA to obtain the true temperature scale on the full-
range of spatial scales. The left panel of Figure 1 shows the end result of
this full process, which has a pixel scale of 0.3″ and FOV as defined by the
interferometer images. The field of view is truncated at the 20% level of the
peak response of the primary (single-dish) beam and has been corrected for the
declining antenna response, which becomes more extreme near the edges.

Due to the fact that ALMA is using the equatorial coordinate system (RA-Dec)
as the output framework for its data sets, as part of the image processing we
performed a coordinate transformation to heliocentric coordinates, to align
the ALMA images with those from other solar instruments. In this procedure,
the center of the field of view (FOV) of each solar image in one scan was
transformed based on the heliocentric coordinates at the mid-time of the scan
to align all of the images to a fixed time. A time dependent rotation was
applied to take into account the rotation between these two frames, leaving
the images aligned with solar north.

Figure 1: (a) ALMA brightness temperature image with contours from HMI
continuum image showing the umbral and penumbral boundary of sunspot AR12470
at 19:02 UT. The box marked with a dashed line shows the FOV of the BBSO GST
H$\alpha$ image in the right panel. (b) GST H$\alpha$ near blue wing (-0.4 Å)
image with the same HMI umbral contour.

To compare the ALMA images with those at other wavelengths, we obtained the UV
and EUV intensity data cube from the Atmospheric Imaging Assembly (AIA) on
board the Solar Dynamics Observatory (SDO) (Lemen et al., 2012) in the
bandpass of 171 Å (Fe IX: the quiet corona) and 304Å (He II: the chromosphere
and transition region). The time cadence of these two data sets is 12 s and
the pixel resolution is 0.6″. The duration of the AIA data we used covers the
full ALMA solar observation from 18:42 UT to 19:48 UT. Images of both
wavelengths were extracted from the full-disk data to fit the FOV of ALMA and
the center coordinates of these images were shifted to correct for solar
rotation. The continuum intensity data product from by the Helioseismic and
Magnetic Imager (HMI) on SDO was used for alignment of between the SDO and GST
images.

The Goode Solar Telescope (GST) operating at BBSO (Cao et al., 2010) observed
this AR nearly simultaneously with ALMA (from 18:50 UT to 20:59 UT) at
H$\alpha$ line-center and off-band wavelengths ($\pm 0.4$ Å and $\pm 0.8$ Å)
using the Visible Imaging Spectrograph (VIS) instrument. The spatial scale of
the GST H$\alpha$ images is 0.034″ per pixel, the VIS bandwidth is 0.08 Å, and
the time cadence was 40 s. Speckle reconstruction is normally performed by
taking 100 frames in rapid succession and applying the algorithm to obtain a
single sharper image. However, this algorithm can produce occasional glitches
in low-light regions such as sunspot umbrae, which are the focus of our study.
Since high spatial resolution is not needed for large-scale features of umbral
oscillations, the unreconstructed H$\alpha$ images were used. The GST
observations cover the entire ALMA observing period $t2$–$t5$ but are missing
the first few minutes during $t1$ set. Since our interest is on the joint
ALMA-GST coverage, we mainly focus on analysing the $t2$–$t5$ ALMA sub-
datasets.

We also used the slit-jaw images in Mg II k band (2796 Å) from the Interface
Region Imaging Spectrograph (IRIS) which is formed in the chromosphere (De
Pontieu et al., 2014) as a comparison to the ALMA temperature oscillation.
IRIS started observations during $t3$ and suffered some cosmic-ray noise
during the first 10 minutes, which can be seen later in the time-distance
plots.

##  3 Analysis

To visualize the oscillatory motion in the sunspot in various wavelengths, we
show in Figure 2 time-distance plots along a line (Figure 2h) that crosses
both umbra and penumbra. The double solid lines in Figure 2h mark the width of
the slit over which the data are summed, with an ALMA map for reference,
overlaid with HMI contours that mark the boundaries of the umbra and penumbra.
The time-distance plots were constructed for each wavelength after co-
alignment, as shown in the other panels of Figure 2. The $y$ coordinate in
each plot corresponds to the $y$ coordinate of the center of the slit in
Figure 2h. Note that the plots are marked with the $y$ coordinate, not
distance along the slit. Multiply the $y$ coordinate by $\sqrt{2}$ for this
roughly $45^{\circ}$ slit to convert to distance.

The horizontal black lines in each panel mark the extent of the umbra. The
ALMA data for $t1$-$t5$ are in Figure 2a, with gaps indicating times of no
data due to calibrations. GST H$\alpha$ center-line and blue-wing data are
shown on the right (Figs 2b,d,f). IRIS 2796 Å slit-jaw data are shown in
Figure 2c. AIA 304 Å and 171 Å data are shown in Figs 2e,f. Non-uniform
chevron-shaped features (Kobanov et al., 2006) that tilt in the same direction
can be seen in Figs 2b–g, suggesting an outward propagation towards the
penumbra. The 3-minute oscillatory motion in these wavelengths also extends
somehow into the penumbra (above or below the solid black lines in Figure
2a–g), but with some merging and increased separation to match the penumbral
5-minute oscillation period. This signature of propagation seems to be far
less apparent in the ALMA data, and the oscillatory signature is only
prominent in the umbra.

Figure 2: Time-distance plots based on the full duration of ALMA band 3
observation, 18:42:33-19:48:31 UT. An average of pixels between the two
diagonal lines marked in panel (h) is used to constructed each time-distance
plot. From (a) to (g): ALMA Band 3; BBSO GST H$\alpha$ center line; IRIS
2796Å; BBSO GST H$\alpha$ near blue wing (-0.4 Å); AIA 304Å; BBSO GST
H$\alpha$ far blue wing (-0.8 Å); AIA 171Å. The horizontal solid lines at
y-offsets 219″ and 234″ mark the umbral boundaries, and the gaps between each
ALMA band 3 time block are represented with white gaps. (h): ALMA image of
sunspot NOAA 12470 with the boundaries of the slit and the contours from HMI
image overlaid to show the umbral boundary.

Our discovery of umbral oscillation in ALMA band 3 is groundbreaking in many
respects. ALMA data provide a unique probe of temperature changes in the solar
atmosphere in a relatively thin layer where the opacity is near unity
(Loukitcheva et al., 2015). It is of interest to compare these oscillations in
temperature with intensity oscillations provided by other instruments that
probe different layers of solar atmosphere. Figure 3 compares the light curves
obtained from a cut at $x=-251.0$″, $y=222.5{\arcsec}$ in the time-distance
plot (Figure 2) for multiple instruments. The location of the cut was chosen
based on the wave patterns in all wavelengths shown in Figure 2 so that it
covers the region of strongest oscillations. The data have been arbitrarily
scaled and offset for better comparison. Though separated by gaps, we can
still see some correspondence between light curves of different wavelengths,
which enable the calculation of phase relations that will be discussed later.

Figure 3: Light curves from an average over pixels between the two slit in
Figure 2h at $x=-251.0$″, $y=222.5$″. The point along the slit was chosen due
to its good representation of the wave pattern. The name of instruments and
observation band are on the left hand label, and the ALMA temperature range is
shown with the scale on the right.

In order to study the oscillation power, we applied the Fourier power spectrum
for each pixel in both the ALMA map and GST image time series. Although the
temporal resolution of ALMA (2 s) and GST ($\sim$40 s) data are quite
different, they both contain sufficient samples to resolve 3-minute
oscillations and therefore can be used to measure the 3-minute oscillation
power. The power spectrum was calculated as follows. For each ALMA (10-minute)
solar scan, a map cube structure was made according to the method described
above so the time series of each pixel can be easily obtained. The GST data
were selected based on start and end time of each ALMA solar scan and made
into map cube as well. Inspired by the method used in Ireland et al. (2015),
the time series of each pixel in the ALMA and GST data cubes were obtained and
normalized by $(I(t))-\langle I(t)\rangle)/\langle I(t)\rangle$, where $I(t)$
is the varying brightness temperature for ALMA or varying intensity for GST,
and $\langle I(t)\rangle$ is the mean over each 10-minute ALMA scan.

A Hanning window was applied to the normalized time series to minimize edge
effects. The full spatial resolution was retained during the process in order
to reach the finest possible detail on the distribution of frequency-related
physical contents as a function of location. We then applied the Fourier
transform to the processed time series to form the power spectrum from each
pixel. Figure 4 shows examples from three pixels in both the ALMA and GST
H$\alpha$ near blue wing (-0.4 Å) data cubes for time $t2$, one in the region
of umbral oscillations, one in the penumbra, and one in the surrounding quiet
Sun. No GST power spectrum could be obtained for the quiet Sun due to limited
spatial coverage (see Figure 1b).

Figure 4: Three examples showing the power spectrum in the sunspot umbra,
penumbra and quiet Sun regions. (a) ALMA image for locating the sample points,
marked with plus signs. The inner white contour is the boundary of the sunspot
umbra while the outer contour is the boundary for the sunspot penumbra. (b)
Power spectrum of GST H$\alpha$ near blue wing (-0.4 Å) data from red point in
the sunspot umbra showing a strong peak in the 3-minute oscillation band
marked with dot-dashed line. (c) The same as (b) for yellow point in the
sunspot penumbra. (d) Power spectrum of ALMA data from the red point in the
sunspot umbra. (e) (f) the same as (d) for sunspot penumbra and quiet Sun
region.

The differences in shape among these power spectra are consistent with what
was seen in the ALMA time-distance plot in Figure 2a. The 3-minute band of the
power spectra, which is the two closest bins to $0.0056$ Hz marked with red
dot-dashed lines, is especially strong for the umbral power spectrum. From
Figure 4, we see that the peak of the umbral power spectrum is one to two
magnitudes larger than the others in both data sets. The umbral spectrum of
ALMA data also shows some evidence for harmonics near $0.015$ Hz, which
strongly suggests that the oscillation deviates from a sinusoid, as we
demonstrate later. The frequency distribution of power is similar for the
penumbra and quiet Sun.

To characterize the power spectrum spatial distribution, for each pixel we
calculated the integrated power in the 3-minute band marked with the red dash-
dot lines in Figure 4. This procedure was carried out for each 10-minute ALMA
scan $t2$ to $t5$. The result for $t2$ is shown in Figure 5, where the
3-minute oscillation power is shown with red contours. Blue contours outline
the umbral and penumbral boundary as seen in the HMI continuum image in panel
(a).

We find that the strong 3-minute oscillation power falls solely inside the
umbral boundary. The oscillations seem to be strongest in the western part of
the umbra, while in the eastern region separated by a small light bridge
structure, the 3-minute oscillations are present but weaker. This conclusion
may be affected by the declining ALMA sensitivity at the edge of the primary
beam, however.

To further compare the power spectra in different regions, we calculated the
mean power for each frequency bin using all the pixels in the umbra, penumbra
and quiet Sun to create the averaged power spectrum of ALMA t2 data for each
region. A combination of two power-law functions, $af^{b}+cf^{d}$, was used in
fitting the power-law tail starting from the frequency bin immediately after
the 3-minute range. The power-law index of the lower frequency term is shown
in Figure 5 panel (b) - (d). The averaged penumbral power spectrum was over-
plotted using red lines on the umbral and quiet spectra for inter-comparison.
The steeper slope shown in panel (b) by the blue dashed line reveals the
influence of the high 3-minute power in umbral region, while the power-law
indexes from the penumbra and quiet Sun are less steep. It is interesting that
the umbral and penumbral high-frequency tail are essentially the same above
0.03 Hz (panel b), while the quiet Sun spectrum in panel (d) shows a notable
excess near 0.03 Hz that contributes to a shallower and extended low-frequency
slope. A red-noise spectrum above about 0.07 Hz is the same for all three
regions, and may reflect residual fluctuations due to the Earth’s atmosphere
that have not been removed by the self-calibration procedure. Note that
although the 5-minute oscillations cannot be well-resolved in the frequency
domain, there is no sign of a peak near 0.0033 Hz in the penumbral spectrum of
Figure 5c, which further confirms what we learned from the time-distance
plots—that the penumbral 5-minute oscillations seen in H$\alpha$ (Figure 2b,d)
are not apparent in the ALMA band 3 data (Figure 2a).

Figure 5: Comparison of ALMA power spectra averaged over umbral, penumbral,
and quiet Sun regions. (a) HMI continuum image at 19:02 UT, overlaid with red
contours showing the 3-min oscillation power and blue contours outlining the
umbra and penumbra from the HMI continuum image. (b) Averaged umbral power
spectrum with blue dashed line showing a region well fit by a power-law. (c)
Same as (b) for the penumbral region. (d) the same as (b) for quiet Sun. The
averaged penumbra power spectrum in panel (c) is over-plotted in red in both
(b) and (d).

Figure 6 shows the distribution of 3-min oscillation power in the $t2$ to $t5$
time ranges (columns) and all bands from ALMA and GST (rows). The rows
represent the spatial evolution of the 3-min power for the wavelength band
marked in the first column. The sunspot umbra is outlined in blue (see Figure
5) while the black contours with gray shading show the 3-min power maps of
each data set. For a given band, the contour levels were set to be the same
for $t2$ to $t5$, which are 10, 30, 50, 70 and 90 percent of the maximum value
of oscillation power over the four power maps.

Figure 6: Spatial distribution and relative power level of 3-min oscillations
in various bands. For this comparison the H$\alpha$ data have been temporally
interpolated and spatially smoothed to match the ALMA resolution. The contours
in ALMA and the H$\alpha$ sub-band images represent 10, 30, 50, 70 and 90
percent of the peak value of the 3-min oscillation power in each wavelength.
Red and green plus symbols mark regions of high and low power, respectively,
whose light curves are plotted in Figure 7.

From Figure 6 we see that the shape of GST H$\alpha$ far red wing (+0.8Å)
3-min power map contours is less organized and weaker than the others. While
the other contours are generally constrained to lie well within the umbra, the
power in the H$\alpha$ far-red wing mainly lies along the umbral boundary. The
relative amplitude of the +0.8Å power is also 5-50 times lower than in other
bands, and hence we do not include it in our subsequent analysis.

For the remaining bands, the 3-min power was in all cases strongest and most
concentrated in $t2$ and reduced to 20-40% of that value in $t3$ to $t5$.
Although it is difficult to tell from the shape and location of these
contours, which suffer from the limited cadence of the H$\alpha$ data, the
overall pattern is consistent with the often seen appearance of new centers of
oscillation power plus the lateral spreading of wave energy from previous
centers (Yurchyshyn et al., 2020).

To determine the phase relationships among ALMA band 3 and GST H$\alpha$ sub-
bands, we plot in Figure 7 the light curves from these wavelengths taken in a
region of high oscillation power (red plus symbol in Figure 6) and for
comparison a neighboring region of low power (green plus symbol in Figure 6).
The selected locations for the high power sample are different from $t2$ to
$t5$, but within each given time block, the same position is used for every
wavelength. The ALMA 2-s cadence data are plotted directly (black curve),
while the lower-cadence H$\alpha$ measurements are plotted in different
symbols and colors, with the solid curves through the points showing the cubic
spline interpolation to the ALMA cadence.

Figure 7: A comparison of light curves from pixels of strong oscillation (left
column) and those from a relatively quiet region (right column) from locations
marked with red and green plus sign in Figure 6. Different bands of ALMA and
GST H$\alpha$ data, which follow the same order top-to-bottom as in Figure 6,
are plotted with different colors as shown by the labels adjacent to the
upper-left panel. The time axis is marked in seconds, sequentially from the
start time of the ALMA observations at 18:42:33UT. The dashed lines in $t3$
left column mark the peak times of these light curves for one particular
oscillation, as discussed in the text. The ALMA maximum and minimum brightness
temperatures are also shown in each left-hand panel, showing a fluctuation of
$\pm 250$-500 K.

To better reveal the phase relationships among the different bands shown in
Figure 7, we show five vertical dashed lines in the left column T3 panel that
indicate the peak positions of the five light curves for one particular
oscillation. The colors match the colors of light curves for the corresponding
wavelengths. These dashed lines reveal a general sequence of the phase
relationship, with H$\alpha$ far blue wing (-0.8 Å or -0.08 nm) peaking first,
then in order of peak time the H$\alpha$ near blue wing (-0.4 Å), ALMA band 3
(100 GHz or 3.0 mm), H$\alpha$ center line and H$\alpha$ near red wing (+0.4
Å). We calculate the delay times based on these peak positions, relative to
the H$\alpha$ far blue wing to be 14 s, 32 s, 64 s and 82 s. This timing
pattern is followed remarkably faithfully in the other individual oscillations
in the left column, with similar time lags. Since the locations of the sample
point in the umbra change in the different periods $t2$ to $t5$, the
relatively stable phase shift between these bands suggests that it is a
characteristic of the propagating waves that should be accounted for in any
physical interpretation.

We also compared the light curves from the strong 3-min oscillation region in
the left column with the ones from a nearby, relatively quiet umbral region
(green plus symbols in each panel of Figure 6) in the right column for
$t2$-$t5$. Unlike the case for the peak oscillation region, where the location
is different for each time window, the quiet region selection is fixed for all
of the time windows. For each time window, the amplitude range was scaled to
match the left column to show the difference between these two regions. One
can see that for $t2$, $t4$, and $t5$ ALMA temperature oscillations are not
apparent when compared to the high-power region, although velocity
fluctuations in the H$\alpha$ bands still persist at a reduced level. In $t3$
the ALMA oscillations do appear in the quiet region, again at a reduced level.
The significant differences in oscillatory power across the umbra, both in
temperature and velocity diagnostics, indicate rapid spatial and temporal
modification of the wave propagation characteristics.

Figure 8: Epoch analysis of the relative phase of the oscillation pattern in
different bands using the light curves from the strong oscillation region in
the left column of Figure 7. Left panel: analysis for each ALMA time window
separately. Each curve is offset vertically for clarity. The ALMA band 3 light
curves (black) were moved to third place to best reflect the phase order. Each
wave in the light curves is distinguished using different line styles (solid,
dashed, and dotted line for the first, second, and third complete wave
respectively). Right panel: an overlay of all waves from the left column.

In order to study the general phase relationships among all bands in the
strong 3-min oscillation region, we performed an epoch analysis where
individual oscillations are overlaid in terms of their phase relative to the
peak of the H$\alpha$ far blue wing ($-$0.8 Å). This was done for individual
time windows $t2$ to $t5$ in the left column of Figure 8, with solid, dashed,
and dotted lines for the first, second, and third oscillations in each time
window, respectively, and for all time windows in the right column. The
oscillations in different wavelengths were offset vertically for clarity, and
the ALMA light curve was moved to the middle position to reflect its ordering
of peak time midway between H$\alpha$ near blue wing ($-$0.4 Å, cyan) and
H$\alpha$ center line (yellow).

The left column of Figure 8 shows that the oscillations in $t2$ are highly
coherent for both the ALMA and GST data, reflecting the stronger oscillation
power during that time window. The peak of the first wave (solid line) in each
band in $t2$ is marked with a symbol to show the progression in phase of the
peaks similar to the vertical dashed lines in Figure 7. The right panel of
Figure 8, overlaying the curves from all the time blocks, shows that this
general pattern persists for all of the oscillations although there is higher
variability in some bands. The H$\alpha$ near blue wing and ALMA oscillations
are particularly coherent for most of the oscillations, and there is a marked
tendency for the drop in ALMA brightness (corresponding to temperature) to be
steeper than its rise.

##  4 Discussion

In order to put the oscillations in temperature detected by ALMA into context
with the wavelength-dependent brightness variations in the H$\alpha$ line, and
thus account for the phase shifts in Figures 7 and 8, a quantitative solar
atmospheric model that exemplifies the sunspot structure is needed. Inspired
by Molnar et al. (2019), we used the RH code (Rybicki & Hummer, 1991, 1992;
Uitenbroek, 2000) to synthesize the chromospheric radiation in 1D, which uses
the Solar Irradiance Physical Modeling (SRPM) as an input for the calculation
(Fontenla et al., 2011). Due to the estimated formation height of ALMA band 3
(Loukitcheva et al., 2015), we used a sunspot umbral model (model S) that
focuses on the photosphere and chromosphere (model index 1006). The RH code
solves combined equations of statistical equilibrium and radiation transfer
for multi-level atoms and molecules under certain input parameters. The
simulation was carried out with non-LTE setup for 6-levels of hydrogen atom
(levels 0 to 5, plus the continuum), while other atomic species were treated
with the LTE assumption. In addition to the default H$\alpha$ output
wavelength grid, we also used the RH code to produce the radiative emission in
the wavelength range from 2.6 mm to 3.4 mm corresponding to the ALMA band 3
observations.

Figure 9: Spectral synthesis results obtained with the RH code. Left panel:
Intensity contribution function (gray shade) vs. height for the H$\alpha$ line
from the FAL S model that simulates a sunspot umbra, overlaid with the
emergent line profile (violet line, scale on the right axis). The wavelengths
are shown relative to H$\alpha$ line center, with four vertical dashed lines
to indicate the H$\alpha$ sub-bands used for the study. The line colors follow
the same scheme as in Figures 7 and 8. Right panel: Contribution function for
the emergent intensity for ALMA band 3 wavelength range from the same FAL S
model, overlaid with the emergent intensity profile (violet line).

The results from the RH simulation of the umbra model are shown in Figure 9
for both wavelength ranges. The contribution functions are illustrated in the
underlying gray scale map while the emergent intensities are plotted in violet
lines. The H$\alpha$ center line position was set to zero, and four dashed
vertical lines indicate the sub-bands used in our study. From this result, we
can see that the contribution function for H$\alpha$ extends over a wide
height range as seen on the left panel. The major part of the center line
intensity comes from a narrow layer in the middle chromosphere ($\sim$ 1100
km) while the continuum is mainly formed in the photosphere. In the ALMA band
3 wavelength range, the formation height seen in the right panel is extended
between 600 km and 1000 km in the chromosphere. Combined with the phase
pattern indicated by the dashed vertical lines in Figure 7, we can conclude
that the H$\alpha$ line wing signatures are most likely dominated by the
velocity of the waves, as observed by Phillis (1975).

Figure 10: A series of interpolations of the H$\alpha$ line profile in phase
order based on the simulation result from RH code as well as the relative
intensities of the four H$\alpha$ sub-bands. The dashed line represents the
original emergent line profile (violet line in Figure 9 left panel), while
each colored line represents Gaussian fitting to the four H$\alpha$ line
intensities at 10 times (marked with 10 vertical dashed lines in Figure 8).
The color bar indicates the time bin. For these fits we kept the continuum and
line-core intensities fixed and fit only line shift and line width. The error
bars plotted for the first Gaussian fit were calculated based on the standard
deviation of the light curves in the nearby quiet region for each of the four
H$\alpha$ sub-bands.

To better understand the influence of the oscillation on the H$\alpha$ line
profile, we use the first complete oscillation in the $t2$ light curves to fit
dynamic H$\alpha$ line profiles. Because no photometric calibration is
available for the GST data, we instead use the normalized values from each
H$\alpha$ sub-band light curve in the top-left panel of Figure 8 to calculate
the relative variations at these four sub-bands at any given time. In this
process, we assume that the line profile at H$\alpha$ off-center bands can not
exceed the continuum, so the light curves are obtained relative to the maximum
values at these wavelengths. For the line center the relative variation is
based on the mean value of its light curve. To anchor the continuum for the
fits, two far off-center wavelengths (+/- 3 Å) were added to represent the
continuum value, we then applied Gaussian fitting to these six points to form
an interpolated line profile for a given time. To further stabilize the fits,
only line width and line center position were allowed to change during the
Gaussian fitting.

The result is shown in Figure 10. Each color in the color bar represents one
of the vertical dashed lines in the top left panel of Figure 8, with numbers
showing the order in time. The original simulation result from the RH code is
also plotted with a dashed line in the figure, and is used to provide the
intensity scale shown on the vertical axis. From the color-coded time bins in
the figure we can form a schematic picture of how the H$\alpha$ line profile
changed with time/oscillation phase—it starts near the original dashed line
profile, shifts to longer wavelengths (blue curves), then rapidly shifts to
shorter wavelengths (green), and finally moves back to the center (red).

Figure 11: Top panel: A comparison between an ALMA light curve (first
oscillation of $t2$) and full width at half maximum (FWHM) of the Gaussian
fitting in Figure 10. Bottom panel: The Doppler velocity corresponding to the
H$\alpha$ center line shift in Figure 10. The upward Doppler velocity is
treated as positive in the figure.

Although these fits are approximate, the behavior of the resulting line full-
width half-maximum (FWHM) and line shift (Doppler velocity) shown in Figure 11
are suggestive. In the upper panel of Figure 11, we plot the temperature
oscillation from ALMA alongside the line width information generated by the
Gaussian fitting. Note that the line width varies in phase with the
temperature variations measured by ALMA, with the greatest width corresponding
to the highest temperature. This is consistent with the finding by Molnar et
al. (2019) that ALMA temperature is correlated with H$\alpha$ line width.

The lower panel of Figure 11 shows the Doppler velocity. The positive Doppler
shift corresponds to the upward propagating waves. The maximum downward
Doppler velocity amplitude is 3.4 km/s while the maximum upward velocity is
3.9 km/s. The asymmetric Doppler velocity profile reveals that the plasma
moves downward slowly with a steady acceleration and then suddenly changes its
direction and moves upward with a much larger acceleration, creating a steep
rise in the velocity diagram.

The clear signature of a sawtooth pattern in the Doppler velocity for this
particular oscillation, which is the largest amplitude one in our
observations, has to be taken with some caution due to the relatively low
cadence ($\sim$40 s) and small number of wavelengths (4) on which it is based.
Even so, we find that both the Doppler velocity magnitude and trend in Figure
11 agree well with those reported by Tziotziou et al. (2006) and Chae &
Litvinenko (2017). This asymmetrical Doppler velocity trend continues in
observations of transition region lines (Tian et al., 2014), with a higher
velocity of order 10 km s-1. Such asymmetries are evidence of steepening or
shock behavior in the atmosphere of sunspot umbrae, which was discovered in
the Doppler shift of Ca II (Lites, 1984) as well as He I (Lites, 1986). In any
case, the asymmetry in the ALMA temperature behavior (slow rise in temperature
followed by a faster fall) is clearly established by our observations with 2-s
cadence, as seen in each of the individual oscillations in Figure 7 and in the
right panel of Figure 8.

The ALMA temperature variation—its magnitude, phase relative to H$\alpha$
Doppler velocity, and its asymmetric shape—all provide new information on the
phenomenon of umbral 3-min oscillations. We have investigated whether these
observations can be explained by current theories of acoustic wave propagation
by performing a simulation from the model of Chae & Goode (2015). We compare
the results of our simulations from their model with the observations, as
shown in Figure 12. In brief outline, the steps we took to produce Figure 12b
are (i) implementing their model in Python, (ii) using it to calculate the
atmospheric parameter perturbations (velocity $v$, number density perturbation
$dn_{e}$, and temperature perturbation $dT$) as a function of height and time,
(iii) modifying the FAL umbral model by adding these perturbations as
fractional changes at multiple time steps, and (iv) running the RH code at
these multiple time steps to simulate the H$\alpha$ line profile and ALMA
brightness temperature. The good agreement of the phase differences between
the observations (Figure 12a) and simulations (Figure 12b) suggest that the
ALMA temperature variations are mainly a consequence of an upward propagating
acoustic wave in the stratified medium of the sunspot umbra.

Figure 12: Comparison of (a) the observed H-alpha intensity in red and blue
wing and ALMA temperature variation with (b) the same quantities from the
hydrodynamic acoustic model of Chae & Goode (2015). Vertical dashed lines mark
positions in phase similar to those in the top left panel of Figure 8

Similar to Figure 10, we plotted the temperature profile based on the
simulated emergent intensity in band 3 wavelengths from the RH code at the 10
phases marked in Figure 12b. The colors from blue to red indicate increasing
phase (or time), as shown in the color bar. The temperature profile starts at
the lowest temperature, rises to the peak, and then falls back. The rapid line
shift that happens between time bins 4 to 7 in Figure 10 corresponds to the
temperature peak in Figure 13, but the earlier peak in the observations
(Figure 12a) relative to the linear model (Figure 12b) suggests a tendency for
the chromosphere to be heated more strongly at this time than described by the
linear theory. The absolute temperature and temperature fluctuation predicted
by the RH code, about 8700 K and $\pm 1100$ K, respectively at 3 mm, are both
somewhat higher than the observations (Figure 7), which are around 8000 K and
$\pm 500$ K. See Loukitcheva et al. (2017) for a detailed comparison of the
brightness temperatures with various chromospheric models for this same
sunspot observed one day earlier.

Figure 13: A series of ALMA temperature in phase order based on the simulation
result from RH code. The dashed line represents the original temperature
profile calculated based on the emergent line profile (purple line) in Figure
8 right panel. The color bar indicates the time bin. The temperature gradient
reflects the fact that shorter wavelength of ALMA band3 sees deeper in the
chromosphere.

##  5 Conclusions

Our study describes the first observations of sunspot oscillations detected by
ALMA. The ALMA observations provide unique measurements of the chromospheric
temperature variation tied to electron temperature via free-free emission, and
so is relatively insensitive to atomic abundances and local non-equilibrium
conditions.

We applied Fourier analysis to the ALMA band 3 data to create spatial maps of
the 3-min power amplitude. The 3-min power maps correlate well to similar maps
for all H$\alpha$ sub-bands except the far red wing (+0.8 Å) where the 3-min
oscillation signal is weak. From the light curves in Figures 7 and 8, we
discovered a relatively constant phase relationship among H$\alpha$ sub-band
intensity and ALMA temperature as well as evidence for a sawtooth pattern in
the ALMA temperature variation, albeit less abrupt than the velocity
transition seen in H$\alpha$. The phase relationship among the different
H$\alpha$ bands sampling different parts of the line profile is consistent
with a shift of the line center during the oscillation as well as a slight
variation in line width. The ALMA temperature and the H$\alpha$ line width
peak at around the same time, but the ALMA temperature shows a more rapid fall
after the peak time.

We used RH code to simulate the H$\alpha$ line profile and the expected ALMA
emission based on the FAL S sunspot atmosphere model. The results are
generally in agreement with observations although the ALMA band 3 brightness
temperature predicted by the FAL S model (8700 K) is some 9% greater than
observed (8000 K). In order to understand the phase relationship of the
observed ALMA temperature variation, we used the linear theory of Chae & Goode
(2015) for the propagation of a hydrodynamic disturbance in a simplified solar
atmosphere to obtain the velocity, density, and pressure perturbations as a
function of height, then applied those perturbations to the FAL S model.
Calculations of the H$\alpha$ and ALMA emission using repeated runs of the RH
code for different times during one oscillation of such a disturbance gave
results that are in reasonable agreement with the model. In particular, the
calculated relative phase relationships among the H$\alpha$ near blue wing
intensity, ALMA band 3 temperature, and H$\alpha$ near red wing intensity are
close to those observed, showing that the ALMA temperature variations are
largely due to atmospheric heating expected for a propagating disturbance.

The symmetric temperature variation from the model is due to the linear theory
used for wave propagation in Chae & Goode (2015) which cannot reproduce the
slight sawtooth pattern observed in both the ALMA band 3 temperatures and
H$\alpha$ Doppler velocity, therefore a study of the non-linear solution for
the wave function is much needed. Based on previous research (Litvinenko &
Chae, 2017; Chae & Litvinenko, 2017), the non-linear effect can provide not
only a physics scenario of the wavefront steepening and shock formation that
is related with the sawtooth pattern, but also an explanation for the 2nd
harmonic seen in the umbra power spectrum in Figure 4. We also note that the
RH code gives an H$\alpha$ line width for the FAL S model that is
significantly narrower than commonly observed (Fricke & Elsässer, 1965;
Wallace et al., 2005). We have done a preliminary study of other models (e.g.
the FAL penumbral model R and the Maltby model (Maltby et al., 1986), which
can correct the H$\alpha$ line-width discrepancy but results in estimated ALMA
temperatures (and heights of formation) that are much too high. Here we simply
note the discrepancy—a more thorough investigation will be addressed in a
future paper now in preparation.

These new observations show the power of ALMA mm-wave observations for
providing entirely new diagnostics of the solar atmosphere. Because of the low
solar activity since these science-verification observations were made, no
additional observations of sunspot umbrae have been possible to date. However,
advances in ALMA capabilities for solar observations, such as higher spatial
resolution, polarization, and additional frequency bands that provide
diagnostics at both lower and higher heights, all offer great promise for
investigations in the near future as solar activity returns. Such observations
should motivate more sophisticated modeling that includes the effect of the
hydrodynamic waves on ionization states as well as the effect of the strong
sunspot magnetic fields.

This paper makes use of the following ALMA data:
ADS/JAO.ALMA#2011.0.000020.SV. ALMA is a partnership of ESO (representing its
member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST
and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the
Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and
NAOJ. BBSO operation is supported by NJIT and US NSF AGS-1821294 grant. GST
operation is partly supported by the Korea Astronomy and Space Science
Institute, the Seoul National University, and the Key Laboratory of Solar
Activities of Chinese Academy of Sciences (CAS) and the Operation, Maintenance
and Upgrading Fund of CAS for Astronomical Telescopes and Facility
Instruments. The National Radio Astronomy Observatory is a facility of the
National Science Foundation operated under cooperative agreement by Associated
Universities, Inc. YC gratefully acknowledges support from the NRAO Student
Observing Support program. The data used in this paper can be obtained from
the following sources: The IRIS data can be obtained from
https://iris.lmsal.com/data.html. GST data can be obtained from the BBSO data
request form http://www.bbso.njit.edu/$\sim$vayur/nst_requests2/. SDO data
(AIA and HMI) can be obtained from http://jsoc.stanford.edu/.

## References

  * Anan et al. (2019) Anan, T., Schad, T. A., Jaeggli, S. A., et al. 2019, ApJ, 882, 161. doi:10.3847/1538-4357/ab357f
  * Bastian et al. (2017) Bastian, T. S., Chintzoglou, G., De Pontieu, B., et al. 2017, ApJ, 845, L19. doi:10.3847/2041-8213/aa844c
  * Beckers & Tallant (1969) Beckers, J. M. & Tallant, P. E. 1969, Sol. Phys., 7, 351
  * Beckers & Schultz (1972) Beckers, J. M. & Schultz, R. B. 1972, Sol. Phys., 27, 61
  * Bogdan & Judge (2006) Bogdan, T. J. & Judge, P. G. 2006, Philosophical Transactions of the Royal Society of London Series A, 364, 313. doi:10.1098/rsta.2005.1701
  * Cao et al. (2010) Cao, W., Gorceix, N., Coulter, R., et al. 2010, Astronomische Nachrichten, 331, 636. doi:10.1002/asna.201011390
  * Centeno et al. (2006) Centeno, R., Collados, M., & Trujillo Bueno, J. 2006, Solar Polarization 4, 358, 465
  * Chae & Goode (2015) Chae, J. & Goode, P. R. 2015, ApJ, 808, 118. doi:10.1088/0004-637X/808/2/118
  * Chae & Litvinenko (2017) Chae, J. & Litvinenko, Y. E. 2017, ApJ, 844, 129. doi:10.3847/1538-4357/aa7be9
  * Christopoulou et al. (2000) Christopoulou, E. B., Georgakilas, A. A., & Koutchmy, S. 2000, A&A, 354, 305
  * Christopoulou et al. (2001) Christopoulou, E. B., Georgakilas, A. A., & Koutchmy, S. 2001, A&A, 375, 617. doi:10.1051/0004-6361:20010887
  * De Pontieu et al. (2014) De Pontieu, B., Title, A. M., Lemen, J. R., et al. 2014, Sol. Phys., 289, 2733. doi:10.1007/s11207-014-0485-y
  * Felipe et al. (2010) Felipe, T., Khomenko, E., & Collados, M. 2010, ApJ, 719, 357. doi:10.1088/0004-637X/719/1/357
  * Felipe et al. (2020) Felipe, T., Kuckein, C., González Manrique, S. J., et al. 2020, ApJ, 900, L29. doi:10.3847/2041-8213/abb1a5
  * Fleck & Deubner (1989) Fleck, B. & Deubner, F.-L. 1989, A&A, 224, 245
  * Fontenla et al. (2011) Fontenla, J. M., Harder, J., Livingston, W., et al. 2011, Journal of Geophysical Research (Atmospheres), 116, D20108. doi:10.1029/2011JD016032
  * Fricke & Elsässer (1965) Fricke, K. & Elsässer, H. 1965, Zeitschrift für Astrophysik, 63, 35
  * Giovanelli (1972) Giovanelli, R. G. 1972, Sol. Phys., 27, 71
  * Ireland et al. (2015) Ireland, J., McAteer, R. T. J., & Inglis, A. R. 2015, ApJ, 798, 1. doi:10.1088/0004-637X/798/1/1
  * Iwai et al. (2017) Iwai, K., Loukitcheva, M., Shimojo, M., et al. 2017, ApJ, 841, L20. doi:10.3847/2041-8213/aa71b5
  * Khomenko & Collados (2015) Khomenko, E. & Collados, M. 2015, Living Reviews in Solar Physics, 12, 6. doi:10.1007/lrsp-2015-6
  * Kobanov et al. (2006) Kobanov, N. I., Kolobov, D. Y., & Makarchik, D. V. 2006, Sol. Phys., 238, 231. doi:10.1007/s11207-006-0160-z
  * Leibacher & Stein (1981) Leibacher, J. W. & Stein, R. F. 1981, NASA Special Publication, 263
  * Lemen et al. (2012) Lemen, J. R., Title, A. M., Akin, D. J., et al. 2012, Sol. Phys., 275, 17. doi:10.1007/s11207-011-9776-8
  * Lites (1984) Lites, B. W. 1984, ApJ, 277, 874. doi:10.1086/161758
  * Lites (1986) Lites, B. W. 1986, ApJ, 301, 1005. doi:10.1086/163964
  * Litvinenko & Chae (2017) Litvinenko, Y. E. & Chae, J. 2017, A&A, 599, A15. doi:10.1051/0004-6361/201629568
  * Löhner-Böttcher et al. (2016) Löhner-Böttcher, J., Bello González, N., & Schmidt, W. 2016, Astronomische Nachrichten, 337, 1040. doi:10.1002/asna.201612430
  * Loukitcheva et al. (2004) Loukitcheva, M., Solanki, S. K., Carlsson, M., et al. 2004, A&A, 419, 747. doi:10.1051/0004-6361:20034159
  * Loukitcheva et al. (2015) Loukitcheva, M., Solanki, S. K., Carlsson, M., et al. 2015, A&A, 575, A15
  * Loukitcheva et al. (2017) Loukitcheva, M. A., Iwai, K., Solanki, S. K., et al. 2017, ApJ, 850, 35. doi:10.3847/1538-4357/aa91cc
  * Maltby et al. (1986) Maltby, P., Avrett, E. H., Carlsson, M., et al. 1986, ApJ, 306, 284. doi:10.1086/164342
  * Maurya et al. (2013) Maurya, R. A., Chae, J., Park, H., et al. 2013, Sol. Phys., 288, 73. doi:10.1007/s11207-013-0286-8
  * Mein & Mein (1976) Mein, N. & Mein, P. 1976, Sol. Phys., 49, 231. doi:10.1007/BF00162447
  * Molnar et al. (2019) Molnar, M. E., Reardon, K. P., Chai, Y., et al. 2019, ApJ, 881, 99. doi:10.3847/1538-4357/ab2ba3
  * Nagashima et al. (2007) Nagashima, K., Sekii, T., Kosovichev, A. G., et al. 2007, PASJ, 59, S631. doi:10.1093/pasj/59.sp3.S631
  * Nindos et al. (2002) Nindos, A., Alissandrakis, C. E., Gelfreikh, G. B., et al. 2002, A&A, 386, 658. doi:10.1051/0004-6361:20020252
  * Nindos et al. (2018) Nindos, A., Alissandrakis, C. E., Bastian, T. S., et al. 2018, A&A, 619, L6. doi:10.1051/0004-6361/201834113
  * Nindos et al. (2021) Nindos, A., Patsourakos, S., Alissandrakis, C. E., et al. 2021, A&A, 652, A92. doi:10.1051/0004-6361/202141241
  * Patsourakos et al. (2020) Patsourakos, S., Alissandrakis, C. E., Nindos, A., et al. 2020, A&A, 634, A86. doi:10.1051/0004-6361/201936618
  * Petry & CASA Development Team (2012) Petry, D. & CASA Development Team 2012, Astronomical Data Analysis Software and Systems XXI, 461, 849
  * Phillis (1975) Phillis, G. L. 1975, Sol. Phys., 41, 71. doi:10.1007/BF00152957
  * Reardon et al. (2008) Reardon, K. P., Lepreti, F., Carbone, V., et al. 2008, ApJ, 683, L207. doi:10.1086/591790
  * Reznikova et al. (2012) Reznikova, V. E., Shibasaki, K., Sych, R. A., et al. 2012, ApJ, 746, 119. doi:10.1088/0004-637X/746/2/119
  * Rybicki & Hummer (1991) Rybicki, G. B. & Hummer, D. G. 1991, A&A, 245, 171
  * Rybicki & Hummer (1992) Rybicki, G. B. & Hummer, D. G. 1992, A&A, 262, 209
  * Shibasaki (2001) Shibasaki, K. 2001, ApJ, 550, 1113. doi:10.1086/319820
  * Shimojo et al. (2017) Shimojo, M., Hudson, H. S., White, S. M., et al. 2017, ApJ, 841, L5. doi:10.3847/2041-8213/aa70e3
  * Skartlien et al. (1994) Skartlien, R., Carlsson, M., & Stein, R. F. 1994, Chromospheric Dynamics, 79
  * Su et al. (2016) Su, J. T., Ji, K. F., Cao, W., et al. 2016, ApJ, 817, 117. doi:10.3847/0004-637X/817/2/117
  * Sych & Nakariakov (2014) Sych, R. & Nakariakov, V. M. 2014, A&A, 569, A72. doi:10.1051/0004-6361/201424049
  * Thomas (1981) Thomas, J. H. 1981, The Physics of Sunspots, 345
  * Tian et al. (2014) Tian, H., DeLuca, E., Reeves, K. K., et al. 2014, ApJ, 786, 137. doi:10.1088/0004-637X/786/2/137
  * Tziotziou et al. (2006) Tziotziou, K., Tsiropoula, G., Mein, N., et al. 2006, A&A, 456, 689. doi:10.1051/0004-6361:20064997
  * Tziotziou et al. (2007) Tziotziou, K., Tsiropoula, G., Mein, N., et al. 2007, A&A, 463, 1153. doi:10.1051/0004-6361:20066412
  * Uexkuell et al. (1983) Uexkuell, M. V., Kneer, F., & Mattig, W. 1983, A&A, 123, 263
  * Uitenbroek (2000) Uitenbroek, H. 2000, ApJ, 536, 481. doi:10.1086/308933
  * Wallace et al. (2005) Wallace, L., Hinkle, K., & Livingston, W. C. 2005, An atlas of sunspot umbral spectra in the visible from 15,000 to 25,500 cm-1 (3920 to 6664 Å) , by Wallace, Lloyd.; Hinkle, Kenneth.; Livingston, W. C. [Tucson, AZ : National Solar Observatory, 2005]
  * Wedemeyer et al. (2016) Wedemeyer, S., Bastian, T., Brajša, R., et al. 2016, Space Sci. Rev., 200, 1
  * White et al. (2017) White, S. M., Iwai, K., Phillips, N. M., et al. 2017, Sol. Phys., 292, 88
  * Yurchyshyn et al. (2015) Yurchyshyn, V., Abramenko, V., & Kilcik, A. 2015, ApJ, 798, 136. doi:10.1088/0004-637X/798/2/136
  * Yurchyshyn et al. (2020) Yurchyshyn, V., Kilcik, A., Şahin, S., et al. 2020, ApJ, 896, 150. doi:10.3847/1538-4357/ab91b8
  * Zhugzhda & Dzhalilov (1982) Zhugzhda, I. D. & Dzhalilov, N. S. 1982, A&A, 112, 16
  * Zhugzhda (1984) Zhugzhda, I. D. 1984, MNRAS, 207, 731. doi:10.1093/mnras/207.4.731


Axel Maas

# Strong isospin breaking in $S\\!p\\!(4)$ gauge theory

Fabian Zierler

###### Abstract

A hidden strongly-interacting sector is a possible candidate for dark matter.
One scenario consistent with astrophysical constraints is a QCD-like theory
with gauge group Sp(4) and two fundamental non-degenerate fermions. We report
exploratory investigations of the light spectrum and decay constants for this
theory using lattice calculations.

##  1 Introduction

The nature of dark matter (DM) is one the most challenging problems in current
particle physics. Recently a class of strongly-interacting, QCD-like gauge
theories with a number of massive fermions has gained some interest as dark
matter candidates [1, 2, 3, 4]. In these theories dark matter is made up of
the pseudo-Goldstones bosons associated with the breaking of the approximate
chiral symmetry and dark matter can self-annihilate through $3\to 2$ processes
in the early universe. In addition non-degenerate fermions are interesting
because they can establish a mass hierarchy in the dark sector. As an example
of such a sector we study here $Sp(4)$ gauge theory with two fundamental non-
degenerate Dirac fermions in isolation. To the best of our knowledge this is
the first time that strong isospin breaking has been studied in a non-$SU(N)$
gauge theory on the lattice. Here we present first, exploratory results for
the masses and decay constant of the pseudoscalar and vector iso-nonsinglet
mesons.

##  2 Setup, global symmetries and observables

The fundamental representation in any symplectic group is pseudo-real. This
entails relations between particles and anti-particles and their Weyl
components. Therefore this theory has a global $U(4)_{I}$ isospin symmetry in
case of vanishing masses. It is broken by the axial anomaly down to
$SU(4)_{I}$. Spontaneous chiral symmetry breaking and/or degenerate fermions
break this symmetry further to $Sp(4)_{I}$. In case of non-degenerate masses
the symmetry is broken down further to $SU(2)_{u}\times SU(2)_{d}$ [5, 6].
These relate the Weyl components of the two dark fermion flavours, the (dark)
up and the (dark) down.

As this theory confines the physical degrees of freedom are hadrons. Due to
the relation between particles and anti-particles and the even number of
colors all hadrons are bosons, and made from an even number of fermions. In
particular, in addition to the usual quark-antiquark mesons additional diquark
(quark-quark and antiquark-antiquark) bound-states appear. It can be shown
that these states have the same correlation functions as their mesonic
counterparts but have flipped parity [7, 8]. We will therefore restrict
ourselves to the study of quark-antiquark pseudoscalar and vector bound-
states. Note, that the corresponding quark-quark states are, however, scalars
and axialvectors.

Most interesting for dark matter physics are the properties of the lightest
states. These are the (pseudo)-Goldstone bosons as well as (likely) the
nonsinglet vector mesons. In the mass degenerate case these quantities have
already been studied in [9, 10].

Because the global symmetry is further broken by the non-degenerate quark
masses, these hadrons will no longer be degenerate: Similar to QCD, the
flavoured and unflavoured mesons will obtain different masses and decay
constants. In the present case the 5 Goldstones decompose into 2 fundamental
representations of each of the SU(2) groups and one singlet. Likwise, the
10-plet of flavoured vector mesons [11] decompose into a 6-plet of flavoured
particles and 4-plet of particles [6].

We will denote the unflavoured pseudoscalar by $\pi^{0}$ and the flavoured
ones as $\pi^{\pm}$ in analogy to QCD. The other pair is made up of quarks,
and remains degenerate to the $\pi$ [6]. Similarly we consider $\rho^{\pm}$
and the unflavoured ones $\rho^{0}$.

##  3 Lattice setup

We use a variant of the HiRep code [12] adapted for $Sp(N)$ gauge theories [9]
to simulate $Sp(4)$ gauge theory using the standard Wilson gauge action. We
include two fundamental Dirac-Wilson-fermions with non-degenerate masses using
the RHMC algorithm.

Figure 1: (left) Goldstone masses in units of the mass at degeneracy for all
three different values of $\beta$. For all three different values the almost
identical behaviour is observed. (right) Unflavoured Goldstone $\pi^{0}$
masses as a function of the spatial lattice extent $L$ for $\beta=7.2$ with
the lightest set of quark masses considered. An extrapolation to infinite
volume using the same fit function as in [10] is included. The value obtained
on the largest lattice agrees within errors with the finite volume
extrapolation.

We have chosen values for the inverse coupling $\beta>6.8$ so that we avoid
the bulk phase transition described in [10]. We have chosen three particular
values of $\beta=6.9,7.05,7.2$. For all these values of the inverse coupling
results for the mass degenerate case have already been obtained on larger
lattices in [10]. We fix the lighter of the two quark masses by requiring an
approximate ratio of $m_{\rho}/m_{\pi}\approx 1.25$ in the degenerate case.
This is suggested by dark matter phenomenology [13], which prefers stable
vector particles and Goldstones of order or above the confinement scale. The
other bare quark mass is then incrementally increased.

We considered lattices up to the size of $L^{3}\times T=14^{3}\times 24$ such
that the finite volume effects are under control. This is demonstrated in
figure 1. Additionally, we show in figure 1 a comparison of the three
different values of $\beta$. It can be seen that even in the case of the
smallest lattice spacing the infinite volume limit and the result obtained on
the largest lattice ($24\times 14^{3}$) agree within errors. We report that in
all ensembles studied the product of $m_{\pi}L$ was always $>5.5$. However, at
very large mass splittings lattice cutoff effects may become relevant, and
will be studied in the future. Details will be provided elsewhere.

Denoting the (dark) quark fields as $u$ and $d$ these hadrons are described by
operators $O_{\Gamma}^{f}$, specifically

$\displaystyle O_{\Gamma}^{\pm}(x,y)$ $\displaystyle=\bar{u}(x)\Gamma d(y)$
(1) $\displaystyle O_{\Gamma}^{0}(x,y)$ $\displaystyle=\bar{u}(x)\Gamma
u(y)-\bar{d}(x)\Gamma d(y).$ (2)

The decay constants are defined through the same parameterization of the
matrix elements as in [10],

$\displaystyle\langle 0|O_{\gamma_{5}\gamma_{\mu}}^{f}|PS\rangle$
$\displaystyle=f_{\pi}^{f}p_{\mu}$ $\displaystyle\langle
0|O_{\gamma_{\mu}}^{f}|V\rangle$
$\displaystyle=f_{\rho}^{f}m_{\rho}^{f}\epsilon_{\mu}$ (3)

where $|PS\rangle$ and $|V\rangle$ denote the ground state of the relevant
meson. Here $\epsilon_{\mu}$ is the polarisation vector which fulfills
$\epsilon_{\mu}\epsilon^{\mu}=1$ and $\epsilon_{\mu}p^{\mu}=0$. The masses are
extracted from the exponential fall-off of the zero-momentum correlator on a
lattice of temporal extent $T$ and spatial extent $L$ at large times $t$

$\displaystyle
C_{O}(\vec{p},t_{x}-t_{y})=\frac{1}{L^{3}}\sum_{\vec{x},\vec{y}}e^{-i\vec{p}(\vec{x}-\vec{y})}\langle
O(\vec{x},t_{x})O^{\dagger}(\vec{y},t_{y})\rangle$ (4) $\displaystyle
C_{O_{\Gamma}}(0,t)\xrightarrow[]{t\to\infty}\frac{|\langle
0|O_{\Gamma}|GS\rangle|^{2}}{2m_{\text{meson}}L^{3}}\left(e^{-m_{\text{meson}}t}+e^{-m_{\text{meson}}(T-t)}\right)$
(5)

where $|GS\rangle$ is the ground state of the meson sourced by the
interpolator $O_{\Gamma}$. In the rest frame the matrix element parameterizing
the pseudoscalar decay constant becomes

$\displaystyle\langle 0|O_{\gamma_{5}\gamma_{0}}|PS\rangle$
$\displaystyle=f_{\pi}m_{\pi}$ (6)

and we can extract all desired information by fitting the correlators at large
$t$ to

$\displaystyle C_{O_{\gamma_{5}\gamma_{0}}}(0,t)$
$\displaystyle\xrightarrow[]{t\to\infty}\frac{(f_{\pi})^{2}m_{\pi}}{2L^{3}}\left(e^{-m_{\pi}t}+e^{-m_{\pi}(T-t)}\right)$
(7) $\displaystyle C_{O_{\gamma_{\mu}}}(0,t)$
$\displaystyle\xrightarrow[]{t\to\infty}\frac{(f_{\rho})^{2}m_{\rho}}{2L^{3}}\left(e^{-m_{\rho}t}+e^{-m_{\rho}(T-t)}\right).$
(8)

For the unflavoured mesons disconnected contributions to the correlator appear
that vanish exactly in the mass-degenerate limit. These are, however,
suppressed by the relatively large quark masses and the small mass difference
$(m_{u}\\!-\\!m_{d})$ and are therefore for the moment neglected.

##  4 Results

Figure 2: Masses of the flavoured and unflavoured Goldstone and iso-nonsinglet
vector mesons masses in lattice units. At sufficiently large fermionic mass
difference the unflavoured vectors get even lighter than the flavoured
Goldstones. Figure 3: Decay constants of the flavoured and unflavoured
Goldstone and iso-nonsinglet vector mesons in lattice units. For large mass
differences of the fermions the decay constants show hints of non-
monotonicity.

The results for the flavoured and unflavoured pseudo-Goldstone and vector
meson masses and decay constants are depicted in figures 2 and 3. The results
are given in lattice units. If desired, a scale can be introduced, e. g. by
fixing the lattice spacing $a$.

We see a clear separation of the flavoured and unflavoured meson states and
can conclude that they are no longer mass-degenerate. For both the pseudo-
Goldstones and the vectors the unflavoured ones are lighter than the others.
Furthermore, at sufficiently large fermion mass difference the unflavoured
vectors get even lighter than the flavoured pseudo-Goldstones. Since one
fermion mass is kept fixed while the other is incrementally increased all
mesons get heavier the further we move away from degeneracy.

In the case of the form factors we again observe that the flavoured and
unflavoured ones are no longer degenerate and the unflavoured form factors are
smaller. In addition we see hints of non-monotonicity in both form factors.

An important observation is that we see clearly qualitatively different
behaviors. Especially, the unflavoured ones do not show a (dominant) linear
increase with the mass of one of the flavours. For the decay constants the
effect is even different between flavoured pseudoscalars and vectors. This
strongly suggests that the composition of the dark hadrons changes
substantially. This can have implications for decay patterns, which will be
studied elsewhere [6]. Note, however, that in our case even at the largest
mass spitting the vector mesons remain stable.

##  5 Summary

We have studied $Sp(4)$ gauge theory with two non-degenerate massive
fundamental Dirac fermions. Four of the five pseudo-Goldstone bosons remain
mass-degenerate for $m_{u}\neq m_{d}$ while the unflavoured one aquires a
lighter mass and is therefore the lightest state in this theory.

We see strong signs for non-linear behavior, which suggests that the flavour
composition of the dark hadrons changes. This can have interesting
consequences for phenomenology. Also, our results show that it is possible to
extract the splitting effects sufficiently well that interpolation even to
tiny mass splittings, as preferred in some dark matter scenarios [14], would
be accessible. A more detailed analysis of the phenomenological implications
as well of lattice details is forthcoming [6].

## Acknowledgments

FZ is supported by the the Austrian Science Fund research teams grant STRONG-
DM (FG1). The computations have been performed on the HPC cluster of the
University of Graz and on the Vienna Scientific Cluster (VSC4). We are
grateful to the authors of [9, 10] for the support with the HiRep code.

## References

  * [1] Y. Hochberg, E. Kuflik, T. Volansky and J.G. Wacker, _Mechanism for Thermal Relic Dark Matter of Strongly Interacting Massive Particles_ , _Phys. Rev. Lett._ 113 (2014) 171301 [1402.5143].
  * [2] Y. Hochberg, E. Kuflik, H. Murayama, T. Volansky and J.G. Wacker, _Model for Thermal Relic Dark Matter of Strongly Interacting Massive Particles_ , _Phys. Rev. Lett._ 115 (2015) 021301 [1411.3727].
  * [3] Y. Hochberg, E. Kuflik and H. Murayama, _SIMP Spectroscopy_ , _JHEP_ 05 (2016) 090 [1512.07917].
  * [4] N. Bernal, X. Chu, S. Kulkarni and J. Pradler, _Self-interacting dark matter without prejudice_ , _Phys. Rev. D_ 101 (2020) 055044 [1912.06681].
  * [5] D.A. Kosower, _SYMMETRY BREAKING PATTERNS IN PSEUDOREAL AND REAL GAUGE THEORIES_ , _Phys. Lett. B_ 144 (1984) 215.
  * [6] S. Kulkarni, A. Maas, S. Mee, M. Nikolic, J. Pradler and F. Zierler in preparation.
  * [7] S. Hands, P. Sitch and J.-I. Skullerud, _Hadron Spectrum in a Two-Colour Baryon-Rich Medium_ , _Phys. Lett. B_ 662 (2008) 405 [0710.1966].
  * [8] R. Lewis, C. Pica and F. Sannino, _Light Asymmetric Dark Matter on the Lattice: SU(2) Technicolor with Two Fundamental Flavors_ , _Phys. Rev. D_ 85 (2012) 014504 [1109.3513].
  * [9] E. Bennett, D.K. Hong, J.-W. Lee, C.J.D. Lin, B. Lucini, M. Piai et al., _Sp(4) gauge theory on the lattice: towards SU(4)/Sp(4) composite Higgs (and beyond)_ , _JHEP_ 03 (2018) 185 [1712.04220].
  * [10] E. Bennett, D.K. Hong, J.-W. Lee, C.J.D. Lin, B. Lucini, M. Piai et al., _Sp(4) gauge theories on the lattice: $N_{f}=2$ dynamical fundamental fermions_, _JHEP_ 12 (2019) 053 [1909.12662].
  * [11] E. Bennett, D.K. Hong, J.-W. Lee, C.-J.D. Lin, B. Lucini, M. Mesiti et al., _$Sp(4)$ gauge theories on the lattice: quenched fundamental and antisymmetric fermions_, _Phys. Rev. D_ 101 (2020) 074516 [1912.06505].
  * [12] L. Del Debbio, A. Patella and C. Pica, _Higher representations on the lattice: Numerical simulations. SU(2) with adjoint fermions_ , _Phys. Rev. D_ 81 (2010) 094503 [0805.2058].
  * [13] A. Berlin, N. Blinov, S. Gori, P. Schuster and N. Toro, _Cosmology and Accelerator Tests of Strongly Interacting Dark Matter_ , _Phys. Rev. D_ 97 (2018) 055033 [1801.05805].
  * [14] N. Bernal, X. Chu and J. Pradler, _Simply split strongly interacting massive particles_ , _Phys. Rev. D_ 95 (2017) 115023 [1702.04906].