White Dwarfs with Infrared Excess from LAMOST Data Release 5

Infrared excess is an important probe of sub-stellar companions and/or debris disks around white dwarfs (WDs). Such systems are still rare for in-depth understanding of their formation and long-term evolution. One of the largest spectroscopic surveys carried out by the Large sky Area Multi-Object fiber Spectroscopic Telescope (LAMOST) recently released more than $3000$ WDs, a significant fraction of which have not undergone excess search. Here we present cross-correlation of LAMOST DR5 WD catalog with the Pan-STARRS, SDSS, UKIDSS, 2MASS, and {\it WISE}. By performing SED (spectral energy distribution) fitting for 846 WDs with $WISE$ detections, we identify 50 candidates with infrared excess, including 7 candidate WD+M dwarf binaries, 31 candidate WD+brown dwarf (BD) binaries and 12 candidate WD+dust disk systems. 8 of the dust disk systems are our new identifications. Utilizing a systematic survey with accurate stellar parameters derived from spectral fitting, our work is an important addition to previous searches for infrared excess from SDSS and {\it Gaia} WDs, and provides a significant ($\gtrsim8\%$) complement to current database of WDs with candidate BD companions and dust disks. The frequencies of WD+BD binaries and WD+dust disk systems are constrained to be $\lesssim3.7\%$ and $\sim1.4\%$, respectively. The properties of candidate dust disk systems are discussed. All of our candidates require follow-up observations for confirmation owing to limited spatial resolution of {\it WISE}.


INTRODUCTION
White dwarfs (WDs) are the final evolutionary stage and remnants of most stars ( 8 − 10M ; Iben et al. 1997). Infrared excess around a WD may come from an M dwarf companion, a brown dwarf (BD) companion, or a debris disk surrounding it. A statistical sample of WD+BD binaries and/or WD+debris disk systems are crucial for testing long-term evolutionary models of planetary systems around stars of 8 − 10M .
BDs are substellar objects of ∼ 13−80 Jupiter masses, too light to sustain nuclear fusion of hydrogen but sufficient for deuterium burning (e.g., Chabrier & Baraffe 2000;Burrows et al. 2011). Observations of BD companions to main-sequence stars have revealed a paucity of separation within ∼ 3 au, namely "BD desert" (e.g., Marcy & Butler 2000;Grether & Lineweaver 2006; Debris disks are believed to form through tidal disruption of asteroids (or comets) perturbed from locations comparable to the Kuiper belt in our solar system (Bonsor et al. 2011) into the Roche lobe of the WD (∼ 1R ; Davidsson 1999) by a remnant planetary system (e.g., Jura 2003;Debes et al. 2012;Veras 2014). Theoretical studies predict that a fraction of planets can survive the red giant phase of the host star (e.g., Villaver & Livio 2007, 2009Nordhaus et al. 2010;Mustill & Villaver 2012), which is supported by discoveries of debris disks (e.g., Gänsicke et al. 2006;Farihi et al. 2009; Barber et al. 2016), disintegrating planetary bodies (e.g., Vanderburg et al. 2015;Xu et al. 2016), and recently planets (Vanderburg et al. 2020;Blackman et al. 2021) orbiting around WDs. Dust grains in the disk will spiral in towards the WD via Poynting-Robert drag (Rafikov 2011), and the sublimated dust will accrete onto the WD's surface due to viscous torques (Rafikov & Garmilla 2012;Veras et al. 2013). This results in a spectroscopically detectable pollution of the otherwise pristine WD atmosphere, since the gravitational settling time is short compared to the WD's cooling age and the metals will rapidly diffuse out of the WD's photosphere (e.g., Koester 2009). Observations suggest that at least 20 − 30 percent of WDs are polluted by metals (e.g., Zuckerman et al. 2003Zuckerman et al. , 2010Koester et al. 2014).
A typical observational feature of debris disks around WDs is an infrared excess over the WD photosphere. The first infrared excess around a WD was detected by Zuckerman & Becklin (1987) and was later interpreted as emission from circumstellar dust by Graham et al. (1990). Afterwards, lots of efforts have been put into searches of infrared excess around WDs (e.g., Debes et al. 2011;Xu et al. 2020;Lai et al. 2021). Currently, about 100 WDs have been identified to host debris disks (e.g., Hoard et al. 2013;Rocchetto et al. 2015;Guo et al. 2015a;Farihi 2016;Rogers et al. 2020;Lai et al. 2021), largely owing to the sensitive infrared observations provided by Spitzer Space Telescope (Werner et al. 2004). However, Spitzer mostly focused its observations on WDs with polluted atmosphere (e.g., Farihi et al. 2009;Xu & Jura 2012) or in a specific range of stellar effective temperature (e.g., Rocchetto et al. 2015;Wilson et al. 2019). The frequency of debris disks around WDs is constrained to be 1 − 4 percent, and for heavily polluted WDs, the occurrent rate of dusty disks is found to exceed 50 percent (e.g., Mullally et al. 2007;Farihi et al. 2009;Barber et al. 2012;Rocchetto et al. 2015;Wilson et al. 2019).
Infrared surveys such as Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) and UKIRT Infrared Deep Sky Survey (UKIDSS; Lawrence et al. 2007) provide opportunities for comprehensive searches of infrared excess around WDs. The first large and deep untargeted search for infrared excess were carried out by Steele et al. (2011) and Girven et al. (2011), via cross-correlating Sloan Digital Sky Survey (SDSS; York et al. 2000) WDs with UKIDSS catalog. The release of SDSS DR7 WD catalog (Kleinman et al. 2013) further triggered excess search via cross-correlation with WISE all-sky survey (Debes et al. 2011;Barber et al. 2014), although source confusion is unavoidable due to the poor spatial resolution (∼ 6 ) of WISE (e.g., Dennihy et al. 2020). More recently, Gaia released its second WD catalog (Gentile Fusillo et al. 2019), providing the largest photometric WD sample for infrared excess search via colors and/or magnitudes (Rebassa-Mansergas et al. 2019;Xu et al. 2020;Lai et al. 2021). These studies identified hundreds of candidates with infrared excess in total, and candidate WD+BD binaries and WD+dust disk systems are fewer ( 100), since only part of the studies classified them further via detailed spectral energy distribution (SED) modeling.
The Large sky Area Multi-Object fiber Spectroscopic Telescope (LAMOST) is a powerful spectroscopic survey operated by National Astronomical Observatories, Chinese Academy of Sciences (Cui et al. 2012). In the past ten years, LAMOST has released thousands of spectroscopically identified WDs with stellar parameters (e.g., effective temperature T eff , surface gravity log g) derived from spectral fitting (Zhao et al. 2013;Zhang et al. 2013;Guo et al. 2015bGuo et al. , 2022, which should in principle be more accurate and reliable than that obtained from SED fitting if the signal-to-noise ratio (S/N) of the spectrum is above some criterion, e.g., S/N> 30 (Guo et al. 2022). The vast majority of LAMOST WDs locate at low Galactic latitudes Deng et al. 2012), providing a complement to SDSS DR7 WD catalog (Kleinman et al. 2013). Recently, LAMOST released its 5th spectroscopically identified WD catalog (Guo et al. 2022). More than 1000 of them have not been searched for infrared excess, offering opportunities for discoveries of new targets of interest, e.g., WDs with a debris disk or a BD companion. Without cut-off in the spectral type or effective temperature of the WDs, this catalog provides an unbiased sample in this point of view, although other selection effects may exist (see discussions in Section 5.3). In this paper, we present systematic search for infrared excess in the LAMOST DR5 WD catalog according to WISE photometry. The poor spatial resolution of WISE make background contamination unavoidable, which is a shortcoming of this work.
The paper is organized as follows. In Section 2, we present photometric data of our sample by crossmatching with various surveys in optical and infrared. Then we describe SED models of different components and search for infrared excess in Section 3. We further check the optical and infrard images of our candidates in Section 4 to exclude possible background contamination or M dwarf companions that are not the focus of this work. Results and discussion are presented in Section 5 and conclusions are given in Section 6.
2. PHOTOMETRIC DATA LAMOST DR5 WD catalog contains 3064 unique sources (see Guo et al. 2022, Table 3). To search for infrared excess around those WDs via SED fitting in Section 3, prior information on stellar parameters (e.g., T eff and log g) is required to simultaneously model the radiation of the WD photosphere as well as an additional component (if any) with the lack of data at wavelengths longer than several mirons. Moreover, Gaia distance is adopted to convert the absolute magnitude to radiative flux for various photometric systems when applying WD cooling model (Bergeron et al. 1995). For these two reasons, we only select targets with T eff , log g, and Gaia distance measurements. Some of the targets have multiple observations, and each spectrum has best-fit values of T eff and log g. We adopt the fitting parameters derived from the highest S/N spectrum. This results in a sample of 1179 WDs (hereafter the parent sample), and all of them were classified as either DA or DB types 1 (Guo et al. 2022).
Then we search for optical and infrared data by crosscorrelating the LAMOST coordinates of the parent sample with various surveys within a radius of 2 (e.g., Debes et al. 2011). Our excess search mainly relies on WISE photometry, with W 1, W 2, W 3, and W 4 bands centered at 3.4, 4.6, 12, and 22 µm, respectively. Crosscorrelating with AllWISE catalog results in 459, 303, 15 and 7 targets, respectively for W 1 − W 4 detections. CatWISE2020 catalog (Marocco et al. 2021) provides an excellent alternative for W 1 and W 2 bands as the detection limit is fainter and the sample size is roughly twice larger than AllWISE. There are 828 and 685 of the parent sample having CatWISE2020 detections in W 1 and W 2 bands, respectively. For W 1 and W 2 photometry, we give priority to CatWISE2020 data, and if is not available, we adopt AllWISE data (the case for 18 targets). Since W 3 and W 4 data are not provided by CatWISE2020, we directly use AllWISE data if available.
We then cross-match the parent sample with the Two-Micron All Sky Survey (2MASS; JHK s ; Skrutskie et al. 2006). Considering the limited resolution of 2MASS, we further cross-match these targets with UKIDSS catalog from Galactic Plane, Galaxy Cluster, and Large Area Survey (YJHK; Lawrence et al. 2007). 293 and 292 of the parent sample are detected in the 2MASS K s and UKIDSS K bands, respectively, and 77 of the targets are detected by both filters. We notice that some of the WDs have inconsistent photometry for 2MASS and UKIDSS data. Since they have similar wavelength coverage and UKIDSS has better spatial resolution, we give priority to UKIDSS data and adopt 2MASS data only if the corresponding band UKIDSS data are not available.
In order to characterize the photospheric radiation of WDs, the parent sample are further cross-matched with Pan-STARRS (g p r p i p z p y p ; Kaiser et al. 2010) catalog. As a consistency check, we also search for their SDSS DR12 detections (ugriz; Alam et al. 2015). This returns 1110 and 895 detections for g p and g bands of the two surveys, respectively, and most of these targets have the other four optical-bands observations. Except for the above-mentioned cases of either Cat-WISE2020 or AllWISE, and either 2MASS or UKIDSS, all of these optical and infrared data are used in the SED fitting if the target is detected in that filter, i.e., the apparent magnitude is given with errors instead of upper limits. We add additional 2% in quadrature to the reported uncertainties in the radiative flux in each photometric band to account for systematic errors.

MODELS AND SED FITTING
In order to search for infrared excess around the WDs and distinguish among possible scenarios for their origin, we adopt T eff and log g measurements of LAMOST WD catalog derived from spectral fitting of Balmer absorption lines (Guo et al. 2022) in the SED fitting, and our models involve (i) photospheric emission of a WD; (ii) the SED of a late-type (M/brown dwarfs) companion; and (iii) radiation spectrum of a dust disk, as detailed below.
(i) Photospheric emission of WDs. For a grid of T eff and log g, WD cooling models 2 (Bergeron et al. 1995) provide an estimate of the mass and cooling age for the WD as well as its synthetic photometry in all optical and infrared bands we studied. For a WD in the LAMOST catalog, we force log g = 7 if it is below 7, and force log g = 9 if it exceeds 9, since 7 and 9 are the lower and upper boundaries of the cooling models. 24 of the parent sample have surface gravity out of this range.
(ii) SED templates of M/brown dwarfs. For late-type dwarf companions, empirical relations have been constructed between the spectral type and SDSS colors as well as absolute J-band magnitude (Hawley et al. 2002). Alternatively, We adopt the colors provided by Best et al. (2018) involving Pan-STARRS, 2MASS, and WISE as a function of the spectral type, and interpolate to get colors of SDSS and UKIDSS filters. Since the absolute J-band magnitude is available only for dwarfs cooler than M6 type, we use the absolute magnitudes provided by Hawley et al. (2002) for all of the M-, L-and T-types dwarfs for consistency.
(iii) Radiation spectra of dust disks. We adopt Jura (2003) model to describe the dust disk emission, which assumes a geometrically flat, optically thick disk illuminated by the central WD and reradiating at infrared. The temperature of each annulus in the disk is determined by the effective temperature (T eff ) of and its distance (R ring ) to the WD, i.e., where R WD is the WD radius. Radiation from each annulus can be approximated as a blackbody spectrum, and its integration over the radius from the inner to outer edges of the disk gives the total disk flux: where i is the inclination of the disk, D is the WD distance to the observer, k B is the Boltzmann's constant, h is the Planck constant, c is the speed of light, and x ≡ hν/(k B T ring ).
For the purpose of SED fitting, the magnitude in each photometric band is converted into the flux density according to the zero points provided by SVO Filter Profile Service (Rodrigo et al. 2012;Rodrigo & Solano 2020). According to the measured T eff and log g for each WD, we interpolate to obtain the absolute magnitude in various photometric bands according to WD cooling models (Bergeron et al. 1995). Then model flux is converted from the absolute magnitude by adopting Gaia distances of the WDs, and we add a free normalization to account for measurement errors of log g and distances. The cooling models are initially fit to the SDSS and Pan-STARRS fluxes of the WD by minimizing where f i,mod is the model flux for the ith photometric band in optical, and f i,obs and σ i,obs are the observed flux and its 1σ error, respectively. Then an infrared excess with respect to the WD cooling model is searched with a significance expressed by the photometric and model uncertainties, i.e., where the model uncertainty is also considered and σ i,mod = 0.05f i,mod is assumed (e.g., Rocchetto et al. 2015;Xu et al. 2020). The criterion we adopt for an excess is χ w > 3 for WISE W 1 band or longward, i.e., the excess is at a > 3σ significance level. Besides that, it should also be satisfied that the source image in this WISE band is unaffected by any known artifacts, i.e., the quality flag (cc flag) of the band is zero.
Once an excess is detected, we add an additional model component, i.e., a late-type dwarf companion or a dust disk. For the case of a dwarf companion, we adopt the WD distance provided by Gaia, and for each spectral type from M1 to T6, χ 2 is calculated, and the lowest χ 2 value corresponds to the best-fit spectral type of the companion.
For the dust disk case, degeneracies exist between the inclination angle and inner disk boundary, especially for the case here that only one or two data points are available for excess modeling. Here we set the inclination and inner boundary of the disk as free parameters. The dust disk can extend to ∼ 100 WD radii. Constraints on the outer disk boundary require photometric data at longer wavelengths, which is not the case here. We therefore take R out = 80 R WD and assume an initial face-on inclination, i.e., i = 0 • (e.g., Debes et al. 2011). The inner disk radii is capped by 50 R WD , and the lower boundary is set by the sublimation temperature (T sub ) of dust grains, as indicated by Equation (1). For silicatedominated dust, T sub ∼ 1200 K, and for pure carbon dust, T sub ∼ 2100 K (Debes et al. 2011). We adopt a more relaxed value of T sub = 3000 K, as it was pointed out that the sublimation temperature of the dust could be much higher due to the high metal vapor pressure in the inner disk (Rafikov & Garmilla 2012).
The companion model is compared with the dust disk model via the reduced-χ 2 (χ 2 /d.o.f ) of SED fitting, and the model with lower χ 2 value is preferred.

IMAGES
The infrared excess found in Section 3 are mainly based on WISE observations, and due to its poor spatial resolution, a significant fraction of the excess could arise from background contamination and/or late-type stellar companions (possibly M dwarfs). LAMOST DR5 WD catalog already identified some WD+M dwarf binaries, while a fraction of M dwarf companions may be too faint to exhibit features in WD spectra and could be missed in the classification. The contaminants radiating in WISE bands are expected to radiate at a little bit shorter wavelengths since the spectrum is continuous. If it is not as cool as BDs, it should be detected in the K band, and possibly z band.
UKIDSS image provides an effective tool to exclude background contamination owing to its superior spatial resolution. However, only a fraction of the WDs with infrared excess have UKIDSS observations. We notice that Pan-STARRS z band data are also good in quality and may provide an alternative for contamination check, although the filters are even shorter in wavelength. As Figure 1 shows, the apparent two point-sources in the UKIDSS K-band image are blended in the WISE W 1 image, and Pan-STARRS z-band can clearly resolve the two point sources. The circle with a radius of 6 can marginally enclose the central source in the WISE image, which is also true for our other candidates with infrared excess. We checked the UKIDSS and Pan-STARRS images of all those WDs with infrared excess selected in Section 3, and those with multiple sources within a radius of 6 are excluded. While WD+BD binaries are not expected to show two point-sources in the K or z band, some of the WDs with an M dwarf companion and having images similar to Figure 1 may be excluded in this process.
For the leftover candidates, we further compare their W ISE images and optical ones using the WISEView tool (http://byw.tools/wiseview), which is useful to look for possible contamination in the field. We remove those candidates with obvious contamination in the W ISE W 1 image.

RESULTS AND DISCUSSION
The SED fitting and image check result in a total of 50 candidates with infrared excess, including 7 candidate WD+M dwarf binaries, 31 candidate WD+BD binaries, and 12 candidate WD+dust disk systems. 37 out of the 50 candidates with infrared excess are new identifications, and 8 candidate dust disk systems are first reported as such systems.
Our results rely upon the accuracy of T eff measurement. An incorrect T eff will result in an inaccurate measurement of WD photospheric radiation, possibly leading to a false or magnified excess in infrared. Moreover, our results may still suffer from background contamination. While UKIDSS and/or Pan-STARRS image check is able to remove contamination mainly radiating at shorter wavelength, those redder contaminants (e.g., red galaxies) may still be missed due to the limited resolution of W ISE. These factors could lead to (very) high χ 2 values in the SED fitting and should be treated with caution. Below we will list all selected candidates regardless of its χ 2 value, and each of them waits for confirmation or refutation via follow-up observations. Table 1 lists the first part of WDs with candidate BD companions. Also listed are the 7 WDs with an M dwarf companion. WD+M dwarf binaries are not the focus of this paper, and our parent sample do not contain WDs spectroscopically classified as WD+M dwarf binaries by LAMOST, and the image check in Section 4 may have excluded some of such systems, leaving only a few here. However, considering that the uncertainties of the spectral type of the companion could be several dex, we also list WD+M dwarf candidates here for reference.

WD+BD binary candidates
The second part of the WD+BD candidates are listed in Table 2. These candidates have comparable reducedχ 2 resulting from the companion and the dust disk models. As examples, Figure 2 shows the SED fitting results for three WD+BD binaries. The best-fit spectral type of the BD companion is L9−9.9, L0−0.9 and L5−5.9, respectively for the left, middle, and right panels.
Two of previously known dust disk systems (WD 2328+107, ID: 470; and WD 0843+516, ID: 505) were misclassified as WD+BD binaries by our SED model fitting, mainly because T eff given in the LAMOST catalog is higher than the best-fit value. An overestimated T eff will shift blueward the model SED of the WD, creating spurious excess at shorter wavelengths, which mimics radiation from a late-type dwarf companion. When the best-fit T eff is adopted, the resulting χ 2 for the companion model and dust disk model are both reduced, and the dust disk model is now preferred. This suggests that the identification of infrared excess and the best-fit J131939.94+111131.3

WISE W1
UKIDSS K Pan-STARRS z Figure 1. A case of confused WISE W 1 (3.4µm; left panel) photometry classified as infrared excess according to SED model fitting but is removed from our candidate list due to contamination visible in the images. UKIDSS K-(1.25µm; middle panel), and Pan-STARRS z-band (right panel) images are presented as comparison. The field of view is 30 × 30 for each cutout, and the blue circle is centered at the LAMOST coordinate of the WD with a radius of 6 . As shown, Pan-STARRS can clearly distinguish the two point sources, though the contaminant is fainter compared to its UKIDSS image.  Figure 2. Examples of SED fitting results for WD+BD candidates. Colored symbols are photometric data from SDSS (magenta diamonds), Pan-STARRS (green triangles), 2MASS (blue squares), UKIDSS (cyan circles), and WISE (red stars). All data have been shown in the figure, and missing symbols in a panel means that there is no photometry in the corresponding survey (same for Figure 4). The flux in W 3 and W 4 bands given with upper limits are not used in the model fitting. The black solid line is the best-fitting model, which is the sum of WD photospheric radiation (blue dashed line) and BD contribution (red dotted line). The best-fit spectral types of the companion are L9−9.9, L0−0.9 and L5−5.9 BDs, respectively for the left, middle, and right panels. The y-axis errors of the data points are too small to be visible.
model are sensitive to the accuracy of T eff measurement. We have moved these two targets to the WD+dust disk classification, and for this reason we also consider the candidates in Table 2 in the following statistical studies on dust disk properties.

WD+dust disk candidates
Our WD+dust disk candidates are listed in Table 3. All these WDs have spectral types of either DA or DB as expected, since this is also the case for our parent sample. Although most of the dust disk candidates found by Debes et al. (2011) have similar stellar spectral types, WDs with known dust disks usually show metal absorption lines in WD spectra and are mostly classified as DAZ or DBZ. However, metal lines could be missed in the WD spectrum with a relatively low S/N. Three of the known WDs with a debris disk (WD 1018+410; ID: 282; WD 0843+516, ID: 505; and Ton 345, ID: 923) classified as DAZ or DBZ in the literature (see Table 1 of Xu et al. 2020) are classified as DA or DB in the LAMOST catalog (Guo et al. 2022). The effective temperatures of these WDs mostly locate at 25000 K, consistent with Note-Columns from left to right represent: (1)-(3) the unique ID (same as the Group ID in Guo et al. 2022) and LAMOST coordinate; (4) G band magnitude given by Gaia DR2; (5)-(9) spectral type, effective temperature, surface gravity, mass and cooling age of the WD given by LAMOST catalog; (10) Gaia distance of the WD; Columns (11)-(12): best-fit spectral type of the companion, and minimum reduced χ 2 . a The WD was reported by Debes et al. (2011) as candidate WD+M dwarf binaries.
b The WD was identified as exhibiting WISE excess by Xu et al. (2020).
the expectation that dust grains may not survive around a too hot WD. Indeed, those published WD+dust disk systems all have similar WD temperatures (see Xu et al. 2020, and references therein).
The cutout images of these candidate WD+dust disk systems are shown in Figure 3. The UKIDSS K-band image is our first choice, and if it is not available, the Pan-STARRS z-band image is shown. The blue circle centered at the LAMOST coordinate of the WD is to illustrate that there is no contamination within a radius of 6 .
The SED fitting results of the dust disk candidates are shown in Figure 4. The disk inclinations of the candidates concentrate at i > 60 • and the inner disk radii cover a broad range from several to tens of WD radii. As the disk is illuminated by the central WD, the disk temperature declines from the inner to outer radii. The blue cutoff of the disk radiation is mainly set by the inner disk boundary and the declining rate at red end is determined by the outer disk boundary. The inclination can affect the normalization, peak, and shape of the disk SED and is difficult to uniquely constrain due to degeneracies with the inner disk radius, especially for those targets with excess only in W 1 band. Therefore, the best-fit parameters of those candidate dust disk systems should be treated as indicative if the dust disk is confirmed. Figure 5 shows mass vs. cooling age distribution of the WDs for the dust disk candidates. Our candidates occupy similar parameter space to the WDs already known to host a dust disk. Including the targets with statistically similar fits between dust disk and companion models in Table 2 (red open squares) does not change the overall distribution.
In Figure 6, we plot distribution of the inner disk radius (R in ) vs. effective temperature of the WDs for  Table 3. The blue circle is centered at the LAMOST coordinate of the WD with a radius of 6 .  Note-Columns (1)-(10): similar to those in Table 1; (11)-(12): best-fit spectral type of the companion, and minimum reduced-χ 2 for the companion model; (13)-(15): best-fit inner radius and inclination of the disk, and minimum reduced-χ 2 for the dust disk model. a The WD was identified to have WISE excess by Xu et al. (2020).
b The WD was confirmed to have Spitzer excess by Lai et al. (2021). Note-Columns (1)-(10): similar to those in Table 1; (11)-(13) best-fit inner radius and inclination of the dust disk, and minimum reduced-χ 2 . a The WD is a previously known dust disk system.
b The WD was reported as a candidate WD+M dwarf binary by Debes et al. (2011).
c The WD was identified with WISE excess by Xu et al. (2020). our candidates. According to Equation (1), we have R in ∝ (T eff /T sub ) 4/3 , where the sublimation temperature T sub of dust grains defines a power-law relation between R in and T eff (black curves). A significant fraction of our candidates locate at T sub = 3000 K, higher than the previously known dust disks and candidate dust disks identified by Debes et al. (2011). This is simply caused by our settings of T sub ≤ 3000 K, which is to account for possible effects that may elevate sublimation temperature of dust grains (Rafikov & Garmilla 2012). Beside that, the distribution of our candidates are in general consistent with the literature.
The mass of the dust disk cannot be constrained if it is optically thick (e.g., Jura 2003). However, if the disk is optically thin, under the assumption of a modified blackbody radiation, the dust disk mass can be roughly estimated via the dust temperature T d and the radiative flux F ν of the disk at a specific wavelength λ, similar to the commonly used method to calculate the dust mass of protoplanetary disks and interstellar dust mass (e.g., Schuller et al. 2009;Utomo et al. 2019): where D is the WD distance, κ(λ) is the dust absorption cross-section per unit mass at λ, and B ν (T d , λ) is the Planck function for temperature T d at λ. We take λ = 2.2 µm and adopt κ(2.2 µm) = 3800 cm 2 g −1 (Shirley et al. 2011). We have assumed multi-color blackbody radiation in the SED models described in Section 3. Here the dust temperature is assumed to be constant across the disk, and can be roughly estimated via the peak radiation of the dust SED according to Wien's displacement law. F ν in Equation (5) is obtained from the model flux radiated by the dust disk, i.e., the red dotted lines in Figure 4. Figure 7 shows the estimated dust disk mass vs. WD cooling age for our candidate WD+dust disk systems. For most of the candidate systems, the dust disk mass locates at ∼ 10 15 − 10 18 g, much lower than the total mass of Saturn's rings (∼ 10 22 g; Esposito 1993). As comparison, the dust disk of G29-38, a well studied dust disk system, was estimated to have a mass of 2×10 17 −10 22 g (the blue diamond in Figure 7), which is sensitive to parameters such as dust grain size and mostly relies on the optically-thin disk assumption (Jura 2003;Reach et al. 2005;Farihi et al. 2008Farihi et al. , 2009. For our candidates, the optically-thin disk assumption here will

Mass (M )
Known disks Debes+11 Our disk candidates Our potential candidates Figure 5. Mass vs. cooling age of the WDs. Red filled squares are our dust disk candidates (in Table 3), and the red open squares are the potential dust disk candidates (in Table 2) with similar χ 2 for the dust disk and the companion models. WDs with known dust disk are shown as the gray triangles, and dust disk candidates identified by Debes et al. (2011) are marked as cyan circles. . Inner disk radius vs. the effective temperature of the WD for (candidate) dust disk systems. The legend is similar to Figure 5. The dotted, dashed, and solid lines are the relations given by dust sublimation temperature of 1200 K, 2100 K and 3000 K, respectively.
also underestimate the dust mass, and the red squares in Figure 7 should be treated as lower limits. Pearson correlation between the two parameters results in a co-10 7 10 8 10 9 t cool (yr)  Table 3 (Table 2). Cyan circles are the dust disk candidates of Debes et al. (2011). All the red squares and cyan circles represent lower limits for the dust mass.
The blue diamond show the lower limit of dust mass given by the literature (see the text). efficient of 0.57 (for the open and closed red squares), implying a positive correlation between the dust mass and WD cooling age. This might tentatively suggest disruption of planetesimals and settling of dust on the disk plane along the central WD's cooling life if such a relation is confirmed in the future. Also shown is the dust disk candidates identified by Debes et al. (2011), where the dust mass is similarly calculated according to Equation (5). The dust mass for their candidates are in general 1 − 2 orders higher than ours. The main reason is that the disk inclination is fixed at i = 0 • for most of their candidates.

The occurrent rate
We found 31 candidate WD+BD binaries from an input of ∼ 3000 WDs in the LAMOST catalog. However, the requirement of W 1 band observations (in Cat-WISE or AllWISE), T eff measurements and Gaia distances left only 846 targets for excess search via SED model fitting. The event rate of WD+BD binaries is therefore 31/846 ≈ 3.7%, about twice higher than value (0.5 − 2 percent) reported in the literature (e.g., Steele et al. 2011;Girven et al. 2011). However, our frequency should be treated as an upper limit. Two of the previously known dust disk systems were originally identified as WD+BD binaries by us, with one listed in Table 1 and the other in Table 2. We expect that a significant fraction of the candidate WD+BD binaries, especially those listed in Table 2, could be dust disk systems if the excess is real. Follow-up observations such as radial velocity measurement of WDs (e.g., Maxted et al. 2006;Steele et al. 2013) or high resolution infrared imaging (e.g., Wilson et al. 2019) may confirm or refute such systems.
The twelve WD+dust disk candidates we identified convert to a frequency of 12/846 ≈ 1.4% for such systems, consistent with 1 − 4 percent reported in the literature (e.g., Rocchetto et al. 2015). If we include the 19 potential dust disk candidates in Table 2 with comparable χ 2 for dust disk and companion models, the frequency of dust disks will be 3.7%. We reaffirm that for some of them, the infrared excess could still be spurious if T eff measurement is inaccurate or could be caused by background contamination. Moreover, the final 846 sample left for SED model fitting may suffer from selection biases, since only WDs of DA or DB types with relatively high S/N spectra have T eff measurements in the LAMOST catalog.

Comparison with the literature
The 50 candidates we found with infrared excess include 4 WDs with previously known dust disks, 5 targets already reported by Xu et al. (2020) as exhibiting un-WISE excess, and 7 sources studied by Xu et al. (2020) yet no excess found by them. Our different results could arise from several factors, e.g., (i) the effective temperature of WDs is obtained by profile-fitting of Balmer absorption lines in LAMOST WD catalog, which could be different from that derived from Gaia's auto fit; (ii) we set the excess criterion as exceeding the WD photospheric radiation at > 3σ level, while they adopt a more significant > 5σ level; (iii) the CatWISE and/or AllWISE photometric data we used are not exactly the same as the unWISE data they adopted. 2 out of the 5 targets commonly identified as exhibiting WISE excess by us were confirmed to have Spitzer excess (Lai et al. 2021).
Besides one WD with known dust disk, 4 of our candidates with infrared excess were also identified by Debes et al. (2011). We have 2 common candidate WD+M dwarf binaries. 2 candidates classified as WD+M dwarf binaries by them are WD+BD and WD+dust disk systems respectively in our classification. One of these two targets were also studied by Xu et al. (2020), who identified it as with no excess, neither in magnitudes nor in colors. An insignificant infrared excess may lead to different classifications by Debes et al. (2011) and us.

CONCLUSIONS
We search for infrared excess around 3064 WDs recently released by LAMOST via SED model fitting of 846 targets with W ISE detections, effective temperature and distance measurement. The WD's effective temperature provided by LAMOST catalog is adopted in the SED fitting, and thus the identification of candidates with infrared excess relies on the accuracy of effective temperature measurement. Possible contamination have been preliminary excluded by visual check of optical and infrared images of these candidates. We identified 50 candidates with infrared excess, including 7 candidate WD+M dwarf binaries, 31 candidate WD+BD binaries, and 12 candidate WD+dust disk systems. The frequencies of WD+BD binaries and WD+dust disk systems are estimated to be 3.7% and ∼ 1.4%, respectively. 8 of these WD+dust disk candidates are first reported as such systems. Our dust disk candidates occupy similar parameter space to previously known dust disk systems in the WD mass vs. cooling age plot. The SED model fitting suggests that the temperatures of the inner dust disks locate at 1200 − 3000 K, generally consistent with that reported in the literature. The dust masses are constrained to have a lower limit of ∼ 10 15 − 10 18 g under the assumption of optically-thin dust disks. Positive correlations are found between the dust mass and WD cooling age. We caution that all these candidates with infrared excess require follow-up infrared imaging or infrared spectroscopy for confirmation.