{'fraction_non_alphanumeric': 0.03336483244893509, 'fraction_numerical': 0.0247328632696082, 'mean_word_length': 4.6532771727963835, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 4, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Variables, periodic variables and contact binaries in WISE

The time-series component of WISE is a valuable resource for the study of stellar variability. We present an analysis of an all-sky sample of ~450,000 AllWISE+NEOWISE infrared light curves of likely variables identified in AllWISE. By computing periodograms of all these sources, we identify ~63,500 periodic variables. Of these, ~44,000 are short-period (P<1 day), near-contact or contact eclipsing binaries, many of which are on the main sequence. We use the periodic and aperiodic variables to test computationally-inexpensive methods of periodic variable classification and identification, utilizing various measures of the probability distribution function of fluxes and of timescales of variability. Combining variability measures from our periodogram and non-parametric analyses with WISE infrared colors and absolute magnitudes, colors and root-mean-square variability from Gaia yields a powerful engine for the identification and classification of periodic variables. Furthermore, we show that the effectiveness of non-parametric methods for the identification of periodic variables is comparable to that of the periodogram but at a much lower computational cost. Future surveys can utilize these methods to accelerate more traditional time-series analyses and to identify evolving sources missed by periodogram-based selections.


INTRODUCTION
Binary systems are relevant for a wide array of astronomical phenomena. They have been used in cosmology as distance indicators (Riess et al. 2011), in stellar astrophysics as proving grounds for precision stellar evolutionary models (Pietrinferni et al. 2004), and can even host planets (Doyle et al. 2011). Binary systems have been linked to transients such as Luminous Red Novae (Tylenda et al. 2011;Kasliwal 2012) and they are thought to serve as progenitors for some of the most fascinating objects in the universe -ultra-compact binaries of white dwarfs, neutron stars and black holes, the tantalizing source population for type Ia supernovae, kilonovae, and gravitational waves (Weisberg et al. 2010;Knigge et al. 2011;Maoz et al. 2014;Postnov & Yungelson 2014;Brown et al. 2016;Belczynski et al. 2016;Smartt et al. 2017;Cowperthwaite et al. 2017;Abbott et al. 2017;Temmink et al. 2020).
The growth in variability datasets has been accompanied by a rise in complexity of methods for variable classification. Traditionally, variable objects are classified based on the similarity of their light curves and colors to known variable prototypes (Gaia Collaboration et al. 2019). Often, a time-series analysis tool such as a periodogram is run to differentiate between peri-odic and aperiodic variables. This step requires a clear understanding of the effects of cadence on period recovery and a careful weighing of the pros and cons of different period-search algorithms. The results of the periodogram-based analysis are often taken, together with measures of light curve morphology and other characteristics of the source (e.g., color) and used as inputs into a classifier (e.g. Debosscher et al. 2007;Richards et al. 2011;Dubath et al. 2011;Richards et al. 2012;Masci et al. 2014;Jayasinghe et al. 2019a,b;Eyer et al. 2019). There are publicly available light curve classifying engines that rely on machine learning techniques (e.g. Kim & Bailer-Jones 2016). Periodograms, light curve fitting, and machine learning classification are potent tools and will continue to be important to the astronomical community. Nevertheless, one downside of these methods is that they are often survey-dependent, difficult to physically interpret, and involve a lengthy learning curve and a lot of computational power in order to implement.
To reduce computational overhead, some previous works have used non-parametric variability measures (Kinemuchi et al. 2006;Palaversa et al. 2013;Drake et al. 2013Drake et al. , 2014aDrake et al. , 2017Torrealba et al. 2015;Findeisen et al. 2015). We build upon this work and develop new, non-parametric light curve analysis techniques and apply these techniques on data from the Wide-Field Infrared Survey Explorer (WISE; Wright et al. 2010;Mainzer et al. 2011). The all-sky and long-term coverage, decreased effects of extinction compared to optical wavelengths, and non-uniform cadence probing a wide range of variability timescales (Hoffman et al. 2012) make WISE a unique probe of Galactic variability. Previously, Chen et al. (2018) used WISE to select ∼50,000 variable candidates of which ∼ 42,000 were binaries. Here, we expand on this work and present an analysis of ∼450,000 light curves of variables selected using AllWISE variability metrics based on r.m.s. flux variations. Our larger period search grid allows us to detect short-period objects that were missed by Chen et al. (2018). We also employ a different period-search algorithm, use data from a more recent release, and cross-match our results with Gaia DR2. In the end, we produce an all-sky sample of ∼ 63,500 periodic variable candidates, of which ∼ 56,000 are binaries.
In Section 2, we discuss the WISE data and our timeseries analysis. In Section 3 we explore the contents of our periodic variable selection and display the results of various cross-matches. In Section 4 we introduce our non-parametric measures and use them to classify periodic variables. In Section 5 we discuss main-sequence (MS) binaries, short-period objects, extra-galactic and young stellar object variability, and the application of non-parametric methods to the identification of periodic variables. We conclude in Section 6. Throughout this paper, WISE magnitudes are quoted on the Vega system, and the following conversions apply: W i AB = W i Vega +∆m i , with ∆m i = (2.699, 3.339, 5.174, 6.620) 1 .

WISE Mission
The Wide-field Infrared Survey Explorer (WISE) was launched in December 2009 and conducted observations of the entire sky in bands centered on 3.4, 4.6, 12 and 22 µm (W1, W2, W3, and W4 respectively) until 2010 when it ran out of coolant (Wright et al. 2010). From 2010 to 2011 it conducted observations in the 3.4 and 4.6 µm bands as a part of the NEOWISE post-cryogenic mission (Mainzer et al. 2011). The spacecraft was in hibernation from the end of the post-cryogenic mission until 2013 when it was reawakened as a part of the NE-OWISE Reactivation (NEOWISE-R) mission (Mainzer et al. 2014). The AllWISE data release includes data from both the original mission and the post-cryogenic mission as of 2013 2 , whereas the NEOWISE Reactivation 2019 Data Release includes all of the data from the time that the spacecraft was awakened from hibernation until 2019.
The WISE spacecraft orbits the Earth with a period of ∼5700 seconds (∼0.066 days=1.6 hours) 3 on a polar orbit -near the dividing line between night and day -and always looks away from the Earth. Every six months, it images the same portion of the sky and obtains at least eight passes on each point of sky due to the partial overlap of the field of view on consecutive orbits. The sources outside the ecliptic are more frequently observed. The cadence of the data is such that every six months there is a collection of ∼10 data points that are each spaced one orbit apart.
The data pre-processing and extraction of the light curves are different for the AllWISE and NEOWISE releases due to differences in the mission during the two phases. AllWISE stacks all the scans, identifies the objects, and then measures the per-scan magnitude of each object at a fixed position. In contrast, NEOWISE identifies objects and measures their photometry in individual scans without stacking. Mainzer et al. (2014) find systematic changes in W1 between AllWISE and NEOWISE-R to be 0.01 magnitude for sources with 8<W1<14 magnitude and on the order of 0.1 magnitude for sources with 14<W1<15 magnitude. Outside this range there are magnitude dependent systematic offsets between AllWISE and NEOWISE data. We limit our sources to 8<W1<15 mag to ensure concordant All-WISE and NEOWISE measurements and to exclude saturated sources (Nikutta et al. 2014;Mainzer et al. 2014). We also explicitly apply a cut to ensure that the difference between the mean magnitude of AllWISE and NE-OWISE be less than 0.15 because the time-series analysis is not reliable in the case of a large magnitude offset. A typical flux measurement in the range 8<W1<15 mag has a error of ∼ 0.025 mag. We incorporate these flux error measurements into our analysis.

WISE Light Curves
WISE reports a measure of the flux variation in each band based on AllWISE single-epoch photometry. Each source is assigned a single-digit var_flg ranging from 0 to 9 in each band, such that the probability that the object's flux does not vary in said band is ∝ 10 −var flag (Hoffman et al. 2012).
We select WISE variable sources having W1 var_flg ≥ 6. We only consider sources having cc_flags == 0 to ensure that there are no imaging artifacts and ext_flag ≤ 1 to ensure that it is not an extended source. After applying these quality cuts, we download ≈ 500,000 light curves from the AllWISE multi-epoch photometry table  and the NEOWISE-R single exposure source table using  a 1 matching distance. To eliminate redundancies and possible extraneous matches in the table, the allwise_cntr listed in the NEOWISE-R single exposure table is mapped to the source_id_mf in the AllWISE multi-epoch photometry table. 4 We exclude sources whose NEOWISE-R data is mapped to more than one AllWISE source. We also exclude AllWISE sources with no corresponding NEOWISE-R data because the number of AllWISE-only observations is not sufficient for our time-series analysis.
For the AllWISE multi-epoch photometry, the following quality cuts were applied: saa_sep>5.0 deg (image outside of the South Atlantic anomaly), moon_masked == 0000 (frame unaffected by light scattered off the moon), and qi_fact > 0.9 (only the highest quality frames) 5 . In addition, points with null photometric measurement uncertainty or null for the reduced χ 2 of the W1 profile-fit are excluded. For NEOWISE analogous cuts are applied 6 and we also exclude points with with null W1 profile fit signal-to-noise ratio. The number of data points reported in our tables is the total remaining after these quality cuts are applied. Example python code for lightcurve download with quality flags is available through Github 7 .

Time Series Analysis
The periodogram is a time-series analysis tool that allows for the location and characterization of periodic signals. For this paper, we use the multi-harmonic analysis of variance (MHAOV) periodogram (Schwarzenberg-Czerny 1996), which has been shown to be highperforming in comparison to other algorithms (Graham et al. 2013). We fit the data with periodic, orthogonal functions and use a statistic, θ, which is the ratio of the squared norm of the model over the squared norm of the residuals (Schwarzenberg-Czerny 2003), to quantify the quality of the fit (Schwarzenberg-Czerny 1998). The fitting procedure is carried out for a grid of test frequencies and a periodogram shows the dependence of the statistic value on the test frequencies. Figure 1 shows an example periodogram in frequency space with a corresponding phase-folded light curve.
For three model parameters, the MHAOV periodogram is statistically equivalent to the classic Lomb-Scargle periodogram (Lomb 1976;Scargle 1982;Ferraz-Mello 1981;Schwarzenberg-Czerny 1999). We use 5 model parameters to allow for better sensitivity to anharmonic oscillations. This corresponds to fitting the data in real space with with a Fourier series of 2 harmonics (Schwarzenberg-Czerny 2003Graham et al. 2013;Lachowicz et al. 2006).
We adopt a frequency grid spacing of 0.0001 days −1 (Graham et al. 2013). We search the frequency range 0.1 to 20 days −1 . This extension beyond the traditional Nyquist limit of ∼ 7.6 days −1 that corresponds to uniform sampling at the WISE orbital cadence is justified when the data is not uniformly spaced (see subsection 2.6 for further discussion).

Periodogram Peak Significance
Upon observing a periodogram peak, we seek to reject the null hypothesis that the light curve in question consists of pure noise and that the highest peak in the periodogram results from the chance alignment of random  errors. Although some inroads have been made toward an analytical understanding of peak significance (e.g. Horne & Baliunas 1986;Koen 1990;Schwarzenberg-Czerny 1999), the results depend in complex ways on the nature of the data and the chosen test frequencies.
We adapt the Monte Carlo method of Frescura et al. (2007Frescura et al. ( , 2008 to quantify peak significance. We generate 12,000 white-noise light curves that have means, flux deviations from mean, observing cadences, and individual data point photometric uncertainties that are representative of our actual sample. To do this, we start with 12,000 randomly-selected WISE light curves. For each light curve, there is an array of time measurements, t i , an array of magnitude measurements, m i , and an array of photometric uncertainties on each data point, σ(m i ) where i ranges from 1 to the number of data points in the light curve, N LC . For each light curve, we calculate the mean of the magnitude distribution, m , and define a measure of its width, w, as the difference between the 84.2 th and 15.8 th percentile. Next, we create a Gaussian distribution with mean m and standard deviation 0.5w and randomly draw N LC values from this distri-bution to create a new array of Gaussian (white) noise data, noise i . By substituting noise i for m i for each light curve, we create 12,000 white-noise light curves.
We run a MHAOV periodogram on each of the whitenoise light curves using the same frequency grid and number of model parameters that were used on the actual data. Figure 2 shows the resultant distribution of maximum statistic values from the periodogram (θ), as a function of the number of observations in the light curve (N LC ). When N LC is low, the light curve contains less information, it is harder to distinguish periodic and aperiodic signals, and it is expected that a higher statistic value is required to credibly reject the null hypothesis. To make Figure 2, we bin the data based on N LC . Every bin has ∼1,500 simulated light curves save the last one on the right which has ∼1,000. In each bin, we study the resultant empirical cumulative distribution function of the calculated maximum statistic values and estimate the 99.9 th percentile. Then we fit a power law model of the form θ cutof f (N LC ) = A(N LC − N 0 ) k to estimate the appropriate cut-off value as a function of N LC to reject the null hypothesis. The best-fit values are A = 154±57, N 0 = 14.9 ± 4.7, and k = −0.344 ± 0.082. For sources that lie above this best-fit line we reject the null hypothesis that the peak results from the chance-alignment of random errors for a white-noise light curve. The bin with the highest number of light-curve data points (i.e. N LC >∼ 300) had the fewest number of sources spread over the largest range. To account for the less-precise estimate of the cumulative distribution function in this region and for the fact that the power law fit approaches zero as N LC increases, we adopt a minimum statistic value of θ = 10 to reject the null hypothesis.

Completeness
The above cut on the maximum statistic seeks to limit the false alarm probability; however, it does not say much about the probability of periodicity itself nor the incidence of truly periodic signals that are passed over (VanderPlas 2018). The completeness, i.e. the proportion of truly periodic signals we expect to detect given our method and data, is a complex function of the nature of periodicity (Schwarzenberg-Czerny 1999). It depends on the magnitude (worse for very faint and very bright objects), amplitude (worse for small amplitude), signal shape and phase (worse for short-duration pulses, especially if they occur in between observing epochs), period (worse for extremely short and long periods), and number of observations (worse for fewer observations).
Completeness is not the main focus of this project, and we do not attempt to characterize completeness across all of these parameters. Instead, we use simulations to gauge our ability to detect contact binaries with nearsinusoidal light curves as a function of signal amplitude and period. We randomly select 100 light curves that have W1 magnitudes between 10 and 14, the range best probed by WISE. We then simulate a sinusoidal signal centered at the W1 magnitude of the light curve and sample it at the original cadence, preserving the individual data point uncertainty associated with each timestamp. In the first simulation, we fix the phase and choose an amplitude characteristic of our recovered periodic variables (0.4 mag). We pick a starting period of 0.051 days and calculate the MHAOV periodogram, phase-fold the result and see if the source would have been classified as periodic. We iterate over all light curves for a given period and then increase the period and repeat until the period exceeds a value of 2 days. We choose periods between 0.05 and 2 days because the majority of our periodic sources are contained in this range. In this interval, aside from the periods that coincide with the WISE observing cadence, the recovery rate is nearly constant at ∼ 77% of the simulated periodic variables. Next, we repeat the same procedure but this time fixing the period at a value of 0.22 days (a value typical of our close binaries) and varying the amplitude between value of 0.05 magnitude and 2 magnitude. For amplitudes of 0.2 mag and above, the distribution is constant with a median recovery rate of 82%. Below 0.2 mag, the recovery rate drops to 14% at a amplitude of 0.05 mag. In summary, we estimate that we can detect ∼ 80% of near-sinusoidal signals with periods between 0.05 and 2 days and peak-to-peak amplitudes above 0.2 magnitude.

Period Uncertainty
After rejecting the null hypothesis that the object is aperodic, the next task is to select the correct period. The observed light curve is a product of the continuous underlying signal and the discrete and unevenly-spaced window function (i.e. observational cadence). Aliasing, or the correlation between frequencies equidistant from one-half of the the inverse sampling rate (Nyquist frequency) is familiar in the case of evenly-spaced sampling. Uniform observations at the WISE satellite period (∼ 95 minutes) would correspond to a Nyquist frequency of ∼ 7.6 days −1 . Deviations from uniformity dampen the effects of aliasing and allow for the detection of periodic components above the traditional Nyquist limit (Eyer & Bartholdi 1999;Koen 2006).
Despite the orthogonality of the multi-harmonic periodogram fit at each individual frequency, the fits on any set of frequencies are generally not independent in the case of irregular sampling (Schwarzenberg-Czerny 1998). This means that the statistic values at different frequencies can be correlated. Sometimes these correlations manifest as separate peaks. Compared to other algorithms, we find the MHAOV periodogram to be effective at damping secondary peaks across the range of probed frequencies. For the purposes of this analysis, we use the apparent period (P app ) as given by the periodogram. In the case of symmetric light curves, such as an eclipsing binary composed of similar stars, P app corresponds to half of the orbital period. Correlations also cause the periodogram peak to have finite width and limit the precision with which the frequency can be determined from the peak. Generally, the uncertainty is a function of signal characteristics, the model used, the temporal baseline, and the quantity and quality of data (see Hartman et al. 2008;Lachowicz et al. 2009;Harding et al. 2013 for examples of different error estimation strategies).
To refine our frequency measurements and to estimate their error, we repeat the periodogram procedure with a refined frequency grid in the vicinity of the main periodogram peak. Our frequency measurement is the position of the likelihood peak on the refined frequency grid. We use the width of the peak on the fine grid and the surrounding background noise level to estimate the frequency error σ f req (Schwarzenberg-Czerny 1991, 1995, 1996. To check this frequency error estimation scheme, we compare with errors derived using emcee (Foreman-Mackey et al. 2013). To have realistic errors for real-istic (non-sinusoidal) light curves, we take a random sample of 100 light curves, phase-fold them, and fit them with a Fourier series using LMFIT (Newville et al. 2014). Then, holding the Fourier coefficients fixed, we fit the light curve in the time domain with emcee using the phase, frequency, and constant offset from zero as parameters. The fitting is done in this way because the WISE cadence made simultaneously fitting the Fourier coefficients, phase, frequency, and constant offset intractable. To get the frequency distribution, we marginalize over the constant and the phase. As seen in Figure 3, the MC-derived errors, σ M C , are correlated with the periodogram-derived errors, σ f req , but are systematically higher by ∼ 25%. In what follows we use a conservative error estimate by multiplying the periodogram-derived error by a factor of 1.25. The period errors are on the order of ∼ 10 −6 -10 −7 days, but they do not capture the error of picking the wrong periodogram peak entirely and thus should be used with caution (VanderPlas 2018).

Selection of Periodic Variables
To select periodic sources, we first require that the maximum statistic of the MHAOV periodgram, θ max , exceeds the cutoff value given by the power-law fit of Figure 2 for the appropriate number of points. Our cuts on false alarm probability are less strict than those of some other surveys that identify periodicity. This offers us the possibility of capturing sources undergoing interesting orbital evolution. For a higher confidence sample, a more restrictive cut on the maximum statistic value should be used. In addition to the cut on the maximum statistic value, we require phase coverage of at least 90%. We estimate phase coverage by dividing the phase-folded light curve into 20 equally sized bins and then calculating the percentage of bins that contain at least 1 data point. The requirement that phase coverage is at least 90% implicitly requires that the sources have at least 18 data points. Finally, we exclude sources whose highest periodogram peak is at 1 or 0.5 times the orbital frequency of the WISE satellite (f peak / ∈ [7.52, 7.68] days −1 and f peak / ∈ [14.65, 15.7] days −1 ). With the above cuts, we find 63,526 periodic variable candidates. Figure 4 shows the distribution of apparent periods and apparent periods over the period errors for these sources. best-neighbor WISE/Gaia cross-match catalog of Marrese et al. (2019). We find 188,043 matches out of a total of 454,103 variable objects (∼ 41.4%) which is similar to the total percentage of AllWISE sources that have a best-neighbor cross-match (39.83% -see Marrese et al. 2019). Of our 63526 periodic variable candidates, 54,756 or ∼ 86.1% have a best-neighbor Gaia cross-match. The higher match rate for periodic variables is due to the fact that the aperiodic variables -many of them stochastically varying young stellar objects -tend to be redder, with a median W1-W3 = 0.916 mag as compared to the periodic variables with W1-W3 = 0.334 mag. As a further check, we cross-match our sample with the WISE young stellar object catalogs of Marton et al. (2016) and find 44,196 matches of which only 289 are flagged as periodic.
We use the Gaia cross-match and considerations of the limitations of WISE to get a sense of the catalog contents and make the case that the sample is dominated by eclipsing binaries. We search for periods between 0.05 and 10 days, but due to the WISE observing cadence, the sensitivity drops for periods above 2 days. As a result, we do not expect many long-period variables in the more luminous part of the color-magnitude diagram. The period sensitivity limit combined with the crowding in the Galactic disk cause us to detect few classical Cepheids. We expect to detect few δ Scuti variables because their variability amplitudes are too low to be reliably detected (Murphy et al. 2019). The photometric sensitivity of WISE also limits our ability to detect variability on the white dwarf sequence.
For periodic objects with a Gaia cross-match, in order to have robust absolute magnitudes and colors, we require: In addition, we restrict phot_bp_rp_excess_factor in accordance with Gaia Collaboration et al. (2018b). After applying these cuts we retain 37,962 sources. Figure 5 shows the Gaia color-absolute magnitude diagram of the WISE periodic variables that satisfy the above quality cuts. They are color-coded by the median apparent period in each color-magnitude bin. Most of the WISE periodic variables are main-sequence stars. On the main sequence, the most common form of periodic stellar variability is eclipsing binaries. For our sample, on average, more massive sources tend to have longer periods. For close binaries, this makes sense because larger stars have a longer limiting period before reaching contact. A group of pulsating RR Lyrae stars can be seen at M G ∼ 1 and BP-RP∼ 0.7 with periods between ∼0.2 days and 1 day (Preston 1964;Kolenberg 2012;Das et al. 2018), significantly different from the surrounding EBs in the color-magnitude diagram. A few periodic sources are located below the main sequence. They may be white dwarf-brown dwarf binaries or cataclysmic variables. Figure 6 shows another view of the Gaia color-absolute magnitude diagram color-coded by the ratio of optical variability to infrared (IR) variability. To calculate the optical flux variability from Gaia, we follow the methods of Hwang et al. (2020) by computing the fractional variability from the Gaia photometric errors and then accounting for instrumental errors. We plot only those sources here that display appreciable variability Optical-to-IR Variability Figure 6. Gaia color-magnitude diagram color-coded by the ratio of optical variability from Gaia G band to infrared (IR) variability from WISE W1. A group of RR Lyrae can be seen in dark blue that vary more in the optical than in the infrared. The majority of the sources are eclipsing binaries and have equivalent optical and infrared variability amplitudes. The black line is the Pleiades main-sequence fit.
in Gaia. The IR r.m.s. fractional flux variability is ∼ 0.4 ln(10) times the standard deviation of the magnitude distribution. The black, dashed line is the Pleiades main-sequence spline fit following Hamer & Schlaufman (2019). Most of the sources have approximately equal optical and IR variability amplitudes, consistent with what is expected for eclipsing binaries. These sources are also offset to the bright side of the main-sequence fit indicating that there is excess flux due to the presence of a second star. A group of RR Lyrae can be seen that vary more in the optical than in the infrared due to the nature of their pulsation (Das et al. 2018).

Comparison with Chen et al. (2018) Catalog
Chen et al. (2018) present a catalog of 50,282 periodic variables from WISE identified with the Lomb-Scargle (Lomb 1976;Scargle 1982) periodogram and classified via light curve fitting. Our work differs from that of Chen et al. (2018) in that we extend the period search grid to shorter period objects, use a different periodogram, involve a Gaia DR2 cross-match, use data from a more recent NEOWISE release, and apply different quality cuts to the initial selection. Chen et al. (2018) also discarded sources with varying periods from AllWISE to NEOWISE and we make no such restriction.
We cross-match our periodic variables with those of Chen et al. (2018) using a matching distance of 1 arcsecond. Increasing the matching distance to 5 arcseconds does not significantly alter the results. Of all the sources with 6 ≤ var_flag ≤ 9 that we analyzed, 41,532 have a cross-match with the 50,282 periodic variables of Chen et al. (2018). Their remaining 8,750 variables were excluded by our initial quality cuts and were never a part of our periodogram analysis. 8,703 of the variables were excluded by our cuts on cc_flags to avoid contamination and confusion and the remaining 47 were excluded by our cuts on ext_flg to remove extended objects. Of the 41,532 variables in common, we mark 38,681 as periodic. We miss some of the sources due to our more stringent cuts on the quality of individual exposures, reducing the number of points in the light curve and preventing us from detecting them as periodic. We find compatible periods for 38,297 (99%) of the crossmatched sources. Sources are said to have compatible periods if they satisfy: The 384 outliers with non-compatible periods correspond to cases where either (a) our more stringent cuts significantly reduce the number of data points and leads to a different period estimate, (b) the detected period is greater than 5 days i.e. in a range poorly probed by WISE, or (c) we measure a period outside the range probed by Chen et al. (2018) indicating that one of the measurements is a harmonic. For ∼88% of the matches, our apparent periods, P app , are one-half the period cited by Chen et al. (2018), P chen18 , which is expected for a sample dominated by close eclipsing binaries where the primary and secondary eclipses are indistinguishable. We define the difference between the half periods ∆P = 0.5P chen18 − P app . Over 94% of the sources for which we recover half periods have | ∆P |< 0.001 days. Chen et al. (2018) classify 33,426 of the cross-matched sources as binaries, 780 as some type of Cepheid, and 2,464 as RR Lyrae. The remainder were assigned an ambiguous classification (Cepheid/Binary, Binary/RR Lyrae, or miscellaneous) due to a lack of characteristic features in their infrared light curves.
Finally, 24,845 of our periodic variable candidates were not in the periodic variable catalog of Chen et al. (2018) because of our expanded period search grid on the short-period end.
We present a catalog of 454,103 WISE variables with a variety of non-parametric measures and periodogram measurements. Non-parametric measures are discussed below in Section 4. We explicitly flag periodic variables and eclipsing binaries (the identification of eclipsing binaries is also discussed in Section 4). The data model is shown in Table 1.

Non-parametric measures of variability
In this section, we explore novel and computationallyinexpensive methods for light-curve analysis. Much information can be gleaned from just the distribution of the observed fluxes, without any analysis of their time dependence. If the cadence is suitable, the distribution of the observed fluxes can be assumed to be randomly drawn from (and therefore to be representative of) the underlying flux probability density function (PDF). The various moments of the observed flux distribution, assumed to be representative of the moments of the underlying PDF, are easy to measure and computationally fast. The first and the second moment -related to the mean magnitude and the r.m.s. variations around the mean -are already directly or indirectly incorporated into our analysis. In particular the second moment of the PDF is related to the var_flg we use as the primary step in our selection. We also compute a non-parametric amplitude (henceforth amplitude), which is the range between the 5th and 95th percentile of the magnitude values.
The third moment (skewness) of the PDF carries a wealth of information. For an object which stays at constant magnitude and occasionally undergoes dimmings due to occultations or eclipses, the PDF is expected to have a tail to fainter magnitudes (which correspond to the fluxes during occultations). For an object that instead undergoes occasional eruptions, we expect a tail to brighter magnitudes. Thus, the 3rd moment of the PDF can help us distinguish eclipsing from eruptive types of variability.
In this paper we introduce a new non-parametric measure, fainter fraction (FF ), which is related to the third moment. Specifically, for each light curve, we calculate the halfway point between the 5th and 95th percentile of magnitude and then calculate the fraction of points with magnitudes greater (fainter) than this halfway point by at least their measurement uncertainty. If FF is greater than 0.5, this means that the object spends most of its time in the faint state, becoming brighter less than half of the time, and it is likely to be of an eruptive type. The use of this measure is demonstrated in Figure 7. We also report the standard skewness in the catalog.
Finally, we introduce another non-parametric measure to describe the timescale of variability. Due to the peculiar cadence of WISE observations -multiple observations within a few days' time, but then not again for about six months -WISE has sensitivity to variability on a 1 day time scale, as well as to long-term variations. We use R, the ratio of the average r.m.s. variability within any day-long stretch to the r.m.s. variability of the entire light curve. Given that WISE has only a limited number of ∼day-long visits to the same position on the sky, this method is computationally inexpensive. For sources with equal signal-to-noise (S/N) of variability, the ratio R is expected to decline from 1 to 0 as we go from objects that vary on day-long timescales to those with only month-or year-long variability. Sources with low S/N of variability have both their short-term and long-term variability in line with the photometric σ.

Identification of Eclipsing Binaries
We identify a few physically-motivated cuts on period and a few of our non-parametric measures that can reliably isolate eclipsing binaries from other types of variability. In Figure 8 we show the kernel density estimate of our periodic variable candidates in the space of fainter fraction and amplitude. The sources cluster in this space and different clusters are characterized by light curves with different morphology. Nearly all of the periodic variables have FF <0.5 indicating that they are occulting. There is an upper sequence at FF of ∼ 0.35 and amplitude between 0.2 and 0.7 that, based on visual inspection, seems to be dominated by near-contact and contact binaries. A lower sequence at FF of ∼ 0.15 and of similar amplitudes has detached binaries with more narrow eclipses. The offshoot from the upper sequence at amplitude ∼ 0.3 mag and FF ∼ 0.25 seems to contain RR Lyrae.
We cross-match our sample with the Gaia RR Lyrae catalog of Clementini et al. (2019) and find 2172 matches. Figure 9 shows the kernel density estimate    This plot shows that the RR Lyrae do indeed cluster in the area of the FF vs. amplitude plot that was suggested by Figure 8. The gray, dotted lines show our initial selection of RR Lyrae. Right: Our period measurements for these sources. As expected for RR Lyrae, the periods are nearly all between 0.2 days and 1 day.
of the cross-matched RR Lyrae in the fainter fractionamplitude space and the period distribution of the RR Lyrae. To identify candidate RR Lyrae among WISE periodic variables, we apply the cut shown by the gray, dotted lines. Specifically, we require that 0.19 < fainter fraction < 0.34 and 0.22 < amplitude < 0.38. These cuts select the entirety of the aforementioned offshoot from the upper sequence and also a portion of the upper sequence itself. To ensure that we are not excluding too many close EBs, we restrict the apparent periods to between 0.25 and 1 days. After applying this cut, both the upper and lower sequences disappear and the RR Lyrae offshoot becomes the dominant feature in the fainter fraction-amplitude space. We choose a lower period bound of 0.25 days because we do not expect many RR Lyrae with periods below 0.25 days −1 and we do not want to exclude a high number of close eclipsing binaries. This figure shows that these three variablesperiod, amplitude, and fainter fraction -provide a means of separating RR Lyrae from EBs without resorting to the color-magnitude diagram. All told, this selection labels 6,632 of the periodic variables as candidate RR Lyrae. Included in this are 1,868 out of the 2,172 (∼ 86%) Gaia RR Lyrae from the crossmatched catalog of Clementini et al. (2019). We exclude the remaining Gaia RR Lyrae from the EB sample as well. In addition, we also cross-match our periodic variables with the Gaia DR2 single-object-study Cepheid catalog (Clementini et al. 2019). As mentioned above, we expect few Cepheids in the sample. The cross-match reveals only 386 Gaia Cepheids and we remove these from the EB sample.
The above cuts remove 1,906 out of 2,464 (∼77%) of the RR Lyrae identified by Chen et al. (2018) while retaining 30,656 out of the 33,426 sources (∼92%) classified by Chen et al. (2018) as EBs. We add these excluded sources back in to our EB sample. We remove the remaining RR Lyrae and also the 780 Cepheids classified by Chen et al. (2018). We are left with a low contamination sample of 56,323 eclipsing binary candidates.

Main-Sequence Binaries
We next turn to the color-magnitude diagram to isolate main-sequence EBs. Of our initial sample of 56,323 EB candidates, 48,902 have a Gaia best-neighbor crossmatch. Applying the quality cuts of Section 3 leaves us with 35,684 sources. To select main-sequence binaries, we start with the Hamer & Schlaufman (2019) Pleiades spline fit shown in Figure 6. Binaries are expected to be found above the main sequence for a wide range of mass ratios (Hurley & Tout 1998). We require that the source have an absolute G-band magnitude above the Pleiades spline fit and within 1.5 magnitudes of the spline fit value for a given color. This leaves 23,079 sources. Figure 10 shows the distribution of apparent period (P app = 0.5P orbital ) versus Gaia BP-RP color for our main-sequence EBs detected with WISE. The distribution is smooth across the entire period range and includes periods not probed by Chen et al. (2018). In blue is a kernel density estimate and the black line represents the theoretical minimum possible apparent period for an equal-mass, contact, main sequence eclipsing binary of a fixed age (1 Gyr) and solar metalicity as a function of mass and color. The calculation of this line Kernel Density Estimate Figure 10. Gaussian kernel density estimate in the space of apparent periods (0.5P orbital ) and Gaia BP-RP color for MS binaries. The black line overlaid with the red-yellow points represent the theoretical minimum apparent period (0.5P orbital ) for a contact, equal-mass, main-sequence binary.
follows Hwang et al. (2020). Briefly, from the PARSEC isochrone (Bressan et al. 2012) for an age of 1 Gyr and solar metalicity we get the stellar radius of an undistorted star (i.e. the Roche lobe volume radius, R L ). Following Eggleton (1983), we use the relationship between R L and the semi-major binary axis, a, to solve for a (R L = 0.38a). Finally, we set the masses of the two constituent stars, M 1 and M 2 , equal to each other and solve for the minimum possible apparent period using where G is the gravitational constant.
Many of our EBs are near this theoretical lower bound indicating that the majority of these systems are contact or near-contact binaries. Interestingly, some of the points are below the black, minimum-period line. These are systems for which the input assumptions break down. They can be systems in which the component stars are not of equal mass, not of identical color, younger than 1 Gyr, or of lower-than-solar metallicity. Some RR Lyrae are located close to the main sequence and may be present on the blue side of the plot (G BP − G RP < 0.6) accounting for the spillover in the bottom left. The dearth of sources hugging the line in the upper-left of the plot does not appear be due to decreased sensitivity to those periods (see subsection 2.5). It is possible that there are physical reasons for the paucity of contact binaries in the corresponding mass range -for example, if the magnetic braking mechanism is responsible for creating contact binaries (Hwang & Zakamska 2019), then it is expected to be inefficient at these masses due to lack of convection at M > 1.3M (Matt et al. 2011).
The peak in the period distribution of our mainsequence EBs is located about at 2P app ≈ 0.34 days. This is higher than the maximum of the period distribution of ∼0.27 days found by Rucinski (2007), but this sample is not volume-limited.

Shortest-Period Objects
There are 146 sources in our eclipsing binary candidates that have apparent periods below 0.1 days. We expect higher contamination and lower period accuracy in this period range due to the limitations of the WISE cadence. We visually inspect the phase-folded light curves and find nine are the result of bad data or an erroneous period measurement and that the remaining 137 represent real periodic signals. Figure 11 shows some example phase-folded light curves for these objects. These sources are prime candidates for comparison with other short-period binary catalogs (e.g Norton et al. 2011;Drake et al. 2014b) and subsequent analysis with the tools of Conroy et al. (2020). We cross-match these sources with SIMBAD with a matching distance of 1 arcsecond. After removing the bad data, we find 42 matches of which 37 were classified in SIMBAD as a specific type of variability. Of these 37, we find that four are pulsating, one is a long-period variable (that also exhibits short-period variability), and one is a Carbon star. The remaining 31 are all classified as some sort of binary. Five of these sources are cataclysmic variables (V* BL Hyi; SDSS J121209.30+013627.7; V* V347 Pav; RX J2218.5+1925; V* LX Ser). We find compatible periods for all of these cataclysmic variables indicating that our period measurements are accurate even in this short-period regime (Ramsay et al. 2004;Schmidt et al. 2005;Thorstensen & Halpern 2009;Avvakumova et al. 2013;Drake et al. 2014a). In addition, five sources are ellipsoidal variables indicating that the periodicity is due to the gravitational distortion of the stars (Morris 1985).

Young Stellar Objects and Extra-Galactic Variables
We compare R values and amplitudes for blazars and young stellar objects (YSOs). We cross-match all vari- In Figure 12 we show the results of this comparison. The distributions for R are largely similar with peaks around ∼ 0.2-0.3 indicating that both YSOs and blazars tend to vary on longer timescales. For amplitude, the YSOs are peaked below 2 magnitude while the blazar distribution is more uniform.

Separation of Periodic and Aperiodic Variables with Non-Parametric Features
Thus far, we have used our non-parametric measures for the classification of periodic variables identified with a periodogram. In this section, we discuss further applications of these measures. In particular, we show that the non-parametric features contain the requisite information to identify periodic sources. Figure 13 shows the distribution of periodic and aperiodic sources for a variety of the non-parametric measures. We additionally show the distribution of W1 signal-to-noise (S/N) ratio, since the quality of the data can influence the determination of periodicity. For many of these measures (especially the ratio of short-term to long-term variability), there is a marked difference between the distributions of periodic and aperiodic variables.
We investigate the potential of discriminating periodic variables from aperiodic variables using non-parametric features alone. Using empirical testing, we select the five most informative features -the W1 magnitude, W1 signal-to-noise ratio, estimated variability amplitude, fainter fraction, and the ratio between long-term and short-term variability. The variability amplitude and standard deviation have similar distributions (Figure 13), but we find that the amplitude is a marginally better discriminator. We apply the quality cuts of Section 2.7 on the W1 magnitude to exclude saturated and very faint sources and to ensure that we are not classifying trivially based on this data quality metric.
We define a two-class classification problem where the 'positive' class is defined as periodic variability, and the 'negative' class corresponds to aperiodic variability. We consider two classification models -the logistic regression and the random forest. In both cases, the general framework is the same -the model coefficients are solved for by 'training' on stars with known classes. The models can then perform predictions on new data, returning a probabilityp that a given star exhibits periodic variability. The logistic regression model assumes a linear relation between the features and the log-odds of a star exhibiting periodic variability (Yu et al. 2011). The random forest assumes no parametric model, instead relying on logical decision trees to map the feature space to the true and false classes (Breiman 2001).
To evaluate our methods, we randomly select 90% of the stars in our sample and train the models on them -we then make predictions on the remaining validation stars (10%) that were left out of training. Repeating our experiment for different random selections (or even different ratios of train/test sample size) does not affect our results. We find that we can correctly identify stars as periodic using only the non-parametric features with 95% accuracy. The random forest model outperforms the logistic regression, likely because a linear model is insufficient to describe the relationships between our features. A key metric for our particular use-case is precision -the probability that a star classified as periodic is genuinely periodic. A higher precision means  Figure 13. Distribution of periodic and aperiodic variables for a variety of non-parametric measures. The peak at 8th magnitude in the W1 magnitude plot is caused by saturation, rather than by physical variability, which is why there is not a corresponding peak in the distribution of periodic variables (Nikutta et al. 2014). We remove these sources as well as very faint sources with the quality cut of Section 2.7. The fainter fraction plot shows that our sample is dominated by periodic variables of the occulting type. The histograms are normalized to have unit area, to account for the fact that our sample contains more aperiodic variables than periodic variables. a lower false positive rate, preventing wasted followup resources. Additionally, the precision better reflects the performance of the model in this application due to the imbalance between classes -there are six times more aperiodic variables than periodic variables. We find that the random forest model trained on the five non-parametric features achieves an average precision of 90% on our test data, significantly better than the logistic regression's 70%. The confusion matrix of the classifier is balanced, i.e. the number of false positives and false negatives is similar. This indicates that the classifier is not biased to arbitrarily favor one class over the other.
An interesting application of these tools is to use the probabilities returned by the classifier to rank interesting candidates. The probabilities themselves can be incorporated into a hierarchical search that uses other prior information to inform the selection. In future large-scale surveys, such a pre-selection will be essential to efficiently allocate computational resources by running the full periodogram analysis on high-confidence periodic candidates first.
The brief demonstration of this section is intended mainly as a proof-of-concept of the information content of these features and has some important caveats. The first caveat is that one of the features (S/N ratio) measures data quality while others deal with real information about the source from the flux PDF and timescales of variability. The discrimination based on data quality is trivial; without good data it is impossible to flag a source as periodic. Therefore, it is important to ensure that the classifier is not biased to predict periodicity when the data quality is high. First, we verify that the S/N ratio is not the most important feature for the classification; by testing different combinations of features we find that several of other measures are more informative, and the majority of the discriminating power comes from the physical metrics and not just the data quality. Including the S/N ratio does improve the classification accuracy by ∼ 5%. We interpret this as the S/N ratio breaking a degeneracy in the variability measures between noisy sources and truly variable sources -objects with low S/N will automatically display variability due to noise, so including the S/N ratio as an explicit feature enables the classifier to incorporate this information. As another check, we run the classification on only high S/N sources, and yield nearly identical performance. Therefore, the classification appears to rely on meaningful information about the source, and does not seem to be biased by data quality.
The second caveat is that this demonstration is limited by the nature of the training sample. For example, due to the WISE cadence, long-period variables are almost completely absent from our training sample and so we do not expect the classifier to be adept at identifying them. Also, the use of periodogram-based classification as the 'ground-truth' for training biases the classifier to be more likely to detect similar types of periodicity as the periodogram, albeit at a much lower computational cost.
That said, the non-parametric approach has some important advantages over the periodogram and can even yield valuable information without reference to a periodogram-based analysis. PDF measures like this are more robust to photometric errors than the periodogram, and are meaningful even in systems with changing periods. These systems are of great interest, but would not normally be detected by the periodogram because implicit in the periodogram approach is the assumption of a constant period.

CONCLUSIONS
In this paper we present an analysis of ∼450,000 WISE variables. The variables are identified based on their r.m.s. variability in the AllWISE survey (Hoffman et al. 2012). We combine their AllWISE and NEOWISE lightcurves and conduct periodogram analysis to identify ∼63,500 periodic variable candidates. To this end, we search for periodicity over a finely spaced grid in frequency between 0.1 and 20 days −1 . For light curves that meet our quality cuts, we model the light curve as a series of periodic, orthogonal functions with five parameters and compute a variance statistic to quantify the quality of the fit. Good phase coverage and a significant value of the statistic at a frequency non-coincident with that of the WISE observing cadence are required to classify an object as periodic.
With our periodic variables in hand, we cross-match with Gaia DR2. We use the distribution of our periodic variables in the space of color-absolute magnitude to get a sense of our catalog contents. We also cross-match with a previous catalog of WISE periodic variables. We find that the sample is dominated by eclipsing binaries.
Next, we turn to less computationally expensive measures of variability. Collapsing the light curve into a flux distribution and measuring non-parametric estimates of its range and second and third moments yields a wealth of information. Together with our non-parametric measure of the timescale of variability (R, the ratio of variability on day timescales to that on month/year timescales), these measures can be used to cheaply and effectively (a) identify periodic variables and (b) classify the type of periodicity.
In terms of classification (b), we show that these interpretable and easily implemented measures provide an effective means of isolating eclipsing binaries from other types of periodic variability. We identify an all-sky sample of ∼56,000 eclipsing binaries in the infrared. Binary formation and evolution is interesting in its own right (Stepien 1995;Fabrycky & Tremaine 2007;Ivanova et al. 2013;Duchêne & Kraus 2013;Borkovits et al. 2016;Moe & Di Stefano 2017;Hwang et al. 2020) and binary systems have been linked to a wide range of astronomical phenomena (e.g. Kasliwal 2012;Belczynski et al. 2016). The majority of these binaries are contact or near-contact making them prime targets for future study.
In terms of identification of periodic variables (a), we demonstrate the high information content of the non-parametric features by using them to identify periodic variables at a much lower computational cost than the traditional, periodogram-based analysis. This type of analysis could be used to speed future studies of periodic variability and, in some cases, bypass the periodogram-based analysis entirely. Furthermore, the non-parametric method overcomes some of the shortcomings of periodograms. Importantly, because it does not implicitly assume a constant signal period as the periodogram does, it offers the possibility of identifying binaries exhibiting orbital evolution on human timescales.
Data Availability: The full catalog of WISE variables, periodic variables, and binaries is available as an electronic supplement to this paper 8 . The data model is listed in Table 1 ACKNOWLEDGMENTS This publication makes use of 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. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/ gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of the cross-match service provided by CDS, Strasbourg. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France. EP was supported in part by the Provost's Undergraduate Research Award at Johns Hopkins University. EP and HCH acknowledge support by the Space@Hopkins seed grant and by NASA ADAP grant. The authors acknowledge support from NASA ADAP grant 80NSSC19K0581. The MHAOV periodogram and the orthogonal, trigonometric function fit are implemented with the time series analysis package of Alex Schwarzenberg-Czerny that is based on a wrapper scheme by Ewald Zietsman. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018.