11institutetext: 1 Max Planck Institute for Astronomy, Königstuhl 17, 69117
Heidelberg, Germany, 11email<EMAIL_ADDRESS>
2 Leiden Observatory, Leiden University, P.O. Box 9513, NL 2300 RA Leiden, The
Netherlands
3 European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching
bei München, Germany
4 Max-Planck-Institut for Extraterrestrial Physics, Gießenbachstrasse 1, 85748
Garching
5 School of Cosmic Physics, Dublin Institute for Advanced Studies, 31
Fitzwilliam Place, Dublin 2, Ireland
6 UK Astronomy Technology Centre, Royal Observatory Edinburgh, Blackford Hill,
Edinburgh EH9 3HJ, UK
7 INAF-Osservatorio Astronomico di Capodimonte, Salita Moiariello 16, I-80131
Napoli, Italy
8 Laboratory for Astrophysics, Leiden Observatory, Leiden University, P.O. Box
9513, NL 2300 RA Leiden, The Netherlands
9 Chalmers University of Technology, Onsala Space Observatory, 439 92 Onsala,
Sweden
10 Instituut voor Sterrenkunde, KU Leuven, Celestijnenlaan 200D, Bus-2410,
3000 Leuven, Belgium
11 Centro de Astrobiologıa (CAB, CSIC-INTA), Carretera de Ajalvir, E-28850
Torrejon de Ardoz, Madrid, Spain
12 DTU Space, Technical University of Denmark. Building 328, Elektrovej, 2800
Kgs. Lyngby, Denmark
13 Dept. of Astrophysics, University of Vienna, Türkenschanzstr. 17, A-1180
Vienna, Austria
14 ETH Zürich, Institute for Particle Physics and Astrophysics, Wolfgang-
Pauli-Str. 27, 8093 Zürich, Switzerland
15 Universite Paris-Saclay, Universite de Paris, CEA, CNRS, AIM, F-91191 Gif-
sur-Yvette, France
16 Department of Astronomy, Oskar Klein Centre; Stockholm University; SE-106
91 Stockholm, Sweden

# JOYS: JWST Observations of Young protoStars

Outflows and accretion in the high-mass star-forming region IRAS 23385+6053

H. Beuther 11 E. F. van Dishoeck 22 L. Tychoniec 33 C. Gieser 44 P.J. Kavanagh
55 G. Perotti 11 M.L. van Gelder 22 P. Klaassen 66 A. Caratti o Garatti 77 L.
Francis 22 W. R. M. Rocha 2288 K. Slavicinska 2288 T. Ray 55 K. Justtanont 99
H. Linnartz 88 C. Waelkens 1010 L. Colina 1111 T. Greve 1212 M. Güdel
1313111414 T. Henning 11 P.-O. Lagage 1515 B. Vandenbussche 1010 G. Östlin
1616 G. Wright 66

(Version of )

###### Abstract

Context. Understanding the earliest stages of star formation, and setting that
into context with the general cycle of matter in the interstellar medium, is a
central aspect of research with the James Webb Space Telescope (JWST).

Aims. The JWST program JOYS (JWST Observations of Young protoStars) aims at
characterizing the physical and chemical properties of young high- and low-
mass star-forming regions, in particular the unique mid-infrared diagnostics
of the warmer gas and solid-state components. We present early results from
the high-mass star formation region IRAS 23385+6053.

Methods. The JOYS program uses the Mid-Infrared Instrument (MIRI) Medium
Resolution Spectrometer (MRS) with its Integral Field Unit (IFU) to
investigate a sample of high- and low-mass star-forming protostellar systems.

Results. The full 5 to 28 $\mu$m MIRI MRS spectrum of IRAS 23385+6053 shows a
plethora of interesting features. While the general spectrum is typical for an
embedded protostar, we see many atomic and molecular gas lines boosted by the
higher spectral resolution and sensitivity compared to previous space
missions. Furthermore, ice and dust absorption features are also present.
Here, we focus on the continuum emission, outflow tracers like the
H2(0–0)S(7), [FeII](${}^{4}F_{9/2}$–${}^{6}D_{9/2}$) and
[NeII](${}^{2}P_{1/2}-^{2}P_{3/2}$) lines as well as the potential accretion
tracer Humphreys $\alpha$ Hi(7–6). The short-wavelength MIRI data resolve two
continuum sources A and B, where mid-infrared source A is associated with the
main mm continuum peak. The combination of mid-infrared and mm data reveals a
young cluster in its making. Combining the mid-infrared outflow tracer H2,
[FeII] and [NeII] with mm SiO data shows a complex interplay of at least three
molecular outflows driven by protostars in the forming cluster. Furthermore,
the Humphreys $\alpha$ line is detected at a 3–4$\sigma$ level towards the
mid-infrared sources A and B. Following Rigliaco et al. (2015), one can
roughly estimate both accretion luminosities and corresponding accretion rates
between $\sim$2.6$\times 10^{-6}$ and $\sim$0.9$\times 10^{-4}$
$M_{\odot}yr^{-1}$. This is discussed in the context of the observed outflow
rates.

Conclusions. The analysis of the MIRI MRS observations for this young high-
mass star-forming region reveals connected outflow and accretion signatures.
Furthermore, they outline the enormous potential of JWST to boost our
understanding of the physical and chemical processes during star formation.

###### Key Words.:

Stars: formation – ISM: clouds – ISM: individual objects: IRAS23385+6053 –
Stars: jets – Stars: massive

##  1 Introduction

Figure 1: Overview of the region around IRAS 23385+6053. The left panel shows
a large-scale overview while the middle and right panels present zoom-ins. The
color-scale in the left and middle panels show the 12 $\mu$m emission from
WISE at an angular resolution of $6.5^{\prime\prime}$ (Wright et al., 2010).
The little squares outline the field of view of the MIRI IFU mosaic in channel
1 (each square $3.7^{\prime\prime}\times 3.7^{\prime\prime}$). The color-scale
in the right panel shows the new MIRI 5.2 $\mu$m continuum emission (see also
Fig. 4). The contours in the middle and right panels present the 1.3 mm
continuum emission (Beuther et al., 2018; Cesaroni et al., 2019). The contours
start at the 5$\sigma$ level of 0.55 mJy beam-1 and continue in 25$\sigma$
steps. The black and red/blue lines in the middle panel show the directions of
outflows identified in Molinari et al. (1998) and Cesaroni et al. (2019),
respectively. The right panel marks the mid-infrared sources A and B as well
as the mm sources labeled ”mm” with the corresponding source letters from
Cesaroni et al. (2019). Mid-infrared source A and mm source mmA2 are spatially
co-located. A linear scale-bar is shown in all panels.

The earliest phases of protostellar evolution are deeply embedded in their
natal cloud cores. Many important physical and chemical processes occur there,
including the infall of the envelope, the formation of disks and outflows, the
main accretion phase and growth of the final star, as well as the chemical
enrichment of the disk directly impacting planet formation. Because of the
very high extinction $A_{V}$ of several hundreds to even thousands of
magnitudes, these deeply embedded phases have to be studied at (mid-)infrared
and (sub)mm wavelengths. While (sub)mm interferometers have allowed detailed
imaging of the cold material in such star-forming cores (for a review see
Motte et al. 2018), mid-infrared continuum and line imaging of the warmer
material at sub-arcsecond spatial resolution has been barely feasible so far.
Only a few space-based studies at lower angular resolution with the Infrared
Space Observatory (ISO) and Spitzer are reported (e.g., Whittet et al. 1996;
van Dishoeck et al. 1998; Gibb et al. 2000; Molinari et al. 2008; An et al.
2011). The most comprehensive ISO legacy study on high-mass protostars is
presented in Gibb et al. (2004).

The situation has now changed dramatically with the advent of the James Webb
Space Telescope (JWST), launched Dec. 25, 2021, just little more than a year
ago. Scientific data delivery began only in July 2022, and we present the
first results from the European MIRI (Mid-Infrared Instrument) Guaranteed Time
Program JOYS (JWST Observations of Young protoStars, PI: E. F. van Dishoeck,
pid 1290).

The JOYS program observes a total of about two dozen low- to high-mass star-
forming regions with the MIRI Medium Resolution Spectrometer (MRS) and its
Integral Field Unit (IFU) to address a broad range of scientific topics.
Evolutionary stages from very young embedded protostars within infrared dark
clouds to Class 0/I protostars as well as high-mass protostellar objects and
hot molecular cores are being covered by the JOYS project. The spatial
resolution of the MIRI instrument of $\sim 0.19^{\prime\prime}$, $\sim
0.58^{\prime\prime}$ and $\sim 0.96^{\prime\prime}$ at 5, 15 and 25 $\mu$m,
respectively, corresponds, at typical distances of the low-mass cores of
$\sim$150 pc, to linear resolution elements of $\sim$29, $\sim$87 and
$\sim$144 au, respectively. For high-mass regions at more typical distances of
$\sim$3 kpc, the linear resolution corresponds roughly to 570, 1740 and 2880
au.

Figure 2: Infrared K-band data are shown in color (taken with the Palomar
60-inch telescope and the NICMOS-3 array, Faustini et al. 2009). The contours
outline the 1.3 mm emission (Beuther et al., 2018; Cesaroni et al., 2019).
Contour levels start at 5$\sigma$ (0.55 mJy beam-1) and continue in 25$\sigma$
steps. The white full and dashed lines show the directions of outflows
identified in Molinari et al. (1998) and Cesaroni et al. (2019), respectively.

We are using a set of key diagnostics in the mid-infrared band between 5 and
28 $\mu$m such as the continuum spectral energy distribution (SED), a range of
atomic and molecular gas lines as well as several solid state bands. More
details about the JOYS project will be presented in a forthcoming publication
by van Dishoeck et al. (in prep.).

Figure 3: Full MIRI spectrum extracted towards the combined sources A and B
(central position: R.A.(J2000.0) 23:40:54.49, Decl. 61:10:27.4) within a
diameter of $2.4^{\prime\prime}$). The color coding separates the different
sub-bands from channel 1 short to channel 4 long. Prominent strong atomic and
molecular lines are labeled (see Table 1). The four lines mainly discussed in
this paper are highlighted by the red arrows at the bottom. The data are not
background-subtracted. The errorbars at 19 and 21 $\mu$m show the approximate
$1\sigma$ rms below and above 20 $\mu$m of $\sim$0.5 and $\sim$2 mJy (note the
log-scale), respectively (Sect. 2).

Here, we focus on analysing and interpreting JWST data of the first source
fully observed within the JOYS program, the high-mass protostellar object IRAS
23385+6053. This high-mass star-forming region was originally identified by
IRAS color selection and association with H2O maser emission (Casoli et al.,
1986; Cesaroni et al., 1988; Wouterloot & Brand, 1989; Palla et al., 1991).
Early interferometric investigations revealed that the region is embedded in a
comparatively dark cloud core surrounded by larger-scale mid-infrared emission
(Molinari et al., 1998). Figure 1 presents an overview of the region at 12
$\mu$m (left panel, from the Wide-field Infrared Survey Explorer – WISE –
mission at $6.5^{\prime\prime}$ resolution, Wright et al. 2010) and 1.3 mm
emission (right panel, Beuther et al. 2018, Cesaroni et al. 2019) showing the
large-scale extended mid-infrared emission that weakens at the positions of
the central bright mm emission. The outline of the 4-field MIRI-IFU coverage
is shown there as well. The $v_{\rm lsr}$ of the region is $-50.2$ km s-1. At
a kinematic distance of 4.9 kpc, the total luminosity measured with the IRAS
mission accounts for $1.6\times 10^{4}$ $L_{\odot}$ (Molinari et al., 1998).
However, higher-spatial-resolution Spitzer/MIPSGAL observations revealed that
only $\sim$3$\times 10^{3}$ $L_{\odot}$ are associated with the central core,
consistent with a high-mass protostar in the main accretion phase (Molinari et
al., 2008). Infrared imaging by Faustini et al. (2009) exhibits an infrared
cluster in the surroundings of the region. A weak K-band source is detected in
their data toward the very center (Fig. 2). An H2 emission-line study at 2.12
$\mu$m revealed shocked H2 emission in the immediate surroundings of the
central mm-core (Wolf-Chase et al., 2012). CH3OH maser emission at 44 and 95
GHz confirm shock emission in IRAS 23385+6053 (Kurtz et al., 2004; Wolf-Chase
et al., 2012). The most recent high-spatial resolution mm observations
($\sim$0.4′′) were conducted within the IRAM NOEMA (Northern Extended
Millimeter Array) large program CORE that studied the physical and chemical
properties of 20 high-mass star-forming regions (Beuther et al., 2018; Gieser
et al., 2021). A detailed case study of IRAS 23385+6053 from the CORE project
is presented in Cesaroni et al. (2019). Based on the line and continuum
emission, Cesaroni et al. (2019) identified six cores, a new outflow candidate
in the northwest-southeast direction (Fig. 1, middle panel) and a tentative
detection of a disk around a $\sim$9 $M_{\odot}$ protostar within the most
massive mm core mmA1 (Fig. 1, right panel).

This first paper in the JOYS-series presents an overview of the overall MIRI
spectrum and then focuses on the continuum emission and a few selected
molecular and atomic lines (H2(0–0)S(7),
[FeII](${}^{4}F_{9/2}$–${}^{6}D_{9/2}$) and Humphreys $\alpha$ Hi(7–6)). The
data outline the general source structure, and we investigate the outflow
properties and potential accretion signatures. Complementary studies will
discuss the overall extended emission of the diverse set of detected atomic
and H2 lines and set them in relation to the molecular gas properties derived
at mm wavelengths (Gieser et al. in prep.), the molecular emission around the
central disk candidate (Francis et al. in prep.) and the ice properties in the
region (Rocha et al. in prep.).

##  2 Observations

IRAS 23385+6053 was observed within the European MIRI GTO program JOYS
(proposal id. 1290) for 2.12 h on August 22, 2022. The entire MIRI bandpass
from $\sim$4.9 to $\sim$28.0 $\mu$m was covered in a 4-field mosaic (Fig. 1)
around the central position of R.A. (J2000.0) 23:40:54.497 and Dec. (J2000.0)
+61:10:27.83. The total exposure time in each the three MIRI spectral filters
was 200 sec, and parallel off-source imaging was conducted in the F1500W
filter. A complementary background field was observed for overall background
subtraction that is needed, e.g., an analysis of the broad-band ice features.
While the full-wavelength spectrum in Fig. 3 shows the non-background
subtracted data, all individual narrow gas spectral line data presented in
this paper were continuum subtracted by fitting a polynomial baseline to
narrow wavelength ranges around these lines.

Figure 4: Continuum images of IRAS 23385+6053. The color-scale presents the
continuum emission extracted from the MIRI IFU data cube in a line-free part
of each spectral setting centered on the wavelength labeled in each panel. The
contours outline the 1.3 mm emission (Beuther et al., 2018; Cesaroni et al.,
2019) starting at the 5$\sigma$ level of 0.55 mJy beam-1 and continuing in
25$\sigma$ steps. In the top-right of each panel, the corresponding spatial
resolution is shown (grey: mid-infrared, line: 1.3 mm). The top-left panel
labels the two main mid-infrared sources A and B as well as marks the main mm
peak mmA1 from Cesaroni et al. (2019) with a green star. A linear scale-bar is
shown as well. Note that mid-infrared source A and B are well separated up to
8 $\mu$m, but not at longer wavelengths.

We reduced the data using a development version of the JWST calibration
pipeline (Bushouse et al., 2022) with v. 11.16.16 and CRDS (Calibration
Reference Data System) context ’jwst_1017.pmap’, respectively. We used this
development version as it includes cosmic ray shower flagging in the jump step
of the Detector1Pipeline.

The pipeline includes fringe correction via Spec2Pipeline using a fringe flat
field derived from spatially extended sources (Mueller et al. in prep). We
further reduced the fringing by applying the residual_fringe step (Kavanagh et
al. in prep.), which is included in the JWST calibration pipeline package but
not switched on by default.

As discussed in Pontoppidan et al. (2022), the accuracy of spacecraft pointing
information for MRS exposures can be affected by guide star catalog errors and
roll uncertainty. We corrected our MRS exposures following the same procedure
as in that paper, i.e., using identified Gaia-
DR3111https://www.cosmos.esa.int/web/gaia/dr3 sources in the MIRI simultaneous
imaging field, leading to a pointing adjustment of $1.6077^{\prime\prime}$ in
Right Ascension and $0.3485^{\prime\prime}$ in Declination. This correction is
still subject to uncertainties in the relative astrometry between the imager
and the MRS, which is expected to be accurate to $\leq 0.1^{\prime\prime}$
(Patapis et al. in prep.).

Finally, we created spectral cubes for each of the 12 MRS sub-bands (four MIRI
channels with each containing a ”short”, ”medium” and ”long” grating setting)
using the level 2 data products with corrected astrometry from which spectra
were extracted. We applied the post-pipeline residual fringe correction tool
to all spectra to account for high frequency fringe residuals, thought to
originate in the MRS dichroics, which are not effectively removed by the
residual_fringe step in the pipeline (Kavanagh et al. in prep). The $1\sigma$
rms of the spectrum presented in Figure 3 is typically around 0.5 mJy below 20
$\mu$m and rises to a few mJy at longer wavelengths.

Figure 5: Outflow images in NOEMA SiO(2-1) (integrated emission, left panel),
JWST<EMAIL_ADDRESS>$\mu$m (middle panel) and
JWST<EMAIL_ADDRESS>$\mu$m (right panel). The red and blue contours show
the red- and blue-shifted NOEMA SiO(2–1) emission in the labeled velocity
regimes ($v_{\rm lsr}=-50.2$ km s-1). Contour levels are from 30 to 90% (step
10%) of the respective integrated peak emission. The magenta contours show the
5.2 $\mu$m continuum emission (steps 20 to 90% by 15%). The green contours in
the left panel outline the 1.3 mm continuum starting at the 20$\sigma$ levels
(1$\sigma$=1.13 mJy beam-1). A linear scale bar is presented in the middle
panel, and the green lines in the right panel outline and label potential
outflows. The two mid-infrared sources are labeled in the right panel as well.
The dashed line in the left panel marks a potentially precessing outflow
corresponding to outflow (I). A green star in all panels marks the main mm
peak position labeled mmA1 in Cesaroni et al. (2019). North and west are
labeled in the right panel. The resolution elements are shown at the bottom
right of each panel, left SiO(2–1) in grey and 1.3 mm continuum as line,
middle and right [FeII] and H2, respectively.

Complementary 1.3 mm continuum data are used from the CORE project (Beuther et
al., 2018; Cesaroni et al., 2019). The angular resolution and sensitivity of
these data are $0.48^{\prime\prime}\times 0.43^{\prime\prime}$ and 0.11 mJy
beam-1, respectively. Furthermore, SiO(2–1) observations are taken from the
CORE follow-up project CORE+ (PI: Caroline Gieser). This project aims to study
the shock and deuterium chemistry in the CORE sample. IRAS 23385+6053 is part
of the CORE+ pilot study (Gieser et al. in prep.). The corresponding angular
resolution and sensitivity of the SiO(2–1) data are $1.4^{\prime\prime}\times
1.2^{\prime\prime}$ and $\sim$0.12 K per 0.8 km s-1, respectively.

##  3 Results

###  3.1 General spectral features

Figure 3 presents the spectrum extracted towards the combined emission from
the two mid-infrared sources A and B (Sect. 3.2 and Fig. 4). We present the
combined spectrum from sources A and B because at wavelengths longward of
$\sim$13 $\mu$m, they cannot be spatially separated anymore (Fig. 4). The
color-coding corresponds to the four MIRI channels with the corresponding
three grating settings each, hence 12 individual sub-bands. It is evident that
JWST can now take high-quality mid-infrared spectra of sources as weak as 10
mJy in short integration times, more than a factor 1000 fainter than previous
data. Several features can directly be identified (see Table 1). Between
$\sim$5 and $\sim$11.5 $\mu$m, the spectrum is dominated by well-known
structures that stem from several ice and solid-state features (e.g., Gibb et
al. 2004; Yang et al. 2022; McClure et al. 2023). In particular, the silicate
band around 10 $\mu$m shows strong absorption, almost reaching full
saturation. Going to longer wavelengths, the general emission rises as
expected from warm dust emission. On top of the continuum emission, ice and
solid-state features as well as many strong emission lines can be identified,
in particular from H2 and atomic lines from [FeII], [ArII], [NiII], [NeII] and
[SI] (all labeled in Fig. 3). In addition to these strong lines, also a few
weaker molecular lines can be identified, in particular those of CO2, C2H2,
CH4 and HCN (Francis et al. in prep.). While the H2 and atomic line emission
is extended, the molecular emission mainly stems from the environment of the
continuum sources. Selected H2 and [FeII] line images are presented in Sect.
3.3, whereas the remaining gas and ice analysis will be discussed in
accompanying papers (Sect. 1).

Table 1: Main spectral features. Feature | $\lambda$ | Feature | $\lambda$
---|---|---|---
| ($\mu$m) |  | ($\mu$m)
H2 S(8) | 5.053 | PAH | 11.238
[FeII] | 5.340 | H2 S(2) | 12.279
H2 S(7) | 5.511 | Hi | 12.370
H2 S(6) | 6.109 | [NeII] | 12.814
[NiII] | 6.636 | CO2 ice/gas | 14.5–15.5
H2 S(5) | 6.910 | H2 S(1) | 17.035
[ArII] | 6.985 | [FeII] | 17.936
CH4 | 7.670 | [FeII] | 24.519
H2 S(4) | 8.025 | [SI] | 25.249
H2 S(3) | 9.665 | [FeII] | 25.988
Figure 6: Channel map of the H2(0–0) S(7) line at 5.511 $\mu$m. The velocity
of each channel is marked in units of km s-1 in the top-left corner of the
panels. The red and blue contours in the 33 km s-1 channel correspond to the
same red- and blue-shifted SiO(2–1) emission as presented in Fig. 5. A linear
scale-bar is shown in that panel as well. The $v_{\rm lsr}$ is $-50.2$ km s-1.

###  3.2 Continuum emission

Figure 4 presents the mid-infrared continuum emission extracted for each of
the twelve MIRI sub-bands averaged over wavelength ranges without significant
line emission. The central wavelengths of each image are marked in all panels.
We show the mid-infrared emission in comparison to the mm continuum emission
from the CORE project (Beuther et al. 2018, Cesaroni et al. 2019, see also
Fig. 1 for more detailed source labels). At the short wavelengths, one clearly
identifies two mid-infrared sources associated with the main mm continuum
core. That is also the area where the K-band emission is detected (Faustini et
al. 2009, Fig. 2). The north-western mid-infrared source A is closer to the
main mm peak position (mmA1 in Cesaroni et al. 2019) whereas the south-eastern
source B is offset from mid-infrared source A by $\sim$0.67′′. At the given
distance of 4.9 kpc, the angular offset between the two mid-infrared sources A
and B corresponds to a projected linear separation of $\sim$3280 au.

While the offset between the mid-infrared sources A and B is well determined,
the absolute positional offsets between the mm and mid-infrared data is less
well constrained. Positions from mm interferometer observations like the NOEMA
CORE data (Beuther et al., 2018) are largely constrained by the phase error,
and should be better than $0.1^{\prime\prime}$ for the presented mm data. The
corresponding 3 mm continuum data do also agree very well with the 1.3 mm data
presented here (Gieser et al. in prep.). In contrast to that, the mid-infrared
positional accuracy depends on the absolute JWST pointing, and the accuracy
with which the data can be improved by using known GAIA positions during the
data reduction. As outlined in Sect. 2, after correlating the parallel imaging
data with GAIA observations, the data were shifted in R.A./Dec. by
$\sim$1.6077${}^{\prime\prime}/0.3485^{\prime\prime}$. In addition to the
alignment with Gaia, this correction needs additional knowledge of the
alignment of the MIRI imager with the MIRI IFU. This alignment should be very
good and accurate to $0.1^{\prime\prime}$ (Sect. 2). This implies that mid-
infrared source A should indeed be very close to the elongation of the main mm
peak (the mm elongation is labeled peak mmA2 in Cesaroni et al. 2019) whereas
the offset of mid-infrared source B to the south-east should be real as well.

Going to longer wavelengths, in the spectral range of the silicate feature
between roughly 9 and 11 $\mu$m, source B almost disappears. This indicates
that the silicate absorption feature toward source B is even deeper and
entirely absorbs the continuum emission, whereas for source A some continuum
emission still remains. A detailed analysis of all solid state features will
be presented in Rocha et al. (in prep.). Moving to even longer wavelengths
beyond the silicate feature, source B re-appears. From $\sim$16 $\mu$m onward,
the further decreasing spatial resolution does not allow us to separate the
two mid-infrared sources anymore, and they merge into a single source at the
long-wavelength end of the MIRI bandpass.

###  3.3 H2, [FeII] and [NeII] line emission

The MIRI spectral line data now allow us to get a better understanding of the
underlying outflow structures in IRAS 23385+6053. While Molinari et al. (1998)
reported one outflow roughly in the north-south direction, Cesaroni et al.
(2019) recently suggested that a second outflow in the northwest-southeast
direction should exist (Fig. 1 right panel). The shocked H2 emission as well
as CH3OH maser emission at 44 and 95 GHz additionally confirm the existence of
shocked gas (Kurtz et al., 2004; Wolf-Chase et al., 2012). Here, we focus on
the high-spatial-resolution H2(0–0) S(7) line at 5.511 $\mu$m and the
[FeII](${}^{4}F_{9/2}$–${}^{6}D_{9/2}$) line at 5.34 $\mu$m. Furthermore, we
compare the results with SiO(2–1) data at 3.6 mm wavelengths obtained with
NOEMA as part of the CORE+ project (Gieser et al. in prep.). Figure 5 presents
a compilation of the three datasets. For the SiO(2–1) data, we show the
integrated emission ($\pm 15$ km s-1 around the $v_{\rm lsr}=-50.2$ km s-1) as
well as separated blue- and red-shifted high-velocity gas emission, whereas
for [FeII] and H2 we present the integrated emission.

Figure 7: Channel map of the Fe[II](${}^{4}F_{9/2}$–${}^{6}D_{9/2}$) line at
5.34 $\mu$m. The velocity of each channel is marked units of km s-1 in the
top-left corner of the panels. A linear scale-bar is shown in the top-left
panel. The $v_{\rm lsr}$ is $-50.2$ km s-1.

Starting with the mm SiO emission, these data are compatible with at least two
outflows. The integrated SiO(2–1) emission (grey-scale in Fig. 5 left) is
consistent with the northwest-southeast outflow ((I) in Fig. 5 right panel)
proposed by Cesaroni et al. (2019). The SiO emission also shows a more
extended integral-shaped structure in the northwest-southeast direction
(dashed line in Fig. 5 left, see discussion in Sect. 4.2). In comparison, the
red- and blue-shifted SiO emission rather points at an outflow more aligned
with the north-south direction ((IIa) & (IIb) in Fig. 5 right panel) as
suggested by Molinari et al. (1998). However, especially the blue-shifted SiO
emission shows extensions towards the east and west that may also be
associated with other outflow structures.

The [FeII] emission (Fig. 5 middle) exhibits a very different structure. It is
dominated by a collimated jet-like outflow in the northeast-southwest
direction ((III) in Fig. 5 right panel). But it also shows some extension
toward the west, spatially associated with the blue-shifted SiO emission.

Last but not least, the H2 data (Fig. 5 right) show emission structures
associated with all the before mentioned outflows. There is strong emission
associated with the northeast-southwest jet (III) traced also in [FeII], but
it shows additional emission towards the northwest, similar to the integrated
SiO emission (I). Furthermore, the H2 data also exhibit strong emission
towards the north (labeled (IIa) and (IIb) in Fig. 5 right panel), and
slightly weaker in the south, that corresponds to the red- and blue-shifted
SiO emission.

In addition to the integrated emission, with a spectral resolving power of
$\sim$3500 in channel 1A (Labiano et al. 2021222https://jwst-
docs.stsci.edu/jwst-mid-infrared-instrument/miri-observing-modes/miri-medium-
resolution-spectroscopy), we can resolve velocity elements down to $\sim$86 km
s-1 (corrected for radial and heliocentric velocities). With an approximate
Nyquist-sampling the channel separation is $\sim$44 km s-1. Figure 6 presents
a channel map for the H2 line at 5.511 $\mu$m. Interestingly, the gas
velocities are indeed so high that we can velocity-resolve the different
components. With the $v_{\rm{lsr}}$ of IRAS 23385+6053 of $-50.2$ km s-1
(Beuther et al., 2018), we are covering a velocity range of roughly $\pm 140$
km s-1. For comparison, we also present the channel map of the [FeII] line at
5.34 $\mu$m in Figure 7. Again, we can resolve blue- and red-shifted gas,
however with a slightly smaller velocity spread of around $\pm 100$ km s-1.

Investigating the velocity structures in a bit more detail, we find that the
H2 emission features towards the north are all redshifted. This is consistent
with the red-shifted SiO emission – although SiO is at lower velocities – that
is also found in the north (Fig. 5). In contrast to that, the strong H2
emission in the northeast is seen in all velocity channels, hence red- and
blue-shifted. If we assign that feature to the northeast-southwest jet (III)
that is also seen nicely in the [FeII] emission (Fig. 5 & 7), this implies
that this jet should be aligned close to the plane of the sky. The outflow (I)
in the northwest-southeast direction can also be identified in the H2 channel
map (Fig. 6), however, barely in blue-shifted gas and more easily at red-
shifted velocities. This velocity structure is not well recovered in the SiO
emission (Fig. 5).

Figure 8: Integrated [FeII](${}^{4}F_{9/2}$–${}^{6}D_{9/2}$) (left) and
[NeII](${}^{2}P_{1/2}-^{2}P_{3/2}$) (right) emission. The black contours show
the 5.2 $\mu$m continuum emission (steps 20 to 90% by 10%). The green contours
outline again the [NeII] emission in $10\sigma$ steps ($1\sigma\sim 0.09$ MJy
sr${}^{-1}\mu$m).

Another potential jet-tracer is the [NeII](${}^{2}P_{1/2}-^{2}P_{3/2}$) line
at 12.814 $\mu$m (e.g., Hollenbach & McKee 1989; Lefloch et al. 2003). Figure
8 presents a comparison of the [NeII] and the [FeII] emission, and one finds
two spatially very different distributions. While the [FeII] emission mainly
traces the northeast-southwest oriented collimated jet-like structure, the
[NeII] emission is rather double-peaked with one peak centered on mid-infrared
source A, whereas the second [NeII] emission peak is offset by $\sim$0.9′′ to
the west. Comparing the [NeII] emission directly to the [FeII], one finds weak
[FeII] emission associated also with the secondary [NeII] peak (left panel of
Fig. 8). This [FeII] emission arises mainly at relatively high projected red-
shifted velocities (see channels at $-$22 and +22 km s-1 in Fig. 7). We note
that no clear velocity structure can be spectrally resolved in the [NeII]
data. The [NeII] emission can stem from EUV and/or X-ray-irradiated surface
layers of accretion disks as well as from jets and outflows (e.g., Hollenbach
& McKee 1989; Lefloch et al. 2003; Glassgold et al. 2007; Pascucci & Sterzik
2009; Guedel et al. 2010). The J-shock models by Hollenbach & McKee (1989)
predict that the [NeII] emission should strongly increase for gas velocities
greater than $\sim$80 km s-1. In addition to this, high densities are also
increasing the [NeII] intensities in these shock models (Fig. 7 in Hollenbach
& McKee 1989). The spatial correspondence of the high-velocity red-shifted
[FeII] and [NeII] emission close to the center of the region, where gas
densities are also highest, is consistent with these J-shock models. In
addition to this, these regions close to the central protostars, traced by the
mid-infrared continuum emission, should also be exposed most to UV radiation.
Hence a combination of fast shocks and UV radiation may favor the [NeII]
emission in that region. The non-detection of compact [NeII] emission in the
rest of the IRAS 23385+6053 region, where partly also high gas velocities
and/or dense gas are observed (e.g., Figs. 6 and 1), confirms that not only
high velocities and dense gas are important, but that the missing UV radiation
further away from the main protostars is likely to play an important role as
well.

In summary, the combination of these datasets reveals at least three outflows,
one in the northwest-southeast direction (I), another in the northeast-
southwest direction (III) and at least one more in the north-south direction
(IIa & IIb). The fact that the H2 emission shows two emission peaks in the
north can be explained by either an expanding outflow cavity of one outflow or
alternatively by two more jet-like structures. Differentiating between these
two scenarios is not possible with the given data. The three, or potentially
even four outflows are marked with white lines in Fig. 5 (right panel).

###  3.4 Humphreys $\alpha$ line emission as an accretion tracer?

One crucial missing parameter in high-mass star formation is the actual
accretion rate onto the protostars. Here, we differentiate the accretion rate
onto the protostar from otherwise more typically reported gas infall rates
within the protostellar envelopes and cores (e.g, Myers et al. 1996; Mardones
et al. 1997; Beuther et al. 2013; Wyrowski et al. 2016; van Dishoeck et al.
2021). While for low-mass T Tauri stars the accretion rate can be measured
with optical or near-infrared lines (e.g., Connelley & Greene 2010; Salyk et
al. 2013; Hartmann et al. 2016; Alcalá et al. 2017), for the deeply embedded
phases in low- and high-mass star formation, one has to resort to longer
wavelengths. Here, we focus on the Humphreys $\alpha$ line Hi(7–6) at 12.37
$\mu$m. Rigliaco et al. (2015) explored this Hi(7–6) line at 12.37 $\mu$m in
Spitzer data of classical T Tauri stars, and they found a tentative
correlation between the integrated Hi(7–6) luminosity and the accretion
luminosity. We now explore this relation also for IRAS 23385+6053. Figure 9
presents the spectrum of the Hi(7–6) line extracted from the position
inbetween mid-infrared sources A and B, and averaged over an aperture with
radius of $2^{\prime\prime}$. While the red line shows a single Gaussian fit,
the blue line presents a two component Gaussian fit. It should be noted that
at this stage we consider this only a tentative detection, even though the
peak flux density of the single-Gaussian fit with a peak flux density of 0.61
mJy reaches a signal-to-noise ratio $S/N$ of 3 (in the channel with highest
flux densities even $S/N\sim 4.1$). The peak flux densities of the two-
component fit (peak flux densities of 0.78 and 0.65 mJy, respectively) have
nominal $S/N$ ratios of 3.8 and 3.2, respectively.

With a nominal central wavelength for the Hi(7–6) line of 12.3719 $\mu$m, the
two components are shifted by $-70$ and $+293$ km s-1, respectively. The full-
width-half-maximum line-widths of the two components are $\sim$73 and
$\sim$436 km s-1, respectively. While the two-component nature of the spectrum
may be attributed to gas associated with the innermost red- and blue-shifted
rotating gas, the line-width of several hundred km s-1 is consistent with gas
at almost free-fall onto the central protostar. The free-fall velocity of gas
around a 9 $M_{\odot}$ protostar (Cesaroni et al., 2019) with a radius of 5
$R_{\odot}$ would be $\sim$829 km s-1, but is likely narrower if the protostar
were significantly bloated (e.g., Hosokawa & Omukai 2009).

Figure 9: Hi(7–6) spectrum at 12.37 $\mu$m centered right in the middle
between sources A and B (Fig. 4), and averaged over an aperture with radius
$2^{\prime\prime}$. The red line presents a single Gaussian fit to the
spectrum whereas the blue line shows a two component Gaussian fit to the
spectrum.

Assuming this to be a real detection, we can derive the integrated flux values
of the one- and two-component fits. Since the integrated values of the
different fits are similar, the following analysis is done for the two-
component fit. However, with the single-component fit one gets comparable
results. The integrated line luminosity of the two-component Gaussian fit is
$\sim 2.3\times 10^{-4}$ $L_{\odot}$. This value lies at the high-luminosity
end of the relation fitted by Rigliaco et al. (2015) to a sample of low-mass
protostars. This then could be interpreted as an indicator that this relation
may also extend to the high-mass regime.

Following equation 1 in Rigliaco et al. (2015), one can convert the line
luminosity to an accretion luminosity of 142 $L_{\odot}$. Assuming now that
all that luminosity is from infalling gas that converts its energy to
radiation, one can infer an estimate of the accretion rate, e.g., Eq. 11.5 in
Stahler & Palla (2005). Taking these numbers at face-value, for an 9
$M_{\odot}$ protostar (Cesaroni et al., 2019) with a radius of 5 $R_{\odot}$
we can estimate an accretion rate of $\sim 2.6\times 10^{-6}$ $M_{\odot}$yr-1.

The above fitted integrated line luminosity is a lower limit since extinction
may weaken the real line emission. We now estimate an extinction correction
using the extinction curve derived by McClure (2009). Based on the 9.7 $\mu$m
silicate absorption feature, Rocha et al. (in prep.) estimate an extinction
$A_{V}\sim 30-40$ mag. Using $A_{V}\sim$30 mag, and following McClure (2009),
that extinction can be converted to a K-band $A_{K}$ and then 12.37 $\mu$m
extinction values of 3.8 and 1.86 mag, respectively. Such a 12.37 $\mu$m
extinction would boost the integrated Hi(7–6) line luminosity by a factor
$\sim$5.56. Using the same approach as above, one can estimate now an
extinction-corrected accretion luminosity of $\sim 5070$ $L_{\odot}$, and an
accretion rate of $\sim 0.9\times 10^{-4}$ $M_{\odot}$yr-1. The latter is
roughly a factor $\sim$36 higher than the estimate without the extinction
correction.

One caveat is that the extinction-corrected accretion luminosity of $\sim
5070$ $L_{\odot}$ is on the same order as the luminosity estimated for the
source of $\sim$3000 $L_{\odot}$ (Molinari et al., 2008). Because of that,
higher accretion luminosities would be unreasonable, and correspondingly, we
do not consider higher extinction values. This accretion luminosity would
imply that most of the bolometric luminosity originates from accretion, thus
the source could be in its main accretion phase. If the approach outlined by
Rigliaco et al. (2015) and followed here is also applicable for young high-
mass protostars, the extinction-corrected accretion rate estimate has to be an
upper limit to the real accretion rate onto the protostar. The accretion rates
will be discussed in depth in Sect. 4.3.

##  4 Discussion

###  4.1 Embedded protostars

What is the nature of the two mid-infrared sources A and B? Are they
individual protostars or are they potentially part of an extended disk-like
structure perpendicular to the northeast-southwest outflow best depicted in
the [FeII] data (Fig. 5)?

The spatial separation between sources A and B of $0.67^{\prime\prime}$
corresponds at the measured distance of 4.9 kpc to a linear separation of
$\sim$3300 au. This seems too large to be consistent with typical disk sizes,
even around high-mass stars (e.g., Kuiper et al. 2010).

Furthermore, we can inspect the full MIRI spectra towards the continuum
sources. Fig. 3 is extracted within a diameter of $2.4^{\prime\prime}$
combining sources A and B. At the shorter wavelengths, where the sources can
still be spatially separated, the main difference between A and B are an even
deeper silicate absorption feature towards source B which is also manifested
in the non-detection of source B in the continuum images at 9.0 and 10.65
$\mu$m (Fig. 4). Comparing the MIRI spectrum to typical spectra from embedded
protostars to evolved disks (e.g., Evans et al. 2003; Gibb et al. 2004), the
MIRI spectrum presented here shows all the features typical for embedded
protostars, from deep silicate absorption features around 10 $\mu$m to the
rising spectral energy distribution towards long wavelengths to several ice
features at shorter wavelengths. Hence, the spectra support the assessment
that sources A and B are two separate embedded protostars.

As outlined in Cesaroni et al. (2019), the IRAS 23385+6053 star-formation
complex is a high-mass cluster in formation. They identify already six
potential star-forming cores from the 1.3 mm continuum emission (Fig. 5 in
Cesaroni et al. 2019, see also the right panel of our Fig. 1). The main mm
peak position – source mmA1 in Cesaroni et al. (2019), see Figs. 1, 4 and 5 –
has no obvious mid-infrared counterpart. However, the main mm peak is
elongated in the southwestern direction (mm source mmA2 in Cesaroni et al.
2019), and this mm-peak mmA2 is likely the counterpart to the mid-infrared
source A identified with JWST. Ahmadi et al. (subm.) analyze the CH3CN mm line
emission around the mm sources mmA1 and mmA2, and in their disk stability
analysis they find very low Toomre Q values near mm peak mmA2 corresponding to
the here identified mid-infrared source A. Ahmadi et al. (subm.) argue that
the low Toomre value could indicate a fragmenting disk-like structure. In that
picture, mid-infrared source A could have potentially formed in that
fragmented disk.

In comparison to that, the mid-infrared source B has no clear mm peak
counterpart and hence can be considered as an independent entity, potentially
of lower mass. Both the mm as well as the mid-infrared observations are
sensitivity limited, hence more sources should likely exist below our
detection limits that could be revealed by deep MIRI imaging.

Although the MIRI spectra are dominated by protostellar envelope features, the
fact that we identify several outflows, and that two of them may likely
emanate from the mid-infrared sources A and B (see the following Sect. 4.2),
the existence of embedded accretion disks is very likely. They are just
difficult to identify because of the limited spatial resolution given the
distance of the source (4.9 kpc), and the spectra being dominated by the
surrounding envelopes. In addition to that, rotational signatures indicative
of an embedded disk around the main mm continuum source mmA1 was already
derived from mm CH3CN data (Cesaroni et al., 2019).

###  4.2 Outflows from mid-infrared to mm wavelengths

The velocity ranges observed by the mid-infrared H2 and [FeII] lines with more
than $\pm 100$ km s-1 are much broader than those observed at mm wavelengths
in the SiO emission of around $\pm 15$ km s-1 (Fig. 5). Furthermore, the
[NeII] emission requires high gas velocities as well (e.g., Hollenbach & McKee
1989). Therefore, the H2, [FeII] and [NeII] lines are tracing genuinely higher
velocity gas than the lines typically observed at mm wavelengths. This implies
that the mid-infrared lines originate directly from the underlying jet,
whereas the SiO emission rather traces outflowing gas at lower velocities.
Since silicon (Si) has to be liberated from grains, some (slow) shocks are
also needed for SiO emission (e.g., Schilke et al. 1997; Anderl et al. 2013).
Hence the SiO emission likely stems from a combination of low-velocity shocks
caused by the jet, as well as entrained outflow gas. This behavior is also
reflected by the excitation energies of the different lines. While the
SiO(2–1) line has an upper level excitation energy $E_{u}/k$ of only 6.3 K
(e.g., Schöier et al. 2005), the presented H2, [FeII] and [NeII] lines at
5.511, 5.34 and 12.814 $\mu$m have much higher excitation energies $E_{u}/k$
$\sim$7202 K, 2694 K and 1123 K,
respectively333https://www.mpe.mpg.de/ir/ISO/linelists/Molecular.html,
https://www.pa.uky.edu/$\sim$peter/newpage/. Hence, the mid-infrared lines
trace ionized as well as neutral shocked material, either from the jet itself
or also from the outflow cavity walls. In contrast, the mm lines, even SiO
that is typically also assumed to be a shock tracer (e.g., Schilke et al.
1997), traces colder components of the jet, where molecules can form at the
internal working surfaces (e.g., Santiago-García et al. 2009). The lesson is
that studying the jets and outflows with a single tracer may miss some
components. In contrast to that, the combined mid-infrared and mm data reveal
a much more complete picture of the dynamical outflows driven from various
embedded protostars. Furthermore, this outflow multiplicity also enforces the
picture of high-mass stars forming exclusively in a clustered mode.

The more extended integral-shaped SiO emission in Fig. 5 corresponds to the
northwest-southeast outflow (I). Could that be the signature of a precessing
outflow or jet? Different jet-precession scenarios are discussed in, e.g.,
Fendt & Zinnecker (1998), and a potential precession reason is the outflow-
driving within a binary or multiple system. Since we are dealing with a
cluster-forming region and several identified protostars (Sect. 4.1),
precession induced by multiplicity may indeed explain the observed bent SiO
emission.

With multiple outflows, it is also important to assess the potential driving
sources for each of them. Inspecting Figure 5 again, the [FeII] jet (III) is
most closely centered on the mid-infrared source A, but neither directly
spatially aligned with mid-infrared source B or mm peak mmA1. Therefore, the
most likely driving source for outflow (III) appears to be the mid-infrared
source A (spatially coincident with the secondary mm-peak mmA2).

Regarding the northwest-southeastern outflow (I), Cesaroni et al. (2019)
proposed a rather straight orientation center on a mm source C (identified in
the H2CO emission, Fig. 1 right panel) at the southwestern edge of the mm
continuum emission (in the middle between the red and blue outflow arrows in
Fig. 1, right panel). Considering the more bent structure we infer from the
SiO and H2 emission (Fig. 5), the driving source is more likely close to the
mid-infrared sources A and B. Since source A is the likely driver of outflow
(III), we propose mid-infrared source B as a candidate driving source for
outflow (I).

Finally, the one (or potentially two) north-south outflows (IIa & IIb) are
spatially associated with mid-infrared sources A and B as well as the main mm-
peak mmA1. Since mid-infrared sources A and B are more likely linked to
outflows (I) and (III), we propose the mm-peak mmA1 to host the driving
source(s) of these north-south outflow(s).

###  4.3 Accretion rates

How do the Humphreys $\alpha$ Hi(7–6) inferred accretion rate estimates
between $2.6\times 10^{-6}$ and $0.9\times 10^{-4}$ $M_{\odot}$yr-1 (Sect.
3.4) compare to outflow rates? Molinari et al. (1998) report integrated
outflow rates for IRAS 23385+6053 of $\geq 10^{-3}$ $M_{\odot}$yr-1 based on
HCO+ data. Based on their lower spatial resolution, these outflow rates likely
combine the different outflows we see in the new data. Assuming momentum
conservation between initial jet and entrained outflow rate, and a velocity
ratio between jet and entrained gas of $\sim$20 (e.g., Beuther et al. 2002),
the jet-mass-flow rate is roughly an order of magnitude lower. Furthermore,
jet models result in ratios between jet-flow rate and accretion rate of
roughly 0.3 (e.g., Tomisaka 1998; Shu et al. 1999). Hence, the measured
outflow rates should correspond to accretion rates on the order 10-4
$M_{\odot}$yr-1. While such accretion rate estimates are almost two orders of
magnitude higher than the non-extinction-corrected Hi(7–6) accretion rate
estimate, they are in the same ballpark as those derived when considering the
extinction correction.

Cesaroni et al. (2019) also estimated accretion rates based on the assumption
that all cores are collapsing in free-fall, and they get accretion rates
between a few times $10^{-4}$ and a few times $10^{-3}$ $M_{\odot}$yr-1. The
corresponding time-scales they estimate from the accretion rates are mostly
below $10^{4}$ yrs. While these accretion rate estimates are again at the
upper end of our Hi(7–6) based estimates, Cesaroni et al. (2019) point out
that their estimated time-scales are too short for typical high-mass star
formation times on the order $10^{5}$ yr. Hence, their estimated time-scales
are lower limits, and correspondingly their estimated accretion rates upper
limits.

Our Hi(7–6) estimated accretion rates appear at the lower end of high-mass
star formation accretion rates if one considers that accretion rates on the
order $10^{-4}-10^{-3}$ $M_{\odot}$yr-1 are needed to form a high-mass star
within a few hundred thousand years (e.g., McKee & Tan 2003). While we cannot
entirely exclude that the main object is not a 9 $M_{\odot}$ protostar but
maybe a multiple system containing lower-mass objects, this seems unlikely
since the luminosity of $3\times 10^{3}$ $L_{\odot}$ (Molinari et al., 2008)
requires a more massive central object (e.g., Mottram et al. 2011; Cesaroni et
al. 2019). The luminosity is also consistent with the estimated $\sim$9
$M_{\odot}$ embedded object based on the rotation curves measured in CH3CN
(Cesaroni et al., 2019). In comparison to that, in the extinction-corrected
accretion estimate, almost the entire luminosity should stem from the
accretion processes.

Furthermore, with a total mass reservoir of the IRAS 23385+6053 star-forming
region of $\sim$510 $M_{\odot}$, and assuming a typical initial mass function
and star formation efficiency, no star above 10 $M_{\odot}$ would be expected
in this region (e.g., Tackenberg et al. 2012). One can estimate, for example,
the maximum stellar mass $m^{*}_{\rm{max}}$ of a cluster following Larson
(2003), Sanhueza et al. (2019) and Morii et al. (subm.):

$m^{*}_{\rm{max}}=15.6\left(\frac{M_{\rm{clump}}}{10^{3}\rm{M}_{\odot}}\frac{\epsilon_{\rm{SFE}}}{0.3}\right)^{0.45}M_{\odot}.$

With a clump mass $M_{\rm{clump}}=510$ $M_{\odot}$ and a star formation
efficiency $\epsilon_{\rm{SFE}}\sim 0.15$, the resulting maximum stellar mass
$m^{*}_{\rm{max}}$ would be $\sim$8.45 $M_{\odot}$, roughly what is found for
IRAS23385+6053 (Cesaroni et al., 2019). Therefore, it could well be that the
main accretion phase is coming to an end and hence we see comparably lower
accretion rates in the Hi(7–6) line. A similar decrease of accretion rates
with time is also seen in low-mass regions (e.g., Evans et al. 2009).

While the outflow rate, and by that the inferred accretion rate, is an
integral measure over the outflow time-scale, the Hi(7–6) inferred accretion
rate should rather correspond to an instantaneous accretion rate onto the star
at the time of observations. Considering that a 9 $M_{\odot}$ protostar has
already formed at the center of IRAS 23385+6053, it may well be that the
outflow-inferred accretion rates correspond to the active accretion phase that
comes to an end now. In that picture, the Hi(7–6) accretion rate gives an
estimate of the still ongoing remaining accretion. In addition to this, (high-
mass) star formation is known to be episodic (e.g., Caratti o Garatti et al.
2017; Kuffmeier et al. 2017; Hunter et al. 2021; Elbakyan et al. 2021), and in
that framework IRAS 23395+6053 may also be in a low-accretion phase at the
time of observations.

Some potential caveats need to be discussed. The original relation between
accretion luminosity and Hi(7–6) line luminosity was inferred for low-mass
protostars, and it is not obvious that this relation has to hold also for the
high-mass regime. The Hi(7–6) line is a recombination line, and while that
should arise in the low-mass regime mainly from accretion processes, high-mass
protostars can also emit UV radiation that could result in recombination line
emission. However, towards IRAS 23385+6053, only extended cm continuum
emission, most likely associated with the extended mid-infrared nebula (Fig.
1), has been detected (Molinari et al., 2002), but no compact cm free-free
emission towards the central protostar could be identified (Molinari et al.,
1998). Therefore, no strong ionizing radiation is emitted from the embedded
high-mass protostar. Furthermore, in that evolutionary stage, high-mass
protostars can also be bloated and then have much lower surface temperatures
(e.g., Hosokawa & Omukai 2009; Hosokawa et al. 2010), again reducing the
potential ionizing radiation. Therefore, from a conceptional point of view, in
the early evolutionary stages of high-mass star formation where the protostar
can be bloated and does not yet emit much ionizing radiation, the relation
between accretion luminosity and Hi(7–6) line luminosity may still be valid.

An additional caveat is that the Hi(7–6) detection is only at a $\sim
3-4\sigma$ level. Other HI emission lines (e.g., Hi(6–5) or Hi(8–7)) are also
in the observed bandpass but remain undetected in IRAS23395+6053. Therefore,
more observations of other high-mass star-forming regions are required to
explore the capability of tracing accretion with the Hi(7–6) and other
recombination lines in more depth. This will become possible, given the fact
that most JOYS targets are still pending in the JWST queue. With a sample of
approximately two dozen regions, sufficient data will become available to
derive statistically unambiguous values. The present work is a first step in
realizing this.

##  5 Conclusions and summary

We present among the first JWST MIRI MRS mid-infrared observations of a young
high-mass star-forming region, IRAS 23385+6053. The spectral coverage between
$\sim$5 and $\sim$28 $\mu$m reveals a plethora of spectral gas lines from
atomic and molecular species. Furthermore, several broader ice-features can be
identified that will be studied in accompanying papers.

Investigating the mid-infrared continuum emission, two mid-infrared sources
are identified at the shortest wavelength with a separation of $\sim$0.67′′
(or $\sim$3280 au). At the long-wavelength end of the spectrum, these two
sources merge into one emission structure because of the lower angular
resolution. While one of the mid-infrared sources (A) is associated with a mm
continuum source (mmA2 in Cesaroni et al. 2019), the other mid-infrared source
B has no obvious mm counterpart. The MIRI spectra of the two sources clearly
show their deeply embedded protostellar nature. Combining the JWST mid-
infrared and previously obtained mm data confirms that we are indeed observing
a cluster in its making.

For the outflow analysis, we focus on two shorter wavelength lines to take
advantage of the high spatial resolution. Combining the outflow-tracing JWST
MIRI H2 at 5.511 $\mu$m and [FeII] at 5.34 $\mu$m lines with the mm SiO(2–1)
emission, we identify at least three independent outflows in the region, and
it is possible to assign individual candidate driving sources to each of them.
Furthermore, although only observed at a spectral resolving power of
$\sim$3500, corresponding to a velocity-resolution of $\sim$86 km s-1, with
given velocity shifts of $\pm 140$ km s-1, we can spatially and spectrally
resolve the red- and blue-shifted high-velocity gas. In particular the red-
shifted high-velocity gas is spatially associated with red-shifted SiO
emission, although the latter is at lower velocities of $\sim$15 km s-1.
Furthermore, the potential additional mid-infrared outflow tracer [NeII] shows
a partially overlapping morphology, with two emission peaks close to the
center of the region. The two [NeII] emission structures can be associated
with high velocity gas also seen in the [FeII] emission from the northwest-
southeast outflow (I), but do not fully reproduce all observed [FeII]
structures. These data show that the higher excited mid-infrared data trace
also the higher-velocity jet-like structures whereas the lower-excited mm
lines preferentially trace lower-velocity jet and entrained outflow gas of the
broader molecular outflow.

Our investigation of a potential accretion rate tracer, the Humphreys $\alpha$
Hi(7–6) line, revealed a 3–4$\sigma$ detection. Using the relation between
integrated Hi(7–6) line luminosity and accretion luminosity, and depending on
extinction-correction, we can estimate accretion luminosities between 142 and
5070 $L_{\odot}$. Assuming furthermore a reasonable protostellar mass and
radius, these luminosities convert to accretion rate estimates between
$\sim$2.6$\times 10^{-6}$ and $\sim$0.9$\times 10^{-4}$ $M_{\odot}$yr-1.
Setting this into context with an accretion rate estimates based on outflow
observations, IRAS 23385+6053 may already be close to the end of the main
accretion phase. However, with data for one region so far and then only a
3–4$\sigma$ detection, these results have to be taken with a bit of caution.
Observations of more high-mass star-forming regions are needed to further
evaluate the possibility of estimating accretion rates with the Hi(7–6) line.
Nevertheless, the present work shows that the Humphreys $\alpha$ line has the
potential to become an important accretion rate tracer with JWST. Future MIRI
observations, also within the JOYS program, will set much tighter constraints
on that.

In summary, these first JWST MIRI MRS observations of a high-mass star-forming
region reveal important new results of the protostellar distribution, the
outflow properties as well as accretion rates in IRAS 23385+6053. The
complementary atomic and molecular gas lines as well as the ice features are
discussed in companion papers. While interesting in themselves, these data
also show the enormous potential for JWST to boost star formation research.

###### Acknowledgements.

The following National and International Funding Agencies funded and supported
the MIRI development: NASA; ESA; Belgian Science Policy Office (BELSPO);
Centre Nationale d’Etudes Spatiales (CNES); Danish National Space Centre;
Deutsches Zentrum fur Luftund Raumfahrt (DLR); Enterprise Ireland; Ministerio
De Economiá y Competividad; Netherlands Research School for Astronomy (NOVA);
Netherlands Organisation for Scientific Research (NWO); Science and Technology
Facilities Council; Swiss Space Office; Swedish National Space Agency; and UK
Space Agency. We thank Fabiana Faustini for providing the infrared data. We
also thank Roy van Boekel for mid-infrared line diagnostic discussions. H.B.
acknowledges support from the Deutsche Forschungsgemeinschaft in the
Collaborative Research Center (SFB 881) “The Milky Way System” (subproject
B1). EvD, MvG, LF, KS, WR and HL acknowledge support from ERC Advanced grant
101019751 MOLDISK, TOP-1 grant 614.001.751 from the Dutch Research Council
(NWO), the Netherlands Research School for Astronomy (NOVA), the Danish
National Research Foundation through the Center of Excellence “InterCat”
(DNRF150), and DFG-grant 325594231, FOR 2634/2. P.J.K. acknowledges financial
support from the Science Foundation Ireland/Irish Research Council Pathway
programme under Grant Number 21/PATH-S/9360. A.C.G. has been supported by
PRIN-INAF MAIN-STREAM 2017 “Protoplanetary disks seen through the eyes of new-
genera-tion instruments” and from PRIN-INAF 2019 “Spectroscopically tracing
the disk dispersal evolution (STRADE)”. K.J. acknowledges the support from the
Swedish National Space Agency (SNSA). T.H. acknowledges support from the
European Research Council under the Horizon 2020 Framework Program via the ERC
Advanced Grant ”Origins” 83 24 28.

## References

  * Alcalá et al. (2017) Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20
  * An et al. (2011) An, D., Ramírez, S. V., Sellgren, K., et al. 2011, ApJ, 736, 133
  * Anderl et al. (2013) Anderl, S., Guillet, V., Pineau des Forêts, G., & Flower, D. R. 2013, A&A, 556, A69
  * Beuther et al. (2013) Beuther, H., Linz, H., & Henning, T. 2013, A&A, 558, A81
  * Beuther et al. (2018) Beuther, H., Mottram, J. C., Ahmadi, A., et al. 2018, A&A, 617, A100
  * Beuther et al. (2002) Beuther, H., Schilke, P., Sridharan, T. K., et al. 2002, A&A, 383, 892
  * Bushouse et al. (2022) Bushouse, H., Eisenhamer, J., Dencheva, N., et al. 2022, JWST Calibration Pipeline
  * Caratti o Garatti et al. (2017) Caratti o Garatti, A., Stecklum, B., Garcia Lopez, R., et al. 2017, Nature Physics, 13, 276
  * Casoli et al. (1986) Casoli, F., Dupraz, C., Gerin, M., Combes, F., & Boulanger, F. 1986, A&A, 169, 281
  * Cesaroni et al. (2019) Cesaroni, R., Beuther, H., Ahmadi, A., et al. 2019, A&A, 627, A68
  * Cesaroni et al. (1988) Cesaroni, R., Palagi, F., Felli, M., et al. 1988, A&AS, 76, 445
  * Connelley & Greene (2010) Connelley, M. S. & Greene, T. P. 2010, AJ, 140, 1214
  * Elbakyan et al. (2021) Elbakyan, V. G., Nayakshin, S., Vorobyov, E. I., Caratti o Garatti, A., & Eislöffel, J. 2021, A&A, 651, L3
  * Evans et al. (2003) Evans, Neal J., I., Allen, L. E., Blake, G. A., et al. 2003, PASP, 115, 965
  * Evans et al. (2009) Evans, N. J., Dunham, M. M., Jørgensen, J. K., et al. 2009, ApJS, 181, 321
  * Faustini et al. (2009) Faustini, F., Molinari, S., Testi, L., & Brand, J. 2009, A&A, 503, 801
  * Fendt & Zinnecker (1998) Fendt, C. & Zinnecker, H. 1998, A&A, 334, 750
  * Gibb et al. (2004) Gibb, E. L., Whittet, D. C. B., Boogert, A. C. A., & Tielens, A. G. G. M. 2004, ApJS, 151, 35
  * Gibb et al. (2000) Gibb, E. L., Whittet, D. C. B., Schutte, W. A., et al. 2000, ApJ, 536, 347
  * Gieser et al. (2021) Gieser, C., Beuther, H., Semenov, D., et al. 2021, A&A, 648, A66
  * Glassgold et al. (2007) Glassgold, A. E., Najita, J. R., & Igea, J. 2007, ApJ, 656, 515
  * Guedel et al. (2010) Guedel, M., Lahuis, F., Briggs, K. R., et al. 2010, ArXiv e-prints [arXiv:1006.2848]
  * Hartmann et al. (2016) Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135
  * Hollenbach & McKee (1989) Hollenbach, D. & McKee, C. F. 1989, ApJ, 342, 306
  * Hosokawa & Omukai (2009) Hosokawa, T. & Omukai, K. 2009, ApJ, 691, 823
  * Hosokawa et al. (2010) Hosokawa, T., Yorke, H. W., & Omukai, K. 2010, ApJ, 721, 478
  * Hunter et al. (2021) Hunter, T. R., Brogan, C. L., De Buizer, J. M., et al. 2021, ApJ, 912, L17
  * Kuffmeier et al. (2017) Kuffmeier, M., Haugbølle, T., & Nordlund, Å. 2017, ApJ, 846, 7
  * Kuiper et al. (2010) Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2010, ApJ, 722, 1556
  * Kurtz et al. (2004) Kurtz, S., Hofner, P., & Álvarez, C. V. 2004, ApJS, 155, 149
  * Labiano et al. (2021) Labiano, A., Argyriou, I., Álvarez-Márquez, J., et al. 2021, A&A, 656, A57
  * Larson (2003) Larson, R. B. 2003, Reports on Progress in Physics, 66, 1651
  * Lefloch et al. (2003) Lefloch, B., Cernicharo, J., Cabrit, S., et al. 2003, ApJ, 590, L41
  * Mardones et al. (1997) Mardones, D., Myers, P. C., Tafalla, M., et al. 1997, ApJ, 489, 719
  * McClure (2009) McClure, M. 2009, ApJ, 693, L81
  * McClure et al. (2023) McClure, M. K., Rocha, W. R. M., Pontoppidan, K. M., et al. 2023, arXiv e-prints, arXiv:2301.09140
  * McKee & Tan (2003) McKee, C. F. & Tan, J. C. 2003, ApJ, 585, 850
  * Molinari et al. (2008) Molinari, S., Faustini, F., Testi, L., et al. 2008, A&A, 487, 1119
  * Molinari et al. (1998) Molinari, S., Testi, L., Brand, J., Cesaroni, R., & Palla, F. 1998, ApJ, 505, L39
  * Molinari et al. (2002) Molinari, S., Testi, L., Rodríguez, L. F., & Zhang, Q. 2002, ApJ, 570, 758
  * Motte et al. (2018) Motte, F., Bontemps, S., & Louvet, F. 2018, ARA&A, 56, 41
  * Mottram et al. (2011) Mottram, J. C., Hoare, M. G., Urquhart, J. S., et al. 2011, A&A, 525, A149
  * Myers et al. (1996) Myers, P. C., Mardones, D., Tafalla, M., Williams, J. P., & Wilner, D. J. 1996, ApJ, 465, L133
  * Palla et al. (1991) Palla, F., Brand, J., Cesaroni, R., Comoretto, G., & Felli, M. 1991, A&A, 246, 249
  * Pascucci & Sterzik (2009) Pascucci, I. & Sterzik, M. 2009, ApJ, 702, 724
  * Pontoppidan et al. (2022) Pontoppidan, K. M., Barrientes, J., Blome, C., et al. 2022, ApJ, 936, L14
  * Rigliaco et al. (2015) Rigliaco, E., Pascucci, I., Duchene, G., et al. 2015, ApJ, 801, 31
  * Salyk et al. (2013) Salyk, C., Herczeg, G. J., Brown, J. M., et al. 2013, ApJ, 769, 21
  * Sanhueza et al. (2019) Sanhueza, P., Contreras, Y., Wu, B., et al. 2019, ApJ, 886, 102
  * Santiago-García et al. (2009) Santiago-García, J., Tafalla, M., Johnstone, D., & Bachiller, R. 2009, A&A, 495, 169
  * Schilke et al. (1997) Schilke, P., Walmsley, C. M., Pineau des Forets, G., & Flower, D. R. 1997, A&A, 321, 293
  * Schöier et al. (2005) Schöier, F. L., van der Tak, F. F. S., van Dishoeck, E. F., & Black, J. H. 2005, A&A, 432, 369
  * Shu et al. (1999) Shu, F. H., Allen, A., Shang, H., Ostriker, E. C., & Li, Z. 1999, in NATO ASIC Proc. 540: The Origin of Stars and Planetary Systems, 193
  * Stahler & Palla (2005) Stahler, S. W. & Palla, F. 2005, The Formation of Stars (ISBN 3-527-40559-3. Wiley-VCH)
  * Tackenberg et al. (2012) Tackenberg, J., Beuther, H., Henning, T., et al. 2012, A&A, 540, A113
  * Tomisaka (1998) Tomisaka, K. 1998, ApJ, 502, L163
  * van Dishoeck et al. (2021) van Dishoeck, E. F., Kristensen, L. E., Mottram, J. C., et al. 2021, A&A, 648, A24
  * van Dishoeck et al. (1998) van Dishoeck, E. F., Wright, C. M., Cernicharo, J., et al. 1998, ApJ, 502, L173
  * Whittet et al. (1996) Whittet, D. C. B., Schutte, W. A., Tielens, A. G. G. M., et al. 1996, A&A, 315, L357
  * Wolf-Chase et al. (2012) Wolf-Chase, G., Smutko, M., Sherman, R., Harper, D. A., & Medford, M. 2012, ApJ, 745, 116
  * Wouterloot & Brand (1989) Wouterloot, J. G. A. & Brand, J. 1989, A&AS, 80, 149
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * Wyrowski et al. (2016) Wyrowski, F., Güsten, R., Menten, K. M., et al. 2016, A&A, 585, A149
  * Yang et al. (2022) Yang, Y.-L., Green, J. D., Pontoppidan, K. M., et al. 2022, ApJ, 941, L13


11institutetext: Institute for Numerical Mathematics, Johannes-Kepler
University Linz, Austria
11email<EMAIL_ADDRESS>11email<EMAIL_ADDRESS>22institutetext: Johann Radon Institute for Computational and Applied
Mathematics, Linz, Austria
22email<EMAIL_ADDRESS>
# A reduced scalar potential approach for magnetostatics avoiding the coenergy

Herbert Egger 1122 Felix Engertsberger 11 Bogdan Radu 22

###### Abstract

The numerical solution of problems in nonlinear magnetostatics is typically
based on a variational formulation in terms of magnetic potentials, the
discretization by finite elements, and iterative solvers like the Newton
method. The vector potential approach aims at minimizing a certain energy
functional and, in three dimensions, requires the use of edge elements and
appropriate gauging conditions. The scalar potential approach, on the other
hand, seeks to maximize the negative coenergy and can be realized by standard
Lagrange finite elements, thus reducing the number of degrees of freedom and
simplifying the implementation. The number of Newton iterations required to
solve the governing nonlinear system, however, has been observed to be usually
higher than for the vector potential formulation. In this paper, we propose a
method that combines the advantages of both approaches, i.e., it requires as
few Newton iterations as the vector potential formulation while involving the
magnetic scalar potential as the primary unknown. We discuss the variational
background of the method, its well-posedness, and its efficient
implementation. Numerical examples are presented for illustration of the
accuracy and the gain in efficiency compared to other approaches.

##  1 Introduction and motivation

The equilibrium distributions of the magnetic fields generated by an imposed
current density $\mathbf{j}_{s}=\operatorname{curl}\mathbf{h}_{s}$ in a
nonlinear magnetically isolated medium satisfy

$\displaystyle\partial_{\mathbf{b}}w(\mathbf{b})+\nabla\psi$
$\displaystyle=\mathbf{h}_{s}\qquad$ $\displaystyle\text{in }\Omega,$ (1)
$\displaystyle\operatorname{div}\mathbf{b}$ $\displaystyle=0\qquad$
$\displaystyle\text{in }\Omega$ (2) $\displaystyle\mathbf{b}\cdot\mathbf{n}$
$\displaystyle=0\qquad$ $\displaystyle\text{on }\partial\Omega.$ (3)

Here $\mathbf{b}$ is the magnetic flux density, $\mathbf{h}_{s}$ the source
field, and $\psi$ the magnetic scalar potential. The magnetic field intensity
$\mathbf{h}=\partial_{b}w(\mathbf{b})$ is obtained as the derivative of the
magnetic energy density $w(\mathbf{b})=w(\mathbf{b},x)$, which describes the
local nonlinear material response [16] and may in general depend on the
spatial position. The first equation thus amounts to the decomposition
$\mathbf{h}=\mathbf{h}_{s}-\nabla\psi$ and Ampere’s law
$\operatorname{curl}\mathbf{h}=\mathbf{j}_{s}$ can be recovered by application
of the $\operatorname{curl}$ operator.

### Vector potential formulation

The magnetic Gauß law (2) can be satisfied by representing
$\mathbf{b}=\operatorname{curl}\mathbf{a}$ in terms of a magnetic vector
potential. The remaining equations can then be restated as

$\displaystyle\operatorname{curl}(\partial_{\mathbf{b}}w(\operatorname{curl}\mathbf{a}))$
$\displaystyle=\mathbf{j}_{s}\qquad$ $\displaystyle\text{in }\Omega$ (4)
$\displaystyle\mathbf{a}\times\mathbf{n}$ $\displaystyle=0\qquad$
$\displaystyle\text{on }\partial\Omega.$ (5)

We tacitly assumed that the boundary $\partial\Omega$ is simply connected and
further applied the $\operatorname{curl}$ operator to (1) as before. Hence (4)
amounts to Ampere’s law stated in terms of the magnetic vector potential. The
weak formulation of this problem is well-suited for analysis and
discretization by $H(\operatorname{curl})$-conforming finite elements [12,
15], and the resulting nonlinear systems can be solved efficiently by Newton-
type methods [7]. Appropriate gauging is required to ensure the uniqueness of
the vector potential on the continuous and the discrete level [1, 10].

### Scalar potential formulation

Let
$w_{*}(\mathbf{h})=\sup_{\mathbf{b}}\\{\langle\mathbf{h},\mathbf{b}\rangle-w(\mathbf{b})\\}$
denote the convex conjugate coenergy density. This allows to recast the
material law $\mathbf{h}=\partial_{\mathbf{b}}w(\mathbf{b})$ in inverse form
$\mathbf{b}=\partial_{\mathbf{h}}w_{*}(\mathbf{h})$. Using
$\mathbf{h}=\mathbf{h}_{s}-\nabla\psi$, which also encodes Ampere’s law, and
invoking the magnetic Gauß law (2) then leads to the nonlinear scalar
potential formulation

$\displaystyle\operatorname{div}(\partial_{\mathbf{h}}w_{*}(\mathbf{h}_{s}-\nabla\psi))$
$\displaystyle=0\qquad$ $\displaystyle\text{in }\Omega$ (6)
$\displaystyle\mathbf{b}\cdot\mathbf{n}$ $\displaystyle=0\qquad$
$\displaystyle\text{on }\partial\Omega.$ (7)

The weak form of this problem is again well-suited for the analysis of the
problem. $H^{1}$-conforming finite elements can now be used for the systematic
discretization and Newton-type methods can be employed for the iterative
solution of the resulting nonlinear systems. We refer to [8, 12] for details
and further information.

######  Remark 1.1.

In three space dimensions, the scalar potential approach achieves the same
approximation order with much fewer degrees of freedom and, additionally,
avoids the need for a complex gauging procedure. The two formulations
mentioned above, however, involve different nonlinear functions and the vector
potential approach seems to require fewer Newton iterations in typical cases
[5].

### Outline of the main results

We start from the mixed formulation (1)–(3), but instead of treating (1) in
strong form, like in previous works [6, 12], we here use the weak form of this
constraint and keep the scalar potential term $\nabla\psi$ in (1) in strong
form. The magnetic Gauß law (2), on the other hand, is treated weakly. Then
$\mathbf{b}$ can be approximated by discontinuous finite elements which allows
to eliminate this variable locally. A key point of our approach is to do this
only after linearization, i.e., in the solution of the linear saddle point
problems arising from the Newton iteration. Overall, we obtain a method that
uses the same linearization $\partial_{\mathbf{b}}w(\mathbf{b})$ of the energy
density as the vector potential formulation, while the linear systems to be
solved in every step of the Newton iteration are of the same structure and can
be solved at the same cost as those of the scalar potential formulation. In
summary, this leads to a modified scalar potential formulation with reduced
computational cost.

##  2 Preliminaries

Let $\Omega\subset\mathbb{R}^{3}$ be a bounded Lipschitz domain and simply
connected. Standard notation for Lebesgue and Sobolev spaces is used
throughout the text; see [13] for details. We write $\langle a,b\rangle$ for
the scalar product in $\mathbb{R}^{3}$ and $\langle
a,b\rangle_{\Omega}=\int_{\Omega}\langle a,b\rangle\,dx$ for the
$L^{2}$-scalar product on $\Omega$. We further assume that
$w:\mathbb{R}^{3}\times\Omega\to\mathbb{R}$ is piecewise smooth in the first
argument and smooth in the second with strongly monotone and Lipschitz
continuous derivative, i.e.,

$\displaystyle\langle\partial_{\mathbf{b}}w(x,\mathbf{b})-\partial_{\mathbf{b}}w(x,\mathbf{b}^{\prime}),\mathbf{b}-\mathbf{b}^{\prime}\rangle$
$\displaystyle\geq\gamma|\mathbf{b}-\mathbf{b}^{\prime}|^{2},$ (8)
$\displaystyle|\partial_{\mathbf{b}}w(x,\mathbf{b})-\partial_{\mathbf{b}}w(x,\mathbf{b}^{\prime})|$
$\displaystyle\leq L|\mathbf{b}-\mathbf{b}^{\prime}|$ (9)

for all $x\in\Omega$ and $\mathbf{b},\mathbf{b}^{\prime}\in\mathbb{R}^{3}$
with uniform constants $L,\gamma>0$. We will not mention the spatial
dependence explicitly in the remainder of the manuscript. As a final
assumption, we require $\mathbf{h}_{s}\in H(\operatorname{curl};\Omega)$ and
set $\mathbf{j}_{s}=\operatorname{curl}\mathbf{h}_{s}$.

### Well-posedness and equivalence

For completeness of the presentation, let us summarize some basic facts about
the well-posedness and mutual connections of the different formulations.

######  Lemma 2.1.

Under the above assumptions, the problem (1)–(3) has a unique solution
$(\mathbf{b},\psi)$ and $\mathbf{h}=\mathbf{h}_{s}-\nabla\psi$ satisfies
$\operatorname{curl}\mathbf{h}=\mathbf{j}_{s}$. Moreover, $\psi$ is also the
unique solution of (6)–(7) and any $\mathbf{a}$ solving (4)–(5) also satisfies
$\mathbf{b}=\operatorname{curl}\mathbf{a}$.

The main arguments to prove existence of a unique solution for the scalar
potential equation can be found in [8]. The remaining details are left to the
reader.

### Discrete mixed problem

Let $\mathcal{T}_{h}$ be a tetrahedral mesh of $\Omega$ satisfying the usual
regularity properties [3, 13].

We will make use of the finite element spaces

$\displaystyle Q_{h}$
$\displaystyle=P_{p-1}(\mathcal{T}_{h})^{d}\quad\text{and}\quad$
$\displaystyle V_{h}$ $\displaystyle=P_{p}(\mathcal{T}_{h})\cap
H^{1}(\Omega)/\mathbb{R},$ (10)

consisting of all vector-valued piecewise polynomials of degree $\leq p-1$ and
all continuous piecewise polynomials of degree $\leq p$ having zero average,
respectively. We further denote by $\langle a,b\rangle_{h}=\sum_{T}\langle
a,b\rangle_{T,h}$ an approximation of the standard $L^{2}$-scalar product
$\langle a,b\rangle_{\Omega}$ obtained by a local quadrature rule applied on
every element. This quadrature rule is assumed to have positive weights and to
integrate polynomials of degree $\leq 2p-2$ exactly.

For the discretization of (1)–(3), we then consider the following method.

######  Problem 2.2.

Find $\mathbf{b}_{h}\in Q_{h}$ and $\psi_{h}\in V_{h}$ such that

$\displaystyle\langle\partial_{\mathbf{b}}w(\mathbf{b}_{h}),\mathbf{b}_{h}^{\prime}\rangle_{h}+\langle\nabla\psi_{h}^{\prime},\mathbf{b}_{h}^{\prime}\rangle_{h}$
$\displaystyle=\langle\mathbf{h}_{s},\mathbf{b}_{h}^{\prime}\rangle_{h}\qquad$
$\displaystyle\forall\mathbf{b}_{h}\in Q_{h}$ (11)
$\displaystyle\langle\mathbf{b}_{h},\nabla\psi_{h}\rangle_{h}$
$\displaystyle=0\qquad$ $\displaystyle\forall\psi_{h}^{\prime}\in V_{h}.$ (12)

With similar arguments as for the continuous problem, we obtain the following.

######  Lemma 2.3.

Under the above assumptions, Problem 2.2 admits a unique solution.

######  Proof 2.4.

Let us note that (11)–(12) corresponds to the first-order optimality condition
for the constrained minimization problem

$\displaystyle\min_{\mathbf{b}_{h}\in Q_{h}}J_{h}(\mathbf{b}_{h}):=\langle
w(\mathbf{b}_{h}),1\rangle_{h}-\langle\mathbf{h}_{s},\mathbf{b}_{h}\rangle_{h}$
(13)
$\displaystyle\qquad\qquad\quad\,\text{s.t.}\quad\langle\mathbf{b}_{h},\nabla\psi_{h}^{\prime}\rangle_{h}=0\quad\forall\psi_{h}^{\prime}\in
V_{h}.$

Under our assumptions, the functional $J_{h}$ to be minimized is strictly
convex and coercive. This guarantees the existence of a unique solution
$\mathbf{b}_{h}\in Q_{h}$. The existence of a unique discrete Lagrange
multiplier $\psi_{h}$ follows from the compatibility of the discretization
spaces; see [2] for related results in the linear case. ∎

### Scalar and vector potential formulations

For later reference let us briefly mention the basic form of the
discretization schemes for the other approaches mentioned in the introduction.
The approximation for the scalar potential problem reads: Find $\psi_{h}\in
V_{h}$ such that

$\displaystyle\langle\partial_{\mathbf{h}}w_{*}(\mathbf{h}_{s}-\nabla\psi_{h}),\nabla\psi_{h}^{\prime}\rangle_{h}$
$\displaystyle=0\qquad\forall\psi_{h}^{\prime}\in V_{h}.$ (14)

The existence of a unique solution again follows by standard arguments [8]. In
our numerical tests, we also compare to a discretization of the vector
potential formulation, which reads: Find $\mathbf{a}_{h}\in W_{h}$ such that

$\displaystyle\langle\partial_{\mathbf{b}}w(\operatorname{curl}\mathbf{a}_{h}),\operatorname{curl}\mathbf{a}_{h}^{\prime}\rangle_{h}$
$\displaystyle=\langle\mathbf{h}_{s},\operatorname{curl}\mathbf{a}_{h}^{\prime}\rangle_{h}\qquad\forall
a_{h}^{\prime}\in W_{h}.$ (15)

Here $W_{h}\subset H_{0}(\operatorname{curl};\Omega)$ is an appropriate sub-
space of edge elements incorporating the required gauging conditions. The
existence of a unique discrete solution can then again be established by well-
known arguments [9].

##  3 Iterative solution algorithm

For the iterative solution of (11)–(12), we consider a Newton method with line
search. Every step of this iteration has the form

$\displaystyle(\mathbf{b}_{h}^{n+1},\psi_{h}^{n+1})$
$\displaystyle=(\mathbf{b}_{h}^{n},\psi_{h}^{n})+\tau^{n}(\delta\mathbf{b}_{h}^{n},\delta\psi_{h}^{n})$
(16)

with increments $(\delta\mathbf{b}_{h},\delta\psi_{h})\in Q_{h}\times V_{h}$
determined by the linearized systems

$\displaystyle\langle\partial_{bb}w(\mathbf{b}_{h}^{n})\delta\mathbf{b}_{h}^{n},\mathbf{b}_{h}^{\prime}\rangle_{h}+\langle\nabla\delta\psi_{h}^{n},\mathbf{b}_{h}^{\prime}\rangle_{h}$
$\displaystyle=-\langle\partial_{b}w(\mathbf{b}_{h}^{n})+\nabla\psi_{h}^{n}-\mathbf{h}_{s},\mathbf{b}_{h}^{\prime}\rangle_{h}$
(17)
$\displaystyle\langle\delta\mathbf{b}_{h}^{n},\nabla\psi_{h}^{\prime}\rangle_{h}$
$\displaystyle=-\langle\mathbf{b}_{h}^{n},\nabla\psi_{h}\rangle_{h}.$ (18)

for all test functions $\mathbf{b}_{h}^{\prime}\in Q_{h}$ and
$\psi_{h}^{\prime}\in V_{h}$. With standard results about linear saddlepoint
problems [2], one can show the following result.

######  Lemma 3.1.

For any given $(\mathbf{b}_{h}^{n},\psi_{h}^{n})\in Q_{h}\times V_{h}$, the
system (17)–(18) has a unique solution
$(\delta\mathbf{b}_{h}^{n},\delta\psi_{h}^{n})\in Q_{h}\times V_{h}$.
Moreover, if $\langle\mathbf{b}_{h}^{n},\nabla\psi_{h}^{\prime}\rangle_{h}=0$
for all $\psi_{h}^{\prime}\in V_{h}$, then one also has
$\langle\delta\mathbf{b}_{h}^{n},\nabla\psi_{h}^{\prime}\rangle_{h}=0$ for all
$\psi_{h}^{\prime}\in V_{h}$.

######  Remark 3.2.

Together with an appropriate line search strategy, one can now establish
global convergence of the iteration (16)–(17) using well-known arguments from
nonlinear optimization. In our numerical tests, we use Armijo backtracking
[11], which allows to guarantee sufficient decay in the energy (13) in every
step. If assuming the initial choice $\mathbf{b}_{h}^{0}=0$, global linear
convergence with a mesh-independent convergence rate follows from the
arguments in [7]. Similar results can also be proven for the other
discretization strategies discussed above.

##  4 Efficient realization of the Newton step

In every step of the proposed iterative method, a linearized discrete saddle
point problem (17)–(18) has to be solved. On the algebraic level, it has the
form

$\displaystyle\begin{pmatrix}A&\,\,B^{\top}\\\
B&0\end{pmatrix}\begin{pmatrix}db\\\ d\psi\end{pmatrix}=\begin{pmatrix}f\\\
g\end{pmatrix}$ (19)

To avoid misunderstanding, we use slightly different notation for the
coordinate vectors $(db,d\psi)$ than for the corresponding functions
$(\delta\mathbf{b}_{h}^{n},\delta\psi_{h}^{n})$. Due to our assumptions, the
matrix $A$ can be shown to be positive definite, and since the finite element
space $Q_{h}$ carries no continuity properties, it has a block-diagonal
structure. One can thus efficiently assemble the Schur complement system

$\displaystyle BA^{-1}B^{\top}d\psi$ $\displaystyle=BA^{-1}f-g$ (20)

whose system matrix $S=BA^{-1}B^{\top}$ corresponds to the stiffness matrix of
a standard finite element discretization for the Poisson problem. In
particular, it has the same properties as that of the discretized and
linearized scalar potential problem (14). In a second step, one can determine
the second component

$\displaystyle db=A^{-1}(f-B^{\top}d\psi),$ (21)

which can again be computed efficiently, since $A^{-1}$ is still block
diagonal. In summary, the Newton-update $(db,d\psi)$ can be computed with a
computational effort comparable to that of the standard scalar potential
formulation.

##  5 Numerical results

To illustrate our theoretical results, we now present some computational tests
in which we compare the proposed methods with the standard scalar- and vector
potential formulations. As a test problem, we use Problem 13 of the Compumag
TEAM suite [4]. The geometry for this problem consists of a coil surrounded by
steel plates; see Figure 1 for a sketch.

Figure 1: Bird’s eye view of the geometry of TEAM problem 13. Pink volumes are
the steel plates; the green area is the coil. The ambient air box is hidden.

#### Material laws.

In the air and the coil, the magnetic response is modeled by the linear
material relation $\mathbf{b}=\mu_{0}\mathbf{h}$, where $\mu_{0}$ is the
vacuum permeability. This corresponds to
$w(\mathbf{b})=\frac{2}{2\mu_{0}}|\mathbf{b}|^{2}$ respectively
$w_{*}(\mathbf{h})=\frac{\mu_{0}}{2}|\mathbf{h}|^{2}$. For the steel plates we
us an isotropic nonlinear relation
$w(\mathbf{b})=\widetilde{w}(|\mathbf{b}|)$. Then
$\widetilde{w}^{\prime}(|\mathbf{b}|)=|\mathbf{h}|$ can be obtained from the
B-H data provided in [4, Figure 3]. The functions
$\widetilde{w}(|\mathbf{b}|)$ and $\widetilde{w}^{\prime\prime}(|\mathbf{b}|)$
which are required to compute $w(\mathbf{b})$ and
$\partial_{\mathbf{bb}}w(\mathbf{b})$ for the nonlinear iteration can then be
obtained by integration resp. differentiation. Validity of the conditions
(8)–(9) can be guaranteed explicitly by this procedure. The setup of the
coenergy density $w_{*}(\mathbf{h})=\widetilde{w}_{*}(|\mathbf{h}|)$ is done
in a similar manner. The solution of our test problem is depicted in Figure 2.

Figure 2: Solution to the magnetostatic problem using our mixed approach.
Black arrows depict the direction of the current in the coil. Scalar value is
the magnitude of the $\mathbf{b}$-field. (computed using Netgen/NGSolve and
displayed using PyVista.

#### Implementation.

We conduct a series of numerical tests using the vector and scalar potential
formulations and the method proposed in this paper. Results for different mesh
sizes $h$ and polynomial degree $p$ are presented. Since singularities occur
at material interfaces, we refine the mesh towards the respective faces and
edges using local mesh sizes $h_{f}=2^{-3}\cdot h$ and $h_{e}=2^{-4}\cdot h$.
The global mesh size is chosen as $h=2^{-k}$, $k=2,3,\ldots$. The
corresponding meshes are generated independently from each other. All tests
were implemented in Netgen/Ngsolve [14]. The same nonlinear solution strategy,
i.e. a Newton method with Armijo line search, was used for all methods with
the same parameters and stopping tolerances. The linear systems arising in
every Newton step were solved with the conjugate gradient method with Jacobi
preconditioner. Again, the same tolerances were used for stopping these
iterations.

#### Numerical results.

In Table 1, we report on the iteration numbers of the nonlinear solver
required for the different methods and for different choices of the excitation
current, the mesh sizes $h$, and the polynomial degrees $p$.

1000 ampere-turns
---
$h$ | $p=1$ | $p=2$ | $p=3$ | $p=4$
$2^{-2}$ | $8\,|\,9\,|\,8$ | $6\,|\,14\,|\,8$ | $7\,|\,13\,|\,8$ | $7\,|\,15\,|\,8$
$2^{-3}$ | $8\,|\,12\,|\,8$ | $6\,|\,16\,|\,7$ | $7\,|\,15\,|\,8$ | $7\,|\,14\,|\,8$
$2^{-4}$ | $9\,|\,15\,|\,7$ | $7\,|\,19\,|\,8$ | $7\,|\,16\,|\,8$ | $7\,|\,15\,|\,8$
$2^{-5}$ | $7\,|\,17\,|\,7$ | $7\,|\,16\,|\,8$ | $7\,|\,20\,|\,8$ | $7\,|\,19\,|\,8$
$2^{-6}$ | $7\,|\,15\,|\,7$ | $7\,|\,23\,|\,7$ | $7\,|\,22\,|\,8$ | $8\,|\,18\,|\,9$

3000 ampere-turns
---
$h$ | $p=1$ | $p=2$ | $p=3$ | $p=4$
$2^{-2}$ | $8\,|\,13\,|\,8$ | $9\,|\,18\,|\,7$ | $8\,|\,16\,|\,9$ | $8\,|\,16\,|\,9$
$2^{-3}$ | $8\,|\,14\,|\,8$ | $9\,|\,13\,|\,8$ | $8\,|\,18\,|\,8$ | $8\,|\,16\,|\,8$
$2^{-4}$ | $8\,|\,12\,|\,8$ | $8\,|\,14\,|\,7$ | $8\,|\,15\,|\,8$ | $8\,|\,20\,|\,8$
$2^{-5}$ | $9\,|\,17\,|\,7$ | $8\,|\,21\,|\,7$ | $8\,|\,27\,|\,8$ | $8\,|\,21\,|\,8$
$2^{-6}$ | $8\,|\,21\,|\,7$ | $8\,|\,21\,|\,8$ | $8\,|\,20\,|\,8$ | $8\,|\,20\,|\,8$

Table 1: Number of iterations for the vector potential formulation, the scalar
potential approach, and the proposed modified scalar potential approach (left
to right).

We observe that in general, the iteration numbers are more or less independent
of the discretization parameters; see [7] for explanation. The two energy-
based methods (left, right) require only about $7$–$8$ iterations to reach the
desired tolerance, while the scalar potential formulation (middle), which is
based on the coenergy functional, requires approximately twice as many. This
is in perfect agreement with the observations made in [5]. Let us note that,
for the same mesh and polynomial order, the vector potential formulation has
about $3$–$6$ times as many degrees of freedom as the corresponding scalar
potential formulation, and in our computations, the solution of the linear
systems, therefore, takes $3$–$6$ times longer. For the proposed modified
scalar potential approach, the linear solves can be done at the same
complexity as for the standard scalar potential formulations. The number of
Newton iterations, and hence also the total computational complexity, is
reduced by about one half.

### Acknowledgement

The authors are grateful for the financial support from the international
FWF/DFG-funded Collaborative Research Centre CREATOR (TRR361/SFB-F90).

## References

  * [1] R. Albanese and G. Rubinacci, Finite element methods for the solution of 3D eddy current problems, Elsevier, 1997, pp. 1–86.
  * [2] D. Boffi, F. Brezzi, and M. Fortin, Mixed finite element methods and applications, Springer, Heidelberg, 2013.
  * [3] P. G. Ciarlet, The finite element method for elliptic problems, vol. 40 of Classics in Applied Mathematics, SIAM, Philadelphia, PA, 2002.
  * [4] Compumag, Team problem 13, www.compumag.org/wp/team/.
  * [5] J. Dular, C. Geuzaine, and B. Vanderheyden, Finite-element formulations for systems with high-temperature superconductors, IEEE Trans. Appl. Supercond., 30 (2020), pp. 1–13.
  * [6] P. Dular, J.-F. Remacle, F. Henrotte, A. Genon, and W. Legros, Magnetostatic and magnetodynamic mixed formulations compared with conventional formulations, IEEE Trans. Magn., 33 (1997), pp. 1302–1305.
  * [7] H. Egger, F. Engertsberger, and B. Radu, Global convergence of iterative solvers for problems of nonlinear magnetostatics, arXiv:2403.18520, (2024).
  * [8] F. Engertsberger, The Scalar Potential Approach in Nonlinear Magnetostatics.  MSc Thesis, JKU Linz, 2023.
  * [9] B. Heise, Analysis of a fully discrete finite element method for a nonlinear magnetic field problem, SIAM J. Numer. Anal., 31 (1994), pp. 745–759.
  * [10] R. Hiptmair, Finite elements in computational electromagnetism, Acta Numerica, 11 (2002), pp. 237–339.
  * [11] C. T. Kelley, Solving nonlinear equations with Newton’s method, Fundamentals of Algorithms, SIAM, 2003.
  * [12] G. Meunier, The Finite Element Method for Electromagnetic Modeling, ISTE Wiley, 2008.
  * [13] P. Monk, Finite element methods for Maxwell’s equations, Oxford University Press, New York, 2003.
  * [14] Netgen/NgSolve.  https://ngsolve.org/.
  * [15] J. C. Nédélec, Mixed finite elements in $\mathbb{R}^{3}$, Numer. Math., 35 (1980), pp. 315–341.
  * [16] P. P. Silvester and R. P. Gupta, Effective computational models for anisotropic soft B-H curves, IEEE Trans. Magn., 27 (1991), pp. 3804–3807.