J-PLUS: Detecting and studying extragalactic globular clusters -- the case of NGC 1023

Extragalactic globular clusters (GCs) are key objects for studying the history of galaxies. The arrival of wide-field surveys such as the Javalambre Photometric Local Universe Survey (J-PLUS) offers new possibilities for the study of these systems. We perform the first study of GCs in J-PLUS to recover information about the history of NGC 1023 taking advantage of wide-field images and 12 filters. We develop the semiautomatic pipeline GCFinder that detects GC candidates in J-PLUS images and can also be adapted to similar surveys. We study the stellar population properties of a sub-sample of GC candidates using spectral energy distribution (SED) fitting. We find 523 GC candidates in NGC 1023, of which $\sim$300 are new. We identify subpopulations of GC candidates, where age and metallicity distributions have multiple peaks. By comparing our results with simulations, we report a possible broad age-metallicity relation, evidence that NGC 1023 experienced accretion events in the past. The dominating age peak is at $10^{10}$ yr. We report a correlation between masses and ages that suggests that massive GC candidates are more likely to survive the turbulent history of the host galaxy. Modeling the light of NGC 1023, we find two spiral-like arms and detect a displacement of the galaxy's photometric center with respect to the outer isophotes and center of GC distribution ($\sim$700 pc and $\sim$1600 pc, respectively), which could be the result of ongoing interaction between NGC 1023 and NGC 1023A. By studying the GC system of NGC 1023 with J-PLUS we showcase the power of multi-band surveys for this kind of study and find evidence of a complex accretion history of the host galaxy.


Introduction
Globular Clusters (GCs) are ubiquitous compact stellar systems found in most galaxies with stellar masses M star > 10 6.8 M ⊙ (e.g.Eadie et al. 2022).Some of these objects are among the oldest objects in the universe (Larsen 2001) with typical ages larger than 10 Gyr (Strader et al. 2005;Chies-Santos et al. 2011).The A well-described property of the globular cluster population of a massive galaxy is its optical color bimodality, showing that there are subpopulations of this class of objects in most massive galaxies (Peng et al. 2006).The bimodality in GC colors is believed to occur due to differences in metallicities.However, age effects and a combination of age and metallicities effects might also play an important role (Brodie & Strader 2006;Lee et al. 2018).It might also be important to take into account non-linear effects in the color-metallicity relations brought by the horizontal-branch morphology in the optical-bands (Richtler 2005;Yoon et al. 2006;Cantiello & Blakeslee 2007;Yoon et al. 2011;Chung et al. 2016;Lee et al. 2018;Villaume et al. 2019;Lee et al. 2020;Kim et al. 2021).Spectroscopic studies have shown that the blue subpopulations of GCs are more metal-poor than the red populations (Beasley et al. 2008;Usher et al. 2012).From chemical evolution models of galaxies as well as from observations, it is known that in dwarf irregular galaxies and low mass galaxies GCs tend to be metal-poor and blue (Lotz et al. 2004).The fact that most galaxies tend to have sub-populations of GCs can be explained by a hierarchical formation: to form a massive system, many small systems are merged throughout time.We note that optical/NIR colors of GC candidates do not have such bimodal distribution in most galaxies, except for NGC 3115, which seems bimodal (Brodie et al. 2012,Cantiello et al. 2014) in any color and metallicity studied (see Cantiello & Blakeslee 2007;Blakeslee et al. 2012;Cantiello et al. 2014;Cho et al. 2016).
Considering how much colors, ages, and metallicities of GCs are important to understanding the assembly of galaxies and their evolution, investigating new photometric bands and colors as well as the interaction of new colors with stellar models and libraries can be interesting to build more detailed spectral energy distributions (SEDs) for these systems.The Javalambre Photometric Local Universe Survey 1 (J-PLUS; Cenarro et al. 2019) operates with a set of 5 broad-band filters based on SDSS (York et al. 2000;Strauss et al. 2002) and 7 narrow-band filters that cover the main stellar indices from 370 to 900 nm ([OII], Ca H+K,D4000, Hδ, Mgb, Hα and CaT).This filter set makes it possible to study GCs with novel colors and more detailed SEDs.Other surveys that can also add new colors to the study of extragalactic GCs are J-PAS (Javalambre Physics of the Accelerating Universe Astrophysical Survey; Benitez et al. 2014) and S-PLUS (Southern Photometric Local Universe Survey; Mendes de Oliveira et al. 2019).J-PAS is composed of 56 narrow-band filters in the optical, while S-PLUS has similar properties to J-PLUS, employing a twin filter system.
To explore the J-PLUS filter set to study extragalactic GCs, we use as a test case the galaxy NGC 1023.NGC 1023 is a SB0 galaxy, located at 11.1 Mpc away (Brodie et al. 2014), with an effective radius of 48 ′′ (Dolfi et al. 2021).This galaxy is located at the center of a small group of galaxies (Tully 1980), and it is currently undergoing a minor merger with NGC 1023 A ( Barbon & Capaccioli 1975;Hart et al. 1980;Capaccioli et al. 1986).NGC 1023 is characterized by a complex and extended HI cloud, whose densest clump is associated with the companion galaxy (Sancisi et al. 1984a;Morganti et al. 2006).NGC 1023 is consistent with being composed of a nearly classical bulge and a fast rotating disk, as extracted from its planetary nebulae system (Noordermeer et al. 2008;Cortesi et al. 2011).Its star cluster system has been explored before in the literature, including spectroscopy (e.g.Larsen 2001;Chies-Santos et al. 2013;Forbes et al. 2014), and it presents complex kinematics, characterized by 1 www.j-plus.esrotation in the inner disk-dominated region, and gradually turning into a pressure supported system in the outer halo-dominated part (Cortesi et al. 2016).The photometric studies trace the GC system up to 8 effective radii (Kartha et al. 2014) and the spectroscopic sample is selected from the photometric sample to maximize the construction of the MOS (Multi-Object Spectroscopy) masks, leading to a non-uniform catalog of GCs ages, metallicities and velocities.
NGC 1023 was observed with the Javalambre Auxiliary Survey Telescope (JAST80).The data was obtained with a Director's Discretionary Time (DDT) proposal observed during science verification.The Globular Cluster Luminosity Function (GCLF) peak is at M V = -7.5 (Harris 2001), which corresponds to V ≈ 22.7 at 11.1 Mpc.J-PLUS reaches g = 21.5 with a signalto-noise ratio S /N = 5 (Cenarro et al. 2019), allowing us to study the brightest part of the GCLF at these distances.On the other hand, the observations of NGC 1023 used in this article were not done with the standard exposure times of J-PLUS and we are able to also detect faint objects.In particular, we detected about 50% of expected objects with magnitudes between the peak of GCLF and 1 sigma below it and no objects between 2 and 3 sigmas below the peak of GCLF.As a result, our sample represents the majority of GC candidates expected for NGC 1023.
In this work, we propose a methodology to detect and select GC candidates from images obtained with JAST80, aiming at exploiting data from J-PLUS, which can be easily adapted for other photometric surveys such as J-PAS and S-PLUS (see also Buzzo et al. 2022 andChies-Santos et al. 2022).With a catalog of GCs, we are set to investigate the stellar population content of GC candidates in NGC 1023 to create an unbiased magnitude limited catalog.This fact makes this galaxy an excellent case of study as we can explore new methodologies and colors, and compare such results to those in literature.
In Section 2 we describe the data used in this article.In Section 3 we present our methodology to detect and select GC candidates, as well as our pipeline GCFinder.In Section 4 we show our results.In Section 5 we discuss our findings.Finally, in Section 6 we present our conclusions.

Data
NGC 1023 was observed in July 2017 through the DDT proposal 1600101 (P.I.Ana Chies Santos) using JAST80 and T80Cam, a the panoramic camera of 9.2k × 9.2k pixels that provides a 2 deg2 field of view (FoV) with a pixel scale of 0.55 arcsec pix −1 (Marín-Franch et al. 2015).
These data are not part of J-PLUS, but they were observed as part of the commissioning period to test the survey capabilities for extragalactic GC science.This galaxy was observed using all filters available at the JAST80 telescope, namely broad-bands u, g, r, i, z and narrow-bands J0378, J0395, J0410, J0430, J0515, J0660, J0861.However, due to problems related to the calibration of the image referring to the J0395 band, we do not use this band in our work.The data are publicly available 2 .The full width at half maximum (FWHM) is presented in Table 1, as well as the exposure times.A color image of NGC 1023 is presented in Figure 1, built using the software Trilogy (Coe et al. 2012).
To perform this study we crop the original wide-field images to a smaller region of approximately 2 000 X 1 500 pixels (≈ 0.1 deg 2 ) around NGC 1023 using IRAF (Tody 1993) tasks.
The images have been reduced with the standard pipeline developed by the Data Processing and Archiving Unit (Unidad de  (Coe et al. 2012).In order to build the color image, the combination of filters (u, J0378, J0410), (J0430, J0515, g) and (J0660, J0861, r, i, z) were used to compose the blue, green and red components, respectively.The FoV is ≈ 0.1 deg 2 .The photometric calibration procedure adopted in J-PLUS (see López-Sanjuan et al. 2019 for details) was not available at the time that our data was processed, and the photometric zeropoints (ZPs) have been obtained separately.ZPs for the bands g, r, i, and z were obtained from stars in the field observed with the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS; Chambers et al. 2016).ZPs for the remaining bands were computed using standard and secondary stars in the field.The bands g and r had the ZPs derived from both procedures, yielding small offsets of 0.087 and 0.069, respectively.Due to the impossibility of inferring the offsets to all bands we do not apply corrections to the magnitudes, using them as provided.These offsets act as adding a small difference between the broad and narrow-bands when SEDs are built, which we expect to cause an increase in χ 2 when performing the SED fitting (see section 4.2).
In the current work, we assume that the field is locate at high galactic latitude (as in the case of J-PLUS) and thus is not significantly affected by extinction.Throughout this work, no correction for line-of-sight extinction was applied.The line-of-sight extinction in the direction of NGC 1023 is estimated to be E(B-V) = 0.052 (following Schlafly & Finkbeiner 2011).

Methodology
We analyze the data presented in Section 2 in two processes, first by compiling a list of GC candidates around NGC 1023, and then by analyzing their stellar population properties.

Detection of candidates with GCFinder
To detect GC candidates in J-PLUS images, we develop the semiautomatic pipeline GCFinder.GCFinder consists of a python code that detects compact sources in a white image (a sum of frames of the 4 broad-bands g, r, i, and z) and performs a selection of GC candidates based on data quality, shape, color and magnitude criteria.The code uses Source Extractor (Bertin & Arnouts 1996) and Montage (Berriman et al. 2004), and run inside a support folder that also contains necessary files to use GCFinder.
A detailed description of GCFinder and technical requirements are given in Appendices A and B. In Appendix C we present the main different methodologies tested for the detection of GC candidates in J-PLUS data before choosing the strategy deployed in GCFinder.
An advantage of GCFinder is that it is a very light code that does not require the modeling of the host galaxy or smoothing filters before the detection of sources.The code is flexible and can be easily adapted to other photometric surveys besides J-PLUS, including the detection of other stellar systems such as Ultra Compact Dwarf Galaxies (UCDs, Phillipps et al. 2001).We run GCFinder on the data described in Section 2 and obtained a set of GC candidates, that are presented in Table 2.We note that a larger number of GC candidates are recovered in redder bands, in agreement with what is expected from old stellar objects.We also note that the S/N of blue bands, if compared to red bands tends to be lower in J-PLUS.
To measure the detection efficiency, i.e. how many GCs the method can recover in the extended light region of the host galaxy's halo, we compare our detections with the catalog from Kartha et al. (2014), used as reference.
The NGC 1023 reference catalog of Kartha et al. (2014) originally had 627 objects including faint sources.We note that the work done in Kartha et al. (2014) use data from MegaCam (Lenzen et al. 2003) at Canada-France-Hawaii Telescope (CHFT), which has a limiting magnitude fainter than T80Cam.In the gband, we are able to detect objects up to a magnitude of around 24 mag and, for comparison purposes, we exclude from the reference catalog objects fainter than the limiting magnitude of our data.
This results in about 200 GC candidates that could be detected within the limiting magnitude of J-PLUS.
Comparing the positions of our GC candidates with the reference catalog (see Figure 2), we retrieve 188 GCs in common which means that GCFinder detected a large fraction (about 93 %) of possible GC candidates.
Based on the numbers presented, we estimate that GCFinder does not detect about 7 % of expected GCs when compared to the literature.
We attribute this difference to the fact that Kartha et al. (2014) detect many globular clusters near the central region of the galaxy since its light central body hampers our ability to detect compact sources, even if modeled (see Appendix C).
On the other hand, we detect 335 new GC candidates with GCFinder due to: 1) new outer halo GCs were identified taking advantage of the wide field of J-PLUS images and 2) reduction in the contamination from Milky Way stars, which was possible from the adopted methodology (see Section A).From these 335 new GC candidates, 48% of them are located beyond the studied region in previous papers such as Kartha et al. (2014) and the rest of them are in the same halo region explored in the literature but are on average fainter than the GC candidates located at same region reported in Kartha et al. (2014) .

Stellar population properties
We obtain stellar population properties for our GC candidates via SED fitting.We use the codes DynBaS and TGASPEX (Magris C. et al. 2015;Mejía-Narváez et al. 2017), adapted to work with the J-PAS, J-PLUS and S-PLUS filter systems as illustrated in González Delgado et al. (2021) for mini-JPAS, where TGASPEX was used.Both DynBaS and TGASPEX are non-parametric fitting codes, in which the star formation history is expressed as an arbitrary superposition of different simple stellar populations.This is the first time that the stellar population properties of GC candidates are derived from J-PLUS data.To investigate the effect of adding narrow-bands to the SED fitting, we perform the fits for 3 different combinations of filters: (a) using only the broad-band filters, (b) using only the narrow-band filters, and (c) using all available filters.
In DynBaS and TGASPEX the best-the fitting solution is obtained by computing the non-negative values of the coefficients x tZ that minimize the merit function used to measure the goodness-of-fit.In Equation 1, F obs λ and f λ,tZ are the observed and model flux in each of the bands, respectively, and σ 2 l is the corresponding uncertainty.The sum is done over all the filters (index λ) and the N models (indices t and Z).N is equal to the number of time steps × the number of metallicities used in the fit.As described by Magris C. et al. (2015), in DynBaS we use N = 3 (hereafter DynBaS3), and minimize χ 2 in Equation 1requiring that the 3 derivatives ∂ ∂x tZ χ 2 = 0.This results in a system of 3 equations with 3 unknowns that we solve using Cramer's rule for all possible combinations of 3 model spectra.The DynBaS3 solution is then the one with the minimum χ 2 , subject to the condition x tZ ≥ 0. In TGASPEX we use the non-negative least squares (NNLS) algorithm (Lawson & Hanson 1974) to find the vector x tZ that minimizes χ 2 .In both codes, we use an outer loop to minimize by the dust attenuation τ V .
We remark that both DynBaS3 and TGASPEX provide independent, deterministic solutions, as opposed to statistical solutions, to the SED fitting problem.The TGASPEX solution, in general, with N ≥ 3, contains the DynBaS3 solution.As has been shown in Magris C. et al. (2015) and Mejía-Narváez et al. (2017), the DynBaS3 and TGASPEX solutions are consistent within errors.González Delgado et al. (2021) show that the TGASPEX solution is consistent with the solutions obtained by other SED fitting codes, including codes that use a fully Bayesian approach.Prieto et al. (2022, in prep) show that the deterministic DynBaS3 and TGASPEX solutions are consistent with those derived following a Bayesian treatment for both methods.In this paper, we opt for the deterministic approach for simplicity and because we know from the cited papers that the Bayesian approach does not add new insight into our problem.
Following González Delgado et al. (2021), we use the frequentist approach of Monte Carlo-ing the input (by adding Gaussian noise with observationally defined amplitudes) and repeating the fit many times (≈ 1000) assuming that the errors in the different bands are uncorrelated, to perform a statistical analysis based on the observed photon-noise and its impact on the results.A probability distribution function (PDF) for each stellar population property is then built by weighting the results from each iteration by the likelihood L ∝ exp(−χ 2 /2).The inferred value for the property of each GC is then obtained directly from the corresponding marginalized PDF.In the end, each population property is characterized by its best value derived directly by DynBaS3 and TGASPEX, and the mean, the median and the percentiles defining the confidence interval in the distribution.The best value is then considered as the best estimate of the unknown true value, and its precision is obtained by averaging the precision determined for each cluster from the PDF built as indicated above.The following properties are given as result: -Total stellar mass (M): Total mass of the stellar population.
It is calculated directly from the mass converted into stars according to our solutions for the GC, reported as log(M/M ⊙ ).-Luminous stellar mass (M ⋆ ): Stellar-mass of the stellar population at present.It is calculated from the mass converted into stars reduced by the mass lost by stars during their evolution, reported as log(M * /M ⊙ ).-Age of the stellar population: We define the mass-weighted logarithmic age (hereafter mass-weighted age) following Cid Fernandes et al. (2013, eq. 9) as where µ tZ is the fraction of mass of the base element with age t and metallicity Z, reported as < log age(yr) > M .Similarly, the light-weighted logarithmic age (hereafter light-weighted age) is defined as where F tZ is the fraction of light in the r filter corresponding to the base element with age t and metallicity Z, reported as < log age(yr) > L .-Metallicity of the stellar population: We define the massweighted logarithmic metallicity (hereafter mass-weighted Z) as reported as < [Z/Z ⊙ ] > M , and the light-weighted logarithmic metallicity (hereafter light-weighted Z) as We define the precision (mean standard deviation, σ j ) for the j th stellar property as where p84 j,i and p16 j,i are, respectively, the percentiles 84 and 16 of the PDF of the j th property for the i th GC candidate.
Example of SEDs chosen randomly are presented in Figure 3 for illustration purposes.

Color distributions
The study of GC systems in massive galaxies has shown that they can have a bimodal optical color distribution (e.g.Larsen et al. 2001;Peng et al. 2006;Kundu & Zepf 2007).At the same time, several studies have shown that the color-metallicity relation is highly non-linear (e.g.Yoon et al. 2006,Kim et al. 2021).
We study the distributions of all possible color-color diagrams with the J-PLUS filter system without mixing broad and narrow-band filters (as broad-and narrow-band filters were calibrated using different methodologies, see Section 2).A total of 25 colors were inspected, 10 using broad-band and 15 colors on narrow-band filters (Figures 4 and 5, respectively).Gaussian mixture modeling (GMM) was performed on the distributions using the Python library Scikit-Learn (sklearn, Pedregosa et al. 2011), following the procedures from Ivezić et al. (2014).We compare the GMM results obtained for one component (black curve) and two components (purple curve), using the Bayesian information criterion (BIC).The BIC makes assumptions about the likelihood that aims to simplify the calculation of the odds ratio and is useful to estimate the statistical significance of clusters found in the data.Therefore we assume that lower BIC values are associated with highly significant clusters, following Ivezić et al. (2014, chapter 5.4).
According to BIC statistics, we find evidence of color bimodality in 17 colors, namely u − g, u − r, u − i, g − r, r − z, i − z, In these cases, the BICs of the GMM with two components have a lower value than the BICs of the GMM with one component.We note that even though the BICs of the GMM with two components are lower than the values for one component in the cases of u − z,r − i, J0410 − J0430 and J0410 − J0861, the differences of the BICs are too small to be conclusive.In Figures 4 and 5 we show BIC values for GMM with one and two components as an example.
De Souza et al. (2017) favors the use of a regularized version of BIC, namely the Integrated Complete Likelihood (ICL).As a sanity check, we repeat the analysis evaluating possible color bimodality using the ICL criterion.According to ICL statistics, we find evidence of color bimodality in 10 colors, namely g In the case of the color J0378 − J0660, the ICL of the GMM with one component has a lower value than the GMM with two components, but the difference is small, hence we consider this case inconclusive.
A table summarising the results from both statistics is presented in Appendix D, Table D.1.

Stellar population properties
In this section, we present the results of the SED fitting performed on a sub-sample of 171 GC candidates.This sub-sample consists of only GC candidates with measured magnitudes in all 0.17 bands.The codes and models used are those described in Section 3.2.The uncertainties of the fits were estimated according to Equation 6 in Section 3.2, extending over the N = 171 GC candidates.In Table 3 we list the resulting values of σ j .The fits reported were performed for A V = 0, the value of <σ(A V )> listed in Table 3 gives an indication of the error of this assumption.
Stellar masses: Figure 6 shows that the distributions of the total mass obtained with TGASPEX and DynBaS3 are consistent, and range from below 10 3 to above 10 6 M⊙, with a peak around 10 5.5 M⊙.Given that GC masses are known in the range 10 4 to 10 6 M⊙ (Brodie & Strader 2006), we interpret the tail towards low mass to be indicative of contaminants (false positives) in our cluster-candidate catalog.Such contaminants represent a negligible fraction, only 5 % of our sample of GC candidates.It is unclear at this point if the more massive systems (stellar mass > 10 6 M⊙) are GCs or UCDs (Phillipps et al. 2001).These lowmass candidates are also intrinsically fainter, which results in lower S/N SEDs.Regarding the distributions retrieved from different filter sets, including the narrow-band filters in the fits has the effect of broadening the distribution of the best mass values.

Ages:
Figure 7 shows the distributions of our fit results for the mass-weighted and light-weighted ages.Old age clusters, with ages ≈ 10 Gyr dominate the distributions for all filter combinations (narrow-band only, broad-band only, broad and narrow-bands combined), with a secondary peak occurring at ages ≈ 9 Gyr.For reasons that remain unclear at this moment, this second intermediate age peak is less pronounced when only the narrow-bands are used in the fit.This may be related to the choice of sub-sample selected to be analyzed via SED fitting, where only the GC candidates with all filters measured were included.DynBaS3 retrieves a higher fraction of old GCs than TGASPEX.We interpret the tail towards youngest ages (log age < 8) as being produced by contaminants present in our candidate catalog.
Metallicities: Figure 8 shows the distributions of our fit results for the best values of mass-weighted Z and light-weighted Z.A clear bimodal distribution of metallicities is seen in most combinations of code and filter set, while the results of light and mass-weighted Z from TGASPEX obtained only with narrowbands show three modes, indicating three different populations.
A stronger tail at low metallicities is derived when only the broad-bands are used.We interpret this result as evidence that the narrow-bands help to constrain the metallicities.It is unclear at this point how the metallicity distributions are bimodal when we only find evidence of color bimodality in part of the studied colors.A possible channel for this could be a non-linear color- As in other galaxies, we interpret the two families of metalpoor and metal-rich GCs to be related to two mechanisms or episodes of star formation (e.g.Brodie & Strader 2006), although we find evidence of three populations when using narrow-bands only and TGASPEX.
Stellar masses vs. ages We explore if there are correlations among the derived parameters.The only correlation identified was between the stellar mass and light-weighted ages, as illustrated in Figure 9.A correlation is also present between stellar mass and mass-weighted ages, albeit less clear.Pfeffer et al. (2018) present globular cluster models in the context of E-MOSAICS project.These models describe the formation, evolution, and the disruption of this class of objects.In their work they find that based on their models and simulations most low-mass clusters were disrupted at redshift 0, therefore they conclude that clusters with higher mass are more likely to survive until the present time, which results in old GC populations having higher characteristic mass when compared with younger GCs.Therefore, we attribute the relation found in this work to be caused by the same processes found in Pfeffer et al. ( 2018), but we also note that we do not calculate ages for all GCs in NGC 1023, therefore our results could be affected by selection effects that are not well characterized.

Specific frequency
The specific frequency (S N ) of the GC population of a galaxy represents the total number of GCs per unit of host galaxy luminosity.Following Kartha et al. (2014), we adopt M V = -21.07± 0.06 and we calculate using the GCLF a N GC = 553 ± 60, which was used to determine the S N .We calculate a S N = 2.1 ± 0.2, which is consistent with the S N = 1.8 ± 0.2 reported in Kartha et al. (2014) and with S N = 1.7 ± 0.3 presented in Yong et al. (2012).Our S N is also consistent with estimations for lenticu- Fig. 5. Color distribution computed using narrow-band filters only.The curves represent the unimodal (black) and bimodal (purple) distributions returned by the GMM analysis.We also show in blue and red the two peaks that compose the bimodal distribution.BICs values for GMM with one and two components are shown as an example.lar galaxies (2 ≤ S N ≤ 6, Kundu & Whitmore 1998;Elmegreen 2000).

Discussion
In this section we discuss the results found in Section 4, to connect the observed properties of the GC system with the evolution of NGC 1023.

The accretion history of NGC 1023
The fact that we can identify bimodal distributions in metallicities from the SED fitting analysis from Section 4.2 could be evidence that there are at least two subpopulations of GCs in the galaxy.
Li & Gnedin (2019) use a novel cluster formation model on a simulated galaxy of the same size as the Milky Way and observed that GC candidates tend to form during major merger events.The merger-induced GC formation scenario has been dis- cussed in various recent articles (Li & Gnedin 2014;Choksi et al. 2018;Choksi & Gnedin 2019).
We can further investigate this scenario by studying agemetallicities relations.In Figure 10 we compare the agemetalliticy relation (hereafter, AMR) obtained from our results to 3 E-MOSAIC simulations (Pfeffer et al. 2018) from Kruijssen et al. (2019a).We choose 3 simulations that have halo masses comparable to the halo mass of NGC1023 adopting the masses from Alabi et al. (2016) and Bílek et al. (2019).
The comparison in Figure 10 is useful to get a handle on the epoch of GC assembly in the NGC 1023.We note that the AMR of MW023 is in better agreement with the AMR of NGC 1023 GC candidates.We also note that our AMR has outliers, similar to the ones found in MW014.Our results then favor the accretion histories of simulations such as MW014 or MW023 and rule out histories such as the one in MW016.The AMR of GC candidates when compared with the simulations indicate that NGC 1023 likely experienced an initial and rapid phase of star formation that might have formed the majority of the GC candidates, since a big amount of GC candidates were formed early in galaxy evolution considering that the age distribution has a dominating peak at ≈ 10 10 yr, agreeing with results from Kruijssen et al. (2019a).
We also note a likely broad AMR, where the GC candidates have a wide range of metallicities.Nevertheless, there is a caveat regarding our interpretation of a broad AMR, as the selection function introduced by GCFinder is not well characterized.As such it remains an open question if the color cuts applied by the pipeline would introduce distortions in the age-metallicities relation.Kruijssen et al. (2019a) find that a wide range of GC metallicities was related to a wide range of progenitor masses.Therefore, we believe our possible broad age-metallicity relation and our wide range of metallicities imply that NGC 1023 experienced mergers and accretion events in the past, resulting in more than one episode of intense star formation.This is in sync with what is been discovered about the formation of the Milky Way.Studies about the formation of our Galaxy are motivated by many surveys created in the last years and generating huge amounts of data.In particular, Gaia survey (Gaia Collaboration et al. 2016a,b;Brown et al. 2018Brown et al. , 2021) ) have been revolutionising our understanding about the Milky Way.Many articles in recent years making use of Gaia were published identifying stars in the Milky Way that are claimed to be accreted from dwarf galaxies that no longer exist.In particular, stars have been claimed to be born in the progenitor galaxies TGASPEX, narrow bands Fig. 9. Correlation between ages and masses derived via SED fitting.A trend where the mean ages of the cluster-candidates increase with stellar masses is seen for all combinations of code and filter set.Left and right illustrate results for DynBaS 3 and TGASPEX respectively.Red triangles represent results obtained using all bands, green circles represent results when using only broad-bands and blue stars represent results obtained from narrow-bands only.
Gaia-Enceladus (Helmi et al. 2018;Belokurov et al. 2018;Das et al. 2020) (which is believed to be the last major merger experienced from the Milky Way), from the Sequoia progenitor (Myeong et al. 2019), from Thamnos 1 and Thamnos 2 (Koppelman et al. 2019), from a structure in the inner Galaxy (Kruijssen et al. 2019b(Kruijssen et al. , 2020;;Horta et al. 2021), among others.We note that there is evidence from kinematic studies (e.g.Romanowsky et al. 2012;Villaume et al. 2019) and simulations (e.g.Muratov & Gnedin 2010;Choksi et al. 2018) about accretion events from dwarf galaxies and that the hierarchical assembly of GC systems is well accepted.In particular, the fact that we only find evidence of color bimodality in some cases is not surprising.GCs trace assembly histories of galaxies and galaxies likely undergo many minor and possibly major mergers throughout their life.In this case, the lack of strong evidence of color bimodality for some colors could be an indicator that more than two subpopulations exist, but we are not able to disentangle them.Puzia et al. (2002) and Blom et al. (2012a,b), for example, found three subpopulations of GCs in the galaxy NGC 4365.

The ongoing interaction with NGC 1023A
NGC 1023 is a barred galaxy in an ongoing interaction with NGC 1023A, a small companion at the outskirt on NGC 1023 at East (Debattista et al. 2002).NGC 1023A was recognised as an independent galaxy by Barbon & Capaccioli (1975) and designated as NGC 1023A by Hart et al. (1980).Capaccioli et al. (1986) classified it as a Magellanic irregular or late-type dwarf galaxy, while Sancisi et al. (1984b) from HI observation found a complex kinematic.On the other hand, Debattista et al. (2002) found a faster bar pattern speed, which would be not compatible with a scenario of a recent formation of the bar by the interaction with NGC 1023A.Nevertheless, the ongoing interaction may have an effect on the overall NGC 1023 structure as well as on the GCs distribution.A small fraction of the GC population of the NGC 1023 system could be associated with NGC 1023A (Cortesi et al. 2016).Due to the morphology and luminosity of NGC 1023A, a definition of its center is challenging, but we estimate that NGC 1023A is at a projected distance of approximately 6.7 kpc from NGC 1023.
Since NGC 1023 is an early type galaxy the tidal effect could have two different dynamics answers depending on how longer or shorter the encounter time is compared to the galaxy's internal cross-time (e.g, Aguilar & White 1985, 1986;Binney & Tremaine 2008).For the outer parts the crossing time could be larger than the encounter time, therefore they may suffer an impulse response.On the other hand, for the inner parts, the crosstime could be smaller than the encounter time, thus they may suffer a typical tidal response.As a consequence, the central parts can be displaced with respect to the outer ones.These offsets can be up to 20% of the observable radius of the galaxy (Lauer 1986(Lauer , 1988;;Davoust & Prugniel 1988;Combes et al. 1995;González-Serrano & Carballo 2000;Mora et al. 2019;Buzzo et al. 2021).Buzzo et al. (2021) show that the nucleus of the lenticular galaxy NGC 3115 has a displacement of 160 arcsec with respect to the outer parts.They interpret it as the result of a recent pericenter passage with its close small companion.To probe if NGC 1023 has a similar feature, we perform a similar photometric analysis following Mora et al. (2019) and Buzzo et al. (2021).We model the isophote contours, in r−band, by using the ELLIPSE task from IRAF (Jedrzejewski 1987), we let free the position angle, ellipticity, and centre of the ellipses.To quantify the offset of the isophotes, we take as reference the photometric center of the galaxy.In Fig. 11, we show the radial profile of the offsets at the top panel, while the outermost isophotes and their respective fitted ellipses with their centers are plotted at the bottom panel.It is clear that from ∼100 arcsec the central part starts to move toward the East-South with respect to the outer parts, the maximum offset is about ∼700 pc.This nuclear displacement is strong evidence that NGC 1023 and NGC 1023A had recently a pericenter passage, just a few hundred million years ago (Combes et al. 1995;Mora et al. 2019).The orientation of the offset could be used as a strong constraint in a numerical simulation of the dynamic encounter of this pair, since the central part of NGC 1023 must have headed into the East-South direction at the pericenter passage.This would limit the family of possible orbits to model the system.(e.g, Combes et al. 1995;Mora et al. 2019).
How this interaction could have affected the distribution of the GCs of NGC1023?To address this question, we calculate the photometric center of the GCs candidates and overlay it on the bottom panel in Fig. 11.We can see that the photometric center of the GCs candidates follows the same direction as the centers of the outer isophotes.The displacement of the nucleus region with respect to the GCs is ∼ 1600 pc.This behavior is expected according to impulse theory (Aguilar & White 1985, 1986;Binney & Tremaine 2008) since the GCs belong to the galactic halo then they have the largest cross times of the galaxy, hence the nuclear displacement should be the largest one.
In addition to this photometric analysis, we study the residual image from the ellipse model, see Fig. 12.The residual image unveils NGC 1023A and the bar (oriented North-East to South-West) of NGC 1023 (Möllenhoff & Heidt 2001;Debattista et al. 2002), besides two possible relic like-spiral arms.The bar has a radius of ∼1100 pc.We note that here we are reporting the presence of these relic spiral-like arms for the first time.It is very plausible that the origin of these structures is also due to the interaction with NGC 1023A, in this case, they would be tidal structures.However, one of the formation mechanisms of lenticular galaxies is gas removal from a spiral galaxy, then these relic spiral-like arms could be a memory of the progenitor galaxy of NGC 1023.These features can also serve as dynamical constraints for a numerical simulation of the encounter.

Summary and Conclusion
In this work, we present the first study about extragalactic globular clusters using J-PLUS data.As a test case, we detect and study the GC system in NGC 1023 with the 12 bands of J-PLUS.To detect the GC candidates we develop GCFinder, a code that can be applied to current and upcoming wide-field multi-band surveys, such as J-PAS and S-PLUS.The pipeline presents good results within the characteristics of the survey for which it was designed.The end product is a code that can be adapted to other photometric surveys and other types of compact stellar systems, such as ultra-compact dwarf galaxies.
With the catalog of GC candidates provided by GCFinder, we perform a study of the stellar population content of a subsample of objects, using SED fitting techniques and photometry from broad and narrow-band filters.To calculate masses, ages, and metallicities of the GC candidates we use the codes DynBaS3 and TGASPEX adapted to work with the J-PLUS filter system.We also carefully model the light of NGC 1023 to investigate possible displacement between the outer isophotes and the distribution center of the GC candidates, which is useful to understand galaxy evolution.
In the following lines we summarize our main findings: -We identify 523 GC candidates in NGC 1023 using GCFinder, being 335 of them not yet reported in the literature.A significant part of these new GC candidates is located in the outer regions of NGC 1023 since we took advantage of the wide field of view of J-PLUS images.We find a specific frequency of S N = 2.1 ± 0.2, which is consistent with estimations for lenticular galaxies in the literature.-We investigate color distributions of the GC candidates, exploring the novel colors provided by J-PLUS.According to BIC statistics, we find evidence of color bimodality in 17 colors -We obtain masses, ages, and metallicities for 171 GC candidates from SED fitting.We find that the peak of the mass distribution is at 10 5.5 M ⊙ .We also find a tail of GC candidates with low masses, which we interpret as likely contaminants in our list of candidates.It is unclear at this point if the more massive systems (stellar mass > 10 6 M⊙) are GCs or UCDs.-The mass-weighted and light-weighted age distributions cover a wide range of ages, with a dominant population at ≈ 10 10 yr.The mass-weighted and light-weighted metallicity distributions of the GC candidates are bimodal in most of the cases (combining different filter sets and codes) and, in a minority of cases, we find three peaks.We note that the inclusion of narrow-band filters helps to constrain the metallicities.These results indicate that there are subpopulations of GCs in NGC 1023 that could exist due to accretion events or different epochs/mechanisms of star formation in the galaxy.-We identify a correlation between light-weighted ages and stellar masses, where older GCs tend to be more massive.This suggests that massive GC candidates in NGC 1023 are more likely to survive the turbulent history of the host galaxy than less massive objects, which is in agreement with the literature on GC systems of galaxies.-The age-metallicity relation is likely broad.A comparison with simulations shows evidence of a likely initial rapid phase of star formation, responsible for the formation of the majority of the GCs.Articles in the literature (e.g.Kruijssen et al. 2019a) found that a wide range of GC metallicities is related to a wide range of masses of progenitor galaxies.Therefore, the broad AMR we find is also evidence of past accretion events experienced by NGC 1023.-We detect that the photometric center has a displacement of ∼700 pc and ∼1600 pc with respect to the outer isophotes and the GC candidate distribution center, respectively.The offsets could be the result of ongoing interaction between NGC 1023 and NGC 1023A.These effects are in excellent agreement with impulse theory (Aguilar & White 1985, 1986;Binney & Tremaine 2008).-The residual map obtained from the photometric model of NGC 1023 unveils two spiral-like arms.These structures are probably due to the NGC 1023 interaction with the NGC 1023A satellite galaxy.
From our main findings, we conclude that it was possible to retrieve new and valuable information about the evolutionary past of NGC 1023 as observed by J-PLUS.The multiple GC populations, relic spiral arms, and displacement between the photometric center of the GC candidates and the isophotal center that we report in this work support a formation of NGC 1023 that involved several minor mergers and group harassment, causing a transformation from spiral to the nowadays lenticular galaxy.Kartha et al. (2014).Top right panel: GCs detected using the Source Extractor method that are also presented in Kartha et al. (2014).Bottom left panel: GCs detected using the median smoothing method that are also found in Kartha et al. (2014).Bottom right panel: Objects detected using ISOFIT and CMODEL method (Ciambur 2015) that are also reported in Kartha et al. (2014).In the background we show the white image of NGC 1023 used for the detection of sources.The FoV is ≈ 0.1 deg 2 .

Fig. 2 .
Fig. 2. In the background we show the white image of NGC 1023 used for the detection of sources.The FoV is ≈ 0.1 deg 2 .Red circles represent GC candidates from the reference catalog (Kartha et al. 2014).Blue circles represent GCs in common between GCFinder selection and the reference catalog.

Fig. 3 .
Fig. 3. Examples of SED fitting for 3 GC candidates.The observed and synthetic SEDs are presented in black and blue, respectively.The model fits of the two codes are very similar, therefore they appear superposed in the figure.

Fig. 4 .
Fig. 4. Color distribution computed using broad-band filters only.The curves represent the unimodal (black) and bimodal (purple) distributions returned by the GMM analysis.We also show in blue and red dashed lines the two peaks that compose the bimodal distribution.BICs values for GMM with one and two components are shown as an example.

Fig. 6 .Fig. 7 .
Fig.6.Distribution of stellar masses obtained from SED fitting.The left panel illustrates results from DynBaS 3, while the right panel shows results of TGASPEX.The distributions show a broad peak with log Stellar Masses between 5 and 6 M ⊙ .Red represents results obtained using all available filters.Green represents results obtained using broad-band filters only.Blue represents results obtained using narrow-band filters only.

Fig. 8 .
Fig.8.Distribution of metallicities obtained from SED fitting.Upper and lower panels illustrate results for DynBaS 3 and TGASPEX respectively.Left and right panels show results for mass-and light-weighted Z, respectively.Red represents results obtained using all available filters.Green represents results obtained using broad-band filters only.Blue represents results obtained using narrow-band filters only.

Fig. 10 .
Fig.10.The AMR for the NGC 1023 GC system based on DynBaS 3 (blue circles) compared to 3 E-MOSAIC simulations: MW014, MW016 and MW023 fromKruijssen et al. (2019a) (grey contours) that have comparable halo masses to NGC 1023.We observe a likely broad AMR of NGC 1023 GC system, where the objects have a broad metallicity distribution and are mainly old.

Fig. 11 .
Fig. 11.NGC 1023 displaced nucleus.Top panel: the offset radius profile with respect to the photometric center of the galaxy.Bottom panel: selected outer isophotes (in black) with their respective fitted ellipses (blue) are overplotted on r-band image.The blue pluses are their respective centers, while the red point is the photometric center of the galaxy and the green point is the photometric center of GCs candidates.The top-right inset is a zoom-in of the inner part to see better the displacements.The blue points in the offset radius profile are the offsets for the selected isophotes.

Fig. 12 .
Fig. 12. Photometric model for NGC 1023.Left panel: The r-band image.Middle panel: the ellipse model.Right panel: the residual map.NGC 1023A is reveled in this map, together with two spiral arms and an inner bar.

Fig. C. 3 .
Fig. C.3.Top left panel: GCs present in the reference catalog fromKartha et al. (2014).Top right panel: GCs detected using the Source Extractor method that are also presented inKartha et al. (2014).Bottom left panel: GCs detected using the median smoothing method that are also found inKartha et al. (2014).Bottom right panel: Objects detected using ISOFIT and CMODEL method(Ciambur 2015) that are also reported inKartha et al. (2014).In the background we show the white image of NGC 1023 used for the detection of sources.The FoV is ≈ 0.1 deg 2 .

Table 1 .
FWHM values and exposure times for NGC 1023 data in each filter.

Table 2 .
Number of GC candidates selected per band after the selection criteria are applied.