{'fraction_non_alphanumeric': 0.03286938302091518, 'fraction_numerical': 0.02574175926965241, 'mean_word_length': 4.484762979683973, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, '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}}
Periodic and Aperoidic Variability in the Molecular Cloud \rho Ophiuchus

Presented are the results of a near-IR photometric survey on 1678 stars in the direction of the \rho Ophiuchus (\rho Oph) star forming region using data from the 2MASS Calibration Database. For each target in this sample, up to 1584 individual J, H and K_s band photometric measurements with a cadence of ~1 day are obtained over 3 observing seasons spanning ~2.5 years; it is the most intensive survey of stars in this region to date. This survey identifies 101 variable stars with \Delta K_s band amplitudes from 0.044 to 2.31 mag and \Delta(J-K_s) color amplitudes ranging from 0.053 to 1.47 mag. Of the 72 young \rho Oph star cluster members included in this survey, 79% are variable; in addition, 22 variable stars are identified as candidate members. Based on the temporal behavior of the K_s time-series, the variability is distinguished as either periodic, long time-scale or irregular. This temporal behavior coupled with the behavior of stellar colors is used to assign a dominant variability mechanism. A new period-searching algorithm finds periodic signals in 32 variable stars with periods between 0.49 to 92 days.


INTRODUCTION
Photometric variability in young stars is an ubiquitous phenomenon. Physical mechanisms causing this variability include, but are not limited to, rotational modulation of starspots, evolution of the circumstellar environment, interstellar extinction, transit events and stellar pulsation. Large sample variability studies such as by Herbst et al. (1994) indicate these mechanisms often operate concurrently resulting in very complex photometric time-series in T Tauri stars. Interpretations are primarily limited because broad band photometry acquired by seeing-limited telescopes are unable to spatially resolve the inner regions around young stars, which are on the order of milliarcseconds for low mass stars in nearby star-forming regions.
High cadence, long temporal baseline photometric surveys, however, can temporally resolve the variations and help identify the dominant mechanism responsible for the variability. This type of survey work has the prospect of identifying peculiar observational phenomenon that may hold clues to unsolved problems in star and planet formation. One example is the physical interpretation fo AA Tau-like variability (Bouvier et al. 2003). This variability presents as discrete drops in flux from a nearly constant "continuum" flux level. The phenomenon is not unqiue to this star; Morales-Calderón et al. (2011) identified 38 stars exhibiting similar variability during an intensive photometric monitoring campaign of young stellar objects (YSOs) in the Orion Nebula Cluster. They find these dips have durations of 1 to a few days decreasing the flux by several tenths of magnitude in the IRAC bands and as much as 1.5 mag in the J band. These dips either appeared only once in 40 days (the survey's temporal baseline) or were periodic with periods between 2 to 14 days. The current physical interpretation for these observations is periodic obscuration of the star by a high latitude "warp" in the circumstelar disk (Bertout 2000;Bouvier et al. 2003). Longer temporal baseline photometric monitoring of additional YSOs will place better constraints on the periodic nature of these events which in turn will confirm/constrain the current physical interpretation.
Intensive photometric monitoring is most useful when it is conducted at a range of wavelengths. Near-IR observations probe the very inner regions from ∼0.01 to 1 AU for low mass stars (Dullemond & Monnier 2010). As young stars are typically rapidly rotating, this region includes both the co-rotation radius and dust sublimation radius. For a median disk lifetime of 5 Myr, this corresponds to millions of dynamical times at these orbital radii (Haisch et al. 2001). Therefore any indication of assymmetries in the circumstellar disk must be driven by some phenomenon (i.e. inclined magnetic dipole, planetary formation). Understanding the physics in this region is of particular interest to protoplanetary formation and migration, as well as how mass accretes onto the star.
ρ Ophiuchus (ρ Oph) makes an excellent laboratory to test the ability of high cadence, long temporal near-IR observations to distinquish between variability mechanisms in young stars. ρ Oph is a dense star-forming region containing a few hundred known young stellar objects (YSOs) with ages ranging from 0.3 to 3 Myr. The region is rich with variable stars; previous surveys having identified more than 100 photometrically variable stars (Greene & Young 1992;Barsony et al. 1997Barsony et al. , 2005Bontemps et al. 2001;Wilking et al. 2005;Alves de Oliveira & Casali 2008). Photometric surveys are limited to the near-to far-IR due to large amounts of visual extinction ranging from 5 to 25 A V in the cloud core (Cambrésy 1999). This complex interstellar environment could itself be responsible for detected photometric variability. Plavchan et al. (2008, hereafter P08) carried out a pilot study of 57 stars in the ρ Oph field using photometry collected by the Two Micron All-Sky Survey Calibration Point Source Working Database (2MASS Cal-PSWDB). That study identified periodic variability in two YSOs given a sample of candidate M stars. The study presented here expands on the initial pilot study performed in P08 and will include the full ρ Oph field data set from the 2MASS Cal-PSWDB to better understand the variability of young stars in this cloud.
In §2, the details concerning the observations and source selection for this survey are discussed. The variability analysis, along with a discussion on advantages to high cadence observations, is presented in §3. The methods used to find both periodic and long time-scale variability are found in §4. In §5, the characteristics of the variability catalog, as well as a discussion of variability mechanisms is presented. This section also discusses 6 stars where two variability mechanisms are estimated to operate concurrently. Finally §6 contains a summary of the findings reported by this work.

Observations
The Two Micron All-Sky Survey (Skrutskie et al. 2006, 2MASS) imaged nearly the entire sky via simultaneous drift scanning in three near-infrared bands (J H K s ) between 1997 and 2001. Observations were taken at the northern Mt. Hopkins Observatory and the southern CTIO facility. Photometric calibration for 2MASS required hourly observations of 35 calibration fields split evenly between the northern and southern hemispheres. Each calibration field is 1 • in length and 8.5 wide. One calibration field lies in the direction of the Ophiuchus constellation. This field is centered at α = 16 h 27 m 15.6 s and δ = -24 • 41 23 (J2000) and covers part of the ρ Oph L1688 cloud core (Bok 1956). These data have an observing cadence of ∼1 epoch per day. A complete observation is comprised of six consecutive 1.3 second scans in declination with a nearly constant right ascension. Each scan is offset by 5 in right ascension to minimize errors from pixel effects. The six scans, or "scan group", are finally co-added to minimize short time-scale and systematic variations. A complete scan group is obtained in approximately 8 minutes (Cutri et al. 2006, §III.2b). The maximum number of scans for a single star is 1584 divided by 6 or 264 scan groups.
Photometry is extracted from the calibration field via the 2MASS Point Source Catalog (2MASS PSC) automated processing system. Details of the system's implementation are described in Cutri et al. (2006); here a brief summary is given. Photometry for sources fainter than J = 9, H = 8.5 and K s = 8 mag, are extracted by profile-fitting. Profile-fitting compares the source flux to a pre-generated point spread function (PSF) via χ 2 minimization. The PSFs are selected from a look-up table with respect to a dimensionless seeing index that is updated regularly during each scan. The seeing index characterizes the atmospheric seeing during specific observations. The library of PSFs is generated by empirically fitting the 50 brightest stars in a single 2MASS calibration scan with a specific average seeing index. This scan is not necessarily of the ρ Oph field, but a calibration field containing a different slice of the sky. An error at the few-percent level may be present in the resulting photometry due to mismatched PSFs arising from rapid seeing variations.
For the few sources brighter than the above cut off magnitudes, photometry are extracted using a 4 fixed aperture corrected using a curve-of-growth. Atmospheric seeing conditions can place as much as 15% of the flux from a point source outside this fixed aperture. A curve-of-growth correction is a constant factor added to measured photometry to simulate measurements taken using an "infinite" aperture. The benefit of this method is avoiding decreased signal-to-noise and potential source confusion arising from large aperture photometry. However, curve-of-growth corrections assume the sources are unresolved single stars that can be approximated by a PSF. Therefore photometry for extended sources (i.e. stars embedded in bright nebular emission) or multiple systems are not properly characterized with this method. All the data scans are compiled in the 2MASS Cal-PSWDB.

Source Identification
The source selection in the ρ Oph field is similar to that described in P08, which is summarized here. A parent sample catalog of 7815 sources is constructed from a co-added deep image of the field (Cutri et al. 2006). For each target in the parent sample, the 2MASS Cal-PSWDB is searched for detections within a 2 matching radius. This radius is several σ larger than the 2MASS Point Source Catalog (PSC) astrometric precision and astrometric bias between the PSWDB and PSC (Zacharias et al. 2005;Skrutskie et al. 2006). This ensures confidence all Cal-PSWDB detections for the parent 7815 sources are found within the PSC astrometric precision.
Of the 7815 stars identified in the parent sample, 1678 stars have a sufficient number of detections for variability and periodic analysis. This sample of 1678 stars is henceforth referred to as the target sample. A 'sufficient number' is defined as stars detected in ≥10% of the observations in either J, H or K s and ≥ 50 detections in the J band. The first constraint ensures a sufficient number of data points for a robust periodogram computation. The 10% limit is an ad hoc limit chosen to reduce the noise present in the variability statistics. The second constraint removes sources near the FOV edges that are not present in most scans. Finally, despite the success of the 2MASS prescription to produce high-quality photometric measurements, occasionally photometry affected by latent image artifacts, spurious detections and poor quality detections still persist in the database. The reader is referred to P08 for a full treatment on how sources with poor photometry are characterized and excluded. Cutri et al. (2006) describes the different varieties of latent image artifacts arising from a number of phenomenon associated with the optical system. These artifacts are identified and removed via visual inspection. Multiple simultaneous detections found within the 2 search radius of a target, which are typically spurious byproducts of the source extraction pipeline, are eliminated. Simultaneous detections are when two (or more) detected sources are identified with a single source in the 2MASS Cal-PSWDB. Secondary detections are typically ∼0.5 -1.5 magnitudes fainter than the primary detection. In addition, they are typically detected in only one passband and only in one of the six scans. Unaccounted for spurious detections can give the appearance of variability and introduce systematic noise into any underlying periodic signals. Photometric measurements with poor spatial fits to the model PSF are also excluded from our analysis. A poor spatial fit occurs when the χ 2 value between the observed stellar profile and a model PSF is > 10. This is flagged as 'E' quality photometry within the Cal-PSWDB. Image saturation, cosmic rays, hot pixels, extended emission or partially resolved doubles could account for this poor quality fit to the photometry (Cutri et al. 2006). Photometry with poor spatial fits are systematically brighter by a few tenths of a magnitude, and this can falsely trigger the identification of variability. Table 1, available only on-line, lists the 1678 stars analyzed for variability. The magnitudes and errors listed are extracted from a co-added image of all calibration scans in that particular band.

Detection and Completeness Limits
For non-variable stars, the photometric measurement uncertainty is characterized by the standard deviation of all photometric measurements in a particular band. P08 showed that this photometric standard deviation as a function of apparent magnitude, for 2MASS photometry, follows the form of two distinct power laws. One power law describes brighter sources, where Poisson statistics dominate the uncertainty, while the second describes the dimmer sources, where the uncertainty is dominated by instrumental noise. The point of intersection between these two power laws, or "break point", designates the survey completeness limit where source detection drops below 100%. This power law model is used to predict the photometric scatter for a star, and any star that has a dispersion significantly (> 5σ) above this is identified as a candidate variable. The model, as a function of apparent magnitude m, is given by the following expression: where a m,l , b m,l , σ a,m,l , and σ b,m,l represent the slope, intercept and respective errors for each fit in each band over magnitude region l. This model is first applied to our sample of 1678 stars using coefficients derived by P08 from the entire 2MASS Calibration Field data set. These coefficients, however, yield a relatively poor fit to the ρ Oph calibration field. The lower noise in the ρ Oph data is attributed to better average seeing conditions during these observations. As a result, the model is re-fit on the ρ Oph data set alone to derive a new set of coefficients. The new coefficients with errors are listed in Table 2. Fig 1 shows the best-fit model along with the observed photometric scatter in each band.
The model yields completeness limits for this survey of 16.63, 15.75 and 15.10 mag in J, H and K s , respectively. These are significantly fainter limits than the 2MASS PSC as a whole, which are 15.8, 15.1 and 14.3 mag in J, H and K s , respectively. The approximate detection limits for this study, found by averaging the apparent magnitudes for the 10 faintest objects meeting our detection criteria, are 17.7, 16.7 and 16.0 mag in J, H and K s respectively.

Selection Criteria for Variability
Numerous surveys have used time-series analysis on multi-wavelength photometry to characterize young star variability (Mathieu et al. 1997;Carpenter et al. 2001Carpenter et al. , 2002Grankin et al. 2007Grankin et al. , 2008Morales-Calderón et al. 2011;Findeisen et al. 2013;Wolk et al. 2013, and   b Range in apparent magnitude; apparent magnitudes of < 6 are excluded from this model. references therein). The methods for identifing stellar variability are nearly as numerous as the variability studies themselves. These include, but are not limited to, the Stetson index, excess photometric dispersion, χ 2 statistic, cross-correlation and Fourier analysis (Stetson 1996;Carpenter et al. 2001;Barsony et al. 2005; Alves de Oliveira & Casali 2008).
Variable stars are identified in this work through 3 complementary methods that are sensitive to different types of variability. A full description of these techniques are presented in P08. The same terminology used in P08 is adopted in this work. Here a summary is presented along with specifics regarding this sample. The first and second methods, "flickering" and "excursive", identify variability in each band individually. The third method uses the Stetson index to identify correlated variability between bands.

Flickering Variability
Flickering variability describes when the star's photometric scatter significantly differs from the predicted scatter. Flickering variability is sensitive to continuous variability, as consistent, substantial variations are needed to significantly increase the observed photometric dispersion above the expected non-variable value. To identify flickering variables, an observed dispersion is calculated for all scan measurements of a star prior to conbining as a scan group. This is then compared to the star's expected dispersion with associated uncertainty, σ, calculated using the noise model described in § 2.2.1 (Eqn 1 and Fig 1). If the observed dispersion exceeds the expected dispersion by more than 5σ, the star is a candidate variable. This search is done separately for each of the 3 bands (J,H,K s ); a star can thus be flagged as a flickering variable in 1, 2 or all 3 bands. Following this criterion, 17 stars flag in only a single band, 23 flag in two bands and 54 flag in all three bands. If variability is intrinsic to the star, the expectation is the flickering will occur in more than one band. Low signal-to-noise photometry (the number of observations, N obs , is typically less than 500) might likely account for the 9 candidate variables that flicker in the K band only. It might also account for the 11 candidate variables that flicker in both the H and K bands. However there is no obvious explanation for the 17 candidate variable stars that flicker only in the J or H, J and H, or J and K bands. The average dispersions for these variable stars in J, H and K s are 0.12 ± 0.46, 0.12 ± 0.43 and 0.11 ± 0.35 mag, respectively. The listed errors are the standard deviation of the average dispersion. These values represent the dispersion intrinsic to the source, or specifically the dispersion after the predicted non-variable measurement dispersion is subtracted in quadrature from the observed dispersion.

Excursive Variability
Excursive variability describes when the average magnitude of a individual scan group is significantly deviant from the mean of all the star's scan groups. Excursive variability is sensitive to short time-scale variations such as a single eclipse event or flare. Excursive candidate variables are identified if the average magnitude for a single scan group exceeds the global mean by more than 5σ, where here σ is the co-added uncertainty in the scan group photometry. As with flickering variability, this search is done separately for each of the 3 bands. From the final variable catalog, 21 stars flag in only a single band, 19 flag in two bands and 41 flag in all three bands. Low signal-to-noise photometry (N obs typically less than 500) might account for the 10 candidate variables that are excursive in the K band only. It might also account for the 12 candidate variables that are excursive in both the H and K bands. However there is no obvious explanation for the 19 candidate variable stars that are excursive only in the J or H, J and H, or J and K bands. The average number of deviant scan groups per star in our variable catalog is 24, 42 and 57 in J, H and K s bands, respectively.

Stetson Index
The Stetson index describes the correlation in a star's photometric variation between different bands. The Stetson index is sensitive to variability whose amplitude is not significantly different between photometric bands. For example, the Stetson index is not sensitive to a strong increase in the Brγ emission line strength that might only affect 1 band. This index has been previously used on other molecular cloud 2MASS variability surveys in Orion A and Chameleon I (Carpenter et al. 2001(Carpenter et al. , 2002. The Stetson index is computed for all 1678 stars; a star is considered a candidate variable if this index is > 0.2. P08 determined this criterion based on 18 of 23 periodic variables, in that work, having indices above this value. The same index is adopted here since the observing methodology is identical in both works. This index is smaller than those adopted for the Orion A (0.55) and Chameleon I (1.00) surveys. The Orion A survey contained 29 epochs over a 36 day temporal baseline and the Chameleon I survey contained 15 epochs over 5 months. The smaller number of observed epochs in each case causes these surveys to be less sensitive to variability and thus in need of a higher index. A Stetson index of zero indicates random noise or no correlation between the photometry in different bands. A positive index indicates correlation between the photometry in two bands. The higher the index, the greater the correlation between the photometry. Using the Stetson index 57 stars flag as variable.

Excluding Seeing Induced Variables
A common way in which a non-variable star is misidentified as variable is from photometric variations caused by changing atmospheric seeing. Both photometric techniques used here (PSF fitting, fixed sums) are susceptible to this, especially in regions that are crowded or where there is bright nebular emission. Seeing estimates, corresponding to the average FWHM for each calibration scan, are provided for the Cal-PSWDB photometry. The typical seeing values range between 2.5" to 2.7" over the entire observing season (Cutri et al. 2006).
The possibility of changes in brightness being correlated with changes in the seeing are first investigated. This is done by computing the Pearson r-correlation statistic for each star, n. The statistic is given by the following: where m is the band, S m,t is the m-band seeing FWHM in arcseconds at epoch t and S is the average seeing in m-band. The separate quantities are summed over all N m,n mband observations for star, n. This statistic spans the range from -1 to 1 with negative values indicating inversely correlated variations and positive values corresponding to directly correlated variations. An inverse correlation means as the seeing worsens the star gets brighter. A direct correlation refers to the opposite effect. Since in Eqn 2, the photometry comparison (numerator) is computed in magnitudes and the photometric standard deviation (denominator) is computed first in flux units then converted to magnitudes, this can result in r values slightly outside the -1 to 1 range. A slight trend exists in the sample of 1678 stars toward an inverse seeing correlation in each band. The average r statistics in J, H and K s are -0.12, -0.11 and -0.05, respectively. Inverse correlation is likely caused by crowded fields where as the seeing worsens, flux from surrounding stars may encroach into the measured star's aperture or spatial profile. While these correlations are not very significant in most cases, it is noted the seeing in one band is slightly correlated with the seeing in another band. This is consistent with multi-band photometry taken simultaneously. To look for correlations between bands, the Pearson index for J and H are plotted in Fig 2. To characterize and flag seeing induced variability, a single seeing test is constructed to provide an estimate of seeing effects on measured photometry. Each correlation statistic (r J , r H , r K ) is considered a component of a single "seeing vector". This vector is rotated and transformed from cartesian to cylindrical coordinates so the z -axis corresponds to r J = r H = r K . This representation causes the seeing correlation to be axisymmetric about the z -axis, thus reducing the characterization of multi-band seeing correlation by one dimension. A "seeing ellipse" is described by where z n is the component of the seeing vector for star n, with standard deviation σ z , along the z -axis. ρ n is the component for star n along the ρ-axis, with the standard deviation σ ρ . Both σ z and σ ρ are determined from the distribution of the ensemble 1678 stars. A candidate variable is flagged as seeing correlated when the seeing vector length is larger than the seeing ellipse for the ensemble. This is the case when the left-hand side of Eqn 3 is greater than unity. This test excludes 19 candidate variables; the variability of these stars is likely solely caused by fluctuations in atmospheric seeing.

Final Variable Catalog
From the target sample of 1678 stars, 101 stars (6%) are identified as variable. These variable stars are referred to as the variable catalog. The variable catalog is listed in Table  3. The full set of light curves, color curves and color-color plots for all variables stars is only available on-line. For inclusion into the variable catalog, a star must not exhibit seeing correlated photometry (see § 3.1.4) and must meet 2 of the 7 variability criteria (see §3.1.1 -3.1.3). In addition, the 2 criteria must be met in different bands or in a single band along with the Stetson criterion. This last condition is imposed in order to prevent identifying variability due to poor quality or spurious photometry that is missed by the previous filters. The amplitudes of variability for stars within the variable catalog span a wide range. The range in ∆K s spans 0.04 to 2.31 mag and ∆(H -K s ) varies from 0.01 to 1.62 mag. The variable catalog contains 47 stars with ∆K s > 0.25 mag and 66 stars with ∆(H -K s ) > 0.1 mag.   Fig 3 contains the co-added calibration field in the direction of ρ Oph; the target sample of 1678 stars and the variable catalog of 101 stars are plotted to show their spatial distribution. It is clear that target stars are not evenly distributed in the field. A demarcation line at δ = -24 • 51' is set as an ad hoc determination of cloud membership. North of this limit is considered "on cloud" while anything south is classified as "off cloud". This demarcation corresponds roughly to where A V = 5 mag (Cambrésy 1999). Comparing the variability north and south of this demarcation, the "on-cloud" variable fraction increases to 15% while the variable fraction for the "field" drops to a mere 1%. This is consistent with the expectation young stars are are more often found spatially close to molecular clouds and are more variable than field stars.

Known Young Stars in the ρ Oph Field
Clues to the formation and evolution of young stars may be revealed by relating the variability to the stars' evolutionary states. As originally proposed by Lada (1987), young stars are classified into four evolutionary stages or classes (Class 0, Class I, Class II and Class III). Class assignment is typically based on photometry through the infrared slope index in the wavelength range from 2 to 25 µm. Class 0 stars represent cloud cores undergoing the initial stages of protostellar collapse. Class I stars are heavily embedded protostars with in-falling material from a circumstellar envelope forming an accretion disk. Class II stars are fully assembled stars with accretion primarily from the circumstellar disk channeled onto the star along magnetic field lines; classical T Tauri (CTTS) stars are another name for Class II stars. The last stage, Class III, represents stars yet to reach the main-sequence with depleted or no accretion disks due to mass accretion onto the star, photo-evaporation, or planet formation. They may nevertheless retain debris disks or disks with depleted inner holes. These stars are also known as weak-lined T Tauri stars (WTTS).
To identify if any of the 1678 sample stars have a previously assigned evolutionary class, the sample is cross-referenced with the ρ Oph L1688 cloud core mid-IR surveys by Bontemps et al. (2001, hereafter B01) and Gutermuth et al. (2009, hereafter G09). The measurements obtained by B01 were taken with the ISO ISOCAM LW2 and LW3 broad band cameras centered on 6.7 µm and 14.7 µm, respectively. The G09 survey obtained measurements in the Spitzer IRAC 3.6, 4.5, 5.8 and 8.0 µm bands complemented by J, H and K s 2MASS data for stars are used when there are no detections in either 5.8 or 8.0 µm. In addition, 24 µm Spitzer MIPS data is also used to verify YSO classifications in cases with high SNR (σ < 0.2 mag) and star luminosity ([24] < 7 mag).
The B01 survey provides YSO classifications for 54 of the 1678 target sample stars, while  (Cambrésy 1999). North of this demarcation contains higher visual extinction.
G09 provides classifications for 58 stars. However, overlapping targets between these surveys results in 40 stars classified by both B01 and G09, yielding YSO classifications for only 72 target sample stars. For 5 stars classified by both B01 and G09, the two surveys disagree on the classification. G09 classifies these 5 stars as belonging in an earlier evolutionary stage by one class than B01 (i.e. WL 22 is classified as Class I by G09 and a Class II by B01). In these cases, the classification by G09 is adopted because of the broader wavelength coverage utilized. Assuming the B01 survey identified all the young stellar objects in the ρ Oph region (425 YSOs), this survey contains ∼17% of these YSOs. Of the 72 stars with YSO classifications, 79% are identified as variable stars. As a function of YSO class, 92% of both Class I (12 of 13) and Class III (11 of 12) are variable stars. The variable fraction decreases to 72% (34 of 47) for Class II stars. The majority (14 stars) of the non-variable YSOs are Class II while ISO-Oph 99 is Class I. All of these stars are located "on cloud". As a YSO evolves in time the median brightness and color variability amplitudes decrease. The median peak-to-trough ∆K s amplitude for Class I, II, and III stars are 0.77, 0.31, and 0.08 mag, respectively. The median peak-to-trough ∆(H -K s ) color amplitudes are 0.81, 0.21, and 0.07 mag for each class respectively.

Advantages of High Cadence Variability Studies
In this section, the advantages of high cadence, long temporal baseline observations in variability studies are investigated. The results of this work are compared to the Alves de Oliveira & Casali (2008, hereafter AC08) survey of the ρ Oph central cloud core. The AC08 survey searched for variability in thousands of target stars within a ∼0.8 deg 2 field of view. These stars were observed in the H and K s bands during 14 epochs spanning May, June and July 2005 and 2006. The magnitudes of target stars fell within 11 to 19 mag in H and 10 to 18 mag in K s . This survey and AC08 have 464 stars in common. The prescription for identifying variables in AC08 is based on χ 2 fitting and cross-correlations between the H and K s photometry. Comparing the number of variables detected from the 464 stars, AC08 identifies 32 (7%) variables while this work identifies 82 (18%). The larger fraction of detected variables by this survey could be attributed to the higher sampling over a longer temporal baseline or from different sensitivities in the adopted variability criteria. To determine which explanation is more probable, histograms of the ∆K s peak-to-trough amplitudes for the variables identified by both this work and AC08 within the joint 464 star sample are computed. Fig  4 contains these histograms as well as the histograms for the ∆(H -K s ) peak-to-trough color amplitudes. It is clear from Fig 4 the fraction of variables with ∆K s < 0.5 mag detected by each survey is nearly identical. The same is true for variables with ∆(H -K s ) < 0.55 mag. Therefore, the higher fraction of variables detected, as compared to AC08, is most likely a consequence of the higher observing cadence. It is worth noting that 7 stars within the joint sample are identified as variable by AC08, but are not in this work. This work identified 5 of these stars as having photometry correlated with seeing. Therefore these stars may have been intrinsically variable within the observing window, however this variability could not be confidently confirmed.
While the detection fraction of low amplitude variables is nearly identical between surveys, the detection fraction of high amplitude variability stars is not. AC08 does not detect variables with ∆K s > 0.7 mag or ∆(H -K s ) > 0.55 mag. This work finds 5.25% of detected variables have ∆K s amplitudes greater than these upper limits. In addition, 6% of detected variables have ∆(H -K s ) color amplitudes greater than these upper limits. Within the 464 star joint sample, 25 stars are identified as variable in both this work and AC08. Strong correlations exist between the difference in amplitudes measured between surveys and the amplitudes measured in this work (see Fig 5). Sparsely sampled photometry will underestimate the amplitude of variability in both magnitude and color.
High cadence, long temporal baseline observations are vital for fully characterizing the variability of young stars. It increases the detection fraction of the survey allowing for more accurate statistics, such as the incidence of variable stars and distribution of variability amplitudes. In addition, this strategy is needed to sample the full amplitude of variability.

New Candidate ρ Oph Members
As photometric variability is an ubiquitous characteristic of young stars, it is a useful tool for assessing youth and potential membership in the ρ Oph star forming region. However, variability alone is not sufficient evidence for identifying potential members and additional constraints are needed, such as spatial location and location on a color-magnitude diagram. Candidate ρ Oph membership is first determined by cross-referencing the final variable catalog with previous surveys to identify previously known ρ Oph members (Strom et al. 1995;Barsony et al. 1997Barsony et al. , 2005Grosso et al. 2000;Ozawa et al. 2005;Wilking et al. 2005;Pillitteri et al. 2010). These are the same surveys used by AC08 to assign membership to their variable stars. This identifies 62 of the 101 variable stars as confirmed members of ρ Oph, which are plotted on a K s versus (H -K s ) color-magnitude diagram in Fig 6. For comparison, the 53 variable stars determined as ρ Oph members by AC08 are also plotted. Eleven stars are identified as ρ Oph in both surveys. A dashed line connects the data for these stars as observed by AC08 and this work.  The difference between amplitudes measured by AC08 and this work. Top: The comparison between measured K s variability. Bottom: The comparison between measured (H -K s ) color variability. In both cases, there is good agreement between the surveys for low amplitude variability. However as the amplitude increases, AC08 underestimates the variability. The dashed line in both plots indicates a difference of zero.
isochrone constructed using NextGen models for masses between 0.02 to 1.4 M at a distance of 129 pc (Baraffe et al. 1998). The distance is the weighted average between previous measurements (Loinard et al. 2008;Mamajek 2008). A star is classified as a new candidate member if it is located "on cloud" (see §3.1.5) and is brighter and redder than the 3 Myr isochrone (see Fig 6). Table 4 contains the 22 stars identified as candidate ρ Oph members from the previously unassociated 39 stars. Candidate member 2MASS J16270597-2428363 is classified as a Class II YSO thereby increasing the likelihood of membership. Follow up spectroscopic observations in the mid-IR for the remaining candidates to determine whether these stars are YSOs will provide additional evidence for membership.

TIME-SERIES ANALYSIS
Characterizing the amplitude, time-scale, and form (e.g. periodic vs. aperiodic) of variability provides valuable insights into the underlying physical mechanism(s) causing the brightness variations. Period-searching algorithms have been very helpful in this regard (e.g. (Lomb 1976;Scargle 1982)). In this section, two separate methods for measuring the time-scales of variability are discussed.

Periodicity Analysis via the Plavchan Algorithm
A novel period-searching algorithm, henceforth called the Plavchan algorithm (PA), is implemented to detect periodicity in identified variable stars. The algorithm described below is a more mature version than the one used in Plavchan et al. (2008). The version of the algorithm is used in the NASA Exoplanet Archive periodogram tool (von Braun et al. 2009;Ramirez et al. 2009). Tens of thousands of test periods are investigated by the PA algorithm with a uniform frequency sampling between 0.1 and 1000 days. For each trial period, P j , the PA starts by generating a phase-folded light curve from the time-series photometry. A phase is defined as the time (t i ) modulo the test period (P j ). This light curve is smoothed via boxcar smoothing with a phase width, p = 0.06. This smoothed light curve is designated as the prior, or reference curve. When the measured photometry for a periodic source is folded to the test period, the photometry is assumed to be approximately continuous and smoothly varying over the phased cycle. The difference between the measured photometry and the prior is computed for every photometric measurement, m i . This difference is compared to the difference between the measured photometry and a "non-variable" straight line, defined by the photometric mean (see Fig 7). A poor fit results when these two differences are equal or nearly equal to each other. A good fit results when the difference between the data and  (Bontemps et al. 2001;Gutermuth et al. 2009) Fig. 6.-Candidate ρ Oph membership. The filled blue circles indicate variables previously identified as ρ Oph members. The green filled circles mark the 22 new candidate ρ Oph members. The blue asterisks indicate the AC08 measured K s and (H -K s ) for variables previously identified as ρ Oph members. The dashed lines connect the same variable as it is detected in this work. The solid black line is a 3 Myr isochrone constructed using NextGen models for masses between 0.02 and 1.4 M for a distance of 129 pc (Baraffe et al. 1998). The red arrow corresponds to a reddening vector A V = 10 mag. the smoothed prior is smallest. This normalization removes the dependence on the absolute value and dispersion in m i . A quality of fit, χ 2 n 0 , is computed by Eqn 4 only over the 40 data points with the poorest fits (n 0 = 40) (i.e. the epochs with the largest difference between the data, m i , and the prior (average) in the denominator (numerator)): where the prior term, m priort , is the mean of m i if m i is within the boxcar smoothing window. The summations in the numerator and denominator in Eqn 4 are over independent sets of poorly fit measurements, since the poorest fit measurements by the prior might not be the same as the measurements that deviate the most from the mean. The best-fits periods have the largest χ 2 n 0 value. In other words, χ 2 n 0 represents the power of the periodic signal. The power indicates, for the PA, the relative improvements of the prior compared to a straight line for a given test period P j .
To evaluate the statistical significance of the power value for a peak period in the periodogram, or in other words to compute a false-alarm probability (FAP), there are several possible quantitative methodologies to arrive at an appropriate probability distribution. The approaches include -one, an analytic derivation from first principles; two, a monte carlo of periodograms generated by randomly swapping measurement values at each epoch; three, the distributions of power values at other periods in the same (adequately sampled) periodogram; and four, the distribution of maximum power values for all sources in an ensemble (mostly non-variable) survey. The first approach is rarely used in the literature, with the noted exception of the Lomb-Scargle periodogram (Scargle 1982). In the case of the Lomb-Scargle periodogram and typical radial velocity surveys, however, systematic errors in the velocity measurements can invalidate the assumptions in the first approach. The second monte carlo approach is often used as a more reliable method for Lomb-Scargle periodograms (Marcy & Butler 1998), and is equally applicable to the PA periodogram. In this section, the third method to evaluate a period's statistical significance is discussed. This third method is readily applicable to most time-series and is the method used in this work for computing the FAP for found periods. In the Appendix, the fourth method is discussed. The fourth method is survey dependent, but provides the insight that the PA periodogram is "well-behaved" with respect to changes in data values, number of observations, and algorithm parameters p and n 0 .
The distribution of power values in an adequately sampled PA periodogram for a nonvariable source is best described by a log-normal distribution. In this instance, adequately sampled means covering a broad dynamic range of periods and sampling the periodogram  BOTTOM: phased to a period of 3.5285 days. This is the most significant period from the periodogram. The dotted line indicates the mean magnitude for this star. The red lines in the middle and bottom panels is the prior generated for each period. Computing the χ 2 40 for the 3.5285 and 16.6672 day periods indicates the power value for the former is ∼9x larger, implying a much larger statistical significance.
at a large number of periods representative of the expected frequency resolution dictated by the cadence . Fig 8 contains the periodogram for the non-periodic star 2MASS J16265576-2508150. The power values vary about a mean value, or a "significance floor". The distribution is slightly asymmetric with a slight bias towards power values greater than the mean, consistent with a normal distribution in log-space. Fig 9 shows this distribution is very similar to the periodogram power value distribution for the boxcar least-squares (BLS) periodogram applied to the same source (Kovács et al. 2002), albeit with a different mean and standard deviation. The BLS periodogram traditionally makes the assumption of a normal distribution for evaluating the statistical significance of a peak period in the distribution of power values from an adequately sampled periodogram. However, again, a log-normal distribution is a more appropriate prescription for the BLS distribution (von Braun et al. 2009;Ramirez et al. 2009). While the assumption of a normal distribution of power values is probably adequate for both algorithms, a normal distribution will ascribe a greater statistical significance (i.e. a smaller FAP) to a peak period than a log-normal distribution. Therefore, the more conservative log-normal distribution is adopted in evaluating the statistical significance of peak periods in both the BLS and PA periodograms.
To determine if a period is statistically significant for a given source is this survey, the log of power values from the PA periodogram are computed, as well as the mean and standard deviation of the log-distribution. Power values that are 5σ outliers in the periodogram are identified as statistically significant periods with low FAP. Each of these significant periods are investigated via visual inspection of the photometry folded to the period in question. Finally the statistical significance of the derived period is confirmed by either the Lomb-Scargle or BLS algorithms, depending on the folded light curve shape. The Lomb-Scargle algorithm is optimized to identify sinusoidal-like periodic variations, while the BLS algorithm is better equipped in identifying eclipse-like periodic variations. Thus, the PA periodogram excels at identifying periodic signatures from both sinusoidal-like and eclipse-like time-series periodic variations (Plavchan et al. 2008). The period error is derived from the 1σ width of a Gaussian fit to the period's peak in the periodogram. In order to avoid confusion in the fit from other peaks, only periods within ±3% of the most significant peak are fit. An upper bound to confident periods is placed at 200 days. Stars are rejected as truly periodic with larger periods since the star will complete at most 3 cycles within the observing baseline. These "periods" are reported as timescales and described in §5.2.1.
From the 101 variables, 32 stars (32%) are identified to exhibit periodic variability with periods ranging from 0.49 to 92 days. Table 5  The PA found two statistically distinct (>20σ) periods for YLW 1C. The time-series folded to the shorter period (5.7792 days) exhibits a sinusoidal-like shape. The time-series folded to the longer period (5.9514 day) exhibits an "eclipse-like" shape where the star periodically dims from a near constant continuum flux. This prompted a search for secondary periods in the other 5 stars that exhibit eclipse-like periodic variability. We found 3 stars (YLW 1C, 2MASS J16272658-2425543, YLW 10C) to vary periodcally at two distinctly different periods; sinusoidal-like variability at one period and eclipse-like variability at the other. Initially the secondary period is not statistically significant; it is only discovered when the time-series of the eclipse event is removed. The PA is run only on the time-series preceding each eclipse ingress and after each eclipse egress. A small number (∼10) of sharp drops outside the eclipse events in the time-series for 2MASS J16272658-2425543 and YLW 10C are also omitted from the PA analysis. Errors in the secondary periods are determined in the same manner as the primary periods.
Since multiple variability mechanisms may be common in variable stars (Herbst et al. 1994;Morales-Calderón et al. 2011), we attempted to search for periodic variability in stars where the variability was complex. For 6 variable stars, the stellar brightess fluxuates about a mean level for one or two consecutive years. During the remaining time, a large amplitude variation is observed lasting longer than 50 days. The PA is run on the nearly constant timeseries omitting the large amplitude variation event. In 2 stars (WL 20W and ISO-Oph 126), the PA found a significant period in the "whitened" time-series. The time-series folded to the appropriate period results in sinusoidal-like variability with an amplitude ∼50% smaller than the large amplitude variation. This larger amplitude variation effectively masked the smaller amplitude periodic signal. For each star, the periodic variability could not be recovered during the large amplitude variation. Fig 10 contains the K s light curves for WL 20W and ISO-Oph 126, as well as the K s light curves folded to the identified periods. For WL 20W, 93 out of 262 scan groups were removed before the PA analysis. For ISO-Oph 126, 149 out of 262 scan groups were removed. Fig 11 shows the periodograms for both stars using the full time-series and the whitened time-series.

Measuring Long Time-Scale Variability
The long temporal baseline of the photometric time-series allows for the analysis of variability on month and year time-scales which are time-scales not well explored for young stars. Long time-scale (>50 days) variability differs from periodic or irregular variability in that the mean flux value may not remain nearly constant from season to season. In addition, MJD The periodograms from the PA only on the photometry not affected by the large amplitude variation. The 2.1026 day period is only seen and significant in the lower periodogram for WL 20. The same is true for the 9.114 day period of ISO-Oph 126. Additionally, this period peak power value is nearly 3x more significant than any power value detected using the complete set of photometry.
the photometry in one season may systematically brighten or dim while remaining constant in the other two seasons. Examples of these two phenomenon in the time-series for WL 20W and ISO-Oph 126 are shown in Fig 10. The intention in this section is to measure the time-scale of the single largest amplitude aperiodic or irregular variation.
Two criteria are used to identify stars exhibiting long time-scale variability. The first criteria is the difference between the photometric mean magnitude from one season to either of the remaining two seasons must be greater than 3σ, where σ is the average photometric error of the data over the entire temporal baseline (see Fig 23, WL 6). The second criteria is that the slope in the photometry in at least one season must be greater than ±5 • . The quality of the line fit determining the slope is assessed by visual inspection. The motivation for the second criterion is illustrated by WL 14 (see Fig 21). An obvious decreasing trend in the photometry is seen in the third season, however the sharp flux drop in the second season causes the mean flux between the two seasons to not satisfy the first criterion. Of the 101 variables, 31 stars (31%) satisfy at least one of these criteria and are designated long time-scale variables (LTVs).
A differencing technique is employed to measure the time-scale over which a LTV changes from one extreme in flux to the other. Fig 12 provides a visual demonstration of this method. In the top panel of Fig 12, a gradual dimming over the entire data set is observed. This global trend is seen in the time-series of 68% of the LTVs. Two different types of variability are believed responsible for the global trend and the long time-scale variation. Removal of the global trends provides an unbiased analysis of the shorter time-scale variation in the time-series superimposed on these trends. The global trend is a sustained, but small amplitude effect superimposed over the time-series including the larger amplitude, long time-scale variation. LTVs with these global trends are split evenly with 50% dimming over time and 50% brightening. The amplitude of the global trends range from 7.5 to 330 millimag/year, with a median value of 26 millimag/year. The median value corresponds to a change in the stellar flux of ∼60 millimag over the temporal baseline.
An accurate time-scale measurement for the largest amplitude variation can be complicated by the presence of small time-scale variability. The middle panel of Fig 12 shows how the light curve is smoothed with a 50 day moving median filter. The length of 50 days is chosen by visual inspection of the smoothed light curves; this timescale suppresses the smaller amplitude, shorter time-scale variability while preserving the shape of the long time-scale variation. The time-scale for the long time-scale variation is set to be the time difference between when the LTV is at one extreme in flux (i.e. brightest state) to the opposite extreme (i.e. dimmest state). This time-scale is determined by subtracting the smoothed magnitude found at time i with the smoothed magnitude found at time j using the following: where M i,j is the ∆mag between time j and time i, m j * is the magnitude at time j, m i * is the magnitude at time i and N obs is the total number of observations. The time between the largest ∆mag is recorded as the time-scale. In many cases, the full time-scale of the variation cannot be measured due to the data sampling. The bottom panel of Fig 13 shows the quantity M i,j as a function of times between measurements i and j. The "landscape" shows multiple peaks each corresponding to various time-scales of variability. The highest peak is only considered as only the time-scale associated with the greatest change in magnitude (i.e. largest M i,j ) is sought. Either extreme flux state may fall within a gap in the photometry or outside the date range of observations. Therefore these time-scales should be treated as lower bounds. The variability time-scales range from 64 to 790 days. Not all LTVs display only one discrete long time-scale variation. ISO-Oph 119 clearly shows two distinct long time-scale variations. For ISO-Oph 119 and similar cases, only the time-scale for the largest amplitude variation is measured. Figs 27 to 30 contain the K s light curves for these LTVs.
Despite observations spanning ∼2.5 years, in most cases it is not possible to conclude whether or not long time-scale variability is periodic. However, 6 LTVs have photometry suggestive of periodic behavior based on visual inspection of the stars' folded light curves correpsonding to periods ranging from 207 to 589 days. The light curves are folded to the most significant period found by the PA. These candidate periodic stars are identified in the first column of Table 6. These sources are not included with the periodic variables as the found periods are greater than the 200 day confidence limit (see §4.1).

DISCUSSION
The observational goal of this study is to measure the amplitudes and timescales of stellar variability, particularly in young stars. This information, in turn, places constraints on the physical mechanisms responsible for the variability. Empirical methods based on correlations between observed magnitudes and color have been employed to characterize stellar variability of young stars (Carpenter et al. 2001(Carpenter et al. , 2002Alves de Oliveira & Casali 2008). These methods consider variability due to rotational modulation of hot or cool starspots, variable extinction, variable mass accretion and structure changes in the circumstellar environment. Cool starspots are believed to be caused by localized magnetic inhibition of convection energy transport. Hot starspots, on the other hand, result from either surface flaring or heating by mass accretion onto the surface along magnetic field lines. Extinction may occur from asymmetries in an accretion disk or even from isolated dense regions of the parent molecular cloud passing through the line of sight. Variable mass accretion rates can cause the star brightness to vary through the clearing of the inner circumstellar disk. In addition, variability may be caused by energy released as material in an accretion disk moves toward a star by viscous processes. Finally, these mechanisms are not mutually exclusive and are often seen to exist simultaneously (Herbst et al. 1994).
Each of the above variability mechanisms can be distinguished based on the temporal nature of the variability and correlations between color variability to stellar brightness. 1 The following set of qualitative observables are developed to classify the observed variability and to connect these variations to physical mechanisms.
• Long-lived cool starspots result in periodic variability with periods consistent with the rotational periods of young stars ( 14 days) (Rebull 2001). This variability is often sinusoidal in shape. At the temperature range of most YSOs, the near-IR wavelength regime samples the Rayleigh-Jeans tail of the stellar energy distribution where the contrast between the starspot and surrounding photosphere is small (e.g. Vrba et al. (1985)). Therefore, the (J -H ) and (H -K s ) colors should remain constant (within photometric errors) as the brightness varies.
• Variability by hot starspots can either result in periodic or irregular variability. Longlived hot starspots caused by accretion onto the stellar surface may result in periodic variability. However it should be noted that accretion induced hot spots may display aperiodic behavior due to a stochastic accretion rate. Variability caused by flares will be aperiodic and will have time-scales on the order of hours to days. As with cool starspots, the period of variability will be consistent with the rotational periods of young stars. In both cases the affected photosphere should be hotter than the surrounding surface resulting in the star becoming bluer as the star brightens (Rodonò & Cutispoto 1988;Panagi & Andrews 1995;Yu & Gan 2006).
• Variable extinction can result in either periodic or long time-scale variability. Variability caused by asymmetries in the inner circumstellar disk, if present, may be periodic with periods from days to weeks. Unlike variability caused by starspots, periodic variable extinction need not appear sinusoidal but present more likely as eclipse-like features. These eclipse-like features are sharp drops or "dips" in the stellar flux with a regularity dependent on the observing cadence. Variability caused by asymmetries in the outer circumstellar disk (> 1 AU) will not be periodic within the temporal baseline of this study due to long period of revolution around the host star. This variability and variable extinction from inter-cloud material can occur on long time-scales, however as the time-scale depends on the system geometry, there is no expectation as to its duration. Variable extinction causes the star to redden as the star dims.
• Variability caused by a variable accretion rate within the circumstellar disk is not expected to be periodic. The time-scale of variability does place constraints on the physics causing this rate change (e.g. disk viscosity, time variable magnetic field) (Armitage 1995;Mahdavi & Kenyon 1998;Lai 1999;Terquem & Papaloizou 2000;Carpenter et al. 2001). During times of lower accretion rates, the inner disk cools and the inner hole becomes larger. This, in turn, decreases the contribution of dust reradiation, particularly in the K s band, to the overall energy budget of the star and circumstellar disk system. Therefore while the total system flux drops, a larger percentage of emitted radiation is from the star causing the system to become bluer as the system dims. However, if the inner circumstellar disk edge is dominated by the dust sublimation temperature, a observationally similar effect will result. In this case, an increased accretion rate raises the star's effective temperature in turn increasing both the distance to the circumstellar disk inner rim and disk vertical height. The result would be that the system would become brighter as it reddens. Both physical scenarios produce a qualitatively identical result to the observed correlation between brightness and color.
In an attempt to identify the dominant variability mechanism, stars in the variable catalog are placed into subclasses based on the observed shape and time-scale of variability. These subclasses are: periodic, long time-scale and irregular. These classifications along with the above criterion identified the likely dominant variability mechanism for 53 of the 101 stars in the variable catalog. The type of variability associated with each star is listed in Table 3 and each sub-class is described in the following subsections. The periodic sub-class accounts for 32% of the variable catalog with the majority (88%) lying "on cloud". Long time-scale variables make up 31% of the variable catalog. All LTVs reside "on cloud". The irregular subclass contains the most members comprising 40% of the variable catalog. Only 68% of irregular variables lie "on cloud". These subclasses are rough descriptions and are by no means mutually exclusive. For instance, WL 20W and ISO-Oph 126 are placed into both the periodic and long time-scale subclasses.
These criteria do not always allow for the dominant variability mechanism to be identified. The main reasons preventing an estimate of the mechanism are: the time-scale/period or color correlation is contrary to the above diagnostics, no dominant amplitude variability is clearly evident, or the photometry in J and H is below the completeness limits in each band resulting in no useful color information. Mechanisms appended with a question mark in Table 3 either possess a marginal color correlation via visual inspection, or the diagnostics did not definitively differentiate between proposed mechanisms.

Periodic Variables
The PA identifies 32 of 101 stars (32%) within the variable catalog as periodic with periods ranging from 0.49 to 92.28 days. Table 4 contains the list of periodic variables. The light curves for certain subsets of periodic variables are very similar in form when phased to the identified period. This allows for periodic variables to be separated into two subcategories: sinusoidal-like and eclipse-like. Assignment to a particular sub-category is based upon visual inspection of the folded light curve in the band with the highest signal-to-noise.

Sinusoidal-like Periodic Variables
Figs 13 to 16 contain the K s light curves for sinusoidal-like periodic variables. This subcategory of sinusoidal-like periodic variables includes the most periodic variables (25 stars) with periods ranging from 0.49 to 25.55 days. The peak-to-trough ∆K s amplitudes range from 0.06 to 1.64 mag, with a median value of 0.29 mag. The peak-to-trough ∆(H -K s ) color amplitudes range from 0.01 to 0.63 mag, with a median value of 0.19 mag. Typically, the light curve folded to the most significant period shows only one sinusoidal cycle. However, four stars (ISO-Oph 100, WL 10, WL 13, YLW 13A) show what could be interpreted as a second cycle at half the frequency (i.e. double the period).
Probable variability mechanisms are identified for these stars by applying the criteria discussed in §5. Correlations are qualitatively examined between the ensemble K s photometry and stellar colors when folded to the star's period. Fig 17 illustrates examples of these interpretations for both cool and hot starspots. For 18 of the sinusoidal-like periodic variables (72%), the variability in the (J -H ) and (H -K s ) colors does not correlate with the K s variability. This favors rotational modulation by cool starspots as the dominant variability mechanism. For 3 sinusoidal-like periodic variables (12%), the (J -H ) and (H -K s ) colors become bluer as the star brightens. This is consistent with the behavior expected from rotational modulation by an accretion induced hot starspot. The (J -H ) and (H -K s ) colors become redder as the star dims for 1 (4%) sinusoidal-like periodic variable. This favors variable extinction as the dominant variability mechanism. Finally for the remaining 3 sinusoidal-like periodic variables (12%), no dominant variability mechanism could be  assigned using the adopted criteria.
All sinusoidal-like periodic variables except 3 (2MASS J1625744-2504017, 2MASS J16271836-2454537 and 2MASS J16272533-2506211) are located "on cloud". This sub-category contains 19 stars with a YSO classification: 3 Class I (25%), 8 Class II (24%) and 8 Class III (73%). 2 The variability mechanism for 2 of the Class I stars is cool starspots, while the mechanism could not be identified for the third. Of the Class II stars, 5 vary due to cool starspots, 1 from an accretion-induced hot starspot and 1 is unknown. The variability of all Class III stars is caused by cool starspots.
Two key points can be made by analyzing the 18 sinusoidal-like periodic variables where cool starspots is the believed variability mechanism. First, cool starspots on young stars persist on preferential longitudes on year timescales. This is evidenced by a lack of phase drift in the folded light curves. This phenomenon of preferential or active longitudes have been associated with a number of chromospherically active stars (e.g. RS CVns, FK Com) (Strassmeier et al. 1988;Zeilik et al. 1988;Henry et al. 1995;Jetsu 1996). Second the variability amplitude due to cool starspots changes on much shorter timescales as evidenced by the significant scatter within the folded light curves. This amplitude variability is likely caused by an evolving starspot covering factor and/or starspot temperature. The covering factor is defined as the area of the observed stellar disk covered by the starspot(s).
ISO-Oph 96, ISO-Oph 133, ISO-Oph 149 and 2MASS J16272533-2506211 differ from the remaining sinusoidal-like periodic variables as the light curves for these 4 stars are asymmetric (i.e. they have a sharp increase in flux then decrease more slowly). The (J -H ) and (H -K s ) color variability are not correlated to the K s variability for the first three asymmetric sinusoidal-like periodic variables. This favors a dominant variability mechanism of rotational modulation by cool starspots. Asymmetric light curves have been observed for both WTTS and chromospherically active dwarf stars (Cutispoto et al. 2001(Cutispoto et al. , 2003Grankin et al. 2008;Frasca et al. 2009). In both cases, the variability is believed to be caused by magnetically generated cool starspots. "Reverse" asymmetric light curves with a slow rise in source flux followed by a steep drop are also observed. Frasca et al. (2009) is able to closely model the RIJH asymmetric light curves of the WTTS V1529 Ori by rotating a stellar surface with two cool starspots of unequal areas separated by ∼130 • in longitude. The size of the leading cool starspot determines if a "forward" or "reverse" asymmetric light curve is seen.
2MASS J16272533-2506211, hereafter designated 'J211', is peculiar due to its unique and difficult to interpret brightness and color variations. Fig 18 contains   The folded K s and color curves for the sinusoidal-like periodic variable ISO-Oph 135. The lack of near-IR color changes with brightness changes favors cool starspots as the variability mechanism. Right: The folded K s and color curves for the sinusoidal-like periodic variable WL 11. The stellar color becomes bluer as the K s photometry becomes brighter. This favors rotational modulation of accretion induced hot starspots as the dominant mechanism. The red line in each plot indicates the mean value.
(H -K s ) photometry for J211 folded to P = 0.485143 ± 0.000050 days. This is the shortest period star among the periodic variables. The peak-to-trough K s amplitude is 0.40 mag and the ∆(H -K s ) color amplitude is 0.30 mag. The (J -H ) color for J211 clearly becomes bluer for a phase duration of ∼0.3 (∼3.5 hrs) centered approximately on the times of maximum brightness. Somewhat surprisingly, however, no similar variation is seen in the (H -K s ) color during the same period. The data is deemed reliable as the J and H photometry are significantly brighter than the survey completeness limits. The variability mechanism is not identified for J211 as this K s -color behavior is inconsistent with any criteria discussed in §5. The shape of the light curve coupled with the short period suggests J211 might be a RR Lyrae variable. However the peculiar color behavior is not expected in these stars.

Eclipse-and Inverse Eclipse-like Periodic Variables
Eclipse-like periodic variables possess photometry containing sharp periodic drops in source flux. The duration of these drops, or possibly "eclipses", is in all cases less than a phase of 0.3 when the photometry is folded to the most significant period. Fig 19 contains the folded light curves for the 6 eclipse-like periodic variables. These eclipse-like periodic variables are assigned to this subclass by visual inspection; it is not possible to confidently state that these sharp changes in photometry are bona fide occultations of star light (i.e. true eclipses).
The periods for the eclipse-like periodic variables range from 2.95 to 8.00 days. The durations of these eclipses range from 5.8 to 12.7 hours. The ∆K s amplitudes range from 0.21 to 0.51, with a median value of 0.31 mag. The ∆(H -K s ) color amplitude range from 0.10 to 0.25 mag, with a median value of 0.11 mag. These amplitudes in both magnitude and color represent the total change in stellar flux and they do not necessarily represent eclipse depths since there is considerable scatter in the out-of-eclipse photometry. The eclipse depths and how they are determined are described below. The variability mechanism for the eclipse is determined in the same manner as with the sinusoidal-like periodic variables. Correlations between the K s photometry and stellar colors (J -H ) and (H -K s ) for the eclipse event are assessed visually and compared with the criteria discussed in §5. For 4 (66%) eclipse-like periodic variables, the (J -H ) and (H -K s ) colors become redder as the star dims. The color correlation coupled with the short, periodic behavior favor extinction, possibly by the inner region of a circumstellar disk, as the dominant variability mechanism causing the eclipse. The variability in the colors during the eclipse for ISO-Oph 102 are not correlated with the K s photometry consistent with variability caused by rotational modulation of cool starspots. Unfortunately, both the J and H photometry for YLW 10C are dimmer than the survey completeness limits and the lack of color information prevents confident identification of a variability mechanism. All 6 eclipse-like periodic variables are classified as YSO Class II (15%) 3 . All stars in this sub-category are located "on cloud" except 2MASS J16271513-2451388. Morales-Calderón et al. (2011, YSOVAR) qualitatively identified a number of similar eclipse-like variables in their mid-IR variability survey of YSOs within the Orion Nebula Cluster. The survey found 38 stars exhibiting brief, sharp drops in stellar flux. These stars are identified as AA Tau or "dipper" variables. The variability mechanism is believed to be due to high latitude warps in the inner accretion disk periodically occulting the star (Bertout 2000;Bouvier et al. 2003). The 4 eclipse-like variables that show evidence of extinction could be considered AA Tau variables, and possibly all 6 eclipse-like systems.
Under the assumption all 6 eclipse-like periodic variables are AA Tau variables, constraints on the spatial location within the circumstellar disk and size of the hypothetical occulter are investigated. Marsh et al. (2010) performed a deep mid-IR imaging survey of the ρ Oph 2MASS Calibration field used in this work. They computed the T ef f , A V and mass for 5 of the 6 eclipse-like periodic varibles. This was done by fitting model spectra to observed SEDs. The SEDs were computed from photometry in the J, H, K, [3.5] and [4.5] bands. The model spectra were obtained using the COND, DUSTY and NextGen models. T ef f and A V are found by minimizing the following equation: where α is a flux scaling factor, f obs λ and f obs λ (T ef f ) are the respective observed and model fluxes at wavelength λ, σ λ is the flux uncertainty, r λ is the absorption at wavelength λ relative to A V and β is a constant penalty parameter. In the Rayleigh-Jeans regime a degeneracy exists between T ef f and A V such that a high temperature star seen through low extinction will have a SED similar to that of a low temperature star seen through high extinction. The parameter β is used to break this degeneracy by penalizing solutions with low values of A V . This parameter was optimized by using the published photometry of a large sample of spectroscopically confirmed brown drawfs. Values of A V have errors between 1 to 2.7 while T ef f is accurate to within 860 K. The COND and DUSTY models then yield a model-unique mass for each star. They used a mass-temperature relationship to derive the mass for hotter stars fit by the NextGen models. They conclude an accuracy in the mass estimate to within For the case of ISO-Oph 106, a 1 Myr isochrone from Siess et al. (2000) is used to determine the T ef f and radius for ISO-Oph 106 by assuming a mass of 0.5 M . The stellar radius for the remaining eclipse-like periodic variables is computed by R = ( L T 4 ef f ) 1 2 , where each quantity is in solar units. The stellar luminosity, L , is computed using the prescription outlined in Natta et al. (2006) except in the case of YLW 10C. The prescription relates the stellar luminosity as a function of J magnitude and extinction A J . The extinction is computed using the (J-H) and (H-K) colors corrected into the CIT system using the ρ Oph extinction law by Kenyon et al. (1998) and the CTTS locus defined by Meyer et al. (1997). The J band photometry for YLW 10C is below the survey completeness limit. The luminosity for this star is found using the 1 Myr isochrone mentioned above and the mass determined by Marsh et al. (2010). The estimated mass, T ef f , A V , and R are in Table 5.
Assuming the occulter has negligible mass and orbits under Keplerian rotation, the occulter's distance from the host star can be computed for each eclipse-like periodic variable. The diameter of the occulter is computed from the duration of the eclipse event and its location in the circumstellar disk. Strictly speaking, the computed diameters are along the orbital path. No assumption concerning the occulter geometry (i.e. spherical, ellipsiod) is made. This diameter is considered a strict lower bound as the disk geometry, occulter impact parameter and occulter opacity are unknown. The eclipse depth, however, is determined by first removing the eclipse-event photometry from the time-series data (see §4.1). The eclipse depth is then determined to be the difference between the maximum magnitude in the eclipse feature relative to a median non-eclipse mean magnitude. This calculation results in occulter distances that range between 1.83 to 7.21 R , with a median value of 3.20 R . The occulter size ranges from 0.58 to 5.07 R , with a median value of 2.05 R . The range in ∆K s eclipse depth is 0.12 to 0.51 mag, with a median value of 0.27 mag. Table 5 summarizes the results of this investigation. Inverse eclipse-like periodic variables are similar to eclipse-like periodic variables, but the "eclipse" is an increase in source flux rather than a decrease. Fig 20 contains the folded K s and color curves for WL 4 and YLW 16A, the 2 inverse eclipse-like periodic variables in the variable catalog. WL 4 is a Class II YSO whose period of variability is 65.61 ± 0.40 days. The peak-to-trough ∆K s amplitude id 0.67 mag and the peak-to-trough ∆(H -K s ) color amplitude is 0.19 mag. The (J -H ) color for WL 4 becomes redder as the star brightens during the inverse eclipse event. However, the (J -H ) color change starts just prior and ends just after the inverse eclipse event. The difference in each case is ∼0.1 in phase or ∼6.6 days. The (H -K s ) color is not correlated with the K s variability.
The period of variability for YLW 16A is longer than WL 4 at 92.28 ± 0.84 days. The amplitudes of variability for YLW 16A, a Class I YSO, are also larger with peak-to-trough ∆K s and ∆(H -K s ) amplitudes of 0.95 and 0.34 mag, respectively. As the J band photometry is dimmer than the survey completeness limits, no reliable (J -H ) color information is available. The (H -K s ) color variability is sinusoidal-like, but is not aligned with the K s variability. Both stars reside "on cloud".
The variability mechanism for both WL 4 and YLW 16A is believed to be related, and similar to the interpretations proposed in separate letters (Plavchan et al. 2008(Plavchan et al. , 2013. Here the proposed variability mechanism is summarized. Both systems contain a visual binary companion detected through high resolution direct imaging (Ratzka et al. 2005;Plavchan et al. 2013). The two visible components for WL 4 are separated by 0.176 and separated by 0.3 in the case of YLW 16A. This corresponds to a projected linear separations of 23 AU and 39 AU, respectively, given a mean distance of 129 pc (see §3.4). In each system, the large amplitude variability is believed to be intrinsic to one of the visible pair. This component is, in turn, hypothesized to be a close binary surrounded by a circumbinary disk; this system is thus a triple system. The influence of the wide companion has caused the plane of the circumbinary disk to be inclined to the orbital plane of the inner binary. The variability results when each component of the inner binary is periodically obscured by the circumbinary disk as the binary orbits around the barycenter. Kusakabe et al. (2005) proposed a similar model to explain the variability for KH-15D.

Long Time-Scale Variables
The largest amplitude variability in long time-scale variables is not observed to be periodic, but show consistent trends, unlike irregular variables. The K s light curves are shown in Fig 21 to    In §4.2, a method for quantifying the time-scale of the extended brightness variations is detailed. These time-scales range from 64 to 790 days, the latter being near the full duration of the observing campaign. The peak-to-trough ∆K s amplitudes for all LTVs range from 0.05 to 2.31 mag, with a median value of 0.29 mag. The peak-to-trough ∆(H -K s ) color amplitudes range from 0.06 to 1.32 mag, with a median value of 0.23 mag. The two most probable mechanisms for these aperiodic variations with time-scales much longer than typical stellar rotation periods are variable extinction and variable mass accretion rates. Fig  25 shows two examples of the change in K s brightness and stellar color caused by these two mechanisms. Extinction causes the star to become redder as the star dims. Changes in the mass accretion rate cause stars to become bluer as the star dims. For 12 LTVs (39%), the (J -H ) and (H -K s ) colors become redder as the star dims favoring variable extinction as the dominant variability mechanism. The (J -H ) and (H -K s ) colors become bluer as the star dims in 11 LTVs (34%) favoring variable mass accretion rates as the dominant variability mechanism. The remaining 7 LTVs (23%), either do not have useful color information because the J or H (or both) photometry is below the survey completeness limits, or the brightness-color correlation does not agree with any of the 4 listed variability criteria. No dominant variability mechanism is assigned to these stars.
One intriguing scenario to explain why LTVs do not seem to favor one variability mechanism over another is the viewing angle. For LTVs where variable accretion is the favored mechanism, the system could be more face-on providing a clearer view of the inner disk hole. Variable extinction due to circumstellar disk asymmetries is more easily seen at higher disk inclinations where the "puffed up" outer disk attentuates the light from the inner disk. As these two mechanisms have opposing brightness-color correlations, systems with no measured correlations may represent intermediate viewing angles. In this case, the measured effects from variable accretion will "cancel" out or confuse the measured effects from variable extinction.
All LTVs are located "on cloud" and 25 stars of the 31 LTVs are classified as a YSO: 7 Class I (58%), 15 Class II (44%) and 3 Class III (27%). 4 The favored variability mechanism in the Class I LTVs is variable extinction for 2 stars, variable mass accretion for 2 stars and unidentified for 3 stars. The variability in the Class II LTVs is consistent with variable extinction in 5 stars, variable mass accretion in 8 stars and is not identified for 2 stars. The K s -color correlation in 2 Class III LTVs favors variable extinction as the dominant variability mechanism while the mechanism for variability in the third Class III LTV is not identified.  This is an example of variability caused by extinction. As the K s magnitude drops the colors become redder. Right: The same light and color curves for the long time-scale variable ISO-Oph 117. This is an example of variability caused by variable mass accretion. As the K s magnitude drops the colors become bluer.
Based upon visual inspection of folded light curves, 6 LTVs are considered candidate periodic variables. These candidate periodic LTVs are denoted in Table 6 and Fig 26 contains their folded K s light curves. The variability time-scales, ranging from 207 to 589 days, for the candidate periodic LTVs are measured using the Lomb-Scargle algorithm. The PA does not find the time-scales found by Lomb-Scargle to be significant. The stars, on average, have higher flux and color amplitude variability than the LTVs taken as a whole. The peak-totrough ∆K s amplitude for these candidate periodic variables range from 0.16 to 1.13 mag, with a median value of 0.77 mag. The peak-to-trough color amplitude range from 0.23 to 1.12 mag, with a median value of 0.82 mag. For 3 candidate periodic LTVs, the (J -H ) and (H -K s ) colors become bluer as the star dims favoring a variable mass accretion rate as the dominant variability mechanism. The (H -K s ) color of ISO-Oph 112 reddens as the star dims. This is consistent with variable extinction as the dominant variability mechanism. A combination of J band photometry below the survey completeness limits and the K s -(H -K s ) color correlation not matching any of the 4 criteria precludes the identification of the dominant variability mechanism for 2 candidate periodic LTVs. All the candidate periodic LTVs are classified as a YSO: 1 Class I star and 5 Class II stars.

Irregular Variables
The variable catalog contains 40 stars (40%) that are clearly variable, but the largest amplitude variability is not periodic or coherent on long time-scales. The K s light curves are located in Figs 27 to 31. Table 8 contains the list of irregular variables.
The ∆K s amplitude range from 0.04 to 1.11 mag, with a median value of 0.14 mag. The ∆(H -K s ) color amplitude range from 0.05 to 0.75 mag with a median value of 0.14 mag. Using the variability criteria, discussed in §5, the primary variability mechanism is only identified for 4 irregular variables (2MASSJ16265861-2446029, ISO-Oph 116, YLW 13B, ISO-Oph 87). The first three 3 exhibit a long time-scale variation in at least one observing season where the star becomes bluer as it dims. This is indicative of variable mass accretion as the variability mechanism. These stars are not considered LTVs as this variability is not the largest amplitude variability in the time-series. The variable accretion for YLW 13B is identified to occur for ∼ 115 days in the second year with a ∆K s ∼0.15 mag and ∆(H -K s ) 0.11 mag (see Fig 30). ISO-Oph 116 varies via variable accretion at least twice (see Fig  28). The first time occurs for ∼170 days in the first year with ∆K s 0.14 mag and ∆(H -K s ) 0.07 mag. The second time occurs for ∼70 days in the third year with ∆K s ∼0.  b Entries with -mark stars with photometry below the survey J band completeness limit.
c Entries with -mark stars with photometry below the survey H band completeness limit.
bluer making the mechanism identification tentative. The YSO classification for this star is unknown. The J band photometry is too dim in 11 irregular variables to identify the variability mechanism. Of these 11 variables, the H band is also too dim in 5 stars.
The variability in ISO-Oph 87 is peculiar due to two flare-like events that occur in the K s photometry on approximately 1400 and 1700 2MASS MJD.. Fig 32 contains the K s , (J -H ) and (H -K s ) photometry for this Class II YSO. If the events are truly related to an increase in stellar activity, the star is expected to become bluer in both (J -H ) and (H -K s ). However, no change is seen in the (J -H ) color and the star reddens in (H -K s ). The first event lasts for ∼10 days with ∆K s 0.20 mag and ∆(H -K s ) 0.13 mag. The second event occurs for ∼6 days and ∆K s 0.18 mag and ∆(H -K s ) 0.16 mag. Two additional, lower amplitude spikes in the K s photometry between 1600 and 1700 2MASS MJD might also be similar flare-like events.
Irregular variables have the smallest percentage (68%) of stars located "on cloud". Only 9 stars in this sub-category are classified as a YSO: 1 Class I (8%), 7 Class II (21%) and 1 Class III (9%). 5 Most (16 stars) of the 22 candidate ρ Oph members are irregular variables. If these candidate members are YSOs then the fraction of YSO irregular variables is comparable to the fraction of periodic and LTV YSOs.

Examples of Multiple Variability Mechanisms
As discussed in the Introduction, variability studies indicate young stars sometimes exhibit complex photometric behavior believed to result from multiple variability mechanisms acting concurrently. For most stars in this survey, only the highest amplitude variability can be confidently characterized. However, a lower amplitude, second type of variability is definitely seen in 7 variable stars. Four of these seven stars (YLW 1C, 2MASS J16272658-2425543, YLW 10C, WL 4) show evidence for two separate, yet statistically distinct periodic variations. The remaining three stars (WL 20W, ISO-Oph 126, WL 15) are periodically variable underneath a higher amplitude, long time-scale variation. The methods used to identify both the primary and secondary variabilities in these stars is discussed in §4.1. The following subsections contain detailed discussions for each of these stars except WL 4 which is described in Plavchan et al. (2008).
• YLW 1C (ISO-Oph 86): This CTTS exhibits both sinusoidal-like and eclipse-like periodic variability at 2 distinct periods; the periods are distinct from each other to a 20σ confidence level. Fig 33 contains   The events become redder as the star brightens rather than becoming bluer as expected for stellar flares.
to the sinusoidal-like period, P = 5.7792 ± 0.0085 days. The peak-to-trough ∆K s amplitude for this variability is 0.14 mag. The variability in both the (J -H ) color and (H -K s ) color is not correlated with the K s variability. This favors variability caused by rotational modulation of a cool starspot(s). Fig 34 contains the K s , (J -H ) and (H -K s ) photometry folded to the eclipse-like period, P = 5.9514 ± 0.0014 days. The ∆K s eclipse depth is 0.29 mag. The (H -K s ) color reddens during the eclipse event consistent with variability caused by extinction. This behavior is not seen in the (J -H ) color.
Since the J band photometry is near the survey completeness limits, the absence of a clear reddening trend may be due to low signal-to-noise in the color curve.
The two periods for YLW 1C is consistent with the interpretation that these events are true occultation events, as proposed for AA Tau (see §5.1.2). The short period, arising from a stellar surface feature(s), traces the stellar rotation rate. The longer period suggests the occultation of the star by an obscuration located just beyond the circumstellar disk co-rotation radius. Following the anaylsis described in §5.1.2, the occulter of YLW 1C is located 7.2 R from the host star and the duration of the occultation is ∼6.4 hours corresponding to a minimum occulter diameter of ∼2.0 R . The reader is reminded this diameter represents the extent within the orbital path and makes no claim on any preferential occulter shape. The eclipse depth is ∆K s = 0.29 mag.
The large occulter diameter argues against the direct detection of a hot protoplanet. However recent imaging results suggest that gas giant planets maybe considerably extended in the mass accretion phase (Quanz et al. 2013;Kraus & Ireland 2012). If true in this case, this would demonstrate the existance of a hot protoplanet with a period of 6 days very near the peak in the period distribution for exoplanets (Wright et al. 2012). Alternatively, the event could be caused by an occultation of a warped portion of a circumstellar disk. This scenario has been proposed to explain the near-to mid-IR variability in LRLL 31 (Flaherty & Muzerolle 2010;Flaherty et al. 2012). While most YSO disk models invoke axisymmetry, objects such as YLW 1C are prompting the creation of more complex models.
• 2MASS J16272658-2425543: This CTTS, designated 'J543' hereafter, is another star exhibiting both sinusoidal-like and eclipse-like variability with two distinctly different periods. Fig 35 contains  is 0.17 mag. Both the (J -H ) and (H -K s ) colors redden during in the eclipse event consistent with variable extinction.
The physical interpretation for the observed variability is identical to that of YLW 1C. The 1.6 day period corresponds to the stellar rotation rate and the 3.0 day period arises from a periodic occultation by an asymmetry in the inner circumstellar disk. The occulter size and distance are 2.8 R and ∼3.0 R . Differing from YLW 1C, the occulter for J543 is located approximately a stellar radius beyond the co-rotation radius.
• YLW 10C (ISO-Oph 122): YLW 10C is the third CTTS where two distinct periods are identified. Fig 37 contains  Given the sinusoidal-like and eclipse-like variability is very similar to both YLW 1C and J543, the same physical interpretation is proposed for this star. However, unlike YLW 1C, the sinusoidal-like variability, presumed to trace the stellar rotation rate, has a longer period by 3.1 hours than the periodic occultations. This places the hypothetical occulter within the co-rotation radius. The size and distance to the occulter are 3.31 R and 3.75 R . The occulter is located within the dust sublimation radius as computed using the formulism of Jura & Turner (1998). This formulism is only an approximation as it does not take into account dust evaporation and condensation rates, grain size, or grain composition.
• WL 20W (YLW 11B, ISO-Oph 126): This CTTS is both periodically variable and variable over a long time-scale. As such, WL 20W is designated both a periodic variable and a LTV. Fig 39 contains the K s , (J -H ) and (H -K s ) photometry for this star. The long time-scale variability begins on ∼1600 2MASS MJD and has a timescale of 122 days. The ∆K s depth is 0.26 mag. Both the (J -H ) and (H -K s ) colors become redder during the long time-scale variation, consistent with variable extinction. The periodic signal is not significant unless the time-series affected by the long time-scale variability is omitted from analysis by the PA (see Fig 11). Fig 40 contains the K s , (J -H ) and (H -K s ) photometry folded to the sinusoidal-like period, P = 2.1026 ± 0.0060 days. The peak-to-trough ∆K s amplitude for the sinusoidal-like variability is 0.19 mag.
Neither the (J -H ) nor the (H -K s ) color is correlated to the K s variability. This favors rotational modulation by cool starspots as the variability mechanism.
The variability of AA Tau has been cited as an explanation for the periodic eclipsing   variables in this survey. However recently it has been discovered that AA Tau is exhibiting a long time-scale dimming on the order of 2-3 mag in the V band (Bouvier et al. 2013). This long time-scale variation is superimposed on top of the periodic variability.
Additionally the system appears to become bluer in this dim state; a phenomenon seen in UX Ori-type variables (Grinin et al. 1991;Herbst et al. 1994). The physical interpretation for UX Ori-type variability is that the star dims due to an asymmetric optically thick occulter beyond the inner circumstellar disk. The bluer color represents a larger contribution of scattered starlight off the occulting material.
This scenario is an alternative explanation than variable mass accretion for the long time-scale variation observed in WL 20W. The eclipse depth corresponds to a 2.9 mag dimming when converted to the V band by using the extinction coefficients given in Cohen et al. (1981); this is consistent with the AA Tau long time-scale variation. The (J -H ) and (H -K s ) colors do become bluer during the long time-scale variation also consistent with the observations of AA Tau. Our overall interpretation for the variability in WL 20W is a central star rotating with a 2 day period that is occulted by a pocket of optically thick material located beyond the inner circumstellar disk.
It is worth noting this star belongs to a triple system that is spatially resolved in the mid-IR. Ressler & Barsony (2001) show the most variable member is, in fact, WL 20S. They classify this source as Class I through SED fitting of mid-IR photometry. They show that WL 20E and WL 20W are nearly constant on decadal timescales whereas the flux of WL 20W increased sixfold in 15 years. As the largest separation between these components is 3.66 , the large aperture size in our work (4 ) includes all three stars and thus cannot rule out the possibility that the measured variability arises from this southern component.
• ISO-Oph 126: This WTTS is similar to WL 20W in that it exhibits both periodic variability and a long time-scale variation. ISO-Oph 126 is also designated both a periodic variable and a LTV. Fig 41 contains the K s , (J -H ) and (H -K s ) photometry for this star. The J band photometry is below the survey completeness limits. Therefore the (J -H ) color is deemed unreliable for analysis. The long time-scale variation dominates the photometry prior to 1400 2MASS MJD with a time-scale of 349 days. The ∆K s depth of this variation from the continuum brightness is 0.10 mag. The (H -K s ) color becomes redder during the long time-scale variation as the star dims. This is consistent with extinction as the variability mechanism. The PA identifies a significant periodic signal when only the portion of the time-series after 1400 2MASS MJD is analyzed (see Fig 11). The K s shows considerable structure during ingress, but is nearly stable during the "eclipse" egress. The red line indicates the mean magnitude of the continuum photometry. The ∼19 day period is clearly evident and appears to continue even through a ∼1 mag drop in K s band flux. The photometry during the larger amplitude flux decrease is too sparse to confidently determine a time-scale. A lack of a trend in the (H -K s ) color during this event highly suggests against extinction except by an opaque occulter. to as caused by cool starspots, to date, is ∆V = 0.63 mag (Strassmeier et al. 1997). This amplitude is nearly 0.3 mag lower than that observed for WL 15. In addition, the contrast between the starspot and surrounding photosphere increases toward bluer wavelengths. Therefore the amplitude will be even larger in the optical. The origin behind the long time-scale variation is equally peculiar. As the (H -K s ) color does not become redder as the star dims, this seems to eliminate extinction as variability mechanism. However, variability by rotational modulation of surface features seems inplausible given the timescale and amplitude of the variation.

SUMMARY
High precision, high cadence J, H, K s photometry is obtained for 7815 stars in the direction of the ρ Oph molecular cloud with a temporal baseline of ∼ 2.5 yrs. Spurious detections, partially resolved doubles or galactic contamination are eliminated from the photometry. The target sample meeting the specifications for time-series variability analysis includes 1678 stars. A seven point variability test is used to identify 101 variable stars, which is 6% of the parent sample. These tests are sensitive to variability on a variety of different time-scales and forms (e.g. sinusoidal, 'eclipse-like', etc.).
Of the 101 stars in the variable catalog, 15% are located "on cloud" while only 1% lie within the "field". Location "on cloud", variability and (H -K s ) colors redder than a 3 Myr isochrone are used to estimate membership in the ρ Oph star forming region. This method has identified 22 stars as previously unknown candidate ρ Oph members.
The effects of observing strategy on variability detection and measured amplitudes is investigated by comparing this work to the ρ Oph variability study performed by Alves de Oliveira & Casali (2008). These two studies have 464 stars in common; AC08 identified 7% as variable stars and this work identifies 18%. The increase in detection fraction is not caused by different sensitivities in the separate variability criteria used in each survey. This work also found both a higher K s and (H -K s ) color amplitude in 25 stars identified as variable in both surveys than measured by AC08. Therefore a high cadence observational strategy will discover more variables within a given set of stars and it will more accurately characterize intrinsically higher amplitude variability.
The K s variability and stellar color behaviors are used to estimate the physical mechanism responsible for the variability. Rotational modulation by long-lived cool starspots is expected to produce colorless, periodic variability. Rotational modulation of long-lived hot starspots (e.g. accretion) is, also, expected to be periodic, while short-lived starspots (e.g. flares) is not. The star becomes bluer as it brightens in both cases. Extinction induced variability is either periodic or exists on long time-scales based on the geometry of the occulter relative to the star. Changes in the mass accretion rate onto the star is not expected to be periodic, but occur on time-scales ranging from days to years. As this rate changes, the star becomes bluer as the star dims.
Identifying periodic variability within the variable catalog is done via a newly improved period-searching algorithm, the Plavchan algorithm. The algorithm tests tens of thousands of periods with uniform frequency sampling between 0.1 to 1000 days. This is done by comparing the observed light curve to a dynamically generated prior. The statistical significance of individual periods is computed via two methods: the distribution of power values at other periods in the same periodogram and the distribution of maximum power values for all sources in an ensemble survey. The Plavchan algorithm finds periodic variability in 32% of the variable catalog with periods ranging from 0.49 to 92 days.
The periodic variables are split into two sub-categories: sinusoidal-like and eclipse/inverse eclipse-like. Sinusoidal-like periodic variability describes a sinusoidal-like change in the observed flux when the time-series is folded to the most significant period. Rotational modulation by cool starspots is believed to be the common variability mechanism in this subcategory. Sinusoidal-like periodic variables are found in each YSO class (3 Class I, 8 Class II, 8 Class III). Eclipse-like periodic variability results in discrete drops, or "dips", in the observed flux when the time-series is folded to the most significant period. Periods range from 2 to 8 days with the duration of these dips lasting less than 30% of one periodic epoch. This sub-category contains 6 stars with a median peak-to-trough ∆K s and ∆(H -K s ) color amplitudes of 0.31 and 0.11 mag, respectively. Variable extinction is the likely mechanism for the eclipse-like variations. All stars in this sub-category are Class II YSOs. The inverse eclipse-like variables, WL 4 and YLW 16A have periods of 65.61 and 92.3 days, respectively. The variability mechanism proposed in both cases is the periodic obscuration of one component in a close binary by a warped circumbinary disk.
In half of the eclipse-like variables (YLW 1C, 2MASS J16272658-2425543, YLW 10C) an additional statistically significant period is identified. This sinusoidal-like periodic variability coupled with the presence of "dips" suggests a rapidly rotating spotted star occulted by a clump of optically thick material in the inner accretion disk. These stars strengthen the interpolation posed to explain the variability of other YSO AA Tau-like variables (Morales-Calderón et al. 2011). The periods corresponding to periodic occultations in YLW 1C and YLW 10C are located near thier respective co-rotation radii. The mechanism driving these occulations could arise from a warped inner circumstellar disk caused by an inclined magnetic dipole, or could be the pre-natal cloud of a forming hot Jupiter.
Long time-scale variables is a variability subclass, containing 31 stars, where the measured flux increases or decreases consistently over months or years. The variability time-scale is measured using a differencing technique and approximates the time between maximum and minimum brightness. The measured time-scales range from 64 to 790 days. The peakto-trough ∆K s amplitudes range from 0.05 to 2.31 mag and the peak-to-trough ∆(H -K s ) color amplitudes range from 0.06 to 1.32 mag. Variable extinction and variable accretion rates are both equally likely to cause long time-scale variability. This subclass contains 25 known YSOs with 7 Class I, 15 Class II and 3 Class III stars.
The time-series of irregular variables are aperiodic and do not vary over any discernible time-scale. This subclass contains more members (40) than either the periodic or long timescale subclasses. The peak-to-trough ∆K s amplitudes range from 0.04 to 1.11 mag and peak-to-trough ∆(H -K s ) color amplitudes range from 0.05 to 0.75 mag. No single dominant variability mechanism explains irregular variability. Only 9 known YSOs (1 Class I, 7 Class II, 1 Class III) are irregular variables.
The CTTS WL 20W and the WTTS ISO-Oph 126 are similar in both have a long timescale variation superimposed onto a periodic signal. In both cases, the physical mechanisms for the variability is consistent with an occultation of a rapidly rotating spotted star by optically thick material outside the inner accretion disk. For WL 20W, the sinusoidal-like variability has a period of 2.1026 days and a peak-to-trough ∆K s amplitude of 0.19 mag. The long time-scale variability has a duration of 122 days with a ∆K s eclipse depth of 0.26 mag. For ISO-Oph 126, the sinusoidal-like variability has a period of 9.114 days and a peakto-trough ∆K s amplitude of 0.06 mag. The long time-scale variability has a duration of 349 days with a ∆K s eclipse depth of 0.10 mag.
The very high amplitude periodic variability measured in the Class I star WL 15 is not consistent with any proposed mechanism. The observed 47 day colorless decrease in K s band brightness of ∼1 mag is also not easily estimated.
From cross-referencing the target sample with two previous surveys, 72 stars have been assigned a YSO classification (13 Class I, 47 Class II, 12 Class III). The variability fraction of these YSOs is 79%. The variability fraction differs according to YSO class with 92% of Class I and Class III stars identified as variable; this fraction drops to 72% for Class II stars. The amplitude of both brightness and color variability are decreasing functions of YSO class. The median peak-to-trough ∆K s amplitude for Class I, II and III stars are 0.77, 0.31 and 0.08 mag, respectively. In addition, the median peak-to-trough ∆(H -K s ) color amplitudes are 0.81, 0.21 and 0.07 mag for each class respectively.
We will now express our thanks to our anonymous referee for their conversation and comments. We would also like to express our thanks to Mary Barsony for her many insightful questions and suggestions. This publication makes use of 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. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. The authors would like to acknowledge the Infrared Processing and Analysis Center's Visiting Graduate Student program for providing the opportunity and funding for this publication.

APPENDIX. Empirical Determination of χ 2 n 0 Dependence on Parameters
In this Appendix an analysis of the dependence of the PA periodogram χ 2 n 0 power values on the number of observations and periodogram parameters n 0 and p is presented. The ensemble survey of mostly non-variable stars is used to carry out this analysis and to present an alternative approach to evaluate the statistical significance of periodogram power values.
First, a random subset of 180 stars is chosen from the survey collection of 1678 stars. These stars are evenly distributed in N obs and K s magnitude. For a given set of parameters p and n 0 , the maximum χ 2 n 0 periodogram power value is computed for the 180 stars. Since other algorithms exist that specialize in finding periodic sources with low numbers of detection (Dworetsky 1983, N obs ∼20), test cases are limited to 0.04 < p and ≤ 0.5 and 12 < n 0 ≤ 250. Fig A-1 shows the dependence of χ 2 n 0 on N obs for the 180 stars with a particular set of p and n 0 . This dependence is somewhat expected -a smaller number of observations can result in an increase in the likelihood for false-positive periodogram peaks.
The distribution of χ 2 n 0 values as a function of N obs is well-described by the functional form: where a, b and c represent real numbers that differ for a given p and n 0 . As the power law index decreases below 1.5 for p < 0.04 and/or n 0 < 12, these ranges are excluded from the analysis. Eqn A-1 is found via trial and error to minimize the residuals when compared to a variety of functional forms tested, rather than from an analytic derivation based upon first principles. The values p and n 0 = 40 -used throughout this paper to identify periodic variablesare an optimal choice of parameters for this survey. They yield the smallest residuals when the 180 test cases are fit in Eqn A-5. Thus, for this survey the maximum peak power in the periodogram for a non-variable star is approximately given by the expression: The validity of the numerical fits (Eqns A-4(a-c)) is verified using an additional 18 test cases with randomly selected values of p and n 0 . The χ 2 n 0 versus N obs distributions are again fit using Eqn A-1, yielding 18 "observed" a, b and c values for each pair of p and n 0 . Predicted values for a, b and c are found by using Eqns A-4(a-c) and compared to the "observed" values. The mean percent errors between the observed and predicted values for the 18 test cases are -1.1±3.4% for f a , 1.1±3.5% for f b and -0.1±1.8% for f c . Table A-1 contains the percent errors for each of the individual 18 test cases. Eqns A-4(a-c) therefore adequately predict the values a, b and c in Eqn A-1 for any set of parameters p and n 0 within the parameter space explored in this survey. Thus, this demonstrates the PA algorithm is reasonably "well-behaved." Eqn A-5 could be applied to different surveys and cadences. However, the particular numerical values in Eqns A-4 and A-6 for a, b and c likely depend on the specific cadence of a survey.
To evaluate the statistical significance of a peak period power value in a periodogram, it does not suffice to identify the expected value for the ensemble survey. The scatter about the expected value is also necessary. This scatter, or standard deviation (σ), of the peak period power values about the expected value depends on N obs in a predictable fashion for this survey (Fig A-3). To characterize this scatter, the scatter for each survey star is grouped into bins as a function of N obs , with a bin size of 25. An average σ is computed for each bin and an empirical fit to this distribution is made, given by Eqn A-7: σ(N obs ) = 2.3790 N obs − 21.6449 + 0.0105 (A-7) Putting it all together, now an affirmative periodicity condition can be defined. The statistical significance of a measured period for a star in this survey is simply a function of the number of observations. Based on visual inspection of star light curves, this condition is defined by the following:  Fig. A-3.-The maximum PA periodogram power value for a star in this survey is found to be well-described by Eqn A-5 (Fig A-1). The scatter about the value predicted by Eqn A-5, σ, is also found to be dependent on N obs as shown in this figure.