The bright end of the exo-Zodi luminosity function: Disk evolution and implications for exo-Earth detectability

We present the first characterisation of the 12um warm dust ("exo-Zodi") luminosity function around Sun-like stars, focussing on the dustiest systems that can be identified by WISE. We detect six new warm dust candidates, five of which have unknown ages. We show that the dustiest old (>Gyr) systems like BD+20 307 are 1 in 10,000 occurrences. Bright warm dust is more common around young (<120Myr) systems, with a ~1% occurrence rate. We show that a two component in situ model where all stars have initially massive warm disks and in which warm debris is also generated at some random time along the stars' main-sequence lifetime, perhaps due to a collision, can explain the observations. However, if all stars only have initially massive warm disks these would not be visible at Gyr ages, and random collisions on the main-sequence are too infrequent to explain the high disk occurrence rate for young stars. That is, neither component can explain the observations on their own. Despite these conclusions, we cannot rule out an alternative model in which comets are scattered from outer regions because the distribution of systems with the appropriate dynamics is unknown. Our in situ model predicts the fraction of stars with exo-Zodi bright enough to cause problems for future exo-Earth imaging attempts is at least ~10%, and is higher for populations of stars younger than a few Gyr. This prediction of ~10% applies to old stars because bright systems like BD+20 307 imply a population of fainter systems that were once bright but are now decaying through fainter levels. Our prediction should be strongly tested by the Large Binocular Telescope Interferometer, providing valuable input for more detailed evolution models. A detection fraction lower than our prediction could indicate that the hot dust in systems like BD+20 307 has a cometary origin due to the quirks of the planetary dynamics.


INTRODUCTION
Perhaps the most important long term goal of humankind is to find and communicate with aliens, the inhabitants of extra-Solar Systems. This goal is a difficult one, with many challenges to overcome before it is technically possible to even image an Earth-like planet around another star. When it is eventually attempted, imaging detection of these planets could still be thwarted by the presence of exo-Zodiacal dust located in or near the terrestrial zone, which with sufficient dust levels can mask the presence of exo-⋆ Email: gkennedy@ast.cam.ac.uk Earths (e.g. Beichman et al. 2006;Absil et al. 2010;Roberge et al. 2012). Therefore, in preparation for future missions that will target specific stars, a critical piece of information is whether those stars harbour dust at a level that precludes the detection of an Earth-like planet.
While such knowledge will be crucial in specific cases, stars searched for warm dust may not be targeted in an unbiased way that yields information on the underlying distribution and its origins. For example, the most promising targets may be deemed to be those without cool Edgeworth-Kuiper belt analogues, or with other specific properties that may unbeknownst to us correlate with the presence of faint warm dust. The results cannot therefore be used to construct the brightness distribution (a.k.a. luminosity function, defined formally at the start of section 3) of warm dust belts, or to make predictions for stars that were not observed or about which no dust was detected. Therefore, an equally important goal is to characterise the luminosity function, with the overall aim of understanding warm dust origin and evolution.
Evolution is of key importance here, because the brightest warm dust sources seen around nearby stars cannot necessarily maintain the same brightness indefinitely. The observed emission comes from small particles that are eventually lost from the system, and the very existence of emission for longer than terrestrial zone orbital timescales-a few years-means that mass in small grains is being lost and replenished, and the small particles necessarily originate in a finite reservoir of larger objects. The question is whether those parent bodies are located at the same distance from the star as the dust (what we call an "in situ" scenario) or at a much more distant location (what we call a "comet delivery" scenario).
The possibility of such different origins means that the rate at which warm dust decays is unknown. For example, if the observed dust has its origins in an outer Edgeworth-Kuiper belt analogue and is scattered inward by planets (i.e. a comet delivery scenario, Nesvorný et al. 2010;Bonsor & Wyatt 2012;, the warm dust could be extremely long lived. Given the lack of evolution seen for main-sequence exo-Kuiper belts around Sunlike stars (e.g. Trilling et al. 2008), dust levels may not appear to decay at all during the stellar main-sequence lifetime. If, however, the dust originates near where it is observed (i.e. in situ), it would decay due to a decreasing number of parent bodies via collisions in a way that is reasonably well understood (e.g. Dominik & Decin 2003;Wyatt et al. 2007a;Löhne et al. 2008). In this scenario it also follows that the brightest warm disks are necessarily the tip of an iceberg of fainter disks; the decay timescale depends on the amount of mass that is present so systems decay as 1/time (e.g. Dominik & Decin 2003), and for every detection of bright warm dust that was created relatively recently, many more exist that are older and fainter, and are slowly grinding themselves to oblivion. Another conclusion from this kind of evolution is that the bright warm dust systems seen around relatively old stars cannot arise from in situ decay that started when the star was born; projecting the 1/time decay back to early times predicts implausibly large initial brightnesses, and delaying the onset of this decay for ∼Gyr timescales within a steady-state scenario requires implausibly large or strong planetesimals (Wyatt et al. 2007a).
These two different origins present very different structural and dynamical pictures of planetary systems, which may manifest differently in the warm dust luminosity function. For example, an observed luminosity function similar to that expected from 1/time decay would be a smoking gun for in situ origin and evolution. However, because the luminosity function in the comet delivery scenario is at least in part set by the (unknown) distribution of systems with the appropriate structural and dynamical properties, and the source regions may themselves be decaying as 1/time, the observation of such a function would not be absolute proof. As a further complication, it may be that the processes that generated the dust in BD+20 307-like systems (i.e. dusty warm systems whose host stars are ∼Gyr old) are completely different to those that result in the faintest systems. It could for example be that they are the result of late (i.e. ∼Gyr) dynamical instabilities , while faint Asteroid belt and Zodiacal cloud-like systems are the decay products of young massive belts that emerged from the protoplanetary disk phase or arise due to comet delivery. Therefore, our work here must be complemented by studies that for example aim to discover whether warm dust is correlated with the presence of cooler dust (e.g. Absil et al. 2013).
Here we use the stars observed by Hipparcos (van Leeuwen 2007) as an unbiased sample to characterise the bright end of the warm dust luminosity function using mid-IR photometry from the WISE mission (Wright et al. 2010). A secondary goal is to discover new bright warm disks, though few discoveries are expected among these stars due to previous IRAS and AKARI studies with a similar goal (e.g. Aumann & Probst 1991;Oudmaijer et al. 1992;Song et al. 2005;Fujiwara et al. 2010Fujiwara et al. , 2012bMelis et al. 2012).
The key difference between those studies and ours is that we count the non-detections, and therefore quantify the rarity of bright warm dust systems. We first outline our sample ( §2) and excess ( §3) selection, and then cull disk candidates from these excesses ( §4). We go on to compare the luminosity function with a simple in situ evolutionary model and make some basic predictions within the context of this model for future surveys for faint warm dust ( §5). Finally, we discuss some of the issues with our model and comet delivery as an alternative (but equally plausible) exo-Zodi origin ( §6).

SAMPLE DEFINITION
Stars are selected from the Hipparcos catalogue (van Leeuwen 2007) using an observational HR diagram. Because warm dust and planets have been discovered over a wide range of spectral types that might be taken as "Sun-like", we use a fairly broad definition of 0.3 < BT − VT < 1.5, in addition only keeping stars with BT < 11 to ensure the Tycho-2 photometric uncertainties are not too large (Høg et al. 2000), particularly for faint cool M types and white dwarves that otherwise creep into the sample. Using the Tycho-2 Spectral Type Catalogue (Wright et al. 2003), the mainsequence spectral types whose most common BT − VT colours are 0.3 and 1.5 are A8 and K5 respectively. While stars have a range of colours at a given spectral type, the bulk of the spectral types in our sample will lie between these bounds. To ensure that we do not select heavily reddened or evolved stars that are wrongly positioned in the HR diagram, we only retain stars within 150pc, and restrict the parallax S/N of stars to be better than 5. Finally, we exclude stars within 2 • of the Galactic plane, as we found that many WISE source extractions, even for relatively bright stars, are unreliable there. With these cuts 27,333 stars remain.
We exclude evolved stars, which may show excesses due to circumstellar material associated with mass loss, by excluding stars with MH p > 5 when BT − VT > 1.1, and with MH p > 25/3(BT − VT − 0.5) when BT − VT 1.1 (3,159 stars). 1 This cut is based purely on by-eye inspection of Figure 1, but as long as it lies in the Hertzsprung gap the sample size is only weakly sensitive to actual position because there are few stars there. The line is broken at 1.1 because a single line does not simultaneously exclude most stars crossing the Hertzsprung gap and include premain-sequence stars such as HD 98800 (the rightmost filled circle at BT −VT = 1.4, Hp = 6). Because they are selected solely using optical photometry and parallax, the 24,174 stars selected by these criteria are unbiased with respect to the presence of warm dust. The Hipparcos HR diagram for the sample is shown in Figure 1, which shows all 27,333 stars within the BT − VT range, the cut used to  exclude evolved stars (dashed line), and those found to have 12µm excesses by the criteria described below (large dots and circles).

Age distribution
It will become clear later that knowing the age distribution of our sample stars is useful. The simplest assumption would be that the ages are uniformly distributed. However, given the tendency for stars to have earlier types than the Sun (Fig 1), this assumption yields an unrealistically high fraction of ∼10Gyr old main-sequence stars (i.e. spectral types later than the Sun). Our sample spans a sufficiently wide range of spectral types that the main-sequence lifetimes vary significantly. For example, the earliest spectral types are late A-types, with ∼Gyr main-sequence lifetimes, while the latest spectral types are late K-types, with mainsequence lifetimes exceeding the age of the Universe. Stars are therefore assumed to have ages distributed uniformly between zero and their spectral type-dependent main-sequence lifetime, with the overall age distribution derived using the distribution of BT − VT colours (i.e. the colours are used as a proxy for the distribution of spectral types, a full description is given in Appendix B). The final age distribution is shown in Figure 2, which we assume throughout. Using a more simplistic assumption of uniformly distributed ages between 0-10Gyr for all stars does not alter any of our conclusions. . Spectral type-weighted age distribution adopted for our sample. The line shows the cumulative fraction (divided by ten). Thus, about 50% of our stars are younger than 2.5Gyr. The dashed line shows the constant relative frequency of 0.1035 and the fourth order polynomial fit used to generate synthetic age populations (see Appendix B).

EXCESS SELECTION
If warm dust is sufficiently bright relative to the host star, it will manifest as a detectable mid-infrared (i.e. 10-20µm) excess relative to the stellar photosphere, which for ∼300K dust peaks near 12µm. The WISE bands are well suited for this task, and because it has about an order of magnitude better sensitivity to 300K dust compared to W4 (see Fig 2 of Kennedy & Wyatt 2012) we use the 12µm W 3 band. 2 Our luminosity function is therefore formally the distribution of 12µm disk-to-star flux density ratios, which we will use in cumulative form. We use the term "luminosity function" largely for conciseness. The disk luminosity of course not only depends on the 12µm flux ratio, but also the disk temperature and stellar properties. If the ratio of the observed to stellar photospheric flux density is R12 = F12/F⋆, and the disk to star flux ratio ξ12 = R12 − 1, then the luminosity function is the fraction of stars with ξ12 above some level.
To derive the 12µm luminosity function, we would ideally have a sample for which the sensitivity to these excesses is the same for all stars. Then the distribution is simply the cumulative excess counts divided by the total sample size. Our use of Hipparcos stars with the above parallax requirement ensures that this ideal is met; all stars in our overall sample have better than 6% photometry at 1σ in W 3, and 98% have better than 2.5% photometry at 1σ. Given that the calibration uncertainty in this band is 4.5% (Jarrett et al. 2011), the fractional photometric sensitivity for all sample stars is essentially the same at about 5% 1σ. Therefore, the disk sensitivity is "calibration limited" for our purposes (see Wyatt 2008), with a 5% contribution from the WISE photometry, and a ∼2% contribution from the stellar photosphere. That is, we are limited to finding 3σ 12µm excesses brighter than about 15% of the photospheric level (flux ratios of R12 1.15). This property means that  Figure 1). The inset shows a smaller range of W 1 − W 3 with a linear count scale. Dark grey bars show all sources, and light grey bars show only sources where the WISE extension flag is not set (i.e. ext flg=0). Curves show Gaussian fits to the overall sample. The dashed line shows the threshold of W 1 − W 3 = 0.1mag for stars to become excess candidates (filled dots in Figure 1). we can select excesses by choosing stars with a red 3.4 − 12µm (W 1 − W 3) colour.
Though our excess limit is a uniform 15% due to the sample stars being bright, this brightness brings other potential problems. Stars brighter than 8.1, 6.7, 3.8, and -0.4mag are saturated in W 1 − 4 respectively, so nearly all stars are saturated in the W1 band to some degree, while only about a hundred are saturated in W3. However, comparison with 2MASS photometry shows that the WISE W1 source extraction is reliable well beyond saturation and shows no biases up to about 4.5mag, and only about 400 stars are brighter than this level. 3 Therefore, saturation actually affects a very small fraction of our stars. However, bright stars are the nearest to Earth and therefore potentially important for future studies if warm dust can be discovered, so we need to be sure that we are not biased against true excesses, i.e. that excesses selected based on W 1 − W 3 are not systematically missed due to some problem associated with W 3. We did not find any signs of W 3-specific issues, but noticed that the S/N of W 1 photometry depends fairly strongly on stellar brightness, with the brightest (i.e. more saturated) stars having lower S/N (i.e. stars brighter than about 7.25mag in VT have S/N< 20). The increased W 1 − W 3 scatter for these stars means some excesses may be missed for those with W 1 that happens to be scattered brighter than the actual brightness. This issue affects ∼1000 stars, and we show in Section 4.4 that we can rule out missing excesses for the worst of these cases. Given that there is no evidence for a problem with using W 1 − W 3 to identify excesses, we now proceed and leave possible issues with the WISE photometry to be identified in the individual source checking stage below.
The W 1 − W 3 distribution of the 24,174 stars is shown in Figure 3. A Gaussian fit centered on W 1 − W 3 = −0.041 with dispersion σ = 0.024 is shown, the fit being consistent with the photometric accuracy described above. The slightly negative average colour arises because our stars are later spectral types than 3 http://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec6 3c.html the reference A0V spectrum used for Vega magnitudes. The logarithmic scale (main plot) shows that systematic effects spread the distribution further for a relatively small number of sources. We therefore looked for properties that appear to predict W 1 − W 3 colours significantly different to -0.041 as indicators of systematic effects. Stars with companion sources at separations of 5-16 ′′ in the Hipparcos catalogue are more likely to have W 1 − W 3 < 0.3, as are WISE-saturated sources. We found that the best indicator of negative W 1 − W 3 was the WISE extension flag (ext flg), which indicates sources not well described by the WISE point spread function, or those near extended sources in the 2MASS Extended Source Catalogue. Aside from a few very bright sources, there is little brightness dependence on W 1 − W 3 colour, so it seems that the bulk of sources with blue WISE colours are due to poor source extractions for confused sources. Five of the warm dust candidates found below have ext flg = 1, but four are already known to host warm dust. For this reason, and because we subsequently check them individually we cull candidates from the full list of 24,174 stars.
The WISE W 3 calibration-limited sensitivity of about 5%, combined with the 2 average W 1 − W 3 colour of -0.04 means that a sensible threshold to use for true 12µm excesses that can be found to be significant is about 0.1mag. While we could choose a slightly lower threshold to include more sources that may turn out to be robust after more detailed analysis, this threshold is also practical in that it limits the number of sources to check by hand to a reasonable number; as we show below the false-positive rate already is fairly high and we expect this rate to increase strongly as the threshold moves closer to -0.04 (e.g. decreasing the threshold from 0.15 to 0.1 roughly doubles the number of sources to check, but only results in 25% more warm dust candidates, and decreasing the threshold to 0.05 would yield 181 more sources to check). Stars later than about K0 have slightly redder W 1 − W 3 than average (tending to about zero colour for BT − VT = 1.5). This trend has the potential to add more false positives, but in practise does not because it only affects a few hundred stars. We therefore identify as promising 12µm excess candidates the 96 mainsequence stars with W 1 − W 3 above 0.1, which are marked as large filled dots in Figure 1. Stars this red may show a 3σ W 3 excess based on the photometric uncertainties of most stars. To study these 96 sources further we use χ 2 minimisation to fit PHOENIX stellar photosphere (SED) models (Brott & Hauschildt 2005) to 2MASS, Hipparcos (Perryman & ESA 1997), Tycho-2 (Høg et al. 2000), AKARI (Ishihara et al. 2010), and WISE W 1 − 2 photometry (Wright et al. 2010) to provide a more accurate prediction of the W 3 flux density to compare with the WISE measurements. The SED fitting method has previously been validated through work done for the Herschel DEBRIS survey (e.g. Kennedy et al. 2012a,b)

WARM DUST CANDIDATES
We now check the 96 stars with candidate excesses for plausibility by inspecting the SEDs and WISE images for each source. The issues that arise are all due to WISE source extraction, but for several different reasons. Some are affected by bright nearby sources and image artefacts, while others are moderately separated binary systems where the WISE source extraction has only measured a single source. In some cases the W 1 photometry is shown to be underestimated, generally due to saturation. The results of the SED fitting alone are not sufficient for recognising spurious excesses and image inspection is generally required. While we inspected all of the images for the few candidates identified here, there are commonly indications in the WISE catalogue-primarily apparent variabilitythat might be used in a larger study. Notes on individual sources are given in Appendix A. We find 25 systems with plausible and significant 12µm excesses, of which all but seven were previously known to host significant levels of warm dust. These seven should be considered as promising candidates and require further characterisation, but are considered here to be real. The 25 are listed in Table 1 and detailed SEDs are shown in Appendix C.
The inclusion or exclusion of targets from the list of 25 in the context of our goal here merits some discussion. Because we wish to find the warm dust luminosity function for main-sequence stars, our sample should include any object that plausibly looks like a debris disk, even if it appears to be an extreme specimen. Though the most extreme disks, such as BD+20 307 are very rare, this rarity is unsurprising, at least in an in situ scenario; brighter (i.e. more massive) disks decay more rapidly due to more frequent collisions, resulting in a decreased detection probability at the peak of their activity.

Protoplanetary disks
We do not include protoplanetary disks in our luminosity function because they represent a qualitatively different phase of evolution to the debris disk phase (for example, models of their decay are not simply power-laws, Clarke et al. 2001). Including protoplanetary disks would require that models of the luminosity function include a prescription for the poorly understood transition from this phase to the debris disk phase, which would be very uncertain and limit our interpretation.
In terms of photometry the main feature that distinguishes protoplanetary disks from debris disks is the breadth of the disk emission spectrum. Protoplanetary disks cover a wide range of radii and extend right down to near the stellar surface, and therefore have near, mid, and far-IR excesses. On the other hand, debris disks are usually well described by blackbodies, sometimes with the addition of spectral features. In some cases (e.g. η Corvi, HD 113766A, Sheret et al. 2004;Wyatt et al. 2005;Morales et al. 2011;Olofsson et al. 2013), debris disks are poorly fit by a single temperature component, suggesting that populations of both warm and cool dust exist, and may argue for a comet delivery scenario for the warm dust with the cool component acting as the comet reservoir (Bonsor & Wyatt 2012). These systems are however easily distinguished from protoplanetary disks because the cool debris disk components always have much lower fractional luminosities than protoplanetary disks (L disk /L⋆ 10 −4 compared to ∼10 −1 , e.g. Trilling et al. 2008).
A more useful discriminant for our purposes is that debris disks never have large near-IR excesses. This lack of near-IR emission is first seen as the protoplanetary disk is dispersed in so called "transition disks" (Skrutskie et al. 1990). We therefore formalise the protoplanetary-debris disk distinction by considering the presence of a near-IR (3.4µm) excess, using the Ks −W 1 colour. While all sources were selected to have a red W 1 − W 3 colour, the addition of a sufficiently large W 1 excess indicates that the source is probably a protoplanetary disk. The Ks − W 1 colours for the sources in Table 1 are all less than 0.2, with three notable exceptions. The first two are the Herbig Ae stars HD 100453 and HD 145718 (Ks − W 1 ≈ 0.9), and the third is a new potential warm dust source, HD 166191 (Ks − W 1 = 0.5), which we discuss further below. The SEDs of all three are shown in Appendix C, and with both near and far-IR excesses are clearly different from the other 22. The detection of significant far-IR excesses by IRAS indicates emission from cold material but we did not use this as a discriminant because it would not have been detectable around all sources, even if a protoplanetary disk was present.
HD 166191 has only recently been noted as a potential warm dust source (Fujiwara et al. 2013), though was reported as having an IR excess from both IRAS and MSX (Oudmaijer et al. 1992;Clarke et al. 2005). 4 In addition, an excess was detected with AKARI-IRC, but not AKARI-FIS (though the upper limit from AKARI-FIS is about 10Jy, so not strongly constraining). An excess was found in the IRAS 60µm band, suggesting that the excess emission extends into the far-infrared. A large excess over a wide range of infrared wavelengths suggests that the emission is at a wide range of temperatures, and therefore that HD 166191 in fact harbours a protoplanetary disk. However, a nearby (1. ′ 7) red source was detected by WISE, MSX, and AKARI-IRC, opening the possibility that because IRAS had poor resolution the far-IR emission in fact comes from the nearby source and not HD 166191. A preliminary conclusion from new Hershel PACS (Pilbratt et al. 2010;Poglitsch et al. 2010) observations is that the disk spectrum is more consistent with a protoplanetary disk; given the WISE photometry, both the 70µm flux density (≈1.7Jy) and the fractional luminosity (≈10%) are larger than would be expected for a warm debris disk. A detailed study of this source will be presented elsewhere.
Finally, though it is not excluded by its Ks − W 1 colour, HD 98800 merits some discussion because whether it is a young debris disk or in the transition between the protoplanetary and debris disk phases is unclear. The disk emission spectrum is well modelled by a blackbody, as are most debris disks. The fractional luminosity is very high at around 10%, meaning that a significant fraction of the star's emission is being reprocessed by the disk. However, the far-IR emission may be decreased compared to that expected from a protoplanetary disk by truncation by the outer binary (Andrews et al. 2010). In addition, the detection of H2 emission suggests that HD 98800 may still harbour some circumstellar gas (Yang et al. 2012), though CO emission has not been detected (Zuckerman et al. 1995;Andrews et al. 2010). Based on the age, blackbody-like emission, and lack of CO detection and near-IR excess we include HD 98800 in our luminosity function, with the expectation that it will decay in a deterministic manner and that it is therefore both a plausible warm dust source and direct progenitor of fainter warm dust systems.

Warm debris disks
Having excluded three sources from the original list of 25, we now consider some properties of the remaining 22. For each source, after fitting a stellar photosphere model we fit a blackbody disk model to any significant infra-red excesses. In general blackbody models are a reasonable match to debris disk emission, though the agreement is poorest for warm dust systems, which commonly have noncontinuum silicate features (e.g. Song et al. 2005;Fujiwara et al. 2010). However, the temperatures derived are representative, and more complicated models are unwarranted because we typically only have 2-3 data points to fit. 5  The warm dust sources have a range of temperatures, which correspond to a range of radial distances. While we are not explicitly looking for habitable-zone dust, it is instructive to compare the temperatures to the width of the "habitable zone", which is naively the radial distance at which the equilibrium temperature is ∼280K (i.e. the temperature of a blackbody at 1AU from the Sun). The width of the habitable zone is of course very uncertain, and is probably at least a factor of two wide in radius (Kasting et al. 1993;Kopparapu et al. 2013). Therefore, the temperature range allowed is at least a factor √ 2 wide, giving a range from approximately 230-320K. In addition, there is considerable uncertainty in the location of dust belts inferred from SED models for several reasons; such as the likely presence of non-continuum spectral features (e.g. all warm dust targets considered transient by Wyatt et al. 2007a, have silicate features). In addition, belts may be a factor ∼2 more distant than inferred from blackbody models due to grain emission inefficiencies at long wavelengths (relative to their sizes, e.g. Booth et al. 2013). However, this factor, which has been derived from cool Kuiper-belt analogues, may not apply to fainter warm dust if larger grains dominate the emission. Finally, the picture of a ring of warm dust is probably over-simplified, particularly at faint levels where Poynting-Robertson drag becomes important. Based on this discussion, we retain all 12µm excess sources in what follows. In any case, removal of the coolest few sources (those below emission from small grains. While this modification is not needed, it ensures that sub-mm fluxes predicted by the models are more realistic than a pure blackbody.

HIP
200K say) would make little difference to our analysis. Temperature is of course worthy of future study, and for example the luminosity function could (with sufficient detections) be extended to include it as a third dimension.
It is clear from Table 1 that the 12µm excess systems are preferentially young, with at least 11 (and possibly 13) being members of the Scorpius-Centaurus association in Lower Centaurus Crux (LCC), Upper Centaurus Lupus (UCL), and Upper Scorpius (US) (i.e. <20Myr old). In addition, HD 98800 belongs to the 8Myr old TW Hydrae association (TWA), HD 165439 was derived to be 11Myr based on isochrone models (Tetzlaff et al. 2011), HD 15407A is a member of the AB Doradus moving group (Melis et al. 2010), and HD 22680 and HD 23157 belong to the 115Myr old Pleiades. Of the 7 stars that are not known to be young, only BD+20 307 is known to be older than about 1Gyr. HD 160959 was suggested to be a member of Upper Centaurus Lupus with 50% probability (i.e. 15Myr old Rizzuto et al. 2012), but also has an isochrone-derived age of 1.4Gyr (Holmberg et al. 2009). Based on the likelihood of stars with excesses to be young, the former age seems more likely to be correct, but merits further study. Similarly, HD 150697 is suggested to be 3.2Gyr old based on isochrones (Holmberg et al. 2009), but with a similar sky location and distance as the ρ Ophiuchi star-forming region (Fajardo-Acosta et al. 2000) requires further study to verify the age.
For a typical main-sequence lifetime of 10Gyr and stars spaced randomly throughout, there is a 1% chance of a star being randomly younger than 100Myr. However, 15 (possibly 18) of the 22 systems are younger than 120Myr, implying that the warm disk occurrence rate for stars this young is at least 5%. Therefore, as with longer wavelengths (i.e. 24µm, Siegler et al. 2007), there appears to be a significant decrease in the incidence of 12µm emission with time. However, this decay is unlikely to be universal; BD+20 307 is at least 1Gyr old, so assuming 1/time decay should have been ∼100 times brighter at 10Myr (i.e. verging on physically impossible given that the fraction of the host star's light that is captured by the dust in this system is currently 3.2%, Weinberger et al. 2011). Therefore, the more reasonable conclusion is that either the dust level in this system is not evolving, as might be expected in a comet delivery scenario, or became very bright some time well after the protoplanetary disk was dispersed, as might be expected for a recent collision .
The stars in Table 1 are biased towards earlier spectral types, though the histogram in Figure 1 shows that this could be explained by more early-type stars in the sample distribution. Based on the histogram there is some evidence that earlier spectral types are more likely to host 12µm excesses. The likely explanation is that later-type stars in the sample are biased somewhat against having younger ages, as it is only the brighter (earlier-type) stars in the nearby young associations that were selected for the Hipparcos sample.

Variability
Given the possibly transient nature of warm dust, particularly in old systems such as BD+20 307, we made a basic search for variability using the WISE data. This search is motivated by the likely orbital timescales being of order years, and the ability of individual collisions to affect the overall dust level (Kenyon & Bromley 2005;Wyatt et al. 2007a;Meng et al. 2012). Indeed, evidence for ∼year timescale variability has been seen in several relatively young warm dust systems (Melis et al. 2012;Meng et al. 2012).
Though stars were only observed over a day or so in a single WISE pass, another was made 6 months later for sources that were in the appropriate RA range. 6 Of the 22 warm dust sources, only 6 WISE did not survive an entire year, so the whole sky was not observed BD+20 307 and HD 103703 show evidence for WISE variability. HD 103703 is marked as variable in W 3 in but not W4. While it may be possible for W 3 to vary while W4 stays constant due to changes in the strength of silicate features or the emission from each being dominated by emission at different radial dust locations, the change of ∼0.1mag (±0.02mag) in only one band is suggestive but requires verification. In contrast, BD+20 307 shows a brightening of 6-8% in both W 3 and W 4 over a six month period (see also Meng et al. 2012), and as shown in Figure 4, these are correlated in the sense that an increase in W 3 brightness is accompanied by an increase at W4. This increase suggests that while the total mass (i.e. including parent bodies) must be decreasing (or at least be constant), effects such as individual collisions, size distribution variations and mineralogical changes can increase the amount of dust as measured by the infrared excess on ∼year long timescales.

Missing warm dust sources?
An apparent omission from our warm dust list is HD 69830. As one of the nearest few dozen Sun-like stars and host to both planets and warm dust, it is of key importance (Beichman et al. 2005). The reason it does not appear here is in fact simple; the excess at 12µm is not large enough to be formally detected with WISE. A detailed SED model that includes all available photometry, including that from the Spitzer Infra-Red Array Camera (Werner et al. 2004;Fazio et al. 2004), shows that the 12µm excess is only at the 5% level, which also agrees with IRAS (see also Figure 2 in Beichman et al. 2011). Another key warm dust source is the F2V star η Corvi, which has W 1 − W 3 of 0.14 and a W 3 excess significance of only 2σ. The "disappearing disk" TYC 8241-2652-1 (Melis et al. 2012) is not included here because it is not in the Hipparcos catalogue. Similarly, well known A-type warm dust sources such as β Pictoris and HD 172555 have significant WISE 12µm excesses, but with BT − VT colours of around 0.2 were not part of the initial sample.
We made a further check for potential warm dust sources within the Unbiased Nearby Stars (UNS) sample, which includes five samples, each having the nearest ∼125 main-sequence stars with spectral types of A, F, G, K, and M (i.e. 125 A-stars, 125 F-stars, etc., Phillips et al. 2010). UNS stars meeting the criteria outlined above are a subset of the larger sample in question here, and SEDs for this sample have been studied in great detail as part of the Herschel DEBRIS Key Programme (e.g. Wyatt et al. 2012;Kennedy et al. 2012a). With the inclusion of IRAS (Moshir et al. 1993) and AKARI (Ishihara et al. 2010) mid-IR photometry, and in many cases Spitzer Infra-Red Spectrograph (IRS, Houck et al. 2004;Lebouteiller et al. 2011) spectra, this subset therefore provides a good check for the sources that are subject to the worst saturation effects. We found no 12µm excesses in this subset that should have resulted in a significant detection with WISE. Figure 5 shows the 12µm luminosity function that results from the 22 unbiased warm dust detections described above. The distribution is generated simply by assuming that these excesses could have been detected around all 24,174 stars, and dividing the cumulative twice (Wright et al. 2010). Here we have not used the so called "3-Band Cryo" data release as any variability found would remain in question due to the change in cooling after the end of the full cryogenic mission. disk to star flux ratio distribution by this number. Thus, while the brightest warm dust systems such as BD+20 307 are known to be rare, we have quantified this rarity to be of order 1 in 10,000. Figure 5 also compares the luminosity function with the limit previously derived for the Kepler field (Kennedy & Wyatt 2012). Stars in the Kepler field are sufficiently faint that the number and distribution of detected excesses is consistent with that expected from chance alignments with background galaxies. The distribution derived here lies slightly above the limit from the Kepler field at the bright end, which most likely arises due to small number variation, though could also be suggesting that the brightest excesses in the Kepler field are in fact real but not very robust against background galaxy confusion.

The 12µm luminosity function
Given the clear bias towards young stars among our excess candidates, the dark lines with added dots in both panels of Figure 6 show our luminosity function split by age (the grey lines are models, described below). "Young" stars are those younger than 120Myr, and "old" stars are those older than 1Gyr (our warm dust sample has no stars between 120Myr and 1Gyr, see also §5.3). This split requires us to use the age distribution shown in Figure 2, for which the number of young stars is 607, and the number of old stars 19,114. The distributions in Figure 6 therefore divide the detections by these total number estimates. Stars with unknown ages are added to each group to illustrate the possible uncertainty in each distribution. The main point here is that old warm dust systems are extremely rare, while young ones are not.

IN SITU EVOLUTION
One powerful advantage of knowing the distribution of warm excesses is that it contains information about how disks evolve. For example, a population of disks whose brightness decays as 1/time that is observed at random times will be distributed with a slope of -1 in Figure 5. A disk spends ten times longer in each successive decade of luminosity than it did in the previous one, so every bright disk that is discovered is the tip of an iceberg of fainter disks that are decaying ever more slowly.
Our basic assumption in now making a model to compare to our derived luminosity function is that warm dust arises from some in situ process, for example a collision or collisions between parent bodies that normally reside where the dust is observed. In the alternative comet delivery scenario, where objects are injected from elsewhere (i.e. larger radial distances), no well informed population evolution model can currently be made because the behaviour depends on the detailed dynamics of each individual system, how those dynamics change over time, and how systems with the appropriate dynamics are distributed. Of course, this difficulty does not mean that comet delivery is a less viable scenario; we discuss it further in Section 6.2.

Model description
Our goal is therefore to use an in situ evolution model and apply it to the observed 12µm luminosity function, testing different scenarios and making predictions for surveys that probe fainter dust. The model described below has been outlined in previous works (Wyatt et al. 2007a,b;Wyatt 2008). The basic premise is that a size distribution of objects is created in catastrophic collisions between a reservoir of the largest objects. The resulting "collisional cascade" size distribution remains roughly constant in shape (Dohnanyi 1969), so the dust level (i.e. the number of smallest objects) is proportional to the number of the largest objects. To connect the mass with the observed dust level requires a relation between the total mass Mtot and the total surface area σtot, which depends on the maximum and minimum object sizes Dmin and Dc, and the slope of the size distribution q, where the number of objects between D and D + dD is n(D) = KD 2−3q (Wyatt et al. 2007b) Mtot = 2.5 × 10 −9 3q − 5 6 − 3q ρσtotDmin 10 9 Dc/Dmin where Mtot is in units of M⊕, ρ is the planetesimal density in kg m −3 , σtot is in AU 2 , Dmin is in µm, and Dc in km. With q = 11/6, the total mass is dominated by the largest objects and the total surface area by the smallest grains. To now connect the surface area with the observed dust level, we assume that the dust emits as a pure blackbody, so where Bν is the Planck function in Jy sr −1 and d is the distance to the star in pc. We approximate the stellar emission as where the stellar luminosity L⋆ is in units of L⊙ and the effective temperature T⋆ is in K. We define ξν as the disk to star flux ratio at some wavelength With equations (1) and (4) and some assumed planetesimal properties we are therefore able to link an excess observed at some wavelength and temperature with the total mass in the disk. Because the size distribution is such that most of the total mass Mtot is contained in the largest objects, the dust level decays at a rate proportional to the large-object collision rate, which has a characteristic timescale t coll (in Myr) (5) 10 -1 10 0 10 1 10 2 F disk /F star at 12µm where dr is the width of the belt at r (both in AU), Q ⋆ D is the object strength (J kg −1 ), and e is the mean planetesimal eccentricity, and M⋆ is the stellar mass. The disk radius can be expressed in terms of temperature with The mass loss rate is dMtot/dt ∝ M 2 tot and the mass decays as (Dominik & Decin 2003) where t is the time since the evolution started and Mtot,0 was the mass available for collisions at that time (with a consequent collision timescale t coll,0 ). This time may be the age of the host star tage, but may be smaller, for example if the belt was generated from a collision event that was the result of a dynamical instability well after the star and planets were formed. An equivalent equation also applies to the disk to star flux ratio Though numerical models find that the evolution can be slower than 1/time (e.g. Löhne et al. 2008;Gáspár et al. 2013), these results are based on disks at large semi-major axes where the largest objects are not in collisional equilibrium. For example, though Löhne et al. (2008) state that their disks decay as t −0.3 to −0.4 , they also note that both the disk and dust masses tend to 1/time decay at times that are sufficiently late that the largest objects have reached collisional equilibrium. Gáspár et al. (2013) find a much slower t −0.08 mass decay for their reference model. However, the slow decay appears to be because the largest (1000km) objects in their size distribution have not started to collide, even at the latest times, and hence the total mass (which is dominated by these objects) is only decaying due to the destruction of smaller objects (the lack of 1000km-object evolution for their reference model can be seen in Fig 1 of Gáspár et al. 2012). Because the warm disks we consider here are very close to their central stars (i.e. a few AU) and the collision rate scales very strongly with radius (eq 5), collisional equilibrium is reached in only a few Myr (which we also confirm a posteriori). Therefore, we consider that while our analytic model is necessarily simplified, a 1/time evolution is justifiably realistic.

Model implementation and interpretation
To implement the model we simply assume some initial disk brightness ξν and that the collision timescale has the same form as equation (5) t coll,0 = C/ξν,0 so C has units of time in Myr and is shorter for disks that are initially more massive, being the collision time for a disk of unity ξν,0 (and some equivalent Mtot,0 that can be calculated using equations 1 and 4). For our model C is the only variable parameter.
To interpret C in terms of physical parameters we combine equations (1), (4), and (5), and rewrite the answer in terms of C, giving C = 7.4 × 10 −12 r 13/3 ρ dr r 1 10 9 Dc Dmin Because the collision timescale depends on the disk mass (or equivalently brightness), C does not depend on the initial disk mass (see Wyatt et al. 2007a), but on a number of other parameters, most of which we can estimate values for. We assume a Sun-like host star with L⋆ = 1L⊙, T⋆ = 5800K, and M⋆ = 1M⊙. We assume that the disk lies at 1AU and hence T disk = 278.3K, and that dr = 1.
We assume a size distribution slope of q = 11/6. Finally, we include the wavelength of observation, λ = 12. With these assumptions equation (10) reduces to The remaining parameters are therefore the maximum and minimum planetesimal sizes, their strengths, and their eccentricities. We may further assume that the minimum grain size is Dmin = 1µm, the approximate size at which grains are blown out by radiation pressure from Solar-type stars. The coefficient for this equation does not change strongly with spectral type, and for example is 5 × 10 −7 for an early F-type star, assuming that the dust remains at the same temperature (i.e. has a larger radius of about 2.5AU) and that the minimum grain size has increased to 4µm due to the increased luminosity.

Excess evolution
In order to test the expectations for disk evolution we use the evolution model described above to construct two simple Monte-Carlo models of 12µm evolution. In the first "initially massive" scenario, all stars have relatively massive and bright warm disks at early times (see §6.1 for further discussion of how this scenario could be interpreted). In this picture the time since the onset of decay is simply the stellar age. All stars have disks regardless of age, but these become progressively fainter for populations of older stars. In the second "random collision" scenario, stars have a single dust creation event at some random time. Though we have a range of spectral types, we set the epoch of this event to be a random time between 0-13Gyr because it is unlikely to be related to stellar evolution (the time distribution of these events is however revisited in §6.1). Our age distribution means that many stars will be observed before such a collision happens (see §4.5), though this detail does not affect the model because it is degenerate in the sense that multiple (instead of single) collisions could be invoked and offset by a shorter collision time to yield the same results. Only those stars that happen to have a collision at a young age and are observed at a young age will have a detectable excess. We assume that stars have the age distribution shown in Figure 2. The initial distribution of flux ratios for both models is assumed to be log-normal with zero mean and unity dispersion in log(F disk /F⋆), so covers the bright end of the luminosity function. We found that the choice of initial distribution does not strongly influence the results, and for example a power-law distribution weighted towards fainter levels but with some bright disks yields very similar results.
We generate disks according to each scenario using two age bins; those that are "young" (<120Myr) and those that are "old" (>1Gyr), using the same numbers of stars noted in §4.5 (607 young stars and 19,114 old stars). We restrict the old stars to be older than 1Gyr because BD+20 307 is ∼1Gyr old, whereas placing the cut at >120Myr results in disk detections in the old group that are not much older than 120Myr, so not as old as the stars we wish to compare the model with. While there are other possible choices of bin locations, the point is that warm excesses around young stars are common, and those around old stars are not, which allows us to distinguish between the models outlined above. In each bin we generate 25 model realisations to illustrate the scatter due to small numbers. The only model parameter is C, which sets how rapidly disks decay, and is varied by hand so that the synthetic distributions have approximately the same level as those observed.
The results from the two models are shown in Figure 6. The left panel shows the initially massive warm belt scenario, while the right panel shows the random collision scenario. Looking first at the initially massive panel, while reasonable agreement is obtained for young stars, no old stars are seen to have large excesses as the dust has decayed significantly by Gyr ages (these disks have F disk /F⋆ ∼ 10 −3 , see 1-2Gyr dotted line for initially massive warm disks in Figure 9). This result echoes the conclusions of Wyatt et al. (2007a), who found that bright debris disks around old stars cannot arise from the in situ evolution considered in this scenario. The observed young-star distribution is usually flatter than most of the models, which might be explained by extra physical processes not included in our model. We note a few of these in section 6.1.
Turning now to the right panel of Figure 6, the random collision scenario reasonably reproduces the luminosity function of warm dust around old stars, but fails to match those seen around young stars. The chance that at least some of these are young, and the possibility that HD 150697 is in fact young, means that a somewhat larger value of C (slower evolution) could be used to find better agreement with BD+20 307 (i.e. shift the model lines upward, thereby making dust from random collisions more common). This slower evolution would however predict about ten times more systems with F disk /F⋆ 0.15 that are not observed. While a few young systems are seen to have excesses due to random collisions (about 1% at F disk /F⋆ ∼ 0.1), these are too infrequent to match the relatively high occurrence rate that is observed. The conclusion is perhaps obvious; if collisions occur randomly over time it is unlikely that a young star will have a collision that is observed soon afterwards. The rapid collisional evolution means that even if they are, the dust levels are unlikely to be extreme. The only significant model parameter is C, which is set to 1Myr for both scenarios and provides reasonable agreement with the observed luminosity function in each case. This parameter need not be the same, as the scenarios could for example have very different sized parent bodies or random velocities. As outlined above in §5.2, for various assumptions C can be interpreted in terms of physical system parameters. Assuming Dmin = 1, equation (11) implies that D 1/2 c (Q ⋆ D ) 5/6 e −5/3 = 3.2 × 10 6 . This value can be compared to previous results for this combination, for example Wyatt et al. (2007b) found 7.4×10 4 for the evolution of Kuiper belt analogues around A-stars, and Kains et al. (2011) found 1.4 × 10 6 for Kuiper belt analogues around Sun-like stars (i.e. both were for disks at much larger radii than those considered here). The latter authors discuss possible reasons for the large difference, which could arise due to real differences between planetesimal populations around A-type and Sun-like stars, or due to differences in the observables such as the 24-70µm colour temperature, which could change due to different grain blowout sizes for example. The point here is that the resulting C is sensible compared to previous results. If it were significantly different the plausibility of the in situ scenario would be questionable.
To create an in situ picture that explains significant levels of warm dust around both young and old stars clearly requires some combination of our two scenarios. We therefore make a simple combined model where systems start with the same initial distribution of warm excesses and also have a single random collision at some point during their lifetimes. This model has the same value of C = 1Myr as before, and is shown in Figure 7. Given that each respective scenario dominates the young and old populations with little effect on the other, that the model is in reasonable agreement with the observed luminosity function is unsurprising. One or two random collisions may be present in the young population of disks, and all old stars have initially massive belts that are now very faint. This model could be tuned a little more by specifying the fraction of systems that are initially massive and/or those that undergo collisions. However, the agreement shown in Figure 7 is good enough to illustrate that such a picture matches the observed distribution well. Therefore, we conclude that the Monte-Carlo evolution model could be a realistic (though not uniquely so) description of the ori- gins of warm dust. This model is of course highly simplified and largely empirical, and we discuss some issues further in Section 6.1.
Though we have shown that it is sensible, the assumption of in situ decay made by the above model may not be correct, with material being delivered to the terrestrial regions from elsewhere, in which case the apparently sensible value of the parameter C is coincidental. In a comet delivery scenario it is likely that the luminosity function is not the same as expected for 1/t evolution (though it could be). One test of our model is therefore to extrapolate to fainter levels, where detections are possible with much more sensitive observations of relatively few stars.

Extrapolation to faint exo-Zodi levels
The level of dust around a "typical" Sun-like star is a key unknown, and crucial information needed for the future goal of imaging Earth-like planets around other stars (e.g. Beichman et al. 2006;Absil et al. 2010;Roberge et al. 2012). There are several avenues for finding the frequency of warm dust at relatively low levels; either to directly detect it around an unbiased sample of stars (e.g. Millan-Gabet et al. 2011), or to make predictions based on the expected evolution of brighter dust. The former approach is clearly preferred as it yields a direct measure, but finding such faint dust is technically difficult and the focus of dedicated mid-IR instruments such as the Bracewell Infrared Nulling Camera (BLINC) at the MMT (Liu et al. 2009), the Keck Interferometer Nuller (KIN, Colavita et al. 2009), and the Large Binocular Telescope Interferometer (LBTI, Hinz 2009). 7 . These instruments are more sensitive to warm dust than photometric methods because their sensitivity is set by their ability to null the starlight and detect the astrophysical flux that "leaks" through the fringes (e.g. Millan-Gabet et al. 2011). The typical scales on which these mid-infrared instruments 10 -4 10 -3 10 -2 10 -1 10 0 10 1 F disk /F star at 12µm are sensitive is 10-100mas so are well suited to the terrestrial regions around the nearest stars. The latter approach of extrapolating results from less sensitive photometric surveys based on evolution models (i.e. the approach here) is easier in the sense that data for bright warm dust systems is readily available, but is of course more uncertain due to the unknown origin and evolution of these systems. A future goal is to fold the results from interferometry back into the modelling process. Figure 8 shows the 12µm luminosity function and one realisation of the combined initially massive and random collision model for 0-13Gyr old stars over a wider disk-star flux ratio range than Figure 6, using our adopted age distribution. The model underpredicts the number of brightest excesses, as this is the typical outcome (as shown in Fig 7), but varies considerably due to the small number of such disks. Below F disk /F⋆ ∼ 0.1 all realisations, and hence the predictions for fainter disks, are very similar.
The right side of the figure is covered by photometric surveys like the one presented here, whereas fainter disks (left side) can only be detected with more sophisticated methods, such as those employed by the KIN and LBTI. We have taken the approximate Solar System 12µm dust level as F disk /F⋆ = 5 × 10 −5 , and the LBTI sensitivity as 20 times this level (Hinz 2009). The predicted limit for exo-Earth detection by a Terrestrial Planet Finder (TPF)like mission is also shown (at ten times the Solar System level). 8 The limit set by the KIN survey (Millan-Gabet et al. 2011), which resulted in no significant detections, is shown, but merits further discussion. Their sample was classified into "high" and "low" dust systems, meaning those with or without detections of cool Kuiper belt analogues. Of the 25 systems observed, two were high dust systems (η Corvi and γ Ophiuchi) and warm dust was confidently detected around the former. η Crv is already well 10 -4 10 -3 10 -2 10 -1 10 0 10 1 F disk /F star at 12µm known to have distinct warm and cool dust belts (Wyatt et al. 2005;Smith et al. 2008), while only cool dust has been detected around γ Oph (Sadakane & Nishida 1986;Su et al. 2008). We use the other 23 low dust stars in calculating the KIN limit, with the caveat that if warm dust levels are positively (or negatively) correlated with cool dust levels, the limit will be slightly too low (or too high). Given that only ∼20% of stars are seen to have cool dust (i.e. many low dust systems could still have relatively high dust levels that could not be detected) this bias will not be very strong. For stars drawn randomly from our assumed age distribution, the predicted disk fraction is about 50% at the LBTI sensitivity. Future work may however show that our assumption of 1/t evolution is not exactly correct; the fraction decreases to 10% if the decay is instead t −1.5 so the prediction is not extremely sensitive to the decay rate. A very slow decay, such as that seen in the reference model of Gáspár et al. (2013) would appear to be ruled out by the KIN upper limit as it would predict many bright warm disks (such slow evolution is not expected for warm disks anyway, see the end of Section 5.1). It is in any case clear that an unbiased LBTI survey of just a few tens of stars would result in a strong test of our prediction, and has the potential for detection of a dozen or so disks if the evolution predicted by our model is correct. Such a result would be a great success for the instrument, but less promising for the study of exo-Earths; another prediction of our model is that roughly only a few tens of percent of stars are amenable to exo-Earth detection if our combined model is correct.
We have already shown that the different scenarios depend strongly on the age of the stellar population that is observed. This effect is shown explicitly in Figure 9, which shows the individual and combined distributions in a series of age bins. Here we have increased the number of stars in each age bin to 10 7 so that the curves are smooth. For both dust creation scenarios, the level of the faintest systems decreases with time as the stars become older. However, because random collisions can happen at any time, the luminosity function for these events is otherwise independent of time.
The initially massive population dominates the distribution for stars younger than about 100Myr, but becomes much less important as time goes on. For stars that are older than ∼500Myr the brightest excesses are already dominated by random collisions if our model is correct.
The solid lines show the combined model, which generally follow the fraction set by the model that dominates at a given disk to star flux ratio. At flux ratios where both components contribute to a similar fraction of systems (i.e. just above the maximum ratio for initially massive belts), the fraction is somewhat greater than expected just for random collisions. This increase arises because for these systems the emission from random collisions is augmented by dust that was initially bright (but is now much fainter).
This comparison of scenarios shows that for the oldest main sequence stars the predicted frequency of dust systems is lower than shown in Figure 8. The faintest dust may still come from those that were initially massive, but the bulk of the population are remnants of random collisions. That is, though the random collisions are very rare and are only bright for a short period of time, their decay through lower levels means that they may in fact dominate the bulk of the warm dust luminosity function for main-sequence stars.

In situ evolution
We have described a simple model of in situ evolution that reproduces the observed distribution of 12µm excesses reasonably well. This model has two components, an initially massive component and a random collision component, both of which are needed to simultaneously explain relatively frequent warm dust around young stars and extremely rare warm dust around ∼Gyr old stars (Fig 7). The predicted luminosity function changes strongly with age due to the different contribution from each component (Fig 9). However, the ∼Gyr old bright warm dust systems imply an underlying population of fainter warm dust that is present, regardless of age. Therefore, if our model is correct, the fraction of stars with warm dust levels sufficient to hinder exo-Earth detection may be a few tens of percent, even if only old stars are targeted.
By using the term "initially massive" for the disk component that is bright at early times, we have implied that it is a pumpedup version of our own Asteroid belt. While such an interpretation may be correct, an alternative possibility is that large amounts of warm dust are generated in the first ∼100Myr as a by-product of terrestrial planet formation (Kenyon & Bromley 2004;Lisse et al. 2008Lisse et al. , 2009Jackson & Wyatt 2012). Such a scenario evolves similarly to an initially massive Asteroid belt analogue, but is different in a few specific ways; the dust levels may initially increase as large planetary embryos that stir the planetesimals to destructive collision velocities form, dust may be visible for longer as mass is released from embryos in planet forming impacts, and planetesimals and dust may also be depleted by dynamical interaction with the embryos. These differences mean that the emission from terrestrial planet formation may not simply decay as 1/time, providing a potential way to distinguish it from initially massive Asteroid belt analogues (e.g. Currie et al. 2008).
One limitation of our model is that processes beyond the basic collisional evolution will modify the model luminosity functions. As noted above, dynamical depletion may follow catastrophic collision events, perhaps with an initially optically thick phase (Jackson & Wyatt 2012). Significant variability and nonlinear decay could also arise at early stages due to changes in dust mineralogy. For example, the silicate peak near 12µm could decrease in strength after a collision event as the smallest grains are removed by radiation pressure. The overall dust level may also vary following the initial creation event as the system settles into dynamical and collisional equilibrium. Such processes might be responsible for the observed luminosity function being somewhat flatter than the models (Figs 6-8).
A more basic limitation is our implementation of random collisions. With only a few confirmed warm dust sources around ∼Gyr old stars there is little leverage to constrain models, so this part of our model is necessarily simple, and not for example linked to the mass remaining from initially massive belts. The largest simplification is that uniform probability of a collision over 0-13Gyr assumed here is unlikely to be the case. Wyatt et al. (2007a) showed that if these collisions are between the largest objects in a continuous size distribution (i.e. just the largest objects in an Asteroid belt analogue), then the probability of detecting a collision drops as 1/time 2 . An alternative possibility is that the collisions are between individual large objects (e.g. planets or planetary embryos) that simply did not collide earlier, or were perhaps subject to some late instability. These objects may have originally been part of a continuous size distribution, but were "stranded" when only a few remained and their chances of destruction became much lower (see Kennedy & Wyatt 2011). However, in these cases collisions are more likely to occur at earlier times, though perhaps with a weaker dependence compared to the large asteroid scenario just described.
That the two "old" stars hosting warm dust, BD+20 307 and possibly HD 150697 are between 1 and 4Gyr old can be used to illustrate that there is a rough constraint on the time dependence of random collisions. This constraint is uncertain because in our sample there are 5 sources of unknown age, and the ages of some others are poorly constrained. However, a 1/time 2 dependence is not favoured because the existence of warm dust around the ∼Gyr old stars implies a significant number should exist around stars with ages from ∼100-1000Myr, which are not observed (see also Melis et al. 2010). The argument is similar for 1/time, but given the uncertainties is harder to rule out. Similarly, for the random collision model with uniform time dependence a few detections are expected in the 4-13Gyr age bin and these are not observed either. Therefore, if the random collision scenario is correct the time dependence of collisions must decrease with time fairly weakly, if at all.

Comet delivery
While we so far have discussed our observational results in the context of an in situ model of evolution, the observed dust may originate elsewhere. Specifically, material may be delivered from more distant regions after interacting with planets, thus replenishing the mass that resides in the warm belt (e.g. Bonsor & Wyatt 2012;. With this kind of evolution, the observed emission need not decay as 1/t, and the luminosity function need not have the slope of -1 seen in Figure 8. Such a scenario may be acting to replenish the Zodiacal cloud in the Solar System. The Zodiacal dust was first suggested to arise from particles dragged radially inwards from the Asteroid belt by Poynting-Robertson drag (Kresak 1980). More recently however, Nesvorný et al. (2010) have shown that most of the Zodiacal cloud as seen by IRAS can be explained by a sufficient level of spontaneous disruptions of Jupiter-family comets (though the exis-tence of dust bands means that Asteroid families still contribute some dust, Dermott et al. 1984). Because Jupiter-family comets are themselves thought to originate in the Edgeworth-Kuiper belt (Levison & Duncan 1997), the location in which the dust is observed (i.e. near Earth) is only linked to the source region by comet dynamics, which are dominated by the Solar System's four giant planets.
In this kind of scenario, the amount of warm dust is plausibly linked to the number of comets in the source region; a change in the number of source objects should be directly reflected in the number scattered inwards, and therefore be reflected in the observed dust level. Therefore, if the source region is in collisional equilibrium at all sizes (i.e. Mtot ∝ 1/t), then the observed dust level should decay as roughly 1/t, with some additional contribution from dynamical depletion (that may be strongest at early times). While the possibility of 1/t decay might appear consistent with the model of section 5, a very different value for the decay timescale C is expected because collisional evolution is much slower at larger radial distances (i.e. at least Gyr instead of Myr, Trilling et al. 2008). Therefore, a naive picture of warm dust occurring in all systems exclusively due an uniformly initially bright level of steady delivery throughout the main-sequence does not work because the brightest systems are not predicted to be rare; their Gyr decay times mean that they should in fact be ubiquitous. The problem with this picture is that comet delivery will not always be initially bright; the level is set by the dynamics of individual systems. It may be that near-constant comet delivery is the dominant warm dust origin and the luminosity function set entirely by the distribution of systems with the appropriate dynamics (Bonsor & Wyatt 2012).
Though we have argued for the slow decay of reservoirs from which comets are delivered, it is also possible that a comet delivery scenario could mimic our initially bright asteroid belt analogues or terrestrial planet formation due to the planetary system dynamics in the first ∼100Myr. A newly formed planetary system will have small bodies scattered throughout, and in addition to those in young Asteroid and Kuiper belt analogues, many will also be on unstable orbits. These unstable objects could contribute to warm dust levels as they are cleared over the first 100Myr or so, thus appearing to decay in the same way as initially massive belts.
Clearly, more work is needed to explore the likelihood of in situ and cometary origins for warm dust. For example, it could be that warm dust originates in some cases from an impact on an object orbiting at ∼1AU by an interloping comet (e.g. Lisse et al. 2012), in which case the material is terrestrial in origin and evolves as such (i.e. in situ), but whose probability of occurrence is more closely related to the comet delivery rate.

Solar System comparison
A basic prediction of our model is that the level of Zodiacal dust seen in the Solar System is roughly a factor of ten fainter than is typical (Fig 8). Given the uncertainties discussed above this predicted rarity is not an issue for our model, but worth some consideration.
If the bulk of our Zodiacal light comes from the Asteroid belt, the lower than predicted dust level could arise in several ways. The initial distribution may simply have more faint belts than we assumed, reflecting the range of outcomes dictated by planetesimal growth in the protoplanetary disk or other processes such as planet migration. For example, in the Solar System the Asteroid belt is thought to have been severely dynamically depleted when Jupiter formed (Petit et al. 2001), and a more general inner Solar System depletion may have occurred if Jupiter was at some point closer to the Sun that it is now (Walsh et al. 2011). This difference could be accounted for by weighting the initial brightness distribution towards lower levels, with a slower disk decay rate to keep the same number of observed warm disks.
If the bulk of our Zodiacal light comes from comet delivery, it may again be that the true initial distribution is different to that assumed here, or that the evolution is not continuous as assumed in our model. For example, it could be that the Solar System's Zodiacal cloud evolution was interrupted by a late dynamical instability that ejected of most of the primordial Kuiper belt (Tsiganis et al. 2005). Removal of most of the source region mass would cause the warm dust level to drop in a way not accounted for by our model. While it is currently unknown if the Kuiper belt has a relatively low mass (e.g. , the relative faintness of the Zodiacal cloud in the comet delivery picture could be explained if such depletion events are relatively uncommon around stars in general. Naturally, many of the possible extensions to our model could involve the effect of planets on the parent planetesimals, whether they be asteroids or comets. Therefore, in addition to looking for correlations between warm and cold dust as a signature of scattering by planets , an important goal is to look for correlations between warm dust and planets. While a tentative positive correlation has been seen between low-mass planets and cool Kuiper belt analogues , no trends have yet been seen for warm dust (e.g. Morales et al. 2012).

SUMMARY AND CONCLUSIONS
We have presented an analysis of the bright end of the 12µm warm dust luminosity function as seen by WISE around nearby Sun-like stars. We report the possible detection of 6 new warm dust candidates, HD 19257, HD 23586, HD 94893, HD 154593, HD 165439, and HD 194931. There is a clear preference for excess emission around young stars, with at least 14/22 systems younger than 120Myr, despite the ∼10Gyr stellar lifetime. However, one of the dustiest systems is ∼1Gyr old, and another may be ∼Gyr old, suggesting that while youth is an important factor for setting the brightness of warm dust, it is not the sole factor.
Using a simple in situ evolution model to interpret the observed luminosity function, we show that neither a picture where all systems have initially massive disks, nor one where all excesses are due to randomly timed collisions, can explain the fact that warm dust is observed around both young and old stars. A simple combination of these two in situ scenarios reproduces the observed luminosity function. However, we cannot rule out that warm dust is in fact dominated by comet delivery. In this scenario comets are scattered from a more distant reservoir, and reliable luminosity function predictions cannot be made because they depend on the distribution of systems with the appropriate dynamics.
We have made simple predictions for the number of exo-Zodi that might be detected by an unbiased LBTI survey. For the expected sensitivity we anticipate a handful of detections and a strong test of our evolution model. Such tests are a crucial step towards understanding the origins and evolution of exo-Zodi, and ultimately the imaging of Earth-like planets.

ACKNOWLEDGEMENTS
We would like to thank the LBTI Science and Instrument teams for fruitful discussions, Fei Dai for preliminary work, and Olivier Absil for a thorough referee report with valuable criticisms and suggestions that improved the manuscript. This work was supported by the European Union through ERC grant number 279973. This research has made use of the following: The NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. Data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. Data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. The SIMBAD database and VizieR catalogue access tools, operated at CDS, Strasbourg, France.

APPENDIX A: NOTES ON INDIVIDUAL SOURCES
Notes for the 96 sources identified as having 12µm excesses based on the criteria outlined in sections 2 and 3. The 25 warm dust systems are shown in bold face. Stars noted as bright and very bright are those where saturation in W1 becomes an issue, meaning flux densities greater than ∼1-10Jy at 3µm. 12. HIP 20612 (CPD-77 172B): About 10" from the primary HIP 20610 (HD 29058). No excess; high W 1 − W 3 caused by low W 1 measurement (which is also marked as variable in the catalogue  Walker & Wolstencroft 1988;Zuckerman & Becklin 1993). Age is about 8Myr. Sometimes classed as a disk in "transition" between the protoplanetary and debris disk phases (Furlan et al. 2007), with the argument perhaps strengthened due to the recent detection of H2 (Yang et al. 2012 (de Zeeuw et al. 1999). The age is therefore about 17Myr. Chen et al. (2011) find a large 70µm flux, suggesting that the 12µm excess is probably the Wien side of cooler (67K) emission. The detailed SED model finds that the 12µm excess is only 2.4σ significant, so this source is not considered as a true warm dust source.
50. HIP 63975 (HD 113766A): Well known warm dust source in Lower Centaurus Crux (de Zeeuw et al. 1999;Oudmaijer et al. 1992;Lisse et al. 2008). The age is therefore about 17Myr. The excess appears around the primary star in this binary system, the companion is about 160AU distant so too distant to affect the warm dust unless the eccentricity is very high (see Smith et al. 2012 Rizzuto et al. (2011Rizzuto et al. ( , 2012 with 69% probability. Hence age is likely 17Myr. Disk temperature is 238K. 55. HIP 64995 (HD 115600): Member of Scorpius-Centaurus in Lower Centaurus Crux (de Zeeuw et al. 1999). The age is therefore about 17Myr. Chen et al. (2005) find a large 70µm flux, suggesting that the 12µm excess is probably the Wien side of cooler (120K) emission. The detailed SED model finds that the 12µm excess is only 1.4σ significant, so this source is not considered a true warm dust source.
56. HIP 65875 (HD 117214): Member of Scorpius-Centaurus in Lower Centaurus Crux (de Zeeuw et al. 1999). The age is therefore about 17Myr. Chen et al. (2011) find a large 70µm flux, suggesting that the 12µm excess is probably the Wien side of cooler (110K) emission. The detailed SED model finds that the 12µm excess is only 2σ significant, so this source is not considered a true warm dust source.
57. HIP 67497 (HD 120326): Member of Scorpius-Centaurus in Lower Centaurus Crux (de Zeeuw et al. 1999). The age is therefore about 17Myr. Chen et al. (2011) find a large 70µm flux, suggesting that the 12µm excess is probably the Wien side of cooler (105K) emission. The detailed SED model finds that the 12µm excess is only 1σ significant, so this source is not considered a true warm dust source.
58. HIP 67970 (HD 121189): Member of Scorpius-Centaurus in Upper Centaurus Lupus (de Zeeuw et al. 1999). The age is therefore about 15Myr. Excess found at 24µm by Chen et al. (2011), but no longer wavelength detection. WISE 12µm excess is only 2.5σ significant, so this source is not considered a true warm dust source. Excess in SED appears reliable. WISE images show possible contamination due to filamentary background at both W3 and W4, though whether the excess is caused by this is hard to tell. The W3 excess is 3.6σ, so the background contribution would not need to be large to cause the excess to be significant. Therefore we do not include this source as a warm dust candidate.
74. HIP 81870 (HD 150697): Noted as a possible warm dust source from IRAS photometry by Fajardo-Acosta et al. (2000), but with the caution that the star lies on the edge of the ρ Oph starforming region. With a distance of 133pc, it may in fact be a young star associated with ρ Oph, but was assigned an age of 3.2Gyr by Holmberg et al. (2009). We find that the IRAS 12µm excess is 2.98σ significant and that the W3 images appear clean, so the excess appears reliable.

APPENDIX B: SAMPLE AGE DISTRIBUTION
While stars of a given spectral type may reasonably be assumed to have random main-sequence ages, our sample spans a sufficiently wide range of spectral types that the main-sequence lifetimes τMS vary significantly. Because our stars are selected based on their Tycho-2 BT − VT colour, to derive the age distribution we want a relation between this colour and τMS.
About 10,000 of our 24,174 stars actually have ages derived by The Geneva-Copenhagen survey of the Solar neighbourhood (Nordström et al. 2004;Holmberg et al. 2009). The ages are derived using isochrones, and therefore have the weakness that stars not sufficiently evolved from the zero-age main-sequence cannot be assigned ages. Indeed, the authors state "that the true age distribution of the full sample cannot be derived directly from the catalogue; careful simulation of the biases operating on the selection of the stars and their age determination will be needed to obtain meaningful results in investigations of this type." Despite this lim- itation however, the catalogue does provide a distribution of ages as a function of colour, shown in Figure B1, which can be used to verify an empirical relation between τMS and BT − VT.
We first plot absolute luminosities as derived by SED fitting against colour for main-sequence stars from the Unbiased Nearby Stars (UNS) sample (Phillips et al. 2010) in Figure B2. The correlation in this HR diagram is reasonably well described as L⋆ ∝ 10 −1.5(B T −V T ) . This relation, combined with τMS ∝ M⋆/L⋆ (i.e. assuming all mass is turned into energy) and L⋆ ∝ M 3.5 ⋆ yields τMS ∝ 10 B T −V T . This dependence is only approximate however. Using it as a starting point, with some experimentation we found that provides a reasonable match to the maximum ages for early-types in Figure B1, while giving an age of 9Gyr at BT − VT = 0.68, corresponding approximately to the Sun. We add a further restriction that ages must be less than 13Gyr.
To derive the age distribution of our sample now requires converting the known distribution of colours. This conversion is done by setting up a series of colour bins, and selecting ages randomly between zero and τMS for the stars in each bin. This process results in the age distribution for our sample shown in Figure 2.
To generate this distribution numerically we make use of a fourth order polynomial fit to the histogram, with coefficients of 0.171124, -0.0528395, 0.00686227, -0.000410586, and 9.29393e-06. We generate a random age in the desired range, and another random number between 0 and 0.1035. If the second random number is lower than the polynomial function evaluated for that age, or if the polynomial is greater than 0.1035, we keep the age. This process is repeated until the desired number of ages is acquired.
Our analysis has been done with both a naive uniform distribution of ages between 0-10Gyr and the more detailed distribution derived above, and our conclusions are the same for either distribution.   Figure C4. 25 warm dust systems cont.