© The Authors 2026

1. Introduction

The stellar mass function (SMF) is a fundamental tool for understanding the properties and evolution of galaxy populations, encompassing both quiescent (Q) and star-forming (SF) systems. It represents the number density of galaxies as a function of their stellar mass, M_★, determined by counting galaxies in stellar mass bins in a given cosmological volume and correcting from biases. The SMF provides critical insights into the evolution of the star formation rate (SFR) over cosmic time, enabling us to track the accumulation of stellar mass across different epochs. Large redshift surveys, such as the Sloan Digital Sky Survey SDSS (York et al. 2000), zCOSMOS (Lilly et al. 2007), and the Galaxy And Mass Assembly (GAMA) survey (Baldry et al. 2010), have significantly advanced our ability to study galaxy populations with unprecedented precision. These surveys facilitate the derivation of the SMF and its evolution across cosmic time, substantially contributing to our understanding of galaxy formation and evolution (see also Muzzin et al. 2013, Peng et al. 2010, and Kelvin et al. 2014). More recent studies, including those from GAMA (Wright et al. 2018; Driver et al. 2022), have further refined and extended these analyses.

The SMF is not only pivotal for observational studies but also serves as a crucial constraint for theoretical models of galaxy formation. In these models, the SMF works as a key diagnostic tool for testing model predictions and refining simulation parameters. A primary challenge in these models is reconciling the halo mass function from simulations with the observed SMF from galaxy surveys. Addressing discrepancies, especially at the low- and high-mass ends, requires iterative refinements through feedback mechanisms, such as supernovae and active galactic nuclei (AGN) feedback, which regulate gas cooling and align the predicted most massive end of the SMF with observed data (Kauffmann et al. 1993; Somerville et al. 2008). Up to z ∼ 7.5, the SMF is well described by a single or double Schechter function (Schechter 1976; Weaver et al. 2023), featuring an exponential decline at a characteristic stellar mass (M^*) and a low-mass slope (α). While the low-mass slope shows only mild evolution, several studies report a moderate but significant increase in M^* from z ∼ 2 to the present day (Marchesini et al. 2009; Tomczak et al. 2014; Adams et al. 2021), likely driven by mass growth through mergers. This suggests that the mechanisms regulating stellar mass assembly vary across cosmic time. In contrast, the normalization of the SMF (Φ^*), especially for SF galaxies, evolves more strongly and reflects the decline in cosmic SFR density (Popesso et al. 2023). Further discussion of these trends can be found in Weaver et al. (2023), Wright et al. (2018).

Despite these advancements, the underlying physical mechanisms governing the SMF, even in the low-redshift Universe, remain poorly understood. Processes such as galaxy–galaxy mergers, stellar and gas kinematics, gas inflows and outflows, and feedback from supernovae and massive black holes are active areas of research. These mechanisms play a crucial role in the transformation of SF galaxies into passive ones, making them a subject of intense investigation in contemporary astrophysics. Studying the SMF separately for SF and Q galaxies helps reveal the physical processes that quench star formation, which are not used to calibrate theoretical models but provide valuable tests of their predictions. Two primary mechanisms are frequently invoked to explain the quenching of star formation: AGN feedback and environmental effects. In the case of AGN feedback, the accretion of matter by a central supermassive black hole injects energy and momentum into the surrounding medium, heating the gas and suppressing star formation. Environmental effects, particularly in high-density regions such as galaxy clusters, can also lead to quenching. Mechanisms such as gas stripping and galaxy interactions can inhibit star formation. While AGN feedback is most efficient in massive galaxies (M_★ >  10^10.5 M_⊙) (Piotrowska et al. 2022), environmental factors primarily influence dwarf galaxies (M_★ < 10⁹ M_⊙; Peng et al. 2010).

In this study, we used data from the Javalambre Photometric Local Universe Survey (J-PLUS¹; Cenarro et al. 2019) third data release (DR3) to estimate the SMF in the local Universe. J-PLUS DR3 covers a sky area of 3284 square degrees and employs a filter set comprising the five broad bands (ugriz) and seven narrow bands, thereby providing improved photometric redshifts (z_phot) and SFRs with respect to previous broadband surveys. A sample of 890 844 galaxies at redshift 0.05 ≤ z ≤ 0.20 and with r ≤ 20 mag was used to estimate the SMF for both SF and Q galaxies. Such a large sample size is critical for conducting robust statistical studies of galaxy populations.

This paper is structured as follows. Section 2 describes the dataset and sample selection. Section 3 outlines the methodology, including quality cuts, spectral energy distribution (SED) fitting, and the 1/V_max method for deriving the SMF in J-PLUS DR3 for Q and SF galaxies. Section 4 presents the results, including the SMF for Q, SF, and all galaxies, along with the Q fraction. The discussion and conclusions are provided in Sects. 5 and 6. Throughout this work, we adopt a Lambda cold dark matter (ΛCDM) cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, h = 0.7, Ω_m = 0.3, and Ω_Λ = 0.7.

2. Data and galaxy properties

2.1. J-PLUS DR3

J-PLUS is carried out at the Observatorio Astrofísico de Javalambre (OAJ²;Cenarro et al. 2014) using the 83 cm Javalambre Auxiliary Survey Telescope (JAST80). The telescope is equipped with a 9.2k × 9.2k pixel camera (Marín-Franch et al. 2015), providing a wide field of view of 2 deg². It features a set of 12 photometric filters, including the five SDSS (u,g,r,i,z) broadband filters and seven medium or narrowband filters. These additional filters target key stellar spectral features: four cover the region around the 4000 Å break (J0378, J0395, J0410, and J0430), one captures the magnesium doublet (J0515), another probes the calcium triplet (J0861), and another is centered on the Hα emission line at rest frame (J0660). These narrow band filters allow more accurate estimates of photometric redshifts and galaxy physical properties. Further details about the J-PLUS observing strategy, data reduction, and general goals can be found in (Cenarro et al. 2019). Consequently, J-PLUS provides broad optical coverage, enabling a variety of studies in stellar astrophysics (e.g., Bonatto et al. 2019; Whitten et al. 2019; Solano et al. 2019; López-Sanjuan et al. 2024a), galaxy evolution across different redshift ranges (e.g., Logroño-García et al. 2019; San Roman et al. 2019), and galaxy clusters (e.g., Molino et al. 2019; Jiménez-Teja et al. 2019). Additionally, it has been used to identify extreme extragalactic emitters (e.g., Spinoso et al. 2020; Lumbreras-Calle et al. 2022).

For this study, we utilized J-PLUS DR3, which was made publicly available in December 2022. J-PLUS DR3 includes photometric information for 1642 pointings, covering a total sky area of 3284 deg². This was reduced to 2881 deg² after masking bright stars, optical artifacts, and overlapping regions. The catalog and different tables used here are publicly accessible on the J-PLUS website³, which contains calibrated photometry (López-Sanjuan et al. 2024b) in several apertures for approximately 47.4 million objects detected in the r band, obtained using SExtractor’s dual-mode.

2.2. J-PLUS DR3 galaxy sample

To select galaxies, we used the point-like probability from table galclass.sglc_prob_star in (López-Sanjuan et al. 2019), estimated by combining the available priors and information in the gri broad bands. Sources with probabilities lower than 0.5 were selected as galaxies.

From the jplus.MagABDualObj table, we retrieved the magnitudes extracted using AUTO aperture (elliptical apertures with semimajor axis equal to twice the Kron radius of the source). The magnitudes of J-PLUS DR3 galaxies were corrected for Galactic extinction using the jplus.MWExtinction table, which contains B − V color excess due to Milky Way dust at the position source. The integrated E(B − V) was estimated from the Schlafly & Finkbeiner (2011) recalibration of the Schlegel et al. (1998) infrared-based dust map. The lower redshift limit z = 0.05 was chosen to minimize the effects of aperture and surface-brightness and to reduce the impact of large-scale local structures on the measured SMF. The photometric redshifts are available in the jplus.PhotoZLephare table. These redshifts were obtained using LePhare (Arnouts & Ilbert 2011). LePhare estimates the photometric redshifts in J-PLUS DR3 by scanning from z_min = 0.0 to z_max = 1.0 in steps of z_step = 0.005. At each redshift, it computes a likelihood based on the χ² of the best-fitting template to the J-PLUS data. The CEFCA_minijpas library used includes 50 synthetic galaxy spectra produced with the Code Investigating GALaxy Emission (CIGALE; Boquien et al. 2019). A prior derived from the VIMOS VLT Deep Survey (VVDS) galaxy counts (Le Fèvre et al. 2005) was applied to obtain the final redshift probability distribution. For details on the configuration, templates, and priors, we refer to Hernán-Caballero et al. (2021).

The ADQL query used to access the J-PLUS DR3 data used in this work is presented in Appendix A, with a total of 890 844 galaxies with extinction-corrected r₀ ≤ 20 mag at 0.05 ≤ z_phot ≤ 0.20 and MASK_FLAGS = 0 (Fig. 1). The selected redshift and magnitude ranges ensured enough signal-to-noise in the J-PLUS photometry and a proper covering of the 4000 Å break with the bluest medium-band filters to derive high-quality physical properties for the galaxies (Sect. 2.3).

Fig. 1. Sky location of J-PLUS DR3 galaxies with r ≤ 20 mag at 0.05 ≤ zphot ≤ 0.20 covering 2881 deg2 (blue points). The black curve depicts the Milky Way.

The accuracy of the z_phot in the sample was estimated by comparing them with the spectroscopic redshifts (z_spec) from SDSS DR12, which has 686 176 counterparts. The difference

$$\begin{array}{c}\hfill \delta z=\frac{z_{\mathrm{phot}}-z_{\mathrm{spec}}}{1+z_{\mathrm{spec}}}\end{array}$$(1)

was computed. A negligible systematic offset of −0.002 and a 0.011 dispersion were measured from a Gaussian fit to the δz distribution. For comparison, the Gaussian distribution compared to the photometric redshifts of Beck et al. (2016) based on ugriz SDSS photometry yields a dispersion of 0.013, i.e., a 20% improvement when J-PLUS DR3 photometry is used.

2.3. Physical properties of galaxies with CIGALE

Spectral energy distribution modeling was performed using CIGALE, with the full configuration provided in Appendix B. For each galaxy, CIGALE was run at a fixed redshift, adopting the LePhare photometric redshift as an input. The star formation history was modeled using the sfhdelayed module, adopting a delayed-τ model with optional recent bursts. The main stellar population spans a range of three τ_main values from 100 to 5000 Myr and five ages from 3000 Myr to 13 Gyr, while the burst component includes three τ_burst values from 5 to 500 Myr, six ages between 5 and 1000 Myr, and four burst mass fractions up to 40%. The spectra of stellar populations were synthesized using the Bruzual & Charlot (2003) models with a Chabrier (2003) initial mass function (IMF) and three metallicities (Z = 0.004, 0.008, and 0.02). Nebular emission was included using the nebular module, with ionization parameters log U = −2.0 and −3.5, gas-phase metallicities Z_gas = 0.014 and 0.022, electron density n_e = 100 cm⁻³, and ionizing photon escape fractions up to 20%. Dust attenuation was modeled with the dustatt_modified_starburst module, following a Calzetti-like law (Calzetti et al. 2000) with five E(B − V) values up to 0.5, no UV bump, and a power-law slope of 0.0 (R_V = 3.1).

The derived parameters used in this study are the stellar mass of the galaxy (in M_⊙ units) and its SFR averaged over the last 10 Myr, in M_⊙ yr⁻¹ units. This SFR was used to cover the typical visibility timescale of the widely used Hα emission line of a star-formation burst. The specific SFR for each galaxy was defined as sSFR = SFR/M_★ [yr⁻¹]. Two examples are presented in Fig. 2. We observe that both fittings are successful, with reduced χ² ≈ 1. The table for the nearly 900 k galaxies analyzed is available in the J-PLUS DR3 repository and in VizieR.

Fig. 2. Examples of SED fitting to J-PLUS 12-band photometry using CIGALE. Panel (a): Q galaxy. Panel (b): SF galaxy. In each panel, the observed fluxes are shown as purple circles with error bars, and the best-fit model spectrum is plotted as a solid black line. The red dots indicate the model-predicted fluxes in each filter. The contributions from the attenuated stellar population (orange line), unattenuated stellar emission (dashed blue line), and nebular emission (green line) are also shown, although difficult to see in the left panel due to their resemblance to the model spectrum. Bottom subpanels: Relative residuals, defined as (Obs − Mod)/Obs, with a horizontal dashed line at zero for reference. A gri composite image from J-PLUS is inset in the bottom right corner of each panel.

The reliability of the derived galaxy properties was tested by comparison with the results from Duarte Puertas et al. (2017). They estimated the stellar masses, SFRs, and specific star formation rates (sSFRs) for a sample of 209 276 SF galaxies with spectroscopic information from SDSS. An empirical correction based on the Calar Alto Legacy Integral Field Area (CALIFA; Sánchez et al. 2012) spectroscopic data to the fiber measurements in SDSS was applied, providing an excellent reference for our photometric measurements. The difference between the Duarte Puertas et al. (2017) and J-PLUS measurements were defined as

$$\begin{array}{c}\hfill \delta \mathcal{X}=log{\mathcal{X}}^{\mathrm{DP}}-log{\mathcal{X}}^{\text{J-PLUS}},\end{array}$$(2)

where the parameters 𝒳 = {M_★, SFR, sSFR} were explored. The distribution of the differences for the 14 063 objects in common was approximated with a Gaussian distribution. The median (μ) and the dispersion (σ) of the Gaussian were used to determine the quality of the J-PLUS measurements.

The stellar masses present a difference of 0.13 dex and a dispersion of 0.2 dex (Fig. 3b). The difference is compatible with typical systematical uncertainties due to the use of different models in stellar mass estimates (∼0.3 dex, Barro et al. 2011). The stellar masses were further tested by applying the Taylor et al. (2011) mass-to-light relation, log(M_*/M_⊙) = 1.15 + 0.70 (g − i)−0.4 M_i, where M_i is the absolute magnitude in the rest-frame i band and (g − i) is the rest-frame color. Both were derived from CIGALE. This relation can be used to estimate the stellar mass with a 0.1 dex precision. We find an offset of −0.08 dex and a remarkable low dispersion of 0.03 (Fig. 4). The comparison with Duarte Puertas et al. (2017) includes the distance uncertainty from the photometric redshifts, while the comparison with Taylor et al. (2011) only accounts for differences in the mass-to-light ratio inference.

Fig. 3. Comparison between SFR (a), stellar mass (b), and sSFR (c) from this study and Duarte Puertas et al. (2017) for the N common sources. The red lines show the best Gaussian fits with median μ and dispersion σ.

Fig. 4. Comparison of stellar masses from this study and Taylor et al. (2011) for the N common sources. The red line shows the best Gaussian fit with median μ and dispersion σ.

The SFR and the sSFR are compared in Figs. 3a and c, respectively. The offsets obtained are −0.01 dex for SFR and −0.15 dex for sSFR, while the dispersions are 0.44 dex and 0.34 dex, respectively. The SFR is measured without bias, and the sSFR offset reflects the difference found in the stellar mass estimate. Because the distance uncertainty affects both the stellar mass and the SFR in the same way, the sSFR dispersion is lower.

The Duarte Puertas et al. (2017) sample is biased toward SF galaxies with r ≲ 18 mag (Appendix C). We divided the sSFR distribution in Fig. 3c by the uncertainty derived from CIGALE. The resulting distribution is Gaussian, with σ = 0.9, close to the expected unity. This implies that CIGALE errors are reliable and overestimated by just 10%. The median of the sSFR uncertainties from CIGALE in the full J-PLUS sample increases from 0.3 dex at 15 < r < 16 mag to 0.45 dex at 19 < r < 20 mag. We conclude that the stellar mass and the sSFR can be properly estimated from J-PLUS optical photometry alone to within a factor of 2–3 at r < 20 mag.

2.4. Quiescent and star-forming galaxies in J-PLUS DR3

The Q and the SF populations were selected using their position in the sSFR versus stellar mass plane. A 2D histogram of the 890 k galaxies under study, weighted by their odds parameter⁴, was used to define the best selection (Fig. 5). Two populations are clearly visible on this diagram, corresponding to SF and Q galaxies. The former peaks at log M_★ = 9.9 dex and logsSFR = −9.5 dex, the latter at log M_★ = 10.6 dex and logsSFR = −11.0 dex. The limit of both populations was determined by identifying the relative minimum of the sSFR histogram with a double Gaussian distribution. We obtained log sSFR_lim = −10.2 dex, with no dependence on stellar mass. In the following, Q galaxies are defined as those with an sSFR lower than sSFR_lim and as SF otherwise. The initial sample was split into 488 369 SF galaxies and 402 475 Q galaxies.

Fig. 5. sSFR vs. stellar mass 2D odds-weighted histogram. The two maxima correspond to SF (higher sSFR values) and Q (lower sSFR values) populations. The limit between both populations is set to log sSFR = −10.2 dex (horizontal dashed line).

Other approaches are possible for splitting the global population into Q and SF galaxies. For example, one can use the rest-frame (U − V) versus (V − J) color-color diagram (Williams et al. 2009) and their ultraviolet and infrared extensions (e.g., Leja et al. 2019). In addition, these diagrams can include dusty SF galaxies in the Q sample (Díaz-García et al. 2019). The analysis performed by Leja et al. (2019) suggests that the usual selection in the (U − V) versus (V − J) color–color diagram is equivalent to a log sSFR_lim = −10.0 yr⁻¹ selection, which is close to our inferred value.

The photometric and derived physical properties of the galaxy sample analyzed in this work were compiled into a catalog, a subset of which is shown in Table F.1. The catalog includes J-PLUS DR3 photometric measurements and stellar population parameters, such as stellar masses, SFRs, and rest-frame colors, estimated using CIGALE. In the next section, the methodology to determine the SMF for the total, SF, and Q populations is explained.

3. Estimation of the stellar mass function

Once we obtained good-quality stellar masses for Q and SF galaxies, we obtained the SMF of J-PLUS DR3. We followed the method described in Pozzetti et al. (2010) for estimating mass completeness, using the 1/V_max technique and bootstrapping to estimate errors.

3.1. Stellar mass completeness

We estimated the stellar mass completeness limit as a function of redshift, M_lim(z), following the method introduced by Pozzetti et al. (2010). This approach computes, for each galaxy, the minimum stellar mass for an apparent magnitude equal to the survey’s limiting magnitude. The resulting function M_lim(z) defines the stellar mass above which the sample is considered complete at each redshift. This correction is necessary in magnitude-limited surveys, which can include low-mass galaxies if they are sufficiently luminous, but can exclude massive galaxies that are too faint to be detected. As a result, the stellar mass completeness depends on both redshift and the distribution of mass-to-light ratios in the sample. Separate completeness limits were therefore derived for different galaxy populations (e.g., SF and Q). To determine the effective limiting magnitude of the J-PLUS DR3 sample, m_lim, we examined the distribution of the r apparent magnitudes in the redshift bins. In each bin, we deselected the faintest 20% of galaxies and retained the brightest among them. This yields a conservative estimate of the survey limit with m_lim = 20.1, consistent with our magnitude cut. Following Pozzetti et al. (2010), we computed for each galaxy i a limiting stellar mass M_lim, i, defined as the mass it would have at its redshift if its apparent magnitude were equal to the survey limit m_lim:

$$\begin{array}{c}\hfill log{M}_{\lim ,i}=log{M}_i+0.4\times (m_i-m_{\lim }),\end{array}$$(3)

where M_i and m_i are the stellar mass and observed magnitude of the galaxy. This relation assumes that luminosity L is proportional to flux F, and that M/L remains constant under rescaling. Next, we divided the sample into redshift bins of width Δz = 0.005, matching the native redshift grid used by LePhare and smaller than the typical scatter of our photometric redshifts, σ_z/(1 + z)≃0.011. In each bin, we selected the faintest 20% of galaxies in the r magnitude and computed their individual M_lim, i values using Eq. (3). From this subsample, we took the 95th percentile of the M_lim, i distribution. This value represents the mass above which 95% of the faintest galaxies would still be included in the sample and was used as the completeness threshold in that redshift bin. This process was repeated in all redshift bins, resulting in a set of completeness values that are fitted with a second-order polynomial. The fitted function M_lim(z) defines the redshift-dependent stellar mass limit, which we computed separately for the total, SF, and Q galaxy populations. The results are shown in Fig. 6. After applying the completeness cut, we excluded all galaxies with M < M_lim(z). For each remaining galaxy, its maximum redshift z_max, i was computed as the redshift at which its stellar mass equals M_lim(z). This was used in the estimation of volume-limited statistics such as V_max, following the procedure illustrated in Fig. 4 of Weigel et al. (2016).

Fig. 6. Stellar mass completeness limit Mlim(z) for the total, Q, and SF galaxy populations (orange, red, and blue dots, respectively). Second-order polynomial fits are shown for the total (dashed orange), Q (dash-dotted red), and SF (dotted blue) populations.

The final mass-complete sample consists of 302 393 SF and 376 171 Q galaxies. We note that the J-PLUS DR3 sample is complete at z ≤ 0.2 for stellar masses above log M_★/M_⊙ = 10.2 dex, hereafter log M_★.

3.2. The 1/V_max method

We used the 1/V_max technique from Schmidt (1968) to correct for Malmquist bias by weighting each galaxy by the maximum volume over which it could be detected. We defined V_max, i as the maximum volume in which a galaxy i at redshift z_i with stellar mass M_★, i could be detected. The number density Φ in mass bin j was computed by summing over all N_j objects in that bin:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_j=\frac{1}{dlog{M}_{★}}{\displaystyle \sum _i^{N_j}}\frac{1}{V_{\max ,i}}[{\mathrm{Mpc}}^{-3}{\mathrm{dex}}^{-1}],\end{array}$$(4)

where dlog M_★ ≈ (log M_★, max − log M_★, min)/n_bin, log M_★, max = 11.5 dex is the maximum considered mass, log M_★, min = 9.0 dex is the minimum considered mass, and n_bin = 10 is the number of bins. The comoving volume V_max, i in a flat universe is (Hogg 1999)

$$\begin{array}{c}\hfill V_{\max ,i}=\frac{4\pi }{3}\frac{{\mathrm{\Omega }}^s}{{\mathrm{\Omega }}^{\mathrm{sky}}}[d_c{({\mathrm{z}}_{\max ,i})}^3-d_{\mathrm{c}}{({\mathrm{z}}_{\min ,i})}^3],\end{array}$$(5)

with Ω^s = 2881 deg² being the effective surface area of J-PLUS DR3 after masking. Here, the surface area of the entire sky is Ω_sky = 41 253 deg², d_c (z) is the comoving distance at redshift z for the assumed cosmology, z_min, i = 0.05, and z_max, i is determined as

$$\begin{array}{c}\hfill z_{\max ,i}=min(z_{\max }^s,z_{\max ,i}^{\mathrm{mass}}).\end{array}$$(6)

We defined $z_{\max }^s=0.20$ as the upper redshift limit of our analysis. For each galaxy i, we computed its individual maximum redshift $z_{\max ,i}^{\mathrm{mass}}$, which corresponds to the redshift at which the galaxy would become undetected due to the stellar mass limit, M_lim(z), as described in Sect. 3.1. This was obtained by inverting the M_lim(z) relation to determine the highest redshift at which the galaxy stellar mass would still be above the limiting mass. This inversion is only meaningful within the mass and redshift ranges where M_lim(z) is defined and monotonic. As a result, galaxies whose stellar masses remain above M_lim(z) throughout the entire redshift range do not yield a finite solution from the inversion. For these cases, we assigned $z_{\max ,i}^{\mathrm{mass}}=z_{\max }^s=0.20$.

3.3. Error budget in the SMFs

Random uncertainties. We estimated the uncertainties in the number density of SMF, Φ_j, via bootstrap resampling. We generated 100 realizations by drawing, with replacement, a sample of the same size as the input sample. For each mass bin, Φ_j was taken as the median across realizations, and the 1σ uncertainty as the standard deviation of the Φ_j distribution.

Systematic uncertainties. To account for systematics between the Q and SF separation, we recomputed the SMFs using sSFR thresholds shifted by ±0.1 dex and took the systematic error in each bin as half of the absolute difference between the two perturbed cases. Additional systematic uncertainties related to stellar population modeling were not explicitly propagated in our SMFs, but are expected to be comparable to those discussed in earlier studies (e.g., Marchesini et al. 2009).

Cosmic variance. Cosmic variance (σ_v) is a significant source of uncertainty in surveys with small sky coverage, such as pencil-beam fields (≲1 deg²), where number-density fluctuations can be large owing to the limited cosmological volume probed (e.g., Díaz-García et al. 2024). In contrast, J-PLUS DR3 covers an effective area of 2881 deg², which should considerably mitigate the impact of cosmic variance on the measured SMF. We estimated σ_v using the framework of Moster et al. (2011), which models cosmic variance based on ΛCDM predictions and incorporates galaxy bias via a halo occupation distribution. From their Table 2, we adopted the root cosmic variance for dark matter, σ_dm = 0.181, obtained for the COSMOS survey, and rescaled it to the J-PLUS DR3 area using their Equation (12), correcting for the difference in survey area (from 1.96 deg² to 2881 deg²). Figure 7 shows that σ_v increases with stellar mass, as expected since massive galaxies are more strongly clustered and typically reside in dense environments. We find that for log(M_★/M_⊙) < 11, σ_v remains below 1%, reaching ∼1% only at log(M_★/M_⊙)∼11. Over the 9.8 < log(M_★/M_⊙) < 11.0 range, cosmic variance is slightly larger than the statistical uncertainties obtained via bootstrapping; nevertheless, it remains small in absolute terms. Because cosmic variance arises from large-scale density fluctuations, it is highly correlated across stellar-mass bins and, to first order, acts as a global normalization uncertainty on the SMF (i.e., primarily on ϕ_★), rather than as independent random noise in each bin (Smith 2012; López-Sanjuan et al. 2017; Díaz-García et al. 2024). We therefore do not add σ_v in quadrature to the per-bin error bars. Rather, we treated it as a small, global systematic uncertainty on the SMF normalization.

Fig. 7. Cosmic variance for this study, estimated with the Moster et al. (2011) prescription (purple circles), compared to the relative errors in the total SMF derived from bootstrapping (cyan squares).

To sum up, our error budget includes: (i) random uncertainties from finite sampling and the 1/V_max weighting captured by the bootstrap resampling; (ii) a systematic component associated with the adopted separation between the Q and SF galaxies, estimated by varying the sSFR threshold by ±0.1 dex; and (iii) an additional, stellar-mass correlated contribution from cosmic variance, σ_v, which primarily affects the overall normalization of the SMF. The final SMF measurements and their total statistical and systematic uncertainties are summarized in Appendix E.

3.4. Schechter function

The observed SMFs can be parameterized with a Schechter (1976) function, given by

$$\begin{array}{c}\hfill \mathrm{\Phi }(M)dM={\mathrm{\Phi }}^{\ast }{\left(\frac{M}{M^{\ast }}\right)}^{\alpha }exp(-\frac{M}{M^{\ast }})dM,\end{array}$$(7)

where M^* is the stellar mass at which the function transitions from a simple power law with slope α at lower masses into an exponential function at higher masses. The normalization Φ^* corresponds to the number density at M_★^*. In our case, for SMFs, it is better to work in log M_★ space and write the Schechter function as

$$\begin{array}{c}\hfill \mathrm{\Phi }dlogM=ln(10){\mathrm{\Phi }}^{\ast }\times e^{-{10}^{logM-log{M}_{★}^{\ast }}}\times {\left({10}^{logM-log{M}_{★}^{\ast }}\right)}^{\alpha +1}dlogM.\end{array}$$(8)

We obtained the Schechter functions from the SMF input stellar masses and stellar mass errors from SED fitting and 1/V_max from the Pozzetti et al. (2010) method. The Obreschkow et al. (2018) code fits the SMFs using the modified maximum-likelihood method, treating the Schechter function as a generative distribution and convolving it with the stellar-mass error distribution to correct for Eddington bias and selection effects. The likelihood was evaluated in log M_★ space using the stellar masses and 1/V_max weights described in Sect. 3.2. Parameter uncertainties were estimated by bootstrap, resampling the galaxy catalog 100 times. Although a double-Schechter form is often adopted for the local SMF (e.g., Baldry et al. 2012; Wright et al. 2018), it does not yield a statistically meaningful improvement over a single Schechter within our fitted mass range. The Bayesian evidence, evaluated over the mass-complete interval, consistently favors the single Schechter for the total, SF, and Q SMFs (Q-SMFs). We therefore adopted a single Schechter function throughout; for more details see Appendix D.

4. Results

4.1. Stellar mass function in J-PLUS DR3

The measured SMF for the total, SF, and Q populations in J-PLUS DR3 at 0.05 ≤ z_phot ≤ 0.2 are presented in Fig. 8a. We find that Q galaxies dominate at log M_★ >  10 dex, with a larger number density of SF galaxies at lower stellar masses.

Fig. 8. (a) Total, Q, and SF-SMFs from this study shown as orange circles, red squares, and blue triangles, respectively, with bootstrapping errors. Schechter fits are shown as dashed orange, dash-dotted red, and dotted blue lines. (b) Total SMF from J-PLUS DR3 (orange circles), compared to Peng et al. (2010) (gray pluses), Xu et al. (2025) (cyan hexagons), Baldry et al. (2012) (purple points), and Díaz-García et al. (2024) (dashed orange curve). GAMA SMFs from Wright et al. (2018) are shown as green, brown, and purple triangles with different orientations. (c) SF-SMF from J-PLUS DR3 (blue triangles), with the blue-shaded region showing total uncertainty (statistical plus systematic). This is compared to Baldry et al. (2012) (dash-dotted black curve), Peng et al. (2010) (cyan pluses), Xu et al. (2025) (light-blue hexagons), and Díaz-García et al. (2024) (dashed light-blue curve). (d) Q-SMF from J-PLUS DR3 (red squares), with the total uncertainty shown as a red-shaded region. This is compared to Baldry et al. (2012) (dash-dotted black curve), Peng et al. (2010) (dark red pluses), Xu et al. (2025) (light-red hexagons), and Díaz-García et al. (2024) (dashed dark-red curve).

The best-fit parameters of the Schechter function are presented in Table 1. A single Schechter function provides a proper description of the data in the stellar mass range sampled, as shown in Fig. 8a. We find that the characteristic mass is higher by 0.4 dex for Q galaxies, $log{M}_{★,\mathrm{Q}}^{\ast }=10.80$ dex, than for the SF population, $log{M}_{★,\mathrm{SF}}^{\ast }=10.43$ dex. This is accompanied by a higher characteristic density for Q galaxies. The faint-end slope is steeper for SF galaxies (α_SF = −1.2) than for Q galaxies (α_Q = −0.7). It is important to note that we also corrected the Eddington bias in the Schechter curve for the high masses of the SF-SMF, where our derived SF-SMF value lies slightly toward higher stellar masses. We observe small errors for all parameters and all populations, with values typically around 0.01. Given the small uncertainties on M_★ (∼0.01–0.02 dex), the ∼0.4 dex difference between the characteristic masses of Q and SF galaxies is significant and confirms that the two populations occupy distinct regions of the SMF. The difference in faint-end slope (Δα ≈ 0.5) is also significant, supporting the picture in which low-mass galaxies are predominantly star forming, while the Q population is increasingly dominated by massive galaxies. Our results imply that Q galaxies account for 45% of the number density in the local Universe at log M_★ >  9 dex, but 75% of the stellar mass density. A direct comparison between Schechter function fits and SMF must account for the strong dependence of the fit on the stellar mass range. Specifically, the faint-end slope, α, cannot be reliably determined in the absence of low-mass galaxies. Although we also fixed α with a modified code based on Obreschkow et al. (2018) to constrain the remaining parameters, the resulting parameters that yield the Schechter function did not satisfactorily match the SMF. For this reason, we ultimately allowed all three Schechter parameters to vary freely.

4.2. Fraction of the quiescent galaxies in J-PLUS DR3

The Q fraction, f_Q, was computed from the binned SMFs as

$$\begin{array}{c}\hfill f_{\mathrm{Q}}=\frac{{\mathrm{\Phi }}_{\mathrm{Q}}}{{\mathrm{\Phi }}_{\mathrm{Q}}+{\mathrm{\Phi }}_{\mathrm{SF}}}.\end{array}$$(9)

The Q fraction in J-PLUS DR3 is presented in Fig. 9. Note that f_Q directly increases with stellar mass. An increase of 40% in the Q fraction per dex is observed, reaching f_Q >  0.95 at log M_★ >  11 dex. This tells us that low-stellar-mass galaxies are dominated by SF galaxies while high-stellar-mass galaxies are dominated by Q ones. This is a clear indication of the existence of an important mass-related quenching mechanism. The errors in the red fraction as a function of stellar mass were computed from the uncertainties in Φ_Q and Φ_SF, and then propagated using Eq. (9). At the high-mass end, the small number of galaxies and the different mass-error distributions of Q and SF populations mean that residual Eddington bias and classification systematics might not be fully captured by our bootstrap errors. As discussed by Weaver et al. (2023), uncertainties in f_Q can be underestimated in this regime.

Fig. 9. Q fraction as a function of stellar mass measured from different surveys. The red dots represent this work. The blue points are from DESI Xu et al. (2025), the black crosses are from Peng et al. (2010), and the green curve is from Díaz-García et al. (2024).

5. Discussion

In this section, we interpret our SMF measurements from J-PLUS DR3 in the context of previous observational studies and theoretical models. We begin by comparing our total, Q, and SF-SMFs with literature results to assess consistency and systematic trends. We then examine the Q fraction as a tracer of galaxy evolution, evaluate the impact of cosmic variance, and compare our findings with predictions from the GAEA semi-analytic model (SAM).

5.1. Stellar mass function

In Fig. 8b, we compare the total SMF from J-PLUS DR3 with several measurements from the literature. The total SMF from the photometric study of Peng et al. (2010) corresponds to slightly higher space densities for all stellar mass bins by ∼0.1 dex. The GAMA survey results from Wright et al. (2018), shown for three redshift bins (0.02 < z < 0.08, 0.08 < z < 0.14, and 0.14 < z < 0.20), also estimate higher space densities with respect to our determination, with typical offsets between 0.05 and 0.15 dex, especially at the high-mass end. In contrast, the SMF from the DESI spectroscopic survey Xu et al. (2025) is in excellent agreement with our measurement, showing differences smaller than 0.05 dex over the full range. Finally, the total SMF from miniJPAS (Bonoli et al. 2021; Díaz-García et al. 2024) is very similar to our total SMF, but at log(M_★/M_⊙) > 11 it shows an excess of galaxies, reaching differences of up to ∼0.5 dex at the highest stellar masses. In Fig. 8c, we compare the star-forming stellar mass function (SF-SMF) from J-PLUS DR3 with previous measurements from Peng et al. (2010), Xu et al. (2025), Baldry et al. (2012), and Díaz-García et al. (2024). Our SF-SMF closely follows the shape and normalization of the SMF from Peng et al. (2010), with a small offset of approximately 0.1 dex. Compared to the DESI measurement from Xu et al. (2025), our SF-SMF lies slightly higher, particularly at log(M_★/M_⊙) > 10.5, where the difference can reach up to ∼0.15 dex. At lower stellar masses, both SMFs agree within the uncertainties. The SF-SMF from Díaz-García et al. (2024) also lies above ours, similarly to Peng et al. (2010) at log(M_★/M_⊙) < 10.5; at higher stellar masses, it predicts a larger number of galaxies. In Fig. 8d, we compare the Q-SMF with results from Peng et al. (2010) and Xu et al. (2025). Both literature SMFs lie slightly above our Q-SMF by up to 0.1 dex, but follow very similar trends. The agreement is particularly good with the DESI Q-SMF across most of the mass range. The Q-SMF from Díaz-García et al. (2024) is below ours, by 0.08 dex, and the rest of the literature for log(M_★/M_⊙) < 11 by 0.1 dex.

It is important to emphasize that once total uncertainties are considered, including statistical errors and systematics arising from the classification of Q and SF galaxies, our measurements remain fully consistent with previous studies. The shaded regions in Figs. 8c and d represent the total uncertainty envelopes, combining both components. These regions are similar in size to the scatter observed across the literature, especially at the low- and high-mass ends. This confirms the robustness of our methodology and supports the reliability of our population classification.

We conclude that the J-PLUS DR3 SMFs are in agreement with previous studies. The offsets between J-PLUS DR3 and previous SMFs (typically ≲0.1–0.15 dex) are comparable to the scatter among literature measurements and can be explained by differences in survey selection, redshift range, and stellar population modeling. Given its wide effective area (2881 deg²), intermediate depth, and 12-band optical photometry with seven narrowband filters, J-PLUS provides a particularly robust low-redshift SMF that we use as a reference for comparison with theoretical models and future surveys.

5.2. Quiescent fraction

In Fig. 9, we present the Q fraction (f_Q) as a function of stellar mass and compare it with two photometric studies: (Peng et al. 2010; Díaz-García et al. 2024) and the spectroscopic DESI results from Xu et al. (2025). The uncertainties in our Q fraction include both the statistical errors from bootstrapping the SMFs of Q and SF galaxies at log(sSFR) = − 10.2, and the systematic variations obtained by shifting this threshold by ±0.1 dex. Our measurement lies between (Peng et al. 2010) and Xu et al. (2025) curves and most closely follows the shape of the photometric result. Compared to Díaz-García et al. (2024), our (f_Q) is higher than theirs for z = 0.1. A notable difference arises in the DESI quenched fraction, which rises more steeply from log(M_★/M_⊙)∼9.5 and reaches a plateau by log(M_★/M_⊙)∼10.5. In contrast, J-PLUS and Peng et al. (2010) show a more gradual transition extending to log(M_★/M ⊙ ) ∼ 11.0. We quantified the agreement via point-by-point residuals after interpolating each curve to the J-PLUS mass grid, $\mathrm{\Delta }{f}_{\mathrm{Q}}(M_{★})\equiv f_{\mathrm{Q}}^{\text{J-PLUS}}-f_{\mathrm{Q}}^{\mathrm{lit}}$. For Xu et al. (2025) (DESI), we obtain RMSE = 0.086 and ⟨Δf_Q⟩= − 0.055; for Peng et al. (2010),RMSE = 0.061 and ⟨Δf_Q⟩= 0.022; and for Díaz-García et al. (2024),RMSE = 0.141 and ⟨Δf_Q⟩= 0.131. Overall, the rising f_Q trend with stellar mass observed in J-PLUS DR3 agrees with previous work and, being independent of SMF normalization, offers a robust probe of mass-dependent galaxy quenching.

5.3. Comparison with theoretical models: The GAEA simulation

The GAlaxy Evolution and Assembly (GAEA) model is a SAM of galaxy formation and evolution developed to run on top of dark matter halo merger trees from cosmological N-body simulations. Specifically, GAEA uses the Millennium Simulation (Springel et al. 2005), a large-volume ΛCDM simulation that traces the hierarchical growth of dark matter halos. The model includes detailed treatments of baryonic processes such as gas cooling, star formation, metal enrichment, stellar and AGN feedback, and environmental quenching. A key strength of GAEA lies in its chemically self-consistent enrichment scheme, which tracks individual elements with time-delayed feedback from type Ia and core-collapse supernovae (De Lucia et al. 2014), along with an updated treatment for stellar feedback (Hirschmann et al. 2016). Star formation is linked to the molecular hydrogen content of the interstellar medium through an H₂-based law (Xie et al. 2017). The most recent version, GAEA2023, introduces updated prescriptions for galaxy quenching and satellite evolution and has been calibrated to match observed SMFs and star formation histories across cosmic time (De Lucia et al. 2024; Fontanot et al. 2017). GAEA thus provides a comprehensive theoretical framework for interpreting statistical galaxy properties in the context of ΛCDM cosmology.

5.3.1. Quiescent and star-forming galaxies in GAEA

To distinguish between Q and SF galaxies in the GAEA simulation, we adopted the sSFR as a classification criterion. GAEA associates an SFR estimate with each model galaxy, and in some cases these values can be very small. To take observational effects into account, we applied the same procedure as in De Lucia et al. (2024) and assigned a new SFR value to each model galaxy with SFR < 10⁻⁴ using Eq. (10), which imposes an upper limit based on the empirical locus of Q galaxies in the SDSS:

$$\begin{array}{c}\hfill log\mathrm{SFR}=0.5\times log{M}_{★}-6.59,\end{array}$$(10)

where M_★ is in units of M_⊙. We applied this value to all GAEA galaxies with SFR < 10⁻⁴ and introduced a lognormal scatter of 0.25 dex to mimic observational uncertainties. In our J-PLUS DR3 data, the separation between Q and SF galaxies is defined by a threshold at log(sSFR) = − 10.2. In contrast, GAEA studies often adopt a more conservative cut at log(sSFR) = − 11.0 to define quenched systems. As shown in Fig. 10, both thresholds lie in the valley between the SF main sequence and the quenched population.

Fig. 10. Specific SFR vs. stellar mass for the GAEA simulation, shown as a 2D histogram. The main sequence of SF galaxies is visible at log(sSFR) >  − 11, while Q galaxies with log(sSFR) <  − 11 form a diagonal stripe, reproducing the locus of Q galaxies in the SDSS (see De Lucia et al. 2024). The dashed contours from our J-PLUS DR3 analysis are overplotted for comparison, as in Fig. 5.

5.3.2. Comparison of the stellar mass function

We constructed SMFs from the GAEA simulation using the modeled stellar masses and sSFR-based population division. Galaxies were binned by stellar mass without applying observational corrections. We selected galaxies at z = 0.11 and log(M_★/M_⊙) > 9. To mimic the J-PLUS selection, we applied an apparent magnitude limit of r < 20 mag and projected the simulation along the line of sight to account for observational geometry. Figure 11a presents the total SMF obtained from our GAEA sample, alongside the J-PLUS DR3 measurements and the GAEA2023 predictions from De Lucia et al. (2024). The two GAEA SMFs agree remarkably well across the intermediate- and high-stellar-mass ranges. The effect of the imposed magnitude limit on the predictions is also visible at low stellar masses. The reduced space densities reflect the loss of galaxies relative to the total GAEA SMF. At the low-mass end, our SMF falls below that of De Lucia et al. (2024) due to the imposed apparent magnitude cut, which excludes faint galaxies and reduces completeness. Compared to J-PLUS DR3, the GAEA SMF lies systematically above the observations at log(M_★/M_⊙) > 10.0, with an offset of up to ∼0.3 dex.

Fig. 11. Comparison of SMFs from J-PLUS DR3 and the GAEA simulation. (a) Total SMF from J-PLUS DR3 (orange circles), our GAEA-based sample using raw snapshot data and J-PLUS-like selection (purple squares), and the GAEA2023 model predictions published in De Lucia et al. (2024) (dotted green curve). We also plot the total SMF from Xu et al. (2025) (dash-dotted blue curve). (b) Q-SMFs derived using different thresholds in the range −11 < log(sSFRlim) <  − 10 (shaded orange-purple curves), compared with J-PLUS DR3 (red squares) and GAEA2023 from De Lucia et al. (2024) (dotted green curve). (c) SF-SMFs from J-PLUS DR3 (blue triangles), with GAEA predictions for varying log(sSFRlim) (shaded orange-purple curves). The dashed black line corresponds to log(sSFRlim) = − 10.4, the value adopted for consistency with J-PLUS DR3, while the solid black line shows the same selection after applying a −0.3 dex stellar mass shift to the SF-SMF.

To investigate the division between Q and SF galaxies, we computed SMFs using six thresholds in the range −11.0 < log(sSFR_lim) <  − 10.0 dex. Figures 11b and c show the resulting Q and SF-SMFs compared with the J-PLUS DR3 measurements and the GAEA2023 SMFs from De Lucia et al. (2024). We find that Q-SMFs are particularly sensitive to the sSFR threshold at lower stellar masses, while the high-mass regime remains more stable. As expected, adopting a more negative sSFR limit increases the number of galaxies classified as SF, leading to a higher SF-SMF and a corresponding decrease in the Q-SMF. This effect is most pronounced below log(M_★/M_⊙)∼10.5, where the population balance is more sensitive to classification. The best match to the J-PLUS DR3 SF-SMF occurs for thresholds near log(sSFR) = − 10.0, while the Q-SMF aligns better with log(sSFR) = − 10.4. The J-PLUS-derived separation value of log(sSFR) = − 10.4 dex thus represents a reasonable compromise for comparing populations in GAEA and J-PLUS. Nevertheless, the GAEA simulation predicts an excess of SF galaxies in the range 9.5 < log(M_★/M_⊙) < 11.0, which drives the total SMF above the observed values. Applying a stellar mass correction of −0.3 dex to the SF population in GAEA improves the agreement with J-PLUS DR3, as shown in Fig. 11c.

We find an offset between GAEA and J-PLUS DR3. While the SMF shapes agree well, the total SMF from GAEA2023 exceeds the observations by up to ∼0.3 dex at log(M_★/M_⊙)∼9.5–10.5. Such differences are within the uncertainties discussed in De Lucia et al. (2024), who report comparable offsets when comparing GAEA2023 to the SMFs from Muzzin et al. (2013) and Weaver et al. (2023).

5.3.3. Comparison of the quiescent fraction

Figure 12 presents the Q fraction as a function of stellar mass, comparing results from J-PLUS DR3 with predictions from the GAEA simulation. We display GAEA f_Q values for multiple sSFR thresholds, as well as the GAEA2023 f_Q values from De Lucia et al. (2024), which remain consistently below the observed f_Q over the full stellar mass range, reflecting the excess of SF galaxies in the simulation. This discrepancy is especially evident at intermediate and high masses. For the threshold log(sSFR_lim) = − 10.4 dex, consistent with the division adopted in our J-PLUS DR3 analysis, the GAEA predictions reproduce the observed Q fraction reasonably well at stellar masses below log(M_★/M_⊙)∼10.5. At higher stellar masses, however, GAEA underpredicts the Q fraction. This indicates that a larger fraction of massive galaxies remain SF in the simulation. This discrepancy persists even after applying a −0.3 dex correction to the stellar masses of SF galaxies. Although the adopted threshold provides a reasonable basis for comparison, the results suggest that quenching is not fully efficient in GAEA at the high-mass end.

Fig. 12. fQ from J-PLUS DR3 (red dots) compared to GAEA predictions for different log(sSFRlim) thresholds (color-coded lines). We also show the GAEA2023 results from De Lucia et al. (2024) (dotted green line). The dashed black line shows the fQ for log(sSFRlim) = − 10.4, and the solid black line applies a −0.3 dex stellar mass correction to the GAEA SF population.

6. Conclusions

We presented the SMF of Q and SF galaxies at 0.05 ≤ z ≤ 0.2 based on high-quality photometric data from J-PLUS DR3. Using a sample of 890 844 galaxies selected with r₀ ≤ 20 mag and 0.05 < z < 0.20, we derived stellar masses and SFRs via SED fitting with CIGALE and separated galaxy populations based on their specific SFRs. The main conclusions of our study are as follows.

Stellar mass functions. We computed the SMFs for total, Q, and SF galaxies using the 1/V_max method, accounting for stellar mass completeness and Eddington bias. All three SMFs are well described by single Schechter functions and are consistent with the literature. Quiescent galaxies dominate at log M_★ >  10 dex, while SF galaxies are more numerous at lower stellar masses. Quiescent galaxies account for 45% of the number density in the local Universe at log M_★ >  9 dex, but 75% of the stellar mass density.

Quiescent fraction. The fraction of Q galaxies increases by 40% per dex in stellar mass, as expected from mass-dependent quenching, reaching f_Q >  0.95 at log M_★ >  11 dex. Our measured Q fraction lies between the values obtained in previous studies, reinforcing the robustness of our classification.

Cosmic variance. Given the large effective area covered by J-PLUS DR3 (2881 deg²), the impact of cosmic variance on our SMF measurements is minor. Using the prescription of Moster et al. (2011), we estimate relative uncertainties due to cosmic variance to be below 1% for log M_★ < 11.5 dex.

Comparison with simulations. We compared our observed SMFs with predictions from the GAEA SAM. GAEA2023 overpredicts the abundance of SF galaxies at 9.5 < log M_★ < 11.0 dex, leading to systematic excess in the total SMF and a corresponding deficit in the Q population. The offset reaches up to ∼0.3 dex and persists despite matching sSFR thresholds and applying mass corrections. These differences are consistent with those reported in De Lucia et al. (2024) and remain within the expected systematic uncertainties.

These results validate the scientific potential of J-PLUS DR3 for low-redshift galaxy evolution studies. The methods and findings presented here can be extended with J-PAS, which will provide deeper and higher-resolution photometry over a wider spectral range. The precision of the derived SMFs and Q fractions, combined with the large survey area, enables statistically robust analyses of environmental effects in future work.

Data availability

The catalog associated with this article, described in Table F.1 is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/709/A157

Acknowledgments

F. D. A. B. acknowledge the funding from Erasmus IBERUS+ Student Mobility for Traineeships, and Programa Becas Ibercaja-CAI Estancias de Investigación for the research stay at INAF to work with G.d.L. and M.F. in the analysis of the GAEA simulations. The authors thank the interesting suggestions of the referee, which improved the paper. J. A. F. O., H. D. S., and A. E. acknowledge the financial support from the Spanish Ministry of Science and Innovation and the European Union – NextGenerationEU through the Recovery and Resilience Facility (RRF) project ICTS-MRR-2021-03-CEFCA. H. D. S. also acknowledges financial support by RyC2022-030469-I grant, funded by MCIU/AEI/10.13039/501100011033 and FSE+. A. L. C. and P. T. R. acknowledge the financial support from the European Union – NextGenerationEU through the RRF program Planes Complementarios con las CCAA de Astrofísica y Física de Altas Energías – LA4. J. V. M. acknowledges financial support by PID2022-136598NB-C32 grant. L. A. D. G. acknowledges financial support from the State Agency for Research of the Spanish MCIU through ’Center of Excellence Severo Ochoa’ award to the Instituto de Astrofísica de Andalucía (CEX2021-001131-S) funded by MCIN/AEI/10.13039/501100011033 and to PID2022-141755NB-I00. Based on observations made with the JAST80 telescope and T80Cam camera for the J-PLUS project at the Observatorio Astrofísico de Javalambre (OAJ), in Teruel, owned, managed, and operated by the Centro de Estudios de Física del Cosmos de Aragón (CEFCA). We acknowledge the OAJ Data Processing and Archiving Unit (UPAD; Cristóbal-Hornillos et al. 2012) for reducing the OAJ data used in this work. Funding for the J-PLUS Project has been provided by the Governments of Spain and Aragón through the Fondo de Inversiones de Teruel; the Aragonese Government through the Research Groups E96, E103, E16_17R, E16_20R, and E16_23R; the Spanish Ministry of Science and Innovation (MCIN/AEI/10.13039/501100011033 y FEDER, Una manera de hacer Europa) with grants PID2021-124918NB-C41, PID2021-124918NB-C42, PID2021-124918NA-C43, and PID2021-124918NB-C44; the Spanish Ministry of Science, Innovation and Universities (MCIU/AEI/FEDER, UE) with grants PGC2018-097585-B-C21 and PGC2018-097585-B-C22; the Spanish Ministry of Economy and Competitiveness (MINECO) under AYA2015-66211-C2-1-P, AYA2015-66211-C2-2, AYA2012-30789, and ICTS-2009-14; and European FEDER funding (FCDD10-4E-867, FCDD13-4E-2685). The Brazilian agencies FINEP, FAPESP, and the National Observatory of Brazil have also contributed to this Project.

References

Appendix A: ADQL of the data

The ADQL query used in this study is:

SELECT 
  obj.TILE_ID,
  obj.number,
  obj.ALPHA_J2000 AS RA,
  obj.DELTA_J2000 AS DEC,
  obj.mask_flags,
  obj.MAG_AUTO AS MAG_AUTO,
  obj.MAG_ERR\_AUTO,
  obj.class_star AS CLASS_STAR,
  lephare.photoz,
  lephare.odds AS ODDS,
  galclass.TILE_ID AS galclass_TILE_ID,
  galclass.number AS galclass_number,
  galclass.sglc_prob_star,
  mwe.TILE_ID AS mwe_TILE_ID,
  mwe.number AS mwe_number,
  mwe.ax AS ax,
  mwe.ax_err AS ax_err

FROM jplus.MagABDualObj AS obj

JOIN jplus.StarGalClass AS galclass
  ON obj.TILE_ID = galclass.TILE_ID
  AND obj.number = galclass.number

JOIN jplus.PhotoZLephare AS lephare
  ON obj.TILE_ID = lephare.TILE_ID
  AND obj.number = lephare.number

JOIN jplus.MWExtinction AS mwe
  ON obj.TILE_ID = mwe.TILE_ID
  AND obj.number = mwe.number

WHERE 
  lephare.photoz >= 0.05 
  AND lephare.photoz < = 0.2025
  AND obj.mask_flags[jplus::rSDSS] = 0
  AND galclass.sglc_prob_star < 0.5
  AND (obj.MAG_AUTO[jplus::rSDSS] - mwe.ax[
        jplus::rSDSS]) <= 20

Appendix B: CIGALE configuration file

Below we summarize the modules and grid of parameters used in the CIGALE configuration file.

#Configuration of the SED creation modules.
[sed_modules_params]
  
  [[sfhdelayed]]
  
    tau_main = 100., 1000., 5000.
    
    age_main = 3000., 5000., 7000., 10000., 
        13000.
    
    tau_burst = 5., 100., 500.

    age_burst = 5.,10.,50., 100., 500., 1000.

    f_burst =  0.001, 0.01, 0.2, 0.4

    sfr_A = 1.0

    normalise = True
 
  
  [[bc03]]
  
    imf = 1

    metallicity =  0.004, 0.008, 0.02

    separation_age = 10
       
  
  [[nebular]]

    logU = -2.0, -3.5

    zgas = 0.014, 0.022

    ne = 100

    f_esc = 0.0, 0.2

    f_dust = 0.0

    lines_width = 300.0

    emission = True
    
  
  [[dustatt_modified_starburst]]

    E_BV_lines = 0.0, 0.1, 0.2, 0.3, 0.5

    E_BV_factor = 0.44

    uv_bump_width = 35.0
 
    uv_bump_amplitude = 0.0

    powerlaw_slope = 0.0

    Ext_law_emission_lines = 1

    Rv = 3.1

Appendix C: Comparison between our sample and the Duarte Puertas et al. (2017) sample

In the first panel of Fig. C.1, the sample of Duarte Puertas et al. (2017) is centered on Hα, and is therefore more focused on blue galaxies. Their sample is also complete down to r_SDSS < 17.7 mag, corresponding to the SDSS main galaxy sample, whereas our sample was selected with r < 20 mag. In the second panel of Fig. C.1, we compare the redshift distributions of the cross-matched sample between our study and Duarte Puertas et al. (2017). The cross-matched sample contains more galaxies at lower redshift, whereas our full sample includes more galaxies at higher redshift, as reflected by the difference in the median values. In the third panel of Fig. C.1, the distribution of the rest-frame color (g − i) is shown. We find that the cross-matched sample is more concentrated toward blue galaxies, whereas our sample shows a clearer bimodality.

Fig. C.1. Distributions of the crossmatch sample (solid blue) and the full sample of our study (dashed red). Left panel: Distribution in r-band magnitude. Middle panel: Distribution of the redshift. Right panel: Distribution of the rest-frame color (g − i)0.

Appendix D: Single versus double Schechter model comparison

We fit Schechter functions with the Obreschkow et al. (2018) modified maximum-likelihood (MML) method implemented in dftools::dffit. In this approach the Schechter function is treated as a generative distribution in log M_★ and is convolved with the stellar-mass error distribution; selection effects are incorporated through the 1/V_max weights described in Sect. 3.2. The single-Schechter function (SSF) has three free parameters (log₁₀ϕ₁,  log₁₀M^★,  α₁), while the double-Schechter function (DSF) adopts a shared characteristic mass M^★ and adds a second component (log₁₀ϕ₂,  α₂), for a total of five free parameters. For completeness, the best-fitting DSF parameters used in this comparison are reported in Table D.1, and a direct visual comparison between the best-fitting SSF and DSF is shown in Fig. D.1.

Fig. D.1. Comparison of best-fitting single- (solid) and double-Schechter (dashed) models over the observed SMFs (points with error bars). The curves show the Eddington-convolved (observed-space) model predictions from the Obreschkow et al. (2018) fitting framework. The vertical dotted lines mark the stellar-mass completeness limits at z = 0.2 (from Fig. 6); bins below these limits are shown as gray open symbols, and the model curves are faded below the completeness limit. For the total and SF-SMFs the double-Schechter curve corresponds to the standard (free) five-parameter optimization, while for the Q-SMF the double-Schechter curve corresponds to the stabilized two-stage (tight) procedure (Stage B) described in the text.

For each sample (total, SF, and Q) we obtain the best-fitting parameters with the MML method and compare the SSF and DSF through the Bayesian evidence returned by the fitting framework (Kass & Raftery 1995; Trotta 2008),

$$\begin{array}{c}\hfill Z\equiv p(D\mid M)={\displaystyle \int p}(D\mid \theta ,M)p(\theta \mid M)d\theta ,\end{array}$$(D.1)

where D denotes the data, M the model, θ the model parameters, p(D ∣ θ, M) the likelihood, and p(θ ∣ M) the prior. The Bayesian evidence thus quantifies the probability of the data under a given model after marginalizing over its parameter space, naturally accounting for both goodness of fit and model complexity. To ensure consistency with the interpretation of Fig. D.1, this comparison is performed over the mass-complete interval for each population. We report the difference in log-evidence as ΔlnZ ≡ lnZ_DSF − lnZ_SSF, where negative values indicate that the SSF is preferred over the DSF. For the total, SF, and Q-SMFs, we obtain ΔlnZ = −2611, −574, and −23, respectively. In all three cases, the Bayesian evidence favors the SSF. As shown in Fig. D.1, the best-fitting SSF and DSF curves are nearly indistinguishable over the fitted mass range.

The MML approach improves on a simple binned estimator by fitting the unbinned stellar-mass distribution directly, while incorporating the redshift-dependent selection function and the convolution with stellar-mass uncertainties, thus correcting for Eddington bias without introducing an arbitrary dependence on binning. In the present analysis, the SSF–DSF comparison is constrained mainly by the mass-complete regime at the upper redshift boundary, that is, above M_lim(z = 0.2) for each population. The lower-mass bins are included in Fig. D.1 to illustrate the observed behavior of the SMF below the mass-complete regime, but they are not uniformly complete across 0.05 ≤ z ≤ 0.2 and therefore do not constrain the SSF–DSF comparison.

The DSF fits used for Fig. D.1 are obtained as follows. For the total and SF-SMFs we use the standard single-stage DSF optimization (free fit), in which all five parameters are varied simultaneously within the MML framework. For the Q-SMF, the faint component is comparatively weakly constrained over the present mass-complete interval and exhibits strong degeneracies between (ϕ₂, α₂, M^★) when all parameters are left free. To ensure a well-posed DSF solution for Q galaxies, we therefore adopt a stabilized two-stage (tight) procedure: the bright end is first anchored with a high-mass prefit (Stage A), and the full mass range is then refitted (Stage B) using a pivot reparameterization of the faint component together with weakly informative priors and an ordering constraint (α₁ >  α₂ + 0.3). This strategy stabilizes the faint component without altering the bright-end description. For the total and Q-SMFs, however, the optimization drives the second component to a negligible normalization (Table D.1), such that the resulting DSF is effectively indistinguishable from the SSF over the fitted mass range. This behavior is illustrated explicitly for the total SMF in Fig. D.2, where the second component remains well below the dominant component and the summed DSF closely follows the SSF. We therefore interpret these cases as indicating that the second component is not meaningfully constrained by the present data, rather than as evidence against a DSF description of the local SMF in general.

Fig. D.2. Double-Schechter decomposition for the total SMF. The observed SMF is shown with points and error bars. The two double-Schechter components are plotted separately (component 1: dashed; component 2: dotted), together with their sum (dash-dotted) and the best-fitting single-Schechter model (solid). The vertical dotted line marks the stellar-mass completeness limit at z = 0.2; bins below this limit are shown as gray open symbols, and model curves are faded below the completeness limit.

Finally, we note that many low-redshift studies find that a DSF provides an improved description of the total SMF when extending to substantially lower stellar masses and/or when a wider dynamic range is available. In our case, the SSF–DSF comparison is effectively anchored to the mass-complete regime at the upper redshift boundary, log M_★ ≥ M_lim(z = 0.2) (Fig. 6), with M_lim(0.2) = 10.2 for the total sample, 10.1 for SF galaxies, and 10.4 for Q galaxies (vertical lines in Fig. D.1). The low-mass bins displayed in Fig. D.1 are not complete across 0.05 ≤ z ≤ 0.2. Within our adopted forward-modeling framework, which includes convolution with the stellar-mass uncertainties, the data in the mass-complete range do not statistically require an additional Schechter component, and the SSF is preferred.

Appendix E: Stellar mass function values

Table E.1 lists the measured binned SMFs for the total, SF, and Q galaxy samples.

Appendix F: Main headers of the catalog

Table F.1 summarizes the main columns of the released catalog together with their units and descriptions.

All Tables

All Figures

	Fig. 1. Sky location of J-PLUS DR3 galaxies with r ≤ 20 mag at 0.05 ≤ zphot ≤ 0.20 covering 2881 deg2 (blue points). The black curve depicts the Milky Way.
In the text

Fig. 2. Examples of SED fitting to J-PLUS 12-band photometry using CIGALE. Panel (a): Q galaxy. Panel (b): SF galaxy. In each panel, the observed fluxes are shown as purple circles with error bars, and the best-fit model spectrum is plotted as a solid black line. The red dots indicate the model-predicted fluxes in each filter. The contributions from the attenuated stellar population (orange line), unattenuated stellar emission (dashed blue line), and nebular emission (green line) are also shown, although difficult to see in the left panel due to their resemblance to the model spectrum. Bottom subpanels: Relative residuals, defined as (Obs − Mod)/Obs, with a horizontal dashed line at zero for reference. A gri composite image from J-PLUS is inset in the bottom right corner of each panel.

In the text

	Fig. 3. Comparison between SFR (a), stellar mass (b), and sSFR (c) from this study and Duarte Puertas et al. (2017) for the N common sources. The red lines show the best Gaussian fits with median μ and dispersion σ.
In the text

	Fig. 4. Comparison of stellar masses from this study and Taylor et al. (2011) for the N common sources. The red line shows the best Gaussian fit with median μ and dispersion σ.
In the text

	Fig. 5. sSFR vs. stellar mass 2D odds-weighted histogram. The two maxima correspond to SF (higher sSFR values) and Q (lower sSFR values) populations. The limit between both populations is set to log sSFR = −10.2 dex (horizontal dashed line).
In the text

	Fig. 6. Stellar mass completeness limit Mlim(z) for the total, Q, and SF galaxy populations (orange, red, and blue dots, respectively). Second-order polynomial fits are shown for the total (dashed orange), Q (dash-dotted red), and SF (dotted blue) populations.
In the text

	Fig. 7. Cosmic variance for this study, estimated with the Moster et al. (2011) prescription (purple circles), compared to the relative errors in the total SMF derived from bootstrapping (cyan squares).
In the text

Fig. 8. (a) Total, Q, and SF-SMFs from this study shown as orange circles, red squares, and blue triangles, respectively, with bootstrapping errors. Schechter fits are shown as dashed orange, dash-dotted red, and dotted blue lines. (b) Total SMF from J-PLUS DR3 (orange circles), compared to Peng et al. (2010) (gray pluses), Xu et al. (2025) (cyan hexagons), Baldry et al. (2012) (purple points), and Díaz-García et al. (2024) (dashed orange curve). GAMA SMFs from Wright et al. (2018) are shown as green, brown, and purple triangles with different orientations. (c) SF-SMF from J-PLUS DR3 (blue triangles), with the blue-shaded region showing total uncertainty (statistical plus systematic). This is compared to Baldry et al. (2012) (dash-dotted black curve), Peng et al. (2010) (cyan pluses), Xu et al. (2025) (light-blue hexagons), and Díaz-García et al. (2024) (dashed light-blue curve). (d) Q-SMF from J-PLUS DR3 (red squares), with the total uncertainty shown as a red-shaded region. This is compared to Baldry et al. (2012) (dash-dotted black curve), Peng et al. (2010) (dark red pluses), Xu et al. (2025) (light-red hexagons), and Díaz-García et al. (2024) (dashed dark-red curve).

In the text

	Fig. 9. Q fraction as a function of stellar mass measured from different surveys. The red dots represent this work. The blue points are from DESI Xu et al. (2025), the black crosses are from Peng et al. (2010), and the green curve is from Díaz-García et al. (2024).
In the text

	Fig. 10. Specific SFR vs. stellar mass for the GAEA simulation, shown as a 2D histogram. The main sequence of SF galaxies is visible at log(sSFR) >  − 11, while Q galaxies with log(sSFR) <  − 11 form a diagonal stripe, reproducing the locus of Q galaxies in the SDSS (see De Lucia et al. 2024). The dashed contours from our J-PLUS DR3 analysis are overplotted for comparison, as in Fig. 5.
In the text

Fig. 11. Comparison of SMFs from J-PLUS DR3 and the GAEA simulation. (a) Total SMF from J-PLUS DR3 (orange circles), our GAEA-based sample using raw snapshot data and J-PLUS-like selection (purple squares), and the GAEA2023 model predictions published in De Lucia et al. (2024) (dotted green curve). We also plot the total SMF from Xu et al. (2025) (dash-dotted blue curve). (b) Q-SMFs derived using different thresholds in the range −11 < log(sSFRlim) <  − 10 (shaded orange-purple curves), compared with J-PLUS DR3 (red squares) and GAEA2023 from De Lucia et al. (2024) (dotted green curve). (c) SF-SMFs from J-PLUS DR3 (blue triangles), with GAEA predictions for varying log(sSFRlim) (shaded orange-purple curves). The dashed black line corresponds to log(sSFRlim) = − 10.4, the value adopted for consistency with J-PLUS DR3, while the solid black line shows the same selection after applying a −0.3 dex stellar mass shift to the SF-SMF.

In the text

	Fig. 12. fQ from J-PLUS DR3 (red dots) compared to GAEA predictions for different log(sSFRlim) thresholds (color-coded lines). We also show the GAEA2023 results from De Lucia et al. (2024) (dotted green line). The dashed black line shows the fQ for log(sSFRlim) = − 10.4, and the solid black line applies a −0.3 dex stellar mass correction to the GAEA SF population.
In the text

	Fig. C.1. Distributions of the crossmatch sample (solid blue) and the full sample of our study (dashed red). Left panel: Distribution in r-band magnitude. Middle panel: Distribution of the redshift. Right panel: Distribution of the rest-frame color (g − i)0.
In the text

Fig. D.1. Comparison of best-fitting single- (solid) and double-Schechter (dashed) models over the observed SMFs (points with error bars). The curves show the Eddington-convolved (observed-space) model predictions from the Obreschkow et al. (2018) fitting framework. The vertical dotted lines mark the stellar-mass completeness limits at z = 0.2 (from Fig. 6); bins below these limits are shown as gray open symbols, and the model curves are faded below the completeness limit. For the total and SF-SMFs the double-Schechter curve corresponds to the standard (free) five-parameter optimization, while for the Q-SMF the double-Schechter curve corresponds to the stabilized two-stage (tight) procedure (Stage B) described in the text.

In the text

Fig. D.2. Double-Schechter decomposition for the total SMF. The observed SMF is shown with points and error bars. The two double-Schechter components are plotted separately (component 1: dashed; component 2: dotted), together with their sum (dash-dotted) and the best-fitting single-Schechter model (solid). The vertical dotted line marks the stellar-mass completeness limit at z = 0.2; bins below this limit are shown as gray open symbols, and model curves are faded below the completeness limit.

In the text