© The Authors 2025

1. Introduction

The bulge in the central region of a disk galaxy is a crucial component for understanding the formation and evolution of disk galaxies. Recent studies have shown the diverse nature of bulges, which can be classified into classical bulges and pseudobulges. Classical bulges, often described as mini-ellipticals at the centers of disks, are believed to form early through violent merger processes, offering vital insights into the early stages of galaxy formation (Eggen et al. 1962; Toomre & Toomre 1972; Gott 1977; Bender et al. 1992; Renzini 1999). In contrast, pseudobulges exhibit disky shapes, relatively young stellar populations, flat luminosity distributions, and rotationally dominated kinematics (Kormendy & Kennicutt 2004; Gadotti 2009; Kormendy 2013). These features suggest that pseudobulges form over long timescales via secular processes, where non-axisymmetric structures such as bars and spiral arms drive gas inflow toward the galaxy center, fueling central star formation (Sellwood & Wilkinson 1993; Kormendy & Kennicutt 2004; Sellwood 2014; Tonini et al. 2016). Distinguishing between classical and pseudobulges is crucial for a deeper understanding of disk galaxy evolution.

In the local Universe, a significant fraction of disk galaxies (∼60%) host bars in their central regions (Eskridge et al. 2000; Menéndez-Delmestre et al. 2007; Erwin 2017), which can influence the measurement of bulge parameters in 2D image decomposition (Laurikainen et al. 2004; Gadotti 2008; Fisher & Drory 2008). Both observations (Lütticke et al. 2000; Laurikainen et al. 2014; Erwin & Debattista 2017; Blaña Díaz et al. 2017; Gadotti et al. 2020) and simulations (Raha et al. 1991; Berentzen et al. 1998; Martinez-Valpuesta & Shlosman 2004; Li & Shen 2015; Fragkoudi et al. 2017) have suggested that the bar consists of two components: an elongated thin bar and a thickened inner structure. In edge-on views, this inner component appears boxy, peanut-shaped or X-shaped, and it is referred to as the boxy/peanut (B/P) bulge, a particular type of pseudobulge. In face-on views, the B/P bulge appears as a lens-like oval structure known as the barlens, which typically extends to about half the length of the bar and has a lower ellipticity (Laurikainen et al. 2011; Athanassoula et al. 2015; Buta et al. 2015; Laurikainen & Salo 2017). Recent studies have shown that the fraction of B/P bulges in barred galaxies in the local Universe ranges from 50–70% (Erwin & Debattista 2017; Kruk et al. 2019; Gadotti et al. 2020). Li et al. (2017) and Erwin & Debattista (2017) further demonstrated that the presence of a B/P bulge is strongly correlated with galaxy stellar mass: The majority (∼80%) of barred galaxies with log(M_⋆/M_⊙) ≥ 10.4 host B/P bulges, compared to only an ∼20% B/P bulge fraction among lower-mass galaxies. The formation mechanisms of the B/P bulge could include the buckling instability of the bar (Raha et al. 1991; Merritt & Sellwood 1994; Martinez-Valpuesta & Shlosman 2004; Smirnov & Sotnikova 2019) or the vertical resonances of the stellar orbits within the bar (Combes & Sanders 1981; Combes et al. 1990; Quillen et al. 2014; Sellwood & Gerhard 2020).

The barlens, which resides in the central region of a galaxy, can easily be mistaken as part of the bulge in face-on views or at moderate inclinations. Previous studies have typically consider three components in disk galaxy image decomposition: the bulge, bar, and disk. Ignoring the barlens component would result in attributing its flux to the bulge, leading to an overestimation of the bulge-to-total (B/T) ratio and Sérsic index (n) (Laurikainen et al. 2018; Erwin et al. 2021). Athanassoula et al. (2015) performed 2D image decomposition on N-body simulations of barred disks to derive the parameters of the barlens. Using generalized elliptical isophotes for all components, they found that the Sérsic index of the barlens is less than two, and it accounts for 0.13–0.33 of the galaxy’s total flux. The barlens ellipticity ranges from 0.7 to 1.0, and its length is 0.4–0.8 times the length of the bar, with only a weak dependence on inclination angle. In a 2D image decomposition of galaxies from the Calar Alto Legacy Integral Field Area (CALIFA) survey, Laurikainen et al. (2018) showed that the barlens contributes approximately 13% of the total galaxy light, with the remaining flux primarily attributed to a possible distinct bulge component that contributes ≤10%. They also found that the barlens typically exhibits an exponential profile (n ≤ 1).

In previous work, the influence of modeling different galaxy structures, including the bar, ring, spiral arms and lens, on the bulge parameters has been extensively studied (Gadotti 2008; Gao & Ho 2017; Chugunov et al. 2024; Marchuk et al. 2024). Athanassoula et al. (2015) have proposed methods to measure the barlens component in 2D image decomposition. However, without gradually increasing the modeling complexity, their modeling strategy could not properly answer this important question, namely the impact of modeling the barlens on the bulge parameters. Built upon the previous efforts in the barlens modeling, here we aim to develop an optimal strategy for modeling the barlens in 2D image decomposition and to quantify its influence on the bulge measurements. We carefully designed a series of decomposition tests that minimize parameter degeneracies between the overlapping structures (i.e., bar, bulge, and barlens). To better constrain the intrinsic parameters of the barlens, we first used mock images from two different simulations: a pure-disk N-body simulation featuring a buckled bar without a preexisting classical bulge and a Milky Way-like galaxy from the Auriga simulation with the bulge built in the previous merger process. Afterward, we applied the same 2D image decomposition strategy to real observational images of two nearby barred galaxies that exhibit the barlens structure. Unlike previous works (e.g., Athanassoula et al. 2015), which primarily emphasized the identification and characterization of the barlens component itself, our primary focus here is to study the influence of the barlens modeling on the bulge parameters. We gradually increased the complexity of the model to better understand the variation of the bulge parameters after the barlens component was included. Moreover, we also investigated the effect of the disk inclination on the structural decomposition results, a factor that has not been thoroughly explored in previous barlens studies.

This paper is organized as follows: In Section 2, we provide an overview of the mock images, observational data, and the decomposition method used to model the barlens component. The results of the mock image decomposition are presented in Section 3, along with applications to real observations. The effects of the inclination angle and other statistical tests on the decomposition results are discussed in Section 4. A summary of our findings is given in Section 5.

2. Data and method

2.1. Mock image creation from simulation

To first understand the intrinsic properties of the barlens, we began with a mock galaxy image generated from a pure-disk galaxy simulation without any preexisting classical bulge (hereafter Model 1). The initial conditions (ICs) were set up following Tepper-Garcia et al. (2021), with a dark matter halo of mass M_h ∼ 1.2 × 10¹² M_⊙, a stellar disk of mass M_d = 4.3 × 10¹⁰ M_⊙, a scale length of 2.5 kpc, and a scale height of 0.3 kpc (resembling a thin disk). Both the dark matter halo and the disk were modeled with 10⁶ particles. The ICs were generated using AGAMA (Vasiliev 2019) and evolved with GADGET4 (Springel et al. 2021) for 5 Gyr. During the evolution of Model 1, a bar formed at ∼1.1 Gyr and buckled into a peanut-shaped bulge at ∼2.3 Gyr.

Next, to evaluate the impact of the barlens on the bulge measurement, we used a simulated galaxy from the Auriga project¹ (Grand et al. 2017, 2024). The Auriga project includes 30 Milky Way-like galaxies simulated within the framework of ΛCDM cosmology using the moving mesh code AREPO (Springel 2010). The simulations account for various physical processes, such as star formation, chemical enrichment, and feedback effects. Further details can be found in Grand et al. (2017). The Auriga project successfully forms a series of barred galaxies. After visually checking the Auriga simulations, we decided to use the galaxy Auriga 6 at z = 0 (hereafter Model 2) to generate mock images, which has a clear bar with a prominent B/P bulge from the edge-on view. The galaxy also features a ring structure surrounding the bar and weak spiral arms in its outer disk.

Similar to previous works (Erwin & Debattista 2013; Athanassoula et al. 2015; Anderson et al. 2022), we created mock images of the simulations based on the stellar mass distribution. The particles were projected onto a 2D grid of pixels, with each pixel value representing the sum of the masses of all the stellar particles within that pixel. To simulate observations from ground-based telescopes, such as the Carnegie-Irvine Galaxy Survey (CGS, Ho et al. 2011; Li et al. 2011), we assumed that the galaxy is located at a distance of 40 Mpc. The simulations were projected in the face-on view (in the X-Y plane), with a pixel scale of approximately 0.26 arcsec pixel⁻¹, or 0.048 kpc pixel⁻¹. The image was then convolved with the point spread function (PSF) from the R-band CGS images². The total flux of the convolved image was then scaled to match the total flux of the R-band images of nearby CGS disk galaxies, with the total apparent magnitude of ∼11.6 mag. Poisson noise and the sky background were then added to the image, with the sky background modeled as a Gaussian distribution, where the mean and standard deviation of the sky pixels were determined based on the sky values of the R-band CGS images.

The mock images of the two simulations in face-on, side-on, and end-on views are shown in Fig. 1. In the side-on view, the B/P bulge is visible, corresponding to the barlens in the face-on view. To fully explore how the disk inclination angle affects the measurements of the barlens and the bulge parameters, we generated mock images for different disk inclinations (i = 30°, 45°, 60°) and various position angles of the bar (PA_bar = 0°, 30°, 60°, and 90°).

Fig. 1. Mock images of Model 1 and Model 2 in the face-on view (top), side-on view (middle), and end-on view (bottom). The barlens and the corresponding B/P bulge can be seen clearly in the face-on view and edge-on views, respectively. The blue ellipses in the top panels show the visually identified outer edge of the barlenses.

2.2. Observational images from the CGS survey

After testing the barlens decomposition method with mock images, we applied the same method to two nearby barred galaxies (NGC 1533 and NGC 7329) in the CGS survey that exhibit clear barlens structures (Li et al. 2017). The R-band images of these galaxies are shown in Fig. 2. The images have a pixel scale of 0.26 arcsec pixel⁻¹, with a median seeing of 1.01 arcsec and a surface brightness depth of 26.4 mag arcsec⁻². In our analysis, we used star-cleaned images and PSF images from Li et al. (2011). The star-cleaned images are free from contamination by foreground stars and background galaxies. In the following 2D image decomposition, an additional mask image is not needed. Comprehensive 2D image decompositions, including the bulge, bar, and ring components, have been carried out for these galaxies by Gao & Ho (2017). However, their work did not take into account the barlens as an additional component, which could potentially influence the bulge properties.

Fig. 2. R-band star-cleaned images of NGC 1533 (left) and NGC 7329 (right) with the barlens component identified in Li et al. (2017). The blue ellipses show the visually identified outer edge of the barlenses.

2.3. Image decomposition strategy

2.3.1. Light profile functions and isophotal shapes

We used GALFIT to perform the 2D image decomposition (Peng et al. 2010). GALFIT provides a variety of analytic functions for modeling light profiles, including Sérsic (Sérsic 1963), exponential, modified Ferrer, Moffat, and King functions. Additionally, the isophotal shape can be modeled using Fourier modes to capture features such as the B/P-shaped isophote, enabling more realistic galaxy models. The weight map used in the fitting was automatically generated by GALFIT.

We adopted the exponential profile to model the disk component. Sérsic profiles were used to model the bulge, bar, and barlens. Extensive experimentation reveals that using a Sérsic profile instead of a modified Ferrer function to model the bar does not significantly influence B/T, Sérsic index, and effective radius (r_e) of the bulge (Gao & Ho 2017). For simplicity, we used the Sérsic profile in this work.

Following Gao & Ho (2017), we used two broken exponential profiles to model the ring or ring-like structures, with one component truncated in the inner region and the other in the outer region. The overlap between these two components creates a bump in the surface brightness profile, corresponding to the ring. For both the mock images and CGS observations, the truncation parameters for the disk component were also set following Gao & Ho (2017). We also tested modeling the ring using an exponential disk combined with a Gaussian ring profile (Erwin 2015); the results, presented in Appendix B, are consistent with those obtained using two broken exponential profiles.

The barlens component becomes B/P-shaped at relatively large inclinations (i ≥ 45°) and for a side-on bar (PA_bar ≥ 60°). To model the barlens in this case, we used the Fourier mode in GALFIT to generate a boxy or X-shaped isophote (Peng et al. 2010), adopting the m = 4 mode.

2.3.2. General procedure of the 2D image decomposition

When performing the 2D image decomposition, the initial parameters of each component are important. Following the fitting strategy in Gao & Ho (2017), we started from the simplest fitting configuration with the bar+disk components (bar+D) for Model 1 and the bulge+disk components (B+D) for Model 2. Then we gradually increased the complexity by adding additional components to the image decomposition. The best-fit parameters of the preceding configurations were used as the initial input parameters for the same component of the subsequent configuration. For Model 1, the second fitting configuration, denoted as bar+disk1+disk2 (bar+D1+D2), employed a pair of broken exponential profiles for the disk to account for the ring. Then the barlens component was added in the third fitting configuration denoted as barlens+bar+disk1+disk2 (bl+bar+D1+D2). In addition, we also considered a fourth fitting configuration, with a single exponential disk rather than the two broken exponential disks, namely, barlens+bar+disk (bl+bar+D). For Model 2, the second fitting configuration included a bar component in addition to the B+D, which is the bulge+bar+disk configuration (B+bar+D). In the third fitting configuration, we again used a pair of broken exponential profiles to model the ring structure at the bar end, and denoted it as bulge+bar+disk1+disk2 (B+bar+D1+D2). In the last fitting configuration, we added the barlens component, denoted as bulge+bar+barlens+disk1+disk2 (B+bar+bl+D1+D2).

For both Model 1 and Model 2, the parameters of the components in different fitting configurations were all set free, except for the barlens, which requires special handling. Following Athanassoula et al. (2015), we visually identified the outer boundary of the barlens from the galaxy morphology that appears as a rounder structure extending to about half of the elongated bar. Moreover, to further confirm the position of the outer edge, we also inspected the 1D surface brightness profile along the minor axis of the bar, where a noticeable change in the light profile slope typically marks the transition between the barlens and the surrounding disk. We then fit an ellipse to the visually identified edge of the barlens, and denote its semi-major axis as a_bl. According to Laurikainen et al. (2018), the effective radius of the barlens is then fixed to approximately 0.4a_bl. In Model 2, the bulge overlaps with both the barlens and the bar, leading to significant parameter degeneracy. To address this, additional constraints were applied to the barlens component. Initially, all parameters were set free in the B+bar+bl+D1+D2 configuration for Model 2 to obtain the initial guess values. The effective radius and axial ratio of the barlens were then constrained based on the visually identified outer edge of the barlens ellipse. Moreover, the Sérsic index of the barlens was fixed to the value obtained from Model 1 (n = 1.0), with alternative values (n = 0.5 and 1.5) tested to assess their impact on the bulge parameters. The final Sérsic index for Model 2 was the one that can result in the smallest residual map. The decomposition process proceeded iteratively. First, the parameters of the bulge, bar, and disk were fixed in turn to break the degeneracy and find reasonable values for each component. Once this is done, all parameters of the disk, bar, and bulge were freed in the B+bar+bl+D1+D2 configuration and the fitting was performed again to obtain the final best-fit parameters for each component³.

In contrast to Athanassoula et al. (2015), in our fitting procedure, the size, axial ratio, and Sérsic index of the barlens were always kept fixed. If the size of the barlens is not constrained, the best-fit value may become unreasonable; for example, the barlens covers nearly the entire length of the bar. Similarly, if the Sérsic index of the barlens is not fixed, it would likely become excessively high due to strong degeneracy with both the bulge and the bar.

3. Results

3.1. Mock image decomposition of Model 1

Model 1 is a pure barred disk galaxy without a preexisting bulge component. Therefore, during the modeling process, we only considered the bar, barlens, and disk components. Figure 3 shows the results of the 2D image decomposition for the face-on mock image of Model 1. The best-fit parameters for the barlens and bar components across different fitting configurations are listed in Table 1. Following Gao & Ho (2017), the errors for the free parameters of the components were calculated using the Monte Carlo method. This involved perturbing the sky background 100 times, with the perturbation drawn from a normal distribution with a standard deviation of 25.13 mag arcsec⁻² (twice the standard deviation of the sky background). The error for each parameter was then determined as the standard deviation of the best-fit parameters from these perturbed images, excluding any failed or unreasonable fitting results. For the fixed parameters (e.g., the barlens axial ratio and effective radius r_e), we measured the geometric parameters of the shape of the barlens 50 times and used the standard deviation of the parameters as the corresponding uncertainties.

Fig. 3. Two-dimensional image decomposition results of the Model 1 mock image in the face-on view. In each column, the top row shows the best-fit 2D model image (left) and the residual image (right). The green and blue ellipses indicate the outlines of the bar and barlens, respectively, with semi-major axis of 2re. The other rows show results of the isophotal analysis, including the radial profiles of ellipticity (second row), the position angle (third row), the surface brightness (forth row), and the residual between data and the best-fit model (bottom row). From left to right, the fitting configuration are bar+D, bar+D1+D2, bl+bar+D1+D2 (free parameters fitting), and bl+bar+D1+D2 (fixing the size and axial ratio of the barlens), respectively (see details in the text). The fifth column (bl+bar+D) shows the result of the fitting with only one exponential disk. The RFF of each configuration is labeled in the model image.

To quantify the improvement in the residuals within the barlens region, we use the residual flux fraction (RFF) parameter (Hoyos et al. 2011) defined as

$$\begin{array}{c}\hfill \mathrm{RFF}={\displaystyle \frac{{\sum }_{i,j\in \mathrm{A}}|I_{i,j}-I_{i,j}^{\mathrm{model}}|-0.8\times {\displaystyle {\sum }_{i,j\in \mathrm{A}}}{\sigma }_{i,j}}{{\sum }_{i,j\in \mathrm{A}}{I}_{i,j}}},\end{array}$$(1)

where the A parameter represents the area of the image. Since we are interested in the residual improvement in the central barlens region, we chose the visually identified area of the barlens for the A parameter (for Model 1, see top left panel of Fig. 1). I_i, j is the pixel value of the mock image and $|I_{i,j}-I_{i,j}^{\mathrm{model}}|$ is the absolute value of the residual at each pixel. σ_i, j is the background standard dispersion. Compared to traditional χ² method, RFF offers a noise-corrected measure of residuals; a perfect model with residuals arising solely from background Gaussian noise should yield RFF ∼ 0.

To fit the mock image, we began with the bar+D configuration, with the results shown in the first column of Fig. 3. The residual image for the bar+D configuration shows relatively large residuals in both the central bar region and the outer disk. This is further highlighted in the bottom panel, which shows the surface brightness profile residuals. In the second column (bar+D1+D2), we use a pair of broken exponential disks to model the ring-like structure outside the bar. This significantly improves the residuals in the disk region (R > 25 arcsec). However, there is little change in either the bar parameters or the residual map in the central bar region (R < 20 arcsec). The RFF parameter decreases slightly from 0.061 to 0.052. The Sérsic index of the bar (n_bar) is approximately one, indicating a near-exponential surface brightness profile in the galaxy’s central region. From the residual map, it is evident that both ends of the bar were not properly modeled, suggesting that the best-fit bar configuration did not adequately represent the overall bar light profile.

After adding the barlens component in the bl+bar+D1+D2 configuration (Fig. 3, third column), there is a significant improvement in both the central region of the residual map and the residual surface brightness profile. The RFF parameter decreases from approximately ∼0.05 to ∼0.03. Additionally, n_bar is reduced from ∼1 to ∼0.5, resulting in a flatter profile. The effective radius of the bar (r_e, bar) increases from around ∼10 arcsec to ∼20 arcsec. The Sérsic index of the barlens (n_bl) is ∼0.8, which is close to an exponential profile and steeper than the bar. The barlens is shorter than the bar, with r_e, bl ∼ 0.4r_e, bar. These results for the Sérsic index and relative length of the barlens are consistent with previous studies (Athanassoula et al. 2015; Laurikainen et al. 2018), confirming the reliability of our decomposition method. This also highlights the necessity of modeling the buckled bar with a short barlens (exponential profile) and a long bar (flattened profile). For the disk, the choice between a single exponential disk and a pair of broken exponential disks has little impact on the barlens parameters. After including a barlens component, the scale length of the inner disk increases, and its central surface brightness becomes fainter, while the parameters of the outer disk remain roughly unchanged. Details of the disk parameters can be found in Appendix C.

In another test, we fixed the size and the axial ratio of the barlens in the bl+bar+D1+D2 configuration (Fig. 3, fourth column), using parameters obtained by fitting an ellipse to the visually identified outer edge of the barlens (blue ellipse in the top left panel of Fig. 1). Specifically, we fixed the effective radius of the barlens (r_e, bl) to 7.25 arcsec (∼0.4a_bl, according to Laurikainen et al. 2018) and set the axial ratio to 0.52. Overall, the best-fit parameters show only small changes, with a slight increase in both n_bl and n_bar. Moreover, RFF remains virtually unaffected by fixing the geometric parameters of the barlens. This test confirms that our fitting strategy, which fixes the barlens parameters to break degeneracies, is reasonable. To further evaluate whether the barlens parameters are influenced by the method used to model the disk, and to ensure that the improvement in the central residual primarily results from adding the barlens component rather than from the more complex disk modeling, we also performed a fit with the bl+bar+D configuration (the fifth column of Fig. 3). The result shows that the choice between an exponential disk and a pair of broken exponential disks has little effect on the barlens parameters. Additionally, RFF remains essentially unchanged compared to the bl+bar+D1+D2 configuration in the fifth column, confirming that the improvement in the central residual is primarily due to the inclusion of the barlens component.

3.2. Mock image decomposition of Model 2

Model 2 is more complicated than Model 1, as it includes a classical bulge formed early in the galaxy’s evolution, followed by bar formation at a later stage. Figure 4 shows the results of 2D image decomposition of the face-on mock image of Model 2. Table 2 lists the best-fit parameters of the bulge, bar, and barlens components for different fitting configurations. We first adopted a simple B+D configuration to fit the mock image. The residual map in the left column of Fig. 4 shows the bar and the ring structure clearly. In the second column, after including a bar, the overall residual decreases significantly. For the bulge, n, r_e and B/T all decrease by ∼50%. Additionally, substituting a single exponential disk with a pair of broken exponential disks in the third column further reduces the residuals in the outer disk, resulting in slightly larger bulge parameters. As expected, the residuals continue to decrease as the complexity of the configuration and the number of components increase. However, within the bar region, there is still room for improvement.

After we added a barlens component (i.e., B+bar+bl+D1+D2 with free parameters in the fourth column of Fig. 4), n, r_e, and B/T of the bulge all decrease, with B/T reducing from 0.2 to 0.05. However, in this case (fit the barlens with free parameters), the best-fit size of the barlens becomes unrealistically small, even smaller than the minor axis of the bar. This is physically unreasonable and consistent with the conclusion of Athanassoula et al. (2015) that free parameter fitting for the barlens can yield unrealistic values. To address this issue, based on the best-fit ellipse of the visually identified barlens shape from the image shown in Fig. 1, we fixed r_e, bl to 6.50 arcsec (∼0.4a_bl) and the barlens axial ratio to 0.79, denoted as B+bar+bl+D1+D2 (geo. fixed). The result in the fifth column of Fig. 4 reveals that fixing the barlens size and axial ratio leads to larger bulge parameters compared to the free-fitting scenario. Specifically, B/T increases to 0.1, with n ∼ 1.7 and r_e ∼ 2.6 arcsec. However, the barlens Sérsic index increases to ∼2, which is larger than the expected values based on Model 1 and other previous works. To further improve the fitting, we performed another test by fixing n_bl to 0.5, 1.0, and 1.5, as suggested by previous studies (Athanassoula et al. 2015; Laurikainen et al. 2018; Erwin et al. 2021). We chose the result for n_bl = 1.0, since it exhibits the lowest RFF (Fig. 4, the sixth column, with a configuration denoted as B+bar+bl+D1+D2 (geo. fixed, n_bl = 1.0). Results for n_bl = 0.5 and 1.5 are presented in Appendix A. Regardless of the value of n_bl (free or fixed to 0.5, 1.0, 1.5), barlens-to-total ratio (bl/T) always remains ∼0.15, approximately 50% of the value in the bulgeless Model 1. Additionally, B/T remains low (∼0.1), which is roughly 50% of the value in the B+bar+D1+D2 configuration. The bulge parameters (n and r_e) are reduced in all cases compared to the B+bar+D1+D2 configuration, with changes of roughly ±0.3 in n and ±0.4 arcsec in r_e for different values of n_bl.

Fig. 4. Two-dimensional image decomposition results of the Model 2 mock image in the face-on view. The red, green, and blue ellipses indicate the outlines of the bulge, bar, and barlens, respectively, with semi-major axis of 2re. Other conventions are similar to Fig. 3. The components used in the fitting configurations are shown on the top.

The residual flux fraction is reduced from ∼0.04 in B+bar+D1+D2 to ∼0.03 in B+bar+bl+D1+D2, confirming that including the barlens component improves the modeling significantly. In principle, the barlens is difficult to model because of the spatial overlap of the barlens with the bar and bulge. Our results confirm that to better model the barlens, it is necessary to fix its geometry parameters and constrain its Sérsic index (also see Athanassoula et al. 2015).

For the parameters of the bar, after using a pair of broken exponential disk to model the ring structure, both r_e, bar and the bar-to-total ratio (bar/T) decrease. The Sérsic index of the bar is ∼0.4, consistent with a relatively flat light profile, contributing to about 0.1 of the total flux. Nonetheless, in the four configurations using a pair of broken exponential disks, the best-fit bar parameters remain consistent regardless of the changes in the configuration complexity or additional constraints imposed on the barlens parameters. Regarding the disk parameters, different from the result in Model 1, when the barlens is included, the parameters of both the inner disk and the outer disk remain roughly unchanged.

3.3. Applications on the CGS images

3.3.1. NGC 1533

As shown in the left panel of Fig. 2, NGC 1533, a disk galaxy that is almost face-on, exhibits a short bar, a clear barlens, and an outer ring structure. Following Gao & Ho (2017), we first fit the data with the B+D, B+bar+D and B+bar+D1+D2 configurations. The best-fit parameters for each component are consistent with those in Gao & Ho (2017). Then we added the barlens component and adopt the B+bar+bl+D1+D2 configuration to further refine the fit. The results are shown in Fig. 5 and listed in Table 3. From the left to right columns of Fig. 5, the residual map improves (with a more smoothed residual pattern) due to that we progressively add more structures in the fitting. In B+bar+D1+D2 (the third column in Fig. 5), a clear ring appears in the central region (∼10 arcsec) of the residual map, suggesting that the central region is not well described by the combination of the bar and the bulge. After adding the barlens component (the fourth column), the ring in the residual map disappears and the RFF decreases from ∼0.02 to ∼0.01.

Fig. 5. Two-dimensional image decomposition results of NGC 1533. Other conventions are similar to Fig. 3. From left to right, configurations are arranged in the following order: B+D, B+bar+D, B+bar+D1+D2 and B+bar+bl+D1+D2 (fixing the Sérsic index, size and axial ratio of the barlens).

The axial ratio and size of the barlens were determined following the method outlined in Section 2. Similar to the decomposition strategy used for Model 2, the Sérsic index of the barlens was fixed to 0.5, 1.0, and 1.5, respectively. For n_bl = 1.5, the corresponding result is shown in the fourth column of Fig. 5 since it has the lowest RFF. In this configuration, the including of the barlens component provides the most significant improvement to the residuals in the central region of the galaxy. bl/T is ∼0.25, slightly larger than that in the mock images of Model 2. After including the barlens, the bulge parameters also change in a manner similar to the simulation results. B/T drops from 0.33 in B+bar+D1+D2 to 0.06 after adding the barlens component in B+bar+bl+D1+D2 (geo. fixed, n_bl = 1.5). Similarly, the effective radius also decrease with the additional barlens introduced. However, the Sérsic index of the bulge increases to ∼2.8, different from the results for Model 2. As for the bar component, it has a flat light profile (n_bar ∼ 0.5) and accounts for ∼0.1 of the total flux in B+bar+D1+D2. After adding the barlens component, the Sérsic index and the effective radius of the bar are slightly increased. For the disk component, with the inclusion of the barlens component, the inner disk scale length continues to increase, while the central surface brightness shows little change. The outer disk is less affected.

Different from the results in Section 3.2, fixing n_bl to different values significantly influence the parameters of each component in NGC 1533. B/T decreases from 0.18 to 0.06 as n_bl increases from 0.5 to 1.5, while the Sérsic index of the bulge increases from ∼2.0 to ∼2.8 accordingly. bl/T averages around 0.17, with the lowest value of ∼0.08 for n_bl = 0.5 and highest value of ∼0.25 for n_bl = 1.5. Indeed, accurately modeling the barlens is challenging. In most cases, the bulge parameters would decrease, especially the size and the light fraction B/T. Results for n_bl = 0.5 and 1.0 can be found in Appendix A.

3.3.2. NGC 7329

As shown in the right panel of Fig. 2, the barred galaxy NGC 7329 is observed at a moderate inclination. We followed the fitting strategy in Gao & Ho (2017) for NGC 7329 with the B+D, B+bar+D and B+bar+D1+D2 configurations. The fitting results are consistent with those in Gao & Ho (2017). The best-fit model image, residual map, and isophotal ellipse profiles are shown in Fig. 6 and listed in Table 4. From the left to right columns, as more components are included (from B+D to B+bar+D1+D2), the residual map is improved significantly with less prominent residual patterns. However, similar to NGC 1533, the B+bar+D1+D2 configuration (Fig. 6, the third column) still shows a prominent positive ring feature in the central region (∼5 arcsec) of the residual map, surrounded by a negative ring. This is also evident in the residual profiles, suggesting that an additional component is needed to better model the central region. After adding the barlens component in the B+bar+bl+D1+D2 configuration, the structures in the central region of the residual map, such as the ring in the third column, become less prominent (Fig. 6, the sixth column).

Fig. 6. Two-dimensional image decomposition results of NGC 7329. Other conventions are similar to Fig. 3. From left to right, configurations are arranged in the following order: B+D, B+bar+D, B+bar+D1+D2 and B+bar+bl+D1+D2 (fixing the Sérsic index, size and axial ratio of the barlens).

Again, the barlens parameters were determined following the method in Section 2 and were fixed in the fitting process. Similarly, n_bl was fixed to 0.5, 1.0 and 1.5. We find that using n_bl = 0.5 produces the flattest residual profile and the smallest RFF among the three cases and is therefore shown in Fig. 6 and Table 4. The bulge parameters remain unaffected by the choice of the barlens Sérsic index. The Sérsic index of the bulge ranges between 1.0 to 1.5 for different n_bl and B/T remains ∼0.11. Additionally, changing n_bl does not significantly affect bl/T, which remains ∼0.07. Results for n_bl = 1.0 and 1.5 are presented in Appendix A.

bl/T equals to 0.07 when n_bl is fixed to 0.5, relatively smaller than that for NGC 1533. By comparing the bulge parameters between the B+bar+D1+D2 and the B+bar+bl+D1+D2, we find that including the barlens will slightly decrease the bulge size (r_e of bulge varies from 3.79 to 2.31), flatten the bulge light profile (n of bulge varies from 1.71 to 1.04), and decrease the bulge light fraction (B/T from 0.19 to 0.11). These variations are consistent with the results observed in the simulated mock images and in NGC 1533.

4. Discussion

4.1. Comparison with previous works

Previous studies have explored methods to model the barlens in 2D image decomposition. Athanassoula et al. (2015) used GALFIT to analyze mock images from N-body simulations of barred galaxies and compared them with observations from the Near-IR S0 Galaxy Survey (NIRS0S, Laurikainen et al. 2011) and Spitzer Survey of Stellar Structure in Galaxies (S⁴G, Sheth et al. 2010) samples. They employed generalized elliptical isophotes to fit all components (Peng et al. 2010), with the Sérsic function applied to model the light profile of the bulge and barlens, and an exponential profile for the disk. The bar was modeled using a modified Ferrer function. Although an iterative process was adopted to refine the fit, their modeling strategy only considers the B+bar+bl+D configuration, rather than incrementally increasing the model complexity. As a result, it was not possible to assess the specific impact of the additional barlens component on the bulge parameters. The barlens size and axial ratio in their work were also fixed based on the visually identified barlens shape. However, in their method, the Sérsic index was left free, different from our method. In addition, they also fix the bar length from visual estimation. The ring structure in the disk was also neglected in their fitting.

Later, Laurikainen et al. (2018) performed 2D image decomposition for 46 barlens galaxies from the CALIFA survey using GALFIT. Similar to Athanassoula et al. (2015), they also employed generalized elliptical isophotes for all components and used the same fitting functions, except for the disk, which was modeled with a Sérsic function. Similar to our fitting procedure, they began with a single component fitting and progressively increased the model complexity, but the ring structure was ignored. They initially visually estimated the size and axial ratio of the barlens, which were then fixed in subsequent fitting process. The size of the barlens was only fixed when necessary, while the axial ratio was kept fixed throughout. The galaxy center and disk orientation were also fixed during their fitting process. However, different from our method, their Sérsic index of the barlens was left free, similar to Athanassoula et al. (2015). This may partly explain the presence of some outliers (with large Sérsic index) in their results.

Recently, Erwin et al. (2021) performs 2D image decomposition on two strongly barred galaxies, NGC 4608 and NGC 4643, using IMFIT (Erwin 2015). In their work, the barlens, referred to as the B/P bulge, was modeled using generalized elliptical isophotes. The ring structure was modeled with an additional Gaussian component, while the surface brightness distribution of the outer flat bar was modeled with the default broken-exponential profile. They started with the B+D configuration and progressively increased model complexity, and all parameters in the fitting process were set free. For the two galaxies in Erwin et al. (2021), and the tests mentioned in their appendix, the barlens Sérsic index ranges from 0.5 to 1.1, similar to the values obtained here.

Although the fitting strategies in these previous works differ from ours, their results are generally consistent with our findings. Athanassoula et al. (2015) demonstrated that the barlens can produce a central peak in the galaxy light profile with a Sérsic index less than two, which aligns with our results for the mock images. The measured bl/T ratios in their models range from 0.13 to 0.43, with an average value of ∼0.2, similar to our results. In Laurikainen et al. (2018), the Sérsic index of the barlens also agrees with our findings. On average, the barlens in their sample has n_bl ∼ 1 and bl/T is ∼10 − 30%. However, there are four galaxies with 2.0 ≤ n_bl ≤ 2.5. These relatively large Sérsic indexes of the barlens highlight the necessity of reasonable constraints on n_bl in the image decomposition

For the bulge component, comparing B/T between different fitting configurations, Laurikainen et al. (2018) finds that B/T decreases from ∼0.3 (B+D), to ∼0.15 (B+bar+D), and to ∼0.06 (B+bar+bl+D), which is consistent with our results. For Model 2, we observe a similar trend, with B/T decreasing from ∼0.45 (B+D), to ∼0.25 (B+bar+D1+D2), and to ∼0.10 (B+bar+bl+D1+D2). Similarly, both the bulge Sérsic index and the effective radius decrease after adding the barlens component. In Erwin et al. (2021), compared to B+D, considering a barlens component leads to a decrease in B/T (from ∼0.5 − 0.7 to ∼0.1) and the Sérsic index (from ∼3 to ∼2). Moreover, Erwin et al. (2021) utilized the IFU data to derive the stellar kinematics of the galaxies and assess the reliability of the bulge type classification based on the Sérsic index obtained from the 2D image decomposition with a barlens component. The IFU kinematics confirms the presence of a classical bulge (n ∼ 2.2) in NGC 4608, a nuclear disk (with a broken exponential distribution) with a compact classical bulge (n ∼ 1.6) in NGC 4643, which is identified in the 2D image decomposition.

4.2. Effect of disk inclination

We now investigate the influence of disk inclination (i) on the modeling of the barlens and the bulge. For Model 1 and Model 2, mock images were generated with i = 30°, 45° and 60°, and PA_bar = 0°, 30°, 60° and 90°. The image decomposition strategy described in Section 2 was applied to these mock images. First, we examined Model 1 at different i and PA_bar to study how the barlens properties change with viewing angle. We then analyzed Model 2 to understand how the bulge parameters vary with different inclination angles.

Figure 7 shows the variation of the Sérsic index and the relative size of the barlens (with respect to the bar effective radius r_e, bl/r_e, bar) for Model 1. The top panel suggests that the best-fit surface brightness profile of the barlens is only slightly affected by the disk inclination, with n_bl decreasing from ∼1 to ∼0.9. In the bottom panel, the relative size of the barlens (r_e, bl/r_e, bar) also shows minimal variation with inclination, decreasing slightly from ∼0.41 to ∼0.35 as the inclination angle increases. In addition, for a given inclination angle, the Sérsic index and the relative barlens size across different PA_bar values are roughly similar, suggesting that the bar position angle is less important here. However, for the model with PA_bar = 0° at i = 60°, the best-fit barlens parameters deviate from those of other PA_bar values. In this case, the major axis of the bar is aligned with the minor axis of the projected disk, producing a less elongated bar structure. This alignment results in a serious degeneracy between different components in the central region, leading to a higher Sérsic index and a larger relative size of the barlens than that for other PA_bar.

Fig. 7. Sérsic index (top panel) and relative size (bottom panel) of the barlens for Model 1 at different disk inclinations and bar position angles. The horizontal dashed line in each panel represents the best-fit nbl or re, bl/re.bar from the face-on mock image of Model 1. As the disk inclination angle increases, both the Sérsic index and the relative barlens size decrease. However, for the model with PAbar = 0°, where the major axis of the bar aligns with the minor axis of the projected disk, nbl and re, bl/re.bar increase at i = 60°.

Next, we analyzed the mock images of Model 2 at different inclination angles with the same fitting strategy as Section 2.3.2. We fixed the Sérsic index of the barlens to 1.0 based on the previous results. We also tested with n_bl fixed to 0.5 and 1.5, and found consistent results with those for n_bl = 1.0. The size and axial ratio of the barlens were determined and fixed using the same method as in Section 3.2 for the face-on view. The best-fit configurations of B+bar+bl+D1+D2 (geo. fixed, n_bl = 1.0) and the corresponding residual maps are shown in Fig. 8, with corresponding isophote profiles (ellipticity, PA, and surface brightness) in Fig. 9. These configurations provide a good representation of both the images and the surface brightness profiles at different inclination angles.

Fig. 8. Best-fit B+bar+bl+D1+D2(fixed) model images and residual maps of the 2D image decomposition for Model 2 at different inclinations and PAbar. From top to bottom rows, the inclination angle increases from 30° to 60°. From left to right columns, PA increases from 0° to 90°. The boxy isothote is used for the barlens at a large inclination angle (i = 45° and 60°, and PAbar = 60° and 90°).

Fig. 9. Isophotal analysis showing the ellipticity, PA, B4 (negative for boxyness, positive for diskyness), surface brightness, and residual profile of the data images and the corresponding best-fit model images in Fig. 8. The top, middle, and bottom rows show the results at i = 30°, 45°, and 60°, respectively. From left to right, PAbar increases from 0° to 90°.

Figure 10 presents a comparison of the bulge Sérsic index and B/T for different fitting configurations for Model 2 at various inclinations. For the bulge component, both B/T and n generally decrease with increasing the complexity of the fitting configurations for different inclinations, similar to the results in Section 3.2 for the face-on mock images of Model 2. Our results show that B/T and n of the bulge component are roughly consistent across different inclination angles and PA_bar values. At all inclination angles, B/T is ∼0.5 with the B+D configuration. Once the bar is introduced in the B+bar+D configuration, B/T drops significantly to ∼0.2. Substituting a simple exponential disk with a pair of broken exponential disks leads to a slight increase in B/T. The inclusion of the barlens further decreases B/T to ∼0.1. The Sérsic index n of the bulge also drops from ∼4 in the B+D configuration to ∼1.4 when both the bar and barlens are included. In conclusion, the bulge parameters are not strongly affected by the inclination, while the inclusion of the bar and barlens components significantly reduces both B/T and n values for the bulge.

Fig. 10. Sérsic index (left column) and B/T (right column) of the bulge of Model 2 at different inclinations and PAbar values for different fitting configurations. From top to bottom, the complexity of the fitting configuration increases.

The comparison of bl/T between different inclinations and PA_bar values are shown in Fig. 11. On average, bl/T is ∼0.11. As the inclination angle increases, the difference of bl/T between different PA_bar values increases. At a higher inclination angle (i = 60°) with PA_bar = 0° and 30°, bl/T decreases. This is likely due to that distinguishing between different substructures becomes more challenging when the bar aligns with the minor axis of the projected disk at higher inclinations. Athanassoula et al. (2015) also examined the inclination effect on bl/T and found a larger variation between different PA_bar values at i = 60° (ranging between ∼0.05 and ∼0.12) than i = 30° (ranging between ∼0.11 and ∼0.13), which is consistent with our results.

Fig. 11. Light fraction of the barlens (bl/T) of Model 2 in the B+bar+bl+D1+D2 (geo. fixed, nbl = 1.0) configuration at different inclinations and PAbar values.

At even larger inclinations (i > 60°), determining the barlens parameters becomes increasingly challenging, primarily due to the overlap between the bulge, bar, barlens, and disk. This overlap leads to a more severe parameter degeneracy, resulting in larger uncertainties in the structural parameters. Additionally, the projection effect can change the intrinsic light profiles and isophotal shapes of the different components compared to the face-on case.

Based on the measurements conducted, we can infer a possible correction for the bulge properties of a moderately inclined galaxy (i ≤ 60°) to recover their intrinsic values in the face-on case. For simple fitting configurations (B+D and B+bar+D), i decides the amount of deviation and PA_bar influences the direction of the deviation (whether the values are over or underestimated). The Sérsic index and B/T deviates by 10–30% from the face-on values, with the deviation increasing with i. It is higher than the face-on value at low PA_bar and lower at high PA_bar. For the B+bar+D1+D2 configuration, inclination still introduces systematic deviations, but PA_bar has a comparatively weaker effect on these parameters. Once again, the Sérsic index and B/T exhibits a 10–30% deviation from the face-on value, with the deviation increasing with i. After including the barlens component in B+bar+bl+D1+D2 (geo. fixed, n_bl = 1.0), at low inclinations (i ≤ 30°), the Sérsic index and B/T is systemically higher for ∼15% across all PA_bar than the face-on value. At higher inclinations, the bulge properties remain ∼15% higher than the face-on value at lower PA_bar; for higher PA_bar, the Sérsic index is 10–20% lower than the face-on value, with the deviation increasing with i, while B/T remains 10% higher than theface-on value.

We also tested image decomposition of the edge-on mock images (i = 90°) for both Model 1 and Model 2. Due to the severe degeneracy of these components and the complex light profiles and irregular isophotal shapes of the disk, bar, and barlens in the edge-on view, the overall residuals are too large for the fittings to be considered acceptable. We decided not to include these results in our analysis. The 2D image decomposition of the edge-on galaxies is difficult and challenging, which should be investigated in more detail in the futurework.

4.3. Boxy/Peanut bulge modeling

As mentioned in Section 2, we used the Fourier mode to model the barlens that becomes the B/P bulge for PA_bar = 60° and 90° at high inclinations (i = 45° and 60°) for both Model 1 and Model 2 (as shown in Fig. 8). At these inclinations, visually determining the axial ratio of the B/P bulge is challenging. After performing the initial free parameter fitting for the highly inclined galaxies, the barlens axial ratio was manually adjusted to better match the B/P bulge shape. The residual maps for the configurations with boxy isophotes of the barlens are compared to those with elliptical isophotes in Fig. 12. As expected, using an elliptical isophote results in a distinct X-shaped feature in the residuals, which is properly modeled by the boxy isophote. For both Model 1 and Model 2, compared to previous results for the barlens, using boxy isophotes reduce RFF in the identified barlens region by ∼30% on average (from ∼0.060 to ∼0.037 for Model 1 and from ∼0.025 to ∼0.018 for Model 2), confirming that our method of modeling the barlens with boxy isophotes at large inclinations is reasonable. For both Model 1 and Model 2, the parameters of the bulge, bar, and barlens remain virtually unchanged before and after adopting the boxy isophote for the barlens.

Fig. 12. Comparison of the residual maps of the configurations with the ellipse isophote (left) and the boxy isophote for the barlens component (right). The top and bottom rows show the results for Model 1 and Model 2, respectively.

4.4. BIC test

The existence of the barlens component has been validated in both simulations and observations (Laurikainen et al. 2014; Athanassoula et al. 2015; Laurikainen & Salo 2017; Li et al. 2017; Erwin & Debattista 2017). In the 2D image decomposition, adding extra components can always improve the fitting. However, it also increases the model complexity to result in potential overfitting issue. It is important to assess whether an additional component is truly necessary.

In the literature, several methods have been proposed to quantify the necessity of adding an additional component. For example, Erwin et al. (2021) used the Akaike information criterion (AIC; Akaike 1974) to quantify the improvement in the results after modeling the barlens in the 2D image decomposition. They observed a dramatic decrease in AIC (larger than 1000) after modeling the barlens, highlighting the reliability of their fitting.

The Bayesian information criterion (BIC) is useful to evaluate whether a new component should be added in the fitting (Schwarz 1978). This criterion evaluates the likelihood of a model while applying a penalty for the number of free parameters in the fit. Compared to AIC, BIC imposes a more stringent penalty for an increase in the number of free parameters. According to Head et al. (2014), Méndez-Abreu et al. (2018), BIC statistic can be written as

$$\begin{array}{c}\hfill \mathrm{BIC}=\frac{{\chi }^2}{A_{\mathrm{psf}}}+kln\left(\frac{N}{A_{\mathrm{psf}}}\right),\end{array}$$(2)

where

$$\begin{array}{c}\hfill {\chi }^2={\displaystyle \sum _{i=1}^N}\frac{{[d_i-f_i]}^2}{{\sigma }_i^2}·\end{array}$$(3)

k is the number of free parameters and N is the pixel number of the image. A_psf is the size area of the PSF at the full width at half maximum (FWHM). d_i is the flux of the i-th pixel in the data image, and f_i is the flux of the i-th pixel in the generated model. σ_i² = d_i + S_i, where S_i is the noise. We note that the BIC statistic cannot be the sole criterion for determining the number of components used in the fitting, as its effectiveness is limited by the complexity of real galaxy structures (Argyle et al. 2018) and may lead to inconsistent results compared to visual inspection. The BIC test should be complemented by other strategies when determining the necessity of including a barlens in the galaxy decomposition, such as visual inspection or isophotal ellipse fitting. The parameters from 2D image decomposition, such as B/T, should be also taken into account when employing the BIC test (Nedkova et al. 2024; Bellstedt et al. 2024).

Table 5 shows the values of BIC for different fitting configurations in different galaxies. When counting the number of free parameters used in calculating BIC, the fixed parameters of the barlens are also included. Configurations with lower BIC values are considered superior. Additionally, it is important to examine the BIC variation between different fitting configurations when new components are included. Typically, a decrease in BIC of more than ∼10 is considered sufficient to indicate that a more complex configuration is preferred over a simpler one (Argyle et al. 2018). As presented in Table 5, adding a bar component and adopting a pair of broken exponential disks can decrease BIC considerably (except for NGC 7329 with B+bar+D). For Model 1, we can see a clear reduction of the BIC by modeling the barlens when comparing the bar+D1+D2 and the bl+bar+D1+D2 configurations. For the two CGS galaxies, we can also see an appreciable decrease in BIC after adding the barlens component, demonstrating the necessity of modeling the barlens. However, for Model 2, the barlens component hardly affects BIC. This is expected, as in Model 2 the residual pattern in the central region of B+bar+bl+D1+D2 (geo. fixed, n_bl = 1.0) does not change significantly compared to B+bar+D1+D2. This is different from the results in the CGS real galaxies, where the residuals in the central region become much less significant. The difference may indicate that the barlens in Model 2 exhibits stronger degeneracy with neighboring structures. Therefore, including a barlens component does not significantly improve the BIC. Model 2 could have more significant light concentration in the central region, to result in higher weights during 2D image decomposition. In this case, it is more important to combine the isophotal shapes and the residualsurface brightness profiles to justify the need for a barlens component for Model 2 (RFF decreases from 0.035 to 0.03 after the barlens is added).

4.5. Potential impacts and future prospects

We have shown that the inclusion of the barlens component in 2D decomposition could potentially affect the bulge parameters, including the bulge Sérsic index, effective radius, and B/T. These changes may have non-negligible implications on the well-established scaling relations of the bulges in the barred galaxies.

For example, the Kormendy relation (Kormendy 1977) links the effective radius and surface brightness of spheroidal systems, with the classical bulges following the relation defined by ellipticals (Kormendy 1985, 1987; Bender et al. 1992), and pseudobulges lying below as low surface brightness outliers (Carollo 1999; Gadotti 2009; Fisher & Drory 2010; Gao et al. 2020). After considering the barlens, the effective radius and surface brightness of the bulge are affected to change its position within the Kormendy relation. The bulge classification (pseudo or classical) could also be affected after including the barlens, since typical classification criteria including the Kormendy relation, Sérsic index, and B/T may be sensitive to the adopted configuration. The impact of varying configuration complexities on bulge classification has not been systematically explored. The modeling of the barlens may also affect the bulge mass–size relation. However, previous studies of composite bulges suggest that the overall impact of different decomposition configurations on this relation is relatively modest (Erwin et al. 2015).

Another possible influence is on the B/T distribution along the Hubble sequence. Since the inclusion of the barlens component can significantly reduce B/T, the B/T difference between the early-type and late-type galaxies may become larger depending on the fraction of the barred galaxies in each morphological types. In addition, B/T in the barred galaxies could be smaller than that in the unbarred galaxies, whereas previous studies have found them to be roughly similar (Gao et al. 2019).

The impact of these above mentioned effects on the scaling relations of the bulges remains unclear and will be examined in our future work with larger observational datasets. Given that our current work is primarily focused on evaluating barlens modeling techniques and establishing the necessity of including the barlens component in the 2D image decomposition, we do not attempt to draw broad conclusions regarding the behavior of these scaling relations.

5. Conclusions

In this paper, we have explored the optimal method for modeling the barlens and the impact of this method on the bulge parameters in galaxy 2D image decomposition, using a pure barred disk N-body simulation (Model 1) and a Milky Way-like disk galaxy from the Auriga simulation (Model 2). We also apply our decomposition technique to two selected galaxies from the CGS and examine the effect of inclination on the modeling of both the bulge and barlens. Moreover, to assess the necessity of including the barlens, we conduct BIC tests on the decomposition results. While our sample is limited, this study serves to demonstrate an effective strategy for modeling the barlens in 2D image decomposition, highlighting the importance of accurately accounting for the barlens component to obtain reliable bulge and disk parameters. A more comprehensive statistical analysis based on a larger sample of CGS galaxies will be presented in future work.

The main results are as follows:

1. The barlens has a near-exponential surface brightness profile with n_bl ∼ 0.9 (based on the results in Model 1). In realistic galaxies with a distinct bulge component (as seen in Model 2 and the two observed galaxies), the barlens accounts for 10–20% of the total galaxy light.

2. In 2D image decomposition, modeling the barlens component helps smooth out the central residual patterns, reduces the RFF, and significantly decrease the BIC. To correctly model the barlens, it is essential to constrain the size, axial ratio, and Sérsic index of the barlens. Without these constraints, the fitting may yield unreasonable parameter. The size and axial ratio can be constrained based on the visually identified outer edge of the barlens in the galaxy image. The Sérsic index of the barlens, according to results from Model 1 and previous studies, should be constrained between 0.5 and 1.5.

3. Ignoring the barlens component in 2D image decomposition can significantly affect the bulge measurements. B/T can be overestimated by approximately 30–60%. When the barlens is correctly modeled, only about 10–20% of the total galaxy light is attributed to the bulge component. Other bulge parameters, such as the Sérsic index and effective radius, are also reduced after modeling the barlens component.

4. Both the bulge and barlens parameters are only weakly affected by the inclination of the galaxy or the position angle of the bar. For the barlens component, the Sérsic index of the barlens, n_bl, shows a slight decrease from ∼1 to ∼0.9, and the relative size of the barlens, r_e, bl/r_e, bar, decreases from ∼0.41 to ∼0.35 as the inclination angle increases. The influence of the PA_bar is even less minimal. For the bulge component, both B/T and the Sérsic index decrease with increasing model complexity for different inclinations. These parameters remain relatively consistent across different inclination angles and PA_bar values. However, as the galaxy becomes more inclined, the parameters deviate from the face-on situation significantly, and the spatial overlap of different structures leads to higher uncertainties in the parameter estimates.

5. When the galaxy has a high inclination angle and a large bar position angle, the barlens should be modeled with boxy isophotes to better capture its shape and improve the central residuals.

Acknowledgments

We thank the referee for the valuable suggestions that helped improve our work. We thank Luis Ho, Martin Bureau, and Chengye Cao for helpful suggestions. This work is supported by the National Natural Science Foundation of China under grant Nos. 12233001, 12533004, by the National Key R & D Program of China under grant No. 2024YFA1611602, by a Shanghai Natural Science Research Grant (24ZR1491200), by the “111” project of the Ministry of Education under grant No. B20019, and by the China Manned Space Project with grant Nos. CMS-CSST-2025-A08, CMS-CSST-2025-A09 and CMS-CSST-2025-A11 and in part by Office of Science and Technology, Shanghai Municipal Government (grant Nos. 24DX1400100, ZJ2023-ZD-001). We thank the sponsorship from Yangyang Development Fund. Yang A. Li acknowledges the China Postdoctoral Science Foundation (No. 2023M742285) and the Postdoctoral Fellowship Program of CPSF (GZC20231611). This work made use of the Gravity Supercomputer at the Department of Astronomy, Shanghai Jiao Tong University.

References

Appendix A: Variation of bulge parameters withdifferent barlens Sérsic index

Fig. A.3. Two-dimensional image decomposition results of fixing different nbl in B+bar+bl+D1+D2 for NGC 7329. Other conventions are similar to Fig. 3. The value of nbl fixed in the fitting is indicated on the model image in each column.

Given that fixing the barlens Sérsic index during the 2D image decomposition is a simplified approach, we examine how fixing the Sérsic index to different values impacts the bulge parameters and explore the potential range of these parameters.

The results of Model 2 are shown in Fig. A.1 and listed in Table A.1. Increasing n_bl from 0.5 to 1.5 leads to an increase in the Sérsic index n_bulge from 1.23 to 1.63. B/T remains largely unchanged (∼0.1), while bl/T stays approximately 0.14 for n_bl = 1.0 and 1.5 but drops to ∼0.09 for n_bl = 0.5. The bar parameters are also not severely affected and remain a flat light profile, with the Sérsic index ∼0.4 and bar/T ∼0.07. Moreover, for B+bar+bl+D1+D2 (geo. fixed, n_bl = 1.0 and 1.5), RFF is ∼0.030, lower than the value in B+bar+D1+D2, indicating the reliability of the fitting results. However, for B+bar+bl+D1+D2 (geo. fixed, n_bl = 0.5), RFF is higher than the that in B+bar+D1+D2, suggesting that this is not a successful fitting configuration. This may indicate that the intrinsic value of n_bl for Model 2 lies between 1.0 and 1.5.

Fig. A.1. Two-dimensional image decomposition results of fixing different nbl in B+bar+bl+D1+D2 for Model 2. Other conventions are similar to Fig. 3. The value of nbl fixed in the fitting is indicated on the model image in each column.

Figure A.2 shows the results of NGC 1533 and the corresponding parameters are listed in Table A.2. Unlike in Model 2, the bulge parameters in NGC 1533 are significantly affected by the variation in n_bl. B/T at n_bl = 0.5 (∼0.18) is three times that at n_bl = 1.5 (∼0.06). Compared to B+bar+D1+D2 (B/T∼0.33), this represents a decrease of 40% to 80%. The bulge Sérsic index ranges from 2.06 to 2.82, decreasing at n_bl = 0.5 and 1.0 and increasing at n_bl = 1.5 compared to B+bar+D1+D2. bl/T shows an opposite trend to B/T, increasing from 0.06 to 0.25 as n_bl increases. The bar component is minimally affected by the different barlens models, in contrast to the bulge and barlens components. Although RFF decreases compared to B+bar+D1+D2 in all three cases, n_bl = 0.5 results in the most unsatisfactory residual patterns in the central region of the residual map, suggesting that the intrinsic n_bl for NGC 1533 is likely between 1.0 and 1.5.

Fig. A.2. Two-dimensional image decomposition results of fixing different nbl in B+bar+bl+D1+D2 for NGC 1533. Other conventions are similar to Fig. 3. The value of nbl fixed in the fitting is indicated on the model image in each column.

Figure A.3 shows the results of NGC 7329 and the corresponding parameters are listed in Table A.3. Similar to Model 2, the bulge Sérsic index increases moderately from 1.04 to 1.52 as n_bl increase from 0.5 to 1.5 and B/T remains essentially unaffected (∼0.11). Moreover, Additionally, variations in n_bl have minimal impact on the barlens parameters, with bl/T remaining around 0.07 and the bar parameters showing negligible change. However, for both n_bl = 1.0 and 1.5, RFF is similar to that in B+bar+D1+D2, suggesting that these fittings do not adequately capture the light distribution in the barlens region. This may indicate that the barlens in NGC 7329 has a relatively flat profile, which can be partly attributed to the inclination effect.

Appendix B: Decomposition with a Gaussian ring

Since modeling the ring structure with two broken exponential disks increases the number of free parameters and may introduce degeneracies between the inner disk and the barlens, we tested an alternative approach using a Gaussian ring combined with an exponential disk to represent the ring structure in Model 2. Specifically, we adopted the configurations B+bar+R+D and B+bar+bl+R+D, where R denotes the Gaussian ring. As GALFIT does not support a Gaussian ring component, we performed the fitting using IMFIT (Erwin 2015), which provides this functionality. To ensure consistency, we employed the same optimization algorithm used in GALFIT, namely the Levenberg–Marquardt algorithm. As in previous fits, we fixed the barlens Sérsic index to 0.5, 1.0, and 1.5, and selected the result with the lowest RFF as the final fit.

The results are shown in Fig. B.1 and listed in Table B.1. Compared to the configurations using broken exponential disks, those employing a Gaussian ring yield similar RFF values and residual maps. Overall, the fitting results are consistent between the two modeling approaches, particularly for the bulge parameters and after the inclusion of the barlens component. When the Gaussian ring is adopted, the bar shows a slightly higher Sérsic index and effective radius. Unlike the case with broken exponential disks, the best-fit Sérsic index for the barlens is 1.5; however, bl/T remains nearly unchanged. Considering that the choice of ring model has no significant impact on the results, we adopt the broken exponential disk in the final analysis to facilitate direct comparison with Gao & Ho (2017).

Fig. B.1. Two-dimensional image decomposition results of using the Gaussian ring profile to model the ring structure in B+bar+R+D and B+bar+bl+R+D for Model 2. Other conventions are similar to Fig. 3.

Appendix C: Best-fit disk parameters from the 2D image decomposition

Table C.1 summarizes the best-fit parameters of the disk component for each model configuration applied to the mock images and the two CGS galaxies analyzed in this study. For configurations involving two broken exponential disks, the parameters of the inner and outer disks are listed separately.

For Model 1, when comparing the bar+D and bl+bar+D configurations, the disk scale length slightly decreases after including the barlens component. However, in the case of bar+D1+D2 versus bl+bar+D1+D2, the scale length of the inner disk increases. Additionally, the inner disk in the broken exponential profile exhibits a flatter slope compared to that of the single exponential disk. The central surface brightness of the inner disk and the parameters of the outer disk remain roughly unchanged. For Model 2, the changes in disk parameters after including the barlens component are similar to those observed in Model 1, except that the scale length of the inner disk is slightly reduced.

The inclusion of the barlens component raises the inner disk scale length more significantly for NGC 1533, whereas all the disk parameters remain roughly unchanged forNGC 7329.

All Tables

All Figures

	Fig. 1. Mock images of Model 1 and Model 2 in the face-on view (top), side-on view (middle), and end-on view (bottom). The barlens and the corresponding B/P bulge can be seen clearly in the face-on view and edge-on views, respectively. The blue ellipses in the top panels show the visually identified outer edge of the barlenses.
In the text

	Fig. 2. R-band star-cleaned images of NGC 1533 (left) and NGC 7329 (right) with the barlens component identified in Li et al. (2017). The blue ellipses show the visually identified outer edge of the barlenses.
In the text

Fig. 3. Two-dimensional image decomposition results of the Model 1 mock image in the face-on view. In each column, the top row shows the best-fit 2D model image (left) and the residual image (right). The green and blue ellipses indicate the outlines of the bar and barlens, respectively, with semi-major axis of 2re. The other rows show results of the isophotal analysis, including the radial profiles of ellipticity (second row), the position angle (third row), the surface brightness (forth row), and the residual between data and the best-fit model (bottom row). From left to right, the fitting configuration are bar+D, bar+D1+D2, bl+bar+D1+D2 (free parameters fitting), and bl+bar+D1+D2 (fixing the size and axial ratio of the barlens), respectively (see details in the text). The fifth column (bl+bar+D) shows the result of the fitting with only one exponential disk. The RFF of each configuration is labeled in the model image.

In the text

	Fig. 4. Two-dimensional image decomposition results of the Model 2 mock image in the face-on view. The red, green, and blue ellipses indicate the outlines of the bulge, bar, and barlens, respectively, with semi-major axis of 2re. Other conventions are similar to Fig. 3. The components used in the fitting configurations are shown on the top.
In the text

	Fig. 5. Two-dimensional image decomposition results of NGC 1533. Other conventions are similar to Fig. 3. From left to right, configurations are arranged in the following order: B+D, B+bar+D, B+bar+D1+D2 and B+bar+bl+D1+D2 (fixing the Sérsic index, size and axial ratio of the barlens).
In the text

	Fig. 6. Two-dimensional image decomposition results of NGC 7329. Other conventions are similar to Fig. 3. From left to right, configurations are arranged in the following order: B+D, B+bar+D, B+bar+D1+D2 and B+bar+bl+D1+D2 (fixing the Sérsic index, size and axial ratio of the barlens).
In the text

Fig. 7. Sérsic index (top panel) and relative size (bottom panel) of the barlens for Model 1 at different disk inclinations and bar position angles. The horizontal dashed line in each panel represents the best-fit nbl or re, bl/re.bar from the face-on mock image of Model 1. As the disk inclination angle increases, both the Sérsic index and the relative barlens size decrease. However, for the model with PAbar = 0°, where the major axis of the bar aligns with the minor axis of the projected disk, nbl and re, bl/re.bar increase at i = 60°.

In the text

	Fig. 8. Best-fit B+bar+bl+D1+D2(fixed) model images and residual maps of the 2D image decomposition for Model 2 at different inclinations and PAbar. From top to bottom rows, the inclination angle increases from 30° to 60°. From left to right columns, PA increases from 0° to 90°. The boxy isothote is used for the barlens at a large inclination angle (i = 45° and 60°, and PAbar = 60° and 90°).
In the text

	Fig. 9. Isophotal analysis showing the ellipticity, PA, B4 (negative for boxyness, positive for diskyness), surface brightness, and residual profile of the data images and the corresponding best-fit model images in Fig. 8. The top, middle, and bottom rows show the results at i = 30°, 45°, and 60°, respectively. From left to right, PAbar increases from 0° to 90°.
In the text

	Fig. 10. Sérsic index (left column) and B/T (right column) of the bulge of Model 2 at different inclinations and PAbar values for different fitting configurations. From top to bottom, the complexity of the fitting configuration increases.
In the text

	Fig. 11. Light fraction of the barlens (bl/T) of Model 2 in the B+bar+bl+D1+D2 (geo. fixed, nbl = 1.0) configuration at different inclinations and PAbar values.
In the text

	Fig. 12. Comparison of the residual maps of the configurations with the ellipse isophote (left) and the boxy isophote for the barlens component (right). The top and bottom rows show the results for Model 1 and Model 2, respectively.
In the text

	Fig. A.3. Two-dimensional image decomposition results of fixing different nbl in B+bar+bl+D1+D2 for NGC 7329. Other conventions are similar to Fig. 3. The value of nbl fixed in the fitting is indicated on the model image in each column.
In the text

	Fig. A.1. Two-dimensional image decomposition results of fixing different nbl in B+bar+bl+D1+D2 for Model 2. Other conventions are similar to Fig. 3. The value of nbl fixed in the fitting is indicated on the model image in each column.
In the text

	Fig. A.2. Two-dimensional image decomposition results of fixing different nbl in B+bar+bl+D1+D2 for NGC 1533. Other conventions are similar to Fig. 3. The value of nbl fixed in the fitting is indicated on the model image in each column.
In the text

	Fig. B.1. Two-dimensional image decomposition results of using the Gaussian ring profile to model the ring structure in B+bar+R+D and B+bar+bl+R+D for Model 2. Other conventions are similar to Fig. 3.
In the text