# Reducing Uncertainty in Sea-level Rise Prediction: A Spatial-variability-
aware Approach

Subhankar Ghosh  Dept. of Computer Science & Eng.
University of Minnesota, Twin Cities
Minneapolis, USA
<EMAIL_ADDRESS>Shuai An  Dept. of Computer Science & Eng.
University of Minnesota, Twin Cities
Minneapolis, USA
<EMAIL_ADDRESS>Arun Sharma  Dept. of Computer Science & Eng.
University of Minnesota, Twin Cities
Minneapolis, USA
<EMAIL_ADDRESS>Jayant Gupta  Dept. of Computer Science & Eng.
University of Minnesota, Twin Cities
Minneapolis, USA
<EMAIL_ADDRESS>Shashi Shekhar  Dept. of Computer Science & Eng.
University of Minnesota, Twin Cities
Minneapolis, USA
<EMAIL_ADDRESS>Aneesh Subramanian  Dept. of Atmospheric & Oceanic Sciences
University of Colorado, Boulder
Boulder, USA
<EMAIL_ADDRESS>
###### Abstract

Given multi-model ensemble climate projections, the goal is to accurately and
reliably predict future sea-level rise while lowering the uncertainty. This
problem is important because sea-level rise affects millions of people in
coastal communities and beyond due to climate change’s impacts on polar ice
sheets and the ocean. This problem is challenging due to spatial variability
and unknowns such as possible tipping points (e.g., collapse of Greenland or
West Antarctic ice-shelf), climate feedback loops (e.g., clouds, permafrost
thawing), future policy decisions, and human actions. Most existing climate
modeling approaches use the same set of weights globally, during either
regression or deep learning to combine different climate projections. Such
approaches are inadequate when different regions require different weighting
schemes for accurate and reliable sea-level rise predictions. This paper
proposes a zonal regression model which addresses spatial variability and
model inter-dependency. Experimental results show more reliable predictions
using the weights learned via this approach on a regional scale.

###### Index Terms:

climate change, sea-level rise, spatial variability, forecasting, machine
learning, regression

##  I Introduction

Sea-level rise is one of the most pressing environmental challenges facing the
world today. In 2021 Intergovernmental Panel on Climate Change (IPCC) [26]
projected that the global sea level is likely to rise by 10 to 32 inches (26 -
82 cm) by 2100, relative to 1986-2005. However, these global projections mask
significant regional variation: sea level is projected to rise more in the
Arctic than in the tropics and more in the mid-latitudes than in the
subtropics. For example, in Chesapeake Bay (near Maryland and Washington DC),
projections based on tide gauges, satellite observations and geo-physical
models indicate that sea level will probably rise between 12-20 inches (30 -
50 cm) by 2050 (starting from 2005) [25]. The need for regional sea-level rise
projections is becoming increasingly urgent. Coastal communities worldwide are
already experiencing the impacts of sea-level rise, such as flooding, erosion,
and saltwater intrusion. These impacts are expected to worsen as the sea level
continues to rise. Regional sea-level rise projections are essential for
coastal planning and adaptation. They can help communities identify areas that
are most at risk and develop strategies to mitigate and adapt to these
impacts.

Given a multi-model ensemble of climate projections from state-of-the-art
climate models of the present, the goal is to accurately and reliably predict
future sea-level rise while lowering the prediction uncertainty. A spatially-
variable weighting scheme assigns different sets of weights to the candidate
models in different regions unlike global climate models (GCMs), which assign
the same set of weights across different regions. Using a spatially variable
weighting scheme ensures that higher priority is given to more reliable models
in a specific region when performing predictions and reducing uncertainty.

Figure 1: Interaction between ice sheets, ocean, and atmosphere on the
Antarctic and Greenland ice sheets. [1]

Predictions from climate models are improved when obtained via multi-model
ensembles. The result is more accurate and produces more reliable predictions
with lower uncertainty. The earliest works generally used an average of the
simulated predictions to obtain the final result. However, this method would
not entirely capture the prediction skill of the models in different regions.
More recent approaches [7, 8] use a global weighting scheme based on
regression or deep learning methods to reduce the discrepancy between
observations and the weighted prediction obtained from the simulations. These
models assign greater weight to more reliable predictions on a global scale
but do not adequately handle the spatial variability arising from regional
factors which affect changes in sea level. They can approximate the thermal
expansion of the oceans but cannot properly account for regional factors like
tsunamis, tides, boundary waves, river runoff, vertical land movement, etc.
The spatial variability arising from such regional factors is a significant
challenge that can lead to different weighting schemes based on the model’s
reliability in that specific region.

This paper explores a zonal regression model to address spatial variability in
sea-level rise prediction. This approach assigns weights to models based on
their prediction skills in a specific region, looking for improved
predictions. Traditional weighting strategies use an inverse of the root mean
square between the observations and the simulations to assign the weights on a
global scale. This paper proposes a spatial-variability-aware regression model
to learn the weights from historical predictions in a region. Empirical
results in Figure 5 using root mean squared errors from the different model
predictions, the baseline method, and the proposed approach display the
effectiveness of the approach in reducing the prediction error.

##  II Background

Global sea levels can be affected by factors like warming oceans and the
addition of freshwater from continental ice (Figure 1). Ocean dynamics can
bring warm waters in contact with glaciers, leading to the decay of ice
sheets. Changes in atmospheric temperature can increase surface temperatures
over continents leading to further loss of ice sheets. This reduction in the
ice mass balance can impact the sea level in different regions of the earth.
Other factors such as tsunamis, tides, boundary waves, river runoff, and
vertical land movement can also impact regional sea levels.

Figure 2 provides a comparison of the significant factors contributing to
changes in the global mean sea level (GMSL). Different satellite systems
measure different GMSL properties. For example, mass-driven changes such as
ice melt and land water storage are measured by satellite altimeters, like
GRACE/GRACE-FO [15]. Argos [16] measures ocean heat content by monitoring
factors such as temperature and salinity changes in the ocean.

Figure 2: Contribution to sea-level rise from thermal expansion and ocean
mass. [6]

##  III Proposed Approach

This paper explores a zonal regression method to represent spatial variability
in weighting multi-model ensemble climate projections. CMIP-6 projections
obtained from the National Center for Atmospheric Research (NCAR) are used as
the simulated data, while Copernicus sea-level anomaly gridded data are used
for the ground truth observation.

Given a set of spatio-temporal tensors representing climate model projections
($P_{r}$) over a region $r$, the proposed weighting strategy aims to find a
set of coefficients ($w_{r}$) for the multi-model ensemble such that the
combination of projections with weights ($\sum_{i}P_{r,i}w_{r,i}$) improves
the accuracy of the sea-level rise prediction while simultaneously reducing
the prediction uncertainty. Here, $i$ represents the model under
consideration.

We formulate the problem as a non-negative least squares problem where the
goal is to find the set of weights $w_{r}$ that minimizes the prediction
error:

$\displaystyle\operatorname*{argmin}_{w_{r,i}}||O_{r}-\sum_{i}P_{r,i}w_{r,i}||,$
(1)

where $w_{r,i}>=0$ and $\sum_{i}w_{r,i}=1$. $O_{r}$ represent the ground truth
observation data over region $r$.

Simultaneously, we also need to reduce the prediction uncertainty:

$\displaystyle minimize\sqrt{\frac{1}{T-1}\sum_{t=1}^{T}(E_{t}-\bar{E})^{2}},$
(2)

where $E_{t}=||O_{r,t}-\sum_{i}P_{r,t,i}w_{r,i}||$. This can be considered as
a multi-objective optimization problem in which we want to minimize both the
prediction error and the prediction uncertainty. One technique for solving
this problem is to use Pareto optimality. The main goal in this technique is
to find solutions such that each quantity (prediction error and prediction
uncertainty) is minimized without affecting the other. Thus we want to find
the Pareto frontier (i.e. set of Pareto efficient solutions) as in Figure 3. A
point $y^{\prime\prime}\in\mathbb{R}^{m}$ is preferred to (dominates) another
point $y^{\prime}\in\mathbb{R}^{m}$ is referred to as $y^{\prime\prime}\succ
y^{\prime}$. In Figure 3, both $A$ and $B$ dominate $C$ in reducing the
prediction error and the uncertainty.

Model similarity can lead to projections biased towards the largest set of
similar models and the underestimation of uncertainties [23]. To address this,
the proposed method down-weights models with a higher covariance in their
predictions. The final weights obtained are a product of the weights obtained
from the regression coefficients and the coefficients if the model is down-
weighted. This zonal regression formulation overcomes the drawbacks of one-
size-fits-all (OSFA) models such as GCMs which do not account for spatial
variability.

##  IV Experimental framework

###  IV-A Experimental Evaluation

A zonal regression model is employed to learn a weighting scheme for the
simulated projections such that the output of the weighted predictions reduces
the discrepancy with the observations. The main experimental goal was to
compare the sea-level rise predictions obtained with the proposed spatial
variability-aware zonal regression approach against the sea-level rise
predictions from a one-size-fits-all method that doesn’t adapt to zones.
Figure 4 shows the overall validation framework. The metric for comparison
with the baseline [7, 8] was solution quality, specifically the root mean
squared error between the combined prediction and the ground truth
observations. The experiments were performed with the CMIP-6 projection models
of sea-level height and the Copernicus satellite data for ground truth.

Figure 3: Pareto frontier for simultaneous reduction of prediction error and
uncertainty.
Figure 4: Experiment Design

###  IV-B Data

This work primarily uses Copernicus [9] sea-level gridded satellite data for
ground truth observations and zos (Sea Surface Height Above Geoid) data from
CMIP-6 (Coupled Model Intercomparison Project 6) for the simulated projections
available from NCAR (National Center for Atmospheric Research).

The sea level satellite observation data [9] from Copernicus provides daily
and monthly global estimates of mean sea level anomaly based on satellite
altimetry measurements from 1993 to 2012. Here sea level anomaly refers to the
height of the sea surface in comparison to the mean sea level at a given time
and region. The dataset is generated using a two-satellite merged
constellation where one satellite is used as a reference while the other is
used to improve accuracy.

The Coupled Model Intercomparison Project (CMIP) [10] is a collaborative
framework for collecting, organizing, and distributing output from multiple
coupled climate models performing common sets of experiments. It is focused on
global climate model (GCM) simulations of the past, current and future
climate. The major components include (1) DECK (Diagnostic, Evaluation, and
Characterization of Klima), which are some common experiments, and CMIP
historical simulations (1850–present); (2) common standards, coordination,
infrastructure, and documentation; (3) ensembles of CMIP-Endorsed Model
Intercomparison Projects (MIPs) specific to a phase of CMIP (CMIP6) that build
upon the DECK and CMIP historical simulations. Some data variables of interest
available from CMIP-6 include air pressure at sea level, air temperature,
surface temperature, sea-level anomaly, surface runoff flux, etc. The
experiments mentioned in this paper focus on sea-level anomaly (zos) data.

Experiments mentioned in this paper were focused on the Eastern North America
(ENA) region due to the availability of more reliable ground truth data in the
ENA region. The experiments were performed on four CMIP-6 models obtained from
NCAR, namely CESM2 [17], CESM2-FV2 [18], CESM2-WACCM [19] and CESM2-WACCM-FV2
[20].

##  V Preliminary results

Figure 5 provides a comparison of the root mean squared errors (RMSE) in sea-
level rise predictions from four simulation models, the baseline method, and
the proposed approach. The RMSE obtained from each model prediction was higher
than the weighted combination of the predictions. In the ENA region, the
proposed method had a median RMSE of $0.35$ while the baseline method [7, 8]
had a median RMSE of $0.36$ across the different time scales. The standard
deviation of the prediction was $0.0177$ for the proposed method while the
standard deviation of the baseline method was $0.0158$. The outliers in a
boxplot represent datapoints that are located outside the whiskers of the box
plot. For example, outliers in Figure 5 are data points $1.5$ times outside
the interquartile range ($IQR$) above the upper quartile ($Q3$) and below the
lower quartile ($Q1$), i.e. ($Q1-1.5*IQR$ or $Q3+1.5*IQR$). This shows that
accounting for spatial variability leads to more reliable predictions.

Figure 5: Root mean squared errors of each model.

##  VI Related work and Discussion

The Fourth National Climate Assessment (2018) mentions a weighting strategy
used by the Intergovernmental Panel on Climate Change for differential
weighting of climate projection models. The approach, based on works by
Sanderson et al. [7, 8], incorporates both skills in climatological
performance and the inter-dependency of models due to factors such as common
parameterizations. A set of weights are obtained via a ranking scheme of the
climate models based on the inverse of the root mean square error (RMSE)
between the observations and the predictions. The final output is a single set
of weights that can be used to find a weighted average of the climate
projections. Inter-dependency between models is determined based on inter-
model RMSE between each pair of models. A model is down-weighted if its inter-
model RMSE with another model is significantly lower than a predetermined
threshold. The final weight for each model is a product of its skill and its
inter-dependency weights. This final set of weights obtained is used for
combining projections irrespective of region and thus the strategy ignores the
spatial variability from various regional physical forcings.

Spatial variability is receiving increasing attention in pattern detection
[12, 13, 27, 28] and prediction studies [11, 29]. Recent approaches [14]
towards predicting regional sea-level rise with machine learning require
explicit input values for the regional physical factors, which might be
difficult to obtain in certain cases. Such bottlenecks can be avoided by
employing the proposed approach.

##  VII Conclusion and Future Work

This paper discussed a new spatial variability-aware zonal regression approach
towards weighting sea-level rise models. Rather than relying on a model’s
global projection skill, this method learns a weighting scheme that accounts
for regional factors in sea-level rise predictions.

Future Work: We plan to extend this work into more sophisticated approaches
such as physics informed models using GWR (geographically weighted regression)
or spatial-variability-aware deep learning methods [11] and work on novel ways
to handle the inter-dependency between the models. A potential future
direction of this work would be to make sea-level rise predictions at a higher
resolution.

Sea-level rise is directly tied to the rise in greenhouse gas emissions.
Therefore, we can improve the quantification of uncertainty in sea-level rise
projections by considering different scenarios that are likely to affect
levels of greenhouse gas emissions in the future. For example, with no human
intervention, scientists project emissions will double by 2100. If all
countries comply with the Paris Agreement, emissions should reach net zero
after 2050. A third scenario assumes countries only meet their current
national commitments for reducing emissions. In our future work, we plan to
explore how such emissions scenarios can be used to reduce uncertainty in sea-
level rise projections.

## Acknowledgment

This material is based upon work supported by the National Science Foundation
under Grants No. 2118285, 2040459, 1901099, and 1916518. We also thank Kim
Koffolt, the iHarp community, and the Spatial Computing Research Group for
valuable comments and refinements.

## References

  * [1] J. Willis and J. Church, “Regional sea-level projection,” Science, vol. 336, pp. 550–551, 2012.
  * [2] Intergovernmental Panel On Climate Change, “Climate change 2007: The physical science basis,” Agenda, vol. 6., no. 07, 2007, pp.333.
  * [3] M. Vermeer and S. Rahmstorf, “Global sea level linked to global temperature,” Proceedings of the national academy of sciences, vol. 106, no. 51, National Acad Sciences, 2009, pp. 21527–21532.
  * [4] W. Pfeffer, J. Harper and S. O’Neel, “Kinematic constraints on glacier contributions to 21st-century sea-level rise,” Science, vol. 321, no. 5894, American Association for the Advancement of Science, pp. 1340–1343, 2008.
  * [5] R. Nicholls et al., “Sea-level rise and its possible impacts given a ‘beyond 4 C world’in the twenty-first century,” Philosophical transactions of the Royal Society A: mathematical, physical and engineering sciences, vol. 369, no. 1934, The Royal Society Publishing, pp. 161–181, 2011.
  * [6] NASA, “Understanding Sea Level,” https://sealevel.nasa.gov/understanding-sea-level/global-sea-level/overview.
  * [7] B. Sanderson, R. Knutti and P. Caldwell, “A representative democracy to reduce interdependency in a multimodel ensemble,” Journal of Climate, vol. 28, no. 13, American Meteorological Society, pp. 5171–5194, 2015.
  * [8] B. Sanderson, M. Wehner and R. Knutti, “Skill and independence weighting for multi-model assessments,” Geoscientific Model Development, vol. 10, no. 6, Copernicus GmbH, pp. 2379–2395, 2017.
  * [9] Copernicus, “Sea level gridded data from satellite observations for the global ocean from 1993 to present,” https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-sea-level-global.
  * [10] V. Eyring et al., “Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization,” Geoscientific Model Development, vol. 9, no. 5, Copernicus GmbH, pp. 1937–1958, 2016.
  * [11] J. Gupta et al., “Spatial variability aware deep neural networks (SVANN): a general approach,” ACM Transactions on Intelligent Systems and Technology (TIST), vol. 12, no. 6, ACM New York, NY, pp. 1–21, 2021.
  * [12] S. Ghosh, J. Gupta, A. Sharma, S. An and S. Shekhar, “Towards geographically robust statistically significant regional colocation pattern detection,” Proceedings of the 5th ACM SIGSPATIAL International Workshop on GeoSpatial Simulation, pp. 11–20, 2022.
  * [13] S. Ghosh, J. Gupta, A. Sharma, S. An and S. Shekhar, “Reducing False Discoveries in Statistically-Significant Regional-Colocation Mining: A Summary of Results,” In 12th International Conference on Geographic Information Science (GIScience 2023)., vol. 277, Schloss Dagstuhl – Leibniz-Zentrum für Informatik, pp. 3:1–3:18, 2023.
  * [14] V. Nieves, C. Radin, G. Camps-Valls, “Predicting regional coastal sea level changes with machine learning,” Scientific Reports, vol. 11, no. 1, Nature Publishing Group UK London, pp. 7650, 2021.
  * [15] NASA Jet Propulsion Lab, “Measuring Earth’s Surface Mass and Water Changes,” https://gracefo.jpl.nasa.gov.
  * [16] Argos (satellite system), “https://www.argos-system.org”.
  * [17] Community Earth System Model 2 (CESM2), “https://www.cesm.ucar.edu/models/cesm2”.
  * [18] Community Earth System Model version 2 - FV2, “https://www.cen.uni-hamburg.de/en/icdc/data/cryosphere/cmip6-sea-ice-area/cesm2-fv2.html”.
  * [19] Whole Atmosphere Community Climate Model, “https://catalogue.ceda.ac.uk/uuid/108b11504ef44f4fb0bb66273e261709”.
  * [20] CESM2-WACCM-FV2, “https://www.wdc-climate.de/ui/cmip6?input=CMIP6.CMIP.NCAR.CESM2-WACCM-FV2.historical”.
  * [21] NASA-GISS-E2, “https://data.giss.nasa.gov/modelE/cmip6/”.
  * [22] AWI Climate Model (AWI-CM), “https://www.awi.de/en/science/climate-sciences/atmospheric-physics/climate-modelling.html”.
  * [23] R. Pathak, H. Dasari, K. Ashok, I Hoteit, “Effects of Observational Uncertainty and Models Similarity on Climate Change Projections,” Research Square Platform LLC, 2023.
  * [24] M. Iturbide, “An update of IPCC climate reference regions for subcontinental analysis of climate model data: definition and aggregated datasets,” Earth System Science Data, vol. 12, no. 4, Copernicus GmbH, pp. 2959–2970, 2020.
  * [25] D. Boesch et al., “Sea-level Rise Projections for Maryland 2023,” UMCES reports on Sea-Level Rise Projections for Maryland, 2023.
  * [26] J. Lee et al., “Future global climate: scenario-based projections and near-term information,” Climate change 2021: The physical science basis. Contribution of working group I to the sixth assessment report of the intergovernmental panel on climate change, Cambridge University Press, pp. 553–672, 2021.
  * [27] A. Sharma, J. Gupta and S. Ghosh, “Towards a tighter bound on possible-rendezvous areas: preliminary results,” In Proceedings of the 30th International Conference on Advances in Geographic Information Systems, pp. 1-11, 2022.
  * [28] M. Farhadloo et al., “SAMCNet: towards a spatially explainable AI approach for classifying MxIF oncology data,” Proceedings of the 28th ACM SIGKDD conference on knowledge discovery and data mining, pp. 2860–2870, 2022.
  * [29] Y. Li et al., “Eco-PiNN: A Physics-informed Neural Network for Eco-toll Estimation,” Proceedings of the 2023 SIAM International Conference on Data Mining (SDM), pp. 838-846, 202.


# Energy-Dependent Analyses of the Gamma-Ray Emission from HESS J1857+026 with
Fermi-LAT

Xiaolei Guo Xi Liu  Yuliang Xin School of Physical Science and Technology,
Southwest Jiaotong University, Chengdu 610031, People’s Republic of China;
<EMAIL_ADDRESS>
###### Abstract

We report the discovery of energy-dependent morphology for the GeV gamma-ray
emission from HESS J1857+026 with more than 13 years of Fermi Large Area
Telescope (LAT) data. The GeV gamma-ray emission from this region is composed
of two extended components. The hard component with an index of $1.74\pm 0.07$
in the energy range of 0.5-500 GeV is spatially coincident with HESS
J1857+026, and its 68% containment radius varies from $\sim 0.44^{\circ}$
below 40 GeV to $\sim 0.30^{\circ}$ above 140 GeV. The hard GeV gamma-ray
spectrum and the energy-dependent morphology of HESS J1857+026 make it favor a
PWN origin, which is associated with the energetic pulsar, PSR J1856+0245. The
soft component with an index of $2.70\pm 0.16$ and another extended gamma-ray
source detected in this region, 4FGL J1857.9+0313e with an index of $2.55\pm
0.07$, are spatially coincidence with two molecular clumps in the northeast
and southwest of HESS J1857+026, which favors the hadronic process, and the
protons could be accelerated by the hypothetical SNR associated with PSR
J1856+0245.

gamma rays: general - gamma rays: ISM - ISM: individual objects (HESS
J1857+026) - radiation mechanisms: non-thermal

##  1 Introduction

More than two hundred of very-high-energy (VHE; $>$ 100 GeV) gamma-ray sources
have been detected with the operations of ground-based Cherenkov telescopes,
such as H.E.S.S. (Aharonian et al., 2004), MAGIC (Aleksić et al., 2016), HAWC
(Abeysekara et al., 2013) and LHAASO (Cao et al., 2023). The gamma-ray
emission could be produced by the hadronic interactions in which $\pi^{0}$
decays into two gamma-ray photons subsequently, via the inverse Compton (IC)
scattering process or through non-thermal bremsstrahlung radiation from high
energy electrons. Most of the VHE gamma-ray sources have been identified to be
pulsar wind nebulae (PWNe), supernova remnants (SNRs), X-ray binaries, TeV
halos and so on. However, some TeV gamma-ray sources are still unidentified,
and the multi-wavelength studies of these sources are crucial for revealing
their nature and probing the origin of cosmic rays.

HESS J1857+026 was first discovered as a VHE gamma-ray source with an
extension of $(0.11\pm 0.08)^{\circ}\times(0.08\pm 0.03)^{\circ}$ (Aharonian
et al., 2008). In H. E. S. S. Collaboration et al. (2018), the morphology of
HESS J1857+026 was described by a two dimensional (2D) Gaussian component with
an approximate size of $0.26^{\circ}\pm 0.06^{\circ}$. Hessels et al. (2008)
discovered an energetic pulsar, PSR J1856+0245, in the direction of HESS
J1857+026, which makes HESS J1857+026 to be a potential PWN candidate. The
period and spin-down luminosity of PSR J1856+0245 are P = 81 ms and $\dot{E}$
= $4.6\times 10^{36}$ erg s-1, respectively, with the characteristic age of
$\tau_{c}$ = 21 kyr. The distance of PSR J1856+0245 was first estimated to be
$\sim$ 9 kpc, derived by the dispersion measure (DM) with electron density
model of Cordes & Lazio (2002). And, an updated distance of $\sim 6.3$ kpc is
given by ATNF Pulsar Catalog
111https://www.atnf.csiro.au/research/pulsar/psrcat/ according to the electron
density model of Yao et al. (2017).

MAGIC carried out follow-up observations and presented energy-dependent
morphology of this region (MAGIC Collaboration et al., 2014). In the energy
range of 300 GeV-1 TeV, the morphology observed by MAGIC is compatible with
H.E.S.S. observation. However, two separate gamma-ray sources, named MAGIC
J1857.2+0263 (hereafter MAG1) and MAGIC J1857.6+0297 (hereafter MAG2), were
detected with the data above 1 TeV. MAG1 is an extended source with the
intrinsic extension of $(0.17\pm 0.03_{\rm stat}\pm 0.02_{\rm
sys})^{\circ}\times(0.06\pm 0.03_{\rm stat}\pm 0.02_{\rm sys})^{\circ}$, while
MAG2 is compatible with a point source. MAGIC Collaboration et al. (2014)
interpreted MAG1 as a PWN powered by PSR J1856+0245, and MAG2 may be
associated with a molecular cloud complex containing an Hii region located at
$\sim$ 3.7 kpc and a possible gas cavity. Another pulsar, PSR J1857+0300, was
discovered in the direction of MAG2 (Lyne et al., 2017). The characteristic
age and spin-down luminosity of PSR J1857+0300 are $\tau_{\rm c}$ $\sim
4.6\times 10^{6}$ yr and $\dot{E}$ $\sim 2.3\times 10^{32}\ {\rm erg}\ {\rm
s^{-1}}$, with the distance of $\sim 6.7$ kpc. Meanwhile, an elliptical
superbubble, was detected with neutral gas observation (Petriella et al.,
2021), which is also spatially coincident with HESS J1857+026. The kinematical
distance of the superbubble is about 5.5 kpc and is close to the DM distance
of PSR J1856+0245 ($\sim 6.3$ kpc). Petriella et al. (2021) concluded that the
TeV emission of HESS J1857+026 originates from the superbubble, and PSR
J1856+0245 is located inside the superbubble. In addition, they found five
molecular components in the velocity interval between 78 and 90 $\rm{km}\
\rm{s}^{-1}$ with ${}^{13}\rm{CO}$(J=1-0) observations, which are probably
associated with the superbubble. And they favor a single gamma-ray source
scenario instead of the superposition of two gamma-ray sources. Petriella et
al. (2021) also performed radio observations at 1.5 GHz and 6.0 GHz with VLA.
Nevertheless, no significant radio emission was detected in this region.

In the energy range of 1-25 TeV, the Water Cherenkov Detector Array (WCDA) of
Large High Altitude Air Shower Observatory (LHAASO) detected an extended
gamma-ray source, 1LHAASO J1857+0245, which is spatially consistent with HESS
J1857+026 (Cao et al., 2023). The spatial template of 1LHAASO J1857+0245 is
described by a 2D Gaussian with $\sigma=0.24^{\circ}\pm 0.04^{\circ}$, and its
gamma-ray spectrum in 1-25 TeV is modeled by a power law with an index of 2.93
$\pm$ 0.07.

The GeV gamma-ray emission from HESS J1857+026 was first detected using a
point source hypothesis by Fermi-LAT (Atwood et al., 2009), while no gamma-ray
pulsation from PSR J1856+0245 was observed (Rousseau et al., 2012). Since no
obvious X-ray emission was detected, only an upper limit of $F_{1-10\ {\rm
keV}}=2\times 10^{-12}$ erg cm-2 s-1 was obtained (Rousseau et al., 2012).
Considering the uncertain origin and the complexity of this region, the
detailed analysis with more GeV observational data will be helpful to
investigate the origin of the gamma-ray emission.

Taking advantage of more than 13 years of Fermi-LAT data, we performed an
energy-dependent analyses of the region around HESS J1857+026 and discussed
the nature of this source. The paper is organized as follows. In Section 2, we
describe the data analysis routines and present our results. In Section 3, we
discuss the radiation mechanisms and the nature of HESS J1857+026. And the
conclusion for this work is presented in Section 4.

##  2 Fermi-LAT Data and Results

The Fermi-LAT Pass 8 data we analyzed are collected from August 4, 2008 to
March 4, 2022 with energies from 500 MeV to 500 GeV. The region of interest
(ROI) is a $10^{\circ}\times 10^{\circ}$ square centered at the TeV gamma-ray
centroid of HESS J1857+026 (R.A. = 284.296∘, Dec. = 2.667∘; H. E. S. S.
Collaboration et al., 2018). To reduce the contamination from the Earth Limb,
events with zenith angle greater than $90^{\circ}$ are eliminated. The
standard analysis software of Fermitools
222http://fermi.gsfc.nasa.gov/ssc/data/analysis/software/ is used with the
instrumental response function (IRF) of P8R3_SOURCE_V3. The models we used to
describe the Galactic and isotropic diffuse emissions include
gll_iem_v07.fits, iso_P8R3_SOURCE_V3_v1.txt and iso_P8R3_SOURCE_V3_PSF3_v1
333http://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html. In
addition, all sources in the incremental version of the fourth Fermi-LAT
source catalog (4FGL-DR3; Abdollahi et al., 2020, 2022) within a radius of 15∘
centered at HESS J1857+026 are included in the model. And the binned maximum
likelihood analysis method with gtlike is applied. During the fitting
procedure, the spectral parameters of all sources located in the ROI are left
free, together with the normalizations of the two diffuse backgrounds.

###  2.1 Spatial Analysis

####  2.1.1 Average Spatial Extension

In the 4FGL-DR3 catalog, the gamma-ray emission of HESS J1857+026 is described
by an uniform disk (named as 4FGL J1857.7+0246e) centered at R.A. =
$284.449^{\circ}$, Dec. = $2.774^{\circ}$ with the 68% containment radius of
$r_{\rm 68}=0.50^{\circ}$, which is given by the analysis of Fermi-LAT
extended Galactic sources (FGES; Ackermann et al., 2017). Meanwhile, there is
a point source named as 4FGL J1857.9+0313c located in the north of the disk,
which has no identified counterpart (Abdollahi et al., 2020). And in the south
of the disk, two point sources named as 4FGL J1857.6+0212 and 4FGL
J1858.3+0209 are identified to be associated with SNR G35.6-0.4 and HESS
J1858+020 (Cui et al., 2021; Zhang et al., 2022). To obtain the spatial
template of HESS J1857+026, we performed the spatial extension analysis in the
energy ranges of 1-3 GeV (low energy) and 10-500 GeV(high energy),
respectively. And in the low energy band, only “PSF3” type (evclass=128 &
evtype=32) data with better angular resolution are selected to reduce the
contamination from nearby sources, while the data with “SOURCE” type
(evclass=128 & evtype=3) are used in the high energy range. After subtracting
the background sources included in the model (except 4FGL J1857.9+0313c in
left panel and 4FGL J1857.7+0246e in the middle and right panel), we created
three $2.5^{\circ}\times 2.5^{\circ}$ TS maps centered at HESS J1857+026 with
the different energy ranges, as shown in Figure 1.

As shown in the left panel of Figure 1, the position of 4FGL J1857.9+0313c
given by 4FGL-DR3 is not coincident with the the gamma-ray peak of this
region. In addition, there are discrepancies between the morphologies of the
GeV emission around HESS J1857+026 and the spatial templates given by FGES
both in the low and high energy bands, as shown in the middle and right panels
of Figure 1. Therefore, we refined the spatial templates of HESS J1857+026 of
the low (hereafter referred to be “Src A”) and high (“Src T”) energies, and
the morphology of 4FGL J1857.9+0313c was also re-analyzed. The uniform disk,
2D Gaussian, as well as the H.E.S.S. image model are tested. Meanwhile, a
point source model with the best-fit coordinate calculated by gtfindsrc is
also applied for Src A and 4FGL J1857.9+0313c to test the spatial extension of
it. And the MAGIC image above 1 TeV is also adopted as the spatial template of
Src T to explore whether or not the gamma-ray emission above 10 GeV is a
superposition of two sources, which is similar to MAG1 and MAG2 (MAGIC
Collaboration et al., 2014). The centroids and extensions of the 2D Gaussian
and uniform disk are fitted by Fermipy
444https://fermipy.readthedocs.io/en/latest/, a PYTHON package that automates
analyses with the Fermi Science Tools (Wood et al., 2017). And the results of
the spatial analysis with the different energy ranges are listed in Table 1.
From Table 1, we can see that both the uniform disk and 2D Gaussian we
analyzed can describe the GeV emission of Src A and 4FGL J1857.9+0313c
(hereafter renamed as 4FGL J1857.9+0313e). And for Src T, a 2D Gaussian model
is the best-fit template instead of the MAGIC image, indicating that the GeV
emission from Src T is in favor of one-source scenario instead of two separate
gamma-ray sources. In the following analyses, the uniform disk is adopted as
the spatial template of Src A and 4FGL J1857.9+0313e, while Src T is described
by the 2D Gaussian model.

Figure 1: $2.5^{\circ}\times 2.5^{\circ}$ TS maps centered on HESS J1857+026.
Each pixel is $0.02^{\circ}$ and the TS maps are smoothed with a Gaussian
kernel of $\sigma=0.04^{\circ}$. The green contours indicate the 12CO($J$=1-0)
emission integrated between 78 and 90 $\rm{km}\ \rm{s}^{-1}$ by FUGIN (Umemoto
et al., 2017). The magenta contours are the H.E.S.S. observation of HESS
J1857+026 (H. E. S. S. Collaboration et al., 2018). The cyan dashed circle is
the $r_{\rm 68}$ of uniform disk of HESS J1857+026 given by FGES (Ackermann et
al., 2017). The red solid circle shows the $r_{\rm 68}$ of 2D Gaussian for
1LHAASO J1857+0245 detected by LHAASO-WCDA (Cao et al., 2023). And the blue
and white crosses are the position of PSR J1856+0245 and PSR 1857+0300,
respectively. Left: TS map for the data below 3 GeV with the diffuse
backgrounds and 4FGL-DR3 sources subtracted except for 4FGL J1857.9+0313e. The
black plus shows the position of 4FGL J1857.9+0313c as a point source given by
4FGL-DR3, and the white circle sign the $r_{\rm 68}$ of best-fit uniform disk
we analyzed here. Middle: TS map for the data below 3 GeV with the diffuse
backgrounds and 4FGL-DR3 sources subtracted except for Src A. The $r_{\rm 68}$
of the uniform disk for SrcA is described by the white circle. Right: TS map
with data above 10 GeV. The white circle marks the $r_{\rm 68}$ of the best-
fit 2D Gaussian model for SrcT.

As shown in the middle panel of Figure 1, the centroid of the gamma-ray
emission from Src A is far from the position of PSR J1856+0245/PSR 1857+0300
and is on the edge of the TeV emission of HESS J1857+026, which suggests that
there is no spatial coincidence between Src A and the pulsars or TeV emission.
For Src T, the extension of the best-fit 2D-Gaussian template is much smaller
than the result given by FGES, which could be attributable to the improvement
of the Galactic diffuse background or the newly detected gamma-ray sources.
Considering the comparable extension and the spatial coincidence between Src T
and H.E.S.S. observation (H. E. S. S. Collaboration et al., 2018), we suggest
that the GeV emission of Src T has the same origin of the TeV emission. In
addition, we found that the spectral index of Src A in the energy range of 1-3
GeV is $\Gamma=2.50\pm 0.34$, while the spectral index of Src T for the data
above 10 GeV is $\Gamma=1.84\pm 0.10$. Considering the distinct GeV spectra
and morphologies, the origins of Src A and Src T are probably different. And
both Src A and Src T are considered in the following analysis.

Table 1: Spatial Properties for the GeV Emission in the Direction of HESS
J1857+026

Energy Range | Spatial Template | Sources | R.A. | Dec. | $r_{\rm 68}\tablenotemark{$\star$}$ |  $\Delta$log L $\dagger$$\dagger$footnotemark:  | d.o.f
---|---|---|---|---|---|---|---
|  |  | (deg) | (deg) | (deg) |  |
1-3 GeV | Point+Disk (4FGL) | 4FGL J1857.9+0313c | $284.478$ | $3.219$ | $-$ | 0 | 10
|  | Src A | $284.449$ | $2.774$ | $0.50$
| Two Point | 4FGL J1857.9+0313c | $284.665$ | $3.260$ | $-$ | 6 | 9
|  | Src A | $284.038$ | $2.589$ | $-$
| Point+Disk | 4FGL J1857.9+0313c | $284.665$ | $3.260$ | $-$ | 26 | 10
|  | Src A | $284.169\pm 0.036$ | $2.504\pm 0.047$ | $0.40\pm 0.03$
| Point+Gaussian | 4FGL J1857.9+0313c | $284.665$ | $3.260$ | $-$ | 25 | 10
|  | Src A | $284.178\pm 0.045$ | $2.503\pm 0.064$ | $0.51\pm 0.07$
| Two Disks | 4FGL J1857.9+0313c | $284.641\pm 0.035$ | $3.274\pm 0.039$ | $0.28\pm 0.04$ | 34 | 11
|  | Src A | $284.079\pm 0.043$ | $2.467\pm 0.041$ | $0.35\pm 0.03$
| Disk+Gaussian | 4FGL J1857.9+0313c | $284.641\pm 0.035$ | $3.274\pm 0.039$ | $0.28\pm 0.04$ | 32 | 11
|  | Src A | $284.028\pm 0.046$ | $2.498\pm 0.041$ | $0.31\pm 0.05$
| Gaussian+Disk | 4FGL J1857.9+0313c | $284.641\pm 0.039$ | $3.273\pm 0.038$ | $0.29\pm 0.06$ | 34 | 11
|  | Src A | $284.079\pm 0.043$ | $2.467\pm 0.041$ | $0.35\pm 0.03$
10-500 GeV | Disk (4FGL) | Src T | $284.449$ | $2.774$ | $0.50$ | $0$ | $5$
| Disk | Src T | $284.361\pm 0.030$ | $2.797\pm 0.029$ | $0.36\pm 0.02$ | $15$ | $5$
| Gaussian | Src T | $284.341\pm 0.031$ | $2.757\pm 0.035$ | $0.40\pm 0.04$ | $20$ | $5$
| H.E.S.S. Image | Src T | $-$ | $-$ | $-$ | $15$ | $2$
| MAGIC Image | Src T | $-$ | $-$ | $-$ | $14$ | $2$

⋆⋆footnotetext: $r_{\rm 68}$ is the 68% containment radius, where $r_{\rm
68}=1.51\sigma$ for 2D Gaussian model and $r_{\rm 68}=0.82\sigma$ for uniform
disk model (Lande et al., 2012).

††footnotetext: Calculated with respect to the spatial model used in 4FGL-DR3.

####  2.1.2 Energy-Dependent Extension Analysis of Src T

To further explore the energy-dependent behavior of Src T, we performed the
extension analyses in the energy ranges of 10-40 GeV, 40-140 GeV and 140-500
GeV, respectively. For Src T, the 2D-Gaussian template is adopted, while the
centroid and extension in each energy band are refitted with Fermipy. The
results with the different energy ranges are listed in Table 2, and the
corresponding TS maps are presented in Figure 2.

The energy-dependent analysis show that the extension of Src T in the energy
range of 10-40 GeV is larger than that of higher energy bands, with $r_{\rm
68}$ varying from $0.44^{\circ}$ below 40 GeV to $\sim 0.30^{\circ}$ above 140
GeV. Similar phenomenons are also observed in the typical PWN HESS J1825-137
at TeV and GeV energies (Aharonian et al., 2006; H. E. S. S. Collaboration et
al., 2019; Principe et al., 2020). In addition, it should be noted that the
gamma-ray emission regions in the energy range of 40-140 GeV and 140 GeV-500
GeV seem to be different, as shown in the middle and right panels of Figure 2.
However, such phenomena could be attributed to the limited statistics of
gamma-ray photons, and more observational data will be helpful to explore the
energy-dependent behavior of Src T.

Table 2: Extension Measurements of Src T in the Different Energy Ranges Energy Range | R.A. | Dec. | $r_{68}$ | TS
---|---|---|---|---
| (deg) | (deg) | (deg) |
10-40 GeV | $284.391\pm 0.056$ | $2.846\pm 0.070$ | $0.44\pm 0.07$ | $53$
40-140 GeV | $284.293\pm 0.049$ | $2.787\pm 0.053$ | $0.31\pm 0.05$ | $52$
140-500 GeV | $284.347\pm 0.062$ | $2.769\pm 0.059$ | $0.30\pm 0.06$ | $40$

Figure 2: $2^{\circ}\times 2^{\circ}$ TS maps centered on HESS J1857+026 in
the energy ranges of 10-40 GeV (left), 40-140 GeV (middle) and 140-500 GeV
(right) . The white circles indicate the $r_{68}$ of 2D Gaussian template for
the different energy bands. The magenta contours are H.E.S.S. observation of
HESS J1857+026 (H. E. S. S. Collaboration et al., 2018). And the position of
PSR J1856+0245 and PSR 1857+0300 are shown as the green and the cyan crosses,
respectively.

###  2.2 Spectral Analysis

The spectral analysis was performed in the energy range of 500 MeV - 500 GeV.
And similar to the spatial analysis, the events of “PSF3” type with a better
angular resolution were selected for the data below 3 GeV. During the analysis
process, the summed likelihood analysis method was adopted, and both Src A and
Src T were included in the model with their spatial templates given in Sect.
2.1.1. The spectra of Src A, Src T and 4FGL J1857.9+0313e are adopted to be
the power law models. And we also tested the GeV spectral curvature for each
source by adopting the log-parabola spectrum, while no significant improvement
was obtained compared with the power law models. The global fitting gives a
hard spectrum for Src T with the photon index of $\Gamma=1.74\pm 0.07$, which
is comparable to the result in Rousseau et al. (2012). And the integral photon
flux is estimated to be $(5.04\pm 1.08)\times 10^{-9}$ ph cm-2 s-1. While the
spectrum of Src A is very soft, and the photon index is fitted to be
$\Gamma=2.73\pm 0.10$, with the integral photon flux of $(1.98\pm 0.19)\times
10^{-8}$ ph cm-2 s-1. And the best-fit photon index of 4FGL J1857.9+0313e is
$\Gamma=2.55\pm 0.07$ with the integral photon flux of $(1.41\pm 0.13)\times
10^{-8}$ ph cm-2 s-1.

To study the spectral energy distributions (SEDs) of these three sources, the
data are divided into nine logarithmically space energy bins. And the summed
likelihood analysis is repeated in each energy bin, with only the
normalizations of sources located in ROI and the diffuse backgrounds in the
model left free, while the spectral parameters are fixed to be the global-
fitting values. For the energy bin with TS value lower than 4.0, an upper
limit with 95% confidence level is calculated. And the SEDs are shown in
Figure 3. The GeV spectrum of Src T could connect with the TeV SED of HESS
J1857+026 smoothly, which suggests the same physical origin.

Figure 3:  Left: The GeV SEDs of SrcA and SrcT are marked by the brown and
black circles, respectively. The red and purple circles indicate the H.E.S.S.
observations by Aharonian et al. (2008) and H. E. S. S. Collaboration et al.
(2018), respectively. And the blue circles show the MAGIC observation in MAGIC
Collaboration et al. (2014). And the cyan butterfly indicates the global power
law spectrum of 1LHAASO J1857+0245 detected by LHAASO-WCDA (Cao et al., 2023).
Right: The GeV SED of 4FGL J1857.9+0313e. The gray solid and dashed lines show
the best-fit power law spectrum and its 1$\sigma$ statistic error. And the
hadronic models with $E_{\rm p,cut}$ = 3 PeV and $E_{\rm p,cut}$ = 1 TeV for
4FGL J1857.9+0313e are shown as the red solid and dashed lines, respectively.

##  3 Discussion

Figure 4: The broadband SED of HESS J1857+026 with the leptonic model and SrcA
with hadronic model. The radio and X-ray flux upper limits are shown as the
green and pink dots (Rousseau et al., 2012; Petriella et al., 2021). The blue
dashed line shows the synchrotron component for HESS J1857+026. The red and
green dotted lines represent IC scattering of different seed photons. And the
black solid curve is the sum of different leptonic radiation components for
HESS J1857+026. The gray solid and dashed lines indicate the hadronic models
with $E_{\rm p,cut}$ = 3 PeV and $E_{\rm p,cut}$ = 1 TeV for SrcA.

The spatial and spectral data analyses above reveal that the diffuse GeV
gamma-ray emission around HESS J1857+026 could be distinguished into two
separate extended components: one of them, namely Src A, has no spatial and
spectral coincidence with HESS J1857+026. And another one, namely Src T, shows
both spatial and spectral consistency with the TeV measurement of HESS
J1857+026, which supports Src T as the GeV counterpart of HESS J1857+026.

Although MAGIC observation revealed two sources in the region of HESS
J1857+026, MAG1 and MAG2, with the events above 1 TeV (MAGIC Collaboration et
al., 2014), there is no evidence show that Src T is composed by two gamma-ray
sources limited by the event statistics and the PSF of Fermi-LAT. Therefore,
we suggest that the GeV emission of HESS J1857+026 originates from a single
gamma-ray source. The right panel of Figure 1 shows that the centroid of GeV
gamma-ray emission is consistent with PSR J1856+0245 associated with MAG1, not
PSR J1857+0300 in the direction of MAG2, which supports the same origin for
HESS J1857+026 and MAG1. And for the possible origin of the gamma-ray emission
from HESS J1857+026, MAGIC Collaboration et al. (2014) suggested a PWN powered
by PSR J1856+0245 with $\dot{E}$ = $4.6\times 10^{36}$ erg s-1, and the spin-
down luminosity of PSR J1857+0300 with $\dot{E}$ $\sim 2.5\times 10^{32}\ {\rm
erg}\ {\rm s^{-1}}$ is too low to power a gamma-ray PWN (Abdo et al., 2013;
Acero et al., 2013). While Petriella et al. (2021) revealed the existence of a
superbubble with the analysis of atomic gas in this region, and suggested the
superbubble origin for the TeV emission from HESS J1857+026. However, the hard
GeV gamma-ray spectrum of HESS J1857+026 makes it to be different from the
typical superbubble with an index of $\sim$ 2.2, e.g. Cygnus Cocoon (Ackermann
et al., 2011; Aharonian et al., 2019), but is similar to the typical PWNe,
e.g. MSH 15-52 (Abdo et al., 2010a), HESS J1825-137 (Grondin et al., 2011). In
addition, the GeV gamma-ray emission of HESS J1857+026 also shows the energy-
dependent morphology, with the emission radius varying from $0.44^{\circ}$
below 40 GeV to $\sim 0.30^{\circ}$ above 140 GeV. And the centroid of the GeV
emission moves towards PSR J1856+0245 with increasing energies. Such
characteristics are also detected in the PWN HESS J1825-137 and HESS J1303-631
(Aharonian et al., 2006; H. E. S. S. Collaboration et al., 2019, 2012;
Principe et al., 2020). All these evidences above support that the gamma-ray
emission of HESS J1857+026 could originate from a PWN associated with PSR
J1856+0245.

For PWNe, the emission from radio to X-rays are normally produced by the
synchrotron emission, whereas the gamma-ray emission are explained by the IC
scattering process (leptonic process). Here, a simple one-zone leptonic model
is applied for HESS J1857+026. The IC background photon fields include the
cosmic microwave background (CMB) and infrared (IR) photon from dust with the
temperature of T $\sim$ 30 K and density of 1 eV cm-3 (Rousseau et al., 2012).
And the distance of HESS J1857+026 is adopted to be 6.3 kpc, derived from the
dispersion measurement of PSR J1856+0245 (Yao et al., 2017). Considering the
absence of radio detection for the PWN, we used the sensitivity of VLA in
Petriella et al. (2021) to calculate the upper limits of radio flux by
assuming that the spatial extension of the radio PWN is same to the that of
the gamma-ray emission with a radius of $0.4^{\circ}$. And the radio upper
limits are estimated to be 9.2$\times$10-14 erg cm-2 s-1 at 1.5 GHz and
1.5$\times$10-13 erg cm-2 s-1 at 5 GHz, respectively. The electron spectrum is
assumed to be a broken power law with an exponential cutoff in the form of
$\frac{dN_{e}}{dE}\propto\frac{(E/E_{br})^{-\alpha_{1}}}{1+(E/E_{br})^{\alpha_{2}-\alpha_{1}}}\rm{exp}\left(-\frac{E}{E_{e,cut}}\right)$
(Grondin et al., 2011; Xin et al., 2018), where $E_{\rm br}$, $E_{\rm e,cut}$,
$\alpha_{1}$ and $\alpha_{2}$ are the break energy, cutoff energy, and indices
of electron spectrum, respectively. The model fitting is performed using the
naima package (Zabalza, 2015).

As shown in Figure 4, the gamma-ray spectrum of HESS J1857+026 can be
reproduced with $\alpha_{1}\sim 2.2$, $\alpha_{2}\sim 3.2$, a break energy
with $E_{\rm br}\sim 3.5$ TeV, and a cutoff energy with $E_{\rm e,cut}\sim 70$
TeV, respectively. The total energy of electrons above 1 GeV, $W_{e}$, is
calculated to be $\sim 1.1\times 10^{49}(d/6.3\ {\rm kpc})^{2}$ erg. The
cooling timescale at break energy is estimated to be $t_{\rm cool}\approx
56({\rm E}/3.5\ {\rm TeV})^{-1}[(U_{\rm ph}+U_{\rm B})/1.66\ {\rm eV}\ {\rm
cm}^{-3}]^{-1}$ kyr, here $U_{\rm ph}$ = 1.26 eV/${\rm cm}^{3}$ and $U_{\rm
B}=B^{2}/8\pi$. And such value is about two times larger than the
characteristic age of PSR J1856+0245, which suggests that the break structure
could be an intrinsic characteristic of injected electronic spectrum, not be
produced by the cooling effect (Gaensler & Slane, 2006). With the electron
spectrum, the radio and X-ray upper limits constrain the magnetic field
strength to be lower than $\sim$4$\mu$G, which is consistent with the typical
values for gamma-ray PWNe (e.g. Grondin et al., 2011, 2013). It should be
noted that the X-ray upper limit was calculated with a radius of $0.1^{\circ}$
by Rousseau et al. (2012). And with the same radius of $0.4^{\circ}$ for the
calculation of radio upper limits, the X-ray flux would be much larger, which
would have no affect to the SED fitting.

Along with the evolution of the PWN into the interstellar medium (ISM), the
energetic electrons could escape and transport dominated by diffusion to form
a detectable halo around the pulsar, which is defined as a pulsar halo. Such
halos were first detected around Geminga and PSR B0656+14 with the TeV gamma-
ray emission Abeysekara et al. (2017). And Di Mauro et al. (2019) also claimed
to detect the corresponding GeV gamma-ray emission around Geminga. Based on
the definition of an electron halo in Giacinti et al. (2020), namely that of
the overdensity of relativistic electrons around pulsar compared with ISM, we
estimate the electronic energy density of HESS J1857+026, $\varepsilon_{\rm
e}$, and compared it with the typical value of ISM, $\varepsilon_{\rm ISM}$ =
0.1 eV cm-3. The gamma-ray emission region of HESS J1857+026 is assumed to be
a sphere with a radius of $\sim 0.40^{\circ}$, which corresponds to a physical
radius of $\sim$ 44 pc for the distance of 6.3 kpc. Using the total energy of
electrons we derived, the electronic energy density of $\varepsilon_{\rm e}$
$\sim 0.60\ {\rm eV}\ {\rm cm}^{-3}$ is obtained, which is much larger than
that of ISM. Therefore, we suggest that the relativistic electrons are still
contained in a region energetically and dynamically dominated by the pulsar.

Figure 5: 12CO($J$=1–0) intensity maps (in the unit of K km s-1) in the
velocity ranges of 50-65 km s-1 (left) and 78-90 km s-1 (right). The green,
white and magenta circles show the GeV gamma-ray extensions of SrcA, HESS
J1857+026, and 4FGL J1857.9+0313e, respectively. The position of PSR
J1856+0245 is marked by the yellow cross. And the cyan contours indicate the
MAGIC gamma-ray flux map above 1 TeV (MAGIC Collaboration et al., 2014).

For the other two extended gamma-ray sources, SrcA and 4FGL J1857.9+0313e,
either of them is not spatially consistent with the identified counterpart in
other wavelengths. By analyzing the 12CO($J$=1-0) line data performed by the
FOREST Unbiased Galactic plane Imaging survey with the Nobeyama 45 m telescope
(FUGIN; Umemoto et al., 2017), two molecular clumps in the southwest (named as
clumpA) and northeast (named as clumpB) of HESS J1857+026 are founded in both
the velocity ranges of 50 - 65 $\rm{km}\ \rm{s}^{-1}$ and 78 - 90 $\rm{km}\
\rm{s}^{-1}$, as shown in Figure 5, which are spatially coincidence with the
GeV gamma-ray emission from SrcA and 4FGL J1857.9+0313e, respectively. The
molecular clouds in the range of 50 - 65 $\rm{km}\ \rm{s}^{-1}$ have been
revealed in MAGIC Collaboration et al. (2014) with corresponding kinetic
distance of $\sim$3.7 kpc. And Petriella et al. (2021) studies the molecular
clouds in the range of 79 - 90 $\rm{km}\ \rm{s}^{-1}$ and derived a kinetic
distance of $\sim$5.5 kpc, which is compatible, within the errors, with the DM
distance of PSR J1856+0245 ($\sim$6.3 kpc). By adopting the value of the
conversion factor of $X_{\mathrm{CO}}=2\times 10^{20}$ cm-2 (K km
s${}^{-1})^{-1}$ (Bolatto et al., 2013), we estimate the total mass contents
of clumpA and clumpB with the different distances, respectively. And the total
masses of clumpA and clumpB are estimated in the region of $0.35\arcdeg$ and
$0.28\arcdeg$ sky integration radii, corresponding to the 68% containment
radii of the extended gamma-ray emission of SrcA and 4FGL J1857.9+0313e.
respectively. For the velocity range of 50 - 65 $\rm{km}\ \rm{s}^{-1}$ with
distance of 3.7 kpc, the total mass contents of clumpA and clumpB are
calculated to be $\sim$ $2.1\times 10^{5}d_{3.7}^{2}M_{\odot}$ and $\sim$
$1.7\times 10^{5}d_{3.7}^{2}M_{\odot}$, corresponding to the average gas
number densities of $\rm{n_{\rm gas,A}}$ = 175 cm-3 and $\rm{n_{\rm gas,B}}$ =
280 cm-3 by assuming a spherical geometry of the gas distribution. For the
velocity range of 79 - 90 $\rm{km}\ \rm{s}^{-1}$ with the compatible distance
of PSR J1856+0245 with 6.3 kpc, the total mass content of clumpA and clumpB
are calculated to be $\sim$ $1.5\times 10^{5}d_{6.3}^{2}M_{\odot}$ and $\sim$
$1.6\times 10^{5}d_{6.3}^{2}M_{\odot}$. And the corresponding average gas
number densities are about $\rm{n_{\rm gas,A}}$ = 27 cm-3 and $\rm{n_{\rm
gas,B}}$ = 50 cm-3, respectively. The spatial coincidence between the extended
gamma-ray emission and the molecular gas suggest the hadronic origin for SrcA
and 4FGL J1857.9+0313e.

In the hadronic scenario, the proton spectrum is assumed to be a single power
law with an exponential cutoff in the form of $\frac{dN_{p}}{dE}\propto
E^{-\gamma}\rm{exp}\left(-\frac{E}{E_{p,cut}}\right)$. And the cutoff energy
of protons cannot be well constrained and was first adopted to be the energy
of the cosmic ray knee with $E_{\rm p,cut}$ = 3 PeV. The hadronic model for
Src A is shown as the gray solid line in Figure 4, and the corresponding
proton spectrum should be much soft with $\gamma\sim 2.8$. The total energy of
protons above 1 GeV is estimated to be $W_{p}\sim 1.6\times 10^{50}(n_{\rm
gas,A}/27\ {\rm cm}^{-3})^{-2}(d/6.3\ {\rm kpc})^{2}$ erg or $W_{p}\sim
1.6\times 10^{48}(n_{\rm gas,A}/175\ {\rm cm}^{-3})^{-2}(d/3.7\ {\rm
kpc})^{2}$ erg. For 4FGL J1857.9+0313e, the spectral index of protons with
$\gamma\sim 2.7$ and total energy of $W_{p}\sim 5.8\times 10^{49}(n_{\rm
gas,B}/50\ {\rm cm}^{-3})^{-2}(d/6.3\ {\rm kpc})^{2}$ erg or $W_{p}\sim
6.4\times 10^{47}(n_{\rm gas,B}/280\ {\rm cm}^{-3})^{-2}(d/3.7\ {\rm
kpc})^{2}$ are needed to explain the GeV gamma-ray emission, which is shown as
the red solid line in the right panel of Figure 3. In addition, considering
the fact that the cut-off energy of protons in the middle-aged SNRs is usually
much lower than PeV (Gelfand et al., 2013; Supán et al., 2023), we decreased
the cutoff energies of protons to estimate the different hadronic models. And
the allowed minimum value of cutoff energy is about 1 TeV for SrcA and 4FGL
J1857.9+0313e, with the total energies of protons of $\sim 1.5\times
10^{50}(n_{\rm gas,A}/27\ {\rm cm}^{-3})^{-2}(d/6.3\ {\rm kpc})^{2}$ erg and
$\sim 5.7\times 10^{49}(n_{\rm gas,B}/50\ {\rm cm}^{-3})^{-2}(d/6.3\ {\rm
kpc})^{2}$ erg, respectively.

For the molecular cloud in 50 -65 $\rm{km}\ \rm{s}^{-1}$ with the distance of
3.7 kpc, there is no candidates as the origin of high energy protons. While
for the molecular clouds in 79 - 90 $\rm{km}\ \rm{s}^{-1}$, the compatible
distance with PSR J1856+0245 suggests that the hypothetical SNR associated
with PSR J1856+0245 could provide enough power assuming that $\sim$10% of the
supernova kinetic energy of $\sim$1051 erg is transferred to the energy of
particles. Moreover, the soft GeV gamma-ray spectra of SrcA and 4FGL
J1857.9+0313e are also similar to that of the old SNRs interacting with
molecular clouds (e.g. IC 443 & W44 Abdo et al., 2010b; Ackermann et al.,
2013). The absence of the detection of associated SNR in this region suggests
that it may have already dissipated into the ambient gas. And further
observations, especially in the radio and X-ray bands, will be crucial to
reveal the physical origin of the gamma-ray emission in this region.

##  4 Conclusion

Using more than 13 years of Fermi-LAT observations, we studied the GeV gamma-
ray emission in the direction of HESS J1857+026, and found that the GeV
emission around HESS J1857+026 is composed of two extended gamma-ray sources:
Src A and Src T. Src T is spatially coincident with HESS J1857+026 and its
hard GeV gamma-ray spectrum could connect with the TeV SED of HESS J1857+026
smoothly, indicating that SrcT could be the GeV counterpart of HESS J1857+026.
In addition, we performed the energy-dependent analyses of the GeV gamma-ray
emission from HESS J1857+026, and its extension decreases towards higher
energies. The energy-dependent morphology and the hard GeV gamma-ray spectrum
of HESS J1857+026 make it favor the PWN origin. And one-zone leptonic model
with a broken power law electronic spectrum can well describe the multi-
wavelength data of HESS J1857+026. SrcA and another extended gamma-ray source,
4FGL J1857.9+0313e with soft GeV gamma-ray spectra, have no identified
counterparts in other wavelength. However, two molecular clumps in the
northeast and southwest of HESS J1857+026 are spatially coincidence with the
GeV gamma-ray emission from SrcA and 4FGL J1857.9+0313e, which suggest the
hadronic process for their gamma-ray emission. With the single power law model
for protons, the GeV spectra of SrcA and 4FGL J1857.9+0313e could also be
explained with the soft proton spectra. And these high energy protons could be
produced by the hypothetical SNR associated with PSR J1856+0245, which may
have already dissipated into the ambient gas.

HESS J1857+026 is one of the peculiar gamma-ray sources, which shows the
energy-dependent morphology. More detailed observations by LHAASO (Cao et al.,
2019) and CTA (Cherenkov Telescope Array Consortium et al., 2019) would be
helpful to explore the particle transport mechanisms, and the future radio and
X-ray observations are crucial to investigate the origin of the gamma-ray
emission in this region.

## Acknowledgments

We acknowledge the use of the Fermi-LAT data provided by the Fermi Science
Support Center. This work is supported by the National Natural Science
Foundation of China under the grants 12103040, 12147208 and U1931204, and the
Natural Science Foundation for Young Scholars of Sichuan Province, China (No.
2022NSFSC1808).

## References

  * Abdo et al. (2010a) Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010a, ApJ, 714, 927, doi: 10.1088/0004-637X/714/1/927
  * Abdo et al. (2010b) —. 2010b, Science, 327, 1103, doi: 10.1126/science.1182787
  * Abdo et al. (2013) Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17, doi: 10.1088/0067-0049/208/2/17
  * Abdollahi et al. (2020) Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33, doi: 10.3847/1538-4365/ab6bcb
  * Abdollahi et al. (2022) Abdollahi, S., Acero, F., Baldini, L., et al. 2022, ApJS, 260, 53, doi: 10.3847/1538-4365/ac6751
  * Abeysekara et al. (2013) Abeysekara, A. U., Alfaro, R., Alvarez, C., et al. 2013, Astroparticle Physics, 50, 26, doi: 10.1016/j.astropartphys.2013.08.002
  * Abeysekara et al. (2017) Abeysekara, A. U., Albert, A., Alfaro, R., et al. 2017, ApJ, 843, 40, doi: 10.3847/1538-4357/aa7556
  * Acero et al. (2013) Acero, F., Ackermann, M., Ajello, M., et al. 2013, ApJ, 773, 77, doi: 10.1088/0004-637X/773/1/77
  * Ackermann et al. (2011) Ackermann, M., Ajello, M., Allafort, A., et al. 2011, Science, 334, 1103, doi: 10.1126/science.1210311
  * Ackermann et al. (2013) —. 2013, Science, 339, 807, doi: 10.1126/science.1231160
  * Ackermann et al. (2017) Ackermann, M., Ajello, M., Baldini, L., et al. 2017, ApJ, 843, 139, doi: 10.3847/1538-4357/aa775a
  * Aharonian et al. (2019) Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, Nature Astronomy, 3, 561, doi: 10.1038/s41550-019-0724-0
  * Aharonian et al. (2004) Aharonian, F., Akhperjanian, A. G., Aye, K. M., et al. 2004, Astroparticle Physics, 22, 109, doi: 10.1016/j.astropartphys.2004.06.006
  * Aharonian et al. (2006) Aharonian, F., Akhperjanian, A. G., Bazer-Bachi, A. R., et al. 2006, A&A, 460, 365, doi: 10.1051/0004-6361:20065546
  * Aharonian et al. (2008) Aharonian, F., Akhperjanian, A. G., Barres de Almeida, U., et al. 2008, A&A, 477, 353, doi: 10.1051/0004-6361:20078516
  * Aleksić et al. (2016) Aleksić, J., Ansoldi, S., Antonelli, L. A., et al. 2016, Astroparticle Physics, 72, 61, doi: 10.1016/j.astropartphys.2015.04.004
  * Atwood et al. (2009) Atwood, W. B., Abdo, A. A., Ackermann, M., et al. 2009, ApJ, 697, 1071, doi: 10.1088/0004-637X/697/2/1071
  * Bolatto et al. (2013) Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207, doi: 10.1146/annurev-astro-082812-140944
  * Cao et al. (2019) Cao, Z., della Volpe, D., Liu, S., et al. 2019, arXiv e-prints, arXiv:1905.02773.  https://arxiv.org/abs/1905.02773
  * Cao et al. (2023) Cao, Z., Aharonian, F., An, Q., et al. 2023, arXiv e-prints, arXiv:2305.17030, doi: 10.48550/arXiv.2305.17030
  * Cherenkov Telescope Array Consortium et al. (2019) Cherenkov Telescope Array Consortium, Acharya, B. S., Agudo, I., et al. 2019, Science with the Cherenkov Telescope Array, doi: 10.1142/10986
  * Cordes & Lazio (2002) Cordes, J. M., & Lazio, T. J. W. 2002, arXiv e-prints, astro.  https://arxiv.org/abs/astro-ph/0207156
  * Cui et al. (2021) Cui, Y., Xin, Y., Liu, S., et al. 2021, A&A, 646, A114, doi: 10.1051/0004-6361/202039041
  * Di Mauro et al. (2019) Di Mauro, M., Manconi, S., & Donato, F. 2019, Phys. Rev. D, 100, 123015, doi: 10.1103/PhysRevD.100.123015
  * Gaensler & Slane (2006) Gaensler, B. M., & Slane, P. O. 2006, ARA&A, 44, 17, doi: 10.1146/annurev.astro.44.051905.092528
  * Gelfand et al. (2013) Gelfand, J. D., Castro, D., Slane, P. O., et al. 2013, ApJ, 777, 148, doi: 10.1088/0004-637X/777/2/148
  * Giacinti et al. (2020) Giacinti, G., Mitchell, A. M. W., López-Coto, R., et al. 2020, A&A, 636, A113, doi: 10.1051/0004-6361/201936505
  * Grondin et al. (2013) Grondin, M. H., Romani, R. W., Lemoine-Goumard, M., et al. 2013, ApJ, 774, 110, doi: 10.1088/0004-637X/774/2/110
  * Grondin et al. (2011) Grondin, M. H., Funk, S., Lemoine-Goumard, M., et al. 2011, ApJ, 738, 42, doi: 10.1088/0004-637X/738/1/42
  * H. E. S. S. Collaboration et al. (2012) H. E. S. S. Collaboration, Abramowski, A., Acero, F., et al. 2012, A&A, 548, A46, doi: 10.1051/0004-6361/201219814
  * H. E. S. S. Collaboration et al. (2018) H. E. S. S. Collaboration, Abdalla, H., Abramowski, A., et al. 2018, A&A, 612, A1, doi: 10.1051/0004-6361/201732098
  * H. E. S. S. Collaboration et al. (2019) H. E. S. S. Collaboration, Abdalla, H., Aharonian, F., et al. 2019, A&A, 621, A116, doi: 10.1051/0004-6361/201834335
  * Hessels et al. (2008) Hessels, J. W. T., Nice, D. J., Gaensler, B. M., et al. 2008, ApJ, 682, L41, doi: 10.1086/590908
  * Lande et al. (2012) Lande, J., Ackermann, M., Allafort, A., et al. 2012, ApJ, 756, 5, doi: 10.1088/0004-637X/756/1/5
  * Lyne et al. (2017) Lyne, A. G., Stappers, B. W., Bogdanov, S., et al. 2017, ApJ, 834, 137, doi: 10.3847/1538-4357/834/2/137
  * MAGIC Collaboration et al. (2014) MAGIC Collaboration, Aleksić, J., Ansoldi, S., et al. 2014, A&A, 571, A96, doi: 10.1051/0004-6361/201423517
  * Petriella et al. (2021) Petriella, A., Duvidovich, L., & Giacani, E. 2021, A&A, 652, A142, doi: 10.1051/0004-6361/202141254
  * Principe et al. (2020) Principe, G., Mitchell, A. M. W., Caroff, S., et al. 2020, A&A, 640, A76, doi: 10.1051/0004-6361/202038375
  * Rousseau et al. (2012) Rousseau, R., Grondin, M. H., Van Etten, A., et al. 2012, A&A, 544, A3, doi: 10.1051/0004-6361/201118685
  * Supán et al. (2023) Supán, L., Castelletti, G., & Lemière, A. 2023, A&A, 679, A22, doi: 10.1051/0004-6361/202347305
  * Umemoto et al. (2017) Umemoto, T., Minamidani, T., Kuno, N., et al. 2017, PASJ, 69, 78, doi: 10.1093/pasj/psx061
  * Wood et al. (2017) Wood, M., Caputo, R., Charles, E., et al. 2017, in International Cosmic Ray Conference, Vol. 301, 35th International Cosmic Ray Conference (ICRC2017), 824.  https://arxiv.org/abs/1707.09551
  * Xin et al. (2018) Xin, Y.-L., Liao, N.-H., Guo, X.-L., et al. 2018, ApJ, 867, 55, doi: 10.3847/1538-4357/aae313
  * Yao et al. (2017) Yao, J. M., Manchester, R. N., & Wang, N. 2017, ApJ, 835, 29, doi: 10.3847/1538-4357/835/1/29
  * Zabalza (2015) Zabalza, V. 2015, in International Cosmic Ray Conference, Vol. 34, 34th International Cosmic Ray Conference (ICRC2015), 922.  https://arxiv.org/abs/1509.03319
  * Zhang et al. (2022) Zhang, X., Chen, Y., Zheng, F.-x., et al. 2022, ApJ, 931, 128, doi: 10.3847/1538-4357/ac6957


#  At First Sight! Zero-Shot Classification of Astronomical Images with Large
Multimodal Models

Dimitrios Tanoglidis Bhuvnesh Jain Department of Physics and Astronomy,
University of Pennsylvania, Philadelphia, PA 19104, USA

###### Abstract

Vision-Language multimodal Models (VLMs) offer the possibility for zero-shot
classification in astronomy: i.e. classification via natural language prompts,
with no training. We investigate two models, GPT-4o and LLaVA-NeXT, for zero-
shot classification of low-surface brightness galaxies and artifacts, as well
as morphological classification of galaxies. We show that with natural
language prompts these models achieved significant accuracy (above 80 percent
typically) without additional training/fine tuning. We discuss areas that
require improvement, especially for LLaVA-NeXT, which is an open source model.
Our findings aim to motivate the astronomical community to consider VLMs as a
powerful tool for both research and pedagogy, with the prospect that future
custom-built or fine-tuned models could perform better.

Computational methods (1965) — Astroinformatics (78) — Interdisciplinary
astronomy(804)

††software: Open GPT-4o API (OpenAI et al., 2023), Transformers (Wolf et al.,
2019), bitsandbytes (Dettmers, 2021)

##  1 Introduction

Multi-modal foundation models have the potential to impact diverse scientific
fields. In astronomy, the vast and varied forms of data offer opportunities.
From images taken at different wavelengths; spectra and time domain
observations; and in some limited cases language labels. Given these potential
opportunities, there have been a few attempts at multi-modal investigations in
astronomy (e.g., Parker et al., 2023; Mishra-Sharma et al., 2024)

Despite the recent impact of deep learning in astronomy, especially in image
classification tasks, a significant challenge remains: the lack of large,
high-quality, labeled datasets necessary for the training of supervised deep
learning models. For each specific task (e.g., morphological classification of
galaxies) a separate model has to be trained using labeled data for that task.
Furthermore, data obtained through one instrument or telescope cannot be
readily used to train a model to perform the same task using data from another
instrument. In such cases, domain adaptation techniques should be employed
(e.g., Ćiprijanović et al., 2022).

Foundation models (Bommasani et al., 2021), on the other hand, are models
trained on broad datasets, that can subsequently be adapted in a range of
downstream tasks, such as classification. The case where the model can
categorize objects or classes that have never been seen before, is referred to
as zero-shot learning. In this research note, we explore the potential of
large, multimodal, Vision-Language models (VLMs; see for example Bordes et al.
(2024)) to perform zero-shot classification of astronomical images, by
leveraging natural language descriptions of the different object classes. Our
results, using generalist models that have not been fine-tuned on astronomical
data, serve as a proof of concept for the potential of VLMs in astronomy.

##  2 METHODS

Figure 1: The prompts used to instruct the multimodal models on the two
classification tasks (upper part), as well as the resulting confusion
matrices, obtained using the GPT-4o (middle part) and LLaVa-NeXT (lower part)
models. The confusion matrices compare the predicted vs. the true labels for
the example images in our test datasets. The numbers in parentheses correspond
to the fraction of images with true label in class $x$ that were predicted to
belong to class $y$.

We select two problems as examples for our study: the classification of Low-
Surface-Brightness galaxies (LSBGs) vs artifacts, as described in Tanoglidis
et al. (2021), and the morphological classification of galaxies into four
categories (smooth and round, smooth and cigar-shaped, edge-on disk, and
unbarred spiral) using the
GalaxyMNIST111https://github.com/mwalmsley/galaxy_mnist dataset (Walmsley et
al., 2022).

We use the state-of-the-art, multimodal, GPT-4o model (OpenAI et al., 2023),
accessed through the OpenAI API. For comparison, we also employ the open-
source, language-vision, LLaVA-Next model (Liu et al., 2024), accessed through
the Transformers (Wolf et al., 2019) library from Hugging Face. Specifically,
we use the version with Mistral-7B (Jiang et al., 2023) as the LLM backbone,
and we perform 4-bit quantization to reduce the memory requirements when
loading the model. Both models can process both textual and image inputs
(prompts) by creating embeddings for the different modalities in a common
embedding space (for this technique see, for example, Radford et al. (2021),
and Mishra-Sharma et al. (2024) for a recent application in astrophysics).

We want to investigate how natural language descriptions of the different
classes can instruct an VLM to accurately classify astronomical images,
without further training. For this reason, we provide the exact prompts used
in our analysis – see the upper part of Fig. 1. For each one of the two
classification problems, we select a sample of 500 images; we run each model
in a loop and we ask it to output only a single number that corresponds to the
class the image belongs. The code we used for our experiments, can be found
at: https://github.com/dtanoglidis/zero-shot-astro.

##  3 RESULTS AND DISCUSSION

The results of the zero-shot classification experiments for the two problems
described in the previous section are presented in Fig. 1. In the middle
(lower) panels we show the confusion matrices (that contrast the predicted vs
the actual labels for the sample of images we classified) from the GPT-4o
(LlaVA-NeXT) model. Classifying LSBGs vs artifacts using GPT-4o reaches an
accuracy of $\sim 83\%$; this can be compared to the $\sim 92\%$ presented in
(Tanoglidis et al., 2021) using custom-trained convolutional neural networks
(CNNs) on a training set of 40k examples. Despite the lower performance, we
emphasize that this was achieved using the prompt descriptions alone, without
any further re-training (or the associated coding and choice of CNN
architecture). The overall accuracy in the morphological classification
problem is lower ($\sim 67\%$), but this is driven mainly by the confusion
between cigar-shaped and edge-on galaxies; visual inspection of some examples
in Fig. 3 shows that it is challenging to distinguish between these two
categories. If these two categories are combined into a single one, the
overall accuracy rises to $\sim 86\%$.

It is important to note that the classification is not a result of
memorization (for this problem, see e.g., Hartmann et al. (2023)) of data
present in the training set: we checked that without providing descriptions of
LSBGs and artifacts, GPT-4o tended to classify every image as belonging to a
galaxy.

Despite the significantly lower performance of LlaVA-NeXT on both problems
(see comparisons between the LLaVA and GPT-4 families at https://llava-
vl.github.io/), we chose to present the results from LLaVA-NeXT since it is an
open-source model and can be fine-tuned (e.g. Mangrulkar et al. (2022)) in the
future for classification tasks on astronomical images and text (see for
example the AstroLLaVA model, currently accessible as a Hugging Face space at:
https://huggingface.co/spaces/universeTBD/astroLLaVA). Another technique that
may improve the performance of VLMs is the inclusion of image examples from
each category in the prompt (few-shot learning).

This research note aims to demonstrate to the astronomy community the power of
multimodal models, particularly Vision-Language Models (VLMs), in addressing
complex classification tasks. While our current results are not yet
competitive with those achieved by supervised models, they highlight the
significant potential of VLMs for various astronomical applications (for a
discussion of different applications and tasks that LVMs can perform, see, for
example Yang et al. (2023)). With ongoing advancements in generative AI, fine-
tuning or developing custom-built models with increased astronomical domain
knowledge, we anticipate that LVMs will play a crucial role in automating and
improving astronomical data analysis in the near future.

Furthermore, applications for pedagogical exercises are promising as well. At
both the high school and undergraduate level, analysis of images from Hubble,
Webb, Subaru and other telescopes can make for appealing exercises. A
multimodal model can provide a powerful first layer to give students easy-to-
use images from full field data, so students can then carry out scientific
analyses such as studies of spiral arms or bars in disk galaxies.

We thank Marc Huertas-Company, Josh Nguyen and Helen Qu for helpful
discussions, Lily Noyes for feedback on explorations of VLMs, and the
UniverseTBD team (https://universetbd.org/) for information about their
upcoming AstroLLaVA model. In this appendix we present example images from the
different classes considered in the classification problems mentioned in the
main text.

Figure 2: Examples of LSBG and artifact images.

Figure 3: Examples of images from the four morphological categories in the
GalaxyMNIST dataset. Notice the visual similarity between the smooth cigar-
shaped and the edge-on disk galaxies.

## References

  * Bommasani et al. (2021) Bommasani, R., Hudson, D. A., Adeli, E., et al. 2021, arXiv e-prints, arXiv:2108.07258, doi: 10.48550/arXiv.2108.07258
  * Bordes et al. (2024) Bordes, F., Yuanzhe Pang, R., Ajay, A., et al. 2024, arXiv e-prints, arXiv:2405.17247, doi: 10.48550/arXiv.2405.17247
  * Ćiprijanović et al. (2022) Ćiprijanović, A., Lewis, A., Pedro, K., et al. 2022, arXiv e-prints, arXiv:2211.00677, doi: 10.48550/arXiv.2211.00677
  * Dettmers (2021) Dettmers, T. 2021, GitHub - TimDettmers/bitsandbytes: Accessible large language models via k-bit quantization for PyTorch. — github.com, https://github.com/TimDettmers/bitsandbytes
  * Hartmann et al. (2023) Hartmann, V., Suri, A., Bindschaedler, V., et al. 2023, arXiv e-prints, arXiv:2310.18362, doi: 10.48550/arXiv.2310.18362
  * Jiang et al. (2023) Jiang, A. Q., Sablayrolles, A., Mensch, A., et al. 2023, arXiv e-prints, arXiv:2310.06825, doi: 10.48550/arXiv.2310.06825
  * Liu et al. (2024) Liu, H., Li, C., Li, Y., et al. 2024, LLaVA-NeXT: Improved reasoning, OCR, and world knowledge.  https://llava-vl.github.io/blog/2024-01-30-llava-next/
  * Mangrulkar et al. (2022) Mangrulkar, S., Gugger, S., Debut, L., et al. 2022, PEFT: State-of-the-art Parameter-Efficient Fine-Tuning methods, https://github.com/huggingface/peft
  * Mishra-Sharma et al. (2024) Mishra-Sharma, S., Song, Y., & Thaler, J. 2024, arXiv e-prints, arXiv:2403.08851, doi: 10.48550/arXiv.2403.08851
  * OpenAI et al. (2023) OpenAI, Achiam, J., Adler, S., et al. 2023, arXiv e-prints, arXiv:2303.08774, doi: 10.48550/arXiv.2303.08774
  * Parker et al. (2023) Parker, L., Lanusse, F., Golkar, S., et al. 2023, arXiv e-prints, arXiv:2310.03024, doi: 10.48550/arXiv.2310.03024
  * Radford et al. (2021) Radford, A., Kim, J. W., Hallacy, C., et al. 2021, arXiv e-prints, arXiv:2103.00020, doi: 10.48550/arXiv.2103.00020
  * Tanoglidis et al. (2021) Tanoglidis, D., Ćiprijanović, A., & Drlica-Wagner, A. 2021, Astronomy and Computing, 35, 100469, doi: 10.1016/j.ascom.2021.100469
  * Walmsley et al. (2022) Walmsley, M., Lintott, C., Géron, T., et al. 2022, MNRAS, 509, 3966, doi: 10.1093/mnras/stab2093
  * Wolf et al. (2019) Wolf, T., Debut, L., Sanh, V., et al. 2019, arXiv e-prints, arXiv:1910.03771, doi: 10.48550/arXiv.1910.03771
  * Yang et al. (2023) Yang, Z., Li, L., Lin, K., et al. 2023, arXiv e-prints, arXiv:2309.17421, doi: 10.48550/arXiv.2309.17421