# Let the Great World Spin: Revealing the Stormy, Turbulent Nature of Young
Giant Exoplanet Analogs with the Spitzer Space Telescope

Johanna M. Vos Department of Astrophysics, American Museum of Natural History,
Central Park West at 79th Street, NY 10024, USA  Jacqueline K. Faherty
Department of Astrophysics, American Museum of Natural History, Central Park
West at 79th Street, NY 10024, USA  Jonathan Gagné Planétarium Rio Tinto
Alcan, Espace pour la Vie, 4801 av. Pierre-de Coubertin, Montréal, Québec,
Canada  Institute for Research on Exoplanets, Université de Montréal, 2900
Boulevard Édouard-Montpetit Montréal, QC Canada H3T 1J4  Mark Marley
University of Arizona Department of Planetary Sciences and Lunar and Planetary
Laboratory  Stanimir Metchev Department of Physics & Astronomy and Centre for
Planetary Science and Exploration, The University of Western Ontario, London,
Ontario N6A 3K7, Canada  John Gizis Department of Physics and Astronomy,
University of Delaware, Newark, DE 19716, USA  Emily L. Rice  Department of
Astrophysics, American Museum of Natural History, Central Park West at 79th
Street, NY 10024, USA  Macaulay Honors College, 35 W. 67th Street, New York,
NY 10023, USA  Kelle Cruz Department of Astrophysics, American Museum of
Natural History, Central Park West at 79th Street, NY 10024, USA  Department
of Physics and Astronomy, Hunter College, City University of New York, 695
Park Avenue, New York, NY 10065, USA  Flatiron Institute, 162 Fifth Avenue,
New York, NY 10010, USA

###### Abstract

We present a survey for photometric variability in young, low-mass brown
dwarfs with the Spitzer Space Telescope. The 23 objects in our sample show
robust signatures of youth and share properties with directly-imaged
exoplanets. We present three new young objects: *2m0349, *2m0951 and *2m0718.
We detect variability in 13 young objects, and find that young brown dwarfs
are highly likely to display variability across the L2–T4 spectral type range.
In contrast, the field dwarf variability occurrence rate drops for spectral
types $>$L9. We examine the variability amplitudes of young objects and find
an enhancement in maximum amplitudes compared to field dwarfs. We speculate
that the observed range of amplitudes within a spectral type may be influenced
by secondary effects such as viewing inclination and/or rotation period. We
combine our new rotation periods with the literature to investigate the
effects of mass on angular momentum evolution. While high mass brown dwarfs
($>30M_{\mathrm{Jup}}$) spin up over time, the same trend is not apparent for
lower mass objects ($<30M_{\mathrm{Jup}}$), likely due to the small number of
measured periods for old, low-mass objects. The rotation periods of companion
brown dwarfs and planetary-mass objects are consistent with those of isolated
objects with similar ages and masses, suggesting similar angular momentum
histories. Within the AB Doradus group, we find a high variability occurrence
rate and evidence for common angular momentum evolution. The results are
encouraging for future variability searches in directly-imaged exoplanets with
facilities such as the James Webb Space Telescope and 30-meter telescopes.

brown dwarfs — exoplanet atmospheres — atmospheric variability

††journal: AAS Journals††facilities: Spitzer Space Telescope, Gaia DR2
(Gaia2016; Gaia2018), Magellan I FIRE (Simcoe2013)††software: BANYAN $\Sigma$
(Gagne2018b), SEDkit (Filippazzo2015), emcee (fm2013), matplotlib
(Hunter2007), celerite2 (Foreman-Mackey2017)

[2MASS J0045$+$16]2m00452MASS J00452143$+$1634446 [2MASS
J0501$-$00]2m05012MASS J05012406$-$0010452 [2MASS J1425$-$36]2m14252MASS
J14252798$-$3650229 [2MASS J0001$+$1535]2m00012MASS J00011217$+$1535355 [2MASS
J0030$-$1450]2m00302MASS J0030300$-$145033 [2MASS J0031$+$5749]2m00312MASS
J00310928$+$5749364 [2MASS J0153$-$6744]2m01532MASS J01531463$-$6744181 [2MASS
J0326$-$2102]2m03262MASS J03264225$-$2102057 [2MASS J0342$-$6817]2m03422MASS
J0342162$-$681732 [2MASS J0349$+$0635]2m03492MASS J03492367$+$0635078 [2MASS
J0355$+$1133]2m03552MASS J03552337$+$1133437 [2MASS J0447$-$1216]2m04472MASS
J04473039$-$1216155 [2MASS J0459$-$2853]2m04592MASS J04590034$-$2853396 [2MASS
J0506$+$5236]2m05062MASS J05065012$+$5236338 [2MASS J0642$+$4101]2m06422MASS
J06420559$+$4101599 [2MASS J0718$-$6415]2m07182MASS J07180871$-$6415310 [2MASS
J0809$+$4434]2m08092MASS J08095903$+$4434216 [2MASS J0951$-$8023]2m09512MASS
J09512690$-$8023553 [2MASS J1551$+$0941]2m15512MASS J15515237$+$0941148 [2MASS
J1647$+$5632]2m16472MASS J16471580$+$5632057 [2MASS J1741$-$4642]2m17412MASS
J17410280$-$4642218 [2MASS J2002$-$0521]2m20022MASS J20025073$-$0521524 [2MASS
J2117$-$2940]2m21172MASS J21171431$-$2940034 [2MASS J2154$-$1055]2m21542MASS
J21543454$-$1055308 [2MASS J2206$+$3301]2m2206+332MASS J22062520$+$3301144
[2MASS J2206$-$4217]2m2206-422MASS J2206450$-$421721 [WISE
J2216$+$1952]2m2216WISE J221628.62$+$195248.1 [2MASS J2322$-$6151]2m23222MASS
J23225299$-$6151275 [2MASS J2343$-$3646]2m23432MASS J23433470$-$3646021 [2MASS
J2139$+$02]2m2139 2MASS J21392676$+$0220226 [SIMP J0136$+$09]s0136 SIMP
J01365662+0933473 [WISE J0047$+$6803]w0047WISE J004701.06$+$680352.1 [2MASS
J2244$+$2043]2m22442MASS J2244316$+$204343 [SDSS 1110$+$0116]s1110SDSS
111010$+$011613

##  1 Introduction

The atmospheres of brown dwarfs and exoplanets with temperatures $<2300$K are
characterised by the presence of condensate clouds, which shape their emergent
colors and spectra. Based on the handful of directly-imaged exoplanets
discovered and investigated to date (e.g. $\beta$ Pictoris b, HR8799 bcde,
2M1207 b, 51 Eri b; Lagrange2010; Marois2008; Marois2010; Chauvin2004;
Macintosh2015), it is clear that in-depth investigations of these planets will
hinge on a thorough understanding of condensate cloud chemistry, and its
relationship with surface gravity, temperature and metallicity.

Time-resolved atmospheric monitoring of brown dwarfs and exoplanets directly
probes atmospheric features such as condensate clouds. In the last decade, our
understanding of clouds in brown dwarf atmospheres has been revolutionized by
a number of infrared studies of rotationally modulated variability. Large
surveys (Buenzli2014; Radigan2014; Wilson2014; Metchev2015; Eriksson2019) have
revealed that variability is common across L and T spectral types
(Metchev2015), with a possible enhancement in variability amplitude and
occurrence rate across the L/T transition (Radigan2014). First results from
dwarf monitoring of Y dwarfs with the Spitzer Space Telescope suggest that
variability is also common at the lowest temperatures (Cushing2016;
Esplin2016; Leggett2016). Multi-wavelength characterization studies have
provided insight into the vertical atmospheric structure of brown dwarfs and
planetary-mass objects, since different wavelengths probe different
atmospheric pressures. For example, comparing variability amplitudes at
different wavelengths provides information on the types of clouds present
(Artigau2009; Radigan2012; Apai2013; Kellogg2017; Biller2018) and phase shifts
between light curves of different wavelengths have revealed multiple cloud
layers (Buenzli2012; Biller2013).

Physical mechanisms that can explain observed spectroscopic and photometric
variations without clouds have also been proposed. Models presented by
Morley2014 for mid-late T dwarfs and Y dwarfs show that hot spots can lead to
variability. Tremblin2020 suggest that variability may be driven by
temperature fluctuations that arise as a result of thermochemical
instabilities. Moreover, brown dwarfs are known to exhibit magnetically-driven
aurorae, which may also drive infrared variations (Hallinan2015; Kao2016;
Pineda2017, e.g.). Spectroscopic variability monitoring will be necessary to
confirm the spectral signatures predicted by each model and to disentangle the
effects of each atmospheric process.

Photometric and spectroscopic monitoring is likely to reveal important
insights into the atmospheres of directly-imaged exoplanets. Indeed, there
have already been a number of variability detections in wide-orbit companions
with masses in the range $\sim 10-30~{}M_{\mathrm{Jup}}$ (Naud2017; Zhou2016;
Zhou2018; Zhou2020). For the closer in directly-imaged exoplanets, there have
been two attempts to detect photometric variability in the HR8799 bc planets
(Apai2016; Biller2021). However, neither study detected variability in the
planets, and placed upper limits of $5-10\%$ on the variability amplitudes.
The James Webb Space Telescope (JWST) will enable unprecedented sensitivity
and stability for future variability observations of directly-imaged
exoplanets such as HR8799 bcde and $\beta$ Pictoris b.

Young, isolated brown dwarfs are emerging as a bridge population to giant
exoplanets. This population is vastly easier to study without the glare of a
bright host star, and exist in greater numbers than the current population of
known directly-imaged exoplanets (Allers2013; Cruz2009; Faherty2016;
Gagne2015c; Liu2016). The colors and spectra of this new population are
remarkably similar to the directly-imaged exoplanets; they exhibit redder
near-infrared colors than higher mass field brown dwarfs with the same
spectral type and appear fainter than their equivalent temperature, more
massive counterparts (Faherty2016; Liu2016). These differences are generally
explained by enhanced atmospheric clouds and/or vertical atmospheric mixing
(e.g. Marley2012; Barman2011), arising from lower surface gravity of the
younger objects. Since surface gravity appears to have a large effect on the
atmospheric and cloud properties of brown dwarfs and directly-imaged
exoplanets, it follows that variability may be distinct for the low-gravity
population.

A number of variability studies have focused on young, low-gravity free-
floating and companion objects, and initial studies of young free-floating
objects suggest that variability may be enhanced compared to field brown
dwarfs. The first detection of cloud-driven variability in a planetary-mass
object for the 23 Myr old $\beta$ Pictoris moving group member PSO J318.5$-$22
showed the largest variability amplitude to date in an L dwarf (Biller2015).
This initial discovery was followed by a notable sample of low-gravity objects
showing high-amplitude variability (Lew2016; Zhou2016; Vos2018; Zhou2020;
Bowler2020). Metchev2015 note a potential enhancement in amplitude in a sample
of five low-gravity L3-L5.5 brown dwarfs at $3.6~{}\mu$m and $4.5~{}\mu$m.
Vos2020 also note an apparent amplitude enhancement for young late-L dwarfs at
$4.5~{}\mu$m based on a compilation of all of the objects monitored with
Spitzer to date. Recently, the high-amplitude variables s0136 and 2m2139 were
identified as members of the $\sim 200$ Myr old Carina-Near moving group
(Gagne2017; ZJZhang2021). There is also evidence for a higher variability
occurrence rate in young brown dwarfs. Vos2019 carried out a large, ground-
based survey for $J$-band variability in isolated. low-gravity dwarfs, finding
that the young population has a variability occurrence rate of
$30^{+16}_{-8}\%$, compared to $11^{+11}_{-4}\%$ of field dwarfs. It is clear
that a large, unbiased survey for variability in young, low-gravity objects is
necessary to test the tentative amplitude and occurrence rate enhancements
reported to date.

In this paper, we present the first large survey for variability in young,
low-gravity brown dwarfs in the mid-infrared, and the most sensitive survey to
date at any wavelength. We describe the sample in Section 2, including the
discovery of three new brown dwarfs with evidence of youth in Section 2.1. We
describe the Spitzer observations, data reduction and light curve analysis in
Sections 3 and 4. We discuss the youth and variability properties of each
target in detail in Section 5. We discuss the variability occurrence rate of
young brown dwarfs in comparison with field brown dwarfs in Section 6, and the
variability amplitudes in Section 7. We put the rotation rates of our young
sample in context with literature data in Section 8, and specifically discuss
the variability properties of AB Doradus members in Section 9. Finally, we lay
out some lessons for variability searches in directly-imaged exoplanets in
Section 10. We summarize our results in Section 11.

Figure 1: Spectral type - color diagram of our sample of young, low-gravity
brown dwarfs (large pink circles) compared to the field population (small grey
circles). Our Spitzer sample consists of 26 L3–T5 brown dwarfs with evidence
of youth and/or low-gravity. The L dwarfs appear redder than the field L dwarf
sequence, which is a signature of their low surface gravity. The field brown
dwarf population is drawn from Best2020.

##  2 An Age-Calibrated Sample of Exoplanet Analogs

We designed our survey to test the variability properties of young, low-
gravity brown dwarfs by comparing our results to the field brown dwarfs
observed by Metchev2015. The similar sample size, observing strategy and
magnitudes of our sample and the Metchev2015 sample allow us to robustly
compare variability properties between them. We present our sample, spectral
types and magnitudes in Table 1 and show our sample in a color-magnitude
diagram in Figure 1. Our sample spans spectral types L2 – T5, and the L dwarfs
in the sample appear redder than the field L dwarf sequence, often a signature
of their low gravity.

###  2.1 Three new brown dwarfs with evidence of youth from the BASS Ultracool
Survey

We present the discovery of three brown dwarfs which display evidence of
youth: *2m0349, *2m0718 and *2m0951. All three objects were discovered as
candidate young objects through the BASS-Ultracool survey (Gagné et al., in
prep) based on their kinematics, which are presented in Table
LABEL:tab:kinematics and discussed in detail in Section 5. Here we present
their near-IR FIRE/Prism spectra and investigate their spectral signatures of
youth.

Figure 2: FIRE/Prism spectra of three new young brown dwarfs presented in this
paper – 2m0349 (L5$\gamma$), 2m0718 (T5) and 2m0951 (L5 (pec)). Spectral
templates for an L5$\gamma$, and L5 and a T5 are overplotted in black. The
spectral resolution ranges from 500 at $J$ to 300 at $K$. These objects were
discovered and classified as young as part of the BASS-Ultracool survey (Gagné
et al., in prep).

####  2.1.1 FIRE/Prism Spectra: Observations and Analysis

We used the FIRE (Simcoe2013) spectrograph on the 6.5 m Baade Magellan
telescope to obtain near-infrared spectra of the three targets to investigate
their spectral properties. Observations were taken on 2018 Feb 01, 2015 Nov 24
and 2016 Jan 22 for 2m0349, 2m0718 and 2m0951 respectively using prism
longslit mode, using ABBA nods in each case. Quartz lamps, ThAr lamps and
telluric standards were taken either before our after each target. Airmass
ranged from 1.2-1.6 for the three observing sequences. All FIRE/Prism data
were reduced and combined using a modified version of the FireHose package as
described in Gagne2015a. We show the $R\sim 400$ FIRE/prism spectra of these
three targets in Figure 2.

####  2.1.2 FIRE/Prism Spectra: Signatures of Youth

We perform a visual spectral classification on each of our new discoveries and
an index-based classification on 2m0349 and 2m0951 to assess the evidence of
youth and low-gravity by comparing our L-type spectra to spectral standard
templates reported by Cruz2018 and T-type spectrum to known T-type spectral
standards.

Based on a visual classification of our FIRE/Prism spectrum of 2m0349, we find
that it is best matched by an L5$\gamma$ spectrum. Using the Allers2013
indices, we find that the spectrum is classified as an L4 vl-g. Thus, the
spectrum provides strong evidence of youth.

We assign 2m0951 a visual classification of L5 (pec). The spectrum is not well
matched by the spectral type templates, particularly in the $H$-band, where
there is a ‘bump’ at $\sim 1.6~{}\mu$m that is not shown in the L5 templates.
The Allers2013 index-based classification assigns an L1.6 int-g spectral type.
The disparity between the two spectral types suggests that the peculiarity of
the spectrum prohibits a robust spectral classification. If we force the
Allers2013 to use an L5 spectral type the gravity classification is int-g.

Visually, 2m0718 is best matched by a T5 spectral type. It is very red,
emitting a higher flux in $K$-band compared to the T5 template. Interestingly,
the spectrum of the planetary-mass AB Doradus member SDSS J111010.01+011613.1
(Gagne2015b) provides an excellent fit to our spectrum. With a spectral type
of T5, we cannot perform the index-based classification of Allers2013 on this
target.

Target Name | Spt Opt | SpT IR | Refs | $J_{\mathrm{2MASS}}$ | $K_{\mathrm{2MASS}}$ | $(J-K)$ | [3.6] maga
---|---|---|---|---|---|---|---
*2m0001 | $\cdots$ | L4$\beta$ | 1 | 15.46 | 13.62 | 1.84 | $12.613\pm 0.004$
*2m0030 | L7 | L4-L6$\beta$ | 2, 3 | 16.28 | 14.48 | 1.80 | $13.365\pm 0.008$
*2m0031 | $\cdots$ | L9 | 4 | 14.95 | 13.22 | 1.74 | $11.956\pm 0.003$
*2m0153 | L2 | L3$\beta$ | 3,6 | 16.41 | 14.42 | 1.99 | $13.349\pm 0.005$
*2m0326 | L5$\beta$ | L5$\beta/\gamma$ | 1,3 | 16.13 | 13.92 | 2.208 | $12.559\pm 0.003$
*2m0342 | L4$\gamma$ | $\cdots$ | 1 | 16.854 | 14.541 | 2.313 | $13.579\pm 0.006$
*2m0349 | $\cdots$ | L5$\gamma$ | TW | 16.742 | 15.077 | 1.665 | $14.181\pm 0.007$
*2m0355 | L5$\gamma$ | L3-L6$\gamma$ | 3, 5 | 14.05 | 11.526 | 2.524 | $10.339\pm 0.002$
*2m0447 | $\cdots$ | T2 | 7 | 16.48 | 15.55 | 0.93 | $11.868\pm 0.002$
*2m0459 | $\cdots$ | L7 | 8 | 17.429 | 15.318 | 2.111 | $13.924\pm 0.006$
*2m0506 | $\cdots$ | T4.5 | 7 | 15.75 | 15.6 | 0.15 | $14.251\pm 0.008$
*2m0642 | $\cdots$ | L9 | 7 | 16.164 | 14.28 | 1.884 | $12.888\pm 0.013$
*2m0718 | $\cdots$ | T5 | TW | 16.62 | 16.69 | -0.07 | $15.440\pm 0.017$
*2m0809 | $\cdots$ | L6 pec | 1 | 16.44 | 14.42 | 2.02 | $13.015\pm 0.006$
*2m0951 | $\cdots$ | L5 pec | TW | 17.098 | 15.28 | 1.82 | $14.554\pm 0.008$
*2m1551 | L4$\gamma$ |  $>$L5$\gamma$ | 3 | 16.32 | 14.31 | 2.009 | $13.243\pm 0.004$
*2m1647 | L7 | L9 pec(red) | 9 | 16.59 | 14.48 | 2.107 | $13.180\pm 0.005$
*2m1741 | $\cdots$ | L6-L8$\gamma$ | 3, 10 | 15.78 | 13.44 | 2.34 | $11.936\pm 0.003$
*2m2002 | L5$\beta$ | L5-7$\gamma$ | 3 | 15.32 | 13.42 | 1.90 | $12.161\pm 0.002$
*2m2117 | $\cdots$ | T0 | 7 | 15.6 | 14.15 | 1.45 | $13.014\pm 0.004$
*2m2154 | L4$\beta$ | L5$\gamma$ | 1,3 | 16.44 | 14.12 | 2.24 | $12.979\pm 0.004$
*2m2206+33 | $\cdots$ | T1.5 | 7 | 16.58 | 15.75 | 0.83 | $14.630\pm 0.011$
*2m2206-42 | L4$\gamma$ | L4$\gamma$ | 3 | 15.56 | 13.61 | 1.95 | $12.539\pm 0.003$
*2m2216 | $\cdots$ | T3 | 7 | 16.59 | 16.28 | 0.31 | $14.182\pm 0.007$
*2m2322 | L2$\gamma$ | L3$\gamma$ | 1,5 | 15.54 | 13.86 | 1.68 | $12.893\pm 0.003$
*2m2343 | $\cdots$ | L3-L6$\gamma$ | 1 | 16.57 | 14.19 | 2.38 | $12.835\pm 0.004$
Table 1: The survey consists of 26 brown dwarfs with spectral and/or kinematic
evidence of youth.

###  2.2 YMG Membership Probabilities of the Sample

For each target in our sample we assess the likelihood of kinematic membership
in nearby young moving groups. Proper motions, radial velocities and
parallaxes used in our kinematic analysis are shown in Table
LABEL:tab:kinematics. We use these as inputs for the BANYAN $\Sigma$ tool
(Gagne2018b), which assesses the probability that a source is a member of a
young moving group based on its observed kinematics. We additionally use the
LACEwING tool (Riedel2017) to check for consistency, but for clarity we only
present the BANYAN $\Sigma$ results. We consider the following moving groups:
TW Hydra (TWA, $10\pm 3~{}$Myr; Bell2015), $\beta$ Pictoris ($\beta$ Pic,
$22\pm 6~{}$Myr; Shkolnik2017), Columba (COL, $42^{+6}_{-4}~{}$Myr; Bell2015),
Tucana-Horologium (THA, $45\pm 4~{}$Myr; Bell2015), Carina (CAR, $45\pm
4~{}$Myr; Bell2015), Argus (ARG, 30–50 Myr; Torres2008), AB Doradus (ABDMG,
110–150 Myr; Luhman2007; Barenfeld2013) and Carina-Near (CARN, $200\pm
50~{}$Myr; Zuckerman2006). Field objects are estimated to have ages $>1$ Gyr.
We show the kinematic information and results from membership tools in Table
LABEL:tab:kinematics. We use the BANYAN $\Sigma$ results to make a final
assessment on each target using the following categories:

  * •

Bona fide member (BF): an object with full kinematic information (parallax,
proper motion, radial velocity), membership probability $>90\%$ and
spectroscopic signs of youth.

  * •

High-likelihood member (HLM): An object with membership probability $>90\%$
but does not have full kinematic information.

  * •

Candidate member (CAN): An object that shows greater probability of moving
group membership than being a field object, with missing kinematic
information.

  * •

Field object (FLD): An object that is kinematically eliminated from falling
into a known nearby moving group based on current kinematic information, or if
the spectrum looks unambiguously field age.

  * •

Young field object (YNG-FLD): An object that is kinematically eliminated from
falling into a known nearby moving group based on current kinematic
information but other observational features point to its youth.

  * •

Ambiguous object (AM): An object that requires updated astrometric precision
because it could either belong to more than one group or it cannot be
differentiated from the field.

As shown in Table LABEL:tab:kinematics, our sample is composed of three bona
fide members (2m0355, 2m1741 and 2m2322), seven high-likelihood members, eight
candidate members, two field objects, four young field objects and two
ambiguous members. We discuss the kinematic and spectral indicators of youth
for each object in detail in Section 5. In Figure 3 we plot the full sample on
a color-magnitude diagram.

Figure 3: Spitzer variability sample shown in a color-magnitude diagram for
objects with measured parallaxes (pink circles) and kinematic parallaxes
(purple diamonds), compared to the larger population of brown dwarfs. For
sources identified as significantly variable in Section 4, their amplitudes
are proportional to the partially transparent symbol size. $0.5\%$ and $1\%$
variability amplitudes are shown in the top left for reference.  We exclude
three sources from our variability analysis which are denoted by black
crosses. 2m0506 and 2m0447 have kinematics that are consistent with field
membership and 2m2117 has highly discrepant parallax measurements. Most of the
sample falls in their expected positions relative to the field brown dwarf
sample. 2m0447 is over-luminous compared to other field brown dwarfs, possibly
due to unresolved binarity. The T5 dwarf 2m0718 is the faintest and coldest
object in the sample.

##  3 Spitzer Observations and Data Reduction

We used the Infrared Array Camera (IRAC; Fazio2004) to observe our targets in
the Channel 1 ($3.6~{}\mu$m) band as part of the Cycle 14 Program: “The Young
and the Restless: Revealing the Turbulent, Cloudy Nature of Young Brown Dwarfs
and Exoplanets” (PID: 14128, PI: J Faherty). The observations were designed
following the recommendations for obtaining high precision photometry from the
Spitzer Science Center, and following the many examples of brown dwarf
variability monitoring in the literature (Metchev2015; Apai2017; Biller2018;
Vos2018; Vos2020; Allers2020; Zhou2020). Science images were obtained in
staring mode, using 12 s exposures for a total monitoring duration of 20 hr.
Science observations were preceded by a 30-minute dithered sequence to remove
the initial slew settling that occurs when acquiring a new target, and
followed by a 10-minute dithered sequence to capture the most accurate dark
image for each observation.

We follow the same data reduction outlined in Vos2020; we measured photometry
from the Basic Calibrated Data (BCD) images produced by the Spitzer Science
Center using pipeline version S19.2. The box_centroider.pro routine was used
to find the centroids of the target and reference stars of similar brightness
in the field of view. Aperture photometry was performed on the target and
reference stars using apertures with radii of $2.0-5.0$ pixels, in steps of
$0.25$. We choose the final aperture size that produces the lowest rms target
light curve. Outliers were identified and rejected from the raw light curves
using a $6\sigma$ clip.

Spitzer/IRAC photometry is known to exhibit a systematic effect due to intra-
pixel sensitivity variations, also known as the pixel phase effect. Intra-
pixel sensitivity variations combined with telescope pointing fluctuations
result in raw photometry that is highly correlated with the $x$ and $y$ sub-
pixel coordinates. We model the pixel phase effect using a cubic function of
the $x$ and $y$ coordinates (Knutson2008; Heinze2015; Metchev2015; Vos2020):

$f(x,y)=P_{0}+P_{1}x+P_{2}y+P_{3}xy+P_{4}x^{2}+P_{5}y^{2}+P_{6}x^{3}+P_{7}y^{3}+P_{8}x^{2}y+P_{9}xy^{2}$
(1)

where $f(x,y)$ represents the measured flux, $P_{i}$ are the fitted
coefficients, and $x$ and $y$ are the sub-pixel coordinates. We correct the
light curves of the target and reference stars using Equation 1, and find that
the fitted coefficients $P_{i}$ are similar for the target and reference
stars. We find that the correction significantly decreases the correlation
between the flux and pixel position for each observation. In Figure 4 we show
an example target light curve before and after correction using the cubic
function listed in Equation 1. This is the final step in producing corrected
light curves.

Figure 4: Raw light curve (top) and corrected light curve for 2m0030. The
combined astrophysical variability and pixel phase model is shown in blue and
the astrophysical variability model is shown in purple.

###  3.1 Determining variability parameters using MCMC

For objects showing variability (the process for identifying variability is
presented in Section 4), we use the Markov Chain Monte Carlo (MCMC) emcee
package (fm2013) to constrain the variability parameters such as rotation
period and variability amplitude. For the MCMC analysis we use 1000 walkers
with 8000 steps, and discard the initial 1000 steps as the burn-in sample. We
check for convergence by visually inspecting the resulting chains for each
parameter to ensure that they are consistent with Gaussian noise.
Additionally, we check that there is no significant difference between the
parameter constraints obtained from the first and second halves of each chain.
Based on these two checks we find that the model converges well for each
variable light curve. We use a sinusoidal model for all of our targets, except
for 2m0030 and 2m0642, whose light curves favor a two-term Fourier function.
We show our measured variability parameters along with their 1$\sigma$
uncertainties in Table 2. Best-fit variability models are overplotted in pink
in Figure 7.

We plot all of our variability detections on a color - magnitude diagram in
Figure 3, where the symbol area is proportional to the amplitude of
variability. This figure does not show any obvious trends between spectral
type and amplitudes and we observe a range of amplitudes throughout the L–T
spectral sequence. We present a more in-depth discussion of variability
amplitudes in Section 7.

##  4 Variability Analysis

Following the approach outlined in Vos2020, we use two independent methods to
identify variable signals in each light curve - a periodogram significance
test and the Bayesian Information Criterion (BIC). We briefly outline both
methods below.

Figure 5: Periodogram of the light curve (red) as well as $x$ and $y$ pixel
positions. For variable targets, we confirm that the periodicity is not
correlated with the periodicity of the $x$ and $y$ pixel coordinates.

###  4.1 Variability Detection with Periodogram Analysis

We calculate the Lomb-Scargle periodograms of our target and reference star
light curves (Scargle1982) to assess the significance of their brightness
fluctuations. For each observation we calculate the $95\%$ and $99\%$
significance thresholds by simulating non-variable light curves from our
observed reference stars. The simulated, non-variable light curves are created
by by randomly rearranging the indices of the reference star light curves,
which produces simulated light curves with Gaussian noise similar to our
observed light curves. We consider a target to be variable only if the
periodogram peak falls above the $99\%$ significance threshold. For variable
targets, we also perform a visual check to ensure that the periodicity of the
object is not correlated with the periodicity of the $x$ and $y$ pixel
coordinates. In Figure 5, we this comparison for the object 2m0718.

###  4.2 Variability Detection Using the Bayesian Information Criterion

We additionally use the BIC to identify significant variability in each
observation (Schwarz1978). The BIC is defined as

$\mathrm{BIC}=-2~{}\mathrm{ln}~{}\mathcal{L}_{\mathrm{max}}+k~{}\mathrm{ln}~{}N$
(2)

where $\mathcal{L}_{\mathrm{max}}$ is the maximum likelihood achievable by the
model, $k$ is the number of parameters in the model and $N$ is the number of
datapoints used in the fit (Schwarz1978).

For each observation, we calculate
$\Delta\mathrm{BIC}=\mathrm{BIC}_{flat}-\mathrm{BIC}_{\mathrm{sin}}$ to assess
whether the variable sinusoidal or non-variable flat model is favored by the
data. The BIC penalizes the sinusoidal model for having additional parameters
compared with the flat model. $\Delta\mathrm{BIC}>0$ indicates that the
sinusoidal model is favored and $\Delta\mathrm{BIC}<0$ indicates that the non-
variable, flat model is favored. A $|\Delta\mathrm{BIC}|$ value between 0 and
6 indicates that one model is positively favored over the other, a value
between 6 and 10 indicates that one model is strongly favored over the other
and values above 10 indicate that one model is very strongly favored over the
other (Schwarz1978). We conservatively consider light curves with
$\Delta\mathrm{BIC}>10$ to be a significant detection of variability for this
survey.

We additionally use the BIC to determine whether any of our light curves favor
a more complex variability model. For each light curve we use the BIC to test
both a sinusoidal model and a two-term Fourier model. Such models have been
used to describe more complex brown dwarf light curves in the literature
(Heinze2013; Yang2016; Vos2018), and may be driven by atmospheres with
multiple atmospheric features and/or evolving cloud structures. We find that a
two-term Fourier series is favored for two of our variable light curves - the
L7 dwarf 2m0030 and the L9 dwarf 2m0642.

###  4.3 Comparison of $\Delta$BIC and Periodogram Methods

In Figure 6 we compare the results of both variability identification methods.
Objects placed in the upper right of the figure are considered variable by
both methods, while objects in the lower left are classified as non-variable
using both methods. We find that in general both methods agree on which light
curves show significant variability as well as the relative significance of
their variability detections. One object, 2m2002, is found to be variable
using the periodogram method but not the $\Delta$BIC test. We conservatively
count it as non-variable for the remainder of the paper. We show our positive
variability detections in Figure 7, and non-detections in Figure 8. In total,
we detect variability in 15 objects, and do not detect variability in 11
objects.

Figure 6: Periodogram power plotted against $\Delta$BIC for all of our survey
targets. Variable objects are shown in pink with non-variable objects are
shown in grey. 17 has been added to the $\Delta$BIC values to include all
targets in the log plot. The dashed lines show the significance thresholds for
each test. Objects placed in the upper left of the figure are considered
variable by both methods, while objects in the lower left are classified as
non-variable using both methods. We find that both methods generally identify
the same variables, except for the object 2m2002, which is significantly
variable according the the periodogram power method but not the $\Delta$ BIC
test.

###  4.4 Modelling Red Noise Contributions Using Gaussian Processes

While our two variability identification methods are in agreement, we must
ensure that the detected variability is astrophysical in nature, and not due
to correlated, or “red”, noise. Red noise can mimic the signatures of
astrophysical variability.

Littlefair2017 presents a framework to determine between astrophysical
variability and red noise using Gaussian Processes (GPs). GPs have recently
been used in the exoplanet and stellar subfields of astronomy to model both
light curve systematics and rotational modulations (e.g. Gibson2012b;
Aigrain2016; Littlefair2017; Angus2017).

For this work, we compare two GP models computed using the celerite2 package
(Foreman-Mackey2017; Foreman-Mackey2018). The first exhibits variability
driven by red noise, as modelled by the Matern-3/2 kernel:

$k(\tau)=\sigma^{2}\left(1+\frac{\sqrt{3}\tau}{\rho}\right)exp\left(-\frac{\sqrt{3}\tau}{\rho}\right)$
(3)

where $\sigma$ and $\rho$ are the GP hyperparameters. The second model
exhibits rotational variability as described by a kernel that is a mixture of
two simple harmonic oscillator kernels given by:

$\displaystyle Q_{1}=\frac{1}{2}+Q_{0}+\delta Q$ (4)
$\displaystyle\omega_{1}=\frac{4\pi Q_{1}}{P\sqrt{4Q_{1}^{2}-1}}$
$\displaystyle S_{1}=\frac{\sigma^{2}}{(1+f)\omega_{1}Q_{1}}$

for the primary term, and

$\displaystyle Q_{2}=\frac{1}{2}+Q_{0}$ (5)
$\displaystyle\omega_{2}=\frac{8\pi Q_{1}}{P\sqrt{4Q_{1}^{2}-1}}$
$\displaystyle S_{2}=\frac{f\sigma^{2}}{(1+f)\omega_{2}Q_{2}}$

for the secondary term. Its hyperparameters are easier to interpret that for
the previous model. $\sigma$ represents the standard deviation of the process,
$P$ represents the primary period of variability, $Q_{0}$ represents the
quality factor for the secondary oscillation, $\delta Q$ is the difference
between quality factors for the first and second modes, and $f$ is the
fractional amplitude of the secondary mode compared to the primary.

For each of our variable light curves we use emcee (fm2013) to fit the GP
model and compare models using the $\Delta BIC$ framework described
previously. For each model we initialize 32 walkers and run for 5000 steps
after discarding a 2000 burn-in sample. Most light curves are fit moderately
well by both models, but the rotation model in particular has a lot of
parameters that are not constrained by the data, particularly for the longer
period rotators. We present the $\Delta BIC_{\mathrm{Mat3/2}-\mathrm{Rot}}$
for objects identified as variable in Section 4.3 in Table 2. We find that all
of the previously identified variables are preferred by the rotation model
over the red noise model, apart from 2m0326, the target identified with to
have a long period. The upward trend is fit by the both models equally well,
but since the Matern-3/2 kernel has fewer paramters it is favored. We conclude
that this target requires longer-term monitoring to identify the cause of its
upward flux trend, and leave it out of the statistical analysis.

Figure 7: Positive variability detections from our survey. The blue points
show 2.5 min cadence. $\Delta$BIC values show that a variable model is favored
in each case. A sinusoidal model is preferred for the majority of our
variables, and a two-term Fourier series is favored for 2m0030 and 2m0642.

Figure 8: Non-detections of variability from our survey. The blue points show
the light curve binned to 2.5 min cadence.

###  4.5 Simultaneous Fit for Intra-Pixel Sensitivity and Variability
Parameters

The pixel-phase effect can in principle be covariant with astrophysical
variability (e.g. Heinze2013). As described in Vos2020, we also perform a
simultaneous fit of the intra-pixel phase effect (Equation 1) and
astrophysical variability. We again use the emcee package (fm2013) to explore
the parameter space. We find that fitting these effects simultaneously
identifies the same variables and non-variables, and results in variability
parameters that are broadly consistent with those presented in Table 2.
However, for the simultaneous fit we consistently find a higher correlation
between the final light curve and the $x,y$ pixel positions compared to the
two-step fitting method; the Kendall $\tau$ correlation coefficient for the
simultaneous fit can be up to an order of magnitude higher than the two-step
method, and occasionally the correlation with pixel position increases
compared to the raw flux. For this reason we favor light curves corrected
using the two-step method in this paper, as the corrected photometry is less
likely to be affected by intra-pixel sensitivity.

Name | SpTa | [3.6] Amp | Period | Periodogram |  $\Delta$BICflat-sinb |  $\Delta BIC_{\mathrm{Mat3/2}-\mathrm{Rot}}$c
---|---|---|---|---|---|---
|  | (%) | (hr) | Power |  |
2m0001 | L4$\beta$ | $0.69\pm 0.04$ | $15.75\pm 0.37$ | 85.4 | 114.3 | 41.3
2m0030 | L7 | $1.52\pm 0.06$ | $4.22\pm 0.02$ | 182.3 | 572.2 | 2.4
2m0031 | L9 | $0.35\pm 0.03$ | $1.64\pm 0.01$ | 40.8 | 73.6 | 23.7
2m0153 | L2 | $0.48\pm 0.07$ | $17.63^{+1.13}_{-0.94}$ | 25.3 | 41.9 | 10.5
2m0326 | L5$\beta$ | $0.34\pm 0.04$ | $32.84^{+4.60}_{-4.62}$ | 40.8 | 73.6 | -4.4
2m0342 | L4$\gamma$ | $0.73\pm 0.07$ | $14.73^{+0.51}_{-0.49}$ | 52.7 | 111.8 | 43.5
2m0349 | L5$\gamma$ | $0.53\pm 0.09$ | $14.62^{+1.08}_{-1.00}$ | 16.8 | 23.4 | 30.1
2m0355 | L5$\gamma$ | $0.26\pm 0.02$ | $9.53\pm 0.19$ | 85.2 | 85.2 | 20.7
2m0447 | T2 | $<0.93$ |  $\cdots$ | 5.2 | -9.1 |  $\cdots$
2m0459 | L7: | $<0.47$ |  $\cdots$ | 4.7 | -7.1 |  $\cdots$
2m0506 | T4.5 | $0.58\pm 0.10$ | $13.77^{+0.82}_{-0.77}$ | 17.1 | 21.7 | 23.1
2m0642 | L9(red) | $2.16\pm 0.16$ | $10.11\pm 0.06$ | 421.1 | 1317.8 | 49.9
2m0718 | T5 | $2.14\pm 0.21$ | $1.080^{+0.004}_{-0.003}$ | 44.9 .2 | 47.5 |
2m0809 | L6p | $0.77\pm 0.06$ | $1.365\pm 0.004$ | 79.1 | 185.3 | 61.2
2m0951 | L5(pec) | $<0.92$ |  $\cdots$ | 2.7 | -12.4 |  $\cdots$
2m1551 | L4$\gamma$ | $<0.39$ |  $\cdots$ | 7.0 | 0.4 |  $\cdots$
2m1647 | L7 | $0.47\pm 0.06$ | $9.234^{+0.25}_{-0.23}$ | 32.0 | 52.7 | 2.8
2m1741 | L6-L8$\gamma$ | $0.35\pm 0.03$ | $15.00^{+0.71}_{-0.57}$ | 45.5 | 96.1 | 32.9
2m2002 | L5$\beta$ | $<0.28$ |  $\cdots$ | 12.1 | -1.3 |  $\cdots$
2m2117 | T0 | $<0.38$ |  $\cdots$ | 4.6 | -6.5 |  $\cdots$
2m2154 | L4$\beta$ | $<0.39$ |  $\cdots$ | 5.3 | -8.2 |  $\cdots$
2m2206+33 | T1.5 | $1.20\pm 0.13$ | $15.91^{+0.62}_{-0.61}$ | 38.4 | 67.8 | 38.5
2m2206-42 | L4$\beta$ | $<0.33$ |  $\cdots$ | 6.8 | 2.9 |  $\cdots$
2m2216 | T3 | $<0.97$ |  $\cdots$ | 9.6 | -15.6 |  $\cdots$
2m2322 | L2$\gamma$ | $<0.41$ |  $\cdots$ | 6.9 | -5.4 |  $\cdots$
2m2343 | L3-L6$\gamma$ | $<0.37$ |  $\cdots$ | 3.6 | -10.9 |  $\cdots$
Table 2: Variability parameters for objects in the survey.

###  4.6 Survey Sensitivity

We estimate the sensitivity of each observation to simulated sinusoidal
variability following the methods outlined in Vos2019; Vos2020. We inject
sinusoidal signals with a ranges of amplitudes and periods into light curves
with Guassian distributed noise similar to that of each target. The simulated
light curves have amplitudes of $0.05-0.5\%$, rotation periods of
$0.5-40~{}$hr and random phase shifts. We then analyze the simulated light
curves using the periodogram analysis described in Section 4 and calculate the
detection probability as the percentage of light curves with a given
variability amplitude and period that produces a periodogram power above the
significance threshold. Due to the range of magnitudes in our sample, the
sensitivity of each observation can change significantly.

For non-variable targets, we determine an upper limit on the variability
amplitude as the 90% contour boundary at 10 hours. Metchev2015 also set their
upper limits assuming 10 hr periods, so our upper limits can be directly
compared with their values. The left panel of Figure 9 shows the median
sensitivity for the survey. Contours represent the rate of detection for
signals of a given amplitude and period. The right panel of Figure 9 shows our
calculated median sensitivity of the Metchev2015 survey for comparison,
calculated in the same way. We note a slight improvement in the sensitivity of
our survey compared to the Metchev2015. Both surveys have similar observation
lengths ($\sim 20$ hr) and the median $3.6~{}\mu$m magnitude for both samples
is 13.1–13.2. The slightly higher sensitivity achieved by our survey is likely
due to the fact that the Metchev2015 survey swapped between the [$3.6~{}\mu$m]
and [$4.5~{}\mu$m] during each monitoring observation, or the differing data
reduction and intra-pixel sensitivity correction methods. It is also possible
that the amplitude of telescope motion for our observations was less than that
of the Metchev2015 survey; the telescope drift during our $\sim 20~{}$hr
observations was consistently below $0.5$ pixels for the observations
presented in this survey, resulting in low-level intra-pixel effects. The
sensitivity plot for each observation is used in Section 6 to estimate the
variability occurrence rates of young and field brown dwarfs.

Figure 9: Median survey sensitivity as a function of variability amplitude and
rotation rate for this survey (left) and the Metchev2015 survey. Contours show
the percentage of recovered signals for each amplitude and period. The
detection rates achieved by both variability surveys are very similar, and
thus can be compared robustly. For non-variable objects, we assign upper
limits on the variability amplitude as the 90% contour at 10 hr.

##  5 Youth and Variability Properties of Inidividual Objects

In this section we discuss the youth and variability properties of each object
in the sample individually.

2m0001  was first identified as a candidate member of the AB Doradus moving
group by Gagne2015c. They classify 2m0001 as a L4$\beta$ spectral type, with
an int-g gravity score (Allers2013). Faherty2016 classified 2m0001 as an
ambiguous member as it lacked a parallax and radial velocity. With updated
proper motion and parallax measurements from Gaia DR2 (Gaia2018), we find that
2m0001 is a high-likelihood member of the AB Doradus moving group (Table
LABEL:tab:kinematics). We note that the proper motions and a parallax
presented in Best2020, which are within $2-3\sigma$ of the Gaia DR2 values,
favor field over AB Doradus membership. Vos2019 monitored 2m0001 in their
$J$-band survey, but did not detect variability above $\sim 2\%$ during their
$3.75$ hr observation. In this work we detect significant $3.6~{}\mu$m
variability with an amplitude of $0.69\pm 0.04\%$ and a period of $15.75\pm
0.37$ hr. Its long period likely prohibits detection in shorter observations –
the sensitivity map for the light curve reported by Vos2019 shows that
sinusoidal signals with amplitudes of $>8-10\%$ would have been detected with
a period $\sim 15$ hr. Its long period is also consistent with other known
variables from the AB Doradus moving group such as WISE J0047+68 ($16.4\pm
0.2$ hr; Vos2018).

2m0030 was given an optical classification of L7 by Kirkpatrick2000, an
infared spectral type of L4-L6$\beta$ by Gagne2015c, and was assigned a fld-g
gravity class by Liu2016. Its kinematic measurements in Table
LABEL:tab:kinematics give it a high probability of membership in Argus
according to BANYAN $\Sigma$. We classify it as a high-likelihood member of
Argus. 2m0030 has been searched extensively for optical and infrared
variability from ground-based telescopes in the literature (Enoch2003;
Clarke2008; Radigan2014). Only Enoch2003 find a marginal detection of
variability in the $i$ band with a period of $~{}1.5$ hr. Clarke2008 and
Radigan2014 report light curves with $1-4\%$ precision and report no evidence
for variability. We detect variability in its light curve, with an amplitude
of $1.46\pm 0.06\%$ and a rotation period of $4.29\pm 0.02$ hr. The light
curve appears stable over 5 rotations, however the non-detections reported in
the literature to date suggest the possibility of an evolving-light curve.
2m0030 is one of two variables whose light curve favors a two-term Fourier fit
over a sinusoidal model.

2m0031/PSO J282.7576+59.5858 was given an L9 spectral type by Best2013, who
note that it is not matched well by any of their spectral templates. They
suggest that the spectrum may be a blended spectrum of multiple objects,
however the spectral indices do not classify it as a binary. Best2020 also
report a parallax of $71\pm 5$ mas. Using the kinematic measurements in Table
LABEL:tab:kinematics, BANYAN $\Sigma$ reports that 2m0031 is a high-likelihood
member (98% probability) of the Carina-Near moving group. 2m0031 has not been
monitored for variability in the past. Our new Spitzer light curve reveals
variability with an amplitude of $0.35\pm 0.03\%$ with a rotation period of
$1.64\pm 0.01$ hr. The light curve shows evidence of light curve evolution
during the 20 hr observations, with the amplitude ranging from $\sim 0.3-1\%$.
The periodogram shows two narrow peaks at $\sim 1.6$ hr, which is suggestive
of differential rotation (Apai2017). This target adds to the small but growing
number of L/T transition brown dwarfs that show light curve evolution on
rotational timescales (e.g. Apai2017; Apai2021).

2m0153 was classified as a L3$\beta$ with a vl-g gravity class by Gagne2015c,
who also identified it as a likely member of the Tucana-Horologium moving
group. Using the kinematics presented in Table LABEL:tab:kinematics, both
BANYAN $\Sigma$ and LACEwING identify it as a likely Tucana-Horologium member,
however parallax and radial velocity measurement are needed to confirm this
object as a member. We classify 2m0153 as a candidate member of Tucana-
Horologium. Our Spitzer light curve shows significant variability with a
$0.48\pm 0.07\%$ amplitude and a period of $17.63^{+1.13}_{-0.94}$ hr.

2m0326 has been identified as a highly probable member of AB Doradus by
Gagne2013 who also assign a L5 $\beta/\gamma$ spectral type, and note that the
spectrum is very red for its spectral type. Cruz2007 suggest that 2m0326 is
likely younger than $500$ Myr based on the strength of its Li absorption. Our
updated kinematic analysis which uses new proper motions from (Marocco2021)
agrees with these findings. We classify 2m0326 as a high-likelihood member of
AB Doradus but a radial velocity is needed to confirm membership. Our Spitzer
monitoring shows a slowly increasing brightness during the observations. The
variability is highly significant, with a periodogram power of $\sim 40$ and a
$\Delta$BIC of $\sim 74$. Our MCMC variability fitting finds a period of
$33\pm 5~{}$hr, however this is extremely uncertain since we did not cover a
full period. The GP analysis presented in Section 4.4 finds that the signal
may be explained by correlated red noise in the data. Whether the variability
is rotationally modulated, related to a longer-term astrophysical process or
caused by systematics is an interesting question, which additional long-term
monitoring will answer.

2m0342 is an L4$\gamma$ with strong membership probability in the Tucana-
Horologium association (Gagne2015c) and Faherty2016 classify it as a high-
likelihood member. Using proper motions from Gagne2015a, we recover similar
results, finding that both BANYAN $\Sigma$ favor membership in Tucana-
Horologium. We classify 2m0342 as a high-likelihood member of Tucana-
Horologium. The light curve of 2m0342 displays significant variability over
the observation. With an estimated period of $14.7\pm 0.5$ hr, we cover just
over one rotation period during our Spitzer observation.

2m0349 is one of three likely young brown dwarfs presented in this work for
the first time. As discussed in Section 2.1, 2m0349 is an L5$\gamma$, with
prominent alkali absorption lines that are suggestive of low-gravity. We use
BANYAN $\Sigma$ to assess membership using measured kinematic information.
Based on sky positions and proper motions from (Marocco2021) (shown in Table
LABEL:tab:kinematics), we find a combined probability of $75\%$ that 2m0349 is
a member of either AB Doradus or $\beta$ Pictoris. Including the radial
velocity measured from a FIRE/Echelle spectrum (Gagné et al. in prep) reduces
the probabilities of AB Doradus and $\beta$ Pictoris membership to $25\%$ and
$19\%$ respectively. While the kinematic analysis doesn’t conclusively
identify moving group membership, the low-gravity signatures in its spectrum
lend support to its likely youth. We consider 2m0349 as an ambiguous member
for the remaining analysis. We detect significant variability in our Spitzer
light curve, with an amplitude of $0.53\pm 0.09\%$ and a period of
$14.62^{+1.1}_{-1.0}$ hr.

2m0355 is one of the prototypical isolated exoplanet analogs. It was
identified as a high-likelihood member of the AB Doradus moving group by
Faherty2013, and confirmed by follow-up measurements and analysis by Liu2016;
Faherty2016. Cruz2009 assign an L5$\gamma$ spectral type. Faherty2013 noted
its extremely red near-infrared colors and underluminosity compared to field
brown dwarfs. These peculiarities were first noted by Metchev2006 for the
young companion HD203030 B (Miles-Paez2017), and are shared by directly-imaged
exoplanets such as 2M1207 b and HR8799 bcd. Our kinematic analysis, which
makes use of Gaia proper motions and parallax confirms that it is a member of
AB Doradus using BANYAN $\Sigma$. We thus confirm 2m0355 as a bona fide member
of AB Doradus. We measure significant variability in its $3.6~{}\mu$m light
curve. We measure a variability amplitude of $0.26\pm 0.02\%$ and a period of
$9.53\pm 0.19$ hr. 2m0355 joins W0047+68 and 2M2244+20 (Vos2018) as variable
bona fide members of the AB Doradus moving group. We specifically focus on the
variability properties of AB Doradus members in Section 9.

2m0447/PSO J071.8769–12.2713 was discovered by Best2015 as a T2 candidate of
the $\beta$ Pictoris moving group. Best2020 measure new proper motions and
parallax for this target which determine that it is likely a field object
using BANYAN $\Sigma$. Proper motions from CatWISE2020 (Marocco2021) also
agree that 2m0447 is likely a field object. We conclude that 2m0447 is likely
a field object and leave it out of the statistical analysis of variability in
young objects in Section 6. In Figure 3 we plot our sample on a color-
magnitude diagram and find that 2m0447 is over-luminous compared to the field
population and speculate that this may be due to an unresolved binary
companion, as suggested by Best2020. Vos2019 monitored this object for
variability in ground-based $J$-band observations, detecting significant
(4.5%) variability during a three hour observation. Our Spitzer light curve
does not display significant variability and we place an upper limit of
$0.9\%$. The variability detection in Vos2019 followed by a non-detection in
this work suggests that 2m0447 may be a variable object whose light curve
evolves over time. Such behavior has been observed in a number of variable L/T
transition brown dwarfs (e.g. Apai2017). Another possibility is that 2m0447
exhibits variability that is only apparent near-IR wavelengths, since
different wavelengths probe different atmospheric layers. .

2m0459 was discovered and classified as an L7 object with very red colors by
Schneider2017. They note spectral features that are consistent with low
surface gravity, but its kinematics did not match well with a young moving
group. Our updated kinematics, which include proper motions from CatWISE2020
(Marocco2021) and a radial velocity from Gagné et al. (in prep) indicates that
2m0459 is a high-likelihood member of Argus. Its red near-IR colors, spectral
signatures of youth and spectral type are remarkably similar to other young,
high-amplitude variables discovered to date such as PSO J318.5$-$22
(Biller2015; Biller2018) VHS 1256 B (Bowler2020; Zhou2020) and WISE 0047+68
(Vos2018; Lew2016). We do not detect variability in the Spitzer light curve of
2m0459, and place an upper limit of $<0.47\%$ on the peak-to-peak variability
amplitude.

2m0506/PSO J076.7092+52.6087 was discovered and classified as a T4.5 spectral
type by Best2015. Based on their measured proper motions, they determined that
2m0506 was a likely member of Argus. We check for moving group membership
using new proper motions and parallax from Best2020, which favor field
membership with a $69\%$ probability. Since this object has a T4.5 spectral
type, we cannot examine the spectrum to look for signatures of youth since
such signatures are not yet known. We classify this object as a field object
and do not include it in the following analysis sections. We detect
variability in 2m0506 during our Spitzer observation and measure an amplitude
of $0.6\pm 0.1\%$ and a rotation period of $13.8\pm 0.8$ hr. Although its
current kinematics favor field membership, its long rotation period is
indicative of youth.

2m0642 was discovered by Mace2013 as an “extremely red” object at the L/T
transition. Best2015 later assign a L9 spectral type based on its $J$-band
profile and the depth of its $2.4~{}\mu$m water absorption feature. Its
extremely red colors strongly suggest the presence of large amounts of dusty
photospheric condensates (Best2015). Best2015 and Gagne2014 both reported that
2m0642 is a likely AB Doradus member based on its proper motions and
positions. Our analysis, which uses proper motions and a new parallax from
Kirkpatrick2020, reports a $86\%$ probability of AB Doradus membership. We
thus consider it a candidate member of AB Doradus. With a measured variability
amplitude of $2.16\pm 0.16\%$, 2m0642 shows the highest [$3.6~{}\mu$m]
amplitude of any L dwarf measured to date, and is comparable to large
amplitude [$4.5~{}\mu$m] variables such as PSO J318.5$-$22 (Biller2018). With
a $\sim 10$ hr period, we observe two full rotations during our Spitzer
monitoring. During those two rotations we see significant evolution in the
light curve, including the emergence of a $\sim 1\%$ dip with a $1$ hr
duration in the light curve at 15 hr. This dip may be due to the emergence of
a Great Red Spot analog (Apai2017), or even a transiting satellite
(Tamburo2019; Limbach2021). Long-term follow-up observations could
differentiate the cause of such dips in the light curve of 2m0642.

2m0718 Using the positions, proper motions and radial velocity shown in Table
LABEL:tab:kinematics, BANYAN $\Sigma$ reports a membership probability of
$85\%$ for the $\beta$ Pictoris moving group. Assuming it is a member of the
$\beta$ Pic moving group, we estimate a mass of $2.6\pm
0.3~{}M_{\mathrm{Jup}}$ using the SED analysis described in Section 8.1.
Coincidentally, with a period of $\sim 1.08$ hr, 2m0718 joins a T7 field brown
dwarf, 2MASS J03480772$-$6022270 (hereafter 2MASS J0348$-$60), as the fastest
rotating brown dwarfs detected to date (Tannock2021). 2MASS J0348$-$60 is
estimated to be spinning at $\sim 45\%$ of its break-up velocity. However,
2m0718 is likely young, and we have classified it as a candidate $\beta$
Pictoris member. Using our estimated mass and radius from Table
LABEL:tab:SEDs, we find that the break-up velocity for 2m0718 is $\sim 2.22$
hr. Thus, if 2m0718 is indeed a $\beta$ Pictoris member it would be spinning
faster than its theoretical break-up velocity. It is further possible that he
observed variability pattern is produced by more than one dominant spatial
feature, as for example seen with the regular positioning of near-equatorial
hot spots on Jupiter (dePater2016). Follow-up high resolution spectroscopy
will allow us to determine whether this object truly rotates so rapidly.
Securing a parallax for 2m0718 and confirming its membership will reveal
whether it is a very low-mass object that is potentially in the process of
breaking up, or whether it is an older object rotating at a more moderate pace
like 2MASS J0348$-$60.

2m0809 was identified by Gagne2015c as an L6 pec (red) brown dwarf, and a
candidate member of the Argus moving group. 2m0809 is classified with an int-g
gravity index following the Allers2013 scheme. Best2020 provide a new parallax
measurement of $42.4\pm 3.6$ mas which differs substantially from the
predicted parallax of $\sim 65$ mas assuming membership in Argus (Gagne2015c).
Using this new parallax, BANYAN $\Sigma$ reports 2m0809 as a field object. Due
to its low-gravity spectral signatures, we consider 2m0809 to be a young field
object. The Spitzer light curve of 2m0809 shows significant variability with
an amplitude of $0.77\pm 0.06\%$ and a period of $1.365\pm 0.004$ hr.

2m0951 is one of the new brown dwarfs presented in this paper. As discussed in
Section 2.1, we assign 2m0951 an L5(pec) spectral type. Using proper motions
from CatWISE2020 (Marocco2021), and a radial velocity from Gagné et al. (in
prep), our kinematics analysis classifies it as a candidate member of Argus,
with a $87\%$ BANYAN $\Sigma$ probability. A parallax measurement is still
necessary to confirm membership. We do not detect variability in the Spitzer
light curve of 2m0951. We place an upper limit of $0.9\%$ on its amplitude.

2m1551 shows clear signatures of youth in its spectrum. It has been classified
as L4$\gamma$ by Faherty2013 and L4 vl-g by Allers2013. However, despite its
clear spectral signatures of youth, 2m1551 is not an obvious member of any
known moving group. Faherty2016 classified it as having ambiguous membership
in young moving groups. Using an updated parallax from Gaia2018 and an updated
radial velocity from Gagné et al. in prep, the BANYAN $\Sigma$ tool reports
that 2m1551 is likely a field object. Due to its strong youth indicators, we
consider that 2m1551 is a young field object. The Spitzer light curve of
2m1551 does not show statistically significant variability during our 20 hour
observation. We place an upper limit of $0.39\%$ on the variability amplitude
for periods of $\sim 10$ hr.

2m1647 was discovered and classified as an L9 (pec) brown dwarf with a very
red spectrum by Kirkpatrick2011. They report a parallax of $116\pm 29$ mas.
More recently, Best2020 report a parallax of $42.9\pm 2.1$ mas, which is not
in agreement with the original parallax measurement from Kirkpatrick2011.
Kirkpatrick2011 identifies 2m1647 as a candidate member of $\beta$ Pictoris or
Argus, with a low probability of field membership. However, the updated
parallax from Best2020 results in a 66% probability that 2m1647 is a field
object, and $34\%$ that it is a Carina-Near member. We classify 2m1647 as a
young field object. We detect significant variability in our Spitzer light
curve, and estimate a variability amplitude of $0.5\pm 0.1\%$ with a period of
$9.2\pm 0.3$ hr.

2m1741 was discovered and classified as L7 (pec) by Schneider2014. They also
identified 2m1741 as a candidate member of either the AB Doradus or the
$\beta$ Pictoris moving group. Faherty2016 classify it as a L6–L8 $\gamma$
based on its near-infrared spectrum. Recently, Kirkpatrick2020 provide a new
parallax measurement, and report a high probability that 2m1741 is a member of
AB Doradus. Our updated analysis, which uses the parallax from Kirkpatrick2020
and a new radial velocity from Gagné et al. (in prep) confirms 2m1741 as a
member of the AB Doradus moving group, with probability of 99% using BANYAN
$\Sigma$. This agrees well with the signatures of youth in its spectrum
(Schneider2014; Faherty2016). We detect significant variability in the Spitzer
light curve of 2m1741, with an amplitude of $0.35\pm 0.03\%$ and a period of
$15.0^{+0.7}_{-0.6}$ hr. Vos2019 previously monitored 2m1741 in $J$-band and
did not detect significant variability during a $\sim 2.5$ hr observation.
However, the sensitivity plot from this observation shows a very small
probability ($\sim 5\%$) of detecting variability with a period of $15$ hr.
Its rotation period matches well with the other late-L AB Doradus bona fide
members with periods of $9-16$ hr (W0047+68, 2M2244+20; Vos2018 and 2m0355;
this work), suggesting a common angular momentum history for these objects. We
discuss 2m1741 in context with other AB Doradus members in Section 9.

2m2002 was classified as an L5-L7$\gamma$ by Gagne2015c, but multiple studies
have been unable to place it in a young moving group (Gagne2015c; Faherty2016;
Vos2019). Our updated kinematics, which includes a new parallax from Gaia DR2
(Gaia2018) similarly do not place it within a known moving group. However, its
spectrum displays clear signature of youth, so we classify it as a young field
object. Vos2019 detect $J$-band variability in this object using UKIRT/WFCAM,
but the 4 hour observation did not yield periodicity. Its Spitzer light curve
does not show significant variability, and we place an upper limit of $0.3\%$
on $3.6~{}\mu$m. Follow-up observations may reveal whether 2m2002 shows
periods of variability followed by ‘quiescent’ phases, such as those objects
reported by Apai2017.

2m2117/PSO J319.3102-29.6682 was discovered and classified as a T0 dwarf by
Best2015. More recently, Best2020 measured a parallax of $52.4\pm 6.8$ mas for
this object, and report a BANYAN $\Sigma$ probability of 99% that is a member
of the $\beta$ Pictoris moving group. However, Kirkpatrick2020 report a
parallax of $76.1\pm 3.5$ mas, which gives a 99% probability that 2m2117 is a
field object. We do not include 2m2117 as a young object in the statistical
analysis of variability Section 6 since its membership is unclear. However, in
Section 8.1 we estimate fundamental parameters for this object for both cases
for completeness. We do not detect significant variability in the Spitzer
light curve of 2m2117, and place an upper limit of $0.38\%$ on its amplitude.

2m2154 was identified as an L5$\beta/\gamma$ candidate member of the Argus
moving group by Gagne2014c; Gagne2015c. Vos2019 later identified it as a
likely member of the Carina-Near moving group. Our updated kinematic analysis
finds that BANYAN $\Sigma$ favors Carina-Near membership, with a membership
probability of $80\%$. We thus classify 2m2154 as a candidate Carina-Near
member. 2m2154 was monitored for $J$-band variability by Vos2019, who find no
evidence for variability during a $\sim 4$ hr observation. Similarly, our
Spitzer light curve does not display any variability during our 20 hour
observation, and we place an upper limit of $0.4\%$ on the variability
amplitude of this object.

2m2206+33/PSO J331.6058+33.0207 was discovered and identified as a T1.5
candidate Argus member by Best2015 based on its kinematics. Our analysis,
which uses the proper motions from Best2015 finds that BANYAN $\Sigma$ favors
Argus ($68\%$). We note that CatWISE2020 (Marocco2021) report discrepant
proper motions for this target which result in a higher probability of field
membership. For this work we classify 2m2206+33 as a candidate Argus member. A
future radial velocity and parallax measurement will help to clarify
membership. We detect significant variability in 2m2206+33 with an amplitude
of $1.2\pm 0.1\%$ and a period of $15.9\pm 0.6$ hr. This is one of the longest
rotation periods measured for any T dwarf, and adds further support that
2m2206+33 may be young.

2m2206-42 was discovered and classified as an L2 dwarf by Kirkpatrick2000.
Kirkpatrick2008 detected lithium in its spectrum, and Gagne2013 report it as a
candidate AB Doradus member. Our updated kinematics analysis, which include
new proper motions and a parallax from Gaia DR2 (Gaia2018), confirms 2m2206-42
as a high-likelihood member of AB Doradus. A radial velocity is needed to
confirm its membership. We do not detect variability in our Spitzer light
curve of 2m2206-42, and place an upper limit of $0.33\%$ on the variability
amplitude.

2m2216/PSO J334.1+19 was discovered and identified as a T3 dwarf by Best2015,
who identify it as a possible comoving companion to the M4 star LSPM
J2216+1952. They find that both components are candidate members of $\beta$
Pictoris. Our independent kinematics analysis, which uses proper motions from
CatWISE2020 (Marocco2021), finds a slightly higher BANYAN $\Sigma$ $\beta$
Pictoris membership probability (85%) and we classify it as a candidate
$\beta$ Pictoris member. We also repeat the BANYAN $\Sigma$ analysis for the
potential comoving primary LSPM J2216+1952 using new Gaia DR2 (Gaia2018)
proper motions and parallaxes, which reports a $99\%$ probability that it is a
field star. Thus, we conclude that 2m2216 is a candidate $\beta$ Pictoris
member, but is likely not a companion to the star LSPM J2216+1952. We do not
detect variability in our Spitzer light curve of 2m2216, and place an upper
limit of $1\%$ on the variability amplitude.

2m2322 is an L2$\gamma$ (Gagne2013) that was identified as a possible comoving
companion the the M5 star 2MASS J23225240-6151114. Faherty2016 classify this
object as a high-likelihood member of Tucana Horologium. Our updated kinematic
analysis combines new proper motions and a parallax from Gaia DR2 (Gaia2018)
with a radial velocity measurement from Faherty2016 to confirm it as a member
of the Tucana Horologium moving group. The primary, which lacks a radial
velocity, shows a $99\%$ membership in Tucana Horologium also. Thus we confirm
that 2m2322 is a member of the Tucana Horologium moving group and consider it
highly likely that it is comoving with the M5 star 2MASS J23225240–6151114.
This system will be the subject of an upcoming paper (Faherty et al. in prep).
We do not detect variability in our Spitzer monitoring of 2m2322 and place a
limit of $0.4\%$ on the maximum variability amplitude for this object.

2m2343 was discovered and classified as an L3-L6$\gamma$ dwarf, with a
spectrum very similar to 2m0355. Faherty2016 identify it as an ambiguous
moving group member, which shows probability of membership in the AB Doradus
and Tucana Horologium moving groups. We use proper motions from CatWISE2020
(Marocco2021) in our updated kinematic analysis. BANYAN $\Sigma$ shows
moderately high probability of membership in AB Doradus. We classify it as a
candidate AB Doradus member. We do not detect significant variability in our
Spitzer light curve of 2m2343, and place an upper limit of $0.4\%$ on its
variability amplitude.

##  6 The Variability Occurrence Rate of Young Brown Dwarfs in the Mid-
infrared

Our variability detections and detection limits can be used to estimate the
fraction of young brown dwarfs that exhibit photometric variability (Vos2019).
The occurrence rate analysis presented in this section is described by
Lafreniere2007; Bonavita2013, and summarized below.

###  6.1 Statistical formalism

We consider variability monitoring observations of $N$ sources enumerated by
$j=1....N$. $f$ is the fraction of objects that exhibit variability with
amplitude and rotation period in the interval
$[a_{\mathrm{min}},a_{\mathrm{max}}]\cap[r_{\mathrm{min}},r_{\mathrm{max}}]$.
For this work we consider amplitude ranges of $0.05-3\%$ and rotation periods
of $0.5-40$ hrs. $p_{j}$ is the probability that such variability would be
detected from each observation. To calculate $p_{j}$ we use the variability
sensitivity maps, $g(a,r)$, shown in Figure 9. We integrate over $g(a,r)$ the
considered range in amplitude and period to obtain the probability of
detecting a variable signal with $a_{\mathrm{min}}<a<a_{\mathrm{max}}$ and
$r_{\mathrm{min}}<r<r_{\mathrm{max}}$ for target $j$:

$p_{j}=\int_{a_{\mathrm{min}}}^{a_{\mathrm{max}}}\int_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}g(a,r)~{}da~{}dr$
(6)

Using this notation, the probability of detecting variability in target $j$ is
the fraction of objects displaying variability multiplied by the probability
that such variability would be detected, or $fp_{j}$. The probability of not
detecting variability is $(1-fp_{j})$. We denote $d_{j}$ the detections made
by the observations such that $d_{j}=1$ for a positive variability detection
for target $j$ and $d_{j}=0$ for a non-detection. Then the probability of
observing a detection or non-detection for a given $f$ is:

$L\big{(}{d_{j}}|f\big{)}=\prod\limits_{j=1}^{N}(1-fp_{j})^{1-d_{j}}(fp_{j})^{d_{j}}$
(7)

According to Bayes’ theorem, the a priori probability density, or prior
distribution $p(f)$, and the likelihood function $L$, can be used to calculate
the posterior distribution $p(f|{d_{j}})$, the probability density updated in
light of the data:

$p\big{(}f|{d_{j}}\big{)}=\frac{L\big{(}{d_{j}}|f\big{)}p(f)}{\int_{0}^{1}L\big{(}{d_{j}}|f\big{)}p(f)\mathrm{d}f}$
(8)

Since there are no independent variability occurrence rate estimates
available, we use an non-informative Jeffrey’s prior. As outlined in Appendix
A, the suitable Jeffrey’s prior is given as:

$J(f)=\sqrt{\sum_{j}\frac{p_{j}}{f(1-fp_{j})}}$ (9)

We calculate the posterior probability distribution, $p(f|{d_{j}})$, which
represents the fraction of variable objects, or variability occurrence rate of
brown dwarfs. We calculate credible intervals for the variability occurrence
rate, $f$, for a given level of credibility, $\alpha$. We follow the method
outlined by Kraft1991 to calculate the confidence intervals
$[f_{\mathrm{min}},f_{\mathrm{max}}]$, i.e. the upper and lower boundaries for
the variability occurrence rate. We calculate the confidence interval by
solving for $[f_{\mathrm{min}},f_{\mathrm{max}}]$ in the equations:

$\alpha=\int_{f_{\mathrm{min}}}^{f_{\mathrm{max}}}p(f|{d_{j}})df$ (10)

and

$p(f_{\mathrm{min}}|d_{j})=p(f_{\mathrm{max}}|d_{j}).$ (11)

For the case of a one-sided distribution, the upper or lower bound of the
credible interval is found by solving:

$\alpha=\int_{f_{\mathrm{min}}}^{1}p(f|{d_{j}})df~{}\mathrm{or}~{}\alpha=\int_{0}^{f_{\mathrm{max}}}p(f|{d_{j}})df$
(12)

We determine the $68\%$ and $95\%$ confidence intervals in this paper.

###  6.2 The Variability Occurrence Rates of L and T Spectral Type Objects

Figure 10: Posterior distribution functions of the variability occurrence rate
of field brown dwarfs (blue) from Metchev2015 and young, low-gravity objects
presented in this work. The shaded areas show the variability occurrence rate
limits for the $95\%$ and $68\%$ confidence levels. Top panels show the
occurrence rates for L dwarfs, the middle panel shows occurrence rates for the
L/T transition and the bottom panel shows the variability occurrence rates for
T dwarfs.

We calculate the posterior distribution functions of the variability
occurrence rates for the low-gravity and field population of brown dwarfs
using the low-gravity sample presented here and the field objects from
Metchev2015. The Metchev2015 sample included eight low-gravity brown dwarfs,
which were added to the low-gravity survey sample presented in this paper and
excluded from the field sample. We discard the magnetically active brown dwarf
2MASSW J0036159+182110 because it was a known variable prior to the
Metchev2015 survey, as well as known binaries. We also discard three targets
from the low-gravity sample presented in this paper: 2m0447, 2m0506 and
2m2117, because their youth and kinematic membership are uncertain, as
discussed in Section 5. We calculate the posterior distribution functions of
the variability occurrence rate for L dwarfs (SpT $<$L9), L/T transition
objects (L9 $\leq$ SpT $\leq$T3.5) and T dwarfs (SpT $>$T3.5). We show the
variability occurrence rates of the young (pink) and field (blue) samples in
Figure 10 and discuss each spectral type bin below.

####  6.2.1 L Dwarf Variability Occurrence Rates

For field L dwarfs, we find a variability occurrence rate of $83-100\%$ (68%
confidence). We find a variability occurrence rate of $81-100\%$ for the low-
gravity L dwarfs, which is consistent with that of the field gravity L dwarfs.
Thus we find no evidence for an enhancement in variability occurrence rate for
the low-gravity L dwarfs. This is in contrast to Vos2019 who find a higher
$J$-band occurrence rate for the low-gravity objects. This apparent
discrepancy could be due to the different wavelength regimes and/or different
sensitivities of each survey. Vos2019 surveyed 30 low-gravity objects using
ground-based $J$-band observations. The $J$-band probes deeper pressure levels
than the Spitzer [$3.6~{}\mu$m] band used in this survey (Buenzli2012). The
$J$-band flux arises from pressures of $\sim 10$ bar while the Spitzer
[$3.6~{}\mu$m] band probes pressures of $\sim 1$ bar. Any difference in the
observed variability properties between these wavelength regimes may suggest
that they are sensitive to different variability mechanisms that take place at
different altitudes. Simultaneous spectroscopic variability monitoring ranging
from the near-infrared to the mid-infrared is necessary to test this
possibility. However, another significant difference is the increased
sensitivity achieved by Spitzer surveys (Metchev2015, this work) compared to
ground-based searches (Radigan2014; Wilson2014; Vos2019; Eriksson2019). The
Spitzer surveys on average probe variability amplitudes that are a factor of
$\sim 10$ smaller than the ground-based surveys mentioned above, and thus
offer a more complete picture of low-level variability. Ground-based surveys
thus may probe the variability occurrence rates of high-amplitude variables,
while the Spitzer surveys may probe a more complete picture of variability at
low to high amplitudes.

####  6.2.2 L/T Transition Variability Occurrence Rates

At the L/T transition (L9–T3.5), field dwarfs show a lower variability
occurrence rate compared to field L dwarf behavior, with a range of $22-64\%$
(68% confidence interval). Metchev2015 report similar findings, using
virtually the same sample (Metchev2015, combined young and field brown dwarfs
for their analysis whereas we removed low-gravity objects from the field dwarf
sample). As discussed by Metchev2015, this is in apparent contrast with the
ground-based $J$-band survey reported by Radigan2014, who find an enhancement
in variability occurrence rate at the L/T transition compared to earlier Ls
and later Ts. The reasons for this difference may be the same reasons outlined
above — the different wavelength regimes and sensitivities for Spitzer and
ground-based near-infrared variability surveys. In contrast to the field
dwarfs, young L/T transition objects maintain their high variability rate in
the range $85-100\%$. These results suggest that there may be an enhanced
variability occurrence rate of low-gravity L/T transition objects compared to
their field brown dwarf counterparts. Interestingly, this tentative
enhancement may also explain the apparent disagreement between the L/T
transition variability occurrence rates between Radigan2014 and Metchev2015
discussed above. Radigan2014 find high-amplitude variability in 4/16 of their
L9–T3.5 sample, resulting in a high occurrence rate of high-amplitude
variables. However, two of these detections, s0136 and 2m2139, have since been
identified as planetary-mass members of the Carina-Near moving group
(Gagne2017; ZJZhang2021). Thus, a reanalysis of the Radigan2014 variability
occurrence rates with updated youth information would be beneficial for
evaluating whether field L/T transition objects have an enhanced variability
rate.

####  6.2.3 T Dwarf Variability Occurrence Rates

For field T dwarfs, the maximum likelihood value drops slightly (46% to 18%),
but the confidence intervals ($3-50\%$) are generally consistent with the L/T
transition objects. Metchev2015 explain the lower rate of variability for
field T dwarfs as due to lower sensitivity to this sample, however since our
calculations take into account the sensitivity of the sample, we do not
believe that the low rate is due to lower sensitivity. For young T dwarfs, our
sample is composed of only two T dwarfs – 2m0718 (variable in this work) and
Ross 458C, which is a known variable (Manjavacas2019), but was observed to be
non-variable by Metchev2015, so is classed as non-variable for our analysis.
The resulting posterior distribution of the variability occurrence rate is
wide, and we find that $>44\%$ of young objects with spectral types greater
than T3.5 are likely to be variable at $68\%$ confidence, $>12\%$ at $95\%$
confidence. With the current sample we cannot robustly compare between the
variability occurrence rates of high and low-gravity objects at spectral types
$>T$4, but as more young T dwarfs are confirmed this will be an interesting
question.

##  7 The Variability Amplitudes of Young Brown Dwarfs at 3.6 micron

Figure 11: [$3.6~{}\mu$m] variability amplitudes of young objects (pink
circles) and field objects (blue triangles) as a function of spectral type.
Sources from this study and those compiled in Vos2020_cat. As is the case for
the field brown dwarfs (Metchev2015), the maximum variability amplitude of the
young objects steadily increases with spectral type from early L to the L/T
transition. It is also notable that the highest amplitude for each bin is
measured for a young object for spectral types L0–L9. Thus, it seems that the
maximum variability amplitudes of young brown dwarfs may be enhanced compared
to field brown dwarfs.

A number of studies have suggested that the variability amplitudes of low-
gravity objects may be enhanced compared to the field brown dwarfs.
(Biller2015; Metchev2015; Vos2020). Clouds in low-gravity atmospheres tend to
form at lower pressures/higher altitudes (Marley2012), and thus higher
variability amplitudes might be expected due to a higher photometric contrast
between regions of thick and thin clouds.

We combine the new variable objects detected in this survey with brown dwarfs
previously monitored for [$3.6~{}\mu$m] variability with Spitzer (compiled in
Vos2020_cat) to investigate whether [$3.6~{}\mu$m] amplitudes are enhanced for
low-gravity brown dwarfs. The T2.5 object 2m2139 has been added to the young
sample since it is a high probability member of the Carina-Near moving group
(ZJZhang2021). Figure 11 shows the [$3.6~{}\mu$m] variability amplitudes of
the young (pink) and field (blue) L2–T5 brown dwarfs. As is the case for the
field brown dwarfs (Metchev2015), the maximum variability amplitude of the
young objects steadily increases with spectral type from early L to the L/T
transition. The range of measured amplitudes in each spectral bin are
generally consistent between the low-gravity and field dwarf populations.
However, it is notable that the highest amplitude for each bin is measured for
a young object for spectral types L0–L9. Thus, it seems that the maximum
variability amplitudes of young brown dwarfs are enhanced compared to field
brown dwarfs. These results suggest that variability amplitudes may be
enhanced for low-gravity objects, but that secondary effects may affect the
observed amplitudes. For example, the viewing inclination (e.g. Vos2017)
likely reduce the observed variability amplitudes for objects viewed closer to
pole-on, and may be responsible for the observed spread in amplitudes at each
spectral type bin. Thus, the maximum variability amplitudes observed in low-
gravity brown dwarfs would be higher than those of the field brown dwarfs,
while both populations display variability down to the lower detection limit.
Measuring the viewing angles of these low-gravity brown dwarfs with maximum
amplitudes would allow us to test this possibility, and will be the focus of a
future paper. Additionally, using three-dimensional atmospheric circulation
models, Tan2021 find that the rotation rate likely affects the size of
atmospheric features in the atmosphere, and thus alters the variability
amplitude driven by these features.

##  8 Rotation Rates of an Age-Calibrated Sample

Identifying variable members of young moving groups is extremely valuable, as
it provides an age-calibrated sample of young brown dwarfs with measured
rotation rates, allowing us to trace angular momentum evolution (e.g.
Scholz2018; Zhou2019; Vos2020). Of particular interest is whether there are
rotational evolution trends related to mass, as is seen in our own solar
system (e.g. Allers2016). We estimate masses for objects whose variability
properties are presented in this paper (as described in the following section)
and gather estimated masses from the literature for the sample presented in
Vos2020.

###  8.1 Estimating Masses with Spectral Energy Distribution Analysis

We analyze the spectral-energy distributions (SEDs) of our young sample in
order to investigate their fundamental parameters. We calculate the
$L_{\mathrm{bol}}$, $T_{\mathrm{eff}}$, radii, masses and log(g) following the
method presented by Filippazzo2015. This method performs a numerical
integration of the combined empirical visible and near-infrared photometry and
spectra as well as WISE and/or Spitzer photometry, accessed through the BDNYC
database (bdnyc_db). We use parallaxes when available, but supplement with
kinematic distances when a target is regarded as a high-likelihood member of a
moving group. We present the fundamental parameters in Table LABEL:tab:SEDs.
Our sample spans temperatures of $\sim 600-1750$ K and model-dependent masses
of $2.5-33~{}M_{\mathrm{Jup}}$, overlapping with the masses and temperatures
of the directly-imaged exoplanets that have been discovered to date (e.g.
Marois2008; Marois2010; Lagrange2010; Bohn2020; Bohn2021).

###  8.2 Rotational Evolution as a Function of Time and Mass

In Figure 12, we show the current census of brown dwarf rotation rates as a
function of age and color-coded by mass. The majority of objects were
originally compiled in Vos2020. The $1-10$ Myr sample was vetted and compiled
by Moore2019 using $J$-band cutoffs to identify brown dwarfs. Scholz2018
calculate masses for the $\sim 2$ Myr Taurus sample, and Scholz2004 provide
masses for the $\sim 4$ Myr $\sigma$ Ori cluster sample. Although Rodriguez-
Ledesma2009 estimate masses for the 1 Myr sample, they are not provided in
their paper. We estimate their masses by comparing their distance and
reddening-corrected $I-$band magnitudes with the Baraffe2015 evolutionary
models, as described by Rodriguez-Ledesma2009. Similarly, Moore2019 do not
provide their estimated masses for the Upper Scorpius region. We estimate
masses for this sample by comparing the targets’ $J$-band magnitudes with the
10 Myr Baraffe2015 models, as described by Moore2019. For the Orion Nebula
Cluster, the Orion Belt region, Upper Sco and Taurus we have added a random
spread in the ages of each cluster member for clarity. The rotation periods of
brown dwarfs at age $>10$ Myr are compiled in Vos2020. We have added the
measured rotation rates of young objects measured in this paper and those that
have since been published by Tannock2021. Masses were either gathered from the
literature (Faherty2016; Filippazzo2015) or estimated in this paper (Table
LABEL:tab:SEDs). As described in Vos2020, we also plot angular momentum
conservation tracks for objects with masses of $10~{}M_{\mathrm{Jup}}$ and
$84~{}M_{\mathrm{Jup}}$ using the evolutionary models of Baraffe2015 as well
as the breakup period for objects with these masses.

Figure 12: Brown dwarf rotation period as a function of age and color-coded by
mass. Stars show bona-fide and high-likelihood members of young moving groups,
circles show candidate members of young moving groups, squares show field
brown dwarfs, which are plotted at 1 Gyr. Upside-down triangles are brown
dwarf cluster members. Triangles show substellar companions with measured
rotation rates. The 20 hr observation duration typically used by brown dwarf
Spitzer variability studies is shown by the dashed line. While the higher mass
brown dwarfs ($>30~{}M_{\mathrm{Jup}}$) can be seen to spin up with age, this
behavior is not observed in the low mass sample ($<30~{}M_{\mathrm{Jup}}$),
most likely because most low-mass objects lie in the intermediate age range.
Companion brown dwarfs and planetary-mass objects at $10-500$ Myr have
rotation periods consistent with those of isolated objects with similar
masses, indicating that they likely have similar formation pathways and
subsequent angular momentum evolution histories.

The majority of brown dwarf rotation rates $>10$ Myr have been measured with
Spitzer, and many surveys including Metchev2015 and this work observe for
$\sim$20 hrs. We show this $\sim$20 hr cutoff as the gray dashed line in
Figure 12. It is notable that we have detected a large number of variables in
the $10-500$ Myr age range with periods close to 20 hr, most notably a cluster
of AB Doradus members with periods of 10–20 hrs. This suggests that by
limiting our observation length, we are only sensitive to the fastest rotators
and may be missing some variable brown dwarfs with rotation rates $>20$ hr.
However, the high variability occurrence rates achieved by large surveys
suggest that the majority of rotation rates have been measured. Until we can
probe longer rotation rates by observing on timescales of many days, it is
important to keep in mind that we may not be measuring the full rotation rate
distribution of brown dwarfs and/or isolated and companion planetary-mass
objects.

We can also split the sample into two mass bins: the low-mass sample
($<30~{}M_{\mathrm{Jup}}$) and the high-mass sample
($>30~{}M_{\mathrm{Jup}}$). It is evident that the higher mass brown dwarf
sample spins faster with age, as expected from angular momentum conservation.
However, the spin up of the low-mass sample is not so evident. Based on the
angular momentum conservation tracks shown in Figure 12, we expect that the
high-gravity brown dwarfs should spin up more than their low-gravity
counterparts. However, we do still expect to see some spin-up for the low-mass
sample if we can probe a large enough age range. In particular, the long
rotation periods of the two older planetary-mass Y dwarfs WISE
J085510.83-071442.5 and WISE J140518.8+553421.3 (Cushing2016; Esplin2016) are
somewhat unexpected, as they are longer than the rotation periods of some
planetary-mass objects at younger ages. However, the $\sim 10$ hr rotation
period of Jupiter suggests that the lowest-mass objects may maintain longer
periods at older ages. A larger survey for variability in Y dwarfs (Cushing et
al. in prep) may help to provide more Y dwarf rotation periods, and reveal the
rotational behavior of low-mass objects at old ages.

Finally, we have included companion objects with measured rotation rates
(Bowler2020; Manjavacas2018; Manjavacas2019; Miles-Paez2019; Zhou2016;
Zhou2018; Zhou2020) in Figure 12 to investigate whether their angular momentum
evolution may be different from the isolated objects focused on in this paper.
Companion brown dwarfs and planetary-mass objects at ages 10–500 Myr have
rotation periods consistent with those of isolated objects with similar mass,
indicating that they likely have similar formation pathways and subsequent
angular momentum evolution histories. A similar conclusion was reached by
Bryan2018 based on rotational velocities of 11 low-mass objects. Based on our
current data, it seems reasonable to use the observed distribution of
planetary-mass brown dwarfs to estimate the expected rotation rates of
directly-imaged exoplanets when planning variability monitoring observations.

##  9 Variability Properties of AB Doradus Members

Within our sample, AB Doradus offers a unique opportunity to study the
variability properties of a sample of coeval objects. Four objects in our
survey are either bona fide or high-likelihood members of the AB Doradus
moving group: 2m0001, 2m0355, 2m1741 and 2m2206-42. While the high-likelihood
member 2m0326 was observed in this survey, its variability results are
somewhat inconclusive so we leave it out of the analysis. Four additional AB
Doradus BF/HLM members have variability monitoring observations in the
literature: *w0047, *2m2244, *2m1425 and *s1110 (Lew2016; Vos2018; Vos2020).

###  9.1 Variability Occurrence Rate of AB Doradus

First we use the methods outlined in Section 6 to determine the rate of
variability for the four AB Doradus members in our sample. We did not include
the other AB Doradus members presented in the literature to avoid biasing our
sample – some targets may have been observed with Spitzer after having been
previously found to exhibit variability while some non-detections may not have
been published. We detect significant variability in 3/4 objects in our
sample, with measured amplitudes from $0.3-0.7\%$. We show the probability
distribution of the variability occurrence rate for AB Doradus member in
Figure 13. We determine that the variability occurrence rate for this sample
is $>70\%$ at $68\%$ confidence and $>40\%$ at $95\%$ confidence. As we
uncover more L and T type members of young moving groups, we can begin to
compare occurrence rates across moving groups, and hence across ages.

Figure 13: Variability occurrence rate of AB Doradus members in our sample. We
determine that the variability occurrence rate for this sample is $>76\%$ at
$68\%$ confidence and $>50\%$ at $95\%$ confidence.

###  9.2 Variability Amplitudes and Rotation in AB Doradus

In Figure 14 we show all ten HLM/BF substellar AB Doradus members that have
been monitored in the Spitzer $[3.6~{}\mu$m] band, including this work and
literature values from Vos2018; Vos2020, plotted on a spectral type color
diagram. We see a large range of [$3.6~{}\mu$m] amplitudes for this sample
from $0.4-3.1\%$. It seems that the observed variability amplitude increases
as the temperature cools, at least from early-L to early-T spectral types.
This behavior is also seen in Figure 11 for both the low-gravity and field-
gravity populations, but it is notable that the amplitudes increase with
decreasing temperature for a coeval population of young objects.

There is a very small range in the measured periods within AB Doradus; all of
the variable objects have periods from $9-17$ hr. This clustering of rotation
periods is also evident in Figure 12 where there is a pile-up of objects in
this region. Since these objects are coeval, with similar masses, the
similarity in their periods suggest that they share common angular momentum
histories. Comparing these variability and rotation trends across other moving
groups will help to reveal the diversity in the atmospheres of coeval objects.
Moreover, combining Spitzer variability information for brown dwarfs with
observations of high- to low-mass stars with the Transiting Exoplanet Survey
Satellite (TESS) will allow us to trace variability properties and rotation
from the highest to lowest mass products of star formation in a coeval sample.

Figure 14: AB Doradus members monitored for variability in the Spitzer
$[3.6~{}\mu$m] band. AB Doradus with detections of variability are shown in
green, where the symbol size is proportional to the variability amplitude and
the rotation rate is shown by color. Non-detections are shown in red. The
field brown sample from Best2020 is shown in gray for context.

##  10 Lessons for variability monitoring of directly-imaged exoplanets

One of the goals of this survey is to establish expectations for detecting
variability in directly-imaged exoplanets. There has been considerable success
with variability searches in a handful of young wide-orbit companions with the
WFC3 instrument on the Hubble Space Telescope (Bowler2020; Lew2019; Zhou2016;
Zhou2019; Zhou2020) as well as some ground-based studies (e.g. Naud2017). Due
to their close proximity to their host star and relative faintness,
variability searches on directly-imaged exoplanets such as HR8799 bcde and
$\beta$ Pic b present considerable challenges. To date, there have been two
ground-based attempts to detect variability in HR8799 bc using the Very Large
Telescope (VLT) (Apai2016; Biller2021), however neither study reached the
sensitivity to detect variability, placing near-infrared upper limits of
$5-10\%$ on the variability amplitudes of the planets. Searching for
variability in these exoplanets will likely require next-generation telecopes
such as the James Webb Space Telescope (JWST) and 30-m class telescopes such
as the European-Extremely Large Telescope (E-ELT), the Thirty Metre Telescope
(TMT) and the Giant Magellan Telescope (GMT), which will provide unprecedented
sensitivity for variability searches in directly-imaged exoplanets. The survey
presented here outlines baseline expectations for the behavior of exoplanets
with similar temperatures, ages and masses.

Our survey suggests that the directly-imaged exoplanets are likely to have a
high variability occurrence rate across the full L and T spectral type regime.
Our results also show that we can expect a large range of variability
amplitudes for the directly-imaged exoplanets with very large amplitudes for a
subset of objects, potentially those with favorable viewing inclinations (e.g.
Vos2017). Comparing the variability behavior of the directly-imaged exoplanets
with the low-mass sample presented in this paper will help to put each
exoplanet in context with a larger sample of variables with similar
temperatures, masses and radii. Our survey also reveals the expected rotation
periods for young, directly-imaged exoplanets. With a median rotation period
of $\sim 12$ hr for young free-floating and companion objects, this work
highlights the importance of long duration observations in order to measure a
full period. For ground-based observations with the E-ELT, TMT and/or GMT
multiple nights of observations may be necessary to capture the periodicity,
however the fastest rotators could be prioritized by measuring $v\sin(i)$
values prior to variability searches (Wang2021). Finally, ongoing variability
surveys of Y dwarfs will be critical in revealing the potential for detecting
variability in the coldest exoplanets (Cushing et al. in prep).

##  11 Discussion and Conclusions

We have presented the most sensitive survey to date for variability in young,
low-gravity, giant planet analogs. We surveyed 26 brown dwarfs with evidence
of youth using the Spitzer Space Telescope at $3.6~{}\mu$m, 23 of which we
classify as young after a detailed study of their kinematic and spectroscopic
evidence of youth. We report the discovery of three new brown dwarfs in this
paper: 2m0349, 2m0951 and 2m0718, all of which show signatures of youth. Our
variability analysis, which makes use of periodogram analysis and the Bayesian
Information Criterion (BIC), identifies 15 variable and 11 non-variable
objects in the full sample. We compare the variability properties of the low-
gravity sample with those of their higher-mass field brown dwarf counterparts
reported in the literature, specifically focusing on the variability
occurrence rates, the variability amplitudes, and the rotation rates.

We determine the variability occurrence rates of young and field brown dwarfs
in three spectral bins: L2$-$L8, L9$-$T3.5 and $\geq$T4. For L dwarfs, we find
that both the young and field objects have variability occurrence rates of
$80-100\%$, i.e. virtually all L dwarfs are likely to be variable at the
$0.05-3\%$ level. Our variability occurrence rate calculations find a
tentative enhancement in the rate of variability in low-gravity objects with
spectral types L9-T3.5 compared to the field dwarf populations with similar
spectral types. The variability occurrence rate of low-gravity objects remains
high from L to T dwarfs, while the occurrence rate appears to drop for field
L/T transition and T dwarfs. The small sample size for young objects with
spectral types $\geq$T4 prevents us from placing meaningful constraints on
their variability occurrence rates. These results are encouraging for future
high-contrast searches for variability in directly-imaged exoplanets with
JWST. While a number of ground-based searches have been carried out for the
HR8799 planets (Apai2016; Biller2021) without detecting variability, JWST will
provide unprecedented sensitivity and precision. The results of these searches
will be directly informed by the variability studies of isolated low-gravity
brown dwarfs such as this one.

We find that the [$3.6~{}\mu$m] variability amplitudes increase steadily with
spectral type from L to early T, as is observed for the field brown dwarfs
(Metchev2015). This similar observed behavior supports our assumptions that
the same variability mechanisms are at play in both samples. We also note that
the maximum variability amplitudes observed for L dwarfs are higher for the
low-gravity objects compared to the field objects. Atmospheric models predict
that clouds are likely to form at higher altitudes in low-gravity atmospheres
(Marley2012), and these higher clouds might result in a larger observed
contrast between more cloudy and less cloudy regions. If this is the case, one
might expect higher variability amplitudes in the low-gravity sample. Indeed,
there are a number of young brown dwarfs displaying extremely large amplitudes
compared to their field brown dwarf counterparts such as 2m2139 ($\sim 26\%$
Radigan2012), VHS 1256 B ($\sim 25\%$ Bowler2020; Zhou2020) and PSO
J318.5$-$22 ($\sim 11\%$ Biller2015). Our work shows that not every young
brown dwarf exhibits such high amplitudes. We speculate that while the higher
contrast from higher altitude clouds may result in higher intrinsic
amplitudes, secondary effects such as viewing inclination (Vos2017) may be
responsible for the range in amplitudes that we observe.

We combine our new rotation rates with those compiled by Vos2020 to
investigate the angular momentum evolution of substellar objects from 1 Myr to
$>1$ Gyr. Firstly, we find that a large number of objects with ages 10–300 Myr
have rotation periods close to 20 hr, which is the observation duration for
many variability searches with Spitzer, including this work. This suggests
that we may be missing the slowest rotators within this sample. Secondly,
while the evolutionary spin-up of brown dwarfs with masses
$>30~{}M_{\mathrm{Jup}}$ is evident, this process is not so clear for the
lower mass brown dwarfs ($<30~{}M_{\mathrm{Jup}}$). Measuring rotation rates
of a larger number of low-mass brown dwarfs at ages $<10$ Myr and $>1$ Gyr
(e.g. Y dwarfs) will help to reveal their spin-up properties over time and as
a function of mass. We also compare the rotation rates of the small number of
low-mass companion objects that have been monitored to date, finding that
their rotation rates are consistent with those of the isolated low-mass brown
dwarfs. These results suggest that we can base future variability monitoring
observations of directly-imaged exoplanets on the rotation rates measured to
date for isolated low-mass brown dwarfs.

AB Doradus offers a unique opportunity to investigate the variability
properties of a sample of coeval objects. Within AB Doradus, we find a
variability occurrence rate of $>70\%$ and $68\%$ confidence. We note that the
amplitudes of AB Doradus members increase with decreasing temperature, as has
been noted for field brown dwarfs (Metchev2015) and low-gravity sources
(Vos2020, this work), but shown here for the first time in a coeval sample. We
also note an apparent clustering of rotation periods for AB Doradus objects,
with all of the measured periods in the range $9-17$ hr. The similarity in
their periods may suggest that these objects share common angular momentum
histories. Comparing these variability trends across other moving groups will
help to reveal how they change with age and mass.

Overall, the results from this survey are encouraging for variability searches
in directly-imaged exoplanets with future facilities such as JWST and 30-meter
telescopes, and the young brown dwarf sample presented here will serve as an
excellent baseline for comparison. Variability monitoring observations with
these facilities will enable unprecedented characterization of the weather
phenomena and angular momentum properties and directly-imaged exoplanets.

##  12 Acknowledgements

The authors would like to thank Azul Ruiz Diaz, Jai Glazer and Sophia Ameneyro
for their help and enthusiasm on this project through the NYC Science Research
Mentoring Program at the American Museum of Natural History. We would like to
thank the anonymous referee for thoughtful suggestions that significantly
improved the paper. We would like to thank Dr Seppo Laine and the Spitzer
Science Center for their expertise and advice on Spitzer data products and
analysis. This work is based on observations made with the Spitzer Space
Telescope, which is operated by the Jet Propulsion Laboratory, California
Institute of Technology under a contract with NASA. J. M. V. acknowledges
support by NSF Award Number 1614527, Spitzer Cycle 14 JPL Research Support
Agreement 1627378 and HST-GO-15924.001-A. J. V. and J.F. acknowledge support
from the Heising-Simons Foundation. This work has made use of data from the
European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia),
processed by the Gaia Data Processing and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. This work has benefited from
The UltracoolSheet, maintained by Will Best, Trent Dupuy, Michael Liu, Rob
Siverd, and Zhoujian Zhang, and developed from compilations by Dupuy2012,
Dupuy2013a, Liu2016, Best2018, and Best2020.

##  Appendix A Jeffrey’s Prior Derivation

The Jeffrey’s prior is a non-informative prior which is invariant under
parameter transformation. Jeffrey’s prior is defined in terms of the Fisher
information:

$J(\theta)\propto\sqrt{I(\theta)}$ (A1)

where the Fisher information $I(\theta)$ is given by:

$I(\theta)=-\mathbb{E}_{\theta}\left[\frac{d^{2}\ell}{d\theta^{2}}\right]$
(A2)

where $\mathbb{E}$ represents the expectation value and $\ell$ is the log
likelihood. For the likelihood presented in Section 6, the Fisher information
is given as:

$I(f)=-\mathbb{E}_{f}\left[\sum_{j}\frac{d_{j}}{f^{2}p_{j}}+\frac{p_{j}(1-d_{j})}{(1-fp_{j})^{2}}\right]$
(A3)

The expectation value of $d_{j}$ is simply the $fp_{j}$, so the Fisher
information can be expressed as:

$I(f)=\sum_{j}\frac{p_{j}}{f(1-fp_{j})}$ (A4)

and Jeffrey’s prior is thus given as

$J(f)=\sqrt{\sum_{j}\frac{p_{j}}{f(1-fp_{j})}}.$ (A5)

This represents a non-informative prior for estimation of the posterior
distribution of the variability occurrence rate, $f$. We plot the prior in
Figure 15.

Figure 15: Jeffrey’s Prior used for estimation of variability occurrence rate
in Section 6.