# A Full Resolution of the 450 $\micron$ Extragalactic Background Light

Qi-Ning Hsu Graduate Institute of Astronomy, National Tsing Hua University,
No. 101, Section 2, Kuang-Fu Road, Hsinchu 30013, Taiwan  Academia Sinica
Institute of Astronomy and Astrophysics (ASIAA), No. 1, Section 4, Roosevelt
Rd., Taipei 10617, Taiwan  L. L. Cowie Institute for Astronomy, University of
Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA  Chian-Chou Chen Academia
Sinica Institute of Astronomy and Astrophysics (ASIAA), No. 1, Section 4,
Roosevelt Rd., Taipei 10617, Taiwan  A. J. Barger Department of Astronomy,
University of Wisconsin-Madison, 475 N. Charter Street, Madison, WI 53706, USA
Department of Physics and Astronomy, University of Hawaii, 2505 Correa Road,
Honolulu, HI 96822, USA  Institute for Astronomy, University of Hawaii, 2680
Woodlawn Drive, Honolulu, HI 96822, USA

###### Abstract

The extragalactic background light (EBL) is the cumulative radiation outside
the Milky Way. The determination of its corresponding primary emitting sources
as well as its total energy level across the entire electromagnetic spectrum
has profound implications for both cosmology and galaxy formation. However,
the detailed origin of the EBL at far-infrared wavelengths, particularly those
close to the peak of the cosmic infrared background, remains unclear. Here we
report the results of our ongoing SCUBA-2 450 $\micron$ survey of 10 massive
galaxy cluster fields. By exploiting the strong gravitational lensing offered
by these clusters, we obtain significant counts down to an unprecedented depth
of $\sim$0.1 mJy at this wavelength, about ten times deeper than that reached
by any other previous survey. The cumulative energy density based on the
counts is 138.1${}^{+23.9}_{-19.3}$ Jy degree-2, or 0.45${}^{+0.08}_{-0.06}$
MJy sr-1. Comparing our measurements to those made by the COBE and Planck
satellites, we find that at this flux density level, the 450 $\micron$ EBL is
entirely resolved by our SCUBA-2 observations. Thus, we find for the first
time that discrete sources produce fully to the 450 $\micron$ EBL, and that
about half of it comes from sources with sub-mJy flux densities. Our deep
number counts provide strong constraints on galaxy formation models.

galaxies: clusters — galaxies: number counts — EBL — submillimeter galaxies

††facilities: JCMT(SCUBA-2)††software: astropy (Astropy Collaboration et al.
2013, 2018), LENSTOOL (Kneib et al., 2011)

##  1 Introduction

The cosmic energy budget, encompassing the total amount of energy radiated
throughout the Universe, is a fundamental aspect of modern astrophysics.
Central to this budget is the extragalactic background light (EBL), which
represents the integrated emission across the electromagnetic spectrum from
astrophysical sources outside the Milky Way (see, e.g., the review by Cooray
2016). Observational measurements of the EBL at different wavelengths allow
one to understand its energy distribution, and, in principle, can provide
insights on the dominant contributors to the EBL (see, e.g., the discussion in
Hill et al. 2018).

One common way to put constraints on the EBL is to perform imaging surveys
using ground-based facilities, which allow one to construct source number
counts and calculate the integrated energy densities (e.g., the early results
in the far-infrared (FIR)/submillimeter by Smail et al. 1997; Barger et al.
1999; Cowie et al. 2002). By comparing with satellite measurements having
lower spatial resolution, which, in principle, take into account any diffuse
emission, one can estimate the energy contributions to the EBL from galaxies
and understand whether galaxies are the dominant contributors. Surveys have
found that the optical and near-infrared EBL primarily originate from the
directly observed star formation in galaxies, while the FIR/submillimeter EBL
mainly comes from the thermal emission from interstellar dust reradiated
starlight (Puget et al., 1996; Fixsen et al., 1998; Dole et al., 2006). These
dusty galaxies, sometimes called submillimeter galaxies (SMGs), are typically
characterized by intense dust emission indicating high rates of star formation
(see, e.g., review articles by Blain et al. 2002; Casey et al. 2014).

The advent of space-based telescopes like the Herschel Space Observatory and
ground-based facilities like the Atacama Large Millimeter/submillimeter Array
(ALMA) has revolutionized our ability to perform either wide-field or deep FIR
surveys, enabling unprecedented studies of galaxy number counts (Oliver et
al., 2010; Berta et al., 2011; Gómez-Guijarro et al., 2022; Cowie et al.,
2023; Chen et al., 2023). These observations, coupled with theoretical models,
provide a comprehensive view of the FIR Universe and its connection to other
astrophysical phenomena (e.g., Lacey et al. 2016; Lagos et al. 2019; Hayward
et al. 2021). However, due to the confusion limits of Herschel and the
inefficient survey capability of ALMA with its small field-of-view, the number
counts at wavelengths close to the peak of the FIR EBL ($\sim 200-300$ µm;
e.g., Odegard et al. 2019) remain limited to the brightest end (e.g., Oliver
et al. 2010).

SCUBA-2, a state-of-the-art submillimeter camera, offers unmatched sensitivity
and high angular resolution at 450 $\micron$. Its optimized design for deep
submillimeter observations makes it an excellent tool for studying dust-rich
galaxies (Holland et al., 2013). By harnessing the sensitivity of SCUBA-2 at
submillimeter wavelengths, we can detect and characterize faint sources at 450
$\micron$, enabling comprehensive studies of galaxy number counts close to the
peak of the FIR EBL and a deeper exploration of the submillimeter Universe
(Geach et al. 2013; Casey et al. 2013; Zavala et al. 2017; Wang et al. 2017;
Lim et al. 2020; Barger et al. 2022; Gao et al. 2024).

Additionally, the technique of gravitational lensing provides a unique means
for studying galaxy number counts at the fainter end (Smail et al., 1997;
Cowie et al., 2002; Knudsen et al., 2008; Johansson et al., 2011; Chen et al.,
2013; Hsu et al., 2016; Cowie et al., 2022). Gravitational lensing occurs when
the gravitational field of a massive object, such as a rich cluster of
galaxies, bends and magnifies the light from distant galaxies. By exploiting
the lensing effect, we can effectively boost the observed flux of background
galaxies, allowing us to probe deeper into the Universe and uncover fainter
sources that would otherwise remain undetected in blank-field surveys.

In this study, we utilize the deep data obtained with SCUBA-2 on 10 massive
galaxy cluster fields to derive robust number counts at 450 $\micron$. Our
analysis properly accounts for selection biases, completeness, and
uncertainties. By constructing the deepest 450 $\mu$m number counts ever, we
aim to unravel the relative contributions of SMGs to the 450 $\mu$m EBL. In
Section 2, we describe our SCUBA-2 data and data reduction. In Section 3, we
describe our methodology, including source extraction, Monte Carlo
simulations, and number counts calculations. In Section 4, we present our
number counts and the integrated energy density. We summarize our results in
Section 5.

##  2 Data and Data Reduction

We retrieved the SCUBA-2 450 $\mu$m data from the CADC archive. We used the
data taken under weather band 1 and band 2 conditions ($\tau_{225~{}\rm
GHz}~{}<$ 0.08) between October 2011 and July 2022 (Cowie et al., 2022). The
scan pattern used for these 10 cluster fields was CV DAISY, which has a
roughly circular field size of $\simeq 6^{\prime}$ in radius. We summarize the
data in Table 1.

We reduced the data following Chen et al. (2013). We used the Dynamic
Iterative Map Maker (DIMM) method in the SMURF package contained in the
STARLINK software (Chapin et al., 2013). This method models individual
components that make up the time-series recorded by the bolometer to produce
science maps. We adopted the “blank field” configuration file, which is
suitable for detecting faint point sources in extragalactic surveys. For each
cluster field, we produced scan maps with a pixel scale of $1^{\prime\prime}$
and then applied the recommended Flux Conversion Factors (FCFs) from Mairs et
al. (2021) to convert the pixel unit from picowatts to Jansky per beam. After
calibration, we used the MOSAIC_JCMT_IMAGES recipe from the Pipeline for
Combing and Analyzing Reduced Data (PICARD) to co-add and mosaic calibrated
scans for each field.

To improve source detection, we applied a matched filter using the
SCUBA2_MATCHED_FILTER recipe in PICARD. This recipe first convolves the map
with a Gaussian to estimate the low spatial frequency noise and then subtracts
it from the original map. We adopted the default $20^{\prime\prime}$ FWHM
value for the Gaussian profile. We verified the flux recovery capability of
SCUBA2_MATCHED_FILTER following Lim et al. (2020), and we adopted a mean
upward correction of $5.3\%$ for the flux loss to our 450 $\mu$m data.

Table 1: SUMMARY of SCUBA-2 450 $\micron$ DATA  Field | RA | Dec | Exposure | Central RMS
---|---|---|---|---
|  |  | (hours) | (mJy)
A370 | 39.9604 | -1.5856 | 29.1 | 2.19
A1689 | 197.8729 | -1.3411 | 27.7 | 2.16
A2390 | 328.3979 | 17.6867 | 39.6 | 2.28
A2744 | 3.5788 | -30.3894 | 24.5 | 3.08
MACS J0416.1-2403 | 64.0349 | -24.0724 | 22.1 | 2.42
MACS J0717.5+3745 | 109.4020 | 37.7564 | 51.5 | 1.51
MACS J1149.5+2223 | 117.3962 | 22.4030 | 35.0 | 1.46
MACS J1423.8+2404 | 215.9486 | 24.0778 | 43.8 | 1.76
MACS J2129.40741 | 322.3592 | -7.6906 | 12.6 | 4.07
RX J1347.5-1145 | 206.8775 | -11.7528 | 20.0 | 2.27

Figure 1: Left: An example SCUBA-2 450 $\micron$ flux density map of the
cluster MACS J1149.5+2223 with a 6′ radius circular footprint. The red circles
on the 450 $\mu$m map show the extracted $>3\sigma$ sources. The orange boxes
mark the sources below 1 mJy after delensing. The white contours represent the
CATS lens model for $z=1.5$ at magnification values of 1.5, 2.0, and 2.5
(moving inwards). Right: Histogram of the SNR values based on the map of MACS
J1149.5+2223. The orange region shows the SNR distribution in the jackknife
map. The blue represents detections in the data map.

##  3 Analysis

###  3.1 Source Extraction

Before source extraction, for each field we generated SCUBA-2 PSF models by
stacking the 10-20 highest signal-to-noise ratio (SNR) source images without
neighboring sources. This method inherently assumes that the typical source
size is much smaller than the beam size of SCUBA-2 in the submillimeter ($\sim
7\farcs 5$ at 450 $\micron$), which is supported by recent ALMA observations
that found typical sizes of subarcseconds (e.g., Simpson et al. 2015; Hodge et
al. 2016; Fujimoto et al. 2017; Gullberg et al. 2018; Tadaki et al. 2020). We
then fitted a double Gaussian profile to model the PSFs, and we used the best-
fit model PSFs for source extraction.

We performed the source extractions as in Hsu et al. (2016). We searched for
the maximum SNR pixel in the central circular region ($6^{\prime}$ radius).
(The pointing centers of the maps are normally close to the cluster
centroids.) We recorded the location and flux density of the pixel, then
subtracted a rescaled PSF centered at this pixel and searched for the next
maximum SNR pixel. Following Hsu et al. (2016), we used a 3$\sigma$ threshold,
which allows us to obtain a better SNR in the number counts at the faint end.
We repeated the extraction process until we reached this threshold. In Figure
1, we show one of our cluster fields as an example for the source extraction.

###  3.2 Delensed Raw Number Counts

We calculated delensed differential number counts at 450 $\micron$ by using
sources with SNR $>3.0$. To compute the demagnified flux densities, we used
the public software LENSTOOL (Kneib et al., 2011) to generate magnification
maps for our cluster fields. Since we do not have redshifts for these sources,
we corrected the flux densities of the sources by adopting estimated median
redshifts of 1.5 based on previous 450 $\micron$ studies (e.g., Chen et al.
2013; Casey et al. 2013). We took our lens models from the LENSTOOL developers
and Hubble Frontier Fields Archive (Ammons et al. 2014; Bradač et al. 2009;
Caminha et al. 2017; Diego et al. 2015; Grillo et al. 2015; Hoag et al. 2016;
Ishigaki et al. 2015; Jullo et al. 2007; Jullo & Kneib 2009; Jauzac et al.
2012, 2014; Johnson et al. 2014; Keeton 2010; Kawamata et al. 2018;
Liesenborgs et al. 2006; Lotz et al. 2017; Merten et al. 2011; McCully et al.
2014; Mohammed et al. 2014; Oguri 2010; Richard et al. 2014; Zitrin et al.
2013).

We determined the magnification factors $\mu_{\rm i}$ of the point sources
from the magnification maps. The demagnified flux density of each source can
be obtained from

$S_{\rm demag,\emph{i}}=\frac{S_{\rm obs,\emph{i}}}{\mu_{i}}\,,$ (1)

where $S_{\rm demag,\emph{i}}$ and $S_{\rm obs,\emph{i}}$ are demagnified and
observed flux densities, respectively. We calculated the effective area
A${}_{\rm eff,\emph{i}}$ of each source on the source plane. We summed over
the pixels whose SNRs were greater than 3$\sigma$ and then converted pixels to
square degrees. We then calculated the delensed raw counts at the $j$–th flux
bin as

$\frac{dN_{\rm raw,j}}{dS_{j}}=\frac{1}{\Delta S_{j}}\sum_{i}^{n}X_{i}$ (2)

and

$X_{i}=\frac{1}{A_{\rm eff,\emph{i}}}\,.$ (3)

Here, Xi represents the number density contribution of each source within that
flux bin. We based the error calculation on Poisson statistics.

###  3.3 Simulations

Corrections for flux boosting, false detections, and incompleteness are needed
in order to obtain the intrinsic number counts. To do this, we ran Monte Carlo
simulations to find the underlying models for our fields. We used the
Schechter function form as our number counts model:

$\rm\frac{dN}{dS}=\Bigl{(}\frac{N_{0}}{S_{0}}\Bigr{)}~{}\Bigl{(}\frac{S}{S_{0}}\Bigr{)}^{\alpha}~{}exp\Bigl{(}-\frac{S}{S_{0}}\Bigr{)}\,.$
(4)

We generated artificial sources whose flux densities were assigned according
to the underlying models for each cluster field. We then randomly distributed
these sources in the source plane. Next, we projected the simulated sources
onto the image plane using LENSTOOL. The outputs of LENSTOOL contain the new
flux densities and positions of the simulated sources in the image plane. We
convolved the simulated sources with the PSF and added them into the jackknife
maps to produce mock observation maps. We produced jackknife maps following
Hsu et al. (2016) by coadding even and odd scans separately. We then
subtracted these two coadded maps and rescaled the value of each pixel by a
factor of $\sqrt{t_{\rm even}}\times\sqrt{t_{\rm odd}}$/(teven \+ todd), where
teven and todd represent the integration times of each pixel from the even and
odd coadded maps, respectively. We also applied matched-filtering to the
jackknife maps, as we did for the real data images.

In order to estimate the true number counts, we adopted an iterative procedure
in our simulations. We generated 15–20 mock maps in each iteration step. We
then performed source extraction and computed the averaged recovered counts of
the mock maps. Next, we corrected the raw counts by using the ratio between
the averaged recovered counts and the raw counts. Finally, we did a $\chi^{2}$
fit to the corrected counts by using the Schechter function. This fit will be
the new counts model for the next iteration. We iterated until the recovered
counts converged with the raw counts to within the 1$\sigma$ uncertainties for
all the flux bins.

Once we obtained the intrinsic number counts models, we used these to produce
500 mock images and then performed the source extractions to create the source
catalogs. Following Gao et al. (2024), we cross-matched the input and output
catalogs within a 1/2 beam FWHM as our search area to find the brightest
counterparts. We considered Soutput/S${}_{\rm input}\leq 3$ a match in this
analysis, which is similar to what has been adopted previously (e.g., Geach et
al. 2017). Following Gao et al. (2024), we estimated the boosting factors,
false detection rates, and completeness of the point sources using a
two–dimensional binning method in our cross-matched catalogs.

###  3.4 Delensed Corrected Number Counts

To derive delensed corrected number counts, we first deboosted the fluxes of
the sources. We then delensed their fluxes and estimated their effective areas
in the source plane. Finally, we derived the delensed corrected number counts
using

$\frac{dN_{\rm corr,j}}{dS_{j}}=\frac{1}{\Delta S_{j}}\sum_{i}^{n}X_{\rm
corr,\emph{i}}$ (5)

and

$X_{\rm corr,\emph{i}}=\frac{1-p_{\rm false,\emph{i}}}{C_{i}A_{\rm
eff,\emph{i}}}\,,$ (6)

where p${}_{\rm false,\emph{i}}$ is the false detection rate, and Ci is the
completeness. We confirmed that the delensed corrected number counts agree
with the intrinsic counts models. We estimated the statistical uncertainties
from the source positions and took these uncertainties into account in the
analyses.

###  3.5 Systematic Uncertainties

Systematic uncertainties related to the modeling of gravitational lensing need
to be taken into account for a proper assessment of the number counts error
budget. In the following, we give our results with estimates of the systematic
uncertainties caused by the redshift distribution of the background lensed
SMGs, the lens models, and the clustering of the source distributions.

In our methodology, we calculated corrected number counts by assuming a median
redshift of 1.5. However, it is expected that the SMGs have a redshift
distribution, which could affect the magnification estimates. This leads to
additional uncertainties in the number counts. To address this, we randomly
assigned redshifts to our sources using the redshift distribution from the
STUDIES survey, which is the deepest 450 $\micron$ blank-field survey (Wang et
al. 2017; Dudzevičiūtė et al. 2021; Gao et al. 2024). We calculated the
corrected number counts again and compared the standard deviations of the
number counts in the randomized redshift sample with the Poisson noise. We
found that the uncertainties caused by the assumptions concerning the redshift
are subdominant, and they are, on average, about 25% of the Poisson
uncertainties. Nevertheless, we include this error budget in the total error
budget.

Table 2: Combined Differential Number Counts at 450 $\micron$

S${}_{450~{}\micron}$ | log10(dN/dS) | $<\mu>^{1}$ | N${}_{\rm total}^{2}$
---|---|---|---
(mJy) | (mJy-1 deg-2) |  |
0.13 | 6.2${}_{-1.4}^{+0.4}$ | 34.4 | 3
0.25 | 5.6${}_{-99.9}^{+0.6}$ | 22.4 | 5
0.45 | 5.1${}_{-0.7}^{+0.3}$ | 14.9 | 12
0.81 | 4.7${}_{-99.9}^{+0.5}$ | 7.4 | 3
1.48 | 4.2${}_{-0.5}^{+0.3}$ | 4.1 | 26
2.70 | 3.5${}_{-0.2}^{+0.2}$ | 2.1 | 108
4.91 | 3.0${}_{-0.1}^{+0.1}$ | 1.6 | 232
8.94 | 2.4${}_{-0.1}^{+0.1}$ | 1.2 | 358
16.28 | 1.5${}_{-0.2}^{+0.2}$ | 1.1 | 228
29.65 | 0.6${}_{-0.8}^{+0.4}$ | 1.2 | 23

11footnotetext: Averaged gravitational magnification

22footnotetext: Total number of sources in each flux bin

To estimate the systematic uncertainties caused by the lens modeling, we
utilized the various lens models provided for the five Hubble Frontier Fields
(A370, A2744, MACS J0416.1-2403, MACS J0717.5+3745, and MACS J1149.5+2223). We
ran the same Monte Carlo simulations as those done on the real data but using
different lens models. We then estimated the systematic uncertainties by
calculating the standard deviations of the corrected counts obtained using the
different lens models. We found the systematic uncertainties to be
subdominant, again about 25% of the Poisson uncertainties. We include this
error budget in the total error budget.

We note that in our Monte Carlo simulations, we randomly distributed the
positions of the artificial sources. This might bias our counts, since we do
not take clustering effects into account when we calculate the boosting
factors, false detection rates, and completeness from the mock maps. To test
whether neglecting clustering could significantly alter our counts results, we
used the empirical catalogs produced by the SIDES simulation (Béthermin et
al., 2017), which inherently include clustering effects, since the simulation
builds upon dark matter lightcones.

Before assessing the clustering effects, we first validated our methodology
for estimating intrinsic counts by performing the following test. We clipped
the original SIDES 2 deg2 simulated map into a set of 25 smaller cutout maps
with a size similar to our SCUBA-2 footprint, and we treated them as different
cluster fields in the source planes. We then lensed these simulated maps to
the image planes using LENSTOOL and used our methodology to find the corrected
counts. We compared the corrected counts to the true counts provided by SIDES.
In this test, we adopted the three-dimensional source positions from the SIDES
catalog instead of having the sources randomly distributed on the sky with an
assigned redshift. We show our results in Figure 2, where the averaged
corrected number counts over the 25 cutout maps are consistent with the true
number counts.

After validating our methodology, we moved on with a similar test. This time
we randomly distributed the source positions to calculate the corrected number
counts using the smaller SIDES cutout maps. We show our results in Figure 2.
While not taking clustering into account can lead to larger uncertainties,
evidently with a larger scatter in the averaged counts, on average, there is
no significant difference between the SIDES counts and the corrected counts.
Our results are therefore consistent with previous studies—either from other
SCUBA-2 450 $\micron$ surveys (Wang et al., 2017), or from Herschel studies
with similar beam sizes but at slightly shorter wavelengths (Béthermin et al.,
2015)—where no significant impact of clustering on the number counts was
found.

Figure 2: Ratios between the corrected number counts and the model counts,
where the model counts are based on the SIDES simulation source catalog, and
the corrected counts were obtained by applying our methodology to the mock
images that were made based on the SIDES source catalog. The top panel shows
the results for the case where the source positions are adopted from SIDES,
which includes clustering. The bottom panel shows the results for the case
where the sky positions of the sources are randomized.

##  4 Results and Discussion

We present our corrected differential number counts for the 10 lensing cluster
fields in the left panel of Figure 3. The solid black curve is the best-fit
Schechter function for our corrected counts, which can be parameterized by
Equation 4 with the following parameters: N0 = 4437.5$\pm$1399.5, S0 =
10.4$\pm$1.7, and $\alpha$=-1.9$\pm$0.1. The uncertainties on each data point
include Poisson noise and the uncertainties from the redshift distribution and
the lens model. Thanks to the powerful effects of strong gravitational
lensing, we have pushed the detection limit down to $\sim$0.1 mJy at 450
$\micron$, a factor of $>$10 improvement over the deepest blank-field counts
(Wang et al. 2017; Gao et al. 2024).

We also calculate the weighted average counts by combining all 10 fields. We
show the results in the right panel of Figure 3, and we provide the
corresponding values in Table 2. Our results are consistent with the previous
measurements shown in Figure 3. However, the various physical or empirical
models (dotted and dashed curves in the right panel) tend to overpredict
source densities at the bright end at 450 $\micron$ ($\gtrsim$1 mJy) by about
10-30%. On the other hand, at the faint end ($\lesssim$1 mJy), our counts are
slightly higher, although not significantly. The physical reasons for the
disagreement between the measurements and the models at the bright end need to
be investigated. Gao et al. (2024) proposed that it could be due to the
mismatch in halo masses between those inferred from clustering measurements
(Lim et al., 2020) and those in the models. Studies of the physical properties
of the 450 $\micron$ sources, such as stellar mass, could potentially shed
more light on this issue.

Figure 3:  Left: Corrected differential number counts for all 10 cluster
fields at 450 $\mu$m. The solid black curve shows the best-fit Schechter
function, while the shaded gray region is the 68% confidence interval. The
uncertainties include Poisson noise, as well as those caused by the underlying
redshift distributions and lens models. Right: Weighted average counts of all
10 cluster fields in filled black circles, along with previous 450 $\micron$
surveys in open symbols (Chen et al. 2013; Hsu et al. 2016; Zavala et al.
2017; Wang et al. 2017; Gao et al. 2024). Model predictions from SIDES
(Béthermin et al., 2017), GALFORM (Cowley et al., 2019), and SHARK (Lagos et
al., 2020) are shown as dotted and dashed curves.

To estimate the contributions of the detected sources to the 450 $\micron$
EBL, we integrated the best-fit Schechter model, which we show as a function
of flux in Figure 4. Above 0.1 mJy where significant counts are obtained, we
find the total energy density at 450 $\mu$m to be 138.1${}^{+23.9}_{-19.3}$ Jy
deg-2, or 0.45${}^{+0.08}_{-0.06}$ MJy sr-1. This corresponds to
103${}^{+18}_{-14}\%$ of the total 450 $\mu$m EBL reported by Odegard et al.
(2019).

Past works have suggested that a broken power law could be a viable
alternative to the Schechter function for the underlying counts models. It can
be described as

$\frac{dN}{dS}=\begin{cases}N_{0}(\frac{S}{S_{0}})^{-\alpha},&if\ S\leq
S_{0}\\\ N_{0}(\frac{S}{S_{0}})^{-\beta},&if\ S>S_{0}\end{cases}$ (7)

We ran further counts analyses based on this model form and found consistent
results compared to those obtained based on Schechter functions. The best-fit
parameters of the broken power law are N0=151.0$\pm$14.6, S0=9.6$\pm$0.4,
$\alpha$=2.0$\pm$0.1, $\beta$=5.6$\pm$1.1. We plot the cumulative energy
density of the best-fit broken power law in Figure 4. The cumulative energy
density based on the broken power law is 126.9${}^{+41.6}_{-41.6}$ Jy deg2.

Our work demonstrates for the first time that discrete sources, are the
dominant contributors to the 450 $\micron$ EBL. Interestingly, about half of
the contribution comes from sources that are fainter than $\sim$1 mJy, below
the typical confusion limit of SCUBA-2 450 $\micron$ images (Gao et al.,
2024).

Noticeably, as shown in Figure 4, the integrated energy density does not
converge when integrating down to 0.14 mJy, suggesting that deeper data are
needed in order to put a tighter constraint on the faint end counts and thus
obtain a converged constraint on the EBL contributions from discrete sources.

Figure 4: Cumulative EBL as a function of flux density at 450 $\mu$m. The
black solid curve was calculated based on our best-fit Schechter function with
the 1$\sigma$ uncertainties in gray shading. The blue solid curve was
calculated based on a broken power law function. The black dotted curve was
calculated by combining our results with the deepest blank-field counts from
Gao et al. (2024). The horizontal red dashed line is the EBL measured by the
COBE$+$Planck satellites (Odegard et al., 2019). The red shaded region shows
the range from COBE estimates (e.g., Puget et al. 1996; Fixsen et al. 1998;
Gispert et al. 2000).

##  5 Summary

In summary, our research represents a significant step forward in the study of
the 450 $\micron$ EBL. Through the innovative combination of gravitational
lensing and the unparalled capability of SCUBA-2, we have achieved a complete
resolution of the 450 $\micron$ EBL and established the dominance of discrete
SMGs as its main contributors. These findings provide a broader understanding
of the FIR/submillimeter regime, the cosmic energy distribution, and the
interplay between galaxies and the diffuse background radiation. Our
measurements could also be helpful for the design of the next generation
submillimeter facilities, such as LST (Kohno et al., 2020) and AtLAST
(Klaassen et al., 2020), which aim to obtain wide blank-field images to a
depth that is similar to what has been reached by this work.

We thank the reviewer for a useful report that has improved the manuscript.
Q.-N.H. and C.-C.C acknowledge support from the National Science and
Technology Council of Taiwan (NSTC 109-2112-M-001-016-MY3 and
111-2112M-001-045-MY3), as well as Academia Sinica through the Career
Development Award (AS-CDA-112-M02). L.L.C acknowledges support from NASA grant
80NSSC22K0483. A.J.B. acknowledges support from a Kellett Mid-Career Award and
a WARF Named Professorship from the University of Wisconsin-Madison Office of
the Vice Chancellor for Research and Graduate Education with funding from the
Wisconsin Alumni Research Foundation. The James Clerk Maxwell Telescope is
operated by the East Asian Observatory on behalf of The National Astronomical
Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics;
the Korea Astronomy and Space Science Institute; the National Astronomical
Research Institute of Thailand; Center for Astronomical Mega-Science (as well
as the National Key R&D Program of China with No. 2017YFA0402700). Additional
funding support is provided by the Science and Technology Facilities Council
of the United Kingdom and participating universities and organizations in the
United Kingdom and Canada. Additional funds for the construction of SCUBA-2
were provided by the Canada Foundation for Innovation. The authors wish to
recognize and acknowledge the very significant cultural role and reverence
that the summit of Maunakea has always had within the indigenous Hawaiian
community. We are most fortunate to have the opportunity to conduct
observations from this mountain. This work utilizes gravitational lensing
models produced by PIs Bradač, Natarajan & Kneib (CATS), Merten & Zitrin,
Sharon, Williams, Keeton, Bernstein and Diego, and the GLAFIC group. This lens
modeling was partially funded by the HST Frontier Fields program conducted by
STScI. STScI is operated by the Association of Universities for Research in
Astronomy, Inc. under NASA contract NAS 5-26555. The lens models were obtained
from the Mikulski Archive for Space Telescopes (MAST).

##  Appendix A Appendix A

In Figure 5 we show the imaging data for the remaining nine fields. The
symbols follow those used in Figure 1.

Figure 5: Flux density maps for the remaining nine fields. The format follows
that adopted in the left panel of Figure 1.

## References

  * Ammons et al. (2014) Ammons, S. M., Wong, K. C., Zabludoff, A. I., & Keeton, C. R. 2014, ApJ, 781, 2, doi: 10.1088/0004-637X/781/1/2
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Barger et al. (2022) Barger, A. J., Cowie, L. L., Blair, A. H., & Jones, L. H. 2022, ApJ, 934, 56, doi: 10.3847/1538-4357/ac67e7
  * Barger et al. (1999) Barger, A. J., Cowie, L. L., & Richards, E. A. 1999, in Astronomical Society of the Pacific Conference Series, Vol. 191, Photometric Redshifts and the Detection of High Redshift Galaxies, ed. R. Weymann, L. Storrie-Lombardi, M. Sawicki, & R. Brunner, 279, doi: 10.48550/arXiv.astro-ph/9907022
  * Berta et al. (2011) Berta, S., Magnelli, B., Nordon, R., et al. 2011, A&A, 532, A49, doi: 10.1051/0004-6361/201116844
  * Béthermin et al. (2015) Béthermin, M., Daddi, E., Magdis, G., et al. 2015, A&A, 573, A113, doi: 10.1051/0004-6361/201425031
  * Béthermin et al. (2017) Béthermin, M., Wu, H.-Y., Lagache, G., et al. 2017, A&A, 607, A89, doi: 10.1051/0004-6361/201730866
  * Blain et al. (2002) Blain, A. W., Smail, I., Ivison, R. J., Kneib, J. P., & Frayer, D. T. 2002, Phys. Rep., 369, 111, doi: 10.1016/S0370-1573(02)00134-5
  * Bradač et al. (2009) Bradač, M., Treu, T., Applegate, D., et al. 2009, ApJ, 706, 1201, doi: 10.1088/0004-637X/706/2/1201
  * Caminha et al. (2017) Caminha, G. B., Grillo, C., Rosati, P., et al. 2017, A&A, 600, A90, doi: 10.1051/0004-6361/201629297
  * Casey et al. (2014) Casey, C. M., Narayanan, D., & Cooray, A. 2014, Phys. Rep., 541, 45, doi: 10.1016/j.physrep.2014.02.009
  * Casey et al. (2013) Casey, C. M., Chen, C.-C., Cowie, L. L., et al. 2013, MNRAS, 436, 1919, doi: 10.1093/mnras/stt1673
  * Chapin et al. (2013) Chapin, E. L., Berry, D. S., Gibb, A. G., et al. 2013, MNRAS, 430, 2545, doi: 10.1093/mnras/stt052
  * Chen et al. (2013) Chen, C.-C., Cowie, L. L., Barger, A. J., et al. 2013, ApJ, 776, 131, doi: 10.1088/0004-637X/776/2/131
  * Chen et al. (2023) Chen, J., Ivison, R. J., Zwaan, M. A., et al. 2023, MNRAS, 518, 1378, doi: 10.1093/mnras/stac2989
  * Cooray (2016) Cooray, A. 2016, Royal Society Open Science, 3, 150555, doi: 10.1098/rsos.150555
  * Cowie et al. (2023) Cowie, L. L., Barger, A. J., & Bauer, F. E. 2023, ApJ, 952, 28, doi: 10.3847/1538-4357/acd763
  * Cowie et al. (2022) Cowie, L. L., Barger, A. J., Bauer, F. E., et al. 2022, ApJ, 939, 5, doi: 10.3847/1538-4357/ac91d2
  * Cowie et al. (2002) Cowie, L. L., Barger, A. J., & Kneib, J. P. 2002, AJ, 123, 2197, doi: 10.1086/339978
  * Cowley et al. (2019) Cowley, W. I., Lacey, C. G., Baugh, C. M., et al. 2019, MNRAS, 487, 3082, doi: 10.1093/mnras/stz1398
  * Diego et al. (2015) Diego, J. M., Broadhurst, T., Benitez, N., et al. 2015, MNRAS, 446, 683, doi: 10.1093/mnras/stu2064
  * Dole et al. (2006) Dole, H., Lagache, G., Puget, J. L., et al. 2006, A&A, 451, 417, doi: 10.1051/0004-6361:20054446
  * Dudzevičiūtė et al. (2021) Dudzevičiūtė, U., Smail, I., Swinbank, A. M., et al. 2021, MNRAS, 500, 942, doi: 10.1093/mnras/staa3285
  * Fixsen et al. (1998) Fixsen, D. J., Dwek, E., Mather, J. C., Bennett, C. L., & Shafer, R. A. 1998, ApJ, 508, 123, doi: 10.1086/306383
  * Fujimoto et al. (2017) Fujimoto, S., Ouchi, M., Shibuya, T., & Nagai, H. 2017, ApJ, 850, 83, doi: 10.3847/1538-4357/aa93e6
  * Gao et al. (2024) Gao, Z.-K., Lim, C.-F., Wang, W.-H., et al. 2024, ApJ, submitted
  * Geach et al. (2013) Geach, J. E., Chapin, E. L., Coppin, K. E. K., et al. 2013, MNRAS, 432, 53, doi: 10.1093/mnras/stt352
  * Geach et al. (2017) Geach, J. E., Dunlop, J. S., Halpern, M., et al. 2017, MNRAS, 465, 1789, doi: 10.1093/mnras/stw2721
  * Gispert et al. (2000) Gispert, R., Lagache, G., & Puget, J. L. 2000, A&A, 360, 1, doi: 10.48550/arXiv.astro-ph/0005554
  * Gómez-Guijarro et al. (2022) Gómez-Guijarro, C., Elbaz, D., Xiao, M., et al. 2022, A&A, 658, A43, doi: 10.1051/0004-6361/202141615
  * Grillo et al. (2015) Grillo, C., Suyu, S. H., Rosati, P., et al. 2015, ApJ, 800, 38, doi: 10.1088/0004-637X/800/1/38
  * Gullberg et al. (2018) Gullberg, B., Swinbank, A. M., Smail, I., et al. 2018, ApJ, 859, 12, doi: 10.3847/1538-4357/aabe8c
  * Hayward et al. (2021) Hayward, C. C., Sparre, M., Chapman, S. C., et al. 2021, MNRAS, 502, 2922, doi: 10.1093/mnras/stab246
  * Hill et al. (2018) Hill, R., Masui, K. W., & Scott, D. 2018, Applied Spectroscopy, 72, 663, doi: 10.1177/0003702818767133
  * Hoag et al. (2016) Hoag, A., Huang, K. H., Treu, T., et al. 2016, ApJ, 831, 182, doi: 10.3847/0004-637X/831/2/182
  * Hodge et al. (2016) Hodge, J. A., Swinbank, A. M., Simpson, J. M., et al. 2016, ApJ, 833, 103, doi: 10.3847/1538-4357/833/1/103
  * Holland et al. (2013) Holland, W. S., Bintley, D., Chapin, E. L., et al. 2013, MNRAS, 430, 2513, doi: 10.1093/mnras/sts612
  * Hsu et al. (2016) Hsu, L.-Y., Cowie, L. L., Chen, C.-C., Barger, A. J., & Wang, W.-H. 2016, ApJ, 829, 25, doi: 10.3847/0004-637X/829/1/25
  * Ishigaki et al. (2015) Ishigaki, M., Kawamata, R., Ouchi, M., et al. 2015, ApJ, 799, 12, doi: 10.1088/0004-637X/799/1/12
  * Jauzac et al. (2012) Jauzac, M., Jullo, E., Kneib, J.-P., et al. 2012, MNRAS, 426, 3369, doi: 10.1111/j.1365-2966.2012.21966.x
  * Jauzac et al. (2014) Jauzac, M., Clément, B., Limousin, M., et al. 2014, MNRAS, 443, 1549, doi: 10.1093/mnras/stu1355
  * Johansson et al. (2011) Johansson, D., Sigurdarson, H., & Horellou, C. 2011, A&A, 527, A117, doi: 10.1051/0004-6361/201016138
  * Johnson et al. (2014) Johnson, T. L., Sharon, K., Bayliss, M. B., et al. 2014, ApJ, 797, 48, doi: 10.1088/0004-637X/797/1/48
  * Jullo & Kneib (2009) Jullo, E., & Kneib, J. P. 2009, MNRAS, 395, 1319, doi: 10.1111/j.1365-2966.2009.14654.x
  * Jullo et al. (2007) Jullo, E., Kneib, J. P., Limousin, M., et al. 2007, New Journal of Physics, 9, 447, doi: 10.1088/1367-2630/9/12/447
  * Kawamata et al. (2018) Kawamata, R., Ishigaki, M., Shimasaku, K., et al. 2018, ApJ, 855, 4, doi: 10.3847/1538-4357/aaa6cf
  * Keeton (2010) Keeton, C. R. 2010, General Relativity and Gravitation, 42, 2151, doi: 10.1007/s10714-010-1041-1
  * Klaassen et al. (2020) Klaassen, P. D., Mroczkowski, T. K., Cicone, C., et al. 2020, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 11445, Ground-based and Airborne Telescopes VIII, ed. H. K. Marshall, J. Spyromilio, & T. Usuda, 114452F, doi: 10.1117/12.2561315
  * Kneib et al. (2011) Kneib, J.-P., Bonnet, H., Golse, G., et al. 2011, LENSTOOL: A Gravitational Lensing Software for Modeling Mass Distribution of Galaxies and Clusters (strong and weak regime), Astrophysics Source Code Library, record ascl:1102.004.  http://ascl.net/1102.004
  * Knudsen et al. (2008) Knudsen, K. K., van der Werf, P. P., & Kneib, J. P. 2008, MNRAS, 384, 1611, doi: 10.1111/j.1365-2966.2007.12820.x
  * Kohno et al. (2020) Kohno, K., Kawabe, R., Tamura, Y., et al. 2020, in Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy X, ed. J. Zmuidzinas & J.-R. Gao (SPIE), doi: 10.1117/12.2561238
  * Lacey et al. (2016) Lacey, C. G., Baugh, C. M., Frenk, C. S., et al. 2016, MNRAS, 462, 3854, doi: 10.1093/mnras/stw1888
  * Lagos et al. (2020) Lagos, C. d. P., da Cunha, E., Robotham, A. S. G., et al. 2020, MNRAS, 499, 1948, doi: 10.1093/mnras/staa2861
  * Lagos et al. (2019) Lagos, C. d. P., Robotham, A. S. G., Trayford, J. W., et al. 2019, MNRAS, 489, 4196, doi: 10.1093/mnras/stz2427
  * Liesenborgs et al. (2006) Liesenborgs, J., De Rijcke, S., & Dejonghe, H. 2006, MNRAS, 367, 1209, doi: 10.1111/j.1365-2966.2006.10040.x
  * Lim et al. (2020) Lim, C.-F., Wang, W.-H., Smail, I., et al. 2020, ApJ, 889, 80, doi: 10.3847/1538-4357/ab607f
  * Lotz et al. (2017) Lotz, J. M., Koekemoer, A., Coe, D., et al. 2017, ApJ, 837, 97, doi: 10.3847/1538-4357/837/1/97
  * Mairs et al. (2021) Mairs, S., Dempsey, J. T., Bell, G. S., et al. 2021, AJ, 162, 191, doi: 10.3847/1538-3881/ac18bf
  * McCully et al. (2014) McCully, C., Keeton, C. R., Wong, K. C., & Zabludoff, A. I. 2014, MNRAS, 443, 3631, doi: 10.1093/mnras/stu1316
  * Merten et al. (2011) Merten, J., Coe, D., Dupke, R., et al. 2011, MNRAS, 417, 333, doi: 10.1111/j.1365-2966.2011.19266.x
  * Mohammed et al. (2014) Mohammed, I., Liesenborgs, J., Saha, P., & Williams, L. L. R. 2014, MNRAS, 439, 2651, doi: 10.1093/mnras/stu124
  * Odegard et al. (2019) Odegard, N., Weiland, J. L., Fixsen, D. J., et al. 2019, ApJ, 877, 40, doi: 10.3847/1538-4357/ab14e8
  * Oguri (2010) Oguri, M. 2010, PASJ, 62, 1017, doi: 10.1093/pasj/62.4.1017
  * Oliver et al. (2010) Oliver, S. J., Wang, L., Smith, A. J., et al. 2010, A&A, 518, L21, doi: 10.1051/0004-6361/201014697
  * Puget et al. (1996) Puget, J. L., Abergel, A., Bernard, J. P., et al. 1996, A&A, 308, L5
  * Richard et al. (2014) Richard, J., Jauzac, M., Limousin, M., et al. 2014, MNRAS, 444, 268, doi: 10.1093/mnras/stu1395
  * Simpson et al. (2015) Simpson, J. M., Smail, I., Swinbank, A. M., et al. 2015, ApJ, 807, 128, doi: 10.1088/0004-637X/807/2/128
  * Smail et al. (1997) Smail, I., Ivison, R. J., & Blain, A. W. 1997, ApJ, 490, L5, doi: 10.1086/311017
  * Tadaki et al. (2020) Tadaki, K.-i., Belli, S., Burkert, A., et al. 2020, ApJ, 901, 74, doi: 10.3847/1538-4357/abaf4a
  * Wang et al. (2017) Wang, W.-H., Lin, W.-C., Lim, C.-F., et al. 2017, ApJ, 850, 37, doi: 10.3847/1538-4357/aa911b
  * Zavala et al. (2017) Zavala, J. A., Aretxaga, I., Geach, J. E., et al. 2017, MNRAS, 464, 3369, doi: 10.1093/mnras/stw2630
  * Zitrin et al. (2013) Zitrin, A., Meneghetti, M., Umetsu, K., et al. 2013, ApJ, 762, L30, doi: 10.1088/2041-8205/762/2/L30


$\displaystyle\inf_{x\in K}\mathbb{P}_{x}(n<\tau_{\partial})\geq\inf_{x\in
K}\mathbb{P}_{x}(n+k_{0}<\tau_{\partial})\geq
c_{1}c\,\mathbb{P}_{\nu}(n<\tau_{\partial})$
$\displaystyle\geq\frac{c_{1}c}{c_{2}}\sup_{y\in
E}\mathbb{P}_{y}(n<\tau_{\partial}).$

This implies (E3) and that $\inf_{x\in
K}\mathbb{P}_{x}(k_{0}<\tau_{\partial})>0$. Hence, (E1) follows from (A1) with
the probability measure $\frac{\nu(\cdot\cap K)}{\nu(K)}$. Moreover, for any
$n$ large enough to have $C\alpha^{n}\leq 1/2$ where the constants $C$ and
$\alpha$ are those of (3.6), we have $\mathbb{P}_{x}(X_{n}\in K\mid
t<\tau_{\partial})\geq\nu_{QSD}(K)-C\alpha^{n}\geq 1/2-\varepsilon>0$ and
hence (E4) is satisfied. The last computation also entails (3.7) with
$n^{\prime}_{4}=n$.

It remains to construct a function $\varphi_{1}$ satisfying (E2). For all
$x\in E$,

$\displaystyle\mathbb{P}_{x}(X_{n}\in E\setminus K\mid n<\tau_{\partial})$
$\displaystyle\leq\nu_{QSD}(E\setminus
K)+C\alpha^{n}\leq\varepsilon+C\alpha^{n}.$

Using (11.14), we deduce that

$\displaystyle\mathbb{P}_{x}(X_{n}\in E\setminus K)$ $\displaystyle\leq
C^{\prime}(\varepsilon+C\alpha^{n})\theta_{0}^{n},$

so that there exists $n_{0}$ large enough such that

$\displaystyle\mathbb{P}_{x}(n_{0}<T_{K}\wedge\tau_{\partial})$
$\displaystyle\leq\frac{1}{3}\theta_{0}^{n_{0}}=\left(\frac{\theta_{0}}{3^{1/n_{0}}}\right)^{n_{0}}.$

From this follows that, for all $k\in\mathbb{N}$ and all $x\in E$,

$\displaystyle\mathbb{P}_{x}(kn_{0}<T_{K}\wedge\tau_{\partial})$
$\displaystyle\leq\left(\frac{\theta_{0}}{3^{1/n_{0}}}\right)^{kn_{0}}.$

In particular, for $\theta_{1}:=\theta_{0}/2^{1/n_{0}}$,

$\displaystyle\varphi_{1}(x):=\mathbb{E}_{x}\left(\theta_{1}^{-T_{K}\wedge\lceil\tau_{\partial}\rceil}\right),\quad\forall
x\in E,$

is a bounded function on $E$ and Lemma 3.3 implies that, for all $x\in E$,

$\displaystyle
P_{1}\varphi_{1}(x)\leq\theta_{1}\varphi_{1}(x)+\|\varphi_{1}\|_{\infty}\mathbbm{1}_{K}(x).$

Since $\theta_{1}<\theta_{0}$, (E2) is proved.

##  12 Proof of the results of Section 4.1

In order to prove Theorem 4.1, we check Condition (F). The goal of Subsection
12.1 is to give the construction of the process $X$ and to check (F0) with
$L=K_{k}$ for some $k\geq 1$. In Subsection 12.2, we explain how (F3) can be
deduced from general Harnack inequalities. Finally, Subsection 12.3 completes
the proof of Theorem 4.1. The proof of Corollary 4.2 is then given in
Subsection 12.4.

###  12.1 Construction of the diffusion process $X$ and Markov property

The goal of this section is to construct a weak solution $X$ to the SDE (4.1)
with absorption out of $D$, and prove that it is Markov and satisfies a strong
Markov property at appropriate stopping times, enough to entail Condition (F0)
for $L=K_{k}$ for some $k\geq 1$. We introduce the natural path space for the
process $X$ as

$\mathcal{D}:=\Biggl{\\{}w:\mathbb{R}_{+}\rightarrow D\cup\\{\partial\\}:\
\forall k\geq 1,\ w\text{ is continuous on }[0,\tau_{k}(w)]\\\ \left.\text{
and }w(t)=\partial,\ \forall t\geq\sup_{k\geq 1}\tau_{k}(w)\right\\},$

where $\tau_{k}(w):=\inf\\{t\geq 0:w_{t}\in D\setminus K_{k}\\}$. Note that
$\mathcal{D}$ contains functions which are not càdlàg since they may not have
a left limit at $\tau_{\partial}-$ and, indeed, it is easy to construct
examples where $X$ is not càdlàg $\mathbb{P}$-a.s.222For example, one may
consider $D$ the open disc of radius 1 centered at 0 in $\mathbb{R}^{2}$,
$\sigma=\text{Id}$ and $b(x)=(-x_{2}\beta(|x|),x_{1}\beta(|x|))$ where
$x=(x_{1},x_{2})\in D$. Decomposing the process in polar coordinates
$(R_{t},\theta_{t}):=(|X_{t}|,\arctan(X^{(1)}_{t}/X^{(2)}_{t}))$, the radius
$R_{t}$ is a 2-dimensional Bessel process, and $X_{t}$ is sent to $\partial$
when $R_{t}$ hits 1 (in a.s. finite time). The angle $\theta_{t}$ is solution
to $d\theta_{t}=R_{t}^{-1}dW_{t}-\beta(R_{t})dt$ before $\tau_{\partial}$, for
some Brownian motion $W$. Hence, if $\beta(r)$ converges sufficiently fast to
$+\infty$ when $r\rightarrow 1$, $\theta_{t}$ a.s. converges to $-\infty$ when
$t\rightarrow\tau_{\partial}-$, so $X$ does not admit a left limit at time
$\tau_{\partial}$. Note also that this definition means that we are looking
for a process $X$ such that

$\displaystyle\tau_{\partial}:=\sup_{k\geq 1}\tau_{D\setminus K_{k}},$

which is the natural definition of $\tau_{\partial}$ when the left limit of
$X$ at time $\tau_{\partial}$ does not exist.

We endow the path space $\mathcal{D}$ with its natural filtration

$\mathcal{F}_{t}=\sigma(w_{s},s\leq t)=\bigvee_{n\geq 1,0\leq
t_{1}<t_{2}<\ldots<t_{n}\leq t}\sigma(w_{t_{1}},w_{t_{2}},\ldots,w_{t_{n}})$

and we follow the usual method which consists in constructing for all $x\in D$
a probability measure $\mathbb{P}_{x}$ on $\mathcal{D}$ and a stochastic
process $(B_{t},t\geq 0)$ on
$\mathcal{D}\times\mathcal{C}(\mathbb{R}_{+},\mathbb{R}^{r})$, such that $B$
is a standard $r$-dimensional Brownian motion under
$\mathbb{P}_{x}\otimes\mathbb{W}^{r}$, where $\mathbb{W}^{r}$ is the
$r$-dimensional Wiener measure and such that $w_{0}=x\ $
$\mathbb{P}_{x}\otimes\mathbb{W}^{r}$-almost surely and the canonical process
$(w_{t},t\geq 0)$ solves the SDE (4.1) for this Brownian motion $B$ on the
time interval $[0,\sup_{k}\tau_{k}(w))$ 333Since $\sigma(x)$ is non-degenerate
for all $x\in D$, the space $\mathcal{C}(\mathbb{R}_{+},\mathbb{R}^{r})$
equipped with the Wiener measure $\mathbb{W}^{r}$ is only used to construct
the Brownian path $B_{t}$ after time $\sup_{k}\tau_{k}(w)$ and could be
omitted for our purpose since we only need to construct the process $B$ up to
time $\sup_{k}\tau_{k}(w)$..

For this construction, we use the fact that $b$ and $\sigma$ can be extended
out of $K_{k}$ to $\mathbb{R}^{d}$ as globally Hölder and bounded functions
$b_{k}$ and $\sigma_{k}$ and such that $\sigma_{k}$ is uniformly elliptic on
$\mathbb{R}^{d}$. Hence (see e.g. [53, Rk. 5.4.30]) the martingale problem is
well-posed for the SDE

$\displaystyle dX^{k}_{t}=b_{k}(X^{k}_{t})dt+\sigma_{k}(X^{k}_{t})dB_{t}.$

Let us denote by $\mathbb{P}^{k}_{x}$ the solution to this martingale problem
for the initial condition $x\in\mathbb{R}^{d}$. This is a probability measure
on $\mathcal{C}:=\mathcal{C}(\mathbb{R}_{+},\mathbb{R}^{d})$, equipped with
its canonical filtration $(\mathcal{G}_{t})_{t\geq 0}$.

For all $k\geq 1$, we define $\tau^{\prime}_{k}(w)=\inf\\{t\geq
0,w_{t}\not\in\text{int}(K_{k})\\}$, where $\text{int}(K_{k})$ is the interior
of $K_{k}$. Since the paths $w\in\mathcal{D}$ or $\mathcal{C}$ are continuous
at time $\tau^{\prime}_{k}$ and $\mathbb{R}^{d}\setminus\text{int}(K_{k})$ is
closed, it is standard to prove that $\tau^{\prime}_{k}$ is a stopping time
for the canonical filtration $(\mathcal{F}_{t})_{t\geq 0}$ on $\mathcal{D}$
and for the canonical filtration $(\mathcal{G}_{t})_{t\geq 0}$ on
$\mathcal{C}$. We define as usual the stopped $\sigma$-fields
$\mathcal{F}_{\tau^{\prime}_{k}}$ and $\mathcal{G}_{\tau^{\prime}_{k}}$, and
we define for all $x\in\text{int}(K_{k})$ the restriction of $\mathbb{P}_{x}$
to $\mathcal{F}_{\tau^{\prime}_{k}}$ as the restriction of
$\mathbb{P}^{k}_{x}$ to $\mathcal{G}_{\tau^{\prime}_{k}}$, where we can
identify the events of the two filtrations since they both concern continuous
parts of the paths. This construction is consistent for $k$ and $k+1$ (meaning
that if $x\in K_{k}$, they give the same probability to events of
$\mathcal{F}_{\tau_{k}}$) by uniqueness of the solutions $\mathbb{P}^{k}_{x}$
and $\mathbb{P}^{k+1}_{x}$ to the above martingale problems. Hence there
exists a unique extension $\mathbb{P}_{x}$ of the above measures to
$\bigvee_{k\geq 1}\mathcal{F}_{\tau^{\prime}_{k}}$. Note that, because of the
specific structure of the path space $\mathcal{D}$, we have

$\bigvee_{k\geq 1}\mathcal{F}_{\tau^{\prime}_{k}}=\mathcal{F}_{\infty}.$
(12.1)

To check this, it suffices to observe that, for all $t\geq 0$ and all
measurable $A\subset D\cup\\{\partial\\}$,

$\displaystyle\\{w_{t}\in A\\}$ $\displaystyle=\\{t<\tau_{\partial},\ w_{t}\in
A\cap D\\}\cup\\{\tau_{\partial}\leq t,\ \partial\in A\\}$
$\displaystyle=\left(\bigcup_{k\geq 1}\\{t<\tau^{\prime}_{k},\ w_{t}\in A\cap
D\\}\right)\cup\left(\bigcap_{k\geq 1}\\{\tau^{\prime}_{k}\leq t,\partial\in
A\\}\right),$ (12.2)

hence $\\{w_{t}\in A\\}\in\bigvee_{k\geq 1}\mathcal{F}_{\tau^{\prime}_{k}}$,
and, proceeding similarly, the same property holds for events of the form
$\\{w_{t_{1}}\in A_{1},\ldots,w_{t_{n}}\in A_{n}\\}$.

We recall (see [53, Section 5.4]) that
$(\mathbb{P}^{k}_{x})_{x\in\mathbb{R}^{d}}$ forms a strong Markov family on
the canonical space $\mathcal{C}$. Our goal is now to prove that the family of
probability measures $(\mathbb{P}_{x})_{x\in D\cup\\{\partial\\}}$, where
$\mathbb{P}_{\partial}$ is defined as the Dirac measure on the constant path
equal to $\partial$, forms a Markov kernel of probability measures, for which
the strong Markov property applies at well-chosen stopping times.

We first need to prove that $(\mathbb{P}_{x})_{x\in D}$ defines a kernel of
probability measures, i.e. that $x\mapsto\mathbb{P}_{x}(\Gamma)$ is measurable
for all events $\Gamma$ of $\mathcal{F}_{\infty}$. We prove it for an event of
the form $\\{w_{t}\in A\\}$, the extension to events of the form
$\\{w_{t_{1}}\in A_{1},\ldots,w_{t_{n}}\in A_{n}\\}$, and hence to all events
of $\mathcal{F}_{\infty}$, being easy. This follows from (12.2):

$\displaystyle\mathbb{P}_{x}(w_{t}\in A)$
$\displaystyle=\lim_{k\rightarrow+\infty}\mathbb{P}_{x}(t<\tau^{\prime}_{k},\
w_{t}\in A\cap D)+\mathbbm{1}_{\partial\in
A}\lim_{k\rightarrow+\infty}\mathbb{P}_{x}(\tau^{\prime}_{k}\leq t)$
$\displaystyle=\lim_{k\rightarrow+\infty}\mathbb{P}^{k+1}_{x}(t<\tau^{\prime}_{k},\
w_{t}\in A\cap D)+\mathbbm{1}_{\partial\in
A}\lim_{k\rightarrow+\infty}\mathbb{P}^{k+1}_{x}(\tau^{\prime}_{k}\leq t).$

Since all the probabilities in the right-hand side are measurable functions of
$x$, so is $x\mapsto\mathbb{P}_{x}(w_{t}\in A)$.

Now, let us prove that $(X_{t},t\geq 0)$ is Markov. It is well-known that this
is implied by the following property: for all $n\geq 1$ and $0\leq
t_{1}\leq\ldots\leq t_{n+1}$ and $A_{1},\ldots,A_{n+1}$ measurable subsets of
$D\cup\\{\partial\\}$,

$\displaystyle\mathbb{P}_{x}(w_{t_{1}}\in A_{1},\ldots,w_{t_{n+1}}\in
A_{n+1})=\mathbb{E}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in
A_{1},\ldots,w_{t_{n}}\in A_{n}}\mathbb{P}_{w_{t_{n}}}(w_{t_{n+1}-t_{n}}\in
A_{n+1})\right].$

We prove this property only for $n=1$. It is easy to extend the proof to all
values of $n\geq 1$. We have

$\mathbb{P}_{x}(w_{t_{1}}\in A_{1},w_{t_{2}}\in
A_{2})=\mathbb{P}_{x}(w_{t_{1}}\in A_{1},w_{t_{2}}\in
A_{2},\tau_{\partial}>t_{2})\\\ +\mathbb{P}_{x}(w_{t_{1}}\in
A_{1},t_{1}<\tau_{\partial}\leq t_{2})\mathbbm{1}_{\partial\in
A_{2}}+\mathbb{P}_{x}(\tau_{\partial}\leq t_{1})\mathbbm{1}_{\partial\in
A_{1}\cap A_{2}}.$

Now, using that $(\mathbb{P}^{k}_{x})_{x\in\mathbb{R}^{d}}$ is a Markov family
for all $k\geq 1$,

$\mathbb{P}_{x}(w_{t_{1}}\in A_{1},w_{t_{2}}\in
A_{2},\tau_{\partial}>t_{2})\\\ \begin{aligned}
&=\lim_{k\rightarrow\infty}\mathbb{P}_{x}(w_{t_{1}}\in A_{1},w_{t_{2}}\in
A_{2},\tau_{k}>t_{2})\\\
&=\lim_{k\rightarrow\infty}\mathbb{P}^{k}_{x}(w_{t_{1}}\in A_{1},w_{t_{2}}\in
A_{2},\tau_{k}>t_{2})\\\
&=\lim_{k\rightarrow\infty}\mathbb{E}^{k}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in
A_{1},t_{1}<\tau_{k}}\mathbb{P}^{k}_{w_{t_{1}}}(w_{t_{2}-t_{1}}\in
A_{2},\tau_{k}>t_{2}-t_{1})\right]\\\
&=\lim_{k\rightarrow\infty}\mathbb{E}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in
A_{1},t_{1}<\tau_{k}}\mathbb{P}_{w_{t_{1}}}(w_{t_{2}-t_{1}}\in
A_{2},\tau_{k}>t_{2}-t_{1})\right]\\\
&=\mathbb{E}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in
A_{1},t_{1}<\tau_{\partial}}\mathbb{P}_{w_{t_{1}}}(w_{t_{2}-t_{1}}\in
A_{2},\tau_{\partial}>t_{2}-t_{1})\right]\end{aligned}$

and similarly

$\displaystyle\mathbb{P}_{x}(w_{t_{1}}\in A_{1},t_{1}<\tau_{\partial}\leq
t_{2})\mathbbm{1}_{\partial\in A_{2}}$
$\displaystyle=\mathbb{E}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in
A_{1},t_{1}<\tau_{\partial}}\mathbb{P}_{w_{t_{1}}}(\tau_{\partial}\leq
t_{2}-t_{1})\right]\mathbbm{1}_{\partial\in A_{2}}$
$\displaystyle=\mathbb{E}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in
A_{1},t_{1}<\tau_{\partial}}\mathbb{P}_{w_{t_{1}}}(\tau_{\partial}\leq
t_{2}-t_{1},w_{t_{2}-t_{1}}\in A_{2})\right].$

Since

$\displaystyle\mathbb{P}_{x}(\tau_{\partial}\leq
t_{1})\mathbbm{1}_{\partial\in A_{1}\cap A_{2}}$
$\displaystyle=\mathbb{E}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in
A_{1},\tau_{\partial}\leq t_{1}}\mathbb{P}_{w_{t_{1}}}(w_{t_{2}-t_{1}}\in
A_{2})\right],$

we have proved that $\mathbb{P}_{x}(w_{t_{1}}\in A_{1},w_{t_{2}}\in
A_{2})=\mathbb{E}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in
A_{1}}\mathbb{P}_{w_{t_{1}}}(w_{t_{2}-t_{1}}\in A_{2})\right]$. This ends the
proof of the Markov property.

To conclude this subsection, let us prove that the strong Markov property
holds for all stopping times $\tau_{F}$ where $F\subset D$ is closed in $D$.
Note that $\tau_{F}$ is indeed a stopping time for the filtration
$\mathcal{F}_{t}$ since
$\tau_{F}=\sup_{k}\tau_{F}\wedge\tau^{\prime}_{k}=\sup_{k}\tau_{(F\cup
D^{c})\cup\text{int}(K_{k})^{c}}$, where the complement is understood in
$\mathbb{R}^{d}$, $(F\cup D^{c})\cup\text{int}(K_{k})^{c}$ is a closed subset
of $\mathbb{R}^{d}$ and all $w\in\mathcal{D}$ is continuous at time
$\tau_{(F\cup D^{c})\cup\text{int}(K_{k})^{c}}$. Let $x\in D$,
$t_{1},t_{2},s\geq 0$ and $A,B\subset D$ be measurable sets. We proceed as
above: first, observe that

$\\{w_{t_{1}}\in A,\ t_{1}<\tau_{F}\leq t_{2},\ w_{\tau_{F}+s}\in B\\}\\\
\begin{aligned} &=\bigcup_{\ell\geq 1}\\{w_{t_{1}}\in A,\ t_{1}<\tau_{F}\leq
t_{2},\ w_{\tau_{F}+s}\in B,\ w_{r}\in K_{\ell}\ \forall
r\in[0,\tau_{F}+s]\\}\\\ &=\bigcup_{\ell\geq 1}\\{w_{t_{1}}\in A,\
t_{1}<\tau_{F}\wedge\tau^{\prime}_{\ell}\leq t_{2},\
w_{\tau_{F}\wedge\tau^{\prime}_{\ell}+s}\in B,\
\tau^{\prime}_{\ell}>\tau_{F}+s\\}.\end{aligned}$

Since $\tau_{F}\wedge\tau^{\prime}_{\ell}$ is a $\mathcal{G}_{t}$-stopping
time on $\mathcal{C}(\mathbb{R}_{+},\mathbb{R}^{d})$ and using the strong
Markov property under $\mathbb{P}^{\ell}$, we deduce that

$\mathbb{P}_{x}(w_{t_{1}}\in A,\ t_{1}<\tau_{F}\leq t_{2},\ w_{\tau_{F}+s}\in
B)\\\ \begin{aligned}
&=\lim_{\ell\rightarrow+\infty}\mathbb{P}^{\ell}_{x}(w_{t_{1}}\in A,\
t_{1}<\tau_{F}\wedge\tau^{\prime}_{\ell}\leq t_{2},\
w_{\tau_{F}\wedge\tau^{\prime}_{\ell}+s}\in B,\
\tau^{\prime}_{\ell}>\tau_{F}+s)\\\
&=\lim_{\ell\rightarrow+\infty}\mathbb{E}^{\ell}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in
A,\ t_{1}<\tau_{F}\wedge\tau^{\prime}_{\ell}\leq
t_{2}}\mathbb{P}^{\ell}_{w_{\tau_{F}\wedge\tau^{\prime}_{\ell}}}(w_{s}\in B,\
s<\tau^{\prime}_{\ell})\right]\\\
&=\lim_{\ell\rightarrow+\infty}\mathbb{E}^{\ell}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in
A,\ t_{1}<\tau_{F}\leq\tau^{\prime}_{\ell}\wedge
t_{2}}\mathbb{P}^{\ell}_{w_{\tau_{F}}}(w_{s}\in B,\
s<\tau^{\prime}_{\ell})\right]\\\
&=\mathbb{E}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in A,\
t_{1}<\tau_{F}\leq\tau_{\partial}\wedge
t_{2}}\mathbb{P}_{w_{\tau_{F}}}(w_{s}\in B,\
s<\tau_{\partial})\right].\end{aligned}$

Similarly,

$\displaystyle\mathbb{P}_{x}(w_{t_{1}}\in A,\ t_{1}<\tau_{F}\leq t_{2},\ $
$\displaystyle w_{\tau_{F}+s}=\partial)$
$\displaystyle=\lim_{\ell\rightarrow+\infty}\mathbb{P}^{\ell}_{x}(w_{t_{1}}\in
A,\ t_{1}<\tau_{F}\leq t_{2}\wedge\tau^{\prime}_{\ell},\
\tau^{\prime}_{\ell}\leq\tau_{F}+s)$
$\displaystyle=\mathbb{E}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in A,\
t_{1}<\tau_{F}\leq
t_{2}\wedge\tau_{\partial}}\mathbb{P}_{w_{\tau_{F}}}(w_{s}=\partial)\right]$

and thus

$\mathbb{P}_{x}(w_{t_{1}}\in A,\ t_{1}<\tau_{F}\leq t_{2},\ w_{\tau_{F}+s}\in
B)=\mathbb{E}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in A,\ t_{1}<\tau_{F}\leq
t_{2}\wedge\tau_{\partial}}\mathbb{P}_{w_{\tau_{F}}}(w_{s}\in B)\right]$

for all $A,B\subset D\cup\\{\partial\\}$ measurable. The previous computation
extends without difficulty to prove

$\mathbb{P}_{x}\left(w_{t_{1}}\in A_{1},\ldots,w_{t_{n}}\in A_{n},\
t_{n}<\tau_{F}\leq t_{n+1},\ w_{\tau_{F}+s_{1}}\in
B_{1},\ldots,w_{\tau_{F}+s_{m}}\in B_{m}\right)\\\
=\mathbb{E}_{x}\left[\mathbbm{1}_{w_{t_{1}}\in A_{1},\ldots,w_{t_{n}}\in
A_{n},\ t_{n}<\tau_{F}\leq t_{n+1}}\mathbb{P}_{w_{\tau_{F}}}(w_{s_{1}}\in
B_{1},\ldots,w_{s_{m}}\in B_{m})\right]$ (12.3)

for all $n,m\geq 1$, $0\leq t_{1}\leq\ldots\leq t_{n+1}$, $0\leq
s_{1}\leq\ldots\leq s_{m}$ and $A_{1},\ldots,A_{n},B_{1},\ldots,B_{m}\subset
D\cup\\{\partial\\}$ measurable. This implies the strong Markov property at
time $\tau_{F}$, in the sense that, for all $k\geq 1$, all $x\in E$ and all
$\Gamma\in\mathcal{F}_{\infty}$,

$\displaystyle\mathbb{P}_{x}\left(w^{\tau_{F}}\in\Gamma\mid\mathcal{H}_{\tau_{F}}\right)=\mathbb{P}_{w_{\tau_{F}}}(\Gamma),\quad\mathbb{P}_{x}\text{-almost
surely},$

where $w^{\tau_{F}}=(w_{\tau_{F}+s},s\geq 0)$ and

$\mathcal{H}_{\tau_{F}}=\sigma\Big{(}\left\\{w_{t_{1}}\in
A_{1},\ldots,w_{t_{n}}\in A_{n},t_{n}<\tau_{F}\leq t_{n+1}\right\\},\\\ 0\leq
t_{1}\leq\ldots\leq t_{n+1},\ A_{1},\ldots,A_{n}\in D\text{
measurable}\Big{)}.$

This form of strong Markov property at time $\tau_{F}$ is enough for our
purpose, since it entails (F0) for $L=K_{k}$ for all $k\geq 1$. It will be
also needed in the next section.

###  12.2 Harnack inequalities

Our goal here is to check Conditions (F1) and (F3) for the diffusion process
constructed above. We will make use of general Harnack inequalities of Krylov
and Safonov [59].

######  Proposition 12.1.

There exist a probability measure $\nu$ on $D$ and a constant $t_{\nu}>0$ such
that, for all $k\geq 1$, there exists a constant $b_{k}>0$ such that

$\displaystyle\mathbb{P}_{x}(X_{t_{\nu}}\in\cdot)\geq b_{k}\nu(\cdot),\
\forall x\in K_{k}.$ (12.4)

Moreover, for all $k\geq 1$ such that $K_{k}$ is non-empty,

$\displaystyle\inf_{t\geq 0}\frac{\inf_{x\in
K_{k}}\mathbb{P}_{x}(t<\tau_{\partial})}{\sup_{x\in
K_{k}}\mathbb{P}_{x}(t<\tau_{\partial})}>0.$ (12.5)

###### Proof.

Consider a bounded measurable function $f:D\rightarrow\mathbb{R}$ with
$\|f\|_{\infty}\leq 1$ and define the application
$u:(t,x)\in\mathbb{R}_{+}\times
E\mapsto\mathbb{E}_{x}[\mathbbm{1}_{t<\tau_{\partial}}f(X_{t})]$. It is proved
in [18] using [59] that, for all $k\geq 1$, there exist two constants
$N_{k}>0$ and $\delta_{k}>0$, which do not depend on $f$ (provided
$\|f\|_{\infty}\leq 1$), such that

$\displaystyle u(\delta_{k}+\delta_{k}^{2},x)\leq
N_{k}u(\delta_{k}+2\delta_{k}^{2},y),\text{ for all $x,y\in K_{k}$ such that
}|x-y|\leq\delta_{k}.$ (12.6)

Note that the proof given in [18] makes use of the following strong Markov
property: for all open ball $B$ such that $B\subset K_{k}$ for some $k\geq 1$,
all $x\in B$, $t\geq 0$ and all measurable
$f:D\cup\\{\partial\\}\rightarrow\mathbb{R}_{+}$,

$\mathbb{E}_{x}\left[f(X_{t})\mathbbm{1}_{\tau_{D\setminus B}\leq
t}\right]=\mathbb{E}_{x}\left[\mathbbm{1}_{\tau_{D\setminus B}\leq
t}\mathchoice{{\mathbb{E}_{X_{\tau_{D\setminus
B}}}\left[f(X_{t-u})\right]\,\smash{\vrule
height=6.00002pt,depth=3.13307pt}}_{\,u=\tau_{D\setminus
B}}}{{\mathbb{E}_{X_{\tau_{D\setminus
B}}}\left[f(X_{t-u})\right]\,\smash{\vrule
height=6.00002pt,depth=3.13307pt}}_{\,u=\tau_{D\setminus
B}}}{{\mathbb{E}_{X_{\tau_{D\setminus
B}}}\left[f(X_{t-u})\right]\,\smash{\vrule
height=4.20001pt,depth=2.6684pt}}_{\,u=\tau_{D\setminus
B}}}{{\mathbb{E}_{X_{\tau_{D\setminus
B}}}\left[f(X_{t-u})\right]\,\smash{\vrule
height=3.0pt,depth=2.6684pt}}_{\,u=\tau_{D\setminus B}}}\right].$

This property follows from (12.3).

Step 1 : Proof of (12.4)

Fix $x_{1}\in D$ and $k_{1}\geq 1$ such that $x_{1}\in\text{int}(K_{k_{1}})$.
Let $\nu$ denote the conditional law
$\mathbb{P}_{x_{1}}(X_{\delta_{k_{1}}+\delta_{k_{1}}^{2}}\in\cdot\mid\delta_{k_{1}}+\delta_{k_{1}}^{2}<\tau_{\partial})$.
Then, for all measurable $A\subset D\cup\\{\partial\\}$, Harnack’s inequality
(12.6) with $f=\mathbbm{1}_{A}$ entails that, for all $x\in D$ such that
$|x-x_{1}|<\delta_{k_{1}}\wedge d(x_{1},D\setminus K_{k_{1}})$,

$\displaystyle\mathbb{P}_{x}(X_{\delta_{k_{1}}+2\delta_{k_{1}}^{2}}\in
A)\geq\frac{\mathbb{P}_{x_{1}}(\delta_{k_{1}}+\delta_{k_{1}}^{2}<\tau_{\partial})}{N_{k_{1}}}\,\nu(A).$

Since the diffusion is locally elliptic and $D$ is connected, for all $k\geq
1$, there exists a constant $d_{k}>0$ such that

$\displaystyle\inf_{x\in K_{k}}\mathbb{P}_{x}(X_{1}\in
B(x_{1},\delta_{k_{1}}\wedge d(x_{1},D\setminus K_{k_{1}}))\geq d_{k}.$

This and Markov’s property entail that, for all $x\in K_{k}$,

$\displaystyle\mathbb{P}_{x}(X_{1+\delta_{k_{1}}+2\delta_{k_{1}}^{2}}\in\cdot)\geq
d_{k}\frac{\mathbb{P}_{x_{1}}(\delta_{k_{1}}+\delta_{k_{1}}^{2}<\tau_{\partial})}{N_{k_{1}}}\,\nu.$

This implies the first part of Proposition 12.1.

Step 2 : Proof of (12.5)

Fix $k\geq 1$ such that $K_{k}$ is non-empty and consider $\ell>k$ such that
$K_{k}$ is included in one connected component of $\text{int}(K_{\ell})$. For
all $t\geq\delta_{\ell}+2\delta_{\ell}^{2}$, the inequality (12.6) applied to
$f(x)=\mathbb{P}_{x}(t-\delta_{\ell}-2\delta_{\ell}^{2}<\tau_{\partial})$ and
the Markov property entail that

$\displaystyle\mathbb{P}_{x}(t-\delta_{\ell}^{2}<\tau_{\partial})\leq
N_{\ell}\mathbb{P}_{y}(t<\tau_{\partial}),\text{ for all $x,y\in K_{\ell}$
such that }|x-y|\leq\delta_{\ell}.$

Since $s\mapsto\mathbb{P}_{x}(s<\tau_{\partial})$ is non-increasing, we deduce
that

$\displaystyle\mathbb{P}_{x}(t<\tau_{\partial})\leq
N_{\ell}\mathbb{P}_{y}(t<\tau_{\partial}),\text{ for all $x,y\in K_{\ell}$
such that }|x-y|\leq\delta_{\ell}.$

Since $K_{k}$ has a finite diameter and is included in a connected component
of $K_{\ell}$, we deduce that there exists $N^{\prime}_{k}$ equal to some
power of $N_{\ell}$ such that, for all
$t\geq\delta_{\ell}+2\delta_{\ell}^{2}$,

$\displaystyle\mathbb{P}_{x}(t<\tau_{\partial})\leq
N^{\prime}_{k}\mathbb{P}_{y}(t<\tau_{\partial}),\text{ for all $x,y\in
K_{k}$.}$

Now, for $t\leq\delta_{\ell}+2\delta_{\ell}^{2}$, we simply use the fact that
$x\mapsto\mathbb{P}_{x}(\delta_{\ell}+2\delta_{\ell}^{2}<\tau_{\partial})$ is
uniformly bounded from below on $K_{k}$ by a constant
$1/N^{\prime\prime}_{k}>0$. In particular,

$\displaystyle\mathbb{P}_{x}(t<\tau_{\partial})\leq 1\leq
N^{\prime\prime}_{k}\mathbb{P}_{y}(\delta_{\ell}+2\delta_{\ell}^{2}<\tau_{\partial})\leq
N^{\prime\prime}_{k}\mathbb{P}_{y}(t<\tau_{\partial}),\text{ for all $x,y\in
K_{k}$.}$

This concludes the proof of Proposition 12.1. ∎

###  12.3 Proof of Theorem 4.1

Our aim is to prove that Condition (F) holds true with $L=K_{k}$ for some
$k\geq 1$. We have already proved (F0), (F1) and (F3) with $L=K_{k}$ for any
$k\geq 1$. Hence we only have to check (F2). Fix
$\rho_{1}\in(\lambda_{0},\lambda_{1})$, $\rho_{2}\in(\lambda_{0},\rho_{1})$
and $p\in(1,\lambda_{1}/\rho_{1})$ and define

$\psi_{1}(x)=\varphi(x)^{1/p},\ \forall x\in D.$

Fix $\rho^{\prime}_{1}\in(\rho_{1},\lambda_{1}/p)$ and

$t_{2}\geq\frac{2s_{1}(C+\lambda_{1})}{\lambda_{1}-p\rho^{\prime}_{1}}\,\vee\,\frac{\log
2}{\rho^{\prime}_{1}-\rho_{1}},$

where the constant $C$ comes from (4.5). Set $L=K_{k_{0}}$ with $k_{0}$ large
enough so that $\nu(K_{k_{0}})>0$ and, using (4.6),

$\mathbb{P}_{x}(s_{1}<\tau_{K_{k_{0}}}\wedge\tau_{\partial})\leq
e^{-(\rho^{\prime}_{1}+C/p)t_{2}}$

for all $x\in D_{0}$.

From the definition of $\lambda_{0}$ and applying the same argument as in Step
2 of the proof of Proposition 12.1 with
$f(x)=\mathbb{P}_{x}(X_{t-\delta_{\ell}-2\delta_{\ell}^{2}}\in L)$ with $\ell$
large enough to have $K_{k_{0}}$ included in one connected component of
$K_{\ell}$, we deduce that

$\displaystyle\liminf_{t\rightarrow+\infty}e^{\rho_{2}t}\inf_{x\in
L}\mathbb{P}_{x}(X_{t}\in L)=+\infty,$

and hence the last line of (F2) is proved with $\gamma_{2}=e^{-\rho_{2}}$.

Let us now check that the first line of Assumption (F2) holds true for all
$x\in D_{0}$ and then for all $x\in D\setminus D_{0}$. For all $x\in D_{0}$,
we have $\psi_{1}(x)\leq\sup_{x\in D_{0}}\\\ varphi^{1/p}(x)<+\infty$, and
hence, for all $t\in[s_{1},t_{2}]$, using Hölder’s inequality and the
definition of $k_{0}$,

$\displaystyle\mathbb{E}_{x}\left(\psi_{1}(X_{t})\mathbbm{1}_{t<\tau_{L}\wedge\tau_{\partial}}\right)$
$\displaystyle\leq\mathbb{E}_{x}\left(\mathbbm{1}_{t<\tau_{\partial}}\varphi(X_{t})\right)^{1/p}\mathbb{P}_{x}(t<\tau_{L}\wedge\tau_{\partial})^{\frac{p-1}{p}}$
$\displaystyle\leq\varphi(x)^{1/p}e^{Ct_{2}/p}\mathbb{P}_{x}(s_{1}<\tau_{L}\wedge\tau_{\partial})^{\frac{p-1}{p}}$
(12.7) $\displaystyle\leq e^{-\rho^{\prime}_{1}t_{2}}\leq
e^{-\rho_{1}t_{2}}\psi_{1}(x).$

To prove (12.7), we used the fact that $\mathcal{L}\varphi\leq C\leq C\varphi$
and Itô’s formula to obtain $P_{t}\varphi\leq e^{Ct}\varphi$. Since this
argument is used repeatedly in the sequel, we give it in details for sake of
completeness. It follows from Itô’s formula that, for all $k\geq 1$,
$\mathbb{P}_{x}$-almost surely,

$\displaystyle
e^{-C\left(t\wedge\tau_{K_{k}^{c}}\right)}\varphi\left(X_{t\wedge\tau_{K_{k}^{c}}}\right)=\varphi(x)$
$\displaystyle+\int_{0}^{t}\mathbbm{1}_{s\leq\tau_{K_{k}^{c}}}e^{-Cs}\left(\mathcal{L}\varphi(X_{s})-C\varphi(X_{s})\right)ds$
$\displaystyle+\int_{0}^{t}\mathbbm{1}_{s\leq\tau_{K_{k}^{c}}}e^{-Cs}\nabla\varphi(X_{s})^{*}\sigma(X_{s})dB_{s}.$

Since $\nabla\varphi(x)$ and $\sigma(x)$ are uniformly bounded on $K_{k}$, the
last term has zero expectation, and thus

$\displaystyle\mathbb{E}_{x}\left[e^{-C\left(t\wedge\tau_{K_{k}^{c}}\right)}\varphi\left(X_{t\wedge\tau_{K_{k}^{c}}}\right)\right]$
$\displaystyle\leq\varphi(x).$

Letting $k\rightarrow+\infty$, we deduce form Fatou’s lemma that

$\displaystyle\mathbb{E}_{x}\left[e^{-Ct}\mathbbm{1}_{t<\tau_{\partial}}\varphi(X_{t})\right]$
$\displaystyle\leq\varphi(x)$ (12.8)

as claimed.

This proves the second line of (F2) for all $x\in D_{0}$ and
$\gamma_{1}=e^{-\rho_{1}}$.

Now, for all $x\in D\setminus D_{0}$, since $D_{0}$ is closed in $D$, it
follows from the strong Markov property (12.3) at time $\tau_{D_{0}}$ that

$\mathbb{E}_{x}\left(\psi_{1}(X_{t_{2}})\mathbbm{1}_{t_{2}<\tau_{L}\wedge\tau_{\partial}}\right)=\mathbb{E}_{x}\left(\mathbbm{1}_{t_{2}-s_{1}<\tau_{L}\wedge\tau_{\partial}\wedge\tau_{D_{0}}}\mathbb{E}_{X_{t_{2}-s_{1}}}\left(\psi_{1}(X_{s_{1}})\mathbbm{1}_{s_{1}<\tau_{L}\wedge\tau_{\partial}}\right)\right)\\\
+\mathbb{E}_{x}\left(\mathbbm{1}_{\tau_{D_{0}}\leq
t_{2}-s_{1}}\mathchoice{{\mathbb{E}_{X_{\tau_{D_{0}}}}\left(\psi_{1}(X_{t_{2}-u})\mathbbm{1}_{t_{2}-u<\tau_{\partial}\wedge\tau_{L}}\right)\,\smash{\vrule
height=6.00002pt,depth=3.5163pt}}_{\,u=\tau_{D_{0}}}}{{\mathbb{E}_{X_{\tau_{D_{0}}}}\left(\psi_{1}(X_{t_{2}-u})\mathbbm{1}_{t_{2}-u<\tau_{\partial}\wedge\tau_{L}}\right)\,\smash{\vrule
height=6.00002pt,depth=3.5163pt}}_{\,u=\tau_{D_{0}}}}{{\mathbb{E}_{X_{\tau_{D_{0}}}}\left(\psi_{1}(X_{t_{2}-u})\mathbbm{1}_{t_{2}-u<\tau_{\partial}\wedge\tau_{L}}\right)\,\smash{\vrule
height=4.20001pt,depth=3.05162pt}}_{\,u=\tau_{D_{0}}}}{{\mathbb{E}_{X_{\tau_{D_{0}}}}\left(\psi_{1}(X_{t_{2}-u})\mathbbm{1}_{t_{2}-u<\tau_{\partial}\wedge\tau_{L}}\right)\,\smash{\vrule
height=3.0pt,depth=3.05162pt}}_{\,u=\tau_{D_{0}}}}\right).$ (12.9)

Using Hölder’s inequality and (12.8), we deduce that, for all $y\in D$,

$\displaystyle\mathbb{E}_{y}\left(\psi_{1}(X_{s_{1}})\mathbbm{1}_{s_{1}<\tau_{L}\wedge\tau_{\partial}}\right)\leq\mathbb{E}_{y}\left(\varphi(X_{s_{1}})\mathbbm{1}_{s_{1}<\tau_{\partial}}\right)^{1/p}\leq
e^{\frac{s_{1}C}{p}}\,\varphi(y)^{1/p}=e^{\frac{s_{1}C}{p}}\,\psi_{1}(y).$

Hence, the first term in the right-hand side of (12.9) satisfies

$\displaystyle\mathbb{E}_{x}\left(\mathbbm{1}_{t_{2}-s_{1}<\tau_{L}\wedge\tau_{\partial}\wedge\tau_{D_{0}}}\mathbb{E}_{X_{t_{2}-s_{1}}}\left(\psi_{1}(X_{s_{1}})\mathbbm{1}_{s_{1}<\tau_{L}\wedge\tau_{\partial}}\right)\right)\leq
e^{\frac{s_{1}C}{p}}\mathbb{E}_{x}\left(\mathbbm{1}_{t_{2}-s_{1}<\tau_{L}\wedge\tau_{\partial}\wedge\tau_{D_{0}}}\psi_{1}(X_{t_{2}-s_{1}})\right).$

As a consequence, using again Hölder’s inequality and applying as above Itô’s
formula using that $\mathcal{L}\varphi(x)\leq-\lambda_{1}\varphi(x)$ for all
$x\notin D_{0}$, one has

$\displaystyle\mathbb{E}_{x}\left(\mathbbm{1}_{t_{2}-s_{1}<\tau_{L}\wedge\tau_{\partial}\wedge\tau_{D_{0}}}\mathbb{E}_{X_{t_{2}-s_{1}}}\left(\psi_{1}(X_{s_{1}})\mathbbm{1}_{s_{1}<\tau_{L}\wedge\tau_{\partial}}\right)\right)$
$\displaystyle\leq
e^{-\lambda_{1}\frac{t_{2}-s_{1}}{p}}e^{\frac{s_{1}C}{p}}\,\varphi(x)^{1/p}$
$\displaystyle\leq
e^{-t_{2}\frac{\rho^{\prime}_{1}+\lambda_{1}/p}{2}}\psi_{1}(x),$

where we used in the last inequality that
$t_{2}\geq\frac{2s_{1}(C+\lambda_{1})}{\lambda_{1}-p\rho^{\prime}_{1}}$.
Moreover, using (12.7), we obtain that the second term in the right-hand side
of (12.9) satisfies

$\mathbb{E}_{x}\left(\mathbbm{1}_{\tau_{D_{0}}\leq
t_{2}-s_{1}}\mathchoice{{\mathbb{E}_{X_{\tau_{D_{0}}}}\left(\psi_{1}(X_{t_{2}-u})\mathbbm{1}_{t_{2}-u<\tau_{\partial}\wedge\tau_{L}}\right)\,\smash{\vrule
height=6.00002pt,depth=3.5163pt}}_{\,u=\tau_{D_{0}}}}{{\mathbb{E}_{X_{\tau_{D_{0}}}}\left(\psi_{1}(X_{t_{2}-u})\mathbbm{1}_{t_{2}-u<\tau_{\partial}\wedge\tau_{L}}\right)\,\smash{\vrule
height=6.00002pt,depth=3.5163pt}}_{\,u=\tau_{D_{0}}}}{{\mathbb{E}_{X_{\tau_{D_{0}}}}\left(\psi_{1}(X_{t_{2}-u})\mathbbm{1}_{t_{2}-u<\tau_{\partial}\wedge\tau_{L}}\right)\,\smash{\vrule
height=4.20001pt,depth=3.05162pt}}_{\,u=\tau_{D_{0}}}}{{\mathbb{E}_{X_{\tau_{D_{0}}}}\left(\psi_{1}(X_{t_{2}-u})\mathbbm{1}_{t_{2}-u<\tau_{\partial}\wedge\tau_{L}}\right)\,\smash{\vrule
height=3.0pt,depth=3.05162pt}}_{\,u=\tau_{D_{0}}}}\right)\\\ \leq
e^{-\rho^{\prime}_{1}t_{2}}\mathbb{P}_{x}(\tau_{D_{0}}\leq t_{2}-s_{1})\leq
e^{-\rho^{\prime}_{1}t_{2}}\psi_{1}(x).$

We finally deduce from (12.9) and from the definition of $L=K_{k_{0}}$ that

$\displaystyle\mathbb{E}_{x}\left(\psi_{1}(X_{t_{2}})\mathbbm{1}_{t_{2}<\tau_{L}\wedge\tau_{\partial}}\right)\leq
2e^{-\rho^{\prime}_{1}t_{2}}\psi_{1}(x)\leq e^{-\rho_{1}t_{2}}\psi_{1}(x),$

where we used that $t_{2}\geq\log 2/(\rho^{\prime}_{1}-\rho_{1})$. This
concludes the proof that the second line of (F2) holds true.

Since $\varphi$ is locally bounded, $\sup_{L}\varphi<\infty$, and hence, using
again (12.8), we deduce that, for all $t\geq 0$,

$\displaystyle\sup_{x\in L}\
\mathbb{E}_{x}(\psi_{1}(X_{t})\mathbbm{1}_{t<\tau_{\partial}})\leq\sup_{x\in
L}\ \mathbb{E}_{x}(\varphi(X_{t})\mathbbm{1}_{t<\tau_{\partial}})\leq
e^{Ct}\,\sup_{x\in L}\varphi(x)<\infty,$

which implies the third line of Assumption (F2).

In addition, because of the local uniform ellipticity of the diffusion $X$,
for all $n_{0}\geq 1$, $\psi_{2}:=\sum_{k=0}^{n_{0}}P_{k}\mathbbm{1}_{L}$ is
uniformly bounded away from zero on all compact subsets of $D$. This and
Theorem 3.5 concludes the proof of Theorem 4.1.

###  12.4 Proof of Corollary 4.2

Using Theorem 3.5, there exists $\lambda^{\prime}_{0}$ such that, for all
$x\in D$,

$\eta(x)=\lim_{t\rightarrow+\infty}e^{\lambda^{\prime}_{0}t}\mathbb{P}_{x}(t<\tau_{\partial}).$

We choose in the definition of $\lambda_{0}$ a ball $B$ such that
$\nu_{QSD}(B)>0$ (recall that $\lambda_{0}$ is independent of the choice of
$B$). Given $x\in D$ such that $\eta(x)>0$,

$\lim_{t\rightarrow+\infty}e^{\lambda^{\prime}_{0}t}\mathbb{P}_{x}(X_{t}\in
B)=\eta(x)\nu_{QSD}(B)\in(0,+\infty).$

Hence, $\lambda_{0}=\lambda^{\prime}_{0}$ and the infimum in the definition of
$\lambda_{0}$ is a minimum. The rest of the properties stated in Corollary 4.2
are direct consequences of Theorem 3.5.

Let us now prove that $\eta$ is $\mathcal{C}^{2}$. First, it follows from [79,
Theorem 7.2.4] that $e^{\lambda_{0}t}\mathbb{P}_{x}(t<\tau_{\partial})$ is
continuous for all $t\geq 0$ (see [18] for a detailed proof). Hence the
uniform convergence in Theorem 2.5 implies that $\eta$ is continuous on $D$.

Now, let $B$ be any non-empty open ball such that $\overline{B}\subset D$. We
consider the following initial-boundary value problem (in the terminology of
[39]) associated to the differential operator $\mathcal{L}$ defined in (4.3)

$\begin{cases}\partial_{t}u(t,x)-\mathcal{L}u(t,x)-\lambda_{0}u(t,x)=0&\text{for
all }(t,x)\in(0,T]\times B,\\\ u(0,x)=\eta(x)&\text{for all }x\in B,\\\
u(t,x)=\eta(x)&\text{for all }(t,x)\in(0,T]\times\partial B.\end{cases}$

Since the coefficients of $\mathcal{L}$ are Hölder and uniformly elliptic in
$\overline{B}$ and since $\eta$ is continuous, we can apply Corollary 1 of
Chapter 3 of [39] to obtain the existence and uniqueness of a solution $u$ to
the above problem, continuous on $[0,T]\times\overline{B}$ and
$\mathcal{C}^{1,2}((0,T]\times B)$. Now, we can apply Itô’s formula to
$e^{\lambda_{0}s}u(T-s,X_{s})$: for all $s\leq\tau_{B^{c}}\wedge T$ and all
$x\in B$, $\mathbb{P}_{x}$-almost surely,

$\displaystyle e^{\lambda_{0}s}u(T-s,X_{s})$
$\displaystyle=u(T,x)+\int_{0}^{s}e^{\lambda_{0}r}\left(-\frac{\partial
u}{\partial t}+\mathcal{L}u+\lambda_{0}u\right)(T-r,X_{r})\,dr$
$\displaystyle+\int_{0}^{s}e^{\lambda_{0}r}\nabla
u(T-r,X_{r})\sigma(X_{r})\,dB_{r}.$

Since $u$ is bounded and continuous on $[0,T]\times\overline{B}$ and $\nabla
u(t,x)$ is locally bounded in $(0,T]\times B$, it follows from standard
localization arguments that

$\displaystyle u(T,x)$
$\displaystyle=\mathbb{E}_{x}\left[e^{\lambda_{0}(T\wedge\tau_{B^{c}})}u(T-(T\wedge\tau_{B^{c}}),X_{T\wedge\tau_{B^{c}}})\right]$
$\displaystyle=\mathbb{E}_{x}\left[e^{\lambda_{0}(T\wedge\tau_{B^{c}})}\eta(X_{T\wedge\tau_{B^{c}}})\right].$

Now, the Markov property and the fact that $P_{t}\eta=e^{-\lambda_{0}t}\eta$
entail that $e^{\lambda_{0}t}\eta(X_{t})$ is a martingale on
$(\mathcal{D},(\mathcal{F}_{t})_{t\geq 0},\mathbb{P}_{x})$, hence

$\eta(x)=\mathbb{E}_{x}\left[e^{\lambda_{0}(T\wedge\tau_{B^{c}})}\eta(X_{T\wedge\tau_{B^{c}}})\right]=u(T,x).$

Therefore, $\eta\in\mathcal{C}^{2}(D)$ and
$\mathcal{L}\eta(x)=-\lambda_{0}\eta(x)$ for all $x\in D$.

## References

  * [1] E. Arjas, E. Nummelin, and R. L. Tweedie.  Semi-Markov processes on a general state space: $\alpha$-theory and quasistationarity.  J. Austral. Math. Soc. Ser. A, 30(2):187–200, 1980/81.
  * [2] K. B. Athreya and P. E. Ney.  Branching processes.  Springer-Verlag, New York, 1972.  Die Grundlehren der mathematischen Wissenschaften, Band 196.
  * [3] R. Azaï s, J.-B. Bardet, A. Génadot, N. Krell, and P.-A. Zitt.  Piecewise deterministic Markov process—recent results.  In Journées MAS 2012, volume 44 of ESAIM Proc., pages 276–290. EDP Sci., Les Ulis, 2014.
  * [4] M. Baudel and N. Berglund.  Spectral theory for random poincaré maps.  SIAM J. Appl. Math., 49(6):4319–4375, 2017.
  * [5] N. Berglund and D. Landon.  Mixed-mode oscillations and interspike interval statistics in the stochastic FitzHugh-Nagumo model.  Nonlinearity, 25(8):2303–2335, 2012.
  * [6] G. Birkhoff.  Extensions of Jentzsch’s theorem.  Trans. Amer. Math. Soc., 85:219–227, 1957.
  * [7] F. M. Buchmann.  Simulation of stopped diffusions.  J. Comput. Phys., 202(2):446–462, 2005.
  * [8] D. L. Burkholder.  Martingale transforms.  Ann. Math. Statist., 37:1494–1504, 1966.
  * [9] P. Cattiaux, P. Collet, A. Lambert, S. Martínez, S. Méléard, and J. San Martín.  Quasi-stationary distributions and diffusion models in population dynamics.  Ann. Probab., 37(5):1926–1969, 2009.
  * [10] P. Cattiaux and S. Méléard.  Competitive or weak cooperative stochastic Lotka-Volterra systems conditioned to non-extinction.  J. Math. Biol., 60(6):797–829, 2010.
  * [11] J. A. Cavender.  Quasi-stationary distributions of birth-and-death processes.  Adv. Appl. Probab., 10(3):570–586, 1978.
  * [12] N. Champagnat, K. Coulibaly-Pasquier, and D. Villemonais.  Exponential convergence to quasi-stationary distribution for multi-dimensional diffusion processes.  ArXiv e-prints, Mar. 2016.
  * [13] N. Champagnat, P. Diaconis, and L. Miclo.  On Dirichlet eigenvectors for neutral two-dimensional Markov chains.  Electron. J. Probab., 17:no. 63, 41, 2012.
  * [14] N. Champagnat and S. Rœlly.  Limit theorems for conditioned multitype Dawson-Watanabe processes and Feller diffusions.  Electron. J. Probab., 13:no. 25, 777–810, 2008.
  * [15] N. Champagnat and D. Villemonais.  Exponential convergence to quasi-stationary distribution and Q-process.  Probab. Theory Related Fields, 164(1):243–283, 2016.
  * [16] N. Champagnat and D. Villemonais.  Population processes with unbounded extinction rate conditioned to non-extinction.  ArXiv e-prints, Nov. 2016.
  * [17] N. Champagnat and D. Villemonais.  Exponential convergence to quasi-stationary distribution for absorbed one-dimensional diffusions with killing.  ALEA Lat. Am. J. Probab. Math. Stat., 14(1):177–199, 2017.
  * [18] N. Champagnat and D. Villemonais.  Lyapunov criteria for uniform convergence of conditional distributions of absorbed Markov processes.  ArXiv e-prints, Apr. 2017.
  * [19] N. Champagnat and D. Villemonais.  Uniform convergence of conditional distributions for absorbed one-dimensional diffusions.  J. Appl. Probab., June 2017.  To appear.
  * [20] N. Champagnat and D. Villemonais.  Uniform convergence of penalized time-inhomogeneous Markov processes.  ESAIM Probab. Stat., Mar. 2018.  To appear.
  * [21] S. D. Chatterji.  An $L^{p}$-convergence theorem.  Ann. Math. Statist., 40:1068–1070, 1969.
  * [22] J.-R. Chazottes, P. Collet, and S. Méléard.  Sharp asymptotics for the quasi-stationary distribution of birth-and-death processes.  Probab. Theory Related Fields, 164(1-2):285–332, 2016.
  * [23] J.-R. Chazottes, P. Collet, and S. Méléard.  On time scales and quasi-stationary distributions for multitype birth-and-death processes.  ArXiv e-prints, Feb. 2017.
  * [24] P. Collet, S. Martínez, and J. San Martín.  Asymptotic laws for one-dimensional diffusions conditioned to nonabsorption.  Ann. Probab., 23(3):1300–1314, 1995.
  * [25] P. Collet, S. Martínez, and J. San Martín.  Quasi-stationary distributions.  Probability and its Applications (New York). Springer, Heidelberg, 2013\.  Markov chains, diffusions and dynamical systems.
  * [26] C. Coron.  Stochastic modeling and eco-evolution of a diploid population.  PhD thesis, Palaiseau, Ecole polytechnique, 2013.
  * [27] C. Coron, S. Méléard, E. Porcher, and A. Robert.  Quantifying the mutational meltdown in diploid populations.  The American Naturalist, 181(5):623–636, 2013.
  * [28] J. N. Darroch and E. Seneta.  On quasi-stationary distributions in absorbing discrete-time finite Markov chains.  J. Appl. Probab., 2:88–100, 1965.
  * [29] J. N. Darroch and E. Seneta.  On quasi-stationary distributions in absorbing continuous-time finite Markov chains.  J. Appl. Probab., 4:192–196, 1967.
  * [30] E. B. Davies and B. Simon.  Ultracontractivity and the heat kernel for Schrödinger operators and Dirichlet Laplacians.  J. Funct. Anal., 59(2):335–395, 1984.
  * [31] P. Del Moral.  Feynman-Kac formulae.  Probability and its Applications (New York). Springer-Verlag, New York, 2004.  Genealogical and interacting particle systems with applications.
  * [32] P. Del Moral.  Mean field simulation for Monte Carlo integration, volume 126 of Monographs on Statistics and Applied Probability.  CRC Press, Boca Raton, FL, 2013.
  * [33] P. Del Moral and D. Villemonais.  Exponential mixing properties for time inhomogeneous diffusion processes with killing.  Bernoulli, 24(2):1010–1032, 2018.
  * [34] M. Faure and S. J. Schreiber.  Quasi-stationary distributions for randomly perturbed dynamical systems.  Ann. Appl. Probab., 24(2):553–598, 2014.
  * [35] P. A. Ferrari, H. Kesten, and S. Martínez.  $R$-positivity, quasi-stationary distributions and ratio limit theorems for a class of probabilistic automata.  Ann. Appl. Probab., 6(2):577–616, 1996.
  * [36] P. A. Ferrari, H. Kesten, S. Martínez, and P. Picco.  Existence of quasi-stationary distributions. A renewal dynamical approach.  Ann. Probab., 23(2):501–521, 1995.
  * [37] P. A. Ferrari, S. Martínez, and P. Picco.  Some properties of quasi-stationary distributions in the birth and death chains: a dynamical approach.  In Instabilities and nonequilibrium structures, III (Valparaíso, 1989), volume 64 of Math. Appl., pages 177–187. Kluwer Acad. Publ., Dordrecht, 1991.
  * [38] P. A. Ferrari, S. Martínez, and P. Picco.  Existence of nontrivial quasi-stationary distributions in the birth-death chain.  Adv. Appl. Probab., 24(4):795–813, 1992.
  * [39] A. Friedman.  Partial differential equations of parabolic type.  Prentice-Hall, Inc., Englewood Cliffs, N.J., 1964.
  * [40] E. Gobet.  Weak approximation of killed diffusion using Euler schemes.  Stochastic Process. Appl., 87(2):167–197, 2000.
  * [41] E. Gobet.  Euler schemes and half-space approximation for the simulation of diffusion in a domain.  ESAIM Probab. Statist., 5:261–297, 2001.
  * [42] G. L. Gong, M. P. Qian, and Z. X. Zhao.  Killed diffusions and their conditioning.  Probab. Theory Related Fields, 80(1):151–167, 1988.
  * [43] P. Good.  The limiting behavior of transient birth and death processes conditioned on survival.  J. Austral. Math. Soc., 8:716–722, 1968.
  * [44] F. Gosselin.  Asymptotic behavior of absorbing Markov chains conditional on nonabsorption for applications in conservation biology.  Ann. Appl. Probab., 11(1):261–284, 2001.
  * [45] M. Hairer.  Convergence of markov processes (lecture notes).  www.hairer.org/notes/Convergence.pdf, 2010.
  * [46] T. E. Harris.  The theory of branching processes.  Die Grundlehren der Mathematischen Wissenschaften, Bd. 119. Springer-Verlag, Berlin; Prentice-Hall, Inc., Englewood Cliffs, N.J., 1963.
  * [47] C. R. Heathcote, E. Seneta, and D. Vere-Jones.  A refinement of two theorems in the theory of branching processes.  Teor. Verojatnost. i Primenen., 12:341–346, 1967.
  * [48] A. Hening and M. Kolb.  Quasistationary distributions for one-dimensional diffusions with singular boundary points.  ArXiv e-prints, Sept. 2014.
  * [49] G. Hinrichs, M. Kolb, and V. Wachtel.  Persistence of one-dimensional AR(1)-sequences.  ArXiv e-prints, Jan. 2018.
  * [50] K. Itô and H. P. McKean, Jr.  Diffusion processes and their sample paths, volume 125 of Die Grundlehren der mathematischen Wissenschaften.  Springer-Verlag, Berlin, 1974.  Second printing, corrected.
  * [51] A. Joffe and F. Spitzer.  On multitype branching processes with $\rho\leq 1$.  J. Math. Anal. Appl., 19:409–430, 1967.
  * [52] O. Kallenberg.  Foundations of modern probability.  Probability and its Applications (New York). Springer-Verlag, New York, second edition, 2002.
  * [53] I. Karatzas and S. E. Shreve.  Brownian motion and stochastic calculus, volume 113 of Graduate Texts in Mathematics.  Springer-Verlag, New York, second edition, 1991.
  * [54] S. Karlin and J. McGregor.  The classification of birth and death processes.  Trans. Amer. Math. Soc., 86:366–400, 1957.
  * [55] M. Kijima and E. Seneta.  Some results for quasi-stationary distributions of birth-death processes.  J. Appl. Probab., 28(3):503–511, 1991.
  * [56] J. F. C. Kingman.  The exponential decay of Markov transition probabilities.  Proc. London Math. Soc. (3), 13:337–358, 1963.
  * [57] R. Knobloch and L. Partzsch.  Uniform conditional ergodicity and intrinsic ultracontractivity.  Potential Anal., 33:107–136, 2010.
  * [58] M. Kolb and D. Steinsaltz.  Quasilimiting behavior for one-dimensional diffusions with killing.  Ann. Probab., 40(1):162–212, 2012.
  * [59] N. V. Krylov and M. V. Safonov.  A property of the solutions of parabolic equations with measurable coefficients.  Izv. Akad. Nauk SSSR Ser. Mat., 44(1):161–175, 239, 1980.
  * [60] A. Lambert.  Quasi-stationary distributions and the continuous-state branching process conditioned to be never extinct.  Electron. J. Probab., 12:no. 14, 420–446, 2007.
  * [61] J. Littin C.  Uniqueness of quasistationary distributions and discrete spectra when $\infty$ is an entrance boundary and 0 is singular.  J. Appl. Probab., 49(3):719–730, 2012.
  * [62] M. Lladser and J. San Martín.  Domain of attraction of the quasi-stationary distributions for the Ornstein-Uhlenbeck process.  J. Appl. Probab., 37(2):511–520, 2000.
  * [63] P. Maillard.  The $\lambda$-invariant measures of subcritical Bienaymé-Galton-Watson processes.  Bernoulli, 24(1):297–315, 2018.
  * [64] P. Mandl.  Spectral theory of semi-groups connected with diffusion processes and its application.  Czechoslovak Math. J., 11 (86):558–569, 1961.
  * [65] R. Mannella.  Absorbing boundaries and optimal stopping in a stochastic differential equation.  Phys. Lett. A, 254(5):257–262, 1999.
  * [66] S. Martínez and J. San Martín.  Classification of killed one-dimensional diffusions.  Ann. Probab., 32(1A):530–552, 2004.
  * [67] S. Martínez, J. San Martín, and D. Villemonais.  Existence and uniqueness of a quasistationary distribution for Markov processes with fast return from infinity.  J. Appl. Probab., 51(3):756–768, 2014.
  * [68] S. Méléard and D. Villemonais.  Quasi-stationary distributions and population processes.  Probab. Surv., 9:340–410, 2012.
  * [69] S. Meyn and R. L. Tweedie.  Markov chains and stochastic stability.  Cambridge University Press, Cambridge, second edition, 2009.  With a prologue by Peter W. Glynn.
  * [70] S. P. Meyn and R. L. Tweedie.  Stability of Markovian processes. III. Foster-Lyapunov criteria for continuous-time processes.  Adv. in Appl. Probab., 25(3):518–548, 1993.
  * [71] Y. Miura.  Ultracontractivity for Markov semigroups and quasi-stationary distributions.  Stoch. Anal. Appl., 32(4):591–601, 2014.
  * [72] Y. Ogura.  Asymptotic behavior of multitype Galton-Watson processes.  J. Math. Kyoto Univ., 15(2):251–302, 1975.
  * [73] R. G. Pinsky.  On the convergence of diffusion processes conditioned to remain in a bounded region for large time to limiting positive recurrent diffusion processes.  Ann. Probab., 13(2):363–378, 1985.
  * [74] R. G. Pinsky.  Explicit and almost explicit spectral calculations for diffusion operators.  J. Funct. Anal., 256(10):3279–3312, 2009.
  * [75] P. E. Protter.  Stochastic integration and differential equations, volume 21 of Applications of Mathematics (New York).  Springer-Verlag, Berlin, second edition, 2004.  Stochastic Modelling and Applied Probability.
  * [76] B. Roynette, P. Vallois, and M. Yor.  Some penalisations of the Wiener measure.  Jpn. J. Math., 1(1):263–290, 2006.
  * [77] E. Seneta and D. Vere-Jones.  On quasi-stationary distributions in discrete-time Markov chains with a denumerable infinity of states.  J. Appl. Probab., 3:403–434, 1966.
  * [78] D. Steinsaltz and S. N. Evans.  Markov mortality models: Implications of quasistationarity and varying initial conditions.  Theo. Pop. Bio., 65(65):319–337, 2004.
  * [79] D. W. Stroock and S. R. S. Varadhan.  Multidimensional diffusion processes.  Classics in Mathematics. Springer-Verlag, Berlin, 2006.  Reprint of the 1997 edition.
  * [80] E. A. van Doorn.  Quasi-stationary distributions and convergence to quasi-stationarity of birth-death processes.  Adv. Appl. Probab., 23(4):683–700, 1991.
  * [81] E. A. van Doorn.  Conditions for the existence of quasi-stationary distributions for birth-death processes with killing.  Stochastic Process. Appl., 122(6):2400–2410, 2012.
  * [82] E. A. van Doorn and P. K. Pollett.  Quasi-stationary distributions for discrete-state models.  European J. Oper. Res., 230(1):1–14, 2013.
  * [83] A. Velleret.  Non-uniform exponential convergence to qsd, with a possibly stabilizing extinction.  In preparation, 2018.
  * [84] D. Vere-Jones.  Ergodic properties of nonnegative matrices. I.  Pacific J. Math., 22:361–386, 1967.
  * [85] D. Villemonais.  Minimal quasi-stationary distribution approximation for a birth and death process.  Electron. J. Probab., 20:no. 30, 18, 2015.
  * [86] J. Wang.  First eigenvalue of one-dimensional diffusion processes.  Electron. Commun. Probab., 14:232–244, 2009.
  * [87] J. Wang.  Sharp bounds for the first eigenvalue of symmetric Markov processes and their applications.  Acta Math. Sin. (Engl. Ser.), 28(10):1995–2010, 2012.
  * [88] A. M. Yaglom.  Certain limit theorems of the theory of branching random processes.  Doklady Akad. Nauk SSSR (N.S.), 56:795–798, 1947.