# Estimating the effective lifetime of the $z\sim 6$ quasar population from
the
composite proximity zone profile

Karna A. Morey MIT Kavli Institute for Astrophysics and Space Research, 77
Massachusetts Ave., Cambridge, MA 02139, USA  Anna-Christina Eilers MIT Kavli
Institute for Astrophysics and Space Research, 77 Massachusetts Ave.,
Cambridge, MA 02139, USA  Frederick B. Davies Lawrence Berkeley National
Laboratory, CA 94720-8139, USA  Max-Planck-Institute for Astronomy, Königstuhl
17, 69117 Heidelberg, Germany  Joseph F. Hennawi Physics Department,
University of California, Santa Barbara, CA 93106-9530, USA  Robert A. Simcoe
MIT Kavli Institute for Astrophysics and Space Research, 77 Massachusetts
Ave., Cambridge, MA 02139, USA

###### Abstract

The lifetime of quasars can be estimated by means of their proximity zone
sizes, which are regions of enhanced flux bluewards of the Lyman-$\alpha$
emission line observed in the rest-frame UV spectra of high-redshift quasars,
because the intergalactic gas has a finite response time to the quasars’
radiation. We estimate the effective lifetime of the high-redshift quasar
population from the composite transmitted flux profile within the proximity
zone region of a sample of $15$ quasars at $5.8\leq z\leq 6.6$ with precise
systemic redshifts, and similar luminosities, i.e. $-27.6\leq
M_{1450}\leq-26.4$, and thus a similar instantaneous ionizing power. We
develop a Bayesian method to infer the effective lifetime from the composite
spectrum, including robust estimates of various sources of uncertainty on the
spectrum. We estimate an effective lifetime of the quasar population as a
whole of $\log_{10}(t_{\rm Q}/{\rm yr})=5.7^{+0.5(+0.8)}_{-0.3(-0.5)}$ given
by the median and $68$th ($95$th) percentile of the posterior probability
distribution. While our result is consistent with previous quasar lifetime
studies, it poses significant challenges on the current model for the growth
of supermassive black holes (SMBHs) located in the center of the quasars’ host
galaxies, which requires that quasar lifetimes are more than an order of
magnitude longer.

Supermassive black holes — Quasars— High-redshift galaxies—Early universe —
method: Astronomy data analysis

††journal: ApJ††thanks: NASA Hubble Fellow

##  1 Introduction

Quasars are powered by accretion onto a central supermassive black hole
(SMBH), additionally, they are the most luminous, non–transient objects in the
known universe and thus they can be detected and observed at very early cosmic
epochs (e.g. Mortlock et al., 2011; Venemans et al., 2013; Wu et al., 2015;
Mazzucchelli et al., 2017; Bañados et al., 2018; Onoue et al., 2019; Yang et
al., 2020; Wang et al., 2020). High–redshift quasars at $z\gtrsim 6$ are of
particular interest for studying the early formation and growth of SMBHs.
Observations show that their SMBHs have masses on the order of $\gtrsim 1$
billion solar masses already at a time when the universe is still very young,
i.e. $\lesssim 1$ Gyr after the Big Bang.

These black holes are believed to grow from an initial seed via accretion
during luminous quasar phases. The duration of the luminous quasar phase,
known as the quasar lifetime $t_{Q}$, is an important parameter to understand
the formation and growth of these SMBHs (Inayoshi et al., 2020). Assuming
exponential growth, i.e.

$\displaystyle M_{\rm BH}(t_{Q})=M_{\rm
seed}\exp\left(\frac{t_{Q}}{t_{S}}\right)$ (1)

with an e-folding or Salpeter time

$\displaystyle t_{S}\simeq 4.5\times
10^{7}\left(\frac{\epsilon}{0.1}\right)\left(\frac{L_{\rm bol}}{\bm{L}_{\rm
edd}}\right)^{-1}~{}\rm yr$ (2)

it requires at least 0.8 Gyr to grow a $10^{9}$ $M_{\odot}$ SMBH from a 100
$M_{\odot}$ seed with a fiducial radiative efficiency of 10% (Shakura &
Sunyaev, 1973) and accretion at the Eddington limit, i.e. $L_{\rm bol}=L_{\rm
edd}$. It remains an open question as to the validity of this simple
exponential growth model and the possibility of additional physics leading to
more complicated light curves (e.g. Di Matteo et al., 2005; Springel et al.,
2005; Hopkins et al., 2005; Novak et al., 2011; Davies et al., 2019a).

At high redshifts, the age of the universe is comparable to the timescale
required for SMBHs to form (Volonteri, 2010). However, measurements of the
timescales of quasar activity have proven to be difficult. At lower redshifts
of $z\sim 2-4$ previous studies have estimated the duty cycle of quasars,
which denotes the total fraction of time compared to the age of the universe
that galaxies spent as luminous quasars and therefore represents an upper
limit on the quasars’ lifetime $t_{Q}$, by comparing the number density of
quasars to the abundance of dark matter halos via clustering studies
(Efstathiou & Rees, 1988; Haiman & Hui, 2001; Martini & Weinberg, 2001;
Martini, 2004). These studies rely on the measurement of the abundance of
quasar host dark matter halos in comparison to the abundance of luminous
quasars to estimate the quasars’ duty cycle (White et al., 2008). However,
these methods have yielded only weakly constrained estimates of $\sim
10^{6}-10^{9}$ years, since clustering estimates are susceptible to high
uncertainty due to variance in parameters governing how quasars populate dark
matter halos (Shen et al., 2009; White et al., 2012; Conroy & White, 2013; Cen
& Safarzadeh, 2015). Other estimates of quasar activity timescales have been
made using an extension of the “Soltan argument” (Soltan, 1982), where the
quasar luminosity function can be compared to quiescent SMBHs observed in
galaxies locally, which have led to an estimate of the quasars’ duty cycle to
between $10^{7}-10^{8}$ years (Yu & Tremaine, 2002).

Furthermore, previous studies have used the quasars’ light echo transverse to
the line–of–sight to measure the time delay between the onset of the quasar
emission and changes in opacity of the intergalactic medium (Adelberger, 2004;
Hennawi & Prochaska, 2007; Schmidt et al., 2017, 2018; Bosman et al., 2020) to
constrain the lifetime of quasars to be $t_{Q}\sim 10^{5}-10^{7}$ years.

Recent studies have shown how the lifetime of quasars $t_{Q}$ can be inferred
by means of the proximity zone region observed in rest-frame UV and optical
spectra of high redshift quasars (Khrykin et al., 2016; Eilers et al., 2017,
2018; Khrykin et al., 2019; Davies et al., 2019b; Eilers et al., 2020). The
proximity zone is a region of enhanced transmitted flux in the vicinity of
quasars that has been ionized by the quasar’s own radiation (Bajtlik et al.,
1988; Haiman & Cen, 2001; Wyithe et al., 2005; Bolton & Haehnelt, 2007; Lidz
et al., 2007; Bolton et al., 2011; Keating et al., 2015). Estimating quasar
lifetimes from the HI proximity zone is a technique that has primarily been
used to estimate $t_{Q}$ of high redshift quasar populations, but recent
studies have used the Helium II line at $z\simeq 3-4$ to measure individual
quasar lifetimes and the quasar lifetime distribution (Khrykin et al., 2021;
Worseck et al., 2021; Khrykin et al., 2019). The intergalactic medium has a
finite response time to the quasar’s radiation and therefore the level of
enhanced flux due to ionized gas around the quasar is highly sensitive to the
lifetime of quasars.

Eilers et al. (2017, 2018) used proximity zones in the H I Ly$\alpha\,$forest
to estimate the lifetime for $z\simeq 6$ quasars. The same method can be
applied to proximity zones in the He II Ly$\alpha\,$forest to estimate the
lifetimes of quasars at low redshift ($z\simeq 3-4$), as was done in Khrykin
et al. (2019). Eilers et al. (2017) found a population of quasars $z\simeq 6$
that show very small proximity zones and therefore are likely to be very young
objects ($t_{Q}\lesssim 10^{5}\mathrm{~{}yr}$).

Assuming a simple “light-bulb” light curve for the quasars, i.e. the quasar
turns on abruptly and emits at a constant luminosity during its entire
lifetime $t_{Q}$, the probability of detecting a quasar age of $t_{Q}^{(i)}$
is $p=t_{Q}^{(i)}/t_{Q}$. Given that these young quasars corresponded to
$\simeq 10\%$ of the population, this suggested a fiducial lifetime of
$t_{Q}=10^{6}$ yr for the quasar population (Eilers et al., 2017, 2020).
Davies et al. (2019a) expanded upon this work by developing a semi-analytic
model of quasar proximity zones, including non-equilibrium interactions
between the IGM and the quasar, using hydrodynamical radiative transfer
simulations to simulate the proximity zone behavior on different timescales.
The observed distribution of proximity zones by Eilers et al. (2017) was
similar to the modeled distribution for $t_{Q}\simeq 10^{6}$ using the
“lightbulb” model. As part of this analysis, Davies et al. (2019a) considered
more general light curves, not just the “lightbulb” model, and found
consistent results for the lifetime, again confirming evidence for
$t_{Q}\simeq 10^{6}$ yr.

Davies et al. (2019b) measured the lifetime of individual high redshift
quasars during the epoch of reionization (EoR), and measured the age of these
quasars to be roughly $\simeq 10^{6}$ yr. By using the neutral IGM as a
counter of the ionizing photons emitted by the quasar, they were able to
constrain the lifetimes of two $z\gtrsim 7$ quasars. By measuring the number
of emitted ionizing photons, they were able to directly measure the radiative
efficiency $\epsilon\simeq 0.09\%$, much smaller than the canonical 10% value
given by Eddington limited accretion. This study notes, however, that if this
behavior is instead caused by UV obscured quasars, these radiative
efficiencies are consistent with a fraction of UV obscured quasars greater
than $82\%$ at redshift $z\simeq 7$. However, further work is needed to
characterize whether this mechanism is indeed radiative inefficiency or UV
obscuration.

These lifetime estimates are more than an order of magnitude shorter than
expected and cause tension with the current model for the growth of SMBHs. A
potential solution could be the presence of very massive initial black hole
seeds (i.e. $M_{\text{seed}}>1000M_{\odot}$), vastly exceeding the mass of
Population III stellar remnants (Lodato & Natarajan, 2006; Visbal et al.,
2014; Habouzit et al., 2016; Schauer et al., 2017). As mentioned above, very
low radiative efficiencies could also explain such short lifetimes, with
radiative efficiency rate $\epsilon\lesssim 0.01-0.001$ previously suggested
by other studies (Volonteri et al., 2015; Trakhtenbrot et al., 2017; Davies et
al., 2019b), or highly UV obscured black hole growth phases (Eilers et al.,
2018; Davies et al., 2019b). Additional estimates of $t_{Q}$, both for
individual quasars and for the entire population, could provide important
constraints on these mechanisms, i.e. provide constraints on the mass of the
initial seeds or the specific radiative efficiency rate $\epsilon$.

In this study, we will estimate the effective lifetime of the high-redshift
quasar population by means of the composite transmitted flux within the
proximity zone region. To this end, we assembled a data set of $15$ quasar
spectra at $z\sim 6$ with a similar absolute magnitude, i.e. a similar
ionizing radiation output. By comparing this stacked proximity zone profile to
outputs from radiative transfer simulations of quasars at different lifetimes,
we constrain the effective lifetime of the quasar population.

This paper is organized as follows: in § 2 we discuss our quasar sample and
the selection criteria. In § 3 we describe our quasar continuum estimation
procedure and the radiative transfer models. In § 4, we discuss our stacking
procedure, as well as the forward modeling of uncertainties, and present our
final result. We summarize our results in § 5.

Throughout this paper, we assume a flat $\Lambda$CDM cosmology of $h=0.685$,
$\Omega_{m}=0.3$, and $\Omega_{\Lambda}=0.7$, which is consistent within the
$1\sigma$ errorbars from Planck Collaboration (2020).

Table 1: Overview of our data sample.

Name | RA [hms] | DEC [dms] | Redshift | $\mathrm{M}_{1450}$ | Telescope/Instrument | SNR | PID | PI
---|---|---|---|---|---|---|---|---
PSOJ011+09 | 00:45:33.566 | +09:01:56.93 | 6.4693 | -26.85 | VLT/X-Shooter | 7 | 0101.B-0272(A) | Eilers
PSOJ036+03 | 02:26:01.873 | +03:02:59.24 | 6.5412 | -27.28 | VLT/X-Shooter | 11 | 0102.A-0154(A)/ 0100.A-0625(A) | D’Odorico
PSOJ065-26 | 04:21:38.050 | -26:57:15.72 | 6.1877 | -27.25 | VLT/X-Shooter | 19 | 098.B-0537(A) | Farina
SDSSJ0842+1218 | 08:42:29.438 | +12:18:50.47 | 6.0763 | -26.91 | VLT/X-Shooter | 26 | 097.B-1070(A) | Farina
ATLAS/PSOJ158-14 | 10:34:46.509 | -14:25:15.86 | 6.0681 | -27.41 | VLT/X-Shooter | 15 | 096.A-0418(B) | Shanks
PSOJ159-02 | 10:36:54.190 | -02:32:37.940 | 6.3809 | -26.59 | VLT/X-Shooter | 9 | 098.B-0537(A) | Farina
SDSSJ1143+3803 | 11:43:38.347 | +38:08:28.823 | 5.8367 | -26.69 | Keck/DEIMOS | 15 | 2017A_U078 | Hennawi
SDSSJ1148+5251 | 11:48:16.64 | +52:51:50.2 | 6.4189 | -26.81 | Keck/ESI | 37 | U43E, U11E, U2E | Becker
|  |  |  |  |  |  | C70E | Djorgovski
|  |  |  |  |  |  | H25aE, H45aE | Cowie
|  |  |  |  |  |  | C97E | Kulkarni
SDSSJ1306+0356 | 13:06:08.259 | +03:56:26.19 | 6.0337 | -26.81 | VLT/X-Shooter | 38 | 084.A-0390(A) | Ryan-Weber
ULASJ1319+0950 | 13:19:11.291 | +09:50:51.49 | 6.1330 | -26.88 | VLT/X-Shooter | 52 | 084.A-0390(A) | Ryan-Weber
SDSSJ1335+3533 | 13:35:50.81 | +35:33:15.82 | 5.9012 | -26.67 | Keck/DEIMOS | 22 | 2017A_U078 | Hennawi
PSOJ217-16 | 14:28:21.371 | -16:02:43.73 | 6.1498 | -26.93 | Magellan/FIRE | 9 |  | Simcoe
CFHQSJ1509-1749 | 15:09:41.781 | -17:49:26.68 | 6.1225 | -27.14 | VLT/X-Shooter | 37 | 085.A-0299(A) | D’Odorico
PSOJ323+12 | 21:32:33.178 | +12:17:55.07 | 6.5872 | -27.06 | VLT/X-Shooter | 14 | 098.B-0537(A) | Farina
PSOJ359-06 | 23:56:32.452 | -06:22:59.31 | 6.1719 | -26.79 | VLT/X-Shooter | 19 | 098.B-0537(A) | Farina

Note. — The sub-mm redshift estimate for PSOJ036+03 was obtained by Venemans
et al. (2015), the redshift estimate for J1319+0950 was obtained by Wang et
al. (2013), the redshift estimate for SDSSJ1335+3533 was obtained by Wang et
al. (2013), the redshift estimates for SDSSJ1143+3808, PSOJ359-06, PSOJ158-14,
and PSOJ011+09 were obtained by Eilers et al. (2020), and all others were
published in Decarli et al. (2018). The SNR values are estimated using a pixel
scale of 0.05 Å/pixel in the rest frame of the quasars.

##  2 Data Sample

For our analysis we select quasars at $5.8<z<6.6$, and an absolute magnitude
at 1450 Å in the rest-frame within a narrow range of $-26.6\lesssim
M_{1450}\lesssim-27.4$. We require that all objects in our sample have
available medium-resolution spectroscopy at optical and possibly also near-
infrared (NIR) wavelengths. Additionally, we select quasars for which precise
redshift estimates ($\Delta v\lesssim 100\,\rm km\,s^{-1}$) based on sub-mm
observations of the [C II] emission line at $158\,\mu$m ($\nu_{\rm
rest}=1900.548$ GHz) or the CO(6–5) emission line at $3\,\rm mm$ ($\nu_{\rm
rest}=691.473\rm\,GHz$) from the Atacama Large Millimetre Array (ALMA) or the
NOrthern Extended Millimeter Array (NOEMA) at the Institute de Radioastronomie
Millimétrique (IRAM) are available, since uncertainties on the quasars’
systemic redshift contribute the largest source of uncertainty to the
proximity zone measurements (Eilers et al., 2017, 2020). Furthermore, we
remove all quasars from the sample that show broad absorption line features
(BALs) or proximate damped Ly$\alpha\,$ systems (pDLAs), which would cause
additional absorption in the proximity zone region and would bias our results.
These selection criteria result in $15$ objects with a mean redshift of
$\langle z\rangle=6.23$, and mean magnitude of $\langle
M_{1450}\rangle=-26.9$. An overview of our data sample is shown in Table 1.

Of the quasars in our sample, eleven were observed with the X-Shooter
instrument at the Very Large Telescope (VLT) (Vernet et al., 2011), two where
observed using the DEep Imaging Multi-Object Spectrograph (DEIMOS) instrument
at the Keck observatory Faber et al. (2003), one was observed using the
Echellette Spectrograph and Imager (ESI) instrument also at the Keck
observatory (Sheinis et al., 2002), and one was observed using the Folded-port
InfraRed Echellette (FIRE) spectrometer instrument at the Magellan telescope
(Simcoe et al., 2013). The individual detector resolutions
$R=\frac{\lambda}{\Delta\lambda}$ are $R\approx 8800$ for the X-Shooter
instrument Vernet et al. (2011), $R\approx 6000$ for the FIRE
instrument(Simcoe et al., 2013), $R\approx 5000$ for the DEIMOS instrument
(Faber et al., 2003), and $R\approx 5400$ for the ESI instrument (Sheinis et
al., 2002).

###  2.1 Data Reduction

All optical and NIR spectroscopic data observed with VLT/X-Shooter and
Keck/DEIMOS were reduced applying the open source python spectroscopic data
reduction package PypeIt (Prochaska et al., 2020). For the reduction we first
subtract the sky emission, which was performed on the 2D images by including
both image differencing between dithered exposures whenever these were
available, and a B-spline fitting procedure (e.g. Bochanski et al., 2009). We
apply an optimal spectrum extraction technique (Horne, 1986) to extract the 1D
spectra. All individual 1D spectra are flux calibrated using the standard
stars $\rm LTT\,3218$ (for spectra observed with VLT/X-Shooter) and $\rm
G191B2B$ or $\rm Feige\,110$ (for spectra taken with Keck/DEIMOS). These
fluxed 1D spectra are stacked and a telluric model and the quasar PCA model
from Davies et al. (2018a) are jointly fitted to the stacked spectra using
telluric model grids from the Line-By-Line Radiative Transfer Model
(LBLRTM111http://rtweb.aer.com/lblrtm.html; Clough et al., 2005; Gullikson et
al., 2014).

Figure 1: Quasar spectra from our sample. We mask the regions of telluric
absorption in grey. All fluxes (black) and noise vectors (grey) have been
inverse-variance smoothed with a 10 pixel filter. We also show the quasar
continuum estimates fitted on the red-side (red) and the blue-side prediction
(blue) from our PCA analysis described in section 3.1. Figure 2: Same as Fig.
1 for the remaining objects in our sample.

We then match the visible (VIS) and NIR regions for the VLT/X-Shooter spectra.
To this end, we compute the mean value of all pixels between 10070 and 10270 Å
in the observed frame, and then scale the NIR spectra according to a constant
factor such that the two means are equal. In this step, we do not normalize
the spectra to unity, but merely scale the NIR spectra to be matched with the
VIS spectra.

For the FIRE spectra, which do not require matching between NIR and visible
parts, we use the IDL-based pipeline
FIREHOSE222http://web.mit.edu/r̃simcoe/www/FIRE/ob_data.htm, developed by
Simcoe et al. (2010). The procedures used for data reduction include cosmic-
ray rejection, flat-fielding, wavelength calibration using OH sky lines
(calibrated to vacuum wavelength), nonlinearity correction, optimal extraction
of one-dimensional spectra Horne (1986), combining individual exposures,
merging multiple echelle orders, and heliocentric corrections. This pipeline
performs sky subtraction via a exposure-specific B-spline model of the sky,
based on the technique developed in Bochanski et al. (2009). We reduce the ESI
spectra using the ESIRedux
pipeline333http://www2.keck.hawaii.edu/inst/esi/ESIRedux/ developed as part of
the XIDL444http://www.ucolick.org/~xavier/IDL// suite of astronomical routines
in the Interactive Data Language (IDL). This pipeline is very similar to
FIREHOSE. More details about modifications to this pipeline for high-redshift
quasars can be found in Eilers et al. (2017).

In the end we mask all spectral regions that are affected by telluric
absorption, namely between 13000 and 15200 Å, 17500 and 20000 Å, and beyond
22500 Å in the observed frame. Figures 1 and 2 show the final spectra of all
$15$ quasars in our sample. We estimate the SNR of the spectra by taking the
median pixel-based ratio of the signal to noise between 1265 and 1295 Å in the
rest frame of the quasar, due to the absence of any emission lines in this
region. We scale the pixel-based SNR for each of the instruments assuming
Gaussian white noise to the X-shooter resolution to compare the SNR
measurements with a common resolution in Table 1.

##  3 Methods

We aim to estimate the lifetime of the quasar population as a whole by means
of their composite proximity zone region. Thus, in order to obtain the stacked
transmission profile within the quasars’ proximity zones, we first have to
normalize the spectra by their continuum emission. Once the quasars are
continuum normalized, we stack them between 1190 and 1218 Å in the rest–frame,
in order to obtain an estimate for the composite flux transmission profile in
the proximity zone region. After stacking the quasars, we compare this stacked
flux profile to 1-d radiative transfer models.

###  3.1 Quasar Continuum Normalization

Figure 3: Continuum normalized spectra (black) along with continuum normalized
noise vectors (grey) within our sample. The continuum normalized spectra are
shown between 1177 and 1222 Å, except for PSOJ217-16, where we only show
between 1190 and 1222 Å, due to detector noise. All spectra have been inverse
variance smoothed with a 3 pixel filter. The grey dashed line marks the
Ly$\alpha\,$line at 1215.67 Å.

To estimate the continuum of the quasar within the proximity zone region, we
fit for the quasar continuum redwards of Ly$\alpha\,$, as this side is
unaffected by IGM absorption. We then use this red-side estimate to predict
the blue-side continuum, which is affected by absorption (Suzuki et al., 2005;
Paris et al., 2011; Davies et al., 2018a; Bosman et al., 2021).

As a first step to estimate the quasar continua we normalize the matched
spectra to unity at $1290\pm 2.5$ Å, for ease of fitting. Then, we sigma clip
outlying values at the $3\sigma$ level from the smoothed spectra (10 pixel
filter), shown in Figures 1 and 2, using the noise vectors from each spectra
as an estimate of the $\sigma$. To estimate the quasar continua, we make use
of the principal component analysis (PCA) method developed by Davies et al.
(2018a), where we use 10 red-side basis spectra (spanning between 1200 and
2850Å for most quasars in the sample) and 6 blue-side basis spectra (spanning
between 1177 and 1280 Å). For SDSSJ1143, SDSSJ1148 and SDSSJ1335 we use 10
truncated red-side basis spectra, instead spanning between 1200 and 1450 Å to
accommodate these quasars where we do not have NIR coverage.

The principal components were determined on a sample of $12,764$ low-redshift
quasars at $z\approx 2$, for which the blue-side continuum is easier to
estimate due to the lower level of absorption by the intergalactic medium.
Davies et al. (2018a) noticed a high degree of correlation or anti-correlation
between the red-side coefficients and blue-side coefficients, allowing the
blue side coefficients to be estimated from the red side. This can be done
using the relationship $b_{i}=\sum_{j=1}^{10}r_{j}X_{ji}$, where the $b_{i}$
and $r_{j}$ are blue side and red side PCA coefficients, respectively, and
$X_{ji}$ is a projection matrix determined by the correlations or anti-
correlations between the blue and red side coefficients from the training
sample.

Using the MCMC affine-invariant ensemble sampler `emcee` (Foreman-Mackey et
al., 2013), we fit each spectrum using $10$ basis spectra redwards of the
Ly$\alpha\,$ emission line spanning between 1220 Å and 2800 Å in the rest
frame, as well as a redshift offset $\delta z$ between the systemic redshift
of the quasar and the best redshift estimate for the PCA. We first take the
median of the posterior distribution to estimate the median red side
coefficients, we then use the projection matrix $X_{ji}$ to estimate the
coefficients for the blue-side continuum. However, unlike Davies et al.
(2018a) who developed this machinery for studying IGM damping wings in $z>7$
quasar spectra (Davies et al., 2018b), we match the blue and red side basis
spectra at $1220$ Å instead of at $1280$ Å (Eilers et al., 2020). This is due
to the fact that all quasars in our sample are at $z\lesssim 6.6$, and thus we
do not expect significant absorption redwards of the Ly$\alpha\,$ emission
line due to IGM damping wings. We set a flat prior on the PCA components to be
the $3\sigma$ boundaries on each of the components determined by Davies et al.
(2018a), and allow a prior $\delta z\in[-0.01,0.03]$ for all quasars except
for PSO J323+12, where we extend the prior to $\delta z\in[-0.05,0.05]$ in
order to allow for a better continuum fit. The quasar spectra and their
estimated continua are shown in Fig. 1 and 2. The continuum normalized spectra
that are obtained when dividing by the estimated continuum emission are shown
for each quasar in Fig. 3.

The quasar continuum normalization procedure introduces additional uncertainty
into our analysis procedure due to the fact that the continuum estimation can
either be an overestimate or underestimate of the true continuum. It can be
seen that the PCA based method underestimates the continuum for SDSSJ1143+3803
and CFHQSJ1509-1749 since the continuum normalized flux is clearly above a
flux of unity close to the Ly$\alpha\,$ line. However, the estimated bias in
the PCA continuum resconstruction is very small, i.e. $<1\%$ (Davies et al.,
2018a; Bosman et al., 2021) and thus the this underestimation is balanced by
the roughly symmetric overestimation of other quasar continua, which can be
less obviously identified as it decreases the continuum normalized flux near
Ly$\alpha$. We address this uncertainty in the quasar continuum estimates in
detail in section 4.1.

###  3.2 Radiative Transfer Simulations

For each quasar in our sample we use a 1-d radiative transfer (RT) code
(Davies et al., 2016) to simulate the effect of ionizing radiation emitted by
the quasar along the line-of-sight (Bolton & Haehnelt, 2007). We take $1000$
skewers originating from the most massive dark matter halos with halo masses
of $4\times 10^{11}\,M_{\odot}\lesssim M_{\rm halo}\lesssim 3\times
10^{12}\,M_{\odot}$ from the cosmological hydrodynamical simulation Nyx with a
box of $100\,{\rm Mpc}\,h^{-1}$ on a side (Almgren et al., 2013; Lukić et al.,
2014). The RT code computes the time–dependent ionization and recombination
rates of six particle species ($e^{-}$, H I, H II, He I, He II, He III), as
well as the heating and cooling of the IGM due to the expansion of the
universe and inverse Compton cooling off the cosmic microwave background.

The Nyx simulation has outputs at $z=5.5$, $6.0$ and $6.5$. We choose the
output closest to the quasar’s emission redshift $z_{\rm em}$ and re-scale the
physical densities by $(1+z_{\rm em})^{3}$. We used the Lusso et al. (2015)
composite composite to relate $M_{1450}$ to the specific flux at the Lyman
limit, and then compute the ionizing flux by assuming a spectral slope of
$\alpha_{\nu}=-1.70$. Previous studies of quasar lifetimes based on proximity
zones have only considered the extent of the proximity zone (with the notable
exception of Davies et al. (2018b, 2019a)), which is defined as the location
at which the continuum-normalized flux smoothed to a resolution of $20$ Å in
the observed frame drops below the 10% flux transmission level (Fan et al.,
2006). Therefore, the precise ionization rate of the ultraviolet background
(UVB) did not play a crucial role, since the total ionization rate was always
dominated by the quasar at the 10% flux transmission level. In this work,
however, we consider the complete flux profile within the vicinity of the
quasar, which extends below the 10% flux transmission level. At these lower
flux transmission levels the ionization rate of the UVB is comparable to the
ionization rate of the quasar and therefore needs to be considered, although
the UVB should not matter too much given that we are assuming uniformity.
Nonetheless, we re-scale the ionization rate of the UVB radiation $\Gamma_{\rm
UVB}$ to match the expected strength at the quasar’s redshift based on Haardt
& Madau (2012).

We then simulate the composite quasar spectrum as a function of quasar
lifetime, assuming a light-bulb light curve, ranging from $\log_{10}(t_{Q}/\rm
yr)=1-8.9$ in steps of $\delta\log_{10}(t_{Q}/\rm yr)=0.1$. Each simulated
spectrum is tuned to that quasar’s particular redshift and luminosity, and the
UVB was calibrated at the particular redshift to match Haardt & Madau (2012).
For each step in the simulation, we compute 1000 different spectra at
different sight lines, i.e. we have a model grid of $15\times 1000\times 80$
different spectra, since there are 15 different quasars, 1000 different
sightlines, and 80 different values of the lifetime, from 1.0 to 8.9 in $\log
t_{Q}$ space. At each time step we stack the simulated spectra by averaging
the flux pixels across all 15 quasars and across all 1000 sight lines, since
each simulation is computed on the same velocity grid. The simulated stacked
spectra for a few selected lifetimes are shown in the left panel of Fig. 4.

Figure 4: Left panel: Proximity zone profile as predicted by radiative
transfer simulations for $\log_{10}(t_{Q}/\rm yr)=$2, 4, 6, and 8. The shaded
regions around each of the mean models show the $1\sigma$ uncertainty range as
described in section 4.2. The grey dashed line marks the Ly$\alpha\,$line at
1215.67 Å. Right panel: Plot of the quasar stack, with a bin size of 0.5 Å,
with the best fit model with $t_{Q}=10^{5.7}$ yr shown in purple.

##  4 Results

###  4.1 Stack of the Quasars’ Proximity Zones

We now stack the continuum–normalized fluxes within the proximity zone region
in order to create a composite flux transmission profile. To this end, we
define a regularly spaced grid between 1190 Å and 1218 Å with a grid spacing
of $0.5$ Å. We then bin the wavelengths and the flux pixels within this new
wavelength grid, i.e. we take an average of all the wavelength pixels and flux
pixels within a particular wavelength bin to obtain the stacked flux profile.
The full composite spectrum within the proximity zone region from the $15$
quasars in our data sample is shown in the right panel of Fig. 4.

###  4.2 Forward-Modeling of Uncertainties

There are two possibilities to estimate uncertainties on the flux pixels
within the stacked proximity zone, i.e. either via bootstrapping of the
continuum-normalized spectra, or via forward modeling the various sources of
uncertainty onto the simulated quasar spectra. Forward modeling the sources of
uncertainty does not give us an estimate of the covariance directly from the
data, rather we use covariance estimated from the forward models in performing
inference of the quasar lifetime $t_{Q}$. This gives a better estimate of the
covariance than estimating it directly from the data, with the only drawback
that the uncertainty on the stack is now model-dependent. Due to the limited
size of our data sample, the covariance matrix obtained via bootstrapping is
very noisy, and thus we forward–model the uncertainties onto the simulated
quasar spectra from the RT simulations by means of forward modeling. We then
estimate the covariance between the different flux bins on the sample of
$10,000$ mock stacks with uncertainties. For each stack, we incorporated four
types of uncertainty: cosmic variance, uncertainty from the PCA continuum
estimation, uncertainty in the systemic redshift, and spectral noise. We
assume that the full uncertainties on the stack can be approximated as a
multivariate Gaussian distribution about a mean model. We will now describe in
detail how these uncertainties are being modeled in the quasar stack. We
repeat this uncertainty analysis for each lifetime, so we obtain $80$
different samples of $10,000$ mock spectra. We then assume that these mock
stacks are from a Gaussian distribution and then estimate the covariance of
this distribution directly using the mock stacks.

The first source of uncertainty that we incorporated is uncertainty from
cosmic variance. Thus, for each quasar in the stack we randomly chose one of
the $1000$ sightlines through the cosmological simulation box, ensuring that
every simulated stack consists of $15$ different sightlines. In order to mimic
the spectral resolution of the data, we convolve each spectrum with a Gaussian
kernel with a standard deviation in velocity space given by
$\sigma_{\text{vel}}=\frac{c}{2.355R}$, where $c$ is the speed of light, $R$
is the resolution of the spectrograph, and the factor of $2.355$ arises from
the full width at half maximum (FWHM) to standard deviation conversion.
Following this, we interpolate the convolved spectrum onto the wavelength grid
of the data for the specific quasar in question. That is, we aim to generate
simulated data for each quasar that has the same wavelength grid as the actual
data, that way we can use the noise from the data to simulate spectral noise.

The Davies et al. (2018a) procedure produces a covariance matrix for the
continuum prediction, which we can directly sample from to create mock
realizations of the uncertainty in the continuum estimate. We then scale the
continuum-normalized simulated spectra for each quasar by a draw from this
distribution. In order to incorporate uncertainties in the systemic redshift
estimate of the quasars, we assume an uncertainty on the systemic redshift of
$\Delta v=100\,\rm km\,s^{-1}$, and randomly draw from a normal distribution
with a mean of $\mu=0$ and the corresponding standard deviation (Eilers et
al., 2020). The simulated quasar spectra are then shifted to the rest-frame
wavelength grid based on this new redshift estimate. Finally, in order to
incorporate spectral noise, we directly sample the noise vectors from the
data, as we have already interpolated the simulated convolved spectra onto the
wavelength grid. To be clear, we are generating sets of 15 quasars each on the
wavelength grid of the data, allowing us to directly sample the noise without
re-scaling, after which we can stack the 15 quasars within a set using our
binning algorithm.

Because the radiative transfer models are dependent on the quasar lifetime, as
an input parameter, we performed this Monte-Carlo sampling of mock stacks for
each value of the lifetime on a grid between 1.0 and 8.9 in $\log_{10}$ years,
with a grid spacing of $\Delta\log_{10}t_{Q}=0.1$. This means that we computed
80 different covariance matrices, for each lifetime dependent model, which we
used when we constructed our likelihood function over the log lifetime grid.
For each lifetime dependent model we use the same set of random draws for the
Monte Carlo sampling to ensure that the likelihood function varies smoothly
with the model parameter, rather than being influenced by the stochasticity of
the draws. We show an example correlation matrix for the model corresponding
to $t_{Q}=10^{5.7}$ yr in Fig. 5.

Figure 5: Correlation matrix for $\log_{10}(t_{Q}/{\rm yr})=5.7$ computed from
Monte-Carlo sampling the distribution of quasar stacks directly from the
simulated quasar spectra from radiative transfer models.

###  4.3 Quasar Lifetime Estimate

We assume that the data is drawn from a distribution described by one of the
MC estimated covariance matrices, centered at the mean model for that
particular lifetime. Additionally, we also assume that each quasar in the
stack has the same lifetime, i.e. our models are made up of stacks of quasars
of the same lifetime (see § 5 for a discussion of this assumption).

Given our assumptions, we construct a lifetime dependent likelihood function

$\begin{split}\mathscr{L}(r_{\lambda}|t_{Q})=&(2\pi)^{-\frac{k}{2}}\text{det}(\Sigma(t_{Q}))^{-1/2}\cdot\exp\bigg{(}-\frac{1}{2}(f_{\lambda,\text{data}}-f_{\lambda,\text{model}}(t_{Q}))^{T}\Sigma(t_{Q})^{-1}(f_{\lambda,\text{data}}-f_{\lambda,\text{model}}(t_{Q}))\bigg{)},\end{split}$
(3)

where $\mathscr{L}(r_{\lambda}|t_{Q})$ is the likelihood of observing the data
given a particular model $t_{Q}$, $k$ is the number of wavelength bins in the
stacked spectrum, and $\Sigma(t_{Q})$ is the model dependent covariance
matrix. However, we only know $r_{\lambda}(t_{Q})$ and $\Sigma(t_{Q})$ at
discrete points in log-lifetime space, since we have outputs of the RT
simulation only at discrete time steps. We calculated the likelihood function
at the discrete points within the $\log_{10}(t_{Q})$ space (i.e. where our
covariance matrices are defined) and interpolated the likelihood function to
estimate the likelihood between the grid-points.

To construct our likelihood we consider wavelengths from $1190$ to $1217$ Å in
the rest-frame. This cutoff at $1190$ Å avoids potential biases in the very
low flux regime, where the flux transmission is completely dominated by the
UVB.

To obtain the posterior probability density function (PDF) for the lifetime,
we use a flat prior for $t_{Q}\in[10,10^{8.9}]$ yr, meaning that the posterior
is merely the normalized version of the likelihood. We normalize the
likelihood function using trapezoidal rule numerical integration to obtain the
result shown in Fig. 6. We obtain a measurement of the effective lifetime of
the quasar population of $\log_{10}(t_{Q}/{\rm
yr})=5.7^{+0.5(+0.8)}_{-0.3(-0.5)}$ from the median and 68th (95th) percentile
of the posterior probability distribution. This best fit model is shown in the
right panel of Fig. 4.

Figure 6: Posterior PDF of the average quasar lifetime $t_{Q}$. Vertical
dashed lines indicate the median (orange), 68% credible interval (red), and
95% credible interval (blue). Figure 7: Modeled stacks of the 15 quasars in
the stack for a mean lifetime of $\log_{10}(t_{Q}/{\rm yr})=5.7$, but for
quasar lifetime distribution widths $\sigma=0,1,2$, along with $1\sigma$
uncertainties on the stack.

##  5 Summary & Discussion

In this paper we estimate the effective lifetime of the quasar population at
$z\sim 6$ by means of the average flux transmission profile in the proximity
zone region around the quasars. To this end, we stack the continuum-normalized
spectra of $15$ quasars at $5.8<z<6.6$ with similar absolute luminosities,
i.e. $-27.4<M_{1450}<-26.6$. For all objects in this data set precise systemic
redshift estimates from sub-mm emission lines are available, since
uncertainties in the systemic redshifts contribute the largest source of
uncertainty to proximity zone estimates (Eilers et al., 2017). This average
composite transmission profile is then compared to simulated transmission
profiles from 1-d RT simulations at different effective quasar lifetimes.

Our best estimate for the effective lifetime of the quasar population as a
whole based on the composite proximity zone profile is $\log_{10}(t_{Q}/{\rm
yr})=5.7^{+0.5(+0.8)}_{-0.3(-0.5)}$, derived from the median of the posterior
probability distribution as well as the $68$th ($95$th) percentile. In order
to estimate the uncertainties on the stacked spectrum we forward-model the
uncertainties arising due to cosmic variance, the PCA continuum-normalization,
the estimate of the quasars’ systemic redshifts, as well as spectral noise
onto the simulated quasar sightlines.

Our measurement of the effective lifetime for the quasar population is
consistent with previous work, which had estimated an average quasar lifetime
of $t_{Q}\sim 10^{6}$ yr (Khrykin et al., 2019; Davies et al., 2019a). Our
results are also consistent with the observed fraction of young quasars within
the early universe with lifetimes of $t_{Q}\lesssim 10^{4}-10^{5}$ yr, which
was determined to be $5\%\lesssim f_{\text{young}}\lesssim 10\%$ (Eilers et
al., 2020, 2021). This is because we expect to see approximately 10% of
quasars to shine for 10% of the lifetime, if the quasars’ ages are uniformly
sampled from $[0,t_{Q}]$.

In reality it is likely that the lifetime of the quasar population can not be
described by a single value, but rather by a distribution of lifetimes with a
mean of $\mu_{t_{Q}}=t_{Q}$ and standard deviation $\sigma_{t_{Q}}$. Because
we do not measure, or even consider the full distribution of quasar lifetimes,
we note that our measurement of the effective lifetime will be comparable to
the average lifetime only in the regime where the variance of the quasar
lifetime distribution is small. Khrykin et al. (2016) showed that a large
variance of the quasar lifetime distribution causes an underestimate of the
lifetime, because quasars with lifetime less than $t_{Q}$ impact the proximity
zone more. Additionally, if we have a very young quasar in our stack, this
could cause our estimate to be an underestimate of the lifetime. This is due
to the fact that the intergalactic gas reaches ionization equilibrium with the
quasar’s radiation for lifetimes longer than the equilibration timescale, i.e.
$t_{Q}>t_{\rm eq}\approx 3\times 10^{4}$ yr. This means that our measurement
of the effective lifetime of the quasar population could potentially
underestimate the mean quasar lifetime, if the width of the distribution is
very wide. In Figure 7 we show modeled stacks of the same mean lifetime
($\log_{10}(t_{Q}/\rm yr)=5.7$) for the overall log-normal distribution of
quasar lifetimes, but with different standard deviations. As we can see our
method is not really sensitive to the width of the quasar lifetime
distribution, unless the width of the lifetime distribution is large, i.e.
$\sigma_{\log_{10}t_{Q}}\gtrsim 2$, in which case our method would
underestimate the average quasar lifetime. However, recent results suggest
only a very modest standard deviation of the distribution of quasar lifetimes,
i.e. $\sigma_{\log_{10}t_{\rm Q}}=0.80^{+0.37}_{-0.27}$ (Khrykin et al.,
2021), suggesting that our measurement of the effective quasar lifetime
actually reflects the average lifetime of the quasar population as a whole.

Another potential source of uncertainty in our measurement is the systematic
bias in the PCA continuum estimation. Davies et al. (2018a) estimate that the
PCA reconstructed continuum overestimates the true quasar continuum on average
by $\sim 1\%$. However, this is a relative error on the flux, i.e. at $\sim
10\%$ flux transmission the estimated error is $0.1\%$ and hence this bias has
only negligible effects on our final results.

Our lifetime measurement is more than an order of magnitude shorter than
expected from exponential growth models for SMBHs (Volonteri, 2010, 2012). A
potential explanation for such short average quasar lifetimes could be
flickering quasar light curves instead of the simple “light-bulb” light
curves, implying that there could be multiple epochs of quasar activity and
concurrent black hole growth. Since the transmitted flux within the proximity
zone around quasars at $z\sim 6$ is only sensitive to the last “on-time” our
composite proximity zone spectrum from this analysis is not sensitive to
potential previous phases of quasar activity.

However, several analysis have shown previously that even flickering quasar
light curves cannot explain the discrepancy between the short average quasar
lifetimes and the time required to explain the SMBH growth with the current
exponential growth model. For instance, Davies et al. (2019b) showed that the
damping wing feature observed in quasar spectra at $z>7$, where the
surrounding IGM still has a very high neutral gas fraction, provides a
constraint on the _total_ number of ionizing UV photons emitted into the IGM
integrated over cosmic time, irrespective of the light-curve shape. Their
analysis shows that the integrated time these quasars have been emitting UV
radiation is still consistent with a light bulb light-curve with $t_{Q}\sim
10^{6}$ years.

At $z\sim 6$, however, the IGM is highly ionized and the quasar spectra do not
exhibit a damping wing. Hence the proximity zones are sensitive only to the
on-time of the last accretion episode, which is only consistent with the total
quasar lifetime for light-bulb light curves. (Davies et al., 2019a) has
analyzed the distribution of $z\sim 6$ proximity zones in the context of
various quasar light curves. While they find a good agreement between the
observed distribution of proximity zones and model predictions when assuming a
simple light-bulb light curve, flickering quasar light curves with multiple
black hole growth phases cannot be excluded for $z\sim 6$ quasars.

Our presented results at $z\sim 6$, the distribution of proximity zones at
$z\sim 6$ (Davies et al., 2019a), the fraction of very young quasars within
the quasar population (Eilers et al., 2020, 2021), as well as the damping wing
analysis at $z>7$ Davies et al. (2019b), are all in good agreement and point
towards quasar lifetimes of $\log_{10}(t_{Q})\simeq 10^{6}$. These
measurements suggests that either a significant fraction of the black hole
growth occurs in UV obscured, dust-enshrouded growth phases in the early
universe, or it could indicate radiatively inefficient mass accretion rates,
such that more of the accreted mass contributes to the growth of the black
holes rather than the quasars’ UV emission (Eilers et al., 2018; Davies et
al., 2019b). In future work, we aim to explore possible mechanisms by which
quasars can accrete sufficient mass while only emitting UV radiation for $\sim
10^{6}$ years on average. Additionally, we will expand this work to measure
the properties of the full quasar lifetime distribution from the individual
measurements, as suggested in (Khrykin et al., 2021). These studies will
provide further insights for our understanding of SMBH formation in the early
universe and whether or how the accretion mechanisms change with cosmic time.

The authors would like to thank the UROP office at MIT for generous funding of
this work.
ACE acknowledges support by NASA through the NASA Hubble Fellowship grant
$\\#$HF2-51434 awarded by the Space Telescope Science Institute, which is
operated by the Association of Universities for Research in Astronomy, Inc.,
for NASA, under contract NAS5-26555.

## References

  * Adelberger (2004) Adelberger, K. L. 2004, The Astrophysical Journal, 612, 706, doi: 10.1086/422804
  * Almgren et al. (2013) Almgren, A. S., Bell, J. B., Lijewski, M. J., Lukić, Z., & Van Andel, E. 2013, The Astrophysical Journal, 765, 39, doi: 10.1088/0004-637X/765/1/39
  * Bajtlik et al. (1988) Bajtlik, S., Duncan, R. C., & Ostriker, J. P. 1988, Astrophysical Journal, 327, 570
  * Bañados et al. (2018) Bañados, E., Venemans, B. P., Mazzucchelli, C., et al. 2018, Nature, 553, 473, doi: 10.1038/nature25180
  * Bochanski et al. (2009) Bochanski, J. J., Hennawi, J. F., Simcoe, R. A., et al. 2009, Publications of the Astronomical Society of the Pacific, 121, 1409, doi: 10.1086/648597
  * Bolton & Haehnelt (2007) Bolton, J. S., & Haehnelt, M. G. 2007, Monthly Notices of the Royal Astronomical Society: Letters, 381, doi: 10.1111/j.1745-3933.2007.00361.x
  * Bolton et al. (2011) Bolton, J. S., Haehnelt, M. G., Warren, S. J., et al. 2011, Monthly Notices of the Royal Astronomical Society: Letters, 416, 70, doi: 10.1111/j.1745-3933.2011.01100.x
  * Bosman et al. (2020) Bosman, S. E. I., Kakiichi, K., Meyer, R. A., et al. 2020, The Astrophysical Journal, 896, 49, doi: 10.3847/1538-4357/ab85cd
  * Bosman et al. (2021) Bosman, S. E. I., Ďurovčíková, D., Davies, F. B., & Eilers, A.-C. 2021, MNRAS, 503, 2077, doi: 10.1093/mnras/stab572
  * Cen & Safarzadeh (2015) Cen, R., & Safarzadeh, M. 2015, Astrophysical Journal Letters, 798, 2, doi: 10.1088/2041-8205/798/2/L38
  * Clough et al. (2005) Clough, S. A., Shephard, M. W., Mlawer, E. J., et al. 2005, Journal of Quantitative Spectroscopy and Radiative Transfer, 91, 233, doi: https://doi.org/10.1016/j.jqsrt.2004.05.058
  * Conroy & White (2013) Conroy, C., & White, M. 2013, Astrophysical Journal, 762, doi: 10.1088/0004-637X/762/2/70
  * Davies et al. (2016) Davies, F. B., Furlanetto, S. R., & McQuinn, M. 2016, Monthly Notices of the Royal Astronomical Society, 457, 3006, doi: 10.1093/mnras/stw055
  * Davies et al. (2019a) Davies, F. B., Hennawi, J. F., & Eilers, A.-C. 2019a, Monthly Notices of the Royal Astronomical Society, 14, 1, doi: 10.1093/mnras/stz3303
  * Davies et al. (2019b) —. 2019b, The Astrophysical Journal, 884, L19, doi: 10.3847/2041-8213/ab42e3
  * Davies et al. (2018a) Davies, F. B., Hennawi, J. F., Bañados, E., et al. 2018a, The Astrophysical Journal, 864, 143, doi: 10.3847/1538-4357/aad7f8
  * Davies et al. (2018b) —. 2018b, The Astrophysical Journal, 864, 142, doi: 10.3847/1538-4357/aad6dc
  * Decarli et al. (2018) Decarli, R., Walter, F., Venemans, B. P., et al. 2018, The Astrophysical Journal, 854, 97, doi: 10.3847/1538-4357/aaa5aa
  * Di Matteo et al. (2005) Di Matteo, T., Springel, V., & Ilernquist, L. 2005, Nature, 433, 604, doi: 10.1038/nature03335
  * Efstathiou & Rees (1988) Efstathiou, G., & Rees, M. J. 1988, Monthly Notices of the Royal Astronomical Society, 230, 5P, doi: 10.1093/mnras/230.1.5P
  * Eilers et al. (2017) Eilers, A.-C., Davies, F. B., Hennawi, J. F., et al. 2017, The Astrophysical Journ, 840, 24, doi: 10.3847/1538-4357/aa6c60
  * Eilers et al. (2018) Eilers, A.-C., Hennawi, J. F., & Davies, F. B. 2018, The Astrophysical Journal, 867, 30, doi: 10.3847/1538-4357/aae081
  * Eilers et al. (2021) Eilers, A.-C., Hennawi, J. F., Davies, F. B., & Simcoe, R. A. 2021, arXiv e-prints, arXiv:2106.04586.  https://arxiv.org/abs/2106.04586
  * Eilers et al. (2020) Eilers, A.-C., Hennawi, J. F., Decarli, R., et al. 2020.  http://arxiv.org/abs/2002.01811
  * Faber et al. (2003) Faber, S. M., Phillips, A. C., Kibrick, R. I., et al. 2003, in Instrument Design and Performance for Optical/Infrared Ground-based Telescopes, ed. M. Iye & A. F. M. Moorwood, Vol. 4841, International Society for Optics and Photonics (SPIE), 1657–1669, doi: 10.1117/12.460346
  * Fan et al. (2006) Fan, X., Strauss, M. A., Becker, R. H., et al. 2006, The Astronomical Journal, 132, 117, doi: 10.1086/504836
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, Publications of the Astronomical Society of the Pacific, 125, 306, doi: 10.1086/670067
  * Gullikson et al. (2014) Gullikson, K., Dodson-Robinson, S., & Kraus, A. 2014, The Astronomical Journal, 148, 53, doi: 10.1088/0004-6256/148/3/53
  * Haardt & Madau (2012) Haardt, F., & Madau, P. 2012, The Astrophysical Journal, 746, 125, doi: 10.1088/0004-637X/746/2/125
  * Habouzit et al. (2016) Habouzit, M., Volonteri, M., Latif, M., Dubois, Y., & Peirani, S. 2016, Monthly Notices of the Royal Astronomical Society, 463, 529, doi: 10.1093/mnras/stw1924
  * Haiman & Cen (2001) Haiman, Z., & Cen, R. 2001, in Astronomical Society of the Pacific Conference Series, ed. M. Umemura & H. Susa, 101
  * Haiman & Hui (2001) Haiman, Z., & Hui, L. 2001, The Astrophysical Journal, 547, 27, doi: 10.1086/318330
  * Hennawi & Prochaska (2007) Hennawi, J. F., & Prochaska, J. X. 2007, The Astrophysical Journal, 655, 735, doi: 10.1086/509770
  * Hopkins et al. (2005) Hopkins, P. F., Hernquist, L., Cox, T. J., et al. 2005, The Astrophysical Journal, 630, 716, doi: 10.1086/432463
  * Horne (1986) Horne, K. 1986, Publications of the Astronomical Society of the Pacific, 98, 609, doi: 10.1086/131801
  * Inayoshi et al. (2020) Inayoshi, K., Visbal, E., & Haiman, Z. 2020, Annual Review of Astronomy and Astrophysics, 58, 27, doi: 10.1146/annurev-astro-120419-014455
  * Keating et al. (2015) Keating, L. C., Haehnelt, M. G., Cantalupo, S., & Puchwein, E. 2015, Monthly Notices of the Royal Astronomical Society, 454, 681, doi: 10.1093/mnras/stv2020
  * Khrykin et al. (2016) Khrykin, I. S., Hennawi, J. F., McQuinn, M., & Worseck, G. 2016, The Astrophysical Journal, 824, 133, doi: 10.3847/0004-637x/824/2/133
  * Khrykin et al. (2019) Khrykin, I. S., Hennawi, J. F., & Worseck, G. 2019, Monthly Notices of the Royal Astronomical Society, 484, 3897, doi: 10.1093/mnras/stz135
  * Khrykin et al. (2021) Khrykin, I. S., Hennawi, J. F., Worseck, G., & Davies, F. B. 2021, Monthly Notices of the Royal Astronomical Society, 505, 649
  * Lidz et al. (2007) Lidz, A., McQuinn, M., Zaldarriaga, M., Hernquist, L., & Dutta, S. 2007, The Astrophysical Journal, 670, 39, doi: 10.1086/521974
  * Lodato & Natarajan (2006) Lodato, G., & Natarajan, P. 2006, Monthly Notices of the Royal Astronomical Society, 371, 1813, doi: 10.1111/j.1365-2966.2006.10801.x
  * Lukić et al. (2014) Lukić, Z., Stark, C. W., Nugent, P., et al. 2014, Monthly Notices of the Royal Astronomical Society, 446, 3697, doi: 10.1093/mnras/stu2377
  * Lusso et al. (2015) Lusso, E., Worseck, G., Hennawi, J. F., et al. 2015, Monthly Notices of the Royal Astronomical Society, 449, 4204, doi: 10.1093/mnras/stv516
  * Martini (2004) Martini, P. 2004, Coevolution of Black Holes and Galaxies (Cambridge University Press), 169
  * Martini & Weinberg (2001) Martini, P., & Weinberg, D. H. 2001, The Astrophysical Journal, 547, 12, doi: 10.1086/318331
  * Mazzucchelli et al. (2017) Mazzucchelli, C., Bañados, E., Venemans, B. P., et al. 2017, The Astrophysical Journal, 849, 91, doi: 10.3847/1538-4357/aa9185
  * Mortlock et al. (2011) Mortlock, D. J., Warren, S. J., Venemans, B. P., et al. 2011, Nature, 474, 616, doi: 10.1038/nature10159
  * Novak et al. (2011) Novak, G. S., Ostriker, J. P., & Ciotti, L. 2011, Astrophysical Journal, 737, doi: 10.1088/0004-637X/737/1/26
  * Onoue et al. (2019) Onoue, M., Kashikawa, N., Matsuoka, Y., et al. 2019, The Astrophysical Journal, 880, 77, doi: 10.3847/1538-4357/ab29e9
  * Paris et al. (2011) Paris, I., Petitjean, P., Rollinde, E., et al. 2011, Astronomy and Astrophysics, 530, 1, doi: 10.1051/0004-6361/201016233
  * Planck Collaboration (2020) Planck Collaboration. 2020, Astronomy & Astrophysics, doi: 10.1051/0004-6361/201833910
  * Prochaska et al. (2020) Prochaska, J. X., Hennawi, J. F., Westfall, K. B., et al. 2020, Journal of Open Source Software, 1, 1, doi: 10.21105/joss.02308
  * Schauer et al. (2017) Schauer, A. T., Regan, J., Glover, S. C., & Klessen, R. S. 2017, Monthly Notices of the Royal Astronomical Society, 471, 4878, doi: 10.1093/MNRAS/STX1915
  * Schmidt et al. (2018) Schmidt, T. M., Hennawi, J. F., Worseck, G., et al. 2018, The Astrophysical Journal, 861, 122, doi: 10.3847/1538-4357/aac8e4
  * Schmidt et al. (2017) Schmidt, T. M., Worseck, G., Hennawi, J. F., Prochaska, J. X., & Crighton, N. H. M. 2017, The Astrophysical Journal, 847, 81, doi: 10.3847/1538-4357/aa83ac
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, Astronomy and Astrophysics, 500, 33, doi: 10.16309/j.cnki.issn.1007-1776.2003.03.004
  * Sheinis et al. (2002) Sheinis, A. I., Bolte, M., Epps, H. W., et al. 2002, Publications of the Astronomical Society of the Pacific, 114, 851, doi: 10.1086/341706
  * Shen et al. (2009) Shen, Y., Strauss, M. A., Ross, N. P., et al. 2009, Astrophysical Journal, 697, 1656, doi: 10.1088/0004-637X/697/2/1656
  * Simcoe et al. (2010) Simcoe, R. A., Burgasser, A. J., Bochanski, J. J., et al. 2010, in Ground-based and Airborne Instrumentation for Astronomy III, ed. I. S. McLean, S. K. Ramsay, & H. Takami, Vol. 7735, International Society for Optics and Photonics (SPIE), 469–480, doi: 10.1117/12.856088
  * Simcoe et al. (2013) Simcoe, R. A., Burgasser, A. J., Schechter, P. L., et al. 2013, Publications of the Astronomical Society of the Pacific, 125, 270, doi: 10.1086/670241
  * Soltan (1982) Soltan, A. 1982, Monthly Notices of the Royal Astronomical Society, 200, 115
  * Springel et al. (2005) Springel, V., White, S. D., Jenkins, A., et al. 2005, Nature, 435, 629, doi: 10.1038/nature03597
  * Suzuki et al. (2005) Suzuki, N., Tytler, D., Kirkman, D., O’Meara, J. M., & Lubin, D. 2005, The Astrophysical Journal, 618, 592, doi: 10.1086/426062
  * Trakhtenbrot et al. (2017) Trakhtenbrot, B., Volonteri, M., & Natarajan, P. 2017, The Astrophysical Journal, 836, L1, doi: 10.3847/2041-8213/836/1/l1
  * Venemans et al. (2015) Venemans, B. P., Walter, F., Zschaechner, L., et al. 2015, The Astrophysical Journal, 816, 37, doi: 10.3847/0004-637x/816/1/37
  * Venemans et al. (2013) Venemans, B. P., Findlay, J. R., Sutherland, W. J., et al. 2013, Astrophysical Journal, 779, doi: 10.1088/0004-637X/779/1/24
  * Vernet et al. (2011) Vernet, J., Dekker, H., D’Odorico, S., et al. 2011, Astronomy & Astrophysics, 536, A105, doi: 10.1051/0004-6361/201117752
  * Visbal et al. (2014) Visbal, E., Haiman, Z., & Bryan, G. L. 2014, Monthly Notices of the Royal Astronomical Society, 445, 1056, doi: 10.1093/mnras/stu1794
  * Volonteri (2010) Volonteri, M. 2010, Astronomy and Astrophysics Review, 18, 279, doi: 10.1007/s00159-010-0029-x
  * Volonteri (2012) —. 2012, Science, 337, 544
  * Volonteri et al. (2015) Volonteri, M., Silk, J., & Dubus, G. 2015, Astrophysical Journal, 804, 1, doi: 10.1088/0004-637X/804/2/148
  * Wang et al. (2020) Wang, F., Davies, F. B., Yang, J., et al. 2020, The Astrophysical Journal, 896, 23, doi: 10.3847/1538-4357/ab8c45
  * Wang et al. (2013) Wang, R., Wagg, J., Carilli, C. L., et al. 2013, Astrophysical Journal, 773, doi: 10.1088/0004-637X/773/1/44
  * White et al. (2008) White, M., Martini, P., & Cohn, J. D. 2008, Monthly Notices of the Royal Astronomical Society, 390, 1179, doi: 10.1111/j.1365-2966.2008.13817.x
  * White et al. (2012) White, M., Myers, A. D., Ross, N. P., et al. 2012, Monthly Notices of the Royal Astronomical Society, 424, 933, doi: 10.1111/j.1365-2966.2012.21251.x
  * Worseck et al. (2021) Worseck, G., Khrykin, I. S., Hennawi, J. F., Prochaska, J. X., & Farina, E. P. 2021, 19, 1.  http://arxiv.org/abs/2101.01196
  * Wu et al. (2015) Wu, X. B., Wang, F., Fan, X., et al. 2015, Nature, 518, 512, doi: 10.1038/nature14241
  * Wyithe et al. (2005) Wyithe, J. S. B., Loeb, A., & Carilli, C. 2005, The Astrophysical Journal, 628, 575, doi: 10.1086/430874
  * Yang et al. (2020) Yang, J., Wang, F., Fan, X., et al. 2020, The Astrophysical Journal, 897, L14, doi: 10.3847/2041-8213/ab9c26
  * Yu & Tremaine (2002) Yu, Q., & Tremaine, S. 2002, Monthly Notices of the Royal Astronomical Society, 335, 965, doi: 10.1046/j.1365-8711.2002.05532.x


# Micro-machined deep silicon atomic vapor cells

S. Dyer  Department of Physics, SUPA, University of Strathclyde, Glasgow G4
0NG, United Kingdom  P. F. Griffin  Department of Physics, SUPA, University of
Strathclyde, Glasgow G4 0NG, United Kingdom  A. S. Arnold  Department of
Physics, SUPA, University of Strathclyde, Glasgow G4 0NG, United Kingdom  F.
Mirando  Kelvin Nanotechnology, University of Glasgow, Glasgow G12 8LS, United
Kingdom  D. P. Burt  Kelvin Nanotechnology, University of Glasgow, Glasgow G12
8LS, United Kingdom  E. Riis  Department of Physics, SUPA, University of
Strathclyde, Glasgow G4 0NG, United Kingdom  J. P. McGilligan
<EMAIL_ADDRESS>Department of Physics, SUPA, University of
Strathclyde, Glasgow G4 0NG, United Kingdom

###### Abstract

Using a simple and cost-effective water jet process, silicon etch depth
limitations are overcome to realize a $6\,$mm deep atomic vapor cell. While
the minimum silicon feature size was limited to a $1.5\,$mm width in these
first generation vapor cells, we successfully demonstrate a two-chamber
geometry by including a $\sim$25 mm meandering channel between the alkali pill
chamber and main interrogation chamber. We evaluate the impact of the channel
conductance on the introduction of alkali vapor density during the pill
activation process, and mitigate glass damage and pill contamination near the
main chamber. Finally, we highlight the improved signal achievable in the
$6\,$mm silicon cell compared to standard $2\,$mm path length silicon vapor
cells.

††preprint: AIP/123-QED

##  I Introduction

The micro-fabrication of atomic sensors has enabled a revolution in the field
deployment and commercialization of chip-scale atomic sensors. [1, 2, 3, 4] A
critical component at the core of many of these chip-scale devices is the
atomic vapor cell, consisting typically of a micro-electro-mechanical systems
(MEMS) based glass-silicon-glass stack bonded under vacuum. [5] Such MEMS
cells have been widely explored for compact sensors, forming hermetic seals
from both anodic and thermo-compression bonds to provide a mass-producible
cell manufacturing solution. [5, 6, 7] While this now mature technology has
enabled the transition to miniature components for atomic sensors, the methods
of manufacturing silicon have typically limited the available optical path
length to around $2\,$mm. These limitations arise due to dry etch phenomena
[8] and challenges associated with long duration wet etch processes. [9]
However, a longer optical path length would be beneficial for atomic sensors
due to the increased atomic absorption and the associated enhancement in
signal-to-noise ratio.

To circumvent the short line-of-sight path length available in micro-
fabricated cells, alternative light routing methods have been explored. A
notable example is a wet-etched elongated silicon cell, with reflective walls
at 54.7∘ from the wafer surface that direct the incident light through the
long axis of the silicon cell, to achieve an optical path length of $7\,$mm.
[10] However, the additional complexity required to achieve optically smooth
walls, as well as the inclusion of glass-etched grating structures, critically
aligned to manipulate the angle of incidence for light routing parallel to the
silicon surface, make this process costly and complex. While novel approaches
to light routing in silicon cells may enable an increased optical path length,
[11, 12] the benefits of simple deep line-of-sight cells would also be
valuable for chip-scale laser cooling, where the 6 mm silicon thickness could
enable a larger optical overlap volume within the cell. [13]

Recent investigations of laser cooling in MEMS cells have reported a
degradation of the achievable trapped atom number due to a limited silicon
frame thickness, which restricts the optical overlap volume within the cell.
[13] A viable solution for both hot and cold atom MEMS cells would be an
increased optical path length in line-of-sight cells.

Here we introduce a micro-machined silicon based vapor cell, manufactured with
water jet cutting to achieve an optical path length of $6\,$mm in a simple
line-of-sight geometry. An alkali pill deposited in the silicon frame provides
a natural Rb density within the cell, with a meandering channel used to reduce
material spray and glass damage in the main chamber of the cell. The step-
change from conventional semiconductor fabrication processes to micro-
machining for silicon-based vapor cell production greatly reduces the
manufacturing costs by mitigating the need for photolithography and etching
processes, while providing an improved optical path length to greatly increase
the capabilities of chip-scale atomic sensors.

##  II Fabrication and set-up

The silicon frame of the cell is manufactured using a high pressure water jet
cutting process. While this method is not limited to the depth restrictions
placed upon micro-fabricated silicon, the minimum feature size is degraded. In
this first generation cell, a minimum feature width of 1.5 mm was realized for
the conductance channel between the pill source and main interrogation
chamber. Non-line-of-sight conductance channels have been shown to be critical
in reducing material spray and glass discoloration in the interrogation region
of the vapor cell. [7, 14, 15] Previously, non-line-of-sight channels have
been dry etched with feature sizes on the order of $200\,\upmu$m$\times
200\,\upmu$m to greatly reduce the conductance into the main chamber, while
permitting the presence of a suitable vapor density. [14] Since the micro-
machined channel’s cross-sectional area is $\sim\times$220 larger than
conventional dry etched channels [14], a long meandering path was machined
between the pill location and main chamber, as shown in Fig. 1 (a). The main
cubic chamber was fabricated with dimensions $(6\,$mm$)^{3}$ to enable wider
beams to be used for an improved signal-to-noise ratio in atomic spectroscopy.

Figure 1:  (a): Illustration of the micro-machined silicon vapor cell. The
cell is composed of a borosilicate glass (BSG) - Si - BSG stack, with total
outer dimensions of ($13\,$mm$)^{2}\times 7\,$mm. The main cubic chamber has a
$(6\,$mm$)^{3}$ volume, connected to an alkali micro-pill dispenser via a
meandering channel with a cross-sectional area of $1.5\,$mm$\times 6\,$mm.
(b)-(c): Scanning electron microscope (SEM) images of the micro-machined
silicon surface.

Following the water jet cutting, the silicon inner wall surface roughness was
evaluated with surface profilometry and scanning electron microscopy (SEM).
SEM images shown in Fig. 1 (b)-(c) were used to qualitatively evaluate the
potential impact of wall roughness on the atomic vapor cell; as an increased
surface area could consume more alkali content and hinder the homogeneity of
applied wall coatings that may be required in later applications. The surface
profilometry measurements determined an arithmetic average wall roughness of
$R_{a}\,=\,5.9\,\upmu$m. While this wall roughness is comparable to recent
work in additively manufactured vacuum chambers [16], it remains orders of
magnitude larger than the wall roughness achievable for DRIE, with short etch
depths ($\leq$100 $\upmu$m) achieving $R_{a}=5\,$nm, [17, 8] and deep etched
silicon ($1.4\,$mm) demonstrating $R_{a}=300\,$nm [18]. In future
applications, chemical smoothing of the inner walls can be applied if the wall
roughness is found to hinder the atomic vapor cell’s long-term performance.
[18]

Following manufacturing, the silicon wafer was anodically bonded to $0.5\,$mm
thick borosilicate glass wafers at both interfaces in a total background
pressure of $\approx 10^{-5}\,$mbar. Prior to bonding, a commercially
available alkali pill dispenser (SAES Rb/AMAX/Pill1-0.6) was deposited in the
dedicated pill chamber of the silicon frame. The alkali pill was selected as
the atomic source due to the relatively clean sourcing compared to azide and
chloride compounds, [19] enabling applications in wavelength referencing. [20,
21] Post bonding, the wafer stack was diced into cells with total dimensions
measuring $13\,$mm$\times 13\,$mm$\times 7\,$mm (L$\times$W$\times$H). To
control the alkali vapor pressure, surface mounted resistors adhered to the
front and back glass surfaces were used to heat the cell. The steady-state
cell temperature was measured with a thermal camera.

A Yb fibre laser at $1070\,$nm was used for pill activation, and focused to
fill the pill surface. A digital single-lens reflex camera was used to extract
images of the cell window damage over time during the pill activation. A
volume Bragg-reflector (VBR) laser at $780\,$nm was fibre-coupled and aligned
through a polarizing beam-splitter (PBS) and the main chamber of the MEMS
cell. The light was retro-reflected for in-situ ‘saturated absorption
spectroscopy’, [22] with the photodiode using a $780\pm 2\,$nm band-pass
filter to avoid saturation from the Yb fibre laser during activation.

##  III Results

Figure 2:  Activation process of the alkali pill with laser heating: prior to
activation (a); during the $10\,$s activation with $3\,$W of $1070\,$nm light
(b); immediately after activation, where glass damage can be seen, as well as
the presence of Rb on the glass surfaces (c). The migration of pill released
material after activation was then observed (d)-(f). Images (g)-(i) show a
cell with no serpentine channel before, during, and after activation,
respectively. (j); Simultaneously measured strongest D2 transition absorption
from the 85Rb $F=3$ ground state in the main chamber during the activation
process for a channel and no-channel cell.

The alkali pill source in each individual cell was laser heated [23] for
$10\,$s with $3\,$W of $1070\,$nm light, focused to fill the $1\,$mm diameter
of the pill surface. The activation process and subsequent material diffusion
is highlighted in Fig. 2 (a)-(f). In this demonstration, the cell was heated
to 60∘C to drive alkali diffusion through the channel and into the main
chamber. Images of a cell with no channel present during activation are shown
in Fig. 2 (g)-(i). The peak absorption evolution of the D2 85Rb $F=3$ ground
state was simultaneously measured (Fig. 2 (j)), from in-situ saturated
absorption spectroscopy during the activation process of both cells types. For
these measurements the laser intensity was 1.6 mW/cm2.

The cell prior to pill activation is seen at $t=0\,$s in (Fig. 2 (a)),
followed by a $10\,$s pulse of laser heating with $3\,$W of $1070\,$nm light,
(shown at $t=5\,$s in Fig. 2 (b)). While it is noted that the pill could be
activated with lower laser powers, the chosen parameters ensured a
reproducible activation between the cells under test. In the short time
following activation, darkening of the glass is observed around the activation
site, which progressively migrates through the conductance channel between
$t=(30-10320)\,$s. Simultaneously, a silver spray of Rb metal can be seen
migrating through the channel, with the growth of atomic absorption plateauing
in the main chamber after $t\sim 3000\,$s. Additionally, an identical cell was
activated without heating during pill activation. In this scenario, the Rb
accumulated near the pill activation site on the rough walls created by the
mechanical micro-machining of the silicon frame. This was overcome by
subsequent heating of the cell to again encourage migration into the main
chamber.

In contrast, a comparable cell without the meandering channel was fabricated,
activated and observed in a similar manner (Fig. 2 (g)-(i)). The absence of
the channel between the pill and interrogation cell results in the activation
process spraying material across the glass and reducing the cell optical
access. Simultaneously, the absorption within the cell saturates much earlier
at $t\sim 400\,$s, as there is no diffusion time constant associated with
reaching the cell interrogation region.

Following the activation process and measurement of alkali vapor density
within the main cell, the cell absorption spectra were analysed, as shown in
Fig. 3. The spectra provided in Fig. 3 (a), was obtained with using a 1/e2
beam waist of $\sim 2\,$mm and peak intensity of $5\,$mW/cm2. For this room
temperature measurement, $\approx$7.5$\%$ peak D2 absorption was measured from
the 85Rb $F=3$ ground states. It is noted that the linewidth of the sub-
Doppler features is $\sim\times 1.4$ broader than comparable measurements in a
glass reference cell over a range of incident intensities. While the cause of
the broadening remains under investigation, the observation of a similar
broadening in other pill based MEMS cells indicates that the cause could be
from non-Rb content released during pill activation.

Fig. 3 (b) illustrates the peak absorption of the 85Rb $F=3$ ground state as a
function of the measured cell temperature. The blue data points and curve show
the measured and theoretically expected absorption at various temperature for
a $6\,$mm cell length respectively. [24]. The $6\,$mm cell data are compared
to theoretical curves for a typical $7\,$cm long glass blown cell (black), and
a $2\,$mm path length micro-fabricated cell (red). For these measurements, a
beam waist of $\sim 2\,$mm and a peak intensity of 0.2 mW/cm2 were used. We
observed good correlation between experimental and theoretical absorption. For
an operation temperature of 47∘C, the $6\,$mm cell is capable of achieving the
same absorption level as a room temperature $7\,$cm cell, while the 2 mm cell
would require heating to 60∘C. Additionally, the optical path length of the
deep silicon machined frame provides a sufficient signal-to-noise ratio for
laser stabilization to the atomic reference cell at room temperature.

Figure 3:  (a): Typical D2 line saturated absorption spectrum (single-shot)
for the $6\,$mm meandering micro-machined cell at room temperature. (b): Peak
absorption from the 85Rb $F=3$ ground state Doppler-broadened profile as a
function of temperature for retro-reflected vapor cells of $7\,$cm (black),
$6\,$mm (blue) and $2\,$mm (red) path-lengths. Blue data points represent
measurements taken with the meandering $6\,$mm vapor cell.

In future work, we will evaluate the performance of the deep micro-machined
cell as a wavelength reference for laser locking in a chip-scale laser cooling
system. Additionally, the ability to fabricate 6 mm thick MEMS cells is
directly transferable to the development of ultra-high-vacuum MEMS cells,
where the increased vacuum volume would benefit the achievable optical overlap
volume and trapped atom number. [25, 13, 26]

##  IV Conclusion

We have demonstrated the suitability of deep micro-machined silicon frames for
the production of MEMS atomic vapor cells. While the first demonstration has a
silicon feature width limitation of 1.5 mm, the inclusion of a meandering
channel has been demonstrated to sufficiently restrict unwanted glass
discoloration from the pill activation process outside of the main chamber,
for unhindered optical interrogation. The increased optical path length was
demonstrated in saturated absorption spectroscopy, with an atomic absorption
significantly larger than is achievable in standard 2 mm MEMS vapor cells. The
manufacturing simplicity and cost reduction of this technique makes this a
highly attractive routine for macro-scale cell fabrication in atomic sensors
where increased optical path lengths are desired.

###### Acknowledgements.

The authors would like to thank R. Boudot for useful conversations and careful
reading of the manuscript. Additionally, we thank B. Lewis and A. Bregazzi for
their reading of the manuscript. The authors acknowledge funding from Innovate
UK project number 50133. J. P. M. gratefully acknowledges funding from a Royal
Academy of Engineering Research Fellowship.

###  IV.1 Conflict of interest

The authors declare they have no competing interest. Approved for Public
Release, Distribution Unlimited.

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Knappe  _et al._ [2004] S. Knappe, V. Shah, P. D. D. Schwindt, L. Hollberg, J. Kitching, L.-A. Liew, and J. Moreland, “A microfabricated atomic clock,” Applied Physics Letters 85, 1460–1462 (2004).
  * Vicarini  _et al._ [2019] R. Vicarini, M. Abdel Hafiz, V. Maurice, N. Passilly, E. Kroemer, L. Ribetto, V. Gaff, C. Gorecki, S. Galliou, and R. Boudot, “Mitigation of temperature-induced light-shift effects in miniaturized atomic clocks,” IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control 66, 1962 – 1967 (2019).
  * Hunter  _et al._ [2018] D. Hunter, R. Jiménez-Martínez, J. Herbsommer, S. Ramaswamy, W. Li, and E. Riis, “Waveform reconstruction with a Cs based free-induction-decay magnetometer,” Opt. Express 26, 30523–30531 (2018).
  * Shah  _et al._ [2007] V. Shah, S. Knappe, P. D. D. Schwindt, and J. Kitching, “Subpicotesla atomic magnetometry with a microfabricated vapour cell,” Nature Photonics 1, 649–652 (2007).
  * Liew  _et al._ [2004] L.-A. Liew, S. Knappe, J. Moreland, H. Robinson, L. Hollberg, and J. Kitching, “Microfabricated alkali atom vapor cells,” Applied Physics Letters 84, 2694–2696 (2004).
  * Kitching [2018] J. Kitching, “Chip-scale atomic devices,” Applied Physics Reviews 5, 031302 (2018).
  * Karlen  _et al._ [2020] S. Karlen, J. Haesler, T. Overstolz, G. Bergonzi, and S. Lecomte, “Sealing of mems atomic vapor cells using Cu-Cu thermocompression bonding,” Journal of Microelectromechanical Systems 29, 95–99 (2020).
  * Laermer  _et al._ [2020] F. Laermer, S. Franssila, L. Sainiemi, and K. Kolari, “Chapter 16 - deep reactive ion etching,” in  _Handbook of Silicon Based MEMS Materials and Technologies (Third Edition)_, Micro and Nano Technologies, edited by M. Tilli, M. Paulasto-Krockel, M. Petzold, H. Theuss, T. Motooka, and V. Lindroos (Elsevier, 2020) third edition ed., pp. 417–446.
  * McGilligan  _et al._ [2022] J. P. McGilligan, K. F. Gallacher, P. F. Griffin, D. J. Paul, A. S. Arnold, and E. Riis, “Invited review: Micro-fabricated components for cold atom sensors,” Submitted to Rev. Sci. Instr. (2022).
  * Chutani  _et al._ [2015] R. Chutani, V. Maurice, N. Passilly, C. Gorecki, R. Boudot, M. Abdel Hafiz, P. Abbé, S. Galliou, J.-Y. Rauch, and E. de Clercq, “Laser light routing in an elongated micromachined vapor cell with diffraction gratings for atomic clock applications,” Scientific Reports 5, 14001 (2015).
  * Nishino, Furuya, and Ono [2021] H. Nishino, Y. Furuya, and T. Ono, “Micro-fabricated vapor cells with sealed Rb atoms by distillation at wafer level and two-step bonding for miniature atomic clocks,” Opt. Express 29, 44316–44321 (2021).
  * Nishino  _et al._ [2019] H. Nishino, M. Hara, Y. Yano, M. Toda, Y. Kanamori, M. Kajita, T. Ido, and T. Ono, “A reflection-type vapor cell using anisotropic etching of silicon for micro atomic clocks,” Applied Physics Express 12, 072012 (2019).
  * Bregazzi  _et al._ [2021a] A. Bregazzi, P. F. Griffin, A. S. Arnold, D. P. Burt, G. Martinez, R. Boudot, J. Kitching, E. Riis, and J. P. McGilligan, “A simple imaging solution for chip-scale laser cooling,” Applied Physics Letters 119, 184002 (2021a).
  * Newman  _et al._ [2019] Z. L. Newman, V. Maurice, T. Drake, J. R. Stone, T. C. Briles, D. T. Spencer, C. Fredrick, Q. Li, D. Westly, B. R. Ilic, B. Shen, M.-G. Suh, K. Y. Yang, C. Johnson, D. M. S. Johnson, L. Hollberg, K. J. Vahala, K. Srinivasan, S. A. Diddams, J. Kitching, S. B. Papp, and M. T. Hummon, “Architecture for the photonic integration of an optical atomic clock,” Optica 6, 680–685 (2019).
  * Maurice  _et al._ [2017] V. Maurice, J. Rutkowski, E. Kroemer, S. Bargiel, N. Passilly, R. Boudot, R. Chutani, S. Galliou, C. Gorecki, L. Mauri, and M. Moraja, “Microfabricated vapor cells for miniature atomic clocks based on post-sealing activated cesium dispensers,” in  _2017 Joint Conference of the European Frequency and Time Forum and IEEE International Frequency Control Symposium (EFTF/IFCS)_ (2017) pp. 636–637.
  * Cooper  _et al._ [2021] N. Cooper, L. Coles, S. Everton, I. Maskery, R. Campion, S. Madkhaly, C. Morley, J. O’Shea, W. Evans, R. Saint, P. Krüger, F. Oručević, C. Tuck, R. Wildman, T. Fromhold, and L. Hackermüller, “Additively manufactured ultra-high vacuum chamber for portable quantum technologies,” Additive Manufacturing 40, 101898 (2021).
  * Chen  _et al._ [2002] K.-S. Chen, A. Ayon, X. Zhang, and S. Spearing, “Effect of process parameters on the surface morphology and mechanical performance of silicon structures after deep reactive ion etching (drie),” Journal of Microelectromechanical Systems 11, 264–275 (2002).
  * Chutani  _et al._ [2014] R. Chutani, M. Hasegawa, V. Maurice, N. Passilly, and C. Gorecki, “Single-step deep reactive ion etching of ultra-deep silicon cavities with smooth sidewalls,” Sensors and Actuators A: Physical 208, 66–72 (2014).
  * Bopp, Maurice, and Kitching [2020] D. G. Bopp, V. M. Maurice, and J. E. Kitching, “Wafer-level fabrication of alkali vapor cells using in-situ atomic deposition,” Journal of Physics: Photonics 3, 015002 (2020).
  * Douahi  _et al._ [2007] A. Douahi, L. Nieradko, J. Beugnot, J. Dziuban, H. Maillote, R. Boudot, S. Guerandel, M. Moraja, C. Gorecki, and V. Giordano, “New vapor cell technology for chip scale atomic clock,” in  _2007 IEEE International Frequency Control Symposium Joint with the 21st European Frequency and Time Forum_ (2007) pp. 58–61.
  * Newman  _et al._ [2021] Z. L. Newman, V. Maurice, C. Fredrick, T. Fortier, H. Leopardi, L. Hollberg, S. A. Diddams, J. Kitching, and M. T. Hummon, “High-performance, compact optical standard,” Opt. Lett. 46, 4702–4705 (2021).
  * Smith and Hughes [2004] D. A. Smith and I. G. Hughes, “The role of hyperfine pumping in multilevel systems exhibiting saturated absorption,” American Journal of Physics 72, 631–637 (2004).
  * Griffin, Weatherill, and Adams [2005] P. F. Griffin, K. J. Weatherill, and C. S. Adams, “Fast switching of alkali atom dispensers using laser-induced heating,” Review of Scientific Instruments 76, 093102 (2005).
  * Siddons  _et al._ [2008] P. Siddons, C. S. Adams, C. Ge, and I. G. Hughes, “Absolute absorption on rubidium D lines: Comparison between theory and experiment,” J. Phys. B 41, 155004 (2008).
  * McGilligan  _et al._ [2020] J. P. McGilligan, K. R. Moore, A. Dellis, G. D. Martinez, E. de Clercq, P. F. Griffin, A. S. Arnold, E. Riis, R. Boudot, and J. Kitching, “Laser cooling in a chip-scale platform,” Applied Physics Letters 117, 054001 (2020).
  * Bregazzi  _et al._ [2021b] A. Bregazzi, S. Dyer, P. F. Griffin, D. P. Burt, A. S. Arnold, E. Riis, and J. P. McGilligan, “Enabling the mass production of a chip-scale laser cooling platform,” in  _Quantum Technology: Driving Commercialisation of an Enabling Science II_, Vol. 11881, edited by M. J. Padgett, K. Bongs, A. Fedrizzi, and A. Politi, International Society for Optics and Photonics (SPIE, 2021) pp. 100 – 110\.