\section{Introduction} \label{sec:intro}
The \emph{Kepler}\xspace mission was instrumental in the discovery of thousands of exoplanets \citep{Borucki2010, Borucki2017}, most of which were found to be closer to their star than Mercury is to our Sun. A closer look at this population of short-period exoplanets revealed two prominent features: the Hot Neptune Desert and the Radius Valley. The Hot Neptune Desert (3$-$10\,$R_\oplus${}; $<$4\,days) is a region in the Gyr-old, period-radius plane with a significant dearth of Neptune-sized planets \citep{Beauge2013}. As Neptune-sized planets should be easier to find in short-period orbits, and many sub-Neptunes have been discovered with longer orbital periods from surveys such as \textit{CoRoT} \citep{baglin2006corot} and \textit{Kepler}, the Hot Neptune Desert appears to not be an observational bias. The Radius Valley, a region with a much lower frequency of planets with radii $\sim$1.8\,$R_\oplus${} than at $\sim$1.3\,$R_\oplus${} (Super-Earths) or at $\sim$2.4\,$R_\oplus${} (sub-Neptunes), was discovered by \citet{fulton2017california} using accurate stellar, and hence planetary, radii from the California-\emph{Kepler}\xspace survey \citep{petigura2017california, johnson2017california}. More recently, the Radius Valley was also detected in \emph{K2}\xspace data \citep{hardegree2020scaling}. 

Both the Hot Neptune Desert and the Radius valley are most likely caused by atmospheric mass loss due to XUV photoevaporation (see, e.g. \citealt{owen2013kepler, owen2017evaporation}) and/or core-powered mass loss (see e.g., \citealt{ginzburg2016super, ginzburg2018core, gupta2019sculpting, gupta2020corecool}). Given that these features are most likely to be evolutionary, it is likely that the primordial population of short-period planets was very different. Uncovering that population requires discovering exoplanets around young stars which can be further used to test the theories of atmopheric mass loss.

The \emph{Kepler}\xspace mission, however, was not able to detect any planets younger than 1\,Gyr since there were only four open clusters (0.7-9\,Gyr) in the \emph{Kepler}\xspace field \citep{meibom2011kepler, meibom2013same}. On the other hand, the \emph{K2}\xspace mission \citep{howell2014k2} discovered several young transiting planets while surveying the Taurus, Hyades, Praesepe and Upper Scorpius (e.g., \citealt{mann2016zodiacal_a, mann2016zodiacal_b, mann2017zodiacal, ciardi2018, rizzuto2017zodiacal, rizzuto2018zodiacal, gaidos2017zodiacal, vanderburg2018zodiacal}). Detailed studies of these planets have already begun to shed light on how these planets differ from their Gyr-old counterparts. For example, the planets in the 23\,Myr-old V1298 Tau system, that currently occupy gaps in the \emph{Kepler}\xspace distribution, is thought to be undergoing photoevaporation and is predicted to ``shrink" or lose their atmospheres with time \citep{david2019four}.





With the Transiting Exoplanet Sky Satellite (TESS) mission \citep{ricker2014transiting}, we now have the unique opportunity to detect planets around stars in more young clusters and associations than those observed with \emph{K2}\xspace. During its primary mission, TESS monitored each strip of the sky (also known as sectors) for approximately 27\,days, which means that it can detect planets with an orbital period of $\sim$14 days (assuming two transits per candidate). Per sector, there are about 200,000 bright nearby stars that were pre-selected for the  two year primary mission to be observed in a short 2-minute cadence mode known as the Candidate Target List (CTL; \citealt{stassun2018tess}). Moreover, TESS obtained images of each sector at a 30-minute cadence which are also known as Full-Frame Images (FFIs). Over a billion stars that were observed in this mode are listed as part of the TESS Input Catalog (TIC), which contains a relatively untapped reservoir of young planetary systems. In fact, while a handful of young planets have already been detected using 2-min cadence data within 200\,pc (e.g.,\citealt{newton2019tess,newton2021tess, rizzuto2020tess,mann2020tess}), only two have been detected thus far using the FFIs \citep{nardiello2020psf,bouma2020cluster}.

Automating the detection of young transiting planets poses several challenges, the biggest one primarily being the removal or modeling of instrumental systematics and stellar variability $-$ the latter of which is more important for young stars due to their high variability. Young stars also tend to have stellar rotation rates that are comparable to the orbital periods of short-period planets that can be detected by TESS \citep[e.g.,][]{Rebull2016}, thereby making it hard to find transiting planets in their light curves. A study led by \citet{nardiello2020psf}, \texttt{PATHOS}, used a PSF approach to extract light curves from TESS Primary Mission FFIs and searched for planets in them. However, most of their clusters were far away ($>$500\,pc) and while they were able to find 35 planets, only two of these were within 200\,pc thereby making precise radial velocity of these targets challenging. Additionally, TESS QuickLook Pipeline\footnote{https://archive.stsci.edu/hlsp/qlp} \citep{huang2020photometry_a, huang2020photometry_b, kunimoto2021quick} only considers stars brighter than T mag = 13.5 and is not optimized for highly variable stars such as the ones found in these clusters. To this effect, we have developed a dedicated pipeline to understand the survey completeness and thus properly calculate the intrinsic planet occurrence rates. 

Here, we present \texttt{pterodactyls}\xspace{}\footnote{https://github.com/rachelbf/pterodactyls} (\textbf{P}ython \textbf{T}ool for \textbf{E}xoplanets: \textbf{R}eally \textbf{O}utstanding \textbf{D}etection and \textbf{A}ssessment of \textbf{C}lose-in \textbf{T}ransits around \textbf{Y}oung \textbf{L}ocal \textbf{S}tars), a pipeline that builds on publicly available and tested tools in order to extract, detrend, search, and vet young transiting planet candidates detected using TESS Primary Mission 30-minute cadence photometry. In Section~\ref{sec:sample}, we discuss how we select the clusters that have been included in this analysis and, Section~\ref{sect:PTERODACTYLS} describes the inner workings of \texttt{pterodactyls}\xspace{}. Finally, in Sections~\ref{sec:results} \& \ref{sec:occ_rates}, we summarize our main results and some preliminary occurrence rates calculations followed by a discussion on our next steps (Section~\ref{sec:summary}).


\section{Sample of Young Stellar Clusters} \label{sec:sample}

\begin{table}[!htb]
    \centering
    \begin{tabular}{|c|c|c|c|c|}
    \hline\hline
    Cluster & Distance (pc) & Age (Myr) & Obs/Total\\
    \hline\hline
    THA & $46^{+8}_{-6}$ & 45$\pm$4 & 201/214\\
    IC~2602 & 146$\pm$5 & $46^{+6}_{-5}$ & 502/504\\
    UCL & 130$\pm$20 & 16$\pm$2 & 719/937\\    
    UMa & $\sim$25 & 414$\pm$23 & 16/17\\
    PiEri & 80$-$226 & 120 & 153/254\\    
    \hline\hline
    \end{tabular}
\caption{{Moving groups and clusters that have been homogeneously searched for transiting planets in this paper. Distance, age, and membership are from
\citet{gagne2018banyan} and \citet{babusiaux2018gaia}.} The last column specifies how many stars were observed during the TESS Primary mission vs. the total number of stars in the cluster.}
\label{table:data}
\end{table}


\begin{figure*}[!htb]
   \includegraphics[width=\textwidth]{pap1_cluster.pdf}
   \caption{Equatorial locations of stars in selected clusters used in this study} 
   \label{fig:cluster}
\end{figure*}

In order to develop and test \texttt{pterodactyls}\xspace, we limited our search to only those young clusters in which TESS, along with followup observations with the \emph{Spitzer}\xspace space telescope \citep{werner2004spitzer}, had already found transiting planets during its primary mission: Tucana-Horologium Association (THA), IC~2602, Upper Centaurus Lupus (UCL), Ursa Major (UMa), and Pisces-Eridani (PiEri). DS TUC~A\,b \citep{newton2019tess}, HIP~67522\,b and c \citep{rizzuto2020tess}, TOI~1726\,b and c \citep{mann2020tess}, and TOI~451\,b, c and d \citep{newton2021tess} were discovered using TESS 2-min cadence data and also had \emph{Spitzer}\xspace data. All of these are confirmed planets except for HIP~67522\,c. In contrast, the planets in IC~2602 were discovered using 30-min cadence data: TOI~837\,b which is a confirmed planet with mass measurements from radial velocity \citep{bouma2020cluster}, and TIC~460950389\,b \citep{nardiello2020psf} which is a Community TESS Object of Interest (CTOI)\footnote{\url{https://exofop.ipac.caltech.edu/tess/view\_ctoi.php}}. 

In summary, a total of ten planets have been discovered so far in the selected clusters: of them eight were discovered in TESS 2-min cadence data and two in 30-min cadence data. Followup observations have confirmed eight of these while HIP~67522\,c and TIC~460950389\,b are only candidate exoplanets.

For the selected clusters and associations, we relied on the memberships provided in the BANYAN $\Sigma$ \citep{gagne2018banyan} and the \textit{Gaia} DR2 open cluster member lists (Gaia Collaboration, \citealt{babusiaux2018gaia}) which gives us a total of 1940 stars out of which 1591 were observed during the TESS Primary Mission. Based on parameters available in the TIC, these stars span a large range in \texttt{Gaia}\xspace G magnitude (1.73-20.45)  and stellar mass (0.11-3.2\,$M_\odot$). Table \ref{table:data} provides a breakdown of the main properties of the moving groups and clusters that were included in this analysis while Figure \ref{fig:cluster} shows the relative positions of the clusters and their respective members that were analyzed in this work. 


\section{Transit Detection Pipeline} \label{sect:PTERODACTYLS}
\begin{figure*}[!htb]
   \includegraphics[width=\textwidth]{flowchart.pdf}
   \caption{Schematic of \texttt{pterodactyls}\xspace. $^\text{*}$Before passing to \texttt{triceratops}\xspace, we conducted a centroid test whenever the Target Pixel File (TPF) was available on MAST.}
   \label{fig:flowchart}
\end{figure*}


\texttt{pterodactyls}\xspace builds on publicly available and tested tools in order to extract the light curves (Section \ref{subsec:extract}), detrend them (Section \ref{subsec:detrend}), search for transit-like signals (Section \ref{subsec:search}), and vet planet candidates (Section \ref{subsec:vet}) in these clusters. 
We further optimized the efficiency of the detrending, search and vetting components of \texttt{pterodactyls}\xspace{} via injection-recovery tests. For the latter, we injected $\sim$1000 transit signals with planets ranging in size from 0.5-10\,$R_\oplus${} and with orbital periods between 0.5-10\,days into a sample of our light curves using the Python-based \texttt{batman}\footnote{\url{https://www.cfa.harvard.edu/\~lkreidberg/batman/}} package \citep{kreidberg2015batman}. We then processed the light curves with \texttt{pterodactyls}\xspace and examined the distribution of planets that were recovered. We used this planet distribution as feedback in order to improve \texttt{pterodactyls}\xspace{}. The specifics of the extensive testing using planet injection-recovery tests are summarized in the respective sub-sections while the main steps of the pipeline are summarized in Figure~\ref{fig:flowchart}.


\subsection{Light Curve Extraction}\label{subsec:extract}
Given that less than 10\% of our targets were part of the pre-selected TESS 2-min cadence targets, we use the primary mission's 30-min cadence photometry extracted from TESS FFIs. Here, we employed \texttt{eleanor}\xspace{}\footnote{http://adina.feinste.in/eleanor/}, an open-source tool developed by \citet{feinstein2019eleanor}. \texttt{eleanor}\xspace{} carries out background subtraction, and removal of spacecraft systematics such as jitter and pointing drift. \texttt{eleanor}\xspace provides four different types of light curves: raw light curves, corrected light curves, light curves with systematics removed using Principle Component Analysis (PCA), and Point Spread Function (PSF) modeled light curves.

In this work, we used PSF-modeled light curves since they better enable the recovery of small signals and minimize the effects of scattered light from the Earth and Moon \citep{feinstein2019eleanor}. \texttt{eleanor}\xspace creates these PSF light curves by modeling the shape of the PSF as a multi-variate Gaussian (major axis, minor axis, and a rotation) for each star at each cadence. Then, \texttt{eleanor}\xspace integrates that Gaussian function to calculate the flux of the target star at each cadence. The quality flags assigned to the FFI data by the TESS team are implemented during the extraction of the light curves by \texttt{eleanor}\xspace (see Table 32 in the TESS Science Data Products Description Document\footnote{\url{https://archive.stsci.edu/files/live/sites/mast/files/home/missions-and-data/active-missions/tess/_documents/EXP-TESS-ARC-ICD-TM-0014-Rev-F.pdf}}). By default, we only keep the data with the bit value 0 i.e., the highest quality data. Once we download the light curves for each sector for a particular target, we ‘stitch’ them together and run the search on all available Primary Mission data for the target.



\subsubsection{Accounting for Systematics}
Any instrumental systematics that affect the shape of the PSF, such as the temperature of the instrument, are well-captured by the PSF model and the flux is conserved across it. However, there were several instrumental effects in these light curves that were not captured during the PSF modeling and remain in the light curve such as mid-sector flux drops caused by data downlink, and offsets in the light curves caused by errors in the uploaded Guidestar tables (see TESS' Data Release Notes\footnote{\url{https://archive.stsci.edu/tess/tess_drn.html}} for a more detailed description of these issues). Since these instrumental effects would  generate spikes in the light curve that affect \texttt{pterodactyls}\xspace' ability to detrend and search for planets, we further optimized these light curves by masking the mid-sector flux drops and re-scaling any offsets. In some cases, \texttt{eleanor}\xspace could not properly fit the PSF flux of the target star possibly due to the target being faint, nearby bright stars, etc. We removed 148 such light curves from our sample.

We also masked a number of cadences in sectors 3, 10, 11, and 12 (see Figure~ \ref{fig:sem}) because we noticed that they triggered an anomalous number of transit-like detections in the search, indicating an issue with the light curves in those regions. We used a test similar to the ``Skye" metric used by \citet{thompson2018planetary} that was also adopted for \emph{K2}\xspace data by \citet{zink2020scaling}. The test measures the number of transit-like signals at each cadence and if the number of signals is $>$3$\sigma$, we masked those cadences before re-running the search.

\begin{figure*}[!htb]
  \includegraphics[width=\textwidth]{SEM_pap1.pdf}
  \caption{Skye Excess Metric for Sectors 3, 10, 11 and 12 respectively. Red dashed line depicts the 3$\sigma$ threshold.}
  \label{fig:sem}
\end{figure*}





\subsection{Light Curve Detrending}\label{subsec:detrend}
In order to detrend the light curves of these young stars, we explored several options presented in \texttt{W\={o}tan}\xspace, an open-source Python package on GitHub\footnote{\url{https://github.com/hippke/wotan/}} developed by \citet{hippke2019wotan}. In their study, \citet{hippke2019wotan} used injection recovery tests and found that, in the presence of significant stellar variation, spline-based methods (where adjacent functional pieces join each other at fixed knots) maximized the fraction of injected planets that were recovered. Following up on this result, we tested the following three different spline techniques (see Appendix~\ref{sec:detrend_test} for an example): i) a robust spline with iterative sigma clipping, ii) a Huber estimator spline \citep{huber1981robust} and, iii) a penalized spline with iterative sigma clipping \citep{eilers1996flexible}. We also tested the commonly used Savitzky-Golay filter \citep{savitzky1964smoothing} and found that it often over-fits for the transit thereby causing an inaccurate transit depth measurement, as pointed out in \citet{hippke2019wotan}. 

We evaluated the detrending based on the low noise or the root mean square (hereafter, RMS) of the light curve, and recovery fraction from the injection-recovery tests. We found that the penalized spline with iterative sigma clipping works best for our sample of young stars as it gives the lowest RMS measured on a 7-hour window and had the highest recovery fraction ($\sim$92\%). We found that, in agreement with \citet{nardiello2020psf}, the brighter stars contribute more flux and have low RMS. We also see a lot of scatter in the faint star region where there is more flux contamination. Furthermore, the ratio of the recovered planetary radii to the injected radii was close to unity, thereby suggesting that we did not over-detrend our light curves.

In carrying out these tests, we also noticed that the edges of light curves are often challenging to detrend due to the lack of data on one side and would often trigger a search. As such, we decided not to consider 0.5\,days worth of data on either edge. We also noticed that the penalized spline uses a cut of 2$\sigma$ i.e., it uses sigma clipping in estimating the trend and does not consider any data that is more than 2$\sigma$ away. While this works well for almost the entire sample, it does not do a good job for light curves with high amplitude ($>$2-3\%) rotational signals where it clips out peaks which leaves periodic residuals that are picked up by the search. We manually identified 20 such light curves and applied a 3$\sigma$ cut to them in order to achieve a better detrending.

\subsubsection{Optimization Based on Stellar Rotation Rates}
A penalized spline is essentially a piecewise polynomial fit in which the routine tries all numbers of splines up to a set maximum, and uses a penalty function to evaluate if a better fit outweighs the added degrees of freedom. The penalized spline provides additional flexibility: unlike other spline methods where the exact number of knots must be chosen manually, the penalized spline identifies the best number of knots by comparing it against the smaller residuals from an improved fit starting from a given maximum number. This feature turns out to be useful for our sample as some young stars are rapid rotators ($<$1\,day) with periodic signals on timescales similar to the orbital period. Using injection-recovery tests, we found that the highest planet detection efficiency could be achieved by limiting the maximum number of splines based on the stellar rotation rate. 

To this effect, we used a Lomb-Scargle periodogram (\citealt{lomb1976least, scargle1982studies}; implemented in \texttt{astropy}) to first calculate the stellar rotation rate and then implement the following step function to determine the maximum number of splines that should be used in order to best detrend any given light curve: 

\begin{itemize}
    \item if $R_{\text{rot}}$ $\leqslant$ 2 days: max splines = 200
    
    \item if $R_{\text{rot}}$ $>$ 2 days and $\leqslant$ 10 days: max splines = 200/$R_{\text{rot}}$
    
    \item if $R_{\text{rot}}$ $\geq$ 10 days: max splines = 25
    
\end{itemize}
If the periodogram was not able to find a period, it defaults to the length of the segment (half a sector since there is a break in the middle). Since that is $\sim$10 days, only a small number of splines (25) is used which is ideal for a light curve of low variability. However, this is a \textit{maximum} number of splines which means that the detrending can still pick the right number of splines even if the rotation period is estimated to be shorter than it really is.

At this stage, we visually inspected all of the light curves and found 80 with variablility trends that were difficult to detrend e.g., extremely fast rotators with a rotational rate $<$0.5\,days\footnote{Based on the stellar rotation rate, we ideally would use $\sim$500 splines to detrend the light curve but we realized that by doing so, we would likely be fitting noise.}. Since we were not able to properly detrend these light curves, we also would not have been able to find planets in them. As such, we removed these light curves from our sample and are left with $\sim$1359 light curves that we used to search for planets.



\subsection{Planet Candidate Search}\label{subsec:search}
We searched the detrended light curves for periodic, transit-like signals with \texttt{transitleastsquares} (hereafter, \texttt{TLS}\xspace; \citealt{hippke2019optimized}). \texttt{TLS}\xspace is optimized for signal detection of small planets by taking into account the effects of limb darkening to better model the ingress and egress shape \citep{manduca1977limb, mandel2002analytic}. \texttt{TLS}\xspace searches for these signals by phase-folding the light curve over a range of epochs, transit durations and orbital periods calculated using a function based on the mass and radius of the host star as detailed in \cite{ofir2014optimizing}. The period search grid is set by \texttt{TLS}\xspace by dividing the length of the light curve by the minimum number of required transits which is set to two.

For the purposes of our work, we defined a Threshhold Crossing Event (TCE) as having a signal-to-noise ratio (SNR) of 7 and a signal detection efficiency (SDE) of 7\footnote{SDE is essentially the SNR of the signal in Fourier space. For more details, see \citet{kovacs2002box}}. Most transit survey studies choose an SDE value between 6 and 10 (see e.g., \citealt{dressing2015occurrence, livingston2018sixty, siverd2012kelt, pope2016transiting, aigrain2016k2sc, feliz2021nemesis}). While a lower SDE would mean higher completeness, a higher SDE would reduce the number of false positives. Our chosen SNR and SDE cuts reduced the number of false positives, which arose from detrending artifacts caused by the significant ($\sim$20\%) number of rapid rotators with periodic signals comparable to stellar rotation rates. Additionally, an SNR and SDE cut of 7 recovers eight out of the ten the previously detected planets (see Section \ref{sub:known}).





\subsection{Candidate Vetting} \label{subsec:vet}
Transiting surveys (like \emph{Kepler}\xspace, \emph{K2}\xspace and TESS) that target a large number of stars face a major challenge while distinguishing false alarms (transit-like events caused by instrumental artifacts) and astrophysical false positives (transit-like events that are astrophysical in origin but are not caused by planets). The process of sorting out false alarms and false positives from real transit signals is known as vetting. In order to filter out any instrumental false alarms that could be due to detrending artifacts, we used criteria inspired by \texttt{RoboVetter} (\emph{Kepler}\xspace; \citealt{thompson2018planetary}) and \texttt{EDI-Vetter} (\emph{K2}\xspace; \citealt{zink2019edivetter, zink2020scaling}) to create four checks in order to address the issues unique to our sample. Only the TCEs that pass the following tests are considered objects of interest: 

\begin{itemize}[leftmargin=*]
\item \textbf{At least two transits with data:} 
We required two transit signals with in-transit data points in a given light curve which allowed us to find planets with orbital periods of up to 13.5\,days around most of the stars (assuming each star only had single sector data). If we required at least three unique transits like \emph{Kepler}\xspace did, we would be limited to transit signals with orbital periods less than 9\,days.

\item \textbf{Orbital period dissimilar to stellar rotation rate:} 
While analysing the results of the search, we often found residuals at orbital periods similar ($\sim$20\%) to the stellar rotation rate (and some of these at half the rotation rates), despite detrending. To combat this, we flag any such TCEs.

\item \textbf{Consistency in individual transit depths:} 
We found that, quite often, false alarms have transit depths with larger variance than real planets. Similar to \emph{Kepler}\xspace's Robovetter Single Event Statistic (SES), we flagged signals with large variations in transit depth as likely false positives. We weed out these signals based on the ratio of the standard deviations of the individual transit depths to the average transit depth, which we set to be greater than 1.5 for a false positive. This cutoff was determined by preserving known planets while keeping the number of false positives low.

\item \textbf{Observed transit duration should be consistent with expected transit duration:} 
On visual inspection, we found that some TCEs had transit durations that were inconsistent with those of planetary transits. For example, stellar variability can lead to dips with a duration longer than a transit, and single occurrence outliers can appear to have a shorter transit duration. In such cases, it is beneficial to compare the observed transit duration to the expected transit duration that is determined using the planetary parameters calculated by \texttt{TLS}\xspace. If the observed transit duration is within 75\% of the expected transit duration, the TCE passes this test. This cut was also set to preserve known planets and reduce the number of false positives.
\end{itemize}

In addition, we also used two checks from \texttt{EDI-Vetter Unplugged} \citep{zink2020scaling} that were initially designed to vet transits in \emph{K2}\xspace data but has now been adapted to vet TESS light curves (for e.g., see \citealt{feliz2021nemesis}).

\begin{itemize}[leftmargin=*]
\item \textbf{Individual Transit Test:} 
Derived from the Marshall test designed by \citet{mullally2016identifying} for \emph{Kepler}\xspace, this test used transit models and those of common non-transiting false positives and fit these to the individual transits. It then used the Bayesian Information Criterion (BIC) to gauge the goodness of these fits, thereby flagging the non-transiting signals that are often assumed to belong to transiting planets.


\item \textbf{Secondary Eclipse Test:} 
This test looked for transit events between the 40-60\% phase range of the folded light curve, allowing for some mild eccentricity variations.  If a meaningful secondary eclipse was identified i.e., the secondary eclipse has a transit depth $>$10\% of the primary transit, this flag was triggered and the transit was labeled a false positive. Secondary eclipses with less than 10\% depth can be associated with reflected light from hot Jupiters.
\end{itemize}


Once we removed the clear false alarms and only have transit-like signals, we were left with 21 \texttt{pterodactyls}\xspace Objects of Interest (hereafter POIs; see Tables \ref{table:known_planets} and \ref{table:binaries}) on which we conducted a centroid test \citep{hedges2021vetting}\footnote{https://pypi.org/project/vetting/} in order to determine whether the transiting signal was indeed coming from the target star (see Appendix~\ref{sec:centroid}). This test compares the distributions of centroids for the data inside and outside the transit. An offset would indicate that the transiting signal was coming from a contaminating source such as a nearby or background star. We chose a probability \textit{p} of 5\% for the Student t-test below which the centroids positions inside and outside transits are not considered to come from the same parent distribution.







All POIs were then passed on to \texttt{triceratops}\xspace\footnote{https://github.com/stevengiacalone/triceratops} \citep{giacalone2020vetting} $-$ a vetting tool designed to distinguish bona fide planets from astrophysical false positives in TESS data. \texttt{triceratops}\xspace calculated the marginal likelihood of each transit-producing scenario using observational constraints such as the transit depth and transit duration from \texttt{TLS}\xspace for each potential host, positions and magnitude of nearby stars from \texttt{Gaia}\xspace \texttt{dr2} \citep{gaiadr2}, priors on planet occurrence and stellar multiplicity rates, and models of planetary transits and eclipsing binaries. These scenarios include transiting planets and eclipsing binaries originating from the target star, known nearby stars, and unknown unresolved stars located near the target on the sky. Lastly, \texttt{triceratops}\xspace{} calculates the Bayes factor between all scenarios involving a transiting planet around the target star and all astrophysical false positive scenarios. Using this Bayes factor, \texttt{triceratops}\xspace can place a statistically robust constraint on the reliability of any given planet candidate detected by TESS. Here, we made a transiting planet probability cut of 50\% in order to preserve known planets.




\section{Fitting for Planet Properties}
Once the POIs have passed all of the vetting tests described in Section \ref{subsec:vet}, we fit the phase-folded light curve\footnote{Since the PSF light curves from \texttt{eleanor}\xspace do not currently have uncertainties, we use the RMS of the detrended light curve as errors while fitting for planet parameters} in order to derive a more accurate planet radius. Since \texttt{TLS} only provided an initial fit to the phase-folded light curve and did not provide uncertainties on the planet radius, we used \texttt{TLS}\xspace' outputs as priors for a more rigorous optimization in order to derive uncertainties using \texttt{EXOTIC} \citep{zellem2020utilizing}. Here, the transit parameters are optimized using the multi-modal nested sampling algorithm called \texttt{UltraNest} \citep{buchner2021ultranest} which is a Bayesian inference tool that uses the Monte Carlo strategy of nested sampling to calculate the Bayesian evidence allowing simultaneous parameter estimation and model selection. A nested sampling algorithm is efficient at probing parameter spaces which could potentially contain multiple modes and pronounced degeneracies in high dimensions; a regime in which the convergence of traditional Markov Chain Monte Carlo (MCMC) techniques becomes incredibly slow (\citealt{Skilling2004}; \citealt{Feroz2008}). This approach has been validated on \textit{Kepler} light curves, TESS data and ground-based measurements \citep{Pearson2018, Pearson2019a, Pearson2019b}.



\section{Results} \label{sec:results}
We processed $\sim$1359 young stellar light curves extracted from TESS FFIs i.e., 30-min cadence data through \texttt{pterodactyls}\xspace (see Appendix~\ref{sec:tics} for the TIC IDs processed in this work). As an illustration, Figure~\ref{fig:pieri_d} shows the different steps that were used in the recovery of the planet TOI~451\,d. We detected 21 POIs that were then passed to \texttt{triceratops}\xspace where their phase-folded light curves were vetted to separate transiting planet signals from astrophysical false positives.

\begin{figure*}[!htb]
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{257605131_d.pdf}
\endminipage\hfill
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{TIC257605131_d_fit.pdf}
\endminipage
\caption{\textbf{Left:} Recovery of TOI~451\,d. Upper panel: PSF flux (in dark gray) and penalized spline trend (in blue). Lower panel: Detrended light curve (PSF flux divided by the trend). The multi-transit signature (in red) can be clearly seen. \textbf{Right:} Upper panel shows Phase-folded light curve of TOI~451\,d (gray data points) along with the best fit model using \texttt{EXOTIC}\xspace (in red). Lower panel shows the residuals in parts per million.}
\label{fig:pieri_d}
\end{figure*}

\subsection{Recovery of Known Exoplanets} \label{sub:known}

\begin{table*}[!htb]
    \centering
    \begin{tabular}{|c|c|c|c|c|c|c|c|} 
    \hline\hline
Planet & Ref. & Cluster & Radius (pub; $R_\oplus$) & Period (pub; days) & Radius (rec; $R_\oplus$) & Period (rec; days) & TP prob \%\\ 
 \hline
DS Tuc A\,b$^{*}$ & a & THA & 5.70$\pm$0.17 & 8.138268$\pm$0.000011 & $5.349^{+0.464}_{-0.458}$ & 8.215$\pm$0.066 & 54.93\\

TIC~460950389\,b & b & IC~2602 & 3.8$\pm$0.2 & 2.8622$\pm$0.0004 & $3.554^{+0.320}_{-0.304}$ & 2.861$\pm$0.005 & 99.69\\

TOI~837\,b & b,c & IC~2602 & $6.9^{+0.6}_{-0.4}$; $8.90^{+0.74}_{-0.71}$ & 8.3252$\pm$0.0006 & $7.224^{+0.410}_{-0.387}$ & 8.325$\pm$0.016 & 97.45\\

HIP~67522\,b$^{*}$ & d & UCL & $9.72^{+0.48}_{-0.47}$ & 6.95993$\pm$0.00034 & $13.400^{+1.354}_{-1.328}$ & 6.893$\pm$0.019 & 86.06\\

HIP~67522\,c$^{*}$ & d & UCL & $8.01^{+0.75}_{-0.71}$ & $54^{+70}_{-24}$ & xx & xx & xx \\

TOI~1726\,b$^{*}$ & e & UMa & 2.15$\pm$0.10 & $7.10793^{+0.00046}_{-0.00034}$ & $2.134^{+0.214}_{-0.213}$ & 7.116$\pm$0.043 & 75.32\\

TOI~1726\,c$^{*}$ & e & UMa & 2.67$\pm$0.12 & 20.5455$\pm$0.0011 & $2.225^{+0.264}_{-0.267}$ & 20.5234$\pm$0.1022& 75.72\\

TOI~451\,b$^{*}$ & f & PiEri & $1.94^{+0.15}_{-0.34}$ & $1.858701^{+0.000027}_{-0.000033}$ & xx & xx & xx\\

TOI~451\,c$^{*}$ & f & PiEri & 3.07$\pm$0.14  & $9.192523^{+0.000064}_{-0.000084}$ & $2.604^{+0.308}_{-0.294}$ & 9.195$\pm$0.036 & 71.11\\

TOI~451\,d$^{*}$ & f & PiEri & 4.03$\pm$0.15 & $16.364981^{+0.000047}_{-0.000049}$ & $3.678^{+0.306}_{-0.294}$ & 16.368$\pm$0.070 & 70.63\\
\hline\hline
\end{tabular}
    \caption{Published and Recovered Properties of planets detected in young clusters. References: a) \citet{newton2019tess};
    b) \citet{nardiello2020psf}; c) \citet{bouma2020cluster};
    d) \citet{rizzuto2020tess}; e) \citet{mann2020tess}; f) \citet{newton2021tess}
    $^{*}$Planet originally detected in 2-min cadence data.}
    \label{table:known_planets}
\end{table*}

In the five clusters that were searched in this work, there are eight confirmed planets and two planet candidates. More specifically, three were single planet systems (DS Tuc~A\,b, TIC~460950389\,b, and TOI~837\,b; the latter two discovered in 30min-cadence data) and seven were part of three multi-planet systems (HIP~67522\,b and c, TOI~1726\,b and c, and TOI\,451\,b, c and d), all discovered in 2-min cadence data (see Section~\ref{sec:sample}). As \texttt{pterodactyls}\xspace, which is automated, is built to search for just one planet per star, the discovery of additional planets in a system is carried out after finding a first planet and masking it out.

Among the 21 POIs detected by \texttt{pterodactyls}\xspace in the first automated run, we find six of the confirmed exoplanets that were previously detected by various teams using TESS (summarized in Table \ref{table:known_planets}; for light curves, see Appendix~\ref{sec:recovery}). \texttt{pterodactyls}\xspace was able to recover the values for the orbital periods and radii for three of these planets that were consistent (within uncertainties) with the literature. This includes two single planets (TIC~460950389\,b, and TOI~837\,b) that were initially discovered in 30-min cadence data, and one multi-planet system (TOI~451\,b, c and d) that was initially discovered in 2-min cadence data. We first recovered TOI~451\,d (4.1\,$R_\oplus$) at an orbital period of 16.36\,days. After masking this transit, we were also able to recover TOI~451\,c (3.1\,$R_\oplus$) at 9.19\,days (see Figures~\ref{fig:pieri_d} and \ref{fig:pieri_c}). However, \texttt{pterodactyls}\xspace was not able to recover TOI~451\,b (1.9\,$R_\oplus$) since it has a low SNR ($\sim$5). We find that all of the confirmed planet and planet candidates that \texttt{pterodactyls}\xspace recovers pass the centroid test.

Two out of the other three recovered transiting systems (out of a total of six as stated above) suffered from sparse cadence (few in-transit data points): DS Tuc~A\,b and HIP~67522\,b. For these, \texttt{pterodactyls}\xspace was able to accurately recover the orbital periods but not the same radii as in the literature. In the case of DS Tuc~A\,b, at first we considerably underestimated the planet's radius value. On further inspection, we noticed that this was caused by very few in-transit data points that was due to the removal of ``bad" data that was flagged using the TESS team's quality flags. In this specific case, we found that inclusion of certain quality flags lead to a better estimation of the planetary radius i.e 5.35\,$R_\oplus$ (see Figure~\ref{fig:dstuc}) which is within errors of the value of 5.70\,$R_\oplus$ that was measured using TESS' 2-min cadence and \emph{Spitzer}\xspace data (see Table~\ref{table:known_planets}). In the case of HIP~67522\,b, the overestimation of the planet radius was severe: we recovered a radius of 13.4\,$R_\oplus$ while the planet's observed radius is 9.72\,$R_\oplus$. This could be due to several factors. First, as can be seen in Figure~\ref{fig:ucl_hip67522}, the middle transit was not recovered (but predicted by \texttt{TLS}\xspace) due to the default edge cutoff for all stars being 0.5\,days. However, the severe overestimation of the radius could also be due to the planet being discovered in the more abundant 2-min cadence data. We tested this by processing the publicly available 2-min cadence light curve through \texttt{pterodactyls}\xspace (see Figure~\ref{fig:hip67522_2mins}) and found that with the increased number of in-transit data points (along with higher RMS compared to 30-min cadence data), we were able to recover a planet radius of 9.8\,$R_\oplus$ which is close to the literature value of 9.72\,$R_\oplus$. We were not able to recover HIP~67522\,c due to it being a single transit event and hence, beyond the scope of \texttt{pterodactyls}\xspace since it requires a minimum of two transits.

Finally, for the TOI~1726 system in UMa, we recovered a transit signal from a 2.48$R_\oplus$\ planet with an orbital period of 11.32\,days. However, the TOI-1726 system has two small planets (b and c) of radii 2.15\,$R_\oplus$ and 2.67\,$R_\oplus$, with orbital periods of 7.11 and 20.55\,days respectively \citep{mann2020tess}. Given that the planet sizes (and hence, transit depths) are similar, the signals were initially detected by \texttt{pterodactyls}\xspace as coming from one planet. Upon further inspection, we realized that \texttt{TLS}\xspace set the maximum orbital period in the grid based on the minimum number of transits. So, for single sector data ($\sim$27\,days) with the minimum number of transits set to 2, the maximum orbital period would be $\sim$13.5\,days which is shorter than the orbital period of TOI~1726\,c. With this in mind, for this case alone, we set the minimum number of transits to one and were able to recover TOI~1726\,c (see Figure~\ref{fig:uma_c}). We then masked this transit and re-processed it through \texttt{pterodactyls}\xspace and recovered TOI~1726\,b as well (see Figure~\ref{fig:uma_b}).


After implementing the aforementioned fixes, we find that all of the recovered orbital periods are the same as those reported in the discovery papers while for the radii, four are the same within 1$\sigma$ and three within 2$\sigma$ (Figure~\ref{fig:stellar_comp}). The recovery of these planets demonstrates the usefulness of \texttt{pterodactyls}\xspace in discovering young planets in 30-min cadence data.



\begin{figure*}[!htb]
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{166527623.pdf}
\endminipage\hfill
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{TIC166527623_b_fit_lowrms.pdf}
\endminipage
\caption{\textbf{Left:} Recovery of HIP 67522\,b. Upper panel: PSF flux (in dark gray) and penalized spline trend (in blue). Lower panel: Detrended light curve (PSF flux divided by the trend). The multi-transit signature (in red) can be clearly seen. \textbf{Right:} Upper panel shows Phase-folded light curve of HIP 67522\,b (gray data points) along with the best fit model using \texttt{EXOTIC}\xspace (in red). Lower panel shows the residuals in parts per million.}
\label{fig:ucl_hip67522}
\end{figure*}

\begin{figure*}[!htb]
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{hip67522b_2mins.pdf}
\endminipage\hfill
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{TIC166527623_b_fit_2mins.pdf}
\endminipage
\caption{\textbf{Left:} Comparison of phase-folded light curve of HIP 67522\,b using 2-min cadence data (orange) and 30-min cadence data (black). \textbf{Right:} Upper panel shows Phase-folded light curve of HIP 67522\,b using 2-min cadence (dark gray points) along with the best fit model using \texttt{EXOTIC}\xspace (in red). Lower panel shows the residuals in parts per million.}
\label{fig:hip67522_2mins}
\end{figure*}


\begin{figure*}[!htb]
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{planet_rad.pdf}
\endminipage\hfill
\minipage{0.5\textwidth}
  \includegraphics[width=\linewidth]{planet_per.pdf}
\endminipage
\caption{\textbf{Left:} Comparison of planetary radii recovered using \texttt{pterodactyls}\xspace (fitted with \texttt{EXOTIC}\xspace) vs. known planet radii for the eight recovered known planets. As can be seen, despite being recovered from 30-min cadence data, the radii are often consistent with those published. In two cases, TOI~1726\,c and TOI~451\,c, we slightly underestimate the planet's radius. We also overestimate the radius of the planet for HIP~6752\,b. \textbf{Right:} Comparison of orbital periods recovered using \texttt{pterodactyls}\xspace vs. known orbital periods for the eight recovered planets to depict consistency between the results of our work and those from the discovery papers. Note: Orbital period errors are extremely small and hence cannot be seen in the plot.}
\label{fig:stellar_comp}
\end{figure*}


\subsection{Astrophysical False Positives Vetting} \label{sub:fps}
Apart from the six confirmed planets, \texttt{pterodactyls}\xspace also detected 15 other POIs (see Table \ref{table:binaries}), including two eclipsing binary signals around two stars in UCL i.e., TIC~129116176 and TIC~148158540. We weren't able to conduct centroid tests on all of these targets since the Target Pixel File on the MAST server (which is used for the centroid test) were only available for three of the candidates (all of the ones in UCL) in Table~\ref{table:binaries} but not for the 12 other candidates (IC2602). As such, we conducted the centroid test whenever the data was available. For the targets, where we could not conduct the centroid test, we visually inspected the eleanor pixel by pixel light curves. While determining the probability of a transiting signal being due an astrophysical false positive or a transiting planet, \texttt{triceratops}\xspace (by default) takes into account planet occurrence rates priors that are based on the Gyr-old exoplanet population (e.g., \citealt{howard2012planet, dong2013fast,petigura2013prevalence,dressing2015occurrence,mulders2015stellar,Mulders2018}). However, if the planet population around young stars is different from that around Gyr-old stars, for example because young planets are larger in size, the probability that \texttt{triceratops}\xspace would calculate is likely to be biased. In order to test the effect of priors on the probabilities of the different scenarios investigated by \texttt{triceratops}\xspace, we also calculated the astrophysical false positive probability using non-informative (flat) priors that were recently implemented in \texttt{triceratops}\xspace.

Out of our remaining 15 POIs, we found that 11 were not affected by the choice of priors and had the same probability of being an astrophysical false positive in both cases. 
In two out of the remaining four, we detected a transiting signal in each each of them (TIC~466579089 and TIC~466035035; see Table \ref{table:binaries}) that had passed all our vetting tests. Using Gyr-old priors, \texttt{triceratops}\xspace calculated the most likely scenario to be a nearby eclipsing binary (NEB). However, when using flat priors, the most likely probability was calculated to be a transiting planet. Given that IC~2602 is an extremely crowded cluster ($\sim$5 stars per TESS pixel), we used \texttt{eleanor}\xspace's pixel level light curves to investigate whether these transits were indeed around our target stars. In doing this, we found that the signal was coming from giant stars (TIC~466579123 and TIC~466035000, respectively) a couple of pixels away and actually contributed more flux to the pixel than our target star. With this information in hand, we ran \texttt{pterodactyls}\xspace on those stars and were able to detect the same signal in the light curve but, this time, \texttt{triceratops}\xspace found that the most likely scenario was that the signal was due to an eclipsing binary around TIC~466579123 and TIC~466035000.

When the false positive probability was calculated using flat priors, the remaining two POIs that were initially calculated to be nearby eclipsing binaries with Gyr-old priors, showed a $>$25\% chance that the signal was either a nearby eclipsing binary or a nearby transiting planet. These POIs have large transits depths that could either be attributed to small stars, brown dwarfs or giant planets since the two populations have comparable radius distributions \citep{chen2017prob}. This, along with the fact that identifying shallow secondary transits can be arduous, makes it harder for statistical validation tools such as \texttt{triceratops}\xspace to properly distinguish between these scenarios \citep{shporer2017three}. In cases such as these, masses obtained via radial velocity follow-up are necessary to differentiate between these scenarios. This illustrates the need to properly take into account flux contamination in young clusters while vetting planet candidates.

\begin{table*}[!htb]
    \centering
    \begin{tabular}{|c|c|c|c|c|c|} 
    \hline\hline
TIC~ID & Cluster & Period (d) & Most likely (Gyr-old priors) & Most likely (Flat priors) & Potential Host Star(s)\\ 
 \hline
390758843 & IC~2602 & 0.890 & 25\% NEB & 25\% NEB & TIC~390758947\\ 

376731139 & IC~2602 & 12.769 & 50\% NEBx2P & 49\% NEBx2P & TIC~376731097\\

911870867 & IC~2602 & 20.018 & 42\% NEB & 40\% NTP; 27\% NEB & TIC~390843430\\

911870847 & IC~2602 & 20.018 & 50\% NEB & 37\%NEB; 30\% NTP & TIC~390843430\\ 

315312718 & IC~2602 & 2.203 & 28\% NEB & 30\% NEB & TIC~315312721; TIC~315312651\\

466035035 & IC~2602 & 1.522 & 99\% EB & 92\% TP$^{x}$ & TIC~466035000\\

466579089 & IC~2602 & 15.910 & 91\% NEB & 44\% NEB; 44\% TP$^{x}$ & TIC~466579123\\

913701079 & IC~2602 & 2.900 & 88\% NEBx2P & 84\% NEBx2P & TIC~465306570\\ 

465306494 & IC~2602 & 2.899 & 86\% NEBx2P & 85\% NEBx2P & TIC~465306570\\

378611172 & IC~2602 & 7.542 & 89\% NEBx2P & 87\% NEBx2P &  TIC~378611218\\

466422383 & IC~2602 & 33.247 & 99\% NEBx2P & 99\% NEBx2P & TIC~466422466\\

390757974 & IC~2602 & 1.705 & 54\% NEB & 54\% NEB & TIC~390758017\\

177631209 & UCL & 2.144 & 32\% EB & 21\% PEB; 21\% NEB & TIC~177631209; TIC~1167076020 \\ 

129116176 & UCL & 4.032 & 74\% EBx2P  & 91\% EBx2P & - \\

148158540 & UCL & 3.700 & 65\% EBx2P & 64\% EBx2P & -\\


\hline\hline
\end{tabular}
    \caption{List of POIs/Potential Binaries detected by \texttt{pterodactyls}\xspace and vetted using \texttt{triceratops}\xspace in both Gyr-old and flat prior mode. \textbf{Key:} TP = Transiting Planet around Target Star; EB = Eclipsing Binary around Target Star; EBx2P = Eclipsing Binary around Target Star with twice the orbital period compared to what was detected by \texttt{TLS}\xspace; NTP = Nearby Transiting Planet; NEB = Nearby Eclipsing Binary; NEBx2P = Nearby Eclipsing Binary with twice the orbital period compared to what was detected by \texttt{TLS}\xspace; PEB: Eclipsing Binary around Primary Star with Unresolved Bound Companion.
    $^x$See discussion in Section \ref{sub:fps} Note: The signals found in TIC~911870867 \& TIC~911870847 as well as TIC~913701079 \& TIC~465306494 have similar orbital periods with the same potential host star suggesting that they might be blended or part of the same multiple system.}
    \label{table:binaries}
\end{table*}


\section{Planet Occurrence Rates in investigated clusters}\label{sec:occ_rates}
\subsection{Detection Efficiency}
An advantage of having an automated planet detection pipeline is that it enables the measurement of the survey detection efficiency which is necessary to compute intrinsic occurrence rates. The detection efficiency of our survey was calculated using injection-recovery tests on a 3$\times$4 grid in $\frac{R_\text{p}}{R_\star}$\  from 0.01 to 0.3  and in orbital period from 0.5 to 27\,days , equally spaced in both $\frac{R_\text{p}}{R_\star}$ and orbital period. The injections are carried out in the $\frac{R_\text{p}}{R_\star}$\ space because most stars in our clusters currently lack stellar properties. We injected a transit signal in each bin in each light curve in our sample while taking into account the flux contribution of each star as calculated by \texttt{triceratops}\xspace. We then passed these injections through \texttt{pterodactyls}\xspace{} and calculated the fraction that were recovered as planets to produce a detection efficiency map per cluster as well as an average map for the entire sample. The average detection efficiency map is presented in Figure~\ref{fig:all_det_eff} and shows that the overall detection efficiency does not exceed 50\% for any bin below 0.1 $\frac{R_\text{p}}{R_\star}$ (i.e., a Jupiter-sized planet assuming it orbits a Sun-like star). Interestingly, all the previously known recovered planets (denoted using blue stars in Figure~\ref{fig:all_det_eff}) are in bins with fairly low values ($<$20\%). This might be due to the intrinsic occurrence of smaller vs. larger planets being higher \citep{mulders2015planet,petigura2018california, mulders2018exoplanet}. 

We also investigated how the detection efficiencies of each cluster contributed to the overall detection efficiency. Figure~\ref{fig:det_eff_comp} illustrates that the detection efficiency of a cluster can be quite different from the average. For example, Pisces-Eridani (top panel) has a significantly larger detection efficiency. IC~2602 (middle panel), on the other hand, has an extremely low detection efficiency. These differences might be due to age (IC~2602 is younger, hence its stars might be more variable) and/or to flux contamination due to stellar crowding, which is significantly more pronounced in IC~2602.


We further investigated the effects of flux contamination on IC~2602's detection efficiency by re-running the injection recovery tests without accounting for flux contamination. In other words, signals are injected into the light curve assuming that all of the flux comes from the target star and hence, there is no transit depth dilution. The detection efficiency map is given in the bottom panel of Figure~\ref{fig:det_eff_comp} and it is similar to that of Pisces Eridani within 10\,days and above 0.06$\frac{R_\text{p}}{R_\star}$. This test demonstrates that flux contamination significantly reduces the detectability of all planets but stellar variability in the youngest associations still plays a role for the smaller planets.

\begin{figure}[!htb]
    \centering
    \includegraphics[width=1.0\linewidth]{All-small-highres-binned.png}
    \caption{Overall detection efficiency of \texttt{pterodactyls}\xspace while taking into account flux contamination. The confirmed planets are depicted using blue stars, except for TOI~1726\,b whose recovery required changing \texttt{pterodactyls}\xspace. Darker bins represent regions of higher detection efficiency. ${n_\star}$ is the total number of stars in our sample while \textit{n} is the total number of injections done. Black box denotes the bin over which the intrinsic occurrence rates were calculated i.e., sub-Neptunes and Neptunes (0.017-0.055 $\frac{R_\text{p}}{R_\star}$\ or 1.8-6\,$R_\oplus$, assuming a solar radius) within 12.5\,days (about half a TESS sector).}
    \label{fig:all_det_eff}
\end{figure}

\begin{figure}[!htb]
\minipage{0.46\textwidth}
  \includegraphics[width=\linewidth]{PiEri-small-highres-binned.png}
\endminipage\hfill
\minipage{0.46\textwidth}
  \includegraphics[width=\linewidth]{IC2602-small-highres-binned.png}
\endminipage\hfill
\minipage{0.46\textwidth}
  \includegraphics[width=\linewidth]{IC2602-small-highres_uncomtam-binned.png}
\endminipage
\caption{A comparison of extreme detection efficiencies to highlight the effect of flux contamination: \textbf{Top: }Detection efficiency of Pisces-Eridani, a moving group with members that are spread out so flux contamination is minimal. \textbf{Middle: }Detection efficiency of IC~2602, a crowded cluster, while taking into account flux contamination while injecting planets. \textbf{Bottom: }Detection efficiency of IC~2602 while assuming no flux contamination while injecting planets.}
\label{fig:det_eff_comp}
\end{figure}

\subsection{Preliminary Calculation of Occurrence Rates}
The intrinsic occurrence rate of planets ($\eta$) can be calculated from the fraction of stars with detected planets in a survey and the survey completeness. Here, we used the inverse detection efficiency method as follows:
\begin{equation}
\eta = \frac{{n_\text{p}}}{n_\star} * \frac{1}{comp_\text{bin}}
\end{equation}
where $comp_\text{bin}$ is the survey completeness evaluated in the defined bin, $n_p$ is the number of detected planets in the bin and $n_\star$ is the number of surveyed stars. The survey completeness is computed by combining the detection efficiency and the geometric transit probability, which is given by
\begin{equation}
f_\text{geo} = \frac{R_\star}{a}
\end{equation}
Here, $R_\star$ is the stellar radius, and $a$ is the average semi-major axis which is calculated using \emph{Kepler}\xspace's third law. The uncertainty on the occurrence rate was calculated from the square root of the number of detected planets in the bin, assuming Poissonian statistics.

In our analysis, we only considered sub-Neptunes and Neptunes (0.017-0.055 $\frac{R_\text{p}}{R_\star}$\ or 1.8-6\,$R_\oplus$, assuming a solar radius) within 12.5\,days (about half a TESS sector)\footnote{Here, we use the $\frac{R_\text{p}}{R_\star}$ fit by \texttt{TLS}\xspace and not the values from \texttt{EXOTIC}\xspace since \texttt{EXOTIC}\xspace is not part of \texttt{pterodactyls}\xspace and is only used to fit for the planetary parameters once the POI is vetted by \texttt{triceratops}\xspace.} because we wanted to understand the primordial population of sub-Neptunes and Neptunes before being stripped off of their atmospheres. With an average detection efficiency of 9\% and a $f_\text{geo}$ of 0.1 (at a geometric mean orbital period of 3.5\,days and assuming a Solar-type star), we compute an occurrence rate of 49$\pm$20\% for our sample of young clusters.

We then use \texttt{epos}\xspace \citep{mulders2018exoplanet} to calculate the \emph{Kepler}\xspace Gyr-old FGK (Sun-like) occurrence in the same bin. \texttt{epos}\xspace{} is a well documented and tested Python code developed by our team\footnote{\url{http://eos-nexus.org/epos/}}, it is available on GitHub\footnote{\url{https://github.com/ GijsMulders/epos}}, and has already been used in several publications to compute occurrence rates as well as to compare planet formation models to the \emph{Kepler}\xspace exoplanetary systems \citep[e.g.,][]{kopparapu2018exoplanet,Pascucci2018,fernandes2019hints,Mulders2019,Pascucci2019,mulders2020earths}. Here, we computed an occurrence rate value of 6.8$\pm$0.3\%, which is much lower than our computed value for young clusters. This could be attributed to several reasons, the major one being that our sample is heavily biased towards clusters with known planets. For example, if we increased the sample of young stars to search by a factor of six, which is about all nearby young stars within $\sim$200\,pc and assumed the same detection efficiency with no new planets, the occurrence would drop to $\sim$8\% which is the same as the \emph{Kepler}\xspace value. Furthermore, since most of the stars in the sample do not have stellar radii, masses or effective temperatures, it is very likely that we are considering non-FGK stars, most likely lower mass stars \citep{henry2006solar,winters2014solar}, while calculating the occurrence rate which would make the comparison to \emph{Kepler}\xspace inaccurate, where the sample is mostly Sun-like stars.  Therefore, in order to better estimate the occurrence rates of planets in young clusters, we will need to expand our search to all nearby clusters, measure the stellar parameters for those stars, and only compute the occurrence for FGK stars in our sample. However, if we continue to see this increased occurrence after we expand our sample to include all nearby clusters and moving groups, it could imply that there is an excess of larger planets at young ages possibly because their atmospheres have not yet been stripped due to photoevaporation and/or core-powered mass loss.


\section{Summary and Future Work}\label{sec:summary}
In this paper, we present \texttt{pterodactyls}\xspace $-$ a Python-based pipeline that is built using publicly available codes, and is specifically designed to extract, detrend, search for, and vet transiting exoplanets in young clusters using TESS' Primary Mission 30-min cadence FFIs. Here, we only consider 5 specific young clusters and moving groups: Tucana-Horologium Association, IC~2602, Upper Centaurus Lupus, Ursa Major, and Pisces Eridani. These clusters were chosen solely due to the presence of eight confirmed planets and two candidates which were used to validate \texttt{pterodactyls}\xspace. Using \texttt{pterodactyls}\xspace,

\begin{itemize}[leftmargin=*]
    \item We were able to recover eight out of ten known transiting planets/candidates. For the two that we were not able to recover (both initially discovered using 2-min cadence data), one planet was too small to detect (PiEri; TOI~451\,b) and the other is a long period, single transit planet (UCL; HIP~67522\,c), which is beyond the capabilities of the pipeline.
    \item We also detected 15 transiting signals, all of which were found to be either eclipsing binaries or nearby eclipsing binaries when their phase-folded light curves were vetted by \texttt{triceratops}\xspace.
    \item While taking into account the stellar flux contribution, we conducted injection-recovery tests in order to measure the detection efficiency of the \texttt{pterodactyls}\xspace. We found that the overall detection efficiency does not exceed 50\% for any bin below 0.1 $\frac{R_\text{p}}{R_\star}$ (i.e., a Jupiter-sized planet orbiting a Sun-like star) and that all the recovered planets are in bins with fairly low values ($<$20\%), hinting that even in young clusters smaller planets are more frequent than giant planets.
    \item For our biased sample, we computed a planet occurrence rate of 49$\pm$20\% which is significantly higher than \emph{Kepler}\xspace's Gyr-old occurrence rates of 6.8$\pm$0.3\% for sub-Neptunes and Neptunes (0.017-0.055 $\frac{R_\text{p}}{R_\star}$\ or 1.8-6\,$R_\oplus$, assuming a Solar radius) within 12.5\,days (about half a TESS sector).
\end{itemize}

Since most of these planets were initially detected in 2-min cadence data, our work shows that 30-min cadence data can be used to detect young transiting planets and thus has the potential for more planet discoveries. The planet occurrence rate we calculated for our sample is higher than that for \emph{Kepler}\xspace's Gyr old stars. However, we realize that this value is heavily biased since we only consider clusters with confirmed/candidate planets and can be further improved upon by using \texttt{pterodactyls}\xspace to search for young ($< 1$\,Gyr) transiting planets in a larger sample of nearby clusters and moving groups ($<$200\,pc) using TESS Primary and Extended Mission FFIs. We also plan on uniformly characterizing all the stars in our sample, so as to only consider FGK stars in future studies. With this young population of transiting, short-period planets, we hope to establish how the radius distribution of transiting exoplanets evolved over time, and therefore provide observational constraints on the mass loss mechanisms of planetary atmospheres.






\software{\texttt{NumPy} \citep{numpy}, \texttt{SciPy} \citep{scipy}, \texttt{Matplotlib} \citep{pyplot}, \texttt{eleanor}\xspace{} \citep{feinstein2019eleanor}, \texttt{W\={o}tan}\xspace{} \citep{hippke2019wotan}, \texttt{transitleastsquares} \citep{hippke2019optimized}, \texttt{vetting} \citep{hedges2021vetting}, \texttt{triceratops}\xspace{} \citep{giacalone2020vetting}, \texttt{EDI-Vetter Unplugged} \citep{zink2019edivetter}, \texttt{EXOTIC}\xspace \citep{zellem2020utilizing}, \texttt{epos}\xspace \citep{Mulders2018}}


\section{Acknowledgements}
R.B.F. would like to thank the following individuals for their expertise, assistance and, invaluable insights throughout the testing of \texttt{pterodactyls}\xspace: Adina D. Feinstein and Benjamin T. Montet (\texttt{eleanor}\xspace), Michael Hippke (\texttt{W\={o}tan}\xspace and \texttt{TLS}\xspace), Robert T. Zellem (\texttt{EXOTIC}\xspace), Christina Hedges (\texttt{vetting}, centroid test), and the Scaling \emph{K2}\xspace team i.e., Jessie L. Christiansen, Sakhee Bhure and Britt Duffy Adkins. G.D.M. acknowledges support from ANID --- Millennium Science Initiative --- ICN12\_009. This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by the NASA's Science Mission Directorate. This material is based upon work supported by the National Aeronautics and Space Administration (NASA) under Agreement No. NNX15AD94G for the program Earths in Other Solar Systems, under Agreement No. 80NSSC21K0593 for the program “Alien Earths”, and Grant No. 80NSSC20K0446 issued through the Astrophysics Data Analysis Program (ADAP). The results reported herein benefited from collaborations and/or information exchange within NASA’s Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA’s Science Mission Directorate. This publication makes use of data products from Exoplanet Watch, a citizen science project managed by NASA’s Jet Propulsion Laboratory on behalf of NASA’s Universe of Learning. This work is supported by NASA under award number NNX16AC65A to the Space Telescope Science Institute.
\clearpage

\bibliographystyle{apj}


\subsubsection*{#1}}

\pagestyle{headings}
\markright{Reference sheet: \texttt{natbib}}

\usepackage{shortvrb}
\MakeShortVerb{\|}

\begin{document}
\thispagestyle{plain}

\newcommand{\textsc{Bib}\TeX}{\textsc{Bib}\TeX}
\newcommand{\texttt{#1}\def\filedate{#2}\def\fileversion{#3}}}{\texttt{#1}\def\filedate{#2}\def\fileversion{#3}}}
\begin{center}{\bfseries\Large
    Reference sheet for \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ usage}\\
    \large(Describing version \fileversion\ from \filedate)
\end{center}

\begin{quote}\slshape
For a more detailed description of the \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ package, \LaTeX\ the
source file \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\texttt{.dtx}.
\end{quote}
\head{Overview}
The \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ package is a reimplementation of the \LaTeX\ |\cite| command,
to work with both author--year and numerical citations. It is compatible with
the standard bibliographic style files, such as \texttt{plain.bst}, as well as
with those for \texttt{harvard}, \texttt{apalike}, \texttt{chicago},
\texttt{astron}, \texttt{authordate}, and of course \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}.

\head{Loading}
Load with |\usepackage[|\emph{options}|]{|\texttt{#1}\def\filedate{#2}\def\fileversion{#3}}|}|. See list of
\emph{options} at the end.

\head{Replacement bibliography styles}
I provide three new \texttt{.bst} files to replace the standard \LaTeX\
numerical ones:
\begin{quote}\ttfamily
 plainnat.bst \qquad abbrvnat.bst \qquad unsrtnat.bst
\end{quote}
\head{Basic commands}
The \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ package has two basic citation commands, |\citet| and
|\citep| for \emph{textual} and \emph{parenthetical} citations, respectively.
There also exist the starred versions |\citet*| and |\citep*| that print
the full author list, and not just the abbreviated one.
All of these may take one or two optional arguments to add some text before
and after the citation.
\begin{quote}
\begin{tabular}{l@{\quad$\Rightarrow$\quad}l}
  |\citet{jon90}| & Jones et al. (1990)\\
  |\citet[chap.~2]{jon90}| & Jones et al. (1990, chap.~2)\\[0.5ex]
  |\citep{jon90}| & (Jones et al., 1990)\\
  |\citep[chap.~2]{jon90}| & (Jones et al., 1990, chap.~2)\\
  |\citep[see][]{jon90}| & (see Jones et al., 1990)\\
  |\citep[see][chap.~2]{jon90}| & (see Jones et al., 1990, chap.~2)\\[0.5ex]
  |\citet*{jon90}| & Jones, Baker, and Williams (1990)\\
  |\citep*{jon90}| & (Jones, Baker, and Williams, 1990)
\end{tabular}
\end{quote}
\head{Multiple citations}
Multiple citations may be made by including more than one
citation key in the |\cite| command argument.
\begin{quote}
\begin{tabular}{l@{\quad$\Rightarrow$\quad}l}
  |\citet{jon90,jam91}| & Jones et al. (1990); James et al. (1991)\\
  |\citep{jon90,jam91}| & (Jones et al., 1990; James et al. 1991)\\
  |\citep{jon90,jon91}| & (Jones et al., 1990, 1991)\\
  |\citep{jon90a,jon90b}| & (Jones et al., 1990a,b)
\end{tabular}
\end{quote}

\head{Numerical mode}
These examples are for author--year citation mode. In numerical mode, the
results are different.
\begin{quote}
\begin{tabular}{l@{\quad$\Rightarrow$\quad}l}
  |\citet{jon90}| & Jones et al. [21]\\
  |\citet[chap.~2]{jon90}| & Jones et al. [21, chap.~2]\\[0.5ex]
  |\citep{jon90}| & [21]\\
  |\citep[chap.~2]{jon90}| & [21, chap.~2]\\
  |\citep[see][]{jon90}| & [see 21]\\
  |\citep[see][chap.~2]{jon90}| & [see 21, chap.~2]\\[0.5ex]
  |\citep{jon90a,jon90b}| & [21, 32]
\end{tabular}
\end{quote}
\head{Suppressed parentheses}
As an alternative form of citation, |\citealt| is the same as |\citet| but
\emph{without parentheses}. Similarly, |\citealp| is |\citep| without
parentheses. Multiple references, notes, and the starred variants
also exist.
\begin{quote}
\begin{tabular}{l@{\quad$\Rightarrow$\quad}l}
  |\citealt{jon90}| & Jones et al.\ 1990\\
  |\citealt*{jon90}| & Jones, Baker, and Williams 1990\\
  |\citealp{jon90}| & Jones et al., 1990\\
  |\citealp*{jon90}| & Jones, Baker, and Williams, 1990\\
  |\citealp{jon90,jam91}| & Jones et al., 1990; James et al., 1991\\
  |\citealp[pg.~32]{jon90}| & Jones et al., 1990, pg.~32\\
  |\citetext{priv.\ comm.}| & (priv.\ comm.)
\end{tabular}
\end{quote}
The |\citetext| command
allows arbitrary text to be placed in the current citation parentheses.
This may be used in combination with |\citealp|.
\head{Partial citations}
In author--year schemes, it is sometimes desirable to be able to refer to
the authors without the year, or vice versa. This is provided with the
extra commands
\begin{quote}
\begin{tabular}{l@{\quad$\Rightarrow$\quad}l}
  |\citeauthor{jon90}| & Jones et al.\\
  |\citeauthor*{jon90}| & Jones, Baker, and Williams\\
  |\citeyear{jon90}|   & 1990\\
  |\citeyearpar{jon90}| & (1990)
\end{tabular}
\end{quote}
\head{Forcing upper cased names}
If the first author's name contains a \textsl{von} part, such as ``della
Robbia'', then |\citet{dRob98}| produces ``della Robbia (1998)'', even at the
beginning of a sentence. One can force the first letter to be in upper case
with the command |\Citet| instead. Other upper case commands also exist.
\begin{quote}
\begin{tabular}{rl@{\quad$\Rightarrow$\quad}l}
  when & |\citet{dRob98}| & della Robbia (1998) \\
  then & |\Citet{dRob98}| & Della Robbia (1998) \\
     &   |\Citep{dRob98}| & (Della Robbia, 1998) \\
     &   |\Citealt{dRob98}| & Della Robbia 1998 \\
     &   |\Citealp{dRob98}| & Della Robbia, 1998 \\
     &   |\Citeauthor{dRob98}| & Della Robbia
\end{tabular}
\end{quote}
These commands also exist in starred versions for full author names.

\head{Citation aliasing}
Sometimes one wants to refer to a reference with a special designation,
rather than by the authors, i.e. as Paper~I, Paper~II. Such aliases can be
defined and used, textual and/or parenthetical with:
\begin{quote}
\begin{tabular}{lcl}
  |\defcitealias{jon90}{Paper~I}|\\
  |\citetalias{jon90}| & $\Rightarrow$ & Paper~I\\
  |\citepalias{jon90}| & $\Rightarrow$ & (Paper~I)
\end{tabular}
\end{quote}
These citation commands function much like |\citet| and |\citep|: they may
take multiple keys in the argument, may contain notes, and are marked as
hyperlinks.
\head{Selecting citation style and punctuation}
Use the command |\bibpunct| with one optional and 6 mandatory arguments:
\begin{enumerate}
\item the opening bracket symbol, default = (
\item the closing bracket symbol, default = )
\item the punctuation between multiple citations, default = ;
\item the letter `n' for numerical style, or `s' for numerical superscript
      style, any other letter for
      author--year, default = author--year;
\item the punctuation that comes between the author names and the year
\item the punctuation that comes between years or numbers when common author
      lists are suppressed (default = ,);
\end{enumerate}

The optional argument is the character preceding a post-note, default is a
comma plus space. In redefining this character, one must include a space if
one is wanted.

Example~1, |\bibpunct{[}{]}{,}{a}{}{;}| changes the output of
\begin{quote}
  |\citep{jon90,jon91,jam92}|
\end{quote}
into [Jones et al. 1990; 1991, James et al. 1992].

Example~2, |\bibpunct[; ]{(}{)}{,}{a}{}{;}| changes the output of
\begin{quote}
  |\citep[and references therein]{jon90}|
\end{quote}
into (Jones et al. 1990; and references therein).

\head{Other formatting options}
Redefine |\bibsection| to the desired sectioning command for introducing
the list of references. This is normally |\section*| or |\chapter*|.

Define |\bibpreamble| to be any text that is to be printed after the heading but
before the actual list of references.

Define |\bibfont| to be a font declaration, e.g.\ |\small| to apply to
the list of references.

Define |\citenumfont| to be a font declaration or command like |\itshape|
or |\textit|.

Redefine |\bibnumfmt| as a command with an argument to format the numbers in
the list of references. The default definition is |[#1]|.

The indentation after the first line of each reference is given by
|\bibhang|; change this with the |\setlength| command.

The vertical spacing between references is set by |\bibsep|; change this with
the |\setlength| command.

\head{Automatic indexing of citations}
If one wishes to have the citations entered in the \texttt{.idx} indexing
file, it is only necessary to issue |\citeindextrue| at any point in the
document. All following |\cite| commands, of all variations, then insert
the corresponding entry to that file. With |\citeindexfalse|, these
entries will no longer be made.

\head{Use with \texttt{chapterbib} package}

The \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ package is compatible with the \texttt{chapterbib} package
which makes it possible to have several bibliographies in one document.

The package makes use of the |\include| command, and each |\include|d file
has its own bibliography.

The order in which the \texttt{chapterbib} and \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ packages are loaded
is unimportant.

The \texttt{chapterbib} package provides an option \texttt{sectionbib}
that puts the bibliography in a |\section*| instead of |\chapter*|,
something that makes sense if there is a bibliography in each chapter.
This option will not work when \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\ is also loaded; instead, add
the option to \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}.

Every |\include|d file must contain its own
|\bibliography| command where the bibliography is to appear. The database
files listed as arguments to this command can be different in each file,
of course. However, what is not so obvious, is that each file must also
contain a |\bibliographystyle| command, \emph{preferably with the same
style  argument}.
\head{Sorting and compressing citations}
Do not use the \texttt{cite} package with \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}; rather use one of the
options \texttt{sort} or \texttt{sort\&compress}.

These also work with author--year citations, making multiple citations appear
in their order in the reference list.

\head{Long author list on first citation}
Use option \texttt{longnamesfirst} to have first citation automatically give
the full list of authors.

Suppress this for certain citations with |\shortcites{|\emph{key-list}|}|,
given before the first citation.

\head{Local configuration}
Any local recoding or definitions can be put in \texttt{#1}\def\filedate{#2}\def\fileversion{#3}}\texttt{.cfg} which
is read in after the main package file.

\head{Options that can be added to \texttt{\char`\\ usepackage}}
\begin{description}
\item[\ttfamily round] (default) for round parentheses;
\item[\ttfamily square] for square brackets;
\item[\ttfamily curly] for curly braces;
\item[\ttfamily angle] for angle brackets;
\item[\ttfamily colon] (default) to separate multiple citations with
     colons;
\item[\ttfamily comma] to use commas as separaters;
\item[\ttfamily authoryear] (default) for author--year citations;
\item[\ttfamily numbers] for numerical citations;
\item[\ttfamily super] for superscripted numerical citations, as in
     \textsl{Nature};
\item[\ttfamily sort] orders multiple citations into the sequence in
     which they appear in the list of references;
\item[\ttfamily sort\&compress] as \texttt{sort} but in addition multiple
     numerical citations are compressed if possible (as 3--6, 15);
\item[\ttfamily longnamesfirst] makes the first citation of any reference
     the equivalent of the starred variant (full author list) and subsequent
     citations normal (abbreviated list);
\item[\ttfamily sectionbib] redefines |\thebibliography| to issue
     |\section*| instead of |\chapter*|; valid only for classes with a
     |\chapter| command; to be used with the \texttt{chapterbib} package;
\item[\ttfamily nonamebreak] keeps all the authors' names in a citation on
     one line; causes overfull hboxes but helps with some \texttt{hyperref}
     problems.
\end{description}
\end{document}
