11institutetext: Leiden Observatory, Leiden University, PO Box 9513, 2300 RA
Leiden, The Netherlands
11email<EMAIL_ADDRESS>22institutetext: European Space Research
and Technology Centre, Keplerlaan 1, 2200 AG Noordwijk, The Netherlands
22email<EMAIL_ADDRESS>33institutetext: Faculty of Aerospace
Engineering, Delft University of Technology, Kluyverweg 1, 2629 HS Delft, The
Netherlands

# Spectral characterisation of the extinction properties of NGC 3603 using
JWST NIRSpec.

Ciarán Rogers 11 Bernhard Brandl Guido de Marchi 11 3 3 22

(Received November 10, 2023; accepted April 15, 2024)

###### Abstract

Context. A necessary ingredient in understanding the star formation history of
a young cluster is knowledge of the extinction towards the region. This has
typically been done by making use of the colour-difference method with
photometry, or similar methods utilising the colour-colour diagram. These
approaches rely on adopting an extinction law with a given total-to-selective
extinction ratio $R(V)$, or determining a value of $R(V)$ through empirical
relationships. They also rely upon accurate spectral classification, reliable
stellar isochrones, and separating field stars from genuine cluster members.

Aims. The colour excess $E(B-V)$ can be independently determined by studying
the decrements of the recombination lines produced by the nebular gas. Having
access to many recombination lines from the same spectral series removes the
need of adopting an extinction curve. Rather, different extinction curves can
be trialled and the most appropriate one selected based on a minimum
$\chi^{2}$ procedure.

Methods. Using the Micro-Shutter Assembly (MSA) on board the Near InfraRed
Spectrograph (NIRSpec), multi-object spectroscopy was performed, yielding 600
nebular spectra from the Galactic massive star formation region NGC 3603. The
recombination line intensity ratios were used to determine independent values
of $E(B-V)$. A series of extinction curves were trialled ranging from $R(V)=2$
to $R(V)=8$. The appropriate value of $R(V)$ was adopted based on the minimum
$\chi^{2}$ procedure.

Results. The extinction characteristics of NGC 3603 are similar to other
Galactic HII regions like Orion, as well as starburst regions such 30 Doradus
in the Large Magellanic Cloud, in that we find a relatively large value of
$R(V)=4.8\pm 1.06$, larger than the Galactic average of $3.1$. We find a
typical value of $E(B-V)=0.64\pm 0.27$, significantly lower than values
determined in previous studies. We also present a stacked nebular spectrum
with a typical continuum S/N = $70$. This spectrum highlights the
recombination lines of the HII region, several s-process elements such as Kr
$III$ and Se $IV$, and molecular $H_{2}$ emission lines. This high S/N
spectrum can act as a helpful template for identifying nebular emission lines.

Conclusions. Using ratios of hydrogen recombination lines, we calculated the
value of $R(V)$, $E(B-V)$ and $A(V)$ for $>200$ lines of sight across NGC
3603. An extinction curve with a typical value $R(V)=4.8\pm 1.06$ is required
to explain the colour excess observed in the nebular spectra. This corresponds
to a typical $E(B-V)=0.64\pm 0.27$. This is significantly lower than what has
been found in previous extinction studies of NGC 3603.

###### Key Words.:

(ISM:) dust, extinction - (ISM:) HII regions

##  1 Introduction

Determining the age, distance, and star formation history of a stellar cluster
relies upon an accurate knowledge of the extinction towards the region.
Extinction is often characterised with the term
$R(V)\equiv\frac{A_{V}}{E(B-V)}$, known as the total-to-selective extinction.
Most regions in the Milky Way follow a typical Galactic extinction law with a
total-to-selective extinction ratio $R(V)=3.1$. However, young regions with
active ongoing star formation tend to exhibit higher values of $R(V)$,
implying that the typical size of dust grains is large in these regions
compared to the Galactic average. Indeed, many studies towards Galactic HII
regions, or starburst galaxies, find systematically higher values of $R(V)$
(e.g. Neckel & Chini 1981; Bautista et al. 1995; De Marchi & Panagia 2014;
Pandey et al. 2000; Calzetti et al. 2000) compared to the Galactic average of
$3.1$. This could be a result of photoevaporation of smaller grains due to
radiation from massive stars and/or the coagulation and growth from small to
large grains, aided by the higher densities in star forming regions (Swamy
1965). Gall et al. (2014) show that supernovae explosions are also highly
efficient means of producing large dust grains.
NGC 3603 is a giant Galactic HII region, $1-3$ million years old, located
$7\pm 1$ kpc away (Pandey et al. 2000; Melena et al. 2008; Pang et al. 2016).
The region hosts an extremely compact OB cluster at its centre, which has made
studying the central region in detail difficult due to limited angular
resolution. With about one-third the mass in stars of 30 Doradus in the Large
Magellanic Cloud (LMC), NGC 3603 is the closest system in the Galaxy to a
genuine starburst region Stolte et al. (2006). It is also the optically
brightest HII region in the galaxy.
Most studies on the star formation history of NGC 3603 rely upon optical
photometry, and make use of either the colour-difference method or two colour
diagrams (TCDs) to determine the extinction (e.g. Moffat 1983; Melnick et al.
1989; Pandey et al. 2000; Sung & Bessell 2004), referred to as (Mof83),
(Mel89), (Pan00), and (Sun04) going forward.
An alternative approach to determine the extinction law and colour excess
towards a region is to study the recombination lines of nebular spectra. By
comparing the ratios of hydrogen emission lines present in nebular spectra to
their theoretical values based on Case B recombination theory Osterbrock &
Ferland (2006), the extinction law and colour excess can be determined, based
on how much the measured ratios deviate from the theoretical ratios. In this
study, we obtained $600$ nebular spectra using the multi-object spectroscopy
(MOS) mode of NIRSpec on board the James Webb Space Telescope (JWST), using
the G235H grating and F170LP filter. This corresponds to a resolving power of
$\frac{\lambda}{\Delta\lambda}=2700$, and a wavelength range of $1.6-3.2\mu
m$. MOS mode utilises the novel Micro-Shutter Assembly (MSA), which consists
of $\sim 250\,000$ individual micro-shutters that can be configured open or
closed depending on the positions of sources in the field of view.
In Sect. 2 we discuss the target selection. In Sect. 3, the data reduction
process is explained. In Sect. 4, the results of our extinction analysis are
presented. We also present a high signal to noise ratio (S/N) stacked nebular
spectrum in Sect. 4.6. In Sect. 5, we discuss our results, assumptions,
uncertainties, and biases, and compare our results to the literature. In Sect.
6 our findings are summarised.

##  2 Selection of targets and observation strategy

The principle aim of the NGC 3603 NIRSpec observations was to gain insight
into the formation processes of the pre-main sequence (PMS) stars in the
region. One hundred stellar sources were selected for observation with NIRSpec
from an initial list of about 10000, originally observed with Hubble Space
Telescope’s (HST) Wide Field Camera 3 (WFC3), see Beccari et al. (2010). The
100 sources observed with JWST had been classified photometrically as either
PMS ($60/100$) or main sequence (MS) ($40/100$), based on $H{\alpha}$ excess
emission. For our JWST observations, each source was placed within a micro-
shutter, with the shutters directly above and below also open, forming a mini-
slit of three micro-shutters. The upper and lower shutters were opened in
order to measure the nebular background, which is crucial in a region like NGC
3603 where the nebular emission is bright. A three nod strategy was employed,
moving the source from the central shutter, to the upper shutter, and finally
to the lower shutter so that source spectra could be measured on different
regions of the detector.
A major constraint on the final source selection was the placement of each
source in a micro-shutter. Sources needed to be relatively isolated, so that
contaminant stars were not observed in the neighbouring background shutters.
Here a contaminant star is simply any star that happens to fall within one of
the open micro-shutters, but is not a target of observation. This phenomenon
was impossible to completely avoid in a region as dense with stars as NGC
3603. The impact of contaminants on our nebular spectra is discussed in the
appendix, section A.
Figure 1 shows a section of NGC 3603. The bright south-east pillar is visible,
along with many projected micro-shutters. Note the sharp drop in nebulosity in
the upper left hand corner of the image.

Figure 1: Region of NGC 3603 with the projected micro-shutter’s angle and
dimensions. This RGB image is composed of High Acuity Wide-field K-band Imager
(HAWK-I) F$212n$ $H_{2}$ (Red), HST $F160W$ (Green), and HST $F168n$ (Blue).
The stellar sources not discussed in this paper lie within the central micro-
shutters, while the upper and lower shutters observe the nebula. The zoomed in
section containing the source $3796$ illustrates the nodding pattern utilised.
JWST nods up and down, moving three neighbouring open micro-shutters north-
east and south-west, observing $6$ nebular regions in total. The position of
the OB cluster is indicated by the blue star and arrows.

The final selection of 100 sources from the original catalogue of 10000 was
primarily motivated by how isolated the sources were, so that contamination
could be kept to a minimum. The final observation strategy consisted of three
MSA configurations, with three nods per configuration. A full discussion of
the MS and PMS stars will follow in Rogers et al. (in prep).
For this study, we focus on the nebular spectra that were obtained. Given the
three nod strategy and two micro-shutters observing the nebular background per
stellar source, six nebular spectra were obtained per source, totalling 600
nebular spectra in total.
Due to the nodding pattern of the telescope, certain nebular regions were
observed multiple times, producing repeat spectra. The pointing accuracy of
JWST is $\leavevmode\nobreak\ 0.1\arcsec$ Lallo & Hartig (2022). Given the
distance to NGC 3603 of $7$ kpc, this small pointing offset results in
sampling regions that are close in angular distance but relatively far apart
($\sim 700$ AU) in terms of linear distance. For this reason we treat the
repeat observations as independent spectra. The typical S/N in the continuum
of each spectrum was $\sim 10$. The recombination lines are quite strong,
having a typical S/N of $\geq 30$, with the strongest lines $\geq 100$.

##  3 Data reduction

###  3.1 NIPS

The uncalibrated data was largely reduced with the ESA Instrument Team’s
pipeline known as NIRSpec Instrument Pipeline Software (NIPS) Alves de
Oliveira et al. (2018). NIPS is a framework for spectral extraction of NIRSpec
data from the count-rate maps, performing all major reduction steps from dark
current and bias subtraction to flat fielding, wavelength and flux calibration
and spectral extraction, with the final product being the 1D extracted
spectrum. One of the final steps before extraction is the rectification of the
spectrum. Since NIRSpec spectra do not run perfectly parallel with the
dispersion direction, as well as bending slightly along the detector, a
rectification step is performed on the so-called irregular 2D spectrum. This
step resamples the 2D spectrum on a linear and uniform grid. The rectified
spectrum, also called the regular 2D spectrum, is a count-rate map consisting
of 3817 pixels in the dispersion direction and seven pixels in the cross-
dispersion direction. The extraction step then simply collapses the 2D
spectrum by summing the seven pixels along each detector column. During the
data reduction process, it became apparent that contamination of the nebular
spectra had occurred. This rendered the affected nebular spectra unusable.
Contamination reduced the number of usable nebular spectra from $600$ to
$471$. A full discussion of this effect is given in appendix A.

##  4 Extinction and colour excess

###  4.1 Measured and theoretical line ratios

It is possible to measure the extinction in a region via the recombination
lines found in the nebular background spectra, (e.g. Calzetti 1997; Calzetti
et al. 1996). Interstellar extinction reddens spectra, extinguishing blue
wavelengths more than red wavelengths. The expected intensity ratio of two
recombination lines for a given electron temperature $T_{e}$ and electron
density $N_{e}$ is known from Case B recombination theory. The ”B” signifies
that the environment is optically thick to all Lyman series photons, that is,
these photons are absorbed. The alternative scenario is Case A, where these
photons are not absorbed Osterbrock & Ferland (2006). Case B is expected to be
valid for HII regions where the density of hydrogen is relatively high ($\geq
100cm^{-3}$). The deviation that the measured recombination line ratios
exhibit compared to Case B is a measure of the extinction that is taking
place.

###  4.2 The choice of recombination lines

A strong Brackett series line in our nebular spectra could not be used in this
extinction analysis - $Br_{10}$. This was because it is located in a region of
the spectrum that is visibly affected by the calibration, and hence we could
not measure a reliable line flux. Figure 2 illustrates this, showing four
nebular spectra from neighbouring micro-shutters.

Figure 2: Example of four nebular spectra exhibiting a raised continuum level
around $Br_{10-4}$ as a result of the NIRSpec calibration. This effect
produces unphysical recombination ratios with respect to Case B recombination,
and so this line has not been included in the remainder of the analysis.

The continuum below the line is visibly raised, as indicated with the dashed
black line tracing the correct continuum. The absolute change in the flux of
$Br_{10}$ is a few $\%$. However, even a small change in the flux of the line
can result in a significant change in the line ratio, which can significantly
affect the resulting value of $E(B-V)$. As such, we have excluded $Br_{10}$
from this analysis.
All of the lines used in this analysis belong to the Brackett series. Using
lines from the same series helps to ensure that the nebular conditions
($T_{e}$ and $N_{e}$) being probed are similar. Ideally all of the lines that
we have used would have the same upper energy level, so that the excitation
potential was the same, meaning that they arise from identical environments.
This was not possible from the lines available to us, and so we opted to rely
on a single hydrogen series, with all lines relatively close together in
wavelength, keeping the physical conditions producing these lines similar. A
recombination line list is given in table 1.

Recombination lines
---
Line |  $\lambda$ ($\mu m$) |  $\lambda^{-1}$ ($\mu m^{-1}$)
$Br_{11}$ |  1.6811 |  0.5948
$Br_{9}$ |  1.8274 |  0.5333
$Br_{8}$ |  1.9445 |  0.5143
$Br_{7}$ |  2.1661 |  0.4617
$Br_{6}$ |  2.6258 |  0.3808
Table 1: Recombination lines used to estimate the extinction towards each
source.

####  4.2.1 Using Pfund recombination lines to probe extinction

There is a rich forest of Pfund recombination lines present in the nebular
spectra. For a number of reasons, it has not been possible to study these
lines to determine the extinction characteristics of NGC 3603. These lines lie
at longer wavelengths than all of the Brackett lines except for $Br_{6}$. This
means the extinction experienced by these lines is weaker, and hence harder to
measure. The lines are also relatively close together in wavelength, making
the differential extinction harder to measure. The lines are also weaker,
resulting in larger uncertainties, which also hinders our determination of
$E(B-V)$. As we shown in section 4.5.1, due to the overall weak differential
extinction that we determine, the Pfund lines are not well suited to measuring
$E(B-V)$ in NGC 3603, though we acknowledge that it would provide a very
interesting secondary probe of extinction towards the region.

####  4.2.2 Measuring the recombination lines

For each extracted nebular spectrum, we used Monte Carlo to estimate the line
flux and uncertainty for the Brackett recombination lines, similar to that of
(Ashraf et al. 2023; Riedel et al. 2017). This approach involved fitting the
recombination lines many times, allowing the lines to vary within their
uncertainties for each iteration. For this study we found that the
recombination lines uncertainties were well sampled after $1000$ iterations.
The flux of each line was determined by calculating the median of the $1000$
flux measurements, and the uncertainty of the flux was the standard deviation
of the $1000$ measurements. We fit the recombination lines with Gaussian
profiles using the non-linear least squares fitting routine ”curve-fit” from
SciPy. In combination with the analytical expression for a Gaussian, this
fitting routine returns the best fitting values of the amplitude, standard
deviation, and central wavelength of the Gaussian. We determined the
equivalent width (EW) of the best fitting Gaussian profile, and converted the
EW to a flux by multiplying by the level of the adjacent continuum. In some
cases the hydrogen recombination lines were partially blended with another,
weaker line. To estimate the impact of these weaker blended lines, we fit a
double Gaussian to both lines and found that the relative flux contribution
from the weaker lines was typically $<1\%$, so we do not consider it to
significantly affect our results. An example of the line profiles of the
Brackett series is shown in figure 3.

Figure 3: Five Brackett series line profiles are shown from a representative
spectrum, plotted in log-scale. The best fitting Gaussian profile is shown as
a red dashed line.

The flux of each hydrogen line was divided by the flux of $Br_{6}$, yielding a
list of line ratios. The uncertainty of each ratio was found by propagating
the error of each flux using the standard error propagation formula. The
$Br_{6}$ line was chosen as the anchor (the line that we compare the other
lines to) because it is present in all of the spectra (in some cases, $Br_{7}$
is located in the NIRSpec detector gap), and is always strongly detected,
resulting in low statistical uncertainties. The measured line ratios then
needed to be compared to the expected ratios from Case B recombination. To do
this, the line ratios of Storey & Hummer (1995) were used for values of
$T_{e}=10000-12500K$ and $n_{e}=100-10000cm^{-3}$. The line ratios are not
particularly sensitive to $T_{e}$ or $N_{e}$ and so our specific choice of
these parameters is not extremely important. Our choice of $T_{e}$ comes from
(Esteban et al. 2004; McLeod et al. 2016), who both find typical electron
temperatures of about $11000K$ in NGC 3603.
Equation 1 gives the relationship between the colour excess, the line ratios,
and the extinction curve,

$\displaystyle
E(B-V)=\frac{-log(r_{obs}/r_{int})}{0.4[\kappa(\lambda_{line1})-\kappa(\lambda_{line2})]}mag,$
(1)

where $r_{obs}$ is the measured ratio of two hydrogen recombination lines and
$r_{int}$ is the theoretical ratio of two hydrogen recombination lines based
on Case B recombination. We refer to the term: $log(r_{obs}/r_{int})$ as the
“decrement” moving forward. The notation of formula 1 is taken directly from
Calzetti (1996), and may lead to some confusion as the letter $r$ has been
used to represent the line ratios, while we use it to represent the extinction
curve in this study. In order to minimise confusion, we use lower case to
differentiate between $r_{obs}$/$r_{int}$ \- the line ratios, and $R(V)$ \-
the extinction curve. The terms $\kappa(\lambda_{line1})$ and
$\kappa(\lambda_{line2})$ are equivalent to $R(\lambda_{line1})$ and
$R(\lambda_{line2})$, that is, the values of the assumed extinction curve at
the wavelengths of the recombination lines in question. We refer to this term
as the “extinction curve differences” moving forward. Having numerous
recombination lines meant that each set of ratios could provide an independent
estimate of $E(B-V)$. The uncertainties of the line fluxes derived from the
Monte Carlo iteration method as well as the uncertainties in $T_{e}$ and
$N_{e}$ were propagated through to each value of $E(B-V)$ using the standard
error propagation formulae.

###  4.3 The choice of extinction curve

If an inappropriate extinction curve is used in equation 1 to determine
$E(B-V)$, different values of $E(B-V)$ will result from different line ratios.
The appropriate extinction curve on the other hand, should produce the same
value of $E(B-V)$ regardless of which line ratio is used, though this may not
hold for lines with very large wavelength separations. This caveat is fully
discussed in section 5. The form of the NIR extinction curve that we have
adopted is that of Fitzpatrick & Massa (2009) (FM09),

$\centering\frac{E(\lambda-V)}{E(B-V)}=[0.349+2.087R(V)]\times\frac{1}{1+(\lambda/0.507)^{\alpha}}-R(V),\@add@centering$
(2)

where $\alpha$ is the exponent that defines the curvature of the extinction
curve, and $R(V)$ has the usual definition of
$R(V)\equiv\frac{A_{V}}{E(B-V)}$. Previous studies of the Galactic extinction
law in the NIR indicated that it could be represented as a universal power law
with $A(\lambda)\propto\lambda^{-\beta}$, where the exponent $\beta$ is
constant for any line of sight (the usually agreed upon value of
$1.6\leq-\beta\leq 1.8$), see (Draine 1989; Sneden et al. 1978; Schultz &
Wiemer 1975). FM09 showed that $\beta$ is not only dependent on the line of
sight, but that the wavelength range used to probe the extinction actually
affects the value of $\beta$, and hence the value of $R(V)$ that is
determined. See section $3.2$ of FM09 for a full discussion. Owing to its
success in determining extinction towards a wide variety of regions, its
flexibility, and simplicity, we have adopted the form of the extinction curve
of FM09 using formula 2. This extinction curve mimics a power law whose
exponent increases with wavelength and can vary for different lines of sight,
while keeping the number of free parameters to just two - $R(V)$ and $\alpha$.
Here the exponent $\alpha$ is substituted for $\beta$ simply to highlight that
it is no longer constant. By using the FM09 extinction curve, which has been
characterised using only broadband photometry, we have assumed that it is
appropriate for our wavelength range. This can be confirmed in the coming
years with JWST spectroscopy.

###  4.4 The $\alpha$-$R(V)$ relation

Photometric extinction studies using high mass stars utilise both an expected
colour excess $E(B-V)$ and expected brightness $M$ for each star. Both of
these values are used in FM09 to determine $\alpha$ and $R(V)$. For nebular
observations, while Case B recombination theory can provide $E(B-V)$ via
equation 1, there is no expected nebular brightness. This prevents us from
directly determining $R(V)$ in equation 2. FM09 found that $\alpha$ and $R(V)$
are related, which allows us to determine $R(V)$ indirectly through $\alpha$.
To do this, we took the values and $1\sigma$ uncertainties for $\alpha$ and
$R(V)$ from FM09, and ran Monte Carlo simulations to fit a line between the
data points, and determine the uncertainty of the fitting parameters. The
$\alpha$ \- $R(V)$ relation is shown in figure 4.

Figure 4: Relationship between $R(V)$ and $\alpha$ determined by Fitzpatrick &
Massa (2009). The uncertainties in $\alpha$ and $R(V)$ are correlated, and the
direction and strength of the error bars reflect that. The orange shaded band
is the 1 $\sigma$ uncertainty of the fit, determined by performing 1000 Monte
Carlo iterations.

The uncertainties of $R(V)$ and $\alpha$ are correlated, as indicated by the
direction of the error bars in figure 4. We took care to only allow the data
points to vary along the directions indicated by their error bars. We
performed the Monte Carlo simulations multiple times, fitting different types
of curves to the data. These included a first- and second-order polynomial, as
well as a power law. We found that these choices did not significantly affect
the final result. The goodness of fits for each curve were very similar, and
the best fitting curve closely resembled a straight line in all three cases.
Hence, we chose to continue with the first-order polynomial. The equation of
the line of best fit is

$\centering\alpha=-0.287(\pm 0.018)\times R(V)+3.309(\pm
0.07).\@add@centering$ (3)

Using this relation, we replaced $\alpha$ in equation 1 with equation 3,
reducing it to a single parameter equation. Having done this, we were now able
to begin fitting different extinction curves to solve for $R(V)$.

###  4.5 Determining $R(V)$ and $E(B-V)$

Following equation 1, in order to calculate $E(B-V)$, the decrement:
$log(r_{obs}/r_{int})$ is divided by the extinction curve difference:
$\kappa(\lambda_{line1})-\kappa(\lambda_{line2})$. In order for each line
ratio to produce the same value of $E(B-V)$, such that the dispersion is zero,
the value of the decrements must be equal to, or be a constant multiple of the
extinction curve differences. As such, by varying $R(V)$, the dispersion in
$E(B-V)$ changes accordingly. By minimising the dispersion, the best fitting
value of $R(V)$ can be found. The process of minimising the dispersion in
$E(B-V)$ is analogous to performing the $\chi^{2}$ statistic test and
minimising $\chi^{2}$, where
$\chi^{2}\equiv\Sigma\frac{(data_{i}-model_{i})^{2}}{model_{i}}$.
Using equation 2, we produced $100$ extinction curves with $R(V)$ going from
$2$ to $8$. To determine the best fitting value of $R(V)$, each decrement was
fit with all $100$ extinction curves, with each extinction curve being
normalised to the level of the first decrement ($Br_{11}/Br_{6}$). By
normalising, any offsets in the y-axis between the extinction curves are
removed, and the difference between the decrements (data) and the extinction
curve (model) are set entirely by the value of the exponent $\alpha$ in
equation 2, which we express in terms of $R(V)$ using equation 3. Through
this, by varying the value of $R(V)$, we could find the value that produced
the minimum $\chi^{2}$ statistic.
To determine $R(V)$ for each decrement, we ran a Monte Carlo simulation,
allowing the value of the decrements to vary within their uncertainties. In
order to fully sample the uncertainties of the decrements, we needed to
generate $10\,000$ realisations per spectrum. This larger number of iterations
compared to section 4.2.2 is due to the decrement uncertainties being larger
than the line flux uncertainties. With each iteration, all $100$ normalised
extinction curves were fit to the current realisation of the four decrements,
calculating the $\chi^{2}$ statistic each time. The minimum $\chi^{2}$ and
corresponding value of $R(V)$ were saved before moving on to the next
realisation. This resulted in $471\times 10\,000$ $\chi^{2}$s. We then
calculated the median $\chi^{2}$ and $R(V)$ for each of the $471$ spectra. For
all $471$ nebular spectra, a typical value of $R(V)=4.9\pm 1.5$ was found.
After careful inspection of the decrements with respect to their best fitting
extinction curves, we chose to filter our sample of 471 lines of sight. This
was based on a combination of low S/N spectra (emission lines with S/N$<3$),
resulting in large line ratio uncertainties, and spectra that produced
unphysical jumps/spikes in the decrement trend. The reason for these
jumps/spikes is still not clear to us, but may originate from unresolved
issues with the NIRSpec calibration. We ordered the spectra in terms of
increasing $\chi^{2}$, and chose to consider the best $200$. In general,
beyond this point the decrements exhibited more unphysical spikes/jumps, with
uncertainties becoming very large, of order $\geq 50\%$. An example of
“reliable” and “unreliable” decrements are shown in figure 5. From the $200$
most reliable decrements, we found the typical value of $R(V)=4.8\pm 1.06$.
The nominal value of $R(V)$ is hardly altered by this filtering procedure,
indicating that we have not introduced biases into $R(V)$ by only considering
certain spectra. We also explored whether the smallest $\chi^{2}$ values
corresponded to the brightest nebular spectra (with the highest S/N), as there
is likely a correlation between brighter spectra and regions with lower
extinction/colour excess. However, we did not find any significant correlation
between the S/N of the spectra and their resulting $\chi^{2}$. As such, we
felt it was justified to exclude the sight lines that increased our
uncertainty in $R(V)$.
The $\chi^{2}$ statistic was used in this analysis as it is a familiar and
readily understood measure of goodness-of-fit. However, the minimisation of
dispersion in $E(B-V)$ has a more immediate physical meaning, in that the
colour excess along the line of sight should be constant, and hence dispersion
should be zero when the appropriate extinction curve is used. We performed the
above analysis a second time, minimising the dispersion in $E(B-V)$ rather
than $\chi^{2}$, and filtered our spectra based on increasing values of
dispersion. By minimising the dispersion in $E(B-V)$, a typical value of
$R(V)=4.97\pm 1.1$ was found, and $198/200$ spectra were overlapping between
the two approaches, showing that the methods are effectively equivalent.

(a) Set of decrements in which a reliable fit with low uncertainty can be
found for the extinction curve with $R(V)=5.0$.

(b) Set of decrements in which both the size of the error bars and the over
all trend of the data points results in a poor fit for any extinction curves.
The best fit is found for $R(V)=2.8$

Figure 5: Two sets of decrements are shown as black solid lines. The
uncertainties of each decrement are indicated with black error bars. Four
normalised extinction curves are also shown. $R(V)=2$ is represented by the
blue, sparsely dotted line. $R(V)=3$ is the magenta dot-dashed line. $R(V)=5$
is the green dashed line, and $R(V)=8$ is the red dotted line. Due to the
relatively small differences in slope between extinction curves with $R(V)=2$
and $R(V)=8$ at these wavelengths, small uncertainties in the line ratios are
required in order to reliably determine the best fitting value of $R(V)$

####  4.5.1 Determining $E(B-V)$

To determine the colour excess for each spectrum, we fit our $200$ spectra
with their best fitting extinction curves. This was done via Monte Carlo as
before, to propagate the uncertainties of the decrements and $R(V)$ to
$E(B-V)$. This produced four values of $E(B-V)$ per spectrum. The median of
these four values was calculated to determine the final colour excess for each
spectrum. The median for all spectra was calculated, representing the typical
value of $E(B-V)$ for our $200$ sight lines. We found a typical value of
$E(B-V)=0.64\pm 0.27$.

####  4.5.2 The total extinction $A(V)$

Given the definition of $R(V)\equiv\frac{A(V)}{E(B-V)}$, having determined
both $R(V)$ and $E(B-V)$, we were able to calculate the total extinction
$A(V)$. We have found $A(V)=3.1\pm 1.5$ mag. This value represents both the
foreground extinction due to the diffuse ISM, as well as the local extinction
happening within NGC 3603. In order to separate these two distinct regimes of
extinction, we have used the Gaia data release 3 catalogue Brown et al. (2021)
of extinction values for $95$ sources located between
$0.12^{\circ}-0.5^{\circ}$ of NGC 3603, at a distance of $6-8$ kpc in order to
determine the foreground extinction due to the diffuse ISM. We have obtained
this data directly from ESA Sky. The visual extinction in the Gaia system,
$A_{0}$, has a wavelength of $0.547\mu m$, just slightly shorter than $A(V)$
at $0.551\mu m$. Due to this small difference in wavelength we consider
$A_{0}$ and $A(V)$ directly comparable. We have found a foreground extinction
of $A_{0}=1.96$ mag. No measurement uncertainties for $A_{0}$ are provided in
the Gaia catalogue. This gives a local extinction value of
$A(V)_{local}=1.1\pm 1.5$ mag. This value is consistent within its
uncertainties with zero. This suggests that all of the extinction that we have
measured towards the region may be occurring in the foreground, and that the
cluster itself does not contribute significantly towards the total extinction.

###  4.6 Stacked nebular spectrum

A high S/N, stacked nebular spectrum has been created by combining and
averaging all of the usable nebular spectra together, shown in figure 7. This
process made use of all $471$ uncontaminated spectra. The spectra were first
normalised, setting the continuum level to unity. They were then combined
together by averaging each wavelength bin. In cases where the flux of a
spectrum in a given wavelength bin deviated by $\geq 10\%$ of the median
value, the wavelength bin for that spectrum was rejected. This removed noisy
pixels, cosmic rays, and other spikes present in the spectra. The result was
an average nebular spectrum for NGC 3603, with a typical continuum S/N = $70$.
In this spectrum, the Paschen, Brackett and Pfund series are prominently in
emission. The Pfund discontinuity can be seen clearly in this spectrum at
$\lambda\geq 2.3\mu m$. He I emission lines are also prominent throughout the
spectrum, along with some weaker HeII emission. The ions Kr $III$ at
$\lambda=2.1986\mu m$ and Se $IV$ at $\lambda=2.287\mu m$, produced by the
s-process in Asymptotic Giant Branch stars are strongly detected in the
nebular spectra. These ions are commonly seen in planetary nebula, (e.g.
Sterling et al. 2017). There are also more than $40$ molecular hydrogen
$H_{2}$ lines present, of which a selection have been annotated in figure 7.
These lines are only seen in a subset of the nebular spectra. Typically these
spectra correspond to regions of enhanced gas density, such as the south-east
pillar shown in figure 1, though the lines are weakly detected in more diffuse
regions as well.
This spectrum has been presented to showcase the rich NIR emission lines
produced by the HII region, and we hope that it will be a useful template for
identifying emission lines from other nebulae/HII regions.

##  5 Discussion

###  5.1 Discussion of assumptions and uncertainties

In this study we have made several assumptions in order to determine $R(V)$
and $E(B-V)$ using formula 2. Here a discussion of those uncertainties is
given.

####  5.1.1 Should $E(B-V)$ be constant for each line ratio?

The basis of our analysis relies upon the correct extinction curve producing a
constant value of $E(B-V)$, regardless of which line ratio is used to
determine it. In Bautista et al. (1995) (B95), the authors measured the
extinction $A(V)$ in Orion, using the nebular hydrogen emission lines from the
Balmer, Paschen and Brackett series. They found that the value of $A(V)$
increased systematically with the mean wavelength of the series, with
$A(V)=1.03\pm 0.05$, $A(V)=1.22\pm 0.10$, and $A(V)=1.86\pm 0.18$ for Balmer,
Paschen and Brackett respectively. The authors suggest that this may be a
result of optical effects due to the complex dust cloud geometry of the
region. We also speculate that it may be a result of higher order spectral
series (eg. Brackett, Pfund, Humphrey) coming preferentially from
hotter/denser regions of the nebula, where the extinction properties may
change. In (B95), the authors go on to say that the increase in $A(V)$ only
becomes significant when the Balmer series is compared to the Brackett series
(difference in $A(V)>3\sigma$). For a single value of $R(V)$, the change in
$A(V)$ must be accompanied by a proportional change in $E(B-V)$. Based on
this, if spectral series separated by large wavelengths are used to probe the
extinction properties of a region, then additional care must be taken to
derive $R(V)$. We feel it is justified in our case to assume a constant
$E(B-V)$ for each line of sight, as we make use of just a single spectral
series.

####  5.1.2 Calzetti equations assumptions

We have used equation 2 to calculate $E(B-V)$, which itself relies upon two
assumptions. These are: (1) the extinction that we are measuring is a purely
foreground, uniform screen of dust, and (2) that our extinction law is
appropriate. We know that (1) may not be completely valid in our case. The
morphology of the gas and dust in NGC 3603 is visibly clumpy and irregular
(though the intra-cluster extinction is consistent with zero from our results,
albeit with large uncertainties). For a uniform foreground dust screen, all
scattered and absorbed light is removed from the line of sight, with none of
the photons being redirected back towards the telescope. The more complex
geometry of dust clouds within NGC 3603 region means that scattered light can
indeed be redirected back towards the telescope, including of course the
emission lines that we use to measure the colour excess. We may therefore be
measuring photons that have experienced scattering by dust, with blue photons
being preferentially scattered compared to red photons, and obtain an
underestimate of the amount of colour excess. The effect of scattering on
extinction and colour excess due to different dust cloud geometries, degree of
clumpiness, and aperture size has been investigated by Krügel (2009) (K09).
One of their final conclusions is that scattering has a non-negligible impact
on colour excess and extinction when the scattering optical depth $\tau_{sca}$
within the telescope aperture is larger than the relative observational error.
The expression for $\tau_{sca}$ is given in (K09): $\tau_{sca}=\frac{1}{2}\phi
Dn\kappa_{sca}$ where $\phi$ is the spatial resolution of the telescope in
radians, $D$ is the distance to the region in $cm$, $n$ is the number density
in $cm^{-3}$ and $\kappa_{sca}$ is the scattering coefficient. Using this
expression we can show that thanks to the exquisite resolution of JWST along
with the relative proximity of NGC 3603 (compared to extragalactic sources),
scattering should not significantly impact our results. For JWST at
$\lambda=2\mu m$, $\phi=0.064\arcsec=3.15\times 10^{-7}$ radians. The distance
to NGC 3603 is $D=2.16\times 10^{22}cm$, the number density is
$10000cm^{-3}\leq n\leq 100cm^{-3}$ and the standard Milky Way value of
$\kappa_{sca}\approx 2.5\times 10^{-22}cm^{-2}$. This gives ${8.4\times
10^{-5}\leq\tau_{sca}\leq 8.4\times 10^{-3}}$. The upper limit is
approximately an order of magnitude smaller than our relative observational
error of $\sim 6\times 10^{-2}$, where our relative observational error refers
to the typical relative uncertainty of each decrement. This indicates that the
effect of scattering does not significantly affect our results. We acknowledge
that without a complete radiative transfer computation, we cannot precisely
quantify the effect of scattered light, but based on the estimate of
$\tau_{sca}$ above, the impact should not be significant for JWST observations
of NGC 3603.
For assumption (2), rather than adopting an extinction curve from the
literature, we have tested a range of extinction curves, and so we are
confident that this second assumption is satisfied in our analysis. A final
point to note on the determination of $E(B-V)$, is that NGC 3603 has a very
low Galactic latitude $\sim-0.5^{\circ}$ Fukui et al. (2013), lying well
within the Galactic plane. As such, it is possible that other Galactic star
clusters may be incidentally located behind NGC 3603, which could affect the
extinction properties that we determine. However, given that NGC 3603 is the
optically brightest HII region in the galaxy, we feel it is a safe assumption
that it is the dominant source of line emission along our lines of sight.

###  5.2 Comparing $E(B-V)$ with the literature

The value of $E(B-V)=0.64\pm 0.27$ that we have found is significantly lower
than other estimates in the literature. In (Mof83), they found $E(B-V)=1.44\pm
0.09$ by spectrally classifying nine sources and using existing photometry to
determine the colour excess. In (Mel89), they found $E(B-V)=1.44\pm 0.13$
within $50\arcsec$ of the central OB cluster. In that study they determine
$E(B-V)$ with photometric data using the pair method, also known as the
colour-difference method or MK method. This technique compares the observed
colour and the theoretical colour of zero age main sequence stars that have
been spectrally classified. They obtained a similar result using their nebular
spectral data using the recombination lines, but critically, they assumed a
Galactic extinction law with $R(V)=3.2$. In (Pan00) they found $E(B-V)=1.41\pm
0.13$ for sources within $50\arcsec$ of the central cluster utilising the P-Q
method from Ogura & Ishida (1981). They derive a value of $R(V)=4.3\pm 0.2$
for NGC 3603. When they included NIR photometry from Eisenhauer et al. (1998)
the value of $R(V)$ increased to $4.2\pm 0.7$, $4.8\pm 0.8$, and $5.2\pm 1.0$,
for the J, H and K filters respectively. The uncertainty was much higher for
these NIR observations, and so (Pan00) adopt their own value of $R(V)=4.3\pm
0.2$ owing to its lower uncertainty. In (Sun04), it was found that
$E(B-V)=1.25$ using the colour-difference method, and using the relation 4
from Guetter & Vrba (1989),

$\centering R(V)=2.45\times\frac{E(V-I)}{E(B-V)},\@add@centering$ (4)

they found a value of $R(V)=3.55\pm 0.12$. They also determined a distance
modulus $V_{0}-M_{V}=14.2\pm 0.2$ for NGC 3603 based on their estimate of the
colour excess. Pang et al. (2016) (Pang16) also found a value of $R(V)$ that
is higher than the typical galactic value. Within the filter system for HST,
they found $R_{f555w}=3.75\pm 0.87$ and a value of $E(f435W-f555W)=1.33\pm
0.12$ mag. They do this using the pair method to determine $E(f435W-f555W)$,
and perform a similar step to determine the total extinction $A_{f555W}$ by
comparing the measured magnitude and theoretical magnitude of zero age main
sequence stars. This allows them to directly calculate $R_{f555W}$ from the
relation $R_{f555W}\equiv\frac{A_{f555W}}{E(f435W-f555W)}$. We note that they
rely upon the distance modulus $V_{0}-M_{V}$ from (Sun04) to determine the
total extinction, but the distance modulus from (Sun04) was determined based
on a different colour excess law than that found by (Pang16). See Melena et
al. (2008) for a full summary of selective extinction values from literature
for NGC 3603.

###  5.3 Spatial sampling of nebular spectra

A possible explanation behind why our value of $E(B-V)$ is so much lower than
what has been found in previous studies may be the spatial distribution of the
nebular spectra in our sample. The extinction studies mentioned in 5.2 make
use of the colours of massive stars that are close to the cluster centre
($\leq 50\arcsec$). Our nebular spectra sample a wider area, with the majority
of our spectra coming from regions $\geq 60\arcsec$ from the cluster core, and
some spectra coming from as far as $\sim 100\arcsec$ from the cluster core.
Many of our lines of sight sample denser regions of the nebula, where large
clouds of gas and dust can be seen. This is qualitatively different from the
inner $50\arcsec$ of the cluster, where the density of gas and dust is much
lower. This is evident in figure 1, in the top left of the image, where the
nebular emission drops off sharply beyond the south-east pillar. This
variation in gas and dust density could lead to different extinction
properties across the region. Indeed, we have found that the best fitting
value of $R(V)$ and resulting $E(B-V)$ varies across the region, with a
$1\sigma$ spread of $\pm 1.25$ and $\pm 0.6$ respectively. As such, we believe
it is likely that our lower typical value of $E(B-V)$ is not contradicting
previous studies, rather it is reflective of the different regions of the
nebula that we are probing.

##  6 Conclusions

Using the NIR nebular spectra of NGC 3603, we have determined a new extinction
law towards the giant HII region, mostly corresponding to regions $\geq
50\arcsec$ from the cluster centre. We have done this by exploiting the strong
Brackett recombination lines present in our nebular spectra obtained with JWST
NIRSpec. By fitting a range of extinction curves to each spectrum, we were
able to home in on the most appropriate extinction curve, which has the form
of FM09’s NIR extinction law. We have found a typical value of $R(V)=4.8\pm
1.06$ for the region as a whole. This corresponds to a colour excess of
$E(B-V)=0.64\pm 0.27$, which is about half the value published in previous
works. Our findings diverge from what has been found previously, but we
believe this is a result of the spatial sampling of the region, and highlights
how spatially variable the extinction law is across the HII region.
We have also created a stacked nebular spectrum by averaging together all of
the usable individual nebular spectra, resulting in a typical continuum S/N =
$70$. This stacked spectrum highlights the emission lines produced by the
nebula of NGC 3603, and can be used as a template for identifying otherwise
hard to detect emission lines in fainter nebulae.

## Acknowledgements

We thank our referee for their suggestions and comments which have greatly
improved both the analysis presented here, as well as the readability and flow
of the paper.

## References

  * Alves de Oliveira et al. (2018) Alves de Oliveira, C., Luetzgendorf, N., Ferruit, P., & Rawle, T. 2018
  * Ashraf et al. (2023) Ashraf, M., Jose, J., Herczeg, G., & Fang, M. 2023, Journal of Astrophysics and Astronomy, 44, 67
  * Bautista et al. (1995) Bautista, M. A., Pogge, R. W., & DePoy, D. 1995, The Astrophysical Journal, 452, 685
  * Beccari et al. (2010) Beccari, G., Spezzi, L., De Marchi, G., et al. 2010, The Astrophysical Journal, 720, 1108
  * Brown et al. (2021) Brown, A. G., Vallenari, A., Prusti, T., et al. 2021, Astronomy & Astrophysics, 649, A1
  * Calzetti (1996) Calzetti, D. 1996, arXiv preprint astro-ph/9610184
  * Calzetti (1997) Calzetti, D. 1997, 408, 403
  * Calzetti et al. (2000) Calzetti, D., Armus, L., Bohlin, R. C., et al. 2000, The Astrophysical Journal, 533, 682
  * Calzetti et al. (1996) Calzetti, D., Kinney, A. L., & Storchi-Bergmann, T. 1996, The astrophysical journal. Chicago. Vol. 458, no. 1, pt. 1 (Feb. 1996), p. 132-135
  * De Marchi & Panagia (2014) De Marchi, G. & Panagia, N. 2014, Monthly Notices of the Royal Astronomical Society, 445, 93
  * Draine (1989) Draine, B. 1989, in Infrared spectroscopy in astronomy, Vol. 290
  * Eisenhauer et al. (1998) Eisenhauer, F., Quirrenbach, A., Zinnecker, H., & Genzel, R. 1998, The Astrophysical Journal, 498, 278
  * Esteban et al. (2004) Esteban, C., García-Rojas, J., Peimbert, M., et al. 2004, The Astrophysical Journal, 618, L95
  * Fitzpatrick & Massa (2009) Fitzpatrick, E. L. & Massa, D. 2009, The Astrophysical Journal, 699, 1209
  * Fukui et al. (2013) Fukui, Y., Ohama, A., Hanaoka, N., et al. 2013, The Astrophysical Journal, 780, 36
  * Gall et al. (2014) Gall, C., Hjorth, J., Watson, D., et al. 2014, Nature, 511, 326
  * Guetter & Vrba (1989) Guetter, H. H. & Vrba, F. J. 1989, The Astronomical Journal, 98, 611
  * Krügel (2009) Krügel, E. 2009, Astronomy & Astrophysics, 493, 385
  * Lallo & Hartig (2022) Lallo, M. & Hartig, G. 2022, Technical Report JWST-STScI-008271, 8271
  * McLeod et al. (2016) McLeod, A., Gritschneder, M., Dale, J., et al. 2016, Monthly Notices of the Royal Astronomical Society, 462, 3537
  * Melena et al. (2008) Melena, N. W., Massey, P., Morrell, N. I., & Zangari, A. M. 2008, The Astronomical Journal, 135, 878
  * Melnick et al. (1989) Melnick, J., Tapia, M., & Terlevich, R. 1989, Astronomy and Astrophysics (ISSN 0004-6361), vol. 213, no. 1-2, April 1989, p. 89-96., 213, 89
  * Moffat (1983) Moffat, A. 1983, Astronomy and Astrophysics, 124, 273
  * Neckel & Chini (1981) Neckel, T. & Chini, R. 1981, Astronomy and Astrophysics Supplement Series, 45, 451
  * Ogura & Ishida (1981) Ogura, K. & Ishida, K. 1981, Publications of the Astronomical Society of Japan, 33, 149
  * Osterbrock & Ferland (2006) Osterbrock, D. E. & Ferland, G. J. 2006, Astrophysics Of Gas Nebulae and Active Galactic Nuclei (University science books)
  * Pandey et al. (2000) Pandey, A. K., Ogura, K., & Sekiguchi, K. 2000, Publications of the Astronomical Society of Japan, 52, 847
  * Pang et al. (2016) Pang, X., Pasquali, A., & Grebel, E. K. 2016, The Astronomical Journal, 151, 23
  * Riedel et al. (2017) Riedel, A. R., Alam, M. K., Rice, E. L., Cruz, K. L., & Henry, T. J. 2017, The Astrophysical Journal, 840, 87
  * Rogers et al. (2022) Rogers, C. R., De Marchi, G., Giardino, G., et al. 2022, 12180, 1425
  * Schultz & Wiemer (1975) Schultz, G. & Wiemer, W. 1975, Astronomy and Astrophysics, vol. 43, no. 1, Sept. 1975, p. 133-139. Research supported by the Deutsche Forschungsgemeinschaft., 43, 133
  * Sneden et al. (1978) Sneden, C., Gehrz, R., Hackwell, J., York, D., & Snow, T. 1978, The Astrophysical Journal, 223, 168
  * Sterling et al. (2017) Sterling, N., Madonna, S., Butler, K., et al. 2017, The Astrophysical Journal, 840, 80
  * Stolte et al. (2006) Stolte, A., Brandner, W., Brandl, B., & Zinnecker, H. 2006, The Astronomical Journal, 132, 253
  * Storey & Hummer (1995) Storey, P. & Hummer, D. 1995, Monthly Notices of the Royal Astronomical Society, 272, 41
  * Sung & Bessell (2004) Sung, H. & Bessell, M. S. 2004, The Astronomical Journal, 127, 1014
  * Swamy (1965) Swamy, K. K. 1965, Publications of the Astronomical Society of the Pacific, 77, 164

##  Appendix A Impact of contaminant stars

Of the $600$ nebular spectra obtained with NIRSpec, $129$ were not usable due
to contamination by an incident star, leaving us with $471$ nebular spectra.
As mentioned in section 2, a limitation on the number of sources that could be
observed in a single MSA configuration came from the intense crowding of NGC
3603. In order to ensure that only well isolated sources were observed, with
clean neighbouring background shutters, a criterion was implemented which
specified that a source could be considered for observation so long as any
other nearby stars were at least five magnitudes less bright in the V band.
This approach was successful in so far that every target star had at least one
clean background shutter. However, the presence of contaminant stars, even
faint ones, has a significant impact on the extracted nebular background
spectrum. Figure 6 shows an extracted spectrum from a nebular background
shutter that contained a contaminant star, along with a spectrum from a clean
neighbouring shutter.

Figure 6: Two nebular spectra, one that has been extracted from a clean
nebular region free of contaminants, and another that is dominated by a
contaminant.

In figure 6 the upper blue spectrum is an example of a contaminated nebular
spectrum. The lower green spectrum is an example of a clean nebular spectrum,
with no contamination. Both of these spectra come from neighbouring micro-
shutters, and so they should look quite similar to each other in terms of over
all brightness and spectral slope. The most obvious difference between the
clean and the contaminated spectra is the shape of the continuum, which for
the contaminated spectrum, takes on a quasi-periodic oscillating profile. It
is not clear to us why this happens, however it is apparent that this occurs
only in cases where there is a contaminant star. This phenomenon was also seen
in previous work with simulated NIRSpec data Rogers et al. (2022). With the
simulated data, contaminants could be added or removed at will, and the same
oscillating pattern disappeared from the nebular spectra once the contaminants
were removed. It is also clear in the 2D spectrum of affected nebular spectra,
where a narrow, bright spectral trace can be seen overlaid on top of the
extended nebular emission.
We have attempted extracting contaminated spectra in different ways to try to
minimise or remove the effect of the contaminant. NIPS allows us to extract
the spectrum in multiple ways. The traditional approach is to simply sum the
flux in each detector column, but it is also possible to take a median value
of the flux in each column, and multiply this value by the number of pixels
per column (in our case for the 2D rectified spectra, this is seven pixels).
The benefit of using the median approach for nebular spectra is that, under
the assumption that the nebula is spatially uniform across the micro-shutter,
the S/N of the extracted spectrum can be improved. This is because there is
some slit loss around the edges of the shutter which can essentially be
corrected by taking the median flux value. In the case where the contaminant
spectrum is located near the upper or lower edge of the spectral trace, the
median pixel value would correspond to the nebular light, excluding the
contaminant flux. In the more common case of bright contaminant stars, whose
spectral trace can cover $\geq 50\%$ of the seven pixels per column, the
stellar flux dominates the spectrum, and taking the median effectively
considers only the stellar flux. For this reason we opted to stick with the
traditional sum method.
Another approach to remove the contaminant flux is to simply not count the
detector rows that contain the stellar flux. Rectifying the spectrum means
that the spectral trace runs perfectly parallel with the dispersion direction,
making it relatively simple to identify which rows contain the contaminant
flux, and skip them. Unfortunately, the resampling of the flux during
rectification does not work perfectly, with some pixels receiving more flux
than they should, and others receiving too little. The result is that
extracting only a portion of the nebular spectrum - the portion that does not
contain contaminant flux, results in its own oscillating continuum due to the
resampling. This is most easily seen in data like ours, with high S/N and high
resolution. Based on this, we opted to simply exclude any nebular spectrum
that contained a contaminant. This brought the number of usable nebular
spectra from $600$ to $471$.

##  Appendix B Stacked nebular spectrum

Figure 7: Stacked spectrum, resulting from the averaging and combining of the
individual nebular spectra.


Minimizing betatron coupling of energy spread and divergence in laser-
wakefield accelerated electron beams

A. Koehler1,2,*, R. Pausch1, J.P. Couperus Cabadağ1, O. Zarini1,2, J.M.
Krämer1, M. Bussmann1, A. Debus1, U. Schramm1,2, A. Irman1

1 Helmholtz-Zentrum Dresden - Rossendorf, Bautzner Landstrasse 400, 01328
Dresden, Germany
2 Technische Universität Dresden, 01062 Dresden, Germany

*<EMAIL_ADDRESS>
###### Abstract

Matched beam loading in laser wakefield acceleration (LWFA), characterizing
the state of flattening of the acceleration electric field along the bunch,
leads to the minimization of energy spread at high bunch charges. Here, we
demonstrate by independently controlling injected charge and acceleration
gradients, using the self-truncated ionization injection scheme, that minimal
energy spread coincides with a reduction of the normalized beam divergence.
With the simultaneous confirmation of a constant beam radius at the plasma
exit, deduced from betatron radiation spectroscopy, we attribute this effect
to the reduction of chromatic betatron decoherence. Thus, beam loaded LWFA
enables highest longitudinal and transverse phase space densities.

The concept of laser wakefield acceleration (LWFA) exploits the ultra-high
accelerating field gradients of up to a few hundred gigavolt-per-meter
generated in the wake of a high-intensity laser pulse as it propagates through
an optically transparent plasma [1, 2]. Electron bunches can thus be
accelerated to GeV energies within acceleration distances of only centimeters
[3]. Beam quality with respect to bunch charge, energy bandwidth, emittance
and pulse-to-pulse stability has improved substantially in the last decade and
is closely linked to a variety of controlled electron injection techniques [4,
5, 6, 7, 8, 9]. Only recently it was demonstrated that laser plasma
accelerators can be tailored for minimum energy spread at highest bunch
charges by reshaping the local accelerating field by matched beam loading [10,
9]. This combination of high charge, essential in the beam loading regime, and
the short bunch duration in the range of $10\text{\,}\mathrm{fs}$ [11, 12, 13,
14] result in high peak-currents exceeding $10\text{\,}\mathrm{k}\mathrm{A}$.
Future applications, such as high-field THz sources [15], laboratory-size
beam-driven plasma accelerators [16, 17, 18] and compact free-electron lasers
[19, 20, 21, 22, 23], will require such compact and further improved high-
brightness electron sources.

For taking full advantage of such high peak currents, not only the beam energy
spread has to be minimized, but also divergence and thus beam emittance [24,
25]. In this letter, we address this point by showing how beam decoherence,
schematically illustrated in Fig. 1, affects the transverse beam quality. In
transverse phase space $(x,p_{x})$ a polychromatic bunch can be divided into
slices representing different energy classes. These then rotate with different
angular frequencies according to the energy-dependent transverse focusing [26,
27]. A phase difference $\Delta\phi>\pi$ leads to full decoherence and thus
the maximum in occupied phase space area, i.e. in beam divergence and size.
Here, we demonstrate that beam loading besides minimizing the final energy
spread also decreases the beam divergence by about $20\text{\,}\mathrm{\char
37\relax}$. Betatron radiation spectroscopy simultaneously confirms a constant
beam radius of less than $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ at the
plasma exit. The small beam size and the minimized divergence improve the
transverse beam quality and thus beam emittance, facilitating beam transport
and novel applications.

Figure 1: Beam decoherence: (a) illustrates the transverse phase space of a
bunch with a finite energy spread consisting of many slices of different
electron energy which rotate in phase space. The beam envelope encircles the
sum of all slices. (b) shows the projection $dQ/dp_{x}$ of the beam on
$p_{x}$, which is typically recorded in experiments and used for measuring the
beam divergence. (c) shows the time-dependence that is caused by the phase
space rotation of the beam envelope. When the slices in (a) span over
$\Delta\phi>\pi$, full decoherence is reached and the modulations in (c)
vanish.

At the matched beam loading condition, a balance between the bunch’s self-
fields and the acceleration field of the wakefield is reached, resulting in a
constant longitudinal acceleration field along the bunch [28, 29], illustrated
in Fig. 4(c) presenting simulation results discussed later in the text. All
electrons within the bunch thus experience the same acceleration field so that
no energy spread is added during the process. Besides acceleration, the
wakefield provides strong linear focusing forces [30]. Off-axis electrons
perform transverse (betatron) oscillations around the beam axis and emit
betatron radiation while being accelerated [31, 32, 33, 34, 2]. Frequency,
amplitude, and phase of the betatron motion strongly depend on the energy of
electrons and their initial off-axis positions. At the end of the plasma
accelerator, the motion of electrons in $(x,p_{x})$ defines the final beam
size and transverse momentum distribution. Betatron radiation enables unique
access to the micron-scale beam size at the plasma exit [35, 33, 36]. Bunches
with finite energy spread experience decoherence of the betatron oscillations.
As a consequence the envelope of a polychromatic beam performs significant
oscillations [26, 27]. The beam loading condition can be effectively used to
control beam decoherence, and thus the final beam emittance.

Figure 2: Experimental setup: The laser is focused on a gas jet and drives a
wakefield. Accelerated electrons are energy analyzed in a magnet spectrometer.
An aluminum foil blocks residual laser light behind the spectrometer. The
scintillator of the betatron profiler intercepts the angular profile of the
betatron radiation. On-axis betatron radiation passes through a hole of the
profiler and is detected by the x-ray camera, which is separated from the
interaction chamber by a beryllium window. An aluminum filter foil attenuates
the betatron flux.

Experiments were performed with the DRACO Ti:Sapphire laser system [37] at the
Helmholtz-Zentrum Dresden – Rossendorf. Fig. 2 shows a schematic of the setup.
Laser pulses of $30\text{\,}\mathrm{f}\mathrm{s}$ full-width at half-maximum
(FWHM) duration with $2.5\text{\,}\mathrm{J}$ energy on target were focused
with an off-axis parabolic mirror (f/20) to a vacuum focal spot size of
$20\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ (FWHM) yielding a
normalized laser intensity of $a_{0}\simeq 2.6$. The Strehl ratio of the
wavefront-corrected beam was measured to be $0.9$ by the laser energy within
the beam waist ($e^{-2}$ of intensity). The laser beam was focused $\sim
1.5\text{\,}\mathrm{mm}$ beyond the entry of a $3\text{\,}\mathrm{mm}$ long de
Laval gas nozzle [38]. The nozzle was operated with He-N2 gas mixtures
containing $0.2\text{\,}\mathrm{\char 37\relax}1.5\text{\,}\mathrm{\char
37\relax}$ of N2. The gas density profile at $1.5\text{\,}\mathrm{mm}$ above
the nozzle exit, i.e., the laser beam axis, was characterized by a dedicated
interferometry setup [38] and consists of a $1.6\text{\,}\mathrm{mm}$ flat top
region. Electrons were injected into the wakefield using self-truncated
ionization injection (STII) [10, 39, 40], which enabled stable and
reproducible shots with high charges of up to $500\text{\,}\mathrm{pC}$
(within FWHM) required for entering the beam loading regime. In order to
measure the energy distribution, a $40\text{\,}\mathrm{c}\mathrm{m}$ long
permanent magnet dipole dispersed accelerated electrons to a set of charge
calibrated scintillator screens (Konica Minolta OG 400) [41] that was read-out
by charge coupled device (CCD) cameras.

A back-illuminated, deep depletion x-ray CCD (Princeton Instruments Pixis-XO
400BR) with
$1340\text{\times}400\text{\,}\mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l}$
recorded betatron radiation emitted from the LWFA process. The camera was
placed inside a dedicated radiation shielded area $12\text{\,}\mathrm{m}$
downstream of the plasma target, covering a solid angle of
$1.5\text{\times}0.8\text{\,}\mathrm{mrad}$. A
$76\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ thick beryllium window
sealed the CCD chip, allowing for cooling and background noise reduction. An
aluminum filter foil of $200\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$
thickness attenuated the betatron flux and enabled single-shot reconstruction
of the betatron spectrum by counting single-photon absorption events [42, 36,
43]. Behind the electron spectrometer, a scintillator screen oriented at
$45\text{\,}\mathrm{\SIUnitSymbolDegree}$ with respect to the beam axis
recorded the angle-resolved betatron profile.

Figure 3: Measured energy spread (a), normalized divergence
$\gamma\sigma_{\theta}$ (b), betatron radius $r_{\beta}\gamma^{-1/4}$ (c) and
calculated betatron phase difference $\Delta\phi$ according to Eq. (1) (d)
obtained from electron bunches with different charges and at three different
plasma densities. Every data point represents the average of up to $15$ shots
with constant experimental parameters.

The data sets presented in this letter were taken during the same experimental
campaign as those of Ref. [10] and follow a similar statistical significance
of up to 15 consecutive shots. The amount of injected charge was tuned by
varying the nitrogen concentration, keeping the laser parameters and plasma
density constant. Due to the beam loading effect, the energy spread of the
accelerated electron beams depends strongly on the injected charge, reaching a
minimum value at around $300\text{\,}\mathrm{p}\mathrm{C}$, as shown in Fig.
3(a). The same trend holds for the energy (see supplementary information
[44]).

In order to investigate the effect of beam loading on the transverse electron
dynamics, the geometric divergence of the bunch $\sigma_{\theta}$ was
extracted from the non-dispersive plane of the electron spectrometer after
having left the plasma accelerator. Fig. 3(b) shows the normalized divergence
$\gamma\sigma_{\theta}$ for various plasma densities where $\gamma$ is the
relativistic Lorentz factor. At around $300\text{\,}\mathrm{p}\mathrm{C}$
where the optimum beam loading occurs [10], a drop of up to
$20\text{\,}\mathrm{\char 37\relax}$ of the normalized beam divergence is
consistently observed for all plasma density sets. This suggests that
accelerator geometry, i.e., cavity size, plasma wavelength and the plasma
down-ramp length, does not play a dominant role in the observed effect.

Additionally, the betatron source size, thus the electron beam radius at the
end of the accelerator, was deduced from the betatron spectrum by modeling
synchrotron radiation [33, 36]. Errors caused by shot-to-shot fluctuation
resulting in jitter of the betatron pointing were avoided by excluding off-
center shots using the transverse betatron profile diagnostic. Only shots
within one data set were selected where the betatron radiation cone was
oriented on-axis. Fig. 3(c) presents the reconstructed betatron source radii
$r_{\beta}\gamma^{-1/4}$. The $\gamma^{-1/4}$ factor normalizes differences in
energy gain [2]. Note that the betatron radius expresses no extremum at
$300\text{\,}\mathrm{p}\mathrm{C}$ but indicates a constant micrometer-small
source size at the plasma accelerator exit for different experimental
parameters. For this spatially compact bunch, the normalized divergence shows
a fairly similar value for the lowest and highest charges. Thus, space charge
does not dominate the beam divergence as the bunch propagates in the drift
space after acceleration. Otherwise, a monotonic increase of the beam
divergence with charge should be visible.

Because space charge cannot explain the reduced beam divergence observed at
matched beam loading, we focus on the influence of beam energy spread on beam
decoherence. For optimum beam loading the relative beam energy spread reaches
a minimum value. The length required for the full decoherence can be estimated
by $l_{dc}=2.35\lambda_{\beta}{\langle E\rangle}/{\Delta E}$ [26, 45]. Here,
$\lambda_{\beta}=\sqrt{2\gamma}\lambda_{p}$ is the betatron wavelength,
$\lambda_{p}=2\pi c/\omega_{p}$ the plasma wavelength, $c$ the speed of light
in vacuum, $\omega_{p}$ the plasma frequency, $\Delta E$ the energy spread
(FWHM), and $\langle E\rangle$ the average energy of the peaked spectrum. For
the observed energy range, $l_{dc}$ is always longer than two millimeters
which is more than the acceleration length $l_{\textrm{acc}}$ of
$1.6\text{\,}\mathrm{mm}$, indicated by simulation [44]. The short
acceleration length suggests that the accelerated electron bunches have not
fully experienced beam decoherence at the end of the plasma channel. Thus, the
difference of the betatron phase between the highest and lowest energy
electrons can quantify the degree of decoherence.

In transverse phase space $(x,p_{x})$, electrons at given energy orbit at a
different frequency compared to electrons at other energies [27]. A group $i$
of electrons with the energy $\gamma_{i}$ gains the betatron phase
$\phi_{i}=\int^{l_{\textrm{acc}}}_{0}\text{d}t\omega_{\beta,i}$ over the
acceleration length, with $\omega_{\beta,i}=2\pi c/\lambda_{\beta,i}$ being
the betatron frequency. With the maximum and minimum phase for $i=1,2$, full
decoherence is equivalent to a phase difference
$\Delta\phi=\phi_{1}-\phi_{2}\geq\pi$. Furthermore, the duration of the
injection process determines the initial transverse phase space. Electrons
injected earlier have already performed phase space rotation while other
electrons are still being trapped. For the case of STII with an injection
duration shorter than $\pi/\omega_{p}$, beam size and transverse momentum
spread can exhibit oscillations in $(x,p_{x})$ [27]. The phase difference
$\Delta\phi$ is given by [27]:

$\Delta\phi\approx\frac{\sqrt{2}}{E_{z,1}}\sqrt{E_{z,1}k_{p}l_{\textrm{acc}}+1}-\frac{\sqrt{2}}{E_{z,2}}\sqrt{E_{z,2}k_{p}l_{\textrm{acc}}+1},$
(1)

where $k_{p}=2\pi/\lambda_{p}$ is the plasma wavenumber. The electric fields
$E_{z,i}$ at the bunch’s head and tail are deduced from the measured maximum
energy gain and energy spread and are normalized to $m_{e}\omega_{p}c/e$,
where $m_{e}$ and $e$ are the electron rest mass and the charge, respectively.
Fig. 3(d) presents the calculated phase difference according to Eq. (1). In
general, the minimum $\Delta\phi$ is in good agreement with the expected
optimum beam loading condition, appearing at slightly lower charges than the
minimum of $\gamma\sigma_{\theta}$. This demonstrates that a low energy spread
leads to the reduction of the beam decoherence resulting in electron beams of
small divergence.

An upper limit for the normalized emittance of
$\epsilon_{n}=\gamma\sigma_{\theta}r_{b}$ can be derived from the
simultaneously recorded beam divergence and radius $r_{b}$. The latter can be
deduced from $r_{\beta}$ by $r_{b}=r_{\beta}/\sqrt{2}$ [46]. Typical values of
the beam size are in the range of $0.56\text{\,}\mathrm{\SIUnitSymbolMicro
m}0.7\text{\,}\mathrm{\SIUnitSymbolMicro m}$ at the accelerator exit. Together
with $\gamma\sigma_{\theta}$ varying from
$4\text{\,}\mathrm{rad}5.6\text{\,}\mathrm{rad}$, the resulting normalized
emittance reaches
$2\text{\,}\mathrm{mm}\text{\,}\mathrm{mrad}4\text{\,}\mathrm{mm}\text{\,}\mathrm{mrad}$
which is comparable with typical values for LWFA when applying ionization
injection [2, 9].

Figure 4: Matched beam loading in simulation: (a) illustrates the focusing
force map with an injected charge of $300\text{\,}\mathrm{p}\mathrm{C}$ at
$z=$1.7\text{\,}\mathrm{mm}$$. The density of trapped and accelerated
electrons from the nitrogen’s K-shell and the focusing force are indicated by
the gray and color scale at the bottom right, respectively. (b) presents two
lineouts of the focusing force in front of the bunch (orange) and at the bunch
(blue) emphasizing its linearity. The black dashed line in (b) indicates the
theoretically predicted value [47]. (c) shows the accelerating fields for
injected charges of $300\text{\,}\mathrm{p}\mathrm{C}$ and
$210\text{\,}\mathrm{p}\mathrm{C}$. The black dotted lines illustrate the
local flattening of the acceleration field for matched beam loading close to
$300\text{\,}\mathrm{p}\mathrm{C}$.

Parallel to the experimental work, three-dimensional particle-in-cell (PIC)
simulations were performed with the PIConGPU code [48, 49] taking realistic
experimental parameters [44]. With these simulations the influence of beam
loading on the focusing fields and transverse phase space can be confirmed.
Fig. 4(a) illustrates the focusing field map from two plasma cavities loaded
with an injected charge of $300\text{\,}\mathrm{pC}$. Note that the focusing
force is extended to smaller $(z-ct)$ because of the large amount of injected
charge, demonstrating the beam loading effect [29]. The blue and orange
lineouts in Fig. 4(b) are sampling the focusing force at different
longitudinal positions, i.e., in the front of the bunch and at the bunch
center, respectively, emphasizing the linearity and independence of space
charge. Fig. 4(c) illustrates the effect of beam loading on the acceleration
field. The weak loading of a small bunch charge
($210\text{\,}\mathrm{p}\mathrm{C}$) maintains the typical linear field
gradient. With strong beam loading ($300\text{\,}\mathrm{p}\mathrm{C}$), the
field becomes constant along the bunch, minimizing the energy spread
accumulated during acceleration. Summarizing, beam loading does not directly
couple the longitudinal and transverse fields of the plasma cavity.

Figure 5: The transverse momentum $p_{x}$ in simulation: (a) indicates an
increase with weak modulations of the transverse momentum of the electrons
during acceleration. (b) shows the betatron oscillation of a macro-particle
subset. The black line shows the semi-analytic model from Ref. [30] assuming
an adiabatic acceleration.

Fig. 5 presents the temporal evolution of the transverse momentum for the
bunch and a subset of $\sim 4000$ macro-particles. As shown in Fig. 5(a), the
entire bunch experiences periodic modulations in $p_{x}$, indicating that the
bunch has not reached full decoherence yet. This confirms the estimate of the
decoherence length. The modulations can be attributed to the betatron
oscillations as denoted for a macro-particles subset in Fig. 5(b). The subset
contains macro-particles with same initial position and momentum. The black
curve represents the theoretically predicted evolution in phase space
according to [30] and agrees well with the simulated data.

In conclusion, we have shown that the experimental observation of the charge
dependent minimum in normalized beam divergence can be explained by incomplete
decoherence of the betatron oscillations in the bunch. The finding is an
essential step towards the generation of electron bunches with small
normalized divergence and energy spread at high peak current. Furthermore,
tuning the energy spread can control the beam divergence via the betatron
decoherence process. The next generation of high-flux light sources will
benefit from less divergent electron sources and beam transport to a
subsequent setup such as an insertion device will be facilitated.

The authors gratefully acknowledge the team from DRACO as well as T. Cowan, M.
Downer, J. Grenzer A. Laso García, K. Steiniger, A. Wagner and R. Zgadzaj for
their support. This work was partly supported by EuCARD-2, funding under Grand
Agreement No. 312453, and by HZDR under programme Matter and Technology. The
PIC simulations were performed at the Center for Information Services and High
Performance Computing (ZIH) at TU Dresden.

## References

  * Esarey  _et al._ [2009] E. Esarey, C. B. Schroeder, and W. P. Leemans, Rev. Mod. Phys. 81, 1229 (2009).
  * Downer  _et al._ [2018] M. C. Downer, R. Zgadzaj, A. Debus, U. Schramm, and M. C. Kaluza, Rev. Mod. Phys. 90, 035002 (2018).
  * Gonsalves  _et al._ [2019] A. J. Gonsalves, K. Nakamura, J. Daniels, C. Benedetti, C. Pieronek, T. C. H. de Raadt, S. Steinke, J. H. Bin, S. S. Bulanov, J. van Tilborg, C. G. R. Geddes, C. B. Schroeder, C. Tóth, E. Esarey, K. Swanson, L. Fan-Chiang, G. Bagdasarov, N. Bobrova, V. Gasilov, G. Korn, P. Sasorov, and W. P. Leemans, Phys. Rev. Lett. 122, 084801 (2019).
  * Geddes  _et al._ [2008] C. G. R. Geddes, K. Nakamura, G. Plateau, C. Toth, E. Cormier-Michel, E. Esarey, C. B. Schroeder, J. Cary, and W. P. Leemans, Physical Review Letters 100, 1 (2008).
  * Osterhoff  _et al._ [2008] J. Osterhoff, A. Popp, Z. Major, B. Marx, T. P. Rowlands-Rees, M. Fuchs, M. Geissler, R. Hörlein, B. Hidding, S. Becker, E. A. Peralta, U. Schramm, F. Grüner, D. Habs, F. Krausz, S. M. Hooker, and S. Karsch, Phys. Rev. Lett. 101, 085002 (2008).
  * Rechatin  _et al._ [2009] C. Rechatin, J. Faure, A. Ben-Ismail, J. Lim, R. Fitour, A. Specka, H. Videau, A. Tafzi, F. Burgy, and V. Malka, Phys. Rev. Lett. 102, 1 (2009).
  * Gonsalves  _et al._ [2011] A. J. Gonsalves, K. Nakamura, C. Lin, D. Panasenko, S. Shiraishi, T. Sokollik, C. Benedetti, C. B. Schroeder, C. G. R. Geddes, J. van Tilborg, J. Osterhoff, E. Esarey, C. Toth, and W. P. Leemans, Nat. Phys. 7, 862 (2011).
  * Buck  _et al._ [2013] A. Buck, J. Wenz, J. Xu, K. Khrennikov, K. Schmid, M. Heigoldt, J. M. Mikhailova, M. Geissler, B. Shen, F. Krausz, S. Karsch, and L. Veisz, Phys. Rev. Lett. 110, 185006 (2013).
  * Irman  _et al._ [2018] A. Irman, J. P. Couperus, A. Debus, A. Köhler, J. M. Krämer, R. Pausch, O. Zarini, and U. Schramm, Plasma Physics and Controlled Fusion 60 (2018), 10.1088/1361-6587/aaaef1.
  * Couperus  _et al._ [2017] J. P. Couperus, R. Pausch, A. Köhler, O. Zarini, J. M. Krämer, M. Garten, A. Huebl, R. Gebhardt, U. Helbig, S. Bock, K. Zeil, A. Debus, M. Bussmann, U. Schramm, and A. Irman, Nature Communications 8, 487 (2017).
  * Glinec  _et al._ [2007] Y. Glinec, J. Faure, A. Norlin, A. Pukhov, and V. Malka, Phys. Rev. Lett. 98, 194801 (2007).
  * Debus  _et al._ [2010] A. D. Debus, M. Bussmann, U. Schramm, R. Sauerbrey, C. D. Murphy, Z. Major, R. Hörlein, L. Veisz, K. Schmid, J. Schreiber, K. Witte, S. P. Jamison, J. G. Gallacher, D. A. Jaroszynski, M. C. Kaluza, B. Hidding, S. Kiselev, R. Heathcote, P. S. Foster, D. Neely, E. J. Divall, C. J. Hooker, J. M. Smith, K. Ertel, A. J. Langley, P. Norreys, J. L. Collier, and S. Karsch, Physical Review Letters 104, 1 (2010).
  * Lundh  _et al._ [2011] O. Lundh, J. Lim, C. Rechatin, L. Ammoura, A. Ben-Ismaïl, X. Davoine, G. Gallot, J.-P. Goddet, E. Lefebvre, V. Malka, and J. Faure, Nature Physics 7, 219 (2011).
  * Zarini  _et al._ [2018] O. Zarini, A. Köhler, J. Couperus, R. Pausch, T. Kurz, S. Schöbel, H. Meißner, M. Bussmann, U. Schramm, A. Irman, and A. Debus, in  _IEEE Advanced Accelerator Concepts Workshop (AAC)_ (IEEE, 2018) pp. 1–5, 978-1-5386-7721-6.
  * Green  _et al._ [2016] B. Green, S. Kovalev, V. Asgekar, G. Geloni, U. Lehnert, T. Golz, M. Kuntzsch, C. Bauer, J. Hauser, J. Voigtlaender, B. Wustmann, I. Koesterke, M. Schwarz, M. Freitag, A. Arnold, J. Teichert, M. Justus, W. Seidel, C. Ilgner, N. Awari, D. Nicoletti, S. Kaiser, Y. Laplace, S. Rajasekaran, L. Zhang, S. Winnerl, H. Schneider, G. Schay, I. Lorincz, A. A. Rauscher, I. Radu, S. Maehrlein, T. H. Kim, J. S. Lee, T. Kampfrath, S. Wall, J. Heberle, A. Malnasi-Csizmadia, A. Steiger, A. S. Müller, M. Helm, U. Schramm, T. Cowan, P. Michel, A. Cavalleri, A. S. Fisher, N. Stojanovic, and M. Gensch, Scientific Reports 6, 22256 (2016).
  * Rosenzweig  _et al._ [1991] J. B. Rosenzweig, B. Breizman, T. Katsouleas, and J. J. Su, Physical Review A 44 (1991), 10.1103/PhysRevA.44.R6189.
  * Hidding  _et al._ [2010] B. Hidding, T. Königstein, J. Osterholz, S. Karsch, O. Willi, and G. Pretzler, Physical Review Letters 104, 1 (2010).
  * Martinez De La Ossa  _et al._ [2013] A. Martinez De La Ossa, J. Grebenyuk, T. Mehrling, L. Schaper, and J. Osterhoff, Physical Review Letters 111 (2013), 10.1103/PhysRevLett.111.245003, 1309.1900 .
  * Schlenvoigt  _et al._ [2008] H.-P. Schlenvoigt, K. Haupt, A. Debus, F. Budde, O. Jäckel, S. Pfotenhauer, H. Schwoerer, E. Rohwer, J. G. Gallacher, E. Brunetti, R. P. Shanks, S. M. Wiggins, and D. A. Jaroszynski, Nature Physics 4, 130 (2008).
  * Fuchs  _et al._ [2009] M. Fuchs, R. Weingartner, A. Popp, Z. Major, S. Becker, J. Osterhoff, I. Cortrie, B. Zeitler, R. Hörlein, G. D. Tsakiris, U. Schramm, T. P. Rowlands-Rees, S. M. Hooker, D. Habs, F. Krausz, S. Karsch, and F. Grüner, Nature Physics 5, 826 (2009).
  * Maier  _et al._ [2012] A. R. Maier, A. Meseck, S. Reiche, C. B. Schroeder, T. Seggebrock, and F. Grüner, Phys. Rev. X 2, 031019 (2012).
  * Huang  _et al._ [2012] Z. Huang, Y. Ding, and C. B. Schroeder, Phys. Rev. Lett. 109, 204801 (2012).
  * Steiniger  _et al._ [2014] K. Steiniger, M. Bussmann, R. Pausch, T. Cowan, A. Irman, A. Jochmann, R. Sauerbrey, U. Schramm, and A. Debus, Journal of Physics B: Atomic, Molecular and Optical Physics 47, 234011 (2014).
  * Eichner  _et al._ [2007] T. Eichner, F. Grüner, S. Becker, M. Fuchs, D. Habs, R. Weingartner, U. Schramm, H. Backe, P. Kunz, and W. Lauth, Phys. Rev. ST Accel. Beams 10, 082401 (2007).
  * André  _et al._ [2018] T. André, I. A. Andriyash, A. Loulergue, M. Labat, E. Roussel, A. Ghaith, M. Khojoyan, C. Thaury, M. Valléau, F. Briquez, F. Marteau, K. Tavakoli, P. N’Gotta, Y. Dietrich, G. Lambert, V. Malka, C. Benabderrahmane, J. Vétéran, L. Chapuis, T. El Ajjouri, M. Sebdaoui, N. Hubert, O. Marcouillé, P. Berteaud, N. Leclercq, M. El Ajjouri, P. Rommeluère, F. Bouvet, J. P. Duval, C. Kitegi, F. Blache, B. Mahieu, S. Corde, J. Gautier, K. Ta Phuoc, J. P. Goddet, A. Lestrade, C. Herbeaux, C. Évain, C. Szwaj, S. Bielawski, A. Tafzi, P. Rousseau, S. Smartsev, F. Polack, D. Dennetière, C. Bourassin-Bouchet, C. De Oliveira, and M.-E. Couprie, Nature Communications 9, 1334 (2018).
  * Mehrling  _et al._ [2012] T. Mehrling, J. Grebenyuk, F. S. Tsung, K. Floettmann, and J. Osterhoff, Physical Review Special Topics - Accelerators and Beams 15, 1 (2012).
  * Xu  _et al._ [2014] X. L. Xu, J. F. Hua, F. Li, C. J. Zhang, L. X. Yan, Y. C. Du, W. H. Huang, H. B. Chen, C. X. Tang, W. Lu, P. Yu, W. An, C. Joshi, and W. B. Mori, Physical Review Letters 112, 1 (2014).
  * Katsouleas  _et al._ [1987] T. Katsouleas, S. Wilks, J. M. Dawson, P. Chen, and J. J. Su, Particle Accelerators 22, 210 (1987).
  * Tzoufras  _et al._ [2008] M. Tzoufras, W. Lu, F. S. Tsung, C. Huang, W. B. Mori, T. Katsouleas, J. Vieira, R. A. Fonseca, and L. O. Silva, Physical Review Letters 101, 1 (2008), 0809.0227 .
  * Khachatryan  _et al._ [2007] A. G. Khachatryan, A. Irman, F. A. Van Goor, and K. J. Boller, Physical Review Special Topics - Accelerators and Beams 10, 1 (2007).
  * Esarey  _et al._ [2002] E. Esarey, B. Shadwick, P. Catravas, and W. Leemans, Phys. Rev. E 65, 056505 (2002).
  * Albert  _et al._ [2008] F. Albert, R. Shah, K. T. Phuoc, R. Fitour, F. Burgy, J.-P. Rousseau, A. Tafzi, D. Douillet, T. Lefrou, and A. Rousse, Phys. Rev. E 77, 056402 (2008).
  * Schnell  _et al._ [2012] M. Schnell, A. Sävert, B. Landgraf, M. Reuter, M. Nicolai, O. Jäckel, C. Peth, T. Thiele, O. Jansen, A. Pukhov, O. Willi, M. C. Kaluza, and C. Spielmann, Phys. Rev. Lett. 108, 1 (2012).
  * Corde  _et al._ [2013] S. Corde, K. Ta Phuoc, G. Lambert, R. Fitour, V. Malka, A. Rousse, A. Beck, and E. Lefebvre, Reviews of Modern Physics 85, 1 (2013).
  * Plateau  _et al._ [2012] G. R. Plateau, C. G. R. Geddes, D. B. Thorn, M. Chen, C. Benedetti, E. Esarey, A. J. Gonsalves, N. H. Matlis, K. Nakamura, C. B. Schroeder, S. Shiraishi, T. Sokollik, J. van Tilborg, C. Toth, S. Trotsenko, T. S. Kim, M. Battaglia, T. Stöhlker, and W. P. Leemans, Phys. Rev. Lett. 109, 064802 (2012).
  * Köhler  _et al._ [2016] A. Köhler, J. P. Couperus, O. Zarini, A. Jochmann, A. Irman, and U. Schramm, Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 829, 265 (2016).
  * Schramm  _et al._ [2017] U. Schramm, M. Bussmann, A. Irman, M. Siebold, K. Zeil, D. Albach, C. Bernert, S. Bock, F. Brack, J. Branco, J. Couperus, T. Cowan, A. Debus, C. Eisenmann, M. Garten, R. Gebhardt, S. Grams, U. Helbig, A. Huebl, T. Kluge, A. Köhler, J. Krämer, S. Kraft, F. Kroll, M. Kuntzsch, U. Lehnert, M. Loeser, J. Metzkes, P. Michel, L. Obst, R. Pausch, M. Rehwald, R. Sauerbrey, H. Schlenvoigt, K. Steiniger, and O. Zarini, Journal of Physics: Conference Series 874, 012028 (2017).
  * Couperus  _et al._ [2015] J. P. Couperus, A. Köhler, T. A. W. Wolterink, A. Jochmann, O. Zarini, H. M. J. Bastiaens, K. J. Boller, A. Irman, and U. Schramm, Nuclear Instruments and Methods in Physics Research, Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 830, 504 (2015).
  * Zeng  _et al._ [2014] M. Zeng, M. Chen, Z. M. Sheng, W. B. Mori, and J. Zhang, Physics of Plasmas 21 (2014), 10.1063/1.4868404, 1312.5825 .
  * Mirzaie  _et al._ [2015] M. Mirzaie, S. Li, M. Zeng, N. A. M. Hafz, M. Chen, G. Y. Li, and Q. J. Zhu, Scientific reports , 1 (2015).
  * Kurz  _et al._ [2018] T. Kurz, J. P. Couperus, J. M. Krämer, H. Ding, S. Kuschel, A. Köhler, O. Zarini, D. Hollatz, D. Schinkel, R. D’Arcy, J.-P. Schwinkendorf, J. Osterhoff, A. Irman, U. Schramm, and S. Karsch, Review of Scientific Instruments 89, 093303 (2018).
  * Jochmann  _et al._ [2013] A. Jochmann, A. Irman, M. Bussmann, J. P. Couperus, T. E. Cowan, A. D. Debus, M. Kuntzsch, K. W. D. Ledingham, U. Lehnert, R. Sauerbrey, H. P. Schlenvoigt, D. Seipt, T. Stöhlker, D. B. Thorn, S. Trotsenko, A. Wagner, and U. Schramm, Phys. Rev. Lett. 111, 114803 (2013).
  * Krämer  _et al._ [2018] J. M. Krämer, A. Jochmann, M. Budde, M. Bussmann, J. P. Couperus, T. E. Cowan, A. Debus, A. Köhler, M. Kuntzsch, A. Laso García, U. Lehnert, P. Michel, R. Pausch, O. Zarini, U. Schramm, and A. Irman, Scientific Reports 8, 1398 (2018).
  * [44] “See supplemental material at [url will be inserted by publisher].” .
  * Michel  _et al._ [2006] P. Michel, C. B. Schroeder, B. A. Shadwick, E. Esarey, and W. P. Leemans, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 74, 1 (2006).
  * Curcio  _et al._ [2017] A. Curcio, M. Anania, F. Bisesto, E. Chiadroni, A. Cianchi, M. Ferrario, F. Filippi, D. Giulietti, A. Marocchino, M. Petrarca, V. Shpakov, and A. Zigler, Physical Review Accelerators and Beams 20, 012801 (2017).
  * Lu  _et al._ [2007] W. Lu, M. Tzoufras, C. Joshi, F. S. Tsung, W. B. Mori, J. Vieira, R. A. Fonseca, and L. O. Silva, Phys. Rev. Spec. Top. - Accel. Beams 10, 1 (2007), 0612227 [physics] .
  * Burau  _et al._ [2010] H. Burau, R. Widera, W. Honig, G. Juckeland, A. Debus, T. Kluge, U. Schramm, T. E. Cowan, R. Sauerbrey, and M. Bussmann, IEEE Transactions on Plasma Science 38, 2831 (2010).
  * Huebl  _et al._ [2018] A. Huebl, R. Widera, B. Worpitz, R. Pausch, H. Burau, M. Garten, S. Starke, A. Grund, A. Debus, A. Matthes, S. Bastrakov, K. Steiniger, I. Göthel, S. Rudat, J. Kelling, and M. Bussmann,  (2018), 10.5281/ZENODO.1478937.