A survey for variable stars with small telescopes: IX -- Evolution of Spot Properties on YSOs in IC5070

We present spot properties on 32 periodic young stellar objects in IC 5070. Long term, $\sim$5 yr, light curves in the $V$, $R$, and $I$-bands are obtained through the HOYS (Hunting Outbursting Young Stars) citizen science project. These are dissected into six months long slices, with 3 months oversampling, to measure 234 sets of amplitudes in all filters. We fit 180 of these with reliable spot solutions. Two thirds of spot solutions are cold spots, the lowest is 2150 K below the stellar temperature. One third are warm spots that are above the stellar temperature by less than $\sim$2000 K. Cold and warm spots have maximum surface coverage values of 40 percent, although only 16 percent of warm spots are above 20 percent surface coverage as opposed to 60 percent of the cold spots. Warm spots are most likely caused by a combination of plages and low density accretion columns, most common on objects without inner disc excess emission in $K-W2$. Five small hot spot solutions have $<3$ percent coverage and are 3000 - 5000 K above the stellar temperature. These are attributed to accretion, and four of them occur on the same object. The majority of our objects are likely to be accreting. However, we observe very few accretion hot spots as either the accretion is not stable on our timescale or the photometry is dominated by other features. We do not identify cyclical spot behaviour on the targets. We additionally identify and discuss a number of objects that have interesting amplitudes, phase changes, or spot properties.


INTRODUCTION
Young Stellar Objects (YSOs) display variability on a wide range of timescales.Rotational modulation caused by surface spots occurs when spots rotate in and out of the observer's line of sight, on the timescale of the rotational period.The fast rotation of young stars (see reviews by Herbst et al. 2007;Bouvier et al. 2014), hence results in photometric variability on the order of days.The amplitude of variation produced by spots depends on the nature and the characteristics of the spot.
Cold spots are regions on the photosphere, analogous to sun spots, that are hundreds of degrees cooler than the stellar temperature.These cause flux variations of generally less than ten percent (e.g.Strassmeier 2009).Cold spots are expected to be present during all stages of the pre-main sequence evolution of stars, whereas hot spots are thought to be exclusive to accreting classical T Tauri (stage 2) stars with full or transition discs (Herbst et al. 1994).
Hot spots are considered to be footprints of accretion with temperatures that are several thousand degrees above the stellar temperature and coverage less than a few percent of the stellar surface (Hartmann et al. 2016;Muzerolle et al. 1998).Emission lines as an indicator of accretion require these high temperatures.However, recent surveys present a more diverse picture with low temperature contrast, high coverage regions (Scholz et al. 2009(Scholz et al. , 2012;;Bozhinova et al. 2016).Espaillat et al. (2021) found low density accretion columns as a possible explanation for these spots, with coverage values up to 20 percent, while Sicilia-Aguilar et al. (2023) found warm regions spatially distinct from the accretion shock on EX Lup and TW Hya.
Multi-filter photometry was used to identify surface spots on T Tauri stars in Bouvier & Bertout (1989), and has been used to study spot properties in many surveys since (e.g.Bouvier et al. 1993Bouvier et al. , 1995;;Carpenter et al. 2001).The strong, active magnetic fields of ★ E-mail: cbah2@kent.ac.uk † HOYS Observer ‡ Beacon Observer § Observadores de Supernovas, https://www.obsn.es/YSOs cause large, complex structures and rotational variation is produced by asymmetric spot distributions.Photometry is therefore limited when compared to the resolution of methods such as Doppler Imaging (Strassmeier 2002;Strassmeier et al. 2003;Donati & Landstreet 2009;Skelly et al. 2009).However, the comparatively low-cost equipment requirements for optical photometry allow for larger, longer surveys and better number statistics.
In Herbert et al. (2023) we developed a flux replacement method to identify surface spots using , , and  band data obtained from the Hunting Outbursting Young Stars (HOYS1 ) citizen science project from an 80 d high cadence section of light curves from objects in IC 5070.We identified 21 objects with cold spots, and six with warm spots with temperatures up to 3000 K above the stellar temperature.In this work we will apply the same method to a sample of YSOs in the same region to characterise their spots over the entire duration of their light curves since the project started in 2014.This will allow us to track the evolution of the spot properties over time and greatly improve our statistics of the distribution of the typical spot temperatures and coverage values on YSOs.
In this paper we detail our data in Sect. 2. In Sect.3 we explain our method to divide the light curves, identify periodic objects and determine amplitudes.The results for the period analysis are presented in Sect. 4. In Sect. 5 we discuss our results for the amplitude, phase, and spot property evolution.A number of objects with properties of particular interest are identified and discussed in Sect.6.

HOYS data and photometry
The photometry data for this work were obtained through the HOYS citizen science project (Froebrich et al. 2018).This uses a world wide network of ≈ 100 amateur, university, and professional telescopes to conduct long term photometric monitoring of 25 nearby ( < 1 kpc) young clusters and star forming regions in optical broadband filters.The aim of this project is to obtain photometry in all filters for all target regions with a cadence of 12 -24 hours over many years.
The data supplied to the project are obtained from a range of instruments.The photometry is calibrated against reference images taken under photometric conditions in , , ,   , and   filters.For simplicity we will refer to them as , , , , and , hereafter.The calibration offsets from instrumental to apparent magnitude for these reference images have been obtained from the Cambridge Photometric Calibration Server2 .To account for colour terms in the calibration due to instrumentation (e.g.slightly non-standard filters) or conditions (e.g.cirrus clouds), an internal calibration is applied to all data following Evitts et al. (2020).We identify non-variable stars in each target region, and use their known magnitudes and colours to determine and correct for the colour terms in each image.This colour calibration also evaluates the photometry uncertainty for each data point and is applied to all HOYS data used in this work.

IC 5070 cluster membership
The star forming region IC 5070 (Pelican Nebula) is part of the W 80 Hii region alongside NGC 7000 (North American Nebula).The entire region is frequently referred to as NAP for the combined region names.The NAP region was found to be at ∼ 796 ± 25 pc in Kuhn et al. (2020), which identified 395 YSOs in six astrometric groups in the region based on Gaia DR2.Almost all objects were found to be less than 3 Myr old, with a majority younger than 1 Myr.Using Gaia DR2 as well, Evitts et al. (2020) determined a distance to the IC 5070 region of 870 +70 −55 pc.In this current work we have identified potential IC 5070 cluster members using Gaia DR3 (Gaia Collaboration et al. 2016Collaboration et al. , 2023) ) astrometry.This is part of a larger work characterising the YSO populations in all HOYS fields (Froebrich et al. 2024).All Gaia sources within 0.6 degrees of the cluster centre were selected.Sources were removed that had parallax values of less than 0.3 mas, a signal-to-noise ration (SNR) of the parallax of less than five, and that were fainter than 18 mag in the Gaia -band.We further removed objects with colours indicative of white dwarfs and cataclysmic variables.Thus, sources with  −  < −0.2 mag, and  −  < 0.0 mag were excluded.Cluster members were initially manually identified through a histogram of the distances, using a bin width of 20 pc.Stars in a distance range of 750 pc to 900 pc were selected.Within that selection, two groups of stars with a coherent proper motion were identified.For each group the median distance and proper motions and their root mean square (rms) scatter were determined.Candidate cluster members were then selected as all stars that are within three standard deviations around these median values.The potential cluster members were assigned to the group they were closest to in proper motion space.A summary of the astrometric properties of the two groups can be found in Table 1.
There were a total of 366 potential cluster members based on the Gaia DR3 astrometry in the two groups of IC 5070.Of those 252 are in group a, the other 114 are in group b.Of these 366 cluster members, we have corresponding HOYS light curves for 131 (84, and 47 in groups a and b, respectively) with at least 100 photometry data points in each of the , , and  filters.The majority of the YSOs with HOYS light curves identified in Froebrich et al. (2021) Table 1.Properties of the two astrometric groups of young stars in the IC 5070 region, as determined in Froebrich et al. (2024).The columns contain the following information: Group: The group name; ,   ,   : median distance, standard error of the median, rms scatter from the median of the cluster members;  /  ,   /  ,   /  : median proper motion in RA/DEC, standard error of the median, rms scatter from the median of the cluster members; N: Number of Gaia DR3 selected potential members; using Gaia EDR3, were again identified as cluster members.Seven of the objects previously analysed are not included in our list.This is due to a parallax or proper motion value in Gaia DR3 that is beyond the selection range applied.Two objects that were not previously investigated are now included in our sample.65 of the 366 potential cluster members have recently been identified as H  emission line stars in Panwar et al. (2023).The overlap between the samples will be examined in more detail in Sect.5.4.

METHODOLOGY
In this section we detail the data analysis used in the paper.To a large extend this work follows our previous work on period identification in HOYS data in Froebrich et al. (2021) and spot property determination from peak to peak amplitude measurements in Herbert et al. (2023).We detail the basic principles again below and highlight changes made to our established procedures for the purpose of this work.

Sliced Light Curves
Our aim is to analyse the properties of surface spots over time.HOYS observations in this field begin in 2014, and the light curves were extracted from our database on October 21st, 2022.We have dissected (sliced) the , , and  light curves into blocks of six months in length, every three months.Thus, consecutive slices of light curve data overlap by three months.The starting date for the first slice was chosen as Feb 14, 2014, 12 noon UT.The following slices are numbered consecutively by integers.The start date for the first slice (indexed as zero) is in the first year of HOYS data available for this field, and this day of the year corresponds to the time of worst observability for the IC 5070 field, given its Right Ascension (i.e. it is at its lowest altitude at midnight).Thus, slices with numbers 3, 7, 11, etc. are centered on February, and hence usually contain the least amount of data.On the other hand, slices numbered 1, 5, 9, etc. are centred on mid August each year and thus generally have the most photometry data points.
For our analyses we require a minimum of photometry data points in each slice.We need at least 50 data points in at least two filters to determine the possible period (see Sect. 3.2), and at least 50 data points in each filter to measure the peak-to-peak amplitudes (see Sect. 3.3).Thus, the first slice with sufficient data for the majority of objects is centred on mid May 2018 (MJD = 58162.5),i.e. slice number 16.The final slice with sufficient data in all filters is slice 33, therefore at most there are 17 slices with sufficient data for spot fitting.Hence, our light curves for the region are typically ∼5 yr in length.The observational cadence of an object, and thus the number of photometry data points available, is affected by its brightness and distance from the centre of the HOYS field.The majority of images in the HOYS database come from amateur telescopes with varying fields of view and limiting magnitude.Thus, generally fainter objects, further from the central part of the field have less data points available for analysis.Note that in all figures where properties are shown as function of time, the data points are plotted at the start date of their respective slice.

Period Identification
Identifying periodic objects in IC 5070 from HOYS data was the focus of study in Froebrich et al. (2021).A double-blind study of period-finding algorithms was conducted on data from an 80 day period in the summer of 2018 (MJD ranging from 58329.5 to 58409.5).During this time, the IC 5070 field was the focus of a HOYS observing campaign, leading to a high cadence data set.Photometry data from this period appears in our slices 17 and 18.The campaign was repeated every August in subsequent years, but for the majority of stars this time remains the part of the light curve data with the highest cadence.
The Froebrich et al. (2021) study concluded that a combination of four period finding algorithms is optimal for completeness in identifying periodic objects in HOYS data.The four methods that were identified are: L1Boot, L1Beta, L2Beta, and the generalised Lomb-Scargle (GLS) periodogram.See the detailed description in Froebrich et al. (2021) and references therein for an in depth discussion of the specifics of each of the methods.The L1Boot, L1Beta, and L2Beta methods are run using the R module RobPer (Thieler et al. 2016).The GLS periodogram is based on Scargle (1982); Zechmeister & Kürster (2009) and we are using its implementation in astropy.
For the purpose of this work, we are making some small changes to the methodology applied in Froebrich et al. (2021).First of all, in that work periodic light curves were identified for all stars in the IC 5070 field.Here, we limit our search to the light curves of stars identified as potential cluster members based on Gaia DR3 astrometry (see Sect. 2.2).Furthermore, we search for significant periods in each of the light curve slices independently.This ensures that for objects which show periodic light curve modulation only at certain times, e.g.due to varying surface spot properties, a period is identified.Thus, a more complete sample of rotational variables of the region is obtained.We made some further minor adjustments to the selection of the best period, compared to Froebrich et al. (2021), which are detailed below.
The four period finding algorithms (L1Boot, L1Beta, L2Beta, and GLS) were applied to the photometry in each slice and filter (, , and ), separately.Only slices which had at least 50 data points in two of these filters were considered.For all four period finding methods periodograms with 1000 test periods were obtained.These test periods were sampled homogeneously in frequency space between 0.5 d and 20 d.This range was chosen as most appropriate based on the typical cadence of our data and the periods identified in Froebrich et al. (2021).
In each periodogram the highest peak was chosen as the candidate period.To ensure periods were real, only candidate periods with a peak periodogram power above 0.15 were considered (for L1Boot, L1Beta and L2Beta).This differs from Froebrich et al. (2021), where a slightly higher threshold of 0.2 was used.The reduction in the detection threshold is justified by the larger amount of data available in the current analysis.Thus, there is more opportunity to identify and remove questionable periods (see below).The GLS equivalent to this peak periodogram power is the False Alarm Probability (FAP).We accepted all periods with a FAP below 0.1 as candidate periods.Despite the geographically distributed nature of the HOYS observations, the typical cadence in the data is approximately one day.Hence, candidate periods identified within one percent of one day were removed.We also removed any candidate periods that were within one percent of 0.5 d and 20 d as these are near the edges of the period parameter space we investigated.
If in a given slice a period was detected matching the above criteria and in at least two filters with a separation of less than 0.05 d, then it was retained as a candidate period for this slice.This generated a list of such candidate periods for each object.In some objects these lists contained periods that were not all identical within the uncertainties.Thus we employed the following methodology to select the most likely correct period in all cases.A histogram of all candidate periods for each object was created with a bin-width of 0.1 d for the entire range of investigated periods (0.5 to 20 d).This resulted in groups of candidate periods in individual or immediately adjacent bins, and distinct groups being separated by one or more bins with no candidate periods within them.Typically these groups only span two adjacent bins (i.e.0.2 d) and the maximum range was 0.4 d.In all cases, the scatter of the periods within each group was below ten percent of the median.
For most objects one of these groups contained the vast majority of candidate periods.We visually inspected phase folded light curves determined with the median period for each group.We found that for phase plots from groups with a small number of periods no clear systematic variability could be visually identified.Thus, only groups with more than five candidate periods were kept.For most objects this left only one single group of candidate periods.In all other cases we selected a dominant group if it had more than 2/3rd of the candidate periods for this object.There were only two objects (18,26) where this was not possible and two almost equal sized groups with different periods remained.In those cases we visually inspected the phase folded light curves and manually selected the most likely group with the most likely correct median period.In both cases the group with the shorter candidate periods was chosen.
Finally, again following the procedure in Froebrich et al. (2021), we determined the final period for each object in the following way.For each filter, a simple Lomb-Scargle periodogram with 1000 homogeneously distributed test periods was determined for the entire light curve within ten percent of the median of the candidate periods in the dominant group.This resulted in a much higher period resolution.The peak in this periodogram was used as the most likely period of this object in each filter.The median of these periods from the , , and  data has been adopted as the final period of the object.Similarly the standard deviation of these three values is used as the final period uncertainty for the object throughout our analysis.These values are listed in Table 2 for all objects.

Peak to peak amplitude identification
In order to analyse the spot property evolution we require peak to peak amplitudes in , , and  over time.This basically follows the procedure established in Herbert et al. (2023).We also use the same notation as in that work.For each periodic object the photometry data slices that were created in Sec.3.1 were phase folded with the period determined in Sec.3.2.A minimum of 50 data points per slice and filter were required for this.The running median in the phase plot was calculated, using a smoothing over 0.1 in phase.The threshold of 50 data points is equivalent to a data point every 3.6 days in the slice to ensure a high signal to noise in the measured peak to peak amplitudes.
Following Herbert et al. (2023) we denote the peak to peak amplitude as Â  , where '' indicates the amplitudes as observed, and  refers to the filter used.It is determined as the difference between the maximum and minimum of the smoothed running median brightness in the phase plot.The associated uncertainties  Â  are determined as the standard error of the mean of the data used to calculate the maximum and minimum median.We also determine the phase position of the maximum of the phase plot in each filter.
In Herbert et al. (2023) we found that a signal to noise ratio (SNR) above three for the peak to peak amplitudes resulted in determined spot properties (spot temperature and coverage) with smaller uncertainties.It further also resulted in more consistent results when adding data from shorter wavelengths into the analysis.Therefore, we discarded all peak to peak amplitudes in slices in which one or more of them had a SNR below three.
Furthermore, we introduce an additional selection criterion considering the phase of the folded light curve.All peak to peak amplitudes in a particular slice were discarded when the phase position of the maximum or minimum in the folded light curve scattered by more than 45 degrees between the three filters.This removed any slices where there was no consistent periodicity between filters.This mostly occurred where the SNR of the peak to peak amplitudes were just above the threshold of three.The phase of each slice shown later on, is calculated relative to a phase zero point.This has been calculated as the median phase value of the maxima in the phase folded light curves across all slices and filters.

Spot property fitting
In Herbert et al. (2023) we established and tested our method for determining spot properties from Â  values.We apply this method to every slice that contains over 50 data points in , , and .We assume that the Â  in each six-month slice is generated by a single, uniform spot with temperature   on an otherwise unspotted surface with temperature  ★ .We assume that the spot rotates completely in and out of the line of sight and therefore consider the spot to be 'in-front' or 'behind' of the star.Using a flux replacement model with PHOENIX synthetic spectra, we model for a given stellar temperature  ★ , 10 6 spots homogeneously sampled in the spot temperature   range of 2000 K ≤ T S ≤ 10000 K and spot coverage  range of the visible surface, 0 ≤  ≤ 0.5.The spectra are convolved with the , , and  filter transmission curves, using the speclite.filters3package.The magnitude difference in each filter between the spotted and un-spotted surface is taken to be the modelled peak to peak amplitude, Â  .We follow the same notation laid out in Herbert et al. (2023).The amplitudes are used together as a set { Â , Â , Â }, and we define that subscript { } denotes the shortest wavelength filter and subsequent filters.In this way, Â{ } = { Â , Â , Â }, which is used throughout our analysis.
The best fitting model is established by finding the minimum  between the modelled amplitude sets Â { } and the observed amplitudes Â { } .This single result will potentially be sensitive to small perturbations in the input peak to peak amplitudes.Thus, the Â { } values are varied within their associated uncertainties for 10000 iterations, and the best fitting models are found for each of them.In the majority of cases there are cold spot and hot spot solutions for the Â { } variations.To remove cases where there is ambiguity between the solutions, we require at least 60% of the 10000 variations to be either hot or cold spot solutions.If such a 60% majority is reached, the median and median absolute deviation (MAD) of the majority are taken as the spot properties and associated uncertainties.The ratio of hot spot to cold spot solutions HS:CS { } , is used as an indicator of reliability for the solution.
The spot fitting requires a known stellar temperature for each star.These were taken, if available, from Fang et al. (2020).The stellar temperatures are rounded to the nearest 50K for the spot fitting.The effect on the spot properties by altering the stellar temperature like this is always less than the statistical uncertainties, for details see Herbert et al. (2023).Where there is no stellar temperature reference available the median temperature of the astrometric group (a or b) the star belongs to (as discussed in Sec.2.2) was used.This is functionally 4100 K for group a, and 3950 K for group b.In Table 2 we identify the astrometric group each star belongs to.The exception to this is object 2, which is significantly brighter at 11 mag in the RP Gaia band.This was identified previously in Herbert et al. (2023) where an effective temperature of 5000 K was used to account for its brightness.In Herbert et al. (2023) we test the effect of changing the stellar temperature with simulated spots.We found that when the stellar temperature was altered by ±200 K, the resultant spot properties varied within the statistical uncertainties.

Rotation Period Distribution
Our data set contains light curves in , , and  for 131 objects.When divided into the six-month slices, 126 objects had sufficient data in one or more slices to search for a period.For 68 objects we detected a period in at least one slice.Of those, 31 objects meet the criteria laid out in Sect.3.2 to determine a final period.We also include an additional object in our sample (19, separate in Table 2) that was removed due to the vicinity of the object to a naked eye star (about three arcminutes away) and the unreliability of the photometry for such objects (see Froebrich et al. (2024)).In the case of object 19, the photometry in the  band is affected, but a good quality period was determined prior to removal from  and  band data, and so is included as part of our period discussion but is not followed up on in spot properties.
From the 32 objects, 22 overlap with the sample published in Froebrich et al. (2021).The objects where we could not determine a final period were also not included in Froebrich et al. (2021), with the exception of V1701 Cyg.This source will be discussed in more detail in Sect.6.1.The Froebrich et al. (2021) sample contained 40 YSOs.Objects in that work may not be included in our sample because of the change to Gaia DR3 astrometry, or because we did not have sufficient available  band data for the object.Our requirement for period searching was the availability of over 50 data points in at least two filters.Thus, we did not investigate objects that had little to no -band coverage as they would not be suitable for our spot fitting analysis.
The final periods and their uncertainties are listed in Table 2 together with other source properties.In this table we list our adopted source ID number (based on the sorted order of their Gaia DR3 ID), the J2000 positions, proper motions, and parallax, the astrometric group they are a member of, the effective temperature used for the spot fitting (from Fang et al. (2020) and discussed in Sect. 3.4), the  Table 2. Target list of all YSOs investigated in this work.For each object we list our ID number, the J2000 coordinates, proper motions, and parallax from Gaia DR3, the astrometric group it belong to, the effective temperature used for spot fitting, the determined period and its uncertainty, its  − 2 colour (discussed in Sect.5.4), the Gaia DR3 ID number, and the ID number of the object in Froebrich et al. (2021).Objects marked with (*) have no effective temperature reference in Fang et al. (2020) and the median of all astrometric group members were used or 5000 K for bright objects (see details in text).final periods and their uncertainties, their  − 2 colour (discussed in Sect.5.4), their Gaia DR3 ID, and the ID numbers for overlapping sources from Froebrich et al. (2021).We show the period distribution of our sample in Fig. 1.In blue a histogram of the periods, with 1 d bins, is displayed.A cumulative distribution function (CDF) is over plotted as a solid black line.The vertical dashed line at  = 5.5 d separates fast and slow rotators.We adopted the same period value for this as in Froebrich et al. (2021) for ease of comparison.The period distribution of our sample is bimodal.This bimodality has been identified in many other samples of YSOs, such as e.g. in ONC, NGC 2264, and the Perseus molecular cloud (Rodríguez-Ledesma et al. 2010;Herbst et al. 2007;Wang et al. 2023).The cause of this bimodality is attributed to star-disc interactions via magnetic fields.Long rotation periods are hence indicating a disc slowing the rotation rate.Conversely, once the disc is dissipated or looses the magnetic connection to the star, the star begins to 'spin up', leading to fast rotators.Of our objects, 23 are found to have  < 5.5 d and 9 have  > 5.5 d.The distribution roughly peaks at ∼3 d and ∼7 d, with a gap in-between.This is, within the statistical noise, comparable to Froebrich et al. (2021), where we found 25 fast and 15 slow rotators.Thus, roughly two thirds of the YSOs in IC 5070 are fast rotators.In Froebrich et al. (2021) a gap in the period distribution in IC 5070 between five and six days was identified.Here we find that Object 20 is the first source that populates this gap with a period of  = 5.52 d.This  object is hence in the transition between the slow and fast rotators in its evolution.

Comparison with literature periods
In Table 2 we show the ID numbers from the Froebrich et al. ( 2021) study for the objects that were previously identified as periodic.The ten objects that do not have such an ID are hence newly identified rotational variables.They were investigated in Froebrich et al. (2021) but no period was found.Note that this earlier work focused only on an 80 d duration, high cadence part of the HOYS IC 5070 light curves in 2018.These data are roughly in the middle of our slice 17, but that slice contains all the data for the 6-month period from MJD = 58162.5− 58345.5.For five of the new periodic sources (object 5, 15, 17, 23, 28) we also do not find any period in slice 17.But they have been identified as periodic in other slices.Hence, most likely the amplitudes of their variability were too low during that time for the periodic signal to be detectable.For the remaining five sources (object 12, 18, 19, 20, 31) we do identify a period in the slice 17 data.This can have two reasons: i) The object has shown slightly stronger (higher amplitude) periodic variability outside the 80-day window investigated in Froebrich et al. (2021), which is included in our slice 17 data.ii) Our slight lowering of the detection threshold in the periodograms (see Sect. 3.2) and the ability to confirm these periods in many other slices.
The periods identified here do not all lie within the uncertainties given in Froebrich et al. (2021).However, the median separation between our periods and the ones in Froebrich et al. ( 2021) is 0.27 percent, with a maximum of 1.45 percent.Thus, there are no significant changes in the determined periods between the two works.
In our work we have used the entire light curve, and are hence able to verify the period with multiple slices.Therefore, the periods we identify here are more accurate/reliable.Bhardwaj et al. (2019) also investigated variability of YSOs in IC 5070.Six objects from their work overlap with objects previously identified in Froebrich et al. (2021), where all but one agree on the period.Our work introduces three additional overlapping objects, namely objects 20, 23, and 28.Object 20 was identified in Bhardwaj et al. (2019) as variable but not periodic.Object 23 was identified in Bhardwaj et al. (2019) with a period of 8.446 d, as opposed to our 8.3115 ± 0.089 d.The earlier value is slightly outside our uncertainty range.Given our longer data set and the identification of this period in seven slices, our period is more reliable.Most noteworthy is object 28 which was identified in Bhardwaj et al. (2019) with a period of 0.518 d.This object, as shown in Table 2, has the longest period in our sample at  = 13.690± 0.262 d.It also has the highest absolute period uncertainty in the sample (object 16 has the highest relative period uncertainty).This is due to the best period measurement from the  band periodogram being quite far from the results using the  and -Band periodograms (see Sect. 3.2).The amplitudes in  and  for this source are much higher than for the  band, which can explain this difference.Furthermore, in our period search we exclude periods within 1% of 0.5 d, as these are on the edge of the investigated parameter space.In order to confirm we did not exclude a period close to 0.5 d as in Bhardwaj et al. (2019), we investigated the original identified candidate periods, before we made these exclusions.Object 28 had no other significant candidate periods outside of the primary group of 15 detections.That group had a median period of 13.620 d, and a standard deviation of 0.098 d.This implies that our period for this source is reliable and the nominal uncertainty is lower than the quoted ±0.2618 d.

Gaia Colour Magnitude Diagram
In Fig. 2 we show the Gaia  −  versus  colour magnitude diagram (based on Gaia DR3) for all 336 potential IC 5070 cluster members (groups a and b).The periodic objects in the sample are overlaid as coloured circles.The symbol size is proportional to the period and the colour indicates the  eff used in the spot property determination (as in Table 2).The periodic variables show the same distribution as the other potential cluster members.We find no significant trend between period and colour.They cover a colour range of 1.8 mag <  −  < 3.0 mag, and have a lower limit of  ∼ 17 mag.This is due to fainter sources having more noisy HOYS photometry, not allowing the detection of low amplitude periodicity.The brightest source is object 2 at  ∼ 12 mag.As discussed in Sec.3.3 we use a temperature of 5000 K for the surface temperature of this source.Object 30 is the hottest object with a temperature of 5500 K.However, it is also the reddest with  −  ∼ 3.0 mag.Fang et al. (2020) provide an extinction of   = 5.0mag for this source.In Herbert et al. (2023) our sample contained three objects similar to this.However, in the updated Gaia DR3 cluster members list, object 30 is the only object that remains.We consider the potential that the object is a contaminant, but it has  = 803.5 pc,   = −1.501,  = −2.871which places it within group a of the cluster.It has a short period with  = 1.4327 ± 0.0001 d, which indicates it has lost its connection with its disc.The  − 2 colour of 0.23 mag indicates a lack of inner disc material.No 3 − 4 colour is available for the object due to a poor signal to noise.
A digitized sky survey (DSS) image of the region is shown in Fig. 3.The positions of all identified periodic objects are over plotted.Object 30 is the only object positioned among one of the dark clouds in the North of the region.It is likely then that this object is strongly reddened by non-local extinction.

RESULTS AND DISCUSSION OF SPOT PROPERTIES
This section will discuss the results from the peak to peak amplitude measurements and the spot property determination for our targets.Three of the 32 objects with periods determined from  and -band data have no fitted Â { } .Object 19 is removed to due to unreliable  band photometry, as mentioned in Sect. 4. The other two (12 and 15) are due to an insufficient number (less than 50) of data points in .

Object Result Summary Plots
All measured and determined properties for our objects are summarised in result summary plots, the first of which, for object 1, is shown in the top left panel in Fig. 4. The plots for all other objects can either be found in the main body of the paper (1,2,3, and 4) if they are discussed in detail, or in the Appendix A. In this section we describe the details of these plots.
The data points in each panel of these figures are plotted at the Modified Julian Date (MJD) of the start of each slice (bottom x-axis), and the slice number is shown at the top x-axis.Above that axis we also note the period for the object.The vertical gray dashed lines indicate one year intervals every start of February.Thus, the slices just left of the gray lines are centred on the times of worst observability of the field and hence usually contain the least amount of data (as discussed in Sect.3.1).
The first (top) panel in these result summary plots (as e.g.Fig. 4) shows the peak to peak amplitude Â  (green), Â  (red), and Â  (black) values as coloured dots and their uncertainties  Â  .The second panel shows the phase position (in degrees) of the maximum brightness in the phase folded light curve in each filter with the same colour coding as in the top panel.The zero point for these phase values has been set as the median of all determined values for each object.The blue line indicates the running median of the phase positions over all filters in the slice and the light blue shaded area indicates the RMS scatter of the points from that running median.
As laid out in Sect.3.3, we require a SNR above three for all amplitude measurements and an RMS scatter below 45 degrees for the minimum/maximum brightness phase positions.Thus, only data in slices where these criteria are satisfied are shown as coloured dots in the top two panels.In slices where one or more of the filters do not have a SNR above three for the peak to peak amplitudes, the 3 uncertainties are shown as light-coloured arrows instead and no phase positions are shown.The first example of this is object 3 (see bottom left panel in Fig. 4).In cases where the amplitudes have the required SNR but the phase position varies too much, the three sigma noise is displayed as light coloured crosses and again no phase positions are shown.This happens e.g. for object 7 (see Fig. A3 in the Appendix).In all slices without any data shown, there are less than 50 data points in at least one of the filters.
The third and fourth (bottom) panels in the result summary plots show the determined spot properties and uncertainties for each slice, if the Â{ } amplitudes and phase positions have the required accuracy.The spot temperatures (in degrees Kelvin) are shown in the third panel.The dashed horizontal line indicates the surface temperature of the star (as listed in Table 2).The bottom panel shows the spot coverage.The displayed error bars represent the MAD uncertainties, as explained in Sect.3.4.The colour coding of the symbols in the bottom two panels follows the scale shown at the right hand side of the figure.It represents the ratio of the hot spot to cold spot solutions (HS:CS { } , see Sect.3.4) during the spot property fitting.Thus, spots that are mostly fitted as cold spots are shown in red, and spots that are mostly classified as hot spots are displayed in blue.

Amplitude and Phase evolution
A quarter of our objects have Â { } values near continuously from the start to the end of their light curve, only missing individual or pairs of Â { } due to insufficient data or low signal to noise.Two objects only contain a single slice with Â { } measurements.The median number of Â { } values is nine, and the maximum is 15.The primary reason for slices missing Â { } values is insufficient data.As mentioned in Sect.3.1, slices centred on February contain the fewest data points, and as such these slices are often missing.A number of objects show longer breaks in Â { } values due to low signal to noise values and have hence the 3  noise shown in their result summary plots.These sources are thus showing longer periods of no or undetectable spots on their surface.We are unable to observe any activity cycles in our sample as this requires longer light curves.
Note that we are detecting the periodic variability due to rotational modulation over the ∼ 180 d length of a light curve slice.

Not showing Â𝑜
{ } values with our required signal to noise during a slice does not mean that the object has become quiescent.Rotational modulation is caused by asymmetry in the spot distribution.Thus, two similar surface features separated by 180 • in longitude, i.e. on opposite 'faces' of the star, will not produce a detectable rotational brightness modulation.Additionally, a significant change of the spot properties (temperature, coverage) or position (phase) during the 180 d in a slice, will not allow us to detect the periodic brightness changes.Such slices will hence usually only have the noise shown in the result summary plots.
Examining the phase position over the length of the light curves provides a further check of the quality of our period determination.In case there is a systematic, linear trend in the phase position over time, then this can be considered as an indicator of a small systematic error in the period of the object.None of our objects shows such a behaviour over the typically 4 -5 yr length of the light curves.Furthermore, the phase position tracks the longitudinal position of the surface feature on the star.Thus, changes in phase can track

Spot Property Distribution
Our previous work (Herbert et al. 2023) investigated spots on YSOs in IC 5070 during an 80 d period based on periodic YSOs identified in Froebrich et al. (2021).We identified 21 stars that had cold spots and six that had warm/hot spots.Thus, the cold spot objects outnumbered the hot spot objects by a factor of 3.5.For about 13 percent (4 stars) of the periodic objects we could either not find a valuable spot solution (2 stars) or they are potential AA Tau-like contaminants (2 stars).
We now consider the determined spot properties from the observed Â { } values for all slices in our periodic objects.Each of the 234 slices with high signal to noise peak to peak amplitudes is considered in the statistical analysis.We note, however, that neighbouring slices do overlap by 3 months.In those 234 amplitude sets we can reliably fit 180 as spots.Of these we identify 125 as cold spot solutions and 55 solutions are warm or hot spots.Thus, cold spot solutions outnumber the warm/hot spot solutions by a factor of about 2.3, slightly smaller than in Herbert et al. (2023).In the remaining 54 slices no reliable spot solution can be found (21) or the best solution is at the edge of the investigated parameter space (33).This fraction of about 23 percent is about twice as high as in Herbert et al. (2023).
The difference between the spot temperature and stellar surface temperature   −  ★ for all our solutions is plotted against the spot coverage in the left hand panel of Fig. 5.The colour coding of the symbols refers to the HS:CS { } ratio.We see that there is a 'gap' in the temperature difference distribution between about two and three thousand Kelvin.We hence refer to spots with a temperature difference above 0 K and below 2500 K as warm spots, and to all spots that have more than 2500 K temperature difference as hot spots.The differences between these groups will be discussed in Sect.5.3.2.In the right hand panel of Fig. 5 we show the distribution of the spot coverage for hot, warm, and cold spots as histograms, as well as cumulative distribution functions (CDFs) for the latter two.

Cold spots
In Herbert et al. (2023) we simulated periodic variability due to extinction from interstellar dust, and found that if we try to fit this as spots, these spots clustered at the lower temperature boundary of the PHOENIX models (2000 K), or at the  = 0.5 boundary for the coverage.We therefore remove spots that are fitted as implausibly large or cold, i.e. within 10 percent of  = 0.5 and   = 2000 K. 27 spots were removed because of their temperature, all but two were within 2.5 percent of the lower temperature limit.Six further spots were removed for being too close to  = 0.5.Although we adopted a ten percent margin near this value, all spots that were removed were within three percent of it.The largest spot remaining in our sample has a coverage of  = 0.408.The clustering around the parameter space boundaries indicates that the peak to peak amplitudes for these objects cannot be fit with our simple spot model.The remaining objects which are considered valid solutions for the spot properties are significantly away from the parameter space borders.
Among the cold spots 90 percent are between 0 K and 1500 K below the stellar surface temperature.The coldest spots with temperatures more than 1500 K below the surface temperature, have a spot coverage of less than  = 0.25.We do not observe cold spots with a coverage above  = 0.25 and more than 1500 K, below the surface temperature, which has several contributing factors.These objects likely have a complex network of surface features, and the model simplifies it to one uniform spot temperature and a fixed surface coverage.The coverage we calculate is the asymmetrical part of the spot distribution, and so it represents a lower limit.The distribution of the cold spot coverage values appears homogeneous from the lower limit of  = 0.06 to  ∼ 0.3, where it tails off.A two sided KS-test indicates a probability of 34.4 percent that the observed coverage distribution in this range is homogeneous.This is well above the 5 percent threshold for rejection of the null hypothesis although not conclusively positive with limited statistics.
Spots are modelled up to 0 K temperature difference compared to the star and zero spot coverage.However, there is a clear observational bias just above and below   −  ★ = 0 K and for very small coverage values.There the temperature difference and/or spot coverage are not sufficient to produce a SNR above three for the peak to peak amplitudes.This is demonstrated in the left panel of Fig. 5, where for the cold spots the minimum coverage is  = 0.061 and the smallest temperature difference is 283 K below the stellar temperature.Thus, we are unable to detect objects above the upper envelope of points in the   − ★ vs.  plot (Fig. 5, left).The lack of objects at larger temperature differences and large coverage values is however real.We note that our simulations in Herbert et al. (2023) have shown that there are small systematic effects in our methodology.These cause the spot properties near the upper envelope to be determined as slightly smaller and at a slightly larger temperature difference.However, these systematic shifts are always smaller than the MAD uncertainties for the spot temperature and coverage.

Warm and hot spots
For spot solutions with temperatures above the surface temperature, the spot parameters cover a wide range of temperatures and sizes.90 percent of these spot solutions have temperatures that are less than ∼2000 K above the stellar temperature.We also find a dearth of spots between 2000 K and 3000 K above the stellar temperature.Due to this large apparent gap in the distribution, we use a spot temperature of 2500 K above the stellar temperature to separate warm and hot spots, and discuss these separately.The number of warm/hot spots in the earlier sample from Herbert et al. ( 2023) was very small.Only six warm/hot spot solutions were found.Five of these were warm spots and one object had a spot temperature of 2650 ± 970 K, i.e was situated in our apparent gap of temperature differences.bf In the future we intend to investigate the YSOs in all HOYS fields and will thus have a much better statistical picture if this apparent gap is indeed real.
We identify five hot spot solutions that have a temperature difference of more than 3000 K above the stellar temperature.The upper limit for the temperature in the spot models is 10000 K.There is no single case of any spot solutions approaching this limit.This is most likely due to the simplicity of our model (it averages the spot temperature over the coverage) and the unavailability of any shorter wavelength (B and U-band) peak to peak amplitudes.All of the hot spots are very small, with a maximum coverage of  = 0.027.This is expected for accretion column footprints (e.g.Hartmann et al. 2016;Muzerolle et al. 1998).Four of the five hot spot solutions occur on a single object (30, see Fig. A21).We will discuss this star in more detail in Sect.6.7.We currently have very limited statistics about the hot spot solutions.
Similar to the cold spots, there is a bias against the detection of small warm spots with low temperature differences compared to the stellar temperature.The smallest temperature difference measured in the warm spot regime is 310 K above the stellar temperature.However, due to the temperatures and filters involved, the regions of undetectable warm spots (defined by the lower envelope in the   −  ★ vs.  plot in the left panel of Fig. 5) are much smaller.Similarly we find again that higher temperature differences between spot and stellar surface usually correspond to lower a coverage.However, small a coverage does not necessarily correspond to a high temperature difference, as a number of warm spots also occupy the low temperature difference and coverage range.Furthermore, there is a much better defined upper envelope in the   −  ★ vs.  plot, above which objects do not seem to exist.
Compared to the hot spot solutions, the warm spot solutions cover a wide range of sizes.The largest warm spot is as large as the largest cold spot.However, unlike the cold spots the distribution in coverage is not homogeneous.The histogram in the right hand panel in Fig. 5 shows a clear peak at lower coverage values.We have performed a two-sided KS-test of the coverage values for the warm spots with coverage values between zero and 0.2.We find a probability of 43.3 percent that the distribution follows a normal distribution in this range.However, we note that there is a clear tail of much larger coverage values and that the low coverage values are influenced by our detection bias.There are eight warm spots (16 percent of total warm spots) that have a coverage above  = 0.2, as opposed to a 60 percent of the cold spots.The largest warm spot in our previous sample from Herbert et al. (2023) was  = 0.125.

Rapidly changing spot properties
If one looks through the result summary plots for all objects (e.g.Fig. 4 or in Appendix A), one finds several examples where a spot solution changes from a warm to a cold spot (or the other way) in neighbouring slices.Usually these cases create single-slice warm spots on otherwise cold spot dominated objects.Note that there is not a single case of three consecutive slices where the middle one has a different classification (warm/cold) from the two others.Here we briefly discuss these changes in the classification from one slice to the next, and how they relate to changes in the HS:CS { } ratio.
In our spot fitting methodology (see Sect. 3.4 and details in Herbert et al. 2023) we vary the measured Â { } values within their uncertainties 10000 times, to estimate the MAD errors for the spot properties, and to determine the HS:CS { } ratio.This ratio indicates the fraction of all solutions that are hot/warm spots compared cold spots.We have introduced a somewhat arbitrary cut-off for these HS:CS { } values between 0.4 and 0.6 for which we consider the spot fitting solutions ambiguous.Examining the relationship of the spot parameters (in particular the temperature) and HS:CS { } over time, we find that in many occasions towards a switch from a cold to a warm spot or visa versa, the HS:CS { } value does approach the boundary for ambiguous solutions.
In our rigorous testing of the spot fitting methodology in Herbert et al. ( 2023) we have shown that a small hot spot can make a significant contribution to the peak to peak amplitudes.In the 180 d long slices used here, the asymmetric part of the spot distribution on the stars can evolve.In particular contributions from accretion hot spots, caused by material flowing along magnetospheric structures which might be short lived compared to the duration of a slice (e.g.Donati et al. 2012), need to be stable throughout the slice to create regular periodic light curves (Kurosawa & Romanova 2014), and hence allowing us to reliably determine their resulting spot properties.Thus, evolving spots will most likely lead to HS:CS { } values changing or moving into the range where no clear solution can be found.Given the very small fraction (5/234, i.e. 2.1 percent) of hot spot solutions, we can conclude that it is highly unlikely that the footprints of the accretion columns are stable over timescales of half a year.We discuss object 30, which seems to be the exception to this in Sect.6.3.On the other hand, the much larger fractions of reliable warm and cold spot solutions found indicate that the physical causes (see Sect. 5.4) for these are typically stable for longer than half a year.
In Fig. 6 we show two examples (object 13 left panel, and object 18 right panel) of these changes for the classification from a warm to a cold spot solution within the duration of one slice.Object 13 shows predominately cold spots, with exceptions in slices 21 and 24.The slices next to them, which contain overlapping data show cold spot solutions.The peak to peak amplitudes surrounding this pattern do not show any corresponding significant change.For slice 21, the HS:CS { } ratio is 0.608, very close to our cut-off.The spot properties for slice 23 are ambiguous since the HS:CS { } ratio is 0.423, and in slice 24 the ratio is 0.795.Thus, during this entire period the surface features on this objects probably evolve on timescales shorter than half a year and thus no clear, indisputable spot solution with our simple model can be found.Similarly, slice 20 in object 18 shows a brief change in classification (warm spot compared to cold spots throughout), accompanied by a HS:CS { } ratio of 0.848.Note that in this slice the peak to peak amplitudes are clearly different, especially in the -band, and that there is no data in the preceding slice 19.Thus, the statistical interpretation of the distribution of the warm and cold spot solutions, needs to take the respective HS:CS { } ratio into consideration.In other words, the more lightly coloured spot solutions in the left panel of Fig. 5 should be considered slightly less reliable than the darker coloured ones.

Relation of Spot and Stellar Properties
In this section we briefly examine how the determined spot properties relate to the stellar properties.We reiterate that these objects are likely to exhibit multiple, potentially evolving, surface features at once.As discussed in Sect.5.2, rotational modulation is caused by the asymmetric part of the spot distribution.In order to represent the spot properties of individual objects, we will refer to dominant spot properties throughout this section.These are defined as the proportion of reliable spot property solutions for an object, that are above (warm/hot) or below (cold) the stellar temperature.

Stellar Rotation Periods and Disc Excess Emission
In Fig. 7 we show how these dominant spot properties relate to the period and  − 2 colours for the objects.The marker colour in the figure indicates the proportion of the spot properties that are warm/hot spots.One (dark blue symbols) indicates that all reliable spot solutions are warm/hot and zero (dark red) means they are all cold.This is not to be confused with the ratio HS:CS { } , which is part of the result of the spot fitting process for an individual slice.
There are three objects that show equal numbers of warm/hot and cold spots.The marker size in the plot is scaled to the number of spot measurements, i.e. larger circles represent sources with more reliable spot solutions.Note that this is not equal to the number of Â { } values, as some of them do result in ambiguous spot solutions for which no spot temperature and coverage can be reliably determined.
In Fig. 7 we show a horizontal dashed line at  −2 = 0.5 mag (Teixeira et al. 2012) which separates sources with a detectable warm inner disc (above) and without (below).The 2 data is obtained in all cases from the WISE All-sky catalogue (Cutri & et al. 2013).For the K-band magnitude we generally use 2MASS (Skrutskie et al. 2006), with the exception of object 17 for which this is not available.We utilise the UKIDSS GPS data (Lucas et al. 2008) for this source.The absence of excess  − 2 emission typically indicates a lack of warm dust in the inner (∼1 AU) disc.The vertical dashed line at 5.5 d in Fig. 7 separates fast (left) from slow (right) rotators.The exact value for this has been adopted for consistency with Froebrich et al. (2021), who found a clear gap in the period distribution in IC 5070 at this value.
We find that 15 objects of our sample do not show a disc excess emission and 14 objects do.This disc fraction of about 50 percent is consistent with the value for this cluster determined in Froebrich et al. (2024), and also with the typical age of about 1 Myr of the objects in the region (Kuhn et al. 2020).Based on this age it is unlikely that the non-disc objects have evolved debris discs.It is more likely that they have cleared their inner discs due to the start of planet formation, and in later stages, photo-evaporation.The objects that have cleared their inner disc are known as transition disc sources.An alternative for the non-disc objects is that they retain the inner disc edge, but the disc contains a cavity in the inner part.These are referred to as pre-transition discs.
Observations at longer wavelengths can be used to probe the presence of disc material further from the star.Unfortunately almost half of our objects do not have a detection with a sufficient SNR in the W3 and W4 bands to characterise the cooler material of the outer disc.There are 17 objects, however, that have a better than 5 detection in all four WISE filters.We determine the slope of the spectral energy distribution (  ) following (Majaess 2013) for all of these.We find that amongst the 17 objects, one is classified as a Class 1 protostar (28), 15 are Class 2 objects, and one (7) is a Class 3 source.Thus, if these objects are representative of the entire sample, we find that 88 percent of our objects are Class 2 sources, in good agreement with the typical age of the stars in the region.
The bimodal period distribution in our sample has been discussed in Sect. 4. A longer period may indicate that the central star still has magnetic ties to the inner disc, slowing rotation.Faster rotators have been spinning up due to a lack of this disc-breaking mechanism (Herbst et al. 2007).The majority (20/29, i.e. 70 percent) of our objects are fast rotators.Of those, a majority of 13/20 (65 percent) show no detectable inner disc.There is a loose positive correlation of the rotation period and disc excess emission, with the exception of a few select objects (14 and 17), which have a longer period without inner disc emission.The very fast rotators ( < 3 d) do not show evidence of an inner disc.

Causes of the Surface Spots
In Fig. 7 we can see that the majority of our objects are dominated by cold spots.This is expected as 70 percent of our reliable spot solutions are cold.The presence of these cold spots on young stars is well understood as being caused by magnetic activity on the stellar surfaces, driven by the rotation of the stars (see reviews by Herbst et al. (2007) and Bouvier et al. (2014)).These spots are known to have a range of temperatures and can reach up to half the visible surface (Strassmeier 1992;Bouvier et al. 1995).
Our sample also includes eight warm/hot spot dominated objects.One of these consists of only a single warm spot measurement.We also find that objects with a higher number of reliable spot solu-tions (larger circles) are more likely to be dominated by cold spots.Investigating the result summary plots for all sources in our sample shows that cold spot dominated objects are more 'consistent' in their behaviour.In other words they are more likely to show cold spot solutions throughout, or for longer continuous periods.On the other hand, warm/hot spot solutions on objects tend to be shorter lived.Only four individual objects in the sample contain four or more warm/hot spot solutions.Thus, a more reliable statistical analysis of these warm/hot spot dominated sources will require a larger sample.
All but one of the warm/hot spot dominated objects have  − 2 colour below 0.5 mag, indicating a lack of an inner disc.Can these mostly warm/hot spot sources be caused by accretion, despite the lack of a detectable inner disc?Accretion on such transitional disc objects has been the source of multiple studies.In similarly aged regions the proportion of accreting transition discs is 70 -80 percent (Kim et al. 2009;Cieza et al. 2010;Merín et al. 2010).Sicilia-Aguilar et al. (2010) found that when transition disc objects were accreting the median accretion rate was comparable to full discs.Similarly, Manara et al. (2014) studied transition disc objects in a similar age range to IC 5070 and found that among the transition disc objects that are still accreting, accretion occurs at a similar rates to full discs.The discs tended to be gas-rich close to the inner edge while being dust depleted.Although, Najita et al. (2015) found that when compared with stellar and disc mass, transition disc objects have lower accretion rates.Furthermore, eleven of our objects are identified in Panwar et al. (2023) as H emission sources.Not all our sources are in the surveyed region.The slit-less spectroscopy data for that study was obtained in December 2012 and January 2013, and so unfortunately are not concurrent and the accretion rate may have varied since 2012/2013.
The above reinforces that many of our hot/warm spot dominated objects are still accreting despite their lack of detectable inner discs.Can the determined hot/warm spot properties be explained by accretion?Hot spots as the footprints of accretion are predicted to have temperatures of 8000-10000 K and a small surface coverage of less than a few percent (Hartmann et al. 2016;Muzerolle et al. 1998).
Recent works have identified a wider range of lower temperatures and higher coverage values (Scholz et al. 2009(Scholz et al. , 2012;;Bozhinova et al. 2016).These lower limits of temperature do, however, not explain the lines identified with spectroscopy which require higher energies (Muzerolle et al. 2003(Muzerolle et al. , 2005)).In Espaillat et al. (2021) the dense line producing region has a coverage of 0.1 percent, but the cooler, extended region has a coverage in the range of 10-20 percent of the stellar surface.When the high density region disappears, the UV emission dramatically decreases but the optical peaks remain.Thus, a large fraction of our warm spot solutions could be explained by low density accretion columns.
Furthermore, Sicilia-Aguilar et al. ( 2023) identified stable accretion hot spots on EX Lup and TW Hya with emission line spectroscopy.The photometry of the objects was dominated by spots that are 500 -2000 K above the stellar surface temperatures.These spots were found to vary in longitudinal position distinctly from the line emission region.In the case of EX Lup, which has a cleaner sinusoidal rotational modulation than TW Hya, the warm spot was sometimes ahead and sometimes behind of the line emission hot spot.Over the course of 55 days the spot reduced in size from 88 percent to 10 percent coverage while simultaneously increasing in temperature from 4700 K to 5200 K.This behaviour is consistent with the rapidly changing spot properties discussed in Sect.5.3.3.
Accretion is not the only possibility for warm spots on young stars.Froebrich et al. (2020) identified warm spots on the active young star V5198 Cyg, an object which shows no evidence of ac-cretion or mass transfer.Chromospheric active regions known as plages or faculae have been identified on stars of varying ages and types, from active main-sequence stars (Strassmeier et al. 1993), to late type stars (Zhang et al. 2014) and commonly on close binaries (Frasca et al. 2008).However, there is limited data or studies of these on young stars.The young solar type star HD 171488, aged 30 -50 Myr, has been observed to have chromospheric plages and spots with a 20 • difference in phase (Frasca et al. 2010).
Several studies have also identified photospheric warm spots on WTTSs (Xiang et al. 2023;Gully-Santiago et al. 2017;Donati et al. 2014).These are also referred to as plages but are distinct from chromospheric plages.Strassmeier et al. (2019) observed warm spots alongside cool spots with opposing polarities, each with a coverage of a few percent, spatially distinct with longitudinal separations of 10 • to 50 • .These are presumed to be caused by mass flow between the regions of opposing polarities.The warm regions in this case would always be accompanied by cool regions.In our methodology, a smaller warm region will dominate the photometric variability in the light curves.This scenario hence may provide an explanation of the behaviour seen in many of our objects (see Sect. 5.3.3),where spot temperatures change (from warm to cold and vice versa) without detectable phase changes.
Photospheric plages have been observed to have coverage values on the smaller side of our warm spot distribution.Donati et al. (2017) observes a coverage of six percent for warm spots but also clarifies that Zeeman Doppler Imaging constrains the lower boundary for coverage.A number of our low coverage warm spots may hence also be explained by photospheric plages.Low density accretion columns as in Espaillat et al. (2021) have an upper limit of 20 percent coverage.However, the very large warm features observed on EX Lup in Sicilia-Aguilar et al. ( 2023) reveal a possible explanation for the small number of larger (  > 0.2) warm spots in our distribution.Object 28, which accounts for five of these eight large warm spots, is discussed in Sect.6.6.

The variable star V1701 Cyg
The variable star V1701 Cyg was identified in Froebrich et al. (2021) (object ID: 4766) with a period of 11.7162 d.This period was not found again in our work.This source is the only object in the entire sample for which this occurred.Slice 17, which contains the majority of the 2018 high cadence data used in Froebrich et al. (2021), resulted in a candidate period of 6.55 d.We found five distinct candidate periods for this object in different slices.The most common being 6.55 d, which was identified in two different slices.The additional four candidate periods were 7. 67 d, 9.69 d, 12.61 d, and 14.58 d.The spot fitting in Herbert et al. (2023) identified V1701 Cyg as a likely AA-Tau type contaminant.A  − 2 colour of 1.6 mag showed the presence of an inner disc, and the slope of the spectral energy distribution (  = −0.62)from WISE classified it as a Class II object.It is therefore possible that the lack of consistent periods indicates dust occultations at varying distances from the central star over time.Hence, while not suited to our analysis here, this object is interesting and should be followed up in detail.

Very Cold Spots or Dust Extinction on Object 2
Object 2 (see result summary in top right panel of Fig. 4) contains ten Â { } measurements.Yet only four reliable spot solutions are found.The remaining six slices were removed as solutions were within 10 K of 2000 K. Four of these slices are 'pairs' with overlapping data.The causes of these very cold spots are hence either dust extinction from the disc, or cold spots below the 2000 K lower limit.
This object was discussed in Sect.4.3, as the brightest object in the sample at G ≈ 12 mag.The stellar temperature of 5000 K was chosen to account for this, and yet the majority of the spot solutions are at a temperature difference of 3000 K below stellar surface.We have considered whether an underestimation of the stellar temperature would lead to the spot solutions changing.In Herbert et al. (2023) we tested the effect of altering stellar temperature by ±200 K, but this altered the spot temperatures by less than the statistical uncertainties.The hottest object in the sample is object 30 (see Sect. 6.7) which has a temperature of 5500 K.It cannot be directly compared with object 2 as it undergoes extreme reddening due to non local extinction.While near a bright rimmed globule (see Fig. 3) it is not located within an obvious dark cloud.Even if the true temperature is 5500 K, then this would still not effect the spot temperatures by more than the typical statistical uncertainties.Object 2 would therefore require a significantly higher temperature in order to increase the temperature of the resultant spot properties.
The object has  − 2 = 0.20 mag and   = −1.39.This classifies it as a Class 2 object with a transitional disc.For dust to produce AA-Tau like variability it has to be close to the co-rotation radius, and the  − 2 value implies that this is not the case.Furthermore, we have shown in Herbert et al. (2023) that dust extinction causing spot solutions to approach the lower temperature limit are caused by achromatic extinction.The Â { } values in object 2, however, are not achromatic.The highest filter with the highest amplitude varies throughout the light curve.In general, Â is expected to be the highest and Â  the lowest.The light curve modulation is strong, with an amplitude (0.2 mag) that is higher than the median of our sample, but the modulation in colour is erratic.
The alternative solution is that the spots can indeed be colder than 2000 K at times.In either case, the source warrants more detailed further investigation.

Phase shifts on Object 3
Object 3, shown in the bottom left panel in Fig. 4 has ∼2.5 yr of detectable rotational variability at the start of the light curve.During that time the phase position remains constant within the measurement uncertainties.Note that the dominant spot temperature changes between warm and cold during that period.This is discussed in Sect.5.3.3.After slice 26 there is break in detectable Â { } values for about 18 months.Two slices in the break have sufficient data but the Â { } values have a SNR below three (due to high noise in the  band).When the Â { } values resume, the amplitudes are reduced and there is a 90 • shift in phase.The higher noise in the gap may be due to more complex spot activity.When the sinusoidal modulation resumes after the gap, a surface spot at a different part of the surface compared to the earlier dominating spot might have formed.

Potential Spot Motion on Object 4
The evolution of object 4 is shown in the bottom right panel in Fig. 4.This source only has four gaps in the Â { } values, mostly in the slices centred on the worst observability of the source.Between slices 17 to 23 the amplitudes increase and decrease slowly, while showing a near constant phase position.This changes from slice 25 onward.A warm spot is found with a phase shift of 90 • compared to the earlier cold spots.Together with the lack of successful spot property fitting for slices 22 and 23, this indicates a more complex situation of no clear single dominating spot for the duration of a slice (half a year).Between slices 25 and 34 the Â { } values increase and then decreases (similar to the first two years).During this time the phase position changes gradually.Again, in some slices no spot properties can be fitted, indicating a constant change in the nature of the dominating spot.However, the longitudinal position of the most asymmetric part of the spot distribution on the surface of the star is systematically moving at a rate of about 90 • per year.

Very Cold Spots on Object 5
Similar to object 2 (Sect.6.2) object 5 has eight spot temperature solutions too close to the lower 2000 K limit.Indeed these two objects account for nearly half of these cases.Object 5 has an effective temperature of 3686 K (Fang et al. 2020),  − 2 = 0.79 mag, and   = −1.02.Thus, it is a Class 2 object with an inner disc.The spot temperatures (albeit low) and coverage values for the reliable solutions do not stand out compared to the values found in the rest of the sample.Thus, similar to object 2, this source might have some of the lowest temperature spots and should be investigated in more detail.

Large Warm Spots on Object 28
Object 28 displays five of the eight warm spot measurements with coverage values over  = 0.2.It has the longest period of the sample at  = 13.69 d and  − 2 = 1.20 mag.It hence has an inner disc, and it is therefore likely that this object is accreting.However, it is not identified as an emission line star in Panwar et al. (2023), but this does not exclude accretion as the data is not contemporary to our observations.It also has   = 0.32 mag, making it the one Class 1 protostar in our sample.
In Herbert et al. (2023) we tested the outcome of simulated variability induced by dust with   = 5 and   = 3.1.This could in principle result in large warm spots.However, this can be excluded for this source as the colour modulation is not in agreement with such extinction.Object 28 has above average peak to peak amplitudes in our sample, with the largest amplitudes in the  band.Coverage values above  = 0.2 are not expected from low density accretion columns (Espaillat et al. 2021).Large warm regions on EX Lup (Sicilia-Aguilar et al. 2023) show that on a stably accreting object it is possible for the photometry to be dominated by large features with low temperature contrast.
Object 28 is one of the two objects in our sample with three consecutive warm spot solutions in slices 21 -23.Slice 24 has insufficient data, and then slices 25 and 26 have two further large warm spot solutions.This period of extreme stability implies large warm features as observed on EX Lup.
The stability seen on object 28 could partly be due to its early evolutionary stage.As pre-main sequence stars develop a radiative core their magnetic fields become increasingly complex, and their dipole field becomes weaker (Gregory et al. 2012;Donati et al. 2011).This object is unlikely to have had sufficient time to develop a radiative core, therefore it is likely to host a strong, simple magnetic field.This supports the stability of the spot observed here.

Stable hot spots on Object 30
Object 30 is discussed in Sect.4, as the hottest star in our sample and with high reddening due to non-local extinction.Both, the period and  − 2 indicate that this object has no inner disc.Object 30 accounts for four of the five hot spot solutions in our sample.These spot solutions have a temperature of more than 3000 K above the stellar temperature.The object has a short period of 1.43272 d and  − 2 = 0.234 mag.It is not included in the emission line star catalogue of Panwar et al. (2023) as it is outside of their search area.
The spot solutions for the object start with two cool spots in slices 17 and 18.This correspond to the high cadence period of the light curve in 2018.The HS:CS { } are 0.27 and 0.34 indicating these solutions are credible but not definitive on the side of cold spots.From slice 21 to the end of the light curve, every Â { } identified leads to a strong warm/hot spot classification, with the lowest HS:CS { } being 0.75.Four of these spot solutions have temperatures above 8900 K and very small coverage values.These align with the expectations for accretion based hot spots.
There are two 'pairs' of hot spots in slices 21/22 and 32/33, with overlapping data.There are two spots on this object which are not 2500 K hotter than the stellar surface and are hence classified as warm spots.These have temperatures of 6700 K and 7300 K, which is not out of the range of accretion based activity.Hence, we can reasonably consider this object to be accreting, and all of the warm/hot spot solutions can be considered accretion based hot spots.Froebrich et al. (2024) identified 366 YSOs in IC 5070 from Gaia DR3 data.Of these, 131 have light curves with at least 100 data points in , , and  from the HOYS citizen science project.We have sliced these light curves into six month long sections, every three months.In this way, half of each slice overlaps with the previous.We searched for a period within each individual slice that contained at least 50 data points in two of the three filters.

CONCLUSIONS
Following the method established in Herbert et al. (2023) we have used the peak to peak amplitudes in , , and  to measure spot properties.Out of the 234 amplitude sets, we find 180 reliable spot solutions.Amplitude sets that are not fit as spots are either near the edge of our parameter space, or create ambiguous solutions between cold and warm spots.The latter is most likely caused by complex spot structures, such as multiple spots with different temperature contrasts, which cannot be fit by our simple model.We are not able to observe any cyclical behaviour in our objects, as this requires longer light curves than currently available.
We measured a wide range of spot properties.Roughly two thirds of all reliable spot solutions are cold.We find the coverage values of cold spots are homogeneously distributed between our lower detection limit of 6 percent and 30 percent, with maximum coverage values of around 40 percent.90 percent of the cold spot solutions are between 0 K and 1500 K below the surface temperature.The minimum temperature difference is ∼ 280 K and the maximum temperature difference is ∼ 2150 K.We do not consider spot solutions that are within ten percent of our lower temperature limit (2000 K), as these measurements are potentially non-spot contamination.Alternatively, the spot temperatures are indeed below the lower temperature limit of our models.We discuss two objects (2, 5) that display the majority of these spot solutions and conclude that in these cases it is likely that the spot temperatures are indeed lower than we are able to fit.
We find that a significant proportion (almost one third) of our spot solutions have temperatures above the stellar temperature, but by less than ∼2000 K.We find no spot temperatures between ∼2000 K and 3000 K above the stellar temperature.Thus, we suggest a limit at 2500 K above the stellar temperature, separating warm and hot spots.Warm spots have a temperature contrast between 0 K and 2500 K and hot spots are more than 2500 K above the stellar temperature.
Warm spots have a wide range of surface coverage values.The maximum coverage is around 40 percent, although 84 percent of warm spots have a coverage of less than 20 percent.The coolest warm spots have a temperature contrast of ∼ 310 K above the stellar temperature.Low temperature contrast, large coverage spots have been observed in other objects (Scholz et al. 2009(Scholz et al. , 2012;;Bozhinova et al. 2016).Accretion hot spots are expected to have temperatures of 8000 -10000 K and coverage values less than a few percent (Hartmann et al. 2016).The line emission associated with accretion requires these high temperatures.However, recent works by Espaillat et al. (2021) and Sicilia-Aguilar et al. ( 2023) have shown that the photometry of stably accreting objects is not necessarily dominated by the small, hot accretion shock region due to extended, low density accretion columns or spatially distinct warm regions.We show that these warm surface regions are common in our sample., , and  band data does not cover the peak of the spectral energy distribution for high temperature accretion shocks, and therefore our methodology is more sensitive to low temperature contrast warm spots.
We find five hot spot solutions that are ∼3000 K above the stellar temperature, with coverage values of less than three percent.These align with expectations for accretion column footprints.Four of these five hot spot solutions are on the same object (30), which is a YSO without a detectable inner disc (in  − 2).We consider this object to be actively accreting over the gap in the inner disc.The rarity (5/180, ∼3 percent) of such hot spot solutions shows that such accretion column footprints are either i) usually not stable on the timescale (6 months) at which we are able to investigate the spot properties, or ii) the photometry is not dominated by the accretion shock region as is discussed with the warm spot regime.
Our entire YSO sample has an inner disc fraction of about 50 percent, based on the  − 2 infrared excess.The objects in our sample are dominated by cold spot variability, which is expected considering that approximately two thirds of our spot solutions are cold.We find that objects with more reliable spot solutions overall tend to have a higher proportion of cold spot solutions.Warm/hot spot dominated objects tend to have fewer reliable spot solutions.This implies that cold spots are longer lived, or more likely to produce rotational variability over longer timescales.Increasing our sample size will allow us to investigate this relationship.
We find eight warm/hot spot dominated objects.Only one of these has a detectable inner disc (28).Accretion is common over an inner disc gap, although whether the rate is equivalent with full disc object is subject to debate, and so these warm/hot spots may be accretion related features.Alternatively, warm photospheric plages have been observed on WTTSs with coverage values of around six percent (Donati et al. 2017).The low coverage warm spots may be accretion related warm regions or plages, whereas the larger regions are more likely accretion related.The single warm/hot spot dominated object with a detectable inner disc is the only Class 1 protostar in our sample.It displays large warm spots in consecutive slices indicating extreme stability in the spot behaviour and therefore accretion.
In the future, we intend to apply our methodology to YSOs in the remaining HOYS fields, significantly increasing the number of investigated YSOs.This will allow for a much better statistical investigation of in particular the warm/hot spot solutions, the relationship between disc presence and spot properties, and the investigation of spot lifetimes.

Figure 1 .
Figure 1.Distribution of measured periods of YSOs in our IC 5070 sample.The solid black line represents the cumulative distribution function.The dashed vertical line placed at a period of 5.5 d separates slow from fast rotators.

Figure 2 .
Figure 2.  −  vs.  colour magnitude diagram.All selected IC 5070 cluster members as shown as black points and periodic objects are in colour mapped to their effective temperature from Fang et al. (2020), or the sample median.The marker size is scaled to the period, as shown in the legend.Overlaid in blue is a 1 Myr PARSEC isochrone(Bressan et al. 2012).

Figure 3 .
Figure 3. DSS Red image, 1.2 degree radius from the centre point.Centre point marked with red cross.Objects in our sample are labelled and color scaled to the Gmag brightness.

Figure 4 .Figure 5 .
Figure 4. Result summary plots for the first four objects in our sample.The data for object 1 are shown in the top left panel, object 2 in the top right panel, object 3 in the bottom left panel, and object 4 in the bottom right panel.In each of the summary plots we show from top to bottom the peak to peak amplitudes, the maximum brightness phase positions, the determined spot temperatures, and the spot coverage against time.For more details see the description in Sect.5.1.

Figure 7 .
Figure 7. Period against  − 2 for our objects.The horizontal dashed line separates objects with a detectable inner disc (above) from sources without (below).The vertical dashed line separates fast and slow rotators.Markers are coloured based on the proportion of warm/hot spot measurements on them, according to the colour scale on the right.Values of zero (dark red) indicate that all spot solutions are cold, values of one (dark blue) mean all spot solutions are warm/hot.The marker size is scaled to number of spot property measurements, as shown in the legend.