# Constraining the geometry of the reflection nebula NGC 2023 with [O i]:
Emission & Absorption

Bhaswati Mookerjea1, Göran Sandell2, Rolf Güsten3, Helmut Wiesemeyer3, Yoko
Okada 4, Karl Jacobs4
1Department of Astronomy & Astrophysics, Tata Institute of Fundamental
Research, Homi Bhabha Road, Mumbai 400005, India
2 Institute for Astronomy, University of Hawai‘i at Manoa, 640 N. Aohoku
Place, Hilo, HI 96720, USA
3 Max Planck Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn,
Germany
4 I. Physikalisches Institut, Universität zu Köln Zülpicher Str. 77, 50937
Köln, Germany
E-mail<EMAIL_ADDRESS>
###### Abstract

We have mapped the NGC 2023 reflection nebula in the 63 and 145 µm transitions
of [O i] and the 158 µm [C ii] spectral lines using the heterodyne receiver
upGREAT on SOFIA. The observations were used to identify the diffuse and dense
components of the PDR traced by the [C ii] and [O i] emission, respectively.
The velocity-resolved observations reveal the presence of a significant column
of low-excitation atomic oxygen, seen in absorption in the [O i] 63 µm
spectra, amounting to about 20–60% of the oxygen column seen in emission in
the [O i] 145 µm spectra. Some self-absorption is also seen in [C ii], but for
the most part it is hardly noticeable. The [C ii] and [O i] 63 µm spectra show
strong red- and blue-shifted wings due to photo evaporation flows especially
in the southeastern and southern part of the reflection nebula, where
comparison with the mid- and high-$J$ CO emission indicates that the C+ region
is expanding into a dense molecular cloud. Using a two-slab toy model the
large-scale self-absorption seen in [O i] 63 µm is readily explained as
originating in foreground low-excitation gas associated with the source.
Similar columns have also been observed recently in other Galactic photon-
dominated-regions (PDRs). These results have two implications: for the
velocity-unresolved extra-galactic observations this could impact the use of
[O i] 63 µm as a tracer of massive star formation and secondly the widespread
self-absorption in [O i] 63 µm leads to underestimate of the column density of
atomic oxygen derived from this tracer and necessitates the use of alternative
indirect methods.

###### keywords:

ISM: clouds – ISM: kinematics and dynamics – submillimetre: ISM – ISM:
structure – stars: formation – ISM:individual (NGC 2023)

††pubyear: 2023††pagerange: Constraining the geometry of the reflection nebula
NGC 2023 with [O i]: Emission & Absorption–9

##  1 Introduction

The fine-structure lines of [C ii] at 158 $\mu$m and [O i] at 63 and 145
$\mu$m are the most important cooling lines in the far-infrared (FIR) and have
long been used to study photon-dominated regions (PDRs). In the rest of the
paper we refer to the [O i] 63 µm and [O i] 145 µm fine-structure lines simply
as [O i] 63 and [O i] 145. PDRs are regions where far-ultraviolet (FUV; 6 eV
$<$ h$\nu$ 13.6 eV) radiation from young massive stars dominates the physics
and the chemistry of the interstellar medium and correspond to the transition
from ionized to molecular gas (Hollenbach & Tielens, 1997). PDRS were first
studied extensively with the Infrared Space Observatory, though the lack of
spectral and spatial resolution was a severe limitation. This changed with the
Kuiper Airborne Observatory (in limited regions) and with the Heterodyne
Instrument for the Far Infrared (HIFI) on Herschel, which had both velocity
resolution and sensitivity to even enable determination of the optical depth
of [C ii] by also observing 13[C ii] (Ossenkopf et al., 2013). After Herschel
this work was taken over by first GREAT on SOFIA and later with the array
receiver upGREAT. These studies have found that the [C ii] emission often can
be moderately optically thick and sometimes significantly self-absorbed
(Mookerjea et al., 2018; Guevara et al., 2020; Graf et al., 2012), which may
underestimate [C ii] column densities by a factor of two to three. The [O i]
63 µm can also be optically thick and self-absorbed, sometimes even more
strongly than [C ii] (see e. g., Mookerjea et al., 2021).

NGC 2023, illuminated by the B2 star HD 37903, is one of the best-studied
reflection nebulae. It is nearby, $\sim$400 pc, and exhibits a strong, nearly
edge-on PDR in which the cavity expands into the surrounding dense molecular
cloud. It has served as a test bed for exploring PDR models (e.g. Kaufman et
al., 2006). It was mapped with GREAT on SOFIA in [C ii] and CO(11–10) with the
single pixel L1/L2 mixer (Sandell et al., 2015), who also mapped the nebula in
13CO(3-2), and various transitions of CO from J=3–2 up to J=7–6.

Here we have additionally mapped the NGC 2023 PDR in the 63 and 145 $\mu$m
fine-structure transitions of [O i] and at the same time obtained a deeper map
of [C ii] using upGREAT on SOFIA. The [O i] observations enable us for the
first time to study the physical conditions in the south eastern and southern
part of the nebula, where the C+ region is expanding into the dense
surrounding molecular cloud as well as estimate the relative contribution of
the dense and diffuse PDR gas to the [C ii] emission.

The newly obtained observations are analyzed in combination with the
previously published maps of (i) 13CO(3–2), CO(6–5), CO(7–6) observed with
APEX with beamsizes of 18$\aas@@fstack{\prime\prime}$5,
9$\aas@@fstack{\prime\prime}$1 and 7$\aas@@fstack{\prime\prime}$7 respectively
and (ii) CO(11–10) observed with SOFIA/GREAT with a beamsize of 23″ (Sandell
et al., 2015).

##  2 Observations

The reflection nebula NGC 2023 was observed on a 90 minute leg with the SOFIA
upGREAT111The German REceiver for Astronomy at Terahertz frequencies (upGREAT)
is a development by the MPI für Radioastronomie and the KOSMA/Universität zu
Köln, in cooperation with the DLR Institut für Optische Sensorsysteme.
receiver (Risacher et al., 2016) in GREAT consortium time on SOFIA flight #668
(Project ID: 83_0731) out of Palmdale, CA, on March 6, 2020. Although NGC 2023
was observed at low flight altitude, 11.6 km (38,000 feet), the observing
conditions were very good with a precipitable water vapor (PWV) around 5–10
µm. The upGREAT was in the Low Frequency Array/High Frequency Array (LFA/HFA)
configuration with both arrays operating in parallel. The V polarization of
the LFA was tuned to [C ii] at 1.9005369 THz while the H polarization was
tuned to the [O i] 145 µm line at 2.06006886 THz. The HFA was tuned to the [O
i] 63 µm at 4,744.77749 THz. We made total power on-the fly (OTF) maps in
classical position switched mode. The reference position was at +390″,$-$70″
relative to HD 37903 ( ($\alpha_{J2000}$: $05^{h}~{}41^{m}~{}38.39^{s}$,
$\delta_{J2000}$: $-02^{\circ}15^{\prime}$
32${}^{\prime\prime}\mskip-7.6mu.\,$5). The OTF map was done in two coverages,
scanning in both x and y. The map was rotated by 45∘. The sampling was done
every 3″ with a sampling time of 0.3 second per dump. This resulted in maps
for the LFA array of $\sim$ 4$\aas@@fstack{\prime}$9 $\times$
3$\aas@@fstack{\prime}$9, while the map size for the HFA array was
3.0′$\times$ 2$\aas@@fstack{\prime}$1, which was enough to cover the SE
quadrant of the reflection nebula, where [O i] was expected to be strong.

The observations were reduced and calibrated by the GREAT team. The GREAT team
also provided beam sizes (14${}^{\prime\prime}\mskip-7.6mu.\,$1 for [C ii],
13${}^{\prime\prime}\mskip-7.6mu.\,$0 for [O i] 145 µm, and
6${}^{\prime\prime}\mskip-7.6mu.\,$3 for [O i] 63 µm) and beam efficiencies
derived from planet observations. The data were corrected for atmospheric
extinction, flat-fielded and calibrated in $T_{\rm mb}$. Further processing of
the data was made with the CLASS222CLASS is part of the Grenoble Image and
Line Data Analysis Software (GILDAS), which is provided and actively developed
by IRAM, and is available at http://www.iram.fr/IRAMFR/GILDAS software.

We also retrieved some high quality [C ii] data from the SOFIA archive, which
extended the [C ii] map to the north. This data set was from project 02_0090
(PI. Els Peeters) and fully calibrated in the SOFIA archive. All data in this
paper are presented in main-beam temperature scale.

Figure 1: Integrated intensity images (in both color and contours) of [C ii]
at 158 µm (Left) and [O i] at 145 µm (Right). The color scale in K km s-1 for
each panel shown on the top. Contour levels (in K km s-1) for [C ii] are 40,
46, 65, 70, 76, 83, 90, 98, 109, 131, 141, 150, 160, 169 and 173. Contour
levels (in K km s-1) for [O i] 145 are 3, 5 to 50 in steps of 5. The
coordinates are shown as offsets relative to the center RA: 05h41m38.39s Dec:
-02∘15′32$\aas@@fstack{\prime\prime}$5, which is also the location of the
illuminating star HD 37903 marked by a star symbol. Positions which are
studied in detail are shown as ’+’. The dashed lines drawn in the left panel
show the directions along which position-velocity diagrams have been studied
in Fig.5. The beamsizes of the [C ii] and [O i] 145 maps are
14$\aas@@fstack{\prime\prime}$1 and 13″ respectively and are shown hatched
circles in each panel.

##  3 Results

###  3.1 Morphology & Kinematics of the Region

Figure 2: Integrated intensity image (color) of integrated intensities of [O
i] 63 µm overlaid with contours of [O i] 145 µm intensities. The color scale
is in K km s-1 and beamsize for the [O i] 63 µm map is shown as hatched
circles in the left bottom corner of the figure. Contour levels (in K km s-1)
for [O i] 145 µm are 3, 5 to 50 in steps of 5. Rest of the details are the
same as in Fig. 1.

Figure 3: Integrated intensity images (color) of CO(7–6), CO(11–10) and
C91$\alpha$ compared with contours of intensities of [O i] 145 µm. The color
scale is in K km s-1. The beamsizes for the CO(7–6), CO(11–10) and C91$\alpha$
maps (in color) are 7$\aas@@fstack{\prime\prime}$7, 23″ and 20″ respectively
and shown on top right corner of each panel as a filled white circle. The data
for the CO(7–6), CO(11–10) emission observed with APEX and GREAT/SOFIA
respectively, are from Sandell et al. (2015) and the C91$\alpha$ map is taken
from Wyrowski et al. (2000). Contour levels (in K km s-1) for [O i] 145 are 3,
5 to 50 in steps of 5. Rest of the details are the same as in Fig. 1.

Sandell et al. (2015) used maps of multiple CO rotational transitions combined
with the [C ii] fine-structure line at 158 µm to probe the overall morphology
of the reflection nebula NGC 2023. These observations showed that the [C ii]
emission traces an expanding ellipsoidal shell-like PDR region region, which
is surrounded by a hot molecular shell. In the southeast, where the [C ii]
region expands into a dense, clumpy molecular cloud ridge, they also detected
emission from high-$J$ CO lines, apparently originating in a thin, hot
molecular shell surrounding the [C ii] emission. These authors found that
there was a clear velocity gradient across the nebula, with the emission being
more blue-shifted in the south and southeast and more red-shifted in the north
and northwest, some of which appeared to be due to an expansion of the nebula.
Sandell et al. also noted that high angular resolution images of vibrationally
excited H2 and PAH emission suggested that the PDR was far from smooth,
exhibiting lumpy ridges and bright filamentary structures. The [C ii] and
high-$J$ CO maps looked smoother, partly because they had insufficient spatial
resolution to see such details. Here we have additionally mapped the fine-
structure transitions of [O i] at 63 and 145 µm, which are good tracers of the
warm and dense regions of the PDR. The upper energy levels of the 63 µm and
145 µm [O i] transitions are 227.7 K and 326.6 K respectively, while formally
the critical densities for these transitions are 5 105 cm-3 and 5 106 cm-3
respectively, for collisions with H2 (Goldsmith, 2019). The [O i] 63 emission,
however, can be very optically thick and self-absorbed, which was seen toward
the S1 PDR in RhoOph (Mookerjea et al., 2021) and the same is true for NGC
2023. However, here we have also observed [O i] 145, which is unaffected by
self-absorption. The [O i] 145 emission is largely optically thin and shows
that the dense PDR has a lumpy and filamentary structure. The [O i] 145 map
looks qualitatively similar to the [C ii] map (Fig. 1), but the [O i] emission
is concentrated to the dense south-eastern and southern part of the PDR and
more structured, while the [C ii] emission is rather smooth. In the north and
northwest the gas densities are too low to collisionally excite [O i]. The [O
i] 63 map, which has a spatial resolution of
6${}^{\prime\prime}\mskip-7.6mu.\,$3, shows that the emission is extremely
patchy over the area where it is detected (Fig. 2). The emission is
essentially completely self-absorbed in the southeast and south, while it is
seen in emission in the east and northeast, where the absorption from the
foreground cloud is less severe. Even here it is strongly affected by
absorption from colder foreground gas giving the integrated line emission a
very patchy and lumpy appearance.

Figure 4: Top panels: Velocity channel maps for [O i] 145 µm in color overlaid
with contours at 3 to 30 K in steps of 2 K. Bottom panels: Velocity channel
maps (in color) for [C ii] overlaid with the [O i] 145 µm map in contours from
3 K to 30 K in steps of 2 K. The color scale for the channel maps is shown at
the top of each set. The position of HD 37903 is shown as a “star". The
coordinates in both panels are offsets in RA and Dec relative to the position
of HD 37903 as in Fig. 1.

We also compare the [O i] 145 map with maps of CO(7–6) and CO(11–10) (Sandell
et al., 2015), which also trace the hot, dense PDR as well as with
C91$\alpha$, which was mapped by Wyrowski et al. (2000) (Fig. 3).

The overall morphology of [O i] 145 matches CO(7–6) and CO(11–10) quite well.
Both appear to trace hot gas in the leading edge of the expanding PDR shell.
That CO(7–6) emission is not seen from the entire PDR where [O i] is detected,
is likely due to the fact that the emission is strong enough only from a
narrow shell and needs to be viewed almost tangentially to have enough column
density to be detected. The emission from the remaining hot gas is likely
below the sensitivity of our observations. C91$\alpha$ shows a different
morphology than CO(11–10), although they both appear to originate in
approximately the same region. However, where C91$\alpha$ is strong, CO(11–10)
is faint, like in the southern part of the PDR shell, where C91$\alpha$ shows
two strong peaks. Both of these peaks are seen in [O i] 145, which appear to
trace the same gas as C91$\alpha$. In the east, where C91$\alpha$ is faint
there is no clear peak in [O i] 145 either. We note that although the [O i]
145 traces the high density ridge and the clumps in the ridge, the emission is
slightly offset compared to CO(11–10) and C91$\alpha$, which appears to be
further outside in the PDR shell, i.e., closer to the cold molecular cloud
(Fig. 3). Based on the detailed analysis of the high-$J$ CO lines, Sandell et
al. (2015) estimated the hot molecular shell to be between 90–120 K and have
densities $n\sim 10^{5}$–106 cm-3.

The velocity-channel map of [O i] 145 emission clearly show that the dense PDR
clumps along the ridge are also at different velocities with the south-western
clumps being more red-shifted than the north-eastern clump (Fig. 4). At v = 9
km s-1 the region northeast of HD 37903 dominates the [O i] 145 and [C ii]
emission. This emission peak, at $\sim$ +27″, +18″, coincides with the
strongest [C ii] emission peak, $\sim$ 55 K, in the whole map. It is also seen
as a distinct peak in CO(6–5), CO(7–6) and CO(11–10), but it is relatively
faint compared to the emission at $\sim$ v = 10.5 km s-1, which dominates the
emission in the southeastern part of the PDR. This suggests that the
northeastern emission peak may have somewhat lower gas densities than the PDR
emission in the ridge.

Figure 5: Position-velocity diagrams for [C ii], 13CO(3–2), CO(4–3), CO(6–5)
and [O i] 145 (from top to bottom in each column) along directions given by
position angles PA of 0, 90∘ and 110∘ as well as along the ridge. The position
angles are measured counterclockwise from north. The colorscales (in K) are
shown next to each panel. All the plots are at their respective original
resolutions.

Optical images show that the NGC 2023 reflection nebula is quite large with a
size of $\sim$ 10′ $\times$ 10′. In the near- and mid-IR it is only about half
the size, $\sim$ 5$\aas@@fstack{\prime}$2 $\times$ 5$\aas@@fstack{\prime}$3
(Mookerjea et al., 2009). This is the regime where the emission is dominated
by the PDR. In the near-IR one can see that HD 37903 illuminates a bright
ridge $\sim$ 120″ – 150″ to the north, which intercepts some of the FUV
radiation. This ridge is dominated by fluorescent H2 emission, see e.g, the
narrowband H2 imaging by Field et al. (1998). The feature Field et al. calls
the IR triangle is quite bright in [C ii] (Fig. 1). The peak of the [C ii]
emission (+27″, +18″) coincides with the head of what Field et al. refer to as
the Seahorse. It is clear that the PDR shell opens up to the north, as well as
in a region east of the IR triangle, Fig. 1. Even in the southeast, where the
C+ region is surrounded by the dense molecular cloud, one can see that FUV
radiation escapes, like the tongue of [C ii] emission sticking out at a
position angle (measured counter-clockwise from North) of $\sim$ 110° (Fig.
1).

The position-velocity diagrams, Fig. 5, provides us with an alternative way to
examine the morphology of the PDR and what the different PDR tracers show us.
There is a clear velocity gradient from south to north with the emission in
the northern part of the PDR being more red-shifted, seen in all the tracers,
even [O i] 145, although the [O i] 145 map only partially covers the PDR
ridge. This suggest that the IR triangle and the northern ridge is on the
backside of the PDR. Near HD 37903 the [C ii] emission shows both blue-and
red-shifted emission, suggesting that the whole C+ region is embedded in the
cloud, because one sees photo evaporation flows both from the front and the
back side of the PDR. The CO(4–3) and CO(6–5) show an outflow just south of HD
37903 and strong blue-and red-shifted emission, where the position-velocity
diagram crosses the outflows from Sellgren C and D. Although [C ii] also shows
strong blue-shifted emission in the south, it is almost certainly due to
strong photo evaporation, as one can see from the velocity-channel maps (Fig.
4). These maps show strong blue-shifted emission over the whole PDR in the
southeast. The east-west position-velocity diagram shows the strong blue-
shifted photo-evaporation flow inside the PDR shell and it appears that some
FUV radiation “leaks” through the shell, i.e., one can see [C ii] emission at
distances $>$ 100″ from HD 37903 at roughly the same velocity, $\sim$ 9.9 km
s-1. Closer to the star the [C ii] emission appears to be somewhat self-
absorbed. [O i] 145 shows blue-shifted emission east of HD 37903, indicating
that the photo evaporating gas is so dense that it is seen even in [O i] 145.
To the west, the [O i] 145 becomes more red-shifted, which is not really
noticeable in [C ii]. The CO lines pick up the blue-shifted outflow slightly
west of HD 37903. The position-velocity diagram at PA 110° shows one of the
regions where [C ii] breaks through the PDR shell. However, it looks about the
same as in the east-west cut. Towards the west-northwest the [C ii] splits up
into two velocity components. However, both [O i] 145 and 13CO(3–2) only show
a single velocity component roughly halfway between the two [C ii] peaks,
indicating that [C ii] must be self-absorbed.

The position-velocity diagram along the ridge cuts through the two strong
C91$\alpha$ peaks. The southern one is faint in [C ii] while it is quite
prominent in [O i] 145. Here the emission is shifted toward $\sim$ 10.7 km
s-1, whereas the clump to the northeast of it is $\sim$ 10.2,km s-1 (Fig. 5).
The velocity of the [O i] 145 emission appears to follow that of C91$\alpha$
(Wyrowski et al., 2000), while these emission clumps are unnoticeable in [C
ii]. The CO lines show strong blue-shifted emission from a molecular outflow
powered by the young star Sellgren D.

###  3.2 Spectral Profiles of [O i] lines

Based on the model proposed by (Sandell et al., 2015) the spectral profiles of
CO and [C ii] detected NGC 2023 consist of several velocity components: the
quiescent cloud seen in low J CO lines and 13CO, the PDR, and line wings from
high velocity molecular outflows seen in the outskirt of the reflection
nebula. The [C ii] spectra also show prominent blue or red-shifted wings from
photo evaporation flows. In order to see which of these components contribute
to the [O i] emissions, we have selected six positions, some of which already
were analyzed by (Sandell et al., 2015). All spectra are smoothed to a common
angular resolution of 15″ (Fig. 6). We have also fitted the [O i] 145 spectra
using Gaussian profiles consisting of one or two velocity components as
appropriate and compared the fitted parameters with those derived for [C ii]
(this work) and CO(6–5), (7–6) and (11–10) transitions (Table 1). We have not
used [O i] 63 µm, because the line is strongly self-absorbed over the whole
area that we have mapped. The [O i] 145 spectra are all single-peaked,
although we do see clear blue-shifted line wings in several positions, like at
40″, -40″ and 23″, -60″, where [C ii] shows very strong blueshifted wings. The
position 23″, -60″, which coincides with the C91$\alpha$ emission peak 2, also
shows blue-shifted emission in [O i] 145, and even in C91$\alpha$. At position
60″, -60″, the [C ii] line is strongly self-absorbed at the velocity of [O i]
145, making the blue-shifted line wing appear stronger than the emission from
the PDR (Table 1). Sandell et al. (2015) interpreted the double-peaked [C ii]
spectra that they saw in their [C ii] map as emission coming from the front
and the backside of the PDR shell, but the optically thin [O i] 145 spectra
show that the dip in the [C ii] spectra is due to self-absorption, not two
separate velocity components. Neither can we distinguish more than one
velocity component in the CO spectra, although there are areas, see Fig. 6 and
Table 1, where the lines appear broadened, most likely due to contribution
both from the front and the back wall of the PDR,

In most of the positions the [O i] 63 spectra show wings and line-widths which
are comparable to the [C ii] spectra, although they are often completely self-
absorbed at the center of the line. The [O i] 145 spectra are much narrower
than the [O i] 63 µm and show profiles which are almost identical to the
CO(7–6) emission, except at position -2″, -78″ where the CO(7–6) shows strong
blue-shifted from the outflow powered by Sellgren’s star C (Sandell, priv.
communication).

Figure 6: Comparison of the [O i] 145 µm spectra at selected positions in the
NGC 2023 region (marked in Fig. 1) with spectra of [C ii], [O i] at 63 µm,
CO(11–10) and CO(7–6). As indicated the spectra are scaled by arbitrary
factors for better visibility. The temperature scale is in $T_{\rm mb}$(K).
The CO(7–6) and CO(11–10) are from Sandell et al. (2015). The CO(11–10)
spectra are at the native resolution of 23″, all other spectra are at a common
resolution of 15 ″.  Table 1: Gaussian fits to spectra at selected positions.

Offset | Tracer | $\int T_{\rm mb}dv$ |  $V_{\rm LSR}$00 |  $\Delta V$00 | $\int T_{\rm mb}dv$ |  $V_{\rm LSR}$00 |  $\Delta V$00
---|---|---|---|---|---|---|---
(″,″) |  | (K km s-1) | (km s-1) | (km s-1) | (K km s-1) | (km s-1) | (km s-1)
(0.0,0.0) | [C ii]a | 113.7$\pm$1.2 | 10.24$\pm$0.02 | 3.22$\pm$0.04 |  |  |
|  13CO(3–2) | 32.0$\pm$0.4 | 10.24$\pm$0.01 | 1.61$\pm$0.02 |  |  |
| [O i] 145 | 27.2$\pm$1.2 | 10.48$\pm$0.04 | 2.16$\pm$0.10 |  |  |
| CO(6–5) | 132.0$\pm$0.5 | 10.25$\pm$0.01 | 2.38$\pm$0.01 |  |  |
| CO(7–6) | 85.0$\pm$0.8 | 10.24$\pm$0.01 | 2.02$\pm$0.02 |  |  |
| CO(11–10) | 4.4$\pm$0.2 | 10.51$\pm$0.05 | 1.91$\pm$0.13 |  |  |
(27,18) | [C ii]b | 159.7$\pm$3.9 | 9.10$\pm$0.02 | 2.02$\pm$0.04 | 51.5$\pm$4.1 | 11.19$\pm$0.06 | 1.97$\pm$0.12
|  13CO(3–2) | 32.9$\pm$0.5 | 10.11$\pm$0.02 | 2.35$\pm$0.04 |  |  |
| [O i] 145 | 43.2$\pm$1.3 | 9.89$\pm$0.04 | 2.50$\pm$0.09 |  |  |
| CO(6–5) | 128.6$\pm$0.3 | 10.05$\pm$0.01 | 2.77$\pm$0.01 |  |  |
| CO(7–6) | 80.3$\pm$1.0 | 10.07$\pm$0.01 | 2.51$\pm$0.03 |  |  |
| CO(11–10) | 7.0$\pm$0.4 | 10.63$\pm$0.06 | 2.10$\pm$0.15 |  |  |
(40,-40) | [C ii] | 72.8$\pm$0.0 | 10.73$\pm$0.01 | 2.80$\pm$0.03 | 79.0$\pm$0.5 | 8.73$\pm$0.0.01 | 3.58$\pm$0.02
|  13CO(3–2) | 53.5$\pm$0.4 | 10.06$\pm$0.01 | 1.44$\pm$0.01 |  |  |
| [O i] 145 | 48.2$\pm$1.6 | 10.60$\pm$0.04 | 2.46$\pm$0.09 | 5.6$\pm$1.4 | 7.71$\pm$0.43 | 3.00$\pm$0.00c
| CO(6–5) | 165.2$\pm$0.3 | 10.30$\pm$0.01 | 2.13$\pm$0.01 |  |  |
| CO(7–6) | 116.2$\pm$0.7 | 10.24$\pm$0.01 | 1.95$\pm$0.01 |  |  |
| CO(11–10) | 28.7$\pm$0.3 | 10.29$\pm$0.81 | 1.61$\pm$0.02 |  |  |
(60,-60) | [C ii] | 39.6$\pm$1.1 | 10.08$\pm$0.04 | 2.20$\pm$0.00c | 15.2$\pm$1.4 | 7.65$\pm$0.09 | 2.40$\pm$0.24
|  13CO(3–2) | 47.7$\pm$0.4 | 10.22$\pm$0.01 | 1.09$\pm$0.01 |  |  |
| [O i] 145 | 10.6$\pm$0.9 | 10.41$\pm$0.05 | 1.13$\pm$0.12 | 5.8$\pm$2.8 | 8.87$\pm$1.04 | 3.00$\pm$0.00c
| CO(6–5) | 133.9$\pm$0.3 | 10.37$\pm$0.01 | 1.58$\pm$0.01 |  |  |
| CO(7–6) | 102.8$\pm$0.7 | 10.35$\pm$0.01 | 1.47$\pm$0.01 |  |  |
| CO(11–10) | 18.1$\pm$0.9 | 10.49$\pm$0.03 | 1.13$\pm$0.07 |  |  |
(-2,-78) | [C ii] | 126.5$\pm$1.1 | 10.17$\pm$0.01 | 2.87$\pm$0.03 |  |  |
|  13CO(3–2) | 55.5$\pm$0.5 | 10.11$\pm$0.01 | 1.98$\pm$0.02 |  |  |
| [O i] 145 | 50.2$\pm$1.1 | 10.77$\pm$0.02 | 1.64$\pm$0.04 |  |  |
| CO(6–5)a | 87.0$\pm$0.3 | 10.81$\pm$0.01 | 1.56$\pm$0.01 |  |  |
| CO(7–6)b | 68.6$\pm$0.8 | 10.68$\pm$0.01 | 1.51$\pm$0.02 |  |  |
| CO(11–10) | 14.2$\pm$1.2 | 10.70$\pm$0.07 | 1.49$\pm$0.15 |  |  |
(23,-60) | [C ii] | 112.4$\pm$2.7 | 10.27$\pm$0.03 | 2.49$\pm$0.09 | 50.9$\pm$1.8 | 8.19$\pm$0.05 | 2.42$\pm$0.17
|  13CO(3–2) | 49.3$\pm$0.9 | 10.00$\pm$0.01 | 1.32$\pm$0.01 |  |  |
| [O i] 145 | 40.7$\pm$3.8 | 10.33$\pm$0.02 | 1.61$\pm$0.08 | 14.7$\pm$3.9 | 9.50$\pm$0.38 | 4.30$\pm$0.65
| CO(6–5) | 120.1$\pm$0.0 | 10.23$\pm$0.01 | 1.97$\pm$0.01 |  |  |
| CO(7–6) | 90.6$\pm$0.8 | 10.21$\pm$0.01 | 1.75$\pm$0.02 |  |  |
| CO(11–10) | 12.6$\pm$0.8 | 10.19$\pm$0.04 | 1.12$\pm$0.10 |  |  |

Note: The CO(7–6) and CO(6–5) spectra have been convolved to a resolution of
15″.
The 13CO(3–2) and CO(11–10) spectra are at their respective original
resolutions of 18$\aas@@fstack{\prime\prime}$5 and 23″.
a Both blue- and red-shifted outflow wings masked out
b Blue-shifted outflow wing masked out
c Not fitted, i.e., kept constant

###  3.3 Column Density of the Atomic Oxygen

Figure 7: The smooth curve (green) shows the fit to the [O i]63 spectrum (red)
obtained by attenuating the scaled (by a factor 2 corresponding to the typical
intensity ratio between the two [O i] lines in temperature units) [O i] 145 µm
spectrum (filled grey histogram) by absorption due to foreground material.
Figure 8: The column density distribution of foreground PDR gas responsible
for the [O i] 63 absorption features in the NGC 2023 region plotted in color
and overlaid with the dense PDR gas traced by [O i] 145 emission. The contour
levels in K km s-1for the [O i] 145 are 3, 5 to 50 in steps of 5. This image
shows that the densest, most opaque area of the foreground cloud is in the
south and southwest, while the foreground extinction is much lower in the east
and in the north.  Table 2: Peak optical depth, corresponding velocity, and
the FWHM of the foreground absorbing gas derived from modeling the [O i] 63 µm
spectra with a two-slab model. Column density of oxygen in foreground cloud
($N_{\rm abs}$) is estimated using equation (1) and in the background PDR
($N_{\rm em}$) estimated using RADEX (Fig. 9).

Position | $\tau_{0}$ | $v_{\rm LSR}$ | $\Delta v$ |  $N_{\rm abs}$(O) |  $N_{\rm em}$(O)a |  $T_{\rm kin}^{b}$(C+) |  $N$(C+) |  $N$(H2)c
---|---|---|---|---|---|---|---|---
|  | km s-1 | km s-1 | (1018 cm-2 | (1018) cm-2) | K | (1018 cm-2) | (1021 cm-2)
(0, 0) | 2.2 | 10.4 | 2.0 | 0.88 | 3.4(2.8) | 74 | 3.4 | 7.8
(27,18) | 2.2 | 10.3 | 2.2 | 0.97 | 6.6(5.5) | 107 | 3.6 | 9.1
(40,-40) | 2.4 | 10.4 | 2.0 | 0.96 | 7.3(6.0) | 80 | 3.8 | 22.9
( 60,-60) | 2.9 | 10.4 | 1.8 | 1.0 | 2.0(1.6) | 60 | 2.4 | 23.9
( -2,-78) | 3.3 | 10.3 | 1.9 | 1.3 | 7.0(5.7) | 84 | 3.2 | 23.0
( 23,-60) | 3.5 | 10.1 | 1.9 | 1.3 | 9.1(7.2) | 81 | 3.8 | 29.0

a Estimated from integrated [O i] 145 intensity (Table 1) assuming $T_{\rm
kin}$ = 120 K and $n$(H2)= 2 105 & 106 cm-3.
The number in brackets corresponds to $N$(O) for $n$(H2)=106 cm-3.
b Estimated from Planck-corrected peak of the [C ii] spectrum as described in
Sec. 4
c Estimated from pixel-by-pixel grey-body fitting of 160, 250, 350 and 500 µm
dust continuum emission using hires, an improved
algorithm for the derivation of high-resolution surface densities from
multiwavelength far-infrared Herschel images (Men’shchikov, 2021).

The [O i] 63 µm emission is strongly self-absorbed over most of the reflection
nebula, whereas the [O i] 145 is always observed in emission with no evidence
for self-absorption. Even though [C ii] and [O i] 63 µm shows bright wings due
to photo-evaporation flows, the photo-evaporating gas is not dense enough to
be seen in [O i] 145 except in a few areas, see e.g., spectra at positions
40″, -40″, and 23″, -60″, where [O i] 145 shows a clear blue-shifted line
wing. In order to derive the distribution of atomic oxygen based on the 63 and
145 µm spectra, we fit the [O i] 63 spectra using a two-slab toy model: the
hot PDR shell (lying on either side of HD37903 along the line-of-sight) is
assumed to be the background layer, while the colder oxygen on the front side
is taken to be the foreground layer. While the gas in this layer, where the
PDR merges into the cold molecular cloud is largely molecular, most of the
oxygen is still atomic. We further assume that the emission from the
background layer is captured by the [O i] 145 emission. In this case the
colder foreground oxygen primarily attenuates the background spectral emission
at 63 µm by the factor $\exp\left[-\tau_{0}\exp\left[-4\ln
2(v-v_{0})^{2}/\Delta v^{2}\right]\right]$. Here $\tau_{0}$, $v_{0}$, and
$\Delta v$ denote the peak optical depth, the velocity at which $\tau_{0}$ of
the foreground cloud occurs and the full width at half maximum (FWHM) of the
absorption profile due to the foreground material. As pointed out earlier, the
[O i] 63 emission arising from a transition involving lower energy levels also
trace the blue- and red-shifted gas in photo-evaporation flows. The [O i] 145
emission, which is only excited in the warm and dense gas, does not show these
broad line wings. The background [O i] 63 emission is modeled by Gaussian
profiles generated from fits to the observed [O i] 145 spectrum scaled by a
factor of 2, which corresponds to the typical intensity ratio between the two
[O i] lines in temperature units and is equal to $I_{63\mu m}/I_{145\mu m}$ =
25 in energy units. Calculations for warm, medium density PDRs with $T_{\rm
kin}$= 100 K and $n$=105 cm-3 suggest a $I_{63\mu m}/I_{145\mu m}$ ratio
between 20 to 33 (Fig.10 Goldsmith, 2019). Assumption of a higher ratio of 33
leads to an increase in the derived oxygen column density by 10–20% of the
value obtained in this work. Figure 7 shows an example of the results of the
fits to the [O i] 63 using the two-slab toy model at the selected positions
and Table 2 presents the corresponding fitted parameters. At the selected
positions, the central velocity of the absorbing foreground cloud lies between
10.1–10.7 km s-1 which matches well with the velocity of the [O i] 145
emission spectra and have FWHM between 1.8–2.4 km s-1. This suggests that the
[O i] 63 absorption occurs in the colder foreground layers of the same PDR gas
that accounts for the emission of [O i] 145. We have performed the fitting of
the [O i] 63 profile over a larger number of positions in the map and obtained
the distribution of the foreground cold absorbing component (Fig. 8).

The fits using the [O i] 145 spectra as the model for the background component
do not reproduce the red- and blue-shifted wings observed in the [O i] 63
spectra, since the gas in the photo-evaporation flows causing these wings is
not dense enough to excite [O i] 145. We estimate the column density of the
absorbing foreground layer by assuming that all the oxygen atoms reside in the
lowest level of their ground state. The column density of all O atoms in the
foreground layer can be then estimated from the center opacity $\tau_{0}$ and
width (in km s-1) of the [O i] 63 line by

$N({\rm O^{0}})=2\times 10^{17}\tau_{0}\Delta v~{}~{}~{}~{}~{}{\rm cm^{-2}}$
(1)

The estimated column densities for the foreground gas ranges from 8 1017 to 2
$\times 10^{18}$ cm-2. We can estimate a lower limit to the extinction if we
assume that all the oxygen is atomic. With a typical gas abundance of
3.2$\times 10^{-4}$ for oxygen and assuming most of the hydrogen to be
molecular we estimate $N$(H2) of the foreground gas to range between 1.2–3.2
1021 cm-2, which translates to a minimum AV of 1.3. However, since one third
of the oxygen is tied up in CO in molecular gas a more realistic estimate for
the AV is 2 mag.

For the estimated column densities, based on non-LTE calculations using RADEX
(Van der Tak et al., 2007), the [O i] 145 line is optically thin over a large
range of temperatures (20–300 K) and volume densities (104–107 cm-3). This is
also consistent with our observations. Since the $T_{\rm mb}$ ratio for the
two [O i] lines depends on the physical conditions, the derived values of the
peak optical depth are subject to the ratio assumed in this work. However, the
fitted central velocity and line widths for the absorption profile are fairly
robust against the assumed scale factor. Furthermore, in the estimate of
$N$(O) we have assumed that all O atoms in the absorbing layer are in the
ground 3P2 level, which is reasonable because there is no hint of absorption
in [O i] 145.

Since the [O i] 63 µm is optically thick it cannot be used to estimate the
column density of oxygen atoms in the background layer. We thus use the
integrated intensities of the [O i] 145 µm spectral lines to estimate $N$(O).
For this we need an estimate of the temperature and density of the emitting
gas. Since the [O i] 145 emission seems to follow the CO(11–10) emission quite
closely, we can use the temperature and density estimated for CO(11–10), for
which Sandell et al. (2015) estimated a kinetic temperature of 120 K and
densities of 2 $10^{5}$–1 $10^{6}$ cm-3. Table 2 presents the column density
of atomic oxygen at the seven selected positions that is required to produce
the observed [O i] 145 emission (Table 1), estimated based on non-LTE
calculations using RADEX (Van der Tak et al., 2007) for these kinetic
temperatures and densities (Table 2). Figure 9 shows a plot of the intensity
of [O i] 145 line as a function of $N$(O) for $T_{\rm k}$ = 120 K and $n$(H2)
= (0.2–1) 106 cm-3.

###  3.4 The morphology of the reflection nebula and the PDR illuminated by HD
37903

Sandell et al. (2015) proposed that the PDR illuminated by HD 37903 formed an
expanding thin, hot shell which is roughly spheroidal with the south-eastern
part being almost spherical because the expansion is slowed down by the dense
surrounding molecular cloud and the shell being much more extended to the
north-west. They interpreted the double peaked [C ii] spectra seen on the
north-western side of the star as originating from the front and backside of
the PDR shell. They argued that the reason for not seeing clear double split
spectra on the southeastern side, was because the expansion was slowed down by
the dense molecular cloud. While their model appears reasonable, we now know
that it has some serious flaws. Even though we do see double peaked [C ii]
spectra north-west of HD 37903, see e.g. Fig. 5, [O i] 145 is single peaked
and centered between the two [C ii] peaks. This confirms that what we see is
self-absorption in [C ii], not two velocity components. The [C ii] lines are
quite broad, so it is probably true that they are broadened by emission from
both the front and the back side of the PDR, but they cannot be separated into
two distinct velocity components. The foreground extinction is definitely
lower in the north and northwest. We do see some blue-shifted [C ii] wings due
to photo-evaporation flows. This is consistent with the cloud emission being
more redshifted in both 13CO(3–2) and [C ii] indicating that it is behind the
PDR. It appears more likely that the PDR is fan-shaped rather than ’egg-
shaped’ and that it is most likely open in the north and north-west allowing
FUV photons to escape, hence creating a large reflection nebula dominating the
emission north and north-west of HD 37903.

##  4 Discussion

Analysis of the observed [O i] 63 and [O i] 145 spectra suggests that the
strong absorption features seen in the [O i] 63 profiles are caused by colder
atomic oxygen lying between the observer and the warm PDR emission excited by
HD 37903. The [C ii] spectra are also affected by self-absorption, but it is
far less severe and barely noticeable over most of the reflection nebula. At
the peak (27″, 18″) of the [C ii] emission (Fig. 6), the self-absorption
results in the shifting of the [C ii] spectra toward blue-shifted velocities
compared to [O i] 145 and the high-$J$ CO lines, which are unaffected by self-
absorption. Even the [C ii] spectrum on HD 37903 appears to be somewhat
affected by self-absorption.

The two-slab toy model, which we used in Section 3.3 to model the observed [O
i] 63 absorption features, shows that the column density of atomic oxygen in
the foreground absorbing gas is quite substantial, $\sim$ 1018 cm-2. This is
because most of the oxygen remains atomic way into the cloud, i.e. from the
hot PDR layer up to about an AV of 10, as already shown by Tielens &
Hollenbach (1985), while most of the carbon becomes molecular, i.e., gets tied
up in CO. The contribution to the oxygen column density from the surface of
the molecular cloud, which only gets ionized by the general interstellar
radiation that is a factor of 104 or 105 times weaker than the PDR illuminated
by HD37903, is negligible. The estimated [O i] column densities depend to some
extent on the “model" for the background emission, which have been constructed
using the [O i] 145 spectra, but as we can see in Table 2 the computed column
densities only depend weakly on the assumed gas densities.

The [O i] 63 spectral profiles seen in NGC 2023 are similar to those observed
in other strongly self-absorbed regions like S1 in RhoOph (Mookerjea et al.,
2021) and toward W3 (Goldsmith et al., 2021), where the absorption was
estimated to be due to colder foreground gas with $N$(O) of 3 1018 and 2–7
$10^{18}$ cm-2 respectively. The method for estimating the background for the
S1 PDR in Rho Oph was similar to what we used in this paper, while for W3,
which was not observed in [O i] 145 it was instead estimated by fitting the [O
i] 63 spectra for positions which showed almost no evidence for absorption. In
the [O i] 63 map in NGC 2023 there are only a few positions which show a
single emission peak without any trace of absorption, thus the approach to use
the [O i] 145 emission which has been measured is more robust.

Since we detect the [13C ii] $F$ = 2–1 line in essentially the entire area
mapped in [C ii], it is possible to check whether the [C ii] emission is
optically thick although the [C ii] map does not go deep enough to measure [C
ii] to [13C ii] ratio in individual positions. We have therefore divided up
the map into smaller regions of about 1 sq. arcmin, each containing about
50–70 spectra. When we did this we found that the ratio of intensities between
[C ii] to [13C ii] within errors is $\sim$ 30–40 after correcting for the
relative intensity of the $F$ = 2–1 line, 0.625. Considering the isotope ratio
12C/13C of 70 for the Orion region, this suggests the [C ii] optical depth to
be $\sim$ 1.5–2, consistent with the value of 2 derived by Sandell et al.
(2015). We obtain an estimate of the minimum kinetic temperature $T_{\rm
kin}$(C+) from the Planck-corrected peak of the [C ii] spectrum assuming a
beam filling factor of unity. The lower limit of $T_{\rm kin}$(C+) thus
estimated for the entire map ranges between 60–92 K. Assuming an optical depth
of 2 and $T_{\rm kin}$(C+), and the integrated [C ii] intensity to estimate
$N$(C+) using the Equation (26) from Goldsmith et al. (2012) modified for
optically thick emission as follows:

$N(C^{+})={\rm 2.91\times 10^{15}\left[1+0.5{\rm
e^{91.25/T_{kin}}}\left(1+\frac{A_{ul}}{C_{ul}}\right)\right]\frac{\tau}{1-e^{-\tau}}\int
T_{\rm mb}d\upsilon}$ (2)

where $A_{\rm ul}$ = 2.3$\times 10^{-6}$ s-1, $T_{\rm kin}$ is the gas kinetic
temperature, the collision rate is $C_{\rm ul}$ = $R_{\rm ul}n$ with $R_{\rm
ul}$ being the collision rate coefficient with H2 or H0, which depends on
$T_{\rm kin}$, and $n$ is the volume density of H. Since the critical density
of the [C ii] transition is $n_{\rm cr}$=3000 cm-3, and since it is likely
that most of the [C ii] detected could be at such densities along with some
emission arising from clumps with densities exceeding 105 cm-3, we assume a
density of 104 cm-3 to estimate the $N$(C+) density distribution of the
region. For this calculation we use the kinetic temperatures of C+ estimated
as above, a density of $n$=104 cm-3, and excitation from C+–H2 collisions,
with $R_{ul}$ = 3.8$\times 10^{-10}$ cm3 s-1. The value of $N$(C+) estimated
in the region show a small range of values between (1.8–4.4)$\times 10^{18}$
cm-2. We note that barring a few exceptions (Okada et al., 2015; Mookerjea et
al., 2021) for the estimate of $N$(C+) often a kinetic temperature of 100–120
K is assumed in contrast to using an estimate from the observed [C ii] peak
temperatures as we have done here. From equation (2) we find that for a
temperature of 60 K the $N$(C+) estimated will be approximately twice the
value that would be estimated for $T_{\rm kin}$=120 K. As seen at the selected
positions, throughout the entire region the ratio of column densities of C+
and O atoms (seen in emission) is consistent with the comparable to the solar
[O]/[C] abundance ratio of 3.5 within the uncertainties (Table 2). This is
also consistent with the structure of the PDR expected where the [O i] 145
emission arises from a dense hot region which also emits in [C ii], although
[C ii] emission additionally arises from the more diffuse and colder PDR gas.

Based on the far-infrared dust continuum maps at 160 to 500 µm we estimate the
molecular hydrogen column density to be around 7$\times 10^{21}$ cm-2, while
to the south, beyond the ridge and into the molecular cloud the values are
typically $\sim 2\times 10^{22}$ cm-2. The estimated H2 and C+ column
densities at the selected positions are consistent with the solar C/H
abundance ratio of 3$\times 10^{-4}$ (Table 2).

##  5 Summary & Conclusions

We have studied the geometry of the NGC 2023 reflection nebula using
spectrally resolved observations of the 63 and 145 µm transitions of [O i]
along with the [C ii] 158 µm transition which enable us to map the three-
dimensional distribution of the photodissociated gas vis-a-vis the location of
the illuminating source HD 37903. Combination of the velocity-resolved [O i]
145 spectra with the [C ii] spectra have allowed us to improve our
understanding of the morphology of the region. We conclude that the PDR is
fan-shaped rather than ’egg-shaped’ and that it is most likely open in the
north and north-west allowing FUV photons to escape, that creates a large
reflection nebula dominating the emission north and north-west of HD 37903.
The [C ii] emission detects PDR gas uniformly distributed with $N$(C+)$\sim
0.9$–$2.2\,10^{18}$ cm-2 suggest the contribution of both diffuse and dense
components, while the [O i] 145 emission detects dense PDR gas with $N$(O)
between (2–10) 1018 cm-2. Additionally, the self-absorbed profile of the [O i]
63 indicate the presence of cold low-excitation atomic oxygen with
$N$(O)$\sim$0.8–1.5 $10^{18}$ cm-2 associated with the nebula and lying in the
foreground between the illuminating source and the observer that is most
pronounced towards the southern and south-western part of the nebula. The
self-absorbed [O i] 63 profiles observed in NGC 2023 together with similar
results in other Galactic PDRs (Mookerjea et al., 2021; Goldsmith et al.,
2021) suggest the need to exercise caution when modelling galactic and
extragalactic [O i] 63 intensities, in particular if spectrally unresolved,
without accompanying observations of the optically thin [O i] 145 line.

## Acknowledgements

BM acknowledges the support of the Department of Atomic Energy, Government of
India, under Project Identification No. RTI 4002. This research has made use
of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI :
10.26093/cds/vizier). The original description of the VizieR service was
published in 2000, A&AS 143, 23"

## Data Availability

The new data observed with SOFIA that is presented here will be available from
the SOFIA data archive (https://irsa.ipac.caltech.edu/applications/sofia/).

## References

  * Field et al. (1998) Field D., Lemaire J. L., Pineau des Forêts G., et al. 1998, A&A, 333, 280
  * Goldsmith (2019) Goldsmith P. F., 2019, ApJ, 887, 54
  * Goldsmith et al. (2012) Goldsmith P. F., Langer W. D., Pineda J. L., Velusamy T., 2012, ApJS, 203, 13
  * Goldsmith et al. (2021) Goldsmith P. F., Langer W. D., Seo Y., Pineda J., Stutzki J., Guevara C., Aladro R., Justen M., 2021, ApJ, 916, 6
  * Graf et al. (2012) Graf U. U., et al., 2012, A&A, 542, L16
  * Guevara et al. (2020) Guevara C., Stutzki J., Ossenkopf-Okada V., et al. 2020, A&A, 636, A16
  * Hollenbach & Tielens (1997) Hollenbach D. J., Tielens A. G. G. M., 1997, ARA&A, 35, 179
  * Kaufman et al. (2006) Kaufman M. J., Wolfire M. G., Hollenbach D. J., 2006, ApJ, 644, 283
  * Lique et al. (2018) Lique F., Kłos J., Alexander M. H., Le Picard S. D., Dagdigian P. J., 2018, MNRAS, 474, 2313
  * Men’shchikov (2021) Men’shchikov A., 2021, A&A, 649, A89
  * Mookerjea et al. (2009) Mookerjea B., Sandell G., Jarrett T. H., 2009, A&A, 507, 1485
  * Mookerjea et al. (2018) Mookerjea B., Sandell G., Vacca W., et al. 2018, A&A, 616, A31
  * Mookerjea et al. (2021) Mookerjea B., Sandell G., Veena V. S., Güsten R., Riquelme D., Wiesemeyer H., Wyrowski F., Mertens M., 2021, A&A, 648, A40
  * Okada et al. (2015) Okada Y., Requena-Torres M. A., Güsten R., Stutzki J., Wiesemeyer H., Pütz P., Ricken O., 2015, A&A, 580, A54
  * Ossenkopf et al. (2013) Ossenkopf V., Röllig M., Neufeld D. A., Pilleri P., Lis D. C., Fuente A., van der Tak F. F. S., Bergin E., 2013, A&A, 550, A57
  * Risacher et al. (2016) Risacher C., et al., 2016, A&A, 595, A34
  * Sandell et al. (2015) Sandell G., Mookerjea B., Güsten R., Requena-Torres M. A., Riquelme D., Okada Y., 2015, A&A, 578, A41
  * Tielens & Hollenbach (1985) Tielens A. G. G. M., Hollenbach D., 1985, ApJ, 291, 722
  * Van der Tak et al. (2007) Van der Tak F. F. S., Black J. H., Schöier F. L., Jansen D. J., van Dishoeck E. F., 2007, A&A, 468, 627
  * Wyrowski et al. (2000) Wyrowski F., Walmsley C. M., Goss W. M., Tielens A. G. G. M., 2000, ApJ, 543, 245

##  Appendix A Result of non-LTE calculations using RADEX

Figure 9 shows the result of the non-LTE radiative transfer model calculations
to estimate $N$(O) assuming $T_{\rm kin}$= 120 K and $n$(H2) = 2 105 and 106
cm-3 performed using RADEX (Van der Tak et al., 2007).

Figure 9: The peak [O i] 145 intensities for $\Delta v=2$ km s-1, $T_{\rm
kin}$=120 K and $n$(H2) = 2 105 (black) and 106 cm-3(red) calculated using the
non-LTE radiative transfer models RADEX (Van der Tak et al., 2007) and using
collision rate coefficients from (Lique et al., 2018).


# Neutron penumbral imaging simulation and reconstruction for Inertial
Confinement Fusion Experiments

Xian-You Wang<EMAIL_ADDRESS>Department of Physics, Chongqing
University, Chongqing 400044, P.R. China  Institute of High Energy Physics,
Chinese Academy of Sciences, Beijing 100049, P.R. China  Zhen-yun Fang
Department of Physics, Chongqing University, Chongqing 400044, P.R. China
Yun-qing Tang<EMAIL_ADDRESS>Department of Physics, Chongqing University,
Chongqing 400044, P.R. China  Institute of theoretical Physics, Chinese
Academy of Sciences, Beijing 100049,P.R. China  Zhi-Cheng Tang  Hong Xiao
Ming Xu  Institute of High Energy Physics, Chinese Academy of Sciences,
Beijing 100049, P.R. China

###### Abstract

Neutron penumbral imaging technique has been successfully used as the
diagnosis method in Inertial Confined Fusion. To help the design of the
imaging systems in the future in CHINA. We construct the Monte carlo imaging
system by Geant4. Use the point spread function from the simulation and decode
algorithm (Lucy-Rechardson algorithm) we got the recovery image.

Inertial confinement fusion (ICF) has attracted much attention from all over
the world. The National Ignition Facility (NIF) programs in America and Laser
MegaJoule (LMJ) in France aim at reaching the inertial confinement fusion of a
deuterium-tritium (DT) filled cryogenic target. China’s ICF program are also
in progress, this paper’s aim is just study the neutron diagnosis system which
will be used in the future ICF program in China.

The process of ICF can be divided into four steps as Fig 1 : First step
generate radiation to drive implosion. There are two approaches to drive
implosion, one is laser direct driven, another is laser energy injected inside
the cavity generates X-rays that heat the DT fuel. Under X-ray or laser
compression, the target density as well as the mean ions’ temperature
increases. Fusion reactions will occur within a hot spot generating 14 MeV
neutrons and alpha particlesICFSim ; ICFSim1 . If the density is sufficient to
trap the energetic alphas or other radiation, these energy will be deposited
on the surrounding cold DT medium. Then the chained thermonuclear reactions
will lead to a combustion in the compressed target. For a well compressed
target, combustion will last long enough to generate more energy than the
injected energy. If ever, drive symmetry is poor or pulse shaping is
imprecise, burning will be low. So we can use the neutrons escaped from the
compressed cores of target to study the compressed process to make sure the
drive be symmetrical.

There are several neutron imaging techniques for inertial confinement fusion
(ICF) experiment to detect the size, shape, and uniformity of the compressed
cores. The pinhole imaging is the easiest way, but the low neutron yield has
an effect on the sensitivity. In NIF experiment a pinhole array is designed to
increase the sensitivitypinhole1 ; pinhole2 ; pinhole3 ; pinhole4 ; pinhole5 .

The another technique is the neutron penumbral imaging (NPI), which have been
widely studied by LMJpenum1 ; penum2 ; penum3 ; penum4 ; penum5 ; penum6  and
other groups penum_other1 ; penum_other2 ; penum_other3 ; penum_other4 ;
penum_other5 ; penum_other6 . Compared with other neutron imaging techniques,
the sensitivity of neutron penumbral imaging is higher. Annular imaging is
nearly the same as penumbral imaging, the only difference of the imaging
process is aperture.

Figure 1: ICF physics process (Radiation compressed, Implosion, Hot spot
creation, Combustion)

In this paper we just based on the simulation model in the referencecqu1 ;
cqu2 ; cqu3 ; cqu4  to design a software package NGPINEX which can be used to
simulate the neutron imaging process and recovery the coded image. The imaging
process simulation have been described in Geant4, and then the recover
algorithm we used Lucy-Richardson methodLR1 ; LR2 , which have been
successfully applied on decoding x-ray ring code imagexrayLR . The whole
program is written in C++, which can be run under the linux system with multi-
core server.

In the program we simulated neutron penumbral imaging and Annular imaging, the
pinhole imaging block is also included in the program. Based on the NIF and
LMJ experiment we design an imaging system aimming to achieve $20\mu m$
resolution. At last we use the Lucy-Richardson method successfully got the
recovered image with $20\mu m$ resolution

##  I The imaging system simulation

The structure of neutron penumbral imaging system is shown in Fig.2. The
source , aperture and detector composed the whole imaging system. The fast
neutrons emitting from a burning target are scattered by the aperture
materials, and form an penumbral image on the detector arrays. In our
simulation the coded aperture is placed at $55mm$ from the source, the
detector is placed at $7m$ from the aperture. The aperture is designed as a
$50mm(h)$ thick $10mm$ diameter(tungsten) cylinder to achieve an effective
absorption of the $14MeV$ neutrons in the opaque parts of the aperture. The
pattern of the aperture is a disk (penumbra) which consists in a biconical
hole, the diameters $(D_{1},D_{2},D_{3})$ of the biconical hole are
$0.3mm,0.38mm,0.535mm$. An enlarged coded image is projected on the neutron
detector which is placed on the target bay floor, at $7895mm$ from the
aperture. The image obtained from the detector is a convolution of the source
spatial distribution with the aperture transfer function.

Figure 2: Penumbral imaging

In the Fig.2 only the red annular of the coded image contains the effective
information of the source shape changing. To achieve the high SNR we can place
a plug inside the hole of an aperture of the same dimensions as the
penumbra(Fig.2) to define a annular, which will be study in another paper.
This pattern can shield a lot of background but also face the alignment
tolerance problem and not enough neutron can be detected.

The neutron detector is composed of $255\times 255$ array of scintillating
fibers Fig.3. Each scintillating fiber can be realized as a sub-detector, the
structure was designed as the left figure of Fig.3. The core part is a tubes
of $460\mu m$ diameter and the material is polystyrene BCF-10 (molecular
formula C6H5CH=CH2, $1.05g/cm^{3}$). 14 MeV neutrons mainly interact by
elastic scattering on hydrogen nuclei inside the fibers. As they lose kinetic
energy, the recoil protons will produce light, a part of the light is guided
through these fibers. Outside the core part is a layer of the clad material,
which is used to avoid the scintillating light interaction in the array of
scintillating fibers. The material of the envelope is Polyvinyl
alcohol(molecular formula CH2=CHOH, $1.26g/cm^{3}$) and the thickness is
$20\mu m$. An optical relay then casts the scintillating light from the end of
the array to an image intensifier tube (IIT) which gates the 14 MeV neutrons
at their arrival time on the detector. The image is then reduced with a
fibered optical taper onto a CCD.

Figure 3: Fiber array

Overall spatial resolution of the system for both the annular and the
penumbral apertures is the same, which is determined by the point spread
functions (PSF) of the coded aperture and detector resolution ($\Delta
s_{detector}$) scaled down the source plane as equation 1. In the paraxial
approximation, the broadening of the apertures point spread function is
identical, determined by the FOV and the aperture to target distance
($L_{0}$).

$\Delta s=\sqrt{\frac{ln(2)\times FOV}{2\times L_{0}\times\mu}+\Delta
s_{detector}^{2}\times\left(\frac{L_{0}}{L_{1}}\right)^{2}}$ (1)

where $\mu$ is the attenuation of the neutrons in tungsten and $L_{1}$ is the
detector to the target distance. With this design, the overall spatial
resolution of the system is about 20 $\mu m$ in source plane.

In the GEANT4 simulation we used the physics list QGSP_BERT_HP which contains
the high precision neutron package (NeutronHP) to simulate the process of
neutrons’ transportation with the energy below 20 MeV down to thermal
energies. The physical process neutron capture, fission, elastic and inelastic
scattering (including absorption) are treated by referring to the ENDF-B VI
cross-section data. Besides, in the refMCNP  has demonstrated the validity of
GEANT4 calculations in neutron generation and transportation which is the same
as MCNP. Moreover, object-oriented programming and highly sophisticated
processing of both electromagnetic interactions and ionization processes
enhance GEANT4 capabilities for the study of mixed fields and complicated
geometries.

In the ICF experiment the neutrons generated by hot spot spread in whole
space. The idealized aperture only allow the neutrons in the solid angle as
Fig.2 to deposit the energy on the detector. But actually because the neutrons
have the high ability to penetrate the start part of the aperture, and also
the neutrons will scatter with the wall of the aperture, these neutrons which
hit the detector will become the background of the imaging. To simulate the
complete imaging process and also save the calculation time we ask the program
only generate random distribution neutron in the cone with $R_{i}$ radius as
Fig.4. As equation 2 $R_{i}$ is determined by $R_{0}$ which is the width of
the detector project on the front face of the aperture.

With this design there are about one half of the events be scattered by the
penumbral aperture.To get the high statistic sample we simulated $1\times
10^{10}$ effective events of penumbral imaging.

Figure 4: The neutron radiated region in the simulation

$\rm R_{i}=\left\\{\begin{array}[]{lll}D_{0}/2,&R_{0}>D_{0}/2\\\
1.3R_{0},&D_{3}/2\leq R_{0}\leq D_{0}/2\\\
D_{3}/2,&R_{0}<D_{3}/2\end{array}\right.$ (2)

##  II Reconstruction of the coded image

At present, there are a lot of reconstruction methods have been developed.
Among them the simplest and quickest method is Wiener filter method, which is
a linear reconstruction method. As we cannot get prior knowledge of the
signal-noise ratio, wiener filter method cannot get good results. Therefore we
can use this method as the reference or use the result as the input of some
nonlinear reconstruction method. In these years there are some nonlinear
reconstruct methods have been introduced to reconstruct penumbral image such
as genetic algorithm and heuristic algorithm proposed by Chen et
al.penum_other1 , and the classical molecular dynamics reconstruction method
proposed by Liu et al.penum_other2 ; penum_other3 ,

In the experiment L. Disdier et al. had successfully used the filtered
autocorrelation techniquefilauto1  to get $20\mu m$ resolution reconstruction
image from the penumbral imaging system with $22.3\mu m$ theoritical
resolution. Filtered autocorrelation technique is different from other method
which don’t need to rely on the simulation PSF in the reconstruction.

Another nonlinear method Lucy-Rechardson algorithmLR1 ; LR2  has been
successfully used on X-ray ring coded imaging in experimentxrayLR . And this
algorithm has been most widely used for restoring images with noise from
counting statistics and the data approximate Poisson statistics. In this
program we first try to use the Lucy-Rechardson method to get reconstructed
image. The process of the imaging can be described as a spatial intensity
distribution $g(x,y)$ convolved by a neutron source spatial intensity
distribution $f(x,y)$ with a point spread function (PSF) $h(x,y)$. Add the
noise $n(x,y)$ the detector image can be expressed as

$g(x,y)=f(x,y)\otimes h(x,y)+n(x,y)$ (3)

Lucy-Rechardson algorithm is an iterative method. Based on the referenceLR1 ;
LR2 , the decode formulation can be expressed as 4. Where PSF(h) and image on
detector(g) come from the simulation, the reconstructed image(f) can be
extracted after k times iterate.

$f_{k+1}=f_{k}(h\odot\frac{g}{h\otimes f_{k}})$ (4)

The PSF for penumbral and annular imaging are presented in Fig. 5
respectively. In Fig. 6 we presented the raw images with a ‘F’ shape neutron
source and the unfolded results with Lucy-Rechardson method. The width of the
‘F’ shape source is $20\mu m$.

Figure 5: Distribution of point spread function(PSF).

Figure 6: ‘F’ source imaging for penumbral aperture(left), unfolded result
(right).

From the Fig. 6 we can see with this design the imaging system can achieve the
aim to reconstruct the neutron source. With Lucy-Rechardson method the
unfolded image of Penumbral aperture is clear. To get the real resolution of
the system we simulate two point source with distance $20\mu m$ imaging on the
detector. In the Fig. 7 present the raw image and unfolded image.

Figure 7: Two point source imaging for penumbral aperture(left), unfolded
result(right).

For the annular imaging reconstruction and detail of the unfolding method we
will present in further study.

##  III Conclusion

We have successfully construct one software NGPINEX which integrated the
imaging system simulation and the reconstruction method. The simulation and
reconstruction result show we can use this imaging system with these
parameters to get $20\mu m$ resolution image of DT implore. To attend the aim
of $5\mu m$ resolution in the future we can use this software to optimize the
system parameter.

Acknowledgments: The authors would like to thank Ming Jiang and Bing-quan Hu
for the help on the simulation model construction. This work was supported in
part by the Fundamental Research Funds for the Central Universities under
Grant No.CDJXS1102209 and the Program for New Century Excellent Talents in
Univresity under Grant No. NCET-10-0882, and by Natural Science Foundation of
China under Grant No.10805082 and No.11075225.

## References

  * (1) PA Bradley, DC Wilson, FJ Swenson, and GL Morgan.  Icf ignition capsule neutron, gamma ray, and high energy x-ray images.  Rev. Sci. Instrum, 74(3, Part 2, SI):1824–1827, MAR 2003.
  * (2) M. J. Moran, S. W. Haan, S. P. Hatchett, N. Izumi, J. A. Koch, R. A. Lerche, and T. W. Phillips.  Energy-resolved neutron imaging for inertial confinement fusion.  Review of Scientific Instruments, 74(3):1701–1704, 2003.
  * (3) Olivier Landoas, L. Disdier, T. Caillaud, I. Thfoin, Jean-Luc Bourgade, B. Rosse, T. Craig Sangster, and Vladimir Yu. Glebov.  High resolution neutron imaging for inertial confinement fusion experiments.  Abstr. Am. Chem. Soc, 238, AUG 16 2009.
  * (4) G. P. Grim, P. A. Bradley, R. D. Day, D. D. Clark, V. E. Fatherley, J. P. Finch, F. P. Garcia, S. A. Jaramillo, A. J. Montoya, G. L. Morgan, J. A. Oertel, T. A. Ortiz, J. R. Payton, P. Pazuchanics, D. W. Schmidt, A. C. Valdez, C. H. Wilde, M. D. Wilke, and D. C. Wilson.  Neutron imaging development for megajoule scale inertial confinement fusion experiments.  In 5th International Conference on Inertial Fusion Sciences and Applications (IFSA2007), volume 112, 2008.
  * (5) Mark D. Wilke, Steven H. Batha, Paul A. Bradley, Robert D. Day, David D. Clark, Valerie E. Fatherley, Joshua P. Finch, Robert A. Gallegos, Felix P. Garcia, Gary P. Grim, Steven A. Jaramillo, Andrew J. Montoya, Michael J. Moran, George L. Morgan, John A. Oertel, Thomas A. Ortiz, Jeremy R. Payton, Peter Pazuchanics, Derek W. Schmidt, Adelaida C. Valdez, Carl H. Wilde, and Doug C. Wilson.  The national ignition facility neutron imaging system.  Rev. Sci. Instrum, 79(10), OCT 2008.
  * (6) D. C. Wilson, G. P. Grim, I. L. Tregillis, M. D. Wilke, M. V. Patel, S. M. Sepke, G. L. Morgan, R. Hatarik, E. N. Loomis, C. H. Wilde, J. A. Oertel, V. E. Fatherley, D. D. Clark, D. N. Fittinghoff, D. E. Bower, M. J. Schmitt, M. M. Marinak, D. H. Munro, F. E. Merrill, M. J. Moran, T. S. F. Wang, C. R. Danly, R. A. Hilko, S. H. Batha, M. Frank, and R. Buckles.  Modeling the national ignition facility neutron imaging system.  Rev. Sci. Instrum, 81(10), OCT 2010.
  * (7) E. N. Loomis, G. P. Grim, C. Wilde, D. C. Wilson, G. Morgan, M. Wilke, I. Tregillis, F. Merrill, D. Clark, J. Finch, D. Fittinghoff, and D. Bower.  Progress toward the development and testing of source reconstruction methods for nif neutron imaging.  Rev. Sci. Instrum, 81(10), OCT 2010.
  * (8) L Disdier, A Rouyer, DC Wilson, A Fedotoff, C Stoeckl, JL Bourgade, VY Glebov, JP Garconnet, and W Seka.  High-resolution neutron imaging of laser imploded dt targets.  Nucl. Instrum. Methods Phys. Res. A, 489(1-3):496–502, AUG 21 2002\.
  * (9) L Disdier, A Rouyer, A Fedotoff, JL Bourgade, FJ Marshall, VY Glebov, and C Stoeckl.  Neutron imaging of icf target plasmas.  Rev. Sci. Instrum, 74(3, Part 2, SI):1832–1836, MAR 2003.
  * (10) L Disdier, RA Lerche, JL Bourgade, and VY Glebov.  Capillary detector with deuterated scintillator for inertial confinement fusion neutron images.  Rev. Sci. Instrum, 75(6):2134–2139, JUN 2004.
  * (11) T. Caillaud, O. Landoas, M. Briat, S. Kime, B. Rosse, I. Thfoin, J. L. Bourgade, L. Disdier, V. Yu. Glebov, F. J. Marshall, and T. C. Sangster.  Development of the large neutron imaging system for inertial confinement fusion experiments.  Rev. Sci. Instrum, 83(3), MAR 2012.
  * (12) L Disdier, A Rouyer, I Lantuejoul, O Landoas, JL Bourgade, TC Sangster, VY Glebov, and RA Lerche.  Inertial confinement fusion neutron images.  Phy. Plasmas, 13(5), MAY 2006.
  * (13) M. C. Ghilea, T. C. Sangster, D. D. Meyerhofer, R. A. Lerche, and L. Disdier.  Aperture tolerances for neutron-imaging systems in inertial confinement fusion.  Rev. Sci. Instrum, 79(2, Part 1), FEB 2008.
  * (14) YW Chen, Z Nakao, K Arakaki, and S Tamura.  Distortion-free reconstruction of neutron penumbral images by a genetic algorithm.  IEEE Trans. Nucl. Sci., 44(3, Part 1):905–910, JUN 1997.
  * (15) YW Chen, Z Nakao, and I Nakamura.  A quantitative evaluation of neutron penumbral imaging with a toroidal-segment aperture.  Ieice Trans. Electron., E80C(2):346–351, FEB 1997.
  * (16) YW Chen, Z Nakao, and K Arakaki.  Genetic algorithms applied to neutron penumbral imaging.  Opt. Rev., 4(1B):209–215, JAN-FEB 1997.
  * (17) Weixin Qian, Shuangxi Qi, Wanli Wang, Jinming Cheng, and Dongbing Liu.  A new nonlinear reconstruction method based on total variation regularization of neutron penumbral imaging.  Rev. Sci. Instrum, 82(9), SEP 2011.
  * (18) Liu Dong-Jian, An Zhu, Tang Chang-Huan, and Zou Lian.  Studies of the aperture design and its pointing precision in fast neutron penumbral imaging.  Chin. Phys. C, 32:235–238, OCT 2008.
  * (19) R. Azuma, S. Nozaki, S. Fujioka, Y. W. Chen, and Y. Namihira.  Heuristic optimization in penumbral image for high resolution reconstructed image.  Rev. Sci. Instrum, 81(10), OCT 2010.
  * (20) Shi-biao Tang, Qing-li Ma, and Ji-Wei Zou.  Numericai simulation on psf of pinhole imaging system for fast neutron.  Nucl. Phys. Rev., 25(2):139–144, JUN 2008.
  * (21) Qingli Ma, Shibiao Tang, Bingquan Hu, Min Jiang, Daimin Zeng, and Zhenyun Fang.  Simulationon point spread function of a pinhole imging system for 6mev neutrons.  Nucl. Tech., 31(1):23–26, JAN 2008.
  * (22) Min JIANG, Zhenyun Fang, Bingquan Hu, Qingli Ma, Wensuo Chen, and Shibiao Tang.  Point spread function of neutron pinhole imaging in high gamma ray background.  High Power Laser and Particle Beams, 20(3):517–520, Mar 2008.
  * (23) Min JIANG, Bingquan Hu, Qingli Ma, Shi-biao Tang, Zhenyun Fang, and Daiming Zeng.  Study on point spread function of neutron pinhole imaging system simulation.  High Power Laser and Particle Beams, 37(2):310–313, Feb 2008.
  * (24) Richards. WH.  Bayesian-based iterative method of image restoration.  J. Opt. Soc. Am. A, 62(1):55–, 1972.
  * (25) L. B. Lucy.  An iterative technique for the rectification of observed distributions.  Astronom. J, 79:745, JUN 1974.
  * (26) Zongqing Zhao, Yongkun Ding, Jianjun Dong, Yidan Haol, Shunchao Wui, Leifeng Caol, and Yikang Pu.  Richardson-lucy method for decoding x-ray ring code image.  Plasma Phys. Control. Fusion, 49(8):1145–1150, AUG 2007.
  * (27) R Lemrani, M Robinson, VA Kudryavtsev, M De Jesus, G Gerbier, and NJC Spooner.  Low-energy neutron propagation in mcnpx and geant4.  Nucl. Instrum. Methods Phys. Res. A, 560(2):454–459, MAY 10 2006\.
  * (28) A Rouyer.  A new simple method for decoding penumbra image: The filtered autocorrelation.  Rev. Sci. Instrum, 74(3, Part 1):1234–1239, MAR 2003.


#  Describing Sen’s Transitivity Condition by Use of Inequality and Equation

Fujun Hou
School of Management and Economics
Beijing Institute of Technology
Beijing, China, 100081
Tel.: +86 10 6891 8960; Fax: +86 10 6891 2483; email<EMAIL_ADDRESS>
###### Abstract

In social choice theory, Sen’s value restriction condition is a sufficiency
condition restricted to individuals’ ordinal preferences so as to obtain a
transitive social preference under the majority decision rule. In this
article, Sen’s transitivity condition is described by use of inequality and
equation. First, for a triple of alternatives, an individual’s preference is
represented by a preference map, whose entries are sets containing the ranking
position or positions derived from the individual’s preference over that
triple of those alternatives. Second, by using the union operation of sets and
the cardinality concept, Sen’s transitivity condition is described by
inequalities. Finally, by using the membership function of sets, Sen’s
transitivity condition is further described by equations.

Keywords: social choice theory, social preference transitivity, Sen’s value
restriction condition, preference map

##  1 Introduction

If the number of alternatives is more than two, Arrow proved that, it is
impossible for the majority decision rule to yield a transitive social
ordering under the requirement that Arrow’s conditions 2-5 are satisfied
(Arrow, 1950, 1951). This famous ”Impossibility Theorem” has a profound impact
on the development of modern social choice theory. To make the majority
decision rule yield a transitive social preference, researchers attempted to
impose some constraints on the individuals’ preferences patterns. Thus there
have appeared many transitivity conditions among which Sen’s value restriction
condition (Sen, 1966) is a representable one since it covers several well-
known sufficiency conditions that have preceded it such as the single-peaked
preferences pattern (Arrow, 1951; Black, 1958), the single-caved preferences
pattern (Inada, 1964), the dichotomous preferences pattern (Inada, 1964), the
two-group separated preferences pattern (Inada, 1964), and the Latin-square-
lessness preferences pattern (Ward, 1965).

Sen’s condition does its work by requiring that in any triple of alternatives
at least one alternative not have the value ’best’, ’worst’, or ’medium’.
Therefore, Sen’s condition is a qualitative one. This article is concerned
with how to describe Sen’s condition in a quantitative manner. In particular,
we will characterize it by use of inequality and equation. Quantitative
description can be a great help for examining Sen’s condition with a computer.

##  2 Preference map and membership preference map

When it comes to the social preference transitivity, Arrow and Sen (among
others) are concerned primarily with a scenario that the individuals’
preferences over a finite alternative set are weak orderings (Arrow, 1951;
Sen, 1970). A weak ordering on a finite set can be and will be here
represented by a preference map. This is crucial for the discussion in this
article.

The preference map representation method assumes that the alternatives in a
tie are tied together occupying common positions, and these positions are
consecutive positive integers. Denote by $A=\\{x_{1},x_{2},\ldots,x_{m}\\}$
the alternative set with $1<m<+\infty$ and by $N=\\{1,2,\ldots,n\\}$ the
individual group with $1<n<+\infty$. The preference map is defined as follows
(Hou, 2015a, 2015b; Hou & Triantaphyllou, 2019).

Definition 1 A sequence $(PM_{i})_{m\times 1}$ is called a preference map (PM)
with respect to a weak order relation $\preceq$ on the alternative set
$A=\\{x_{1},x_{2},\ldots,x_{m}\\}$ such that

$PM_{i}=\\{|\xi_{i}|+1,|\xi_{i}|+2,\ldots,|\xi_{i}|+|\eta_{i}|\\},$ $None$

where the notation $|\centerdot|$ represents the cardinality of a set;
$\xi_{i}$ is the predominance set of alternative $x_{i}$, i.e.,
$\xi_{i}=\\{x_{k}\mid x_{k}\in A,x_{i}\prec x_{k}\\}$; and $\eta_{i}$ is the
indifference set of alternative $x_{i}$, i.e., $\eta_{i}=\\{x_{k}\mid x_{k}\in
A,x_{i}\sim x_{k}\\}$.

For illustrative purpose, we consider an example.

Example 1 Suppose that three voters provide their preferences over an
alternative triple $\\{x,y,z\\}$ as follows.

$\begin{array}[]{ll}\text{Voter 1}:&x\succ y\succ z,\\\ \text{Voter 2}:&x\sim
y\succ z,\\\ \text{Voter 3}:&x\sim y\sim z.\end{array}$

According to Def.1, the voters’ preferences can be represented by the
following preference maps

$\begin{array}[]{c}{\hskip 8.5359pt\begin{array}[]{cccc}\text{Voter 1}:&\hskip
5.69046pt\text{Voter 2}:&\hskip 11.38092pt\text{Voter 3}:\end{array}}\\\
{\begin{array}[]{cc}{\begin{array}[]{c}x\\\ y\\\
z\end{array}}&{\left[\begin{array}[]{c}\\{1\\}\\\ \\{2\\}\\\
\\{3\\}\end{array}\right],\left[\begin{array}[]{c}\\{1,2\\}\\\ \\{1,2\\}\\\
\\{3\\}\end{array}\right],\left[\begin{array}[]{c}\\{1,2,3\\}\\\
\\{1,2,3\\}\\\ \\{1,2,3\\}\end{array}\right],}\end{array}}\end{array}$

from which one can see the alternatives’ ranking positions corresponding to
the voters’ preferences. For instance, voter 1 would like to rank alternative
$x$ at position 1, while voter 2 would like to rank alternative $x$ at
position 1 or position 2, and positions 1-3 are all agreed by voter 3 with
regard to alternative $x$’s rank order. According to Sen (1966), by the way,
voter 3’s preference in Example 1 is a sort of unconcerned preference over the
triple $\\{x,y,z\\}$. A unconcerned individual is a person who is indifferent
between all the alternatives (Sen, 1966), and accordingly, the preference of a
unconcerned person is called unconcerned preference. In contrast, a concerned
individual is a person who is not indifferent between all the alternatives
(Sen, 1966).

Because the preference map entries are sets rather than numbers, a preference
map will have a one-to-one membership correspondence which is none other than
a 0-1 matrix as defined below.

Definition 2 For a preference map $(PM_{i})_{m\times 1}$, its membership
preference map (mPM) is a $m\times m$ 0-1 matrix $(mPM_{i,j})_{m\times m}$
whose elements are defined as

${mPM_{i,j}}=\begin{cases}1,&\text{if $j\in PM_{i}$},\\\
0,&\text{otherwise}.\end{cases}$ $None$

When the previous data are used, then the corresponding mPMs of the PMs in
Example 1 are as follows:

$\begin{array}[]{c}{\hskip 8.5359pt\begin{array}[]{cccc}\text{Voter 1}:&\hskip
5.69046pt\text{Voter 2}:&\hskip 11.38092pt\text{Voter 3}:\end{array}}\\\
\begin{array}[]{c}{\hskip 22.76228pt\begin{array}[]{cccc}\\{1\\}\hskip
2.27626pt\\{2\\}\hskip 2.27626pt\\{3\\}&\hskip 2.84544pt\\{1\\}\hskip
2.27626pt\\{2\\}\hskip 2.27626pt\\{3\\}&\hskip 2.84544pt\\{1\\}\hskip
2.27626pt\\{2\\}\hskip 2.27626pt\\{3\\}\end{array}}\\\
{\begin{array}[]{cc}{\begin{array}[]{c}x\\\ y\\\
z\end{array}}&{\left[\begin{array}[]{ccc}1&0&0\\\ 0&1&0\\\
0&0&1\end{array}\right],\left[\begin{array}[]{ccc}1&1&0\\\ 1&1&0\\\
0&0&1\end{array}\right],\left[\begin{array}[]{ccc}1&1&1\\\ 1&1&1\\\
1&1&1\end{array}\right],}\end{array}}\end{array}\end{array}$

where the left column indicates the alternatives, the first row indicates the
voters, the second row indicates the ranking positions, and a 0-1 value
indicates whether a certain alternative can or cannot be ranked at a certain
position by a certain voter.

##  3 Sen’s condition in inequalities

We start in this section by introducing Sen’s transitivity condition and some
underlying concepts. Thereafter, we describe Sen’s transitivity condition in
inequalities and illustrate an example.

Theorem 1 (Sen’s transitivity condition) The method of majority decision is a
social welfare function satisfying Arrow’s Conditions 2-5 and the consistency
condition for any number of alternatives, provided the preferences of
concerned individuals over every triple of alternatives is Value-Restricted,
and the number of concerned individuals for every triple is odd.

Two concepts that are used in Theorem 1 are defined by Sen (1966) as follows.

  * $\diamond$

Value of alternative: The value of an alternative in a triple for an
individual preference ordering is its characteristic of being ”best,” ”worst,”
or ”medium.”

  * $\diamond$

Value-Restricted preference pattern: A set of individual preferences over a
triple of alternatives such that there exist one alternative and one value
with the characteristic that the alternative never has that value in any
individual’s preference ordering, is called a Value-Restricted Preference
pattern over that triple for those individuals.

Before presenting our result, we discuss the key aspects of our analysis.

An alternative in a triple being a Value-Restricted one means that this
alternative is excluded from being assigned at least one position from
$\\{1,2,3\\}$. From Definition 1 we know that an entry of a PM is a set
containing an alternative’s possible ranking position or positions. When the
individuals are concerned ones over a triple of alternatives and their
preferences over this triple are represented by PMs, if an alternative in the
triple is a Value-Restricted one, then the union set of the concerned
individuals’ PMs’ entries corresponding to this alternative will not include
at least one element from $\\{1,2,3\\}$, i.e., in this case the cardinality
number of the specified union set is less than 3.

The above analysis indicates that, Sen’s condition can be stated in an
alternative way. This is described explicitly next as Theorem 1’.

Theorem 1’ Assume that the alternative set is
$A=\\{x_{1},x_{2},\ldots,x_{m}\\}$ where $2<m<+\infty$, the individual set is
$N=\\{1,2,\ldots,n\\}$ where $2<n<+\infty$, and the individuals’ preferences
over the alternative set are weak orderings. Denote by
$(x_{t_{1}},x_{t_{2}},x_{t_{3}})$ the triples of alternatives, where
$t=1,2,\ldots,\binom{m}{3}$. Denote by $V_{t}$ the maximal set including
concerned individuals who are specified in terms of their preferences over
that triple $(x_{t_{1}},x_{t_{2}},x_{t_{3}})$, and represent their preferences
over that triple by PMs. Specifically, for triple
$(x_{t_{1}},x_{t_{2}},x_{t_{3}})$ the concerned individuals’ PMs are denoted
by

$\begin{array}[]{c}{\hskip 8.5359pt\begin{array}[]{cccc}\text{Concerned
individual k's PM with respect to triple
}(x_{t_{1}},x_{t_{2}},x_{t_{3}}):\end{array}}\\\
{\begin{array}[]{cc}{\begin{array}[]{c}x_{t_{1}}\\\ x_{t_{2}}\\\
x_{t_{3}}\end{array}}&{\left[\begin{array}[]{c}PM_{t_{1}}^{(k)}\\\
PM_{t_{2}}^{(k)}\\\
PM_{t_{3}}^{(k)}\end{array}\right]}=\left[PM_{t_{i}}^{(k)}\right]_{3\times
1},k\in V_{t}.\end{array}}\end{array}$

The method of majority decision is a social welfare function satisfying
Arrow’s Conditions 2-5 and the consistency condition for any number of
alternatives, provided every triple of alternatives satisfy

$\exists
t_{i}\left(\big{|}\cup_{k=1}^{|V_{t}|}PM_{t_{i}}^{(k)}\big{|}<3\right),|V_{t}|\text{
is odd};t=1,2,\ldots,\binom{m}{3}.$ $None$

To illustrate, we consider an example used by Sen (1966).

Example 2 Five voters provide the following preference orderings over four
alternatives $\\{w,x,y,z\\}$:

$\begin{array}[]{ll}\text{Voter 1}:&w\sim x\succ y\succ z,\\\ \text{Voter
2}:&x\sim w\succ z\succ y,\\\ \text{Voter 3}:&z\sim x\succ y\succ w,\\\
\text{Voter 4}:&z\succ y\sim x\succ w,\\\ \text{Voter 5}:&z\succ y\succ x\succ
w.\end{array}$

There are 4 possible alternative triples and they are all Value-Restricted,
and thus the majority decision yields a transitive social preference ordering.

Here we examine the triples with Inequality (3). With regard to triple
$(w,x,y)$, the voters’ preferences are represented by PMs:

$\begin{array}[]{c}{\hskip 28.45274pt\begin{array}[]{cccccc}\text{Voter
1}:&\hskip 5.69046pt\text{Voter 2}:&\hskip 2.84544pt\text{Voter 3}:&\hskip
0.56917pt\text{Voter 4}:&\hskip 0.56917pt\text{Voter 5}:\end{array}}\\\
{\begin{array}[]{cc}{\begin{array}[]{c}w\\\ x\\\
y\end{array}}&{\left[\begin{array}[]{c}\\{1,2\\}\\\ \\{1,2\\}\\\
\\{3\\}\end{array}\right],\left[\begin{array}[]{c}\\{1,2\\}\\\ \\{1,2\\}\\\
\\{3\\}\end{array}\right],\left[\begin{array}[]{c}\\{3\\}\\\ \\{1\\}\\\
\\{2\\}\end{array}\right],\left[\begin{array}[]{c}\\{3\\}\\\ \\{1,2\\}\\\
\\{1,2\\}\end{array}\right],\left[\begin{array}[]{c}\\{3\\}\\\ \\{2\\}\\\
\\{1\\}\end{array}\right].}\end{array}}\end{array}$

The union set of the second entries of the above PMs is $\\{1,2\\}$, whose
cardinality number is less than 3. Hence Inequality (3) is satisfied by the
triple $(w,x,y)$. Similarly, one can easily check that Inequality (3) is
fulfilled by the other three triples.

##  4 Sen’s condition in equations

In this section, Sen’s transitivity condition will be described in such a way
that the membership preference map (mPM) is employed to characterize the
Value-Restricted requirement.

In Theorem 1’, the concerned individual’s preferences over a triple of
alternatives are represented by PMs. From Definition 2 we know that a PM has a
one-to-one correspondence of mPM. A very natural way is to further represent
the PMs in Theorem 1’ by mPMs and the thoerem will then transform into another
form. Actually, since an element of a mPM indicates whether an alternative can
be ranked at a certain position or not, thus if a certain alternative in a
triple is Value-Restricted, then none of the concerned individuals would like
to rank the certain alternative at a certain position, i.e., in this case the
0-1 element that corresponds to the certain alternative and the certain
ranking position of each concerned individual’s mPM is 0. As a consequence of
this analysis, the problem of examining inequalities for a transitive social
ordering is equivalent to checking a set of equations. This is described
explicitly next as Theorem 1”.

Theorem 1” Assume that the alternative set is
$A=\\{x_{1},x_{2},\ldots,x_{m}\\}$ where $2<m<+\infty$, the individual set is
$N=\\{1,2,\ldots,n\\}$ where $2<n<+\infty$, and the individuals’ preferences
over the alternative set are weak orderings. Denote by
$(x_{t_{1}},x_{t_{2}},x_{t_{3}})$ the triples of alternatives, where
$t=1,2,\ldots,\binom{m}{3}$. Denote by $V_{t}$ the maximal set including
concerned individuals who are specified in terms of their preferences over
that triple $(x_{t_{1}},x_{t_{2}},x_{t_{3}})$, and represent their preferences
over that triple by PMs. Specifically, for triple
$(x_{t_{1}},x_{t_{2}},x_{t_{3}})$ the concerned individuals’ PMs are denoted
by

$\begin{array}[]{c}{\hskip 8.5359pt\begin{array}[]{cccc}\text{Concerned
individual k's PM with respect to triple
}(x_{t_{1}},x_{t_{2}},x_{t_{3}}):\end{array}}\\\
{\begin{array}[]{cc}{\begin{array}[]{c}x_{t_{1}}\\\ x_{t_{2}}\\\
x_{t_{3}}\end{array}}&{\left[\begin{array}[]{c}PM_{t_{1}}^{(k)}\\\
PM_{t_{2}}^{(k)}\\\
PM_{t_{3}}^{(k)}\end{array}\right]}=\left[PM_{t_{i}}^{(k)}\right]_{3\times
1},k\in V_{t}.\end{array}}\end{array}$

The concerned individual k’s PM is further represented by a mPM:

$\begin{array}[]{c}{\hskip 8.5359pt\begin{array}[]{cc}\text{Concerned
individual k's mPM with respect to triple
}(x_{t_{1}},x_{t_{2}},x_{t_{3}}):\end{array}}\\\
\begin{array}[]{c}{\begin{array}[]{cc}\\{1\\}\hskip 36.98866pt\\{2\\}\hskip
36.98866pt\\{3\\}\end{array}}\\\
{\begin{array}[]{cc}{\begin{array}[]{c}x_{t_{1}}\\\ x_{t_{2}}\\\
x_{t_{3}}\end{array}}&{\left[\begin{array}[]{ccc}mPM_{t_{1,1}}^{(k)}&mPM_{t_{1,2}}^{(k)}&mPM_{t_{1,3}}^{(k)}\\\
mPM_{t_{2,1}}^{(k)}&mPM_{t_{2,2}}^{(k)}&mPM_{t_{2,3}}^{(k)}\\\
mPM_{t_{3,1}}^{(k)}&mPM_{t_{3,2}}^{(k)}&mPM_{t_{3,3}}^{(k)}\end{array}\right]_{3\times
3},}\end{array}k\in V_{t}.}\end{array}\end{array}$

The method of majority decision is a social welfare function satisfying
Arrow’s Conditions 2-5 and the consistency condition for any number of
alternatives, provided every triple of alternatives satisfy

$\min_{i}\min_{j}\left\\{\sum_{k=1}^{|V_{t}|}mPM_{t_{i,j}}^{(k)}\right\\}=0,|V_{t}|\text{
is odd};t=1,2,\ldots,\binom{m}{3}.$ $None$

To illustrate, we reconsider Example 2.

With regard to triple $(w,x,y)$, the voters’ mPMs are:

$\begin{array}[]{c}{\hskip 17.07182pt\begin{array}[]{cccccc}\text{Voter
1}:&\hskip 17.07182pt\text{Voter 2}:&\hskip 14.22636pt\text{Voter 3}:&\hskip
14.22636pt\text{Voter 4}:&\hskip 14.22636pt\text{Voter 5}:\end{array}}\\\
\begin{array}[]{c}{\hskip 22.76228pt\begin{array}[]{cccccc}\\{1\\}\hskip
2.27626pt\\{2\\}\hskip 2.27626pt\\{3\\}&\hskip 2.84544pt\\{1\\}\hskip
2.27626pt\\{2\\}\hskip 2.27626pt\\{3\\}&\hskip 2.84544pt\\{1\\}\hskip
2.27626pt\\{2\\}\hskip 2.27626pt\\{3\\}&\hskip 2.84544pt\\{1\\}\hskip
2.27626pt\\{2\\}\hskip 2.27626pt\\{3\\}&\hskip 2.84544pt\\{1\\}\hskip
2.27626pt\\{2\\}\hskip 2.27626pt\\{3\\}\end{array}}\\\
{\begin{array}[]{cc}{\begin{array}[]{c}w\\\ x\\\
y\end{array}}&{\left[\begin{array}[]{ccc}1&1&0\\\ 1&1&0\\\
0&0&1\end{array}\right],\left[\begin{array}[]{ccc}1&1&0\\\ 1&1&0\\\
0&0&1\end{array}\right],\left[\begin{array}[]{ccc}0&0&1\\\ 1&0&0\\\
0&1&0\end{array}\right],\left[\begin{array}[]{ccc}0&0&1\\\ 1&1&0\\\
1&1&0\end{array}\right],\left[\begin{array}[]{ccc}0&0&1\\\ 0&1&0\\\
1&0&0\end{array}\right].}\end{array}}\end{array}\end{array}$

The sum of the above five 0-1 matrices is:

$\left[\begin{array}[]{ccc}2&2&3\\\ 4&4&0\\\ 2&2&2\end{array}\right],$

from which one can see that Eq.(4) is fulfilled in that the $(2,3)th$ element
of the sum matrix is 0. Likewise, it is easy to check that Eq.(4) is fulfilled
by the other three triples.

##  5 Concluding remarks

Sen’s transitivity condition is a sort of qualitative description. This
article attempted to describe it in a quantitative manner. Our way to do this
was to represent the concerned individual’s preferences over a triple of
alternatives by preference maps (PMs) and membership preference maps (mPMs).
Each of these presented a theorem establishing a quantitative description of
Sen’s value restriction condition constrained on concerned individuals’
preferences for the majority decision to yield a transitive social preference
ordering. The practical significance of our results lies in the testing of
Sen’s transitivity condition with a computer. This is particularly true when,
as presented in Theorem 1”, the condition is constructed based on 0-1
matrices.

## References

  * [1]
  * [2]
  * [3] Arrow K J. A difficulty in the concept of social welfare. Journal of political economy, 1950, 58(4): 328-346.
  * [4] Arrow K J. Social Choice and Individual Values. New York: Wiley, 1951.
  * [5] Black D. The Theory of Committees and Elections, Cambridge: Cambridge University Press, 1958.
  * [6] Hou F. A consensus gap indicator and its application to group decision making. Group Decision and Negotiation, 2015a,24(3), 415–428.
  * [7] Hou F. The prametric-based GDM selection procedure under linguistic assessments. Fuzzy Systems (FUZZ-IEEE), 2015 IEEE International Conference on. IEEE, 2015b: 1-8.
  * [8] Hou F, Triantaphyllou E. An iterative approach for achieving consensus when ranking a finite set of alternatives by a group of experts. European Journal of Operational Research, 2019, 275(2), 570-579.
  * [9] Inada K. A note on the simple majority decision rule. Econometrica, 1964, 525-531.
  * [10] Sen A K. A possibility theorem on majority decisions. Econometrica: Journal of the Econometric Society, 1966, 491-499.
  * [11] Sen A K. Collective Choice and Social Welfare. San Francisco: Holden Day, 1970; republished Amsterdam: Elsevier North-Holland, 1979.
  * [12] Ward B. Majority voting and alternative forms of public enterprise. The Public Economy of Urban Communities, 1965, 112-126.