Galaxy–galaxy strong lensing simulation with GPU acceleration across surveys and multi-bands

1 Introduction

In the Stage-IV sky image surveys, approximately O ( 1 0 5 ) galaxy–galaxy strong lenses (GGLs) are expected to be discovered (Collett, 2015; Acevedo Barroso et al., 2025) by ground- or space-based telescopes, such as the Legacy Survey of Space and Time (LSST) (Ivezić et al., 2019) and Euclid (Euclid Collaboration: Mellier et al., 2025). These expected samples will enable investigations into galaxy structures (Koopmans et al., 2006), the mass function and projected cumulative mass fraction of substructure (Vegetti and Koopmans, 2009), time-delay measurement for the Hubble constant (e.g., Suyu et al., 2010), the details of the magnified high-z lensed galaxy (e.g., Brammer et al., 2012), and cosmological parameters with high statistical significance (Li et al., 2024). In particular, GGLs serve as unique probes for constraining dark matter properties (e.g., Tortora et al., 2010; Sonnenfeld et al., 2015).

However, the large number of candidates hidden within the vast volume of survey data poses new challenges for lens finding and modeling. High-efficiency lens-finding and modeling tools are required to manage the influx of images from current and future surveys. In recent years, various machine learning (ML)-based tools have been developed to tackle these challenges, pioneered by early works using convolutional neural networks (CNNs; Hezaveh et al., 2017; Perreault Levasseur et al., 2017).

For lens finding, current advanced tools include classification CNNs (e.g., Petrillo et al., 2017; Lanusse et al., 2018; Li et al., 2021; Nagam et al., 2025; Storfer et al., 2026), self-attention-based networks (e.g., Thuruthipilly et al., 2022; Thuruthipilly et al., 2024; Parlange et al., 2026), and U-Net-based or ensemble methods (e.g., Nagam et al., 2024; Inchausti et al., 2026; Zhong et al., 2025). These tools achieve superior accuracy and efficiency compared to traditional methods applied in past surveys.

For lens modeling, approaches have evolved from direct parameter estimation via regression CNN (e.g., Schuldt et al., 2023) to more complex architectures, such as Bayesian neural networks (BNNs; e.g., Gentile et al., 2023; Euclid Collaboration: Busillo et al., 2026), differentiable probabilistic models (e.g., Gu et al., 2022; Cao et al., 2025), and simulation-based inference with neural density estimation (e.g., Legin et al., 2023) or normalizing flows (e.g., Mishra-Sharma and Yang, 2022).

Despite these advancements, ML-based methods require a vast amount of training data tailored to specific survey characteristics. For the joint analysis of multi-band data across different surveys (e.g., combining LSST and Euclid), it is essential to simulate identical lens systems that are physically consistent across different instrument response functions. Furthermore, reliable lens mass estimation imposes strict requirements on photometric redshift accuracy, which requires a simulation capable of faithfully reproducing the multi-band color relations of galaxies.

This study is a direct extension and further development of our previous work presented in Zhong et al. (2025). We developed the GGSL-UNet architecture and demonstrated its capability to extract two-dimensional structural information from multi-band images (Zhong et al., 2025). Building on the deblending framework of GGSL-UNet, the current work introduces a fully GPU-accelerated, tensor-based lens simulation pipeline. While Zhong et al. (2025) focused on the neural network’s architecture and downstream deblending performance, this manuscript advances the methodology by achieving full physical consistency through synthetic spectral energy distributions (SEDs) and expanding the application to joint, cross-survey modeling across Kilo-Degree Survey (KiDS), LSST, and Euclid specifications.

To address these needs, we develop a highly efficient, tensor-based pipeline that leverages GPU acceleration. In this framework, the simulated lens encodes color information from synthesized galaxy spectral energy distributions (SEDs) at the given redshift, enabling the generation of high-fidelity, multi-band training sets for cross-survey analysis.

The article is organized as follows: Section 2 briefly introduces the simulation model implemented using a GPU framework; Section 3 introduces the different survey configurations adopted in the simulation; Section 4 compares with other simulations and the application of the deblending technique; and Section 5 presents our conclusions.

2 Lens image simulations

The simulation of a GGLs image is, in principle, performed by providing an expected source flux and a deflector mass model, and then applying Poisson electron statistics on the charge-coupled device (CCD) camera. Meanwhile, additional noise, such as readout noise, is superposed as Gaussian fluctuations. In this work, we implement a fully differentiable lens simulation pipeline using PyTorch tensors (Paszke et al., 2019). This allows deep learning algorithms and GPU acceleration to be directly integrated into the modeling procedure, constrains lens parameters using gradient-based modeling techniques, or directly searches for the minimal-loss optimal lens modeling parameters via error backpropagation.

2.1 Lens model

The simulated lens system consists of four components: the foreground deflector light, the mass model, the background source light model, and the shear. This simulation is numerical and based on a light ray originating from a 2D tensor array of fixed-coordinate sub-grids in the lens plane. This tensor array will be directly or indirectly substituted into the above-mentioned part. The whole simulation is then tensorialized, meaning the numerical units and operators are implemented using PyTorch tensors and operators at the GPU level, enabling gradients to be propagated to the physical parameter level. The source light is modeled by the Sérsic profile (Sérsic, 1963; Ciotti and Bertin, 1999; MacArthur et al., 2003), while the projected 2-dimensional mass distribution is modeled by the cored isothermal sphere (CIS) profile (Narayan and Bartelmann, 1996) as follows: where I 0 is the central surface brightness amplitude, θ ⃗ represents the radius vector of the 2D angular position in the image plane, θ h is the effective (half-light) radius, n is the Sérsic index, b n ≈ 2 n − 1 / 3, θ E is the Einstein radius, and s is the core radius. For CIS, | θ ⃗ | represents the elliptical coordinates | θ ⃗ | 2 = q x 2 + y 2 / q, where q is the axis ratio. The light model of the lens, as well as the sky flux in units of electron (or analog-digital unit (ADU)) flux, reflects the intrinsic properties of the electron counts in the CCD. The foreground and background light models are independent, with different Sérsic parameters, denoted as F f and F b, respectively. Those parameters are sampled from the distributions in Table 1, where we pick up the typical elliptical galaxy with the redshift range of 0.0–0.6 as the lens, and the star-forming galaxy with the redshift range of 0.8–2.4 as the source, both of which mostly have a Sérsic index between 1.0 (1.5) and 6.0, and the half-light radius range around 0 . 5 ″ to 3 . 0 ″ ( 2 . 0 ″ ) at redshift equal to 0.5. Table 1 provides further details. The lens/deflector center ( x 0 , y 0 ) is randomly selected within −0.5 to 0.5 arcsec of the image center. After the lens center is determined, the source center ( x 1 , y 1 ) is randomly selected within −0.4 to 0.4 arcsec relative to the lens center. The deflection angle to the background source light is α, calculated by the lensing potential Ψ and the shear potential Γ (Narayan and Bartelmann, 1996) via: where the κ ̃ is the Fourier transformation of κ, F − 1 is the inverted Fourier transformation, ( x , y ) represents the coordinates, and ( γ 1 , γ 2 ) represents two components of the external shear. The lensed background light is obtained by remapping the source using the transformed coordinate grid, F b ( θ ⃗ − α ⃗ ). The foreground and background systems will be rotated and shifted independently to introduce the additional location parameters listed in Table 1. Finally, the total expected flux (in units of electrons) from the lens on the detector is the sum of the lens light and the lensed source light as follows:where the c f / b is the scaling factor that controls the strength of the signal of the modeled lens or source, which can be obtained by setting an expected magnitude.

F θ ⃗ = I 0 exp − b n | θ ⃗ | / θ h 1 / n , (1)

κ θ ⃗ = θ E 2 | θ ⃗ | 2 − s 2 , (2)

Ψ = F − 1 − 2 κ ̃ / k 2 , (3)

Γ = 1 2 x 2 − y 2 γ 1 + 2 x y γ 2 , (4)

α ⃗ = ∇ ⃗ Ψ + Γ , (5)

F l θ ⃗ = c f F f θ ⃗ + c b F b θ ⃗ − α ⃗ , (6)

TABLE 1. The lens parameters in the r band. A similar table can be found in Li et al. (2021).

The ranges represent sampling a value from flat uniform distributions.

The half-light radius is set to be measured at the reference redshift z 0 = 0.5, then linearly scaled by foreground redshift, z f, or background source redshift, z b.

For some parameters (e.g., ( x 2 , y 2 ) and ( x 0 , y 0 )), we will take the same values.

In summary, the entire lens modeling pipeline on the GPU begins with an initial 2D tensor grid of coordinates, with channels corresponding to different photometric bands (e.g., for KiDS, four channels corresponding to the u , g , r , i bands), which is then substituted into Equations 1–6 above using custom-defined tensor operations. Compared with lenstronomy running on a CPU, the GPU-based simulation is more than three orders of magnitude faster (∼ 30 s / i m a g e / c o r e vs. 0.02 s / i m a g e / G P U)⁴, for images with the same resolution of 64 × 64 pixels and a pixel scale of 0 . 2 ″.

2.2 Redshift and multi-band model

To simulate realistic multi-band observations, we adopt the r band as the reference for the structural parameters. For other bands, the light models of the foreground and background sources are allowed to vary within perturbations less than 20% in 3 σ in the half-light radius θ h and Sérsic index n to mimic wavelength-dependent morphological variations (Vulcani et al., 2014). We adopt the Gaussian distribution to cover all situations and avoid imposing a prior constraint on the trained models. These perturbations Δ are sampled as vectors from a normal distribution N described in Equations 7–9: where the vectors Δ , Θ h , n represent the parameters across different bands (e.g., 4-dim vector for KiDS simulation), while θ h r and n r are the values in the r band. I is an identity matrix. While the light profiles vary across bands, the parameters of center location, primary angle, and axis ratio constrained from the r-band mass model remain fixed for a given system. Given the redshifts of the foreground and background sources, we will synthesize their SEDs, shift them to the given redshifts, and sample the amplitudes in the different wavelength bands using the filters according to Equations 10–12. The magnitude in other specific bands is derived by multiplying the redshifted SED by the corresponding filter transmission curve. where the T ( ω ) is the multi-band filters, and A r is the amplitude of SED in the r band. We assign redshifts to the foreground and background galaxies and synthesize their SEDs using the flexible stellar population synthesis (FSPS; Conroy et al., 2009; Conroy and Gunn, 2010), assuming a Chabrier 2003 initial mass function (IMF). The star-formation history (SFH) is modeled with a τ-model, in which the star-formation rate (SFR) follows an exponential decline, SFR ∝ e − t / τ for simplicity, characterized by the timescale τ.

Δ ∼ N 0 , 0.2 I , (7)

Θ h = θ h r 1 + Δ , (8)

n = n r 1 + Δ , (9)

A = ∫ S E D ω T ω d ω , (10)

Δ m = − 2.5 ⁡ log 10 A / A r , (11)

m = m r + Δ m , (12)

In this article, we adopt distinct stellar population models for the foreground and background SED to reflect the typical GGL configuration: a “red” passive foreground galaxy (e.g., an early-type galaxy lens in Bolton et al., 2006) and a “blue” star-forming background source (e.g., Li et al., 2021; Euclid Collaboration: Walmsley et al., 2026). The foreground SED assumes a quiescent galaxy, while the background SED assumes an active star-forming galaxy.

For the quiescent galaxy, the stellar age is set to the age of the Universe at the given redshift, assuming a cosmology from Planck Collaboration et al. (2020). For simplicity, we assume the low-redshift quiescent galaxy with an age of approximately 8–10 Gyr (Gallazzi et al., 2005), such that its post-quenching SED is effectively similar to those of stellar populations as old as the age of the Universe. The parameter log 10 τ is sampled from [−1.0, 0.0] (rapidly quenched within 1 Gyr; Carnall et al., 2018). The logarithmic metallicity relative to the solar value, log 10 ( Z / Z ⊙ ), is sampled from the fixed grid [−0.25, 0.00, 0.25] (metal-rich; see Peng et al., 2015). Star formation is required to cease earlier than 5 Gyr, and the dust attenuation is kept low, with E B − V drawn from [0.0, 0.2] (Falco et al., 1999), which is assumed to be dust-poor.

For the star-forming galaxy, the stellar age is sampled from 1 Gyr up to the age of the Universe. The parameter log 10 τ is sampled within [0.0, 1.0] (adopted from empirical e-folding time; e.g., Asari et al., 2007; Ciesla et al., 2017); the metallicity, log 10 ( Z / Z ⊙ ), is sampled from the fixed grid [−1.00, −0.75, −0.50, −0.25, 0.00, 0.25] (which includes metal-poor regimes; see Peng et al., 2015); the dust reddening E B − V is sampled from [0.0, 0.6] (dust-rich; see Calzetti, 2001; Xiao et al., 2012); and nebular emission lines are included in the SED (primarily H α, H β, O II, and O III; see Byler et al., 2017). We ensure that the background redshift is larger than the foreground redshift. More complex SED models can also be adopted, including those based on directly observed spectra.

2.3 Noise and observation model

The observational simulation incorporates instrumental effects, including point spread function (PSF), sky background, dark current, readout noise, thermal noise, and other sources. The counted electrons from the sources and the sky generally obey Poisson statistics, and the thermal and readout noise usually obey Gaussian statistics. The noise model is separated into two parts: one is the Poissonian statistical electron (CCD electrons associated with photon counts) from the lens signal, and the sky is denoted as N e, and the other is the contributions of read-out noise and thermal noise, which are denoted as N σ.

The PSF is modeled as a Gaussian kernel defined by full width at half maximum (FWHM) using Equations 13 and 14. We note that the PSF varies not only with band but also with atmospheric turbulence. Here, the Gaussian kernels approximate the typical value across different bands, using a fixed FWHM as a baseline. The future FWHM (Full Width at Half Maximum) will be treated as a learnable parameter within the modeling framework, with GPU acceleration, to adapt to input image quality, especially for ground-based observations. The expected number of electrons collected by the CCD pixel at the corresponding position is given by the convolution of the flux with the PSF, plus the sky background contribution, as shown by Equation 15:where T is the exposure time, and the sky flux, F s, is set as a constant dependent on the limiting magnitude of the survey.

P S F ∝ exp − Θ ⃗ 2 2 σ 2 , (13)

σ = F W H M 2 2 ⁡ ln ⁡ 2 . (14)

N ̄ e = F l ⊗ P S F + F s ⋅ T , (15)

The received total electrons are sampled from the Poissonian distribution given an expectation of source and sky electrons plus a Gaussian fluctuation. After the sky subtraction, the final observed sky-subtracted electrons will then be converted to ADU in the CCD. On balance, we model the photon noise using a Poisson distribution and add a Gaussian component to represent readout and thermal noise described using Equations 16–19: The strength of N e is proportional to T, σ 2 is the total noise variance per pixel, I l is the final image in ADU, G is the gain, and we assume the sky background is subtracted in the final product. In this simulation, we ignore cosmic rays, artificial light, and optical “ghosts” (Euclid Collaboration: Paterson et al., 2026) for simplicity.

N e ∼ P o i s s o n N ̄ e , (16)

N σ ∼ G a u s s i a n 0 , N ̄ σ , (17)

σ 2 = N ̄ e + N ̄ σ 2 , (18)

I l = N e + N σ − F s T ⋅ G . (19)

In general, a survey’s exposure time is fixed. To obtain an image with the expected magnitudes of foreground and background, we adjust the flux amplitude using c f / b. The magnitude is calculated by the total number of electrons from the source following Equation 20:where A is the area of the aperture, and Z P is the zero points of the specified survey. The source strength, c f / b, can be calculated if given a specified magnitude. The summation is over all the pixels in the cut-out area. The signal-to-noise (SNR) of the foreground or background is defined as Equation 21:The sky flux is calculated by the limiting magnitude, which is the electrons integrated in the 5 σ level in a given aperture (2 arcsec in KiDS, Hildebrandt et al. (2017)), following Equation 22:where the m l i m is the limiting magnitude, A is the area of the aperture, and we ignore the read-out noise.

m f / b = − 2.5 ⁡ log 10 ∑ A c f / b ⋅ F f / b ⋅ G + Z P , (20)

S N R f / b = ∑ A F f / b ⋅ T ⋅ G ∑ A σ . (21)

F s = 1 0 − 0.8 ( m l i m − Z P ) 25 A , (22)

In a real survey, the co-added image is a stack of several to tens of individual exposure images of the same target to enhance its SNR. In our simulation, we accordingly increase the effective exposure time T to n T to mimic the stacked-enhancement SNR situation for simplicity, where n denotes the number of exposure images. As a consequence, the stacked process will also average the readout and thermal noise, reducing them to 1 / n under ideal conditions.

2.4 Simulation pipeline and numerical accuracy

Algorithm 1 Forward Simulation Pipeline: Initialized with a 2D coordinate grid tensor, followed by lens modeling, SED rescaling, noise injection, and sky subtraction.

The complete forward simulation pipeline for the lensed image is described in Algorithm 1, where each stage (1–8) is computed sequentially from the results of the preceding stage. Our simulation only considers the science-level sci images, in which the PSF, theoretical electron-counting Poisson noise, readout Gaussian noise, and sky subtraction are accounted for, while the CCD or optical models are beyond the purview of the lens’s science application and are tied to each survey’s pipelines. The final image is thus a deterministic function of the coordinate grid θ ⃗ and the physical parameters listed in Tables 1 and 2. The closed-form expression of the simulation can be written as I ( θ ⃗ ∣ θ param ). When the coordinate grid and parameters are tensorized, the entire closed-form expression becomes vectorized and differentiable. Using a fixed grid, we can apply variational optimization to obtain the best-fit solution, which is efficiently achieved via backpropagation with respect to the parameters θ param.

TABLE 2. Parameters of the noise and observation model.

Because the KiDS flux is normalized, we use a fixed zero point as the demonstration.

For comparison with KiDS, we used the 10-year coadded image, which has the highest SNR.

Those limit magnitudes are for one exposure image, while the 10y LSST depths are deeper.

Zero points of LSST are from https://smtn-002.lsst.io.

The pixel scale of the near-infrared (NIR) band should be 0.3; however, in the GPU tensor, the coordinate grid should be identical for the gradient calculation, so we take a uniform 0.1 for each band.

The optical and NIR bands of Euclid are I E (550–900 nm), Y E (920–1146 nm), J E (1146–1372 nm), and H E (1372–2000 nm).

Thereby, to assess the effect of grid resolution on the resulting image, we generate simulations using oversampling factors of 1 ×, 5 ×, and 10 × and adopt the 20 × oversampled grid as the reference against which the root mean square (RMS) is computed. After downsampling to the KiDS resolution, Figure 1 shows example multi-band images produced at each oversampling factor. Averaged over 100 simulated lenses with identical parameters, the RMS relative to the 20 × reference is 1677.5, 324.9, and 203.4 for the 1 ×, 5 ×, and 10 × grids, respectively, indicating that the simulation converges and becomes numerically stable as the grid resolution increases.

3 Surveys

3.1 KiDS

The KiDS (de Jong et al., 2013, KiDS) is a prominent Stage-III optical imaging survey primarily designed for weak lensing cosmology. This survey covers a footprint of 1,347 d e g 2, overlapping with the VISTA Kilo-degree Infrared Galaxy Survey (VIKING; Edge et al., 2013), providing nine bands of optical to near-infrared (NIR) photometrical (u, g, r, i, Z, Y, J, H, K s) images of the target, with a pixel scale of 0.2 arcsec/pixel.

In this work, we adopt the r band as the reference band, which typically exhibits the best seeing conditions (PSF FWHM= 0.7 0 ″). We simulate the four visible band (u, g, r, i) images with the corresponding pixel scale, PSF, exposure time, depth, and zero point consistent with KiDS Data Release (Hildebrandt et al., 2017; Wright et al., 2019; Wright et al., 2024). The specific simulation parameters are detailed in Table 2, and corresponding simulation examples for different surveys are shown in Figure 2.

3.2 LSST

The LSST (Ivezić et al., 2019), conducted by the Vera C. Rubin Observatory, defines the ground-based Stage-IV era. It is based on the instrumental specifications of the Simonyi Survey Telescope, an 8.4 m class telescope with a three-mirror anastigmat optical design. This configuration is notable for providing sharp images over an exceptionally wide 9.6 d e g 2 field of view, which is a feature that influences the spatial variation of the PSF across the focal plane. It adopts an effective light-collecting area equivalent to that of a 6.42-m single-mirror telescope.

The LSST camera (LSSTCam) is a 3.2-gigapixel mosaic of 189 4K × 4K CCDs with a nominal pixel scale of 0.2 arcsec per pixel. Our simulation generates images for the six broadband LSST filters (u, g, r, i, z, and y), which collectively span the wavelength range from 320 nm to 1060 nm. To generate simulated strong lensing images consistent with LSST, we simulate the observational characteristics of the 10-year coadded data from the main Wide-Fast-Deep (WFD) survey. The noise properties of the simulated images are calibrated to match the final coadded depth of the 10-year WFD survey, which is planned to cover approximately 18,000 deg 2 of the southern sky. The noise level in each band is set to correspond to the official 5σ point-source limiting magnitudes for the 10-year coadded dataset (see Table 2).

3.3 Euclid

The European Space Agency’s Euclid mission (Euclid Collaboration: Mellier et al., 2025) is designed to replicate the capabilities of its two primary instruments: the visible instrument (VIS) and the near-infrared spectrometer and photometer (NISP).

The VIS instrument provides high-resolution imaging through a broad optical filter, the I E band, which covers 550–900 nm. Its focal plane applies a mosaic of 36 CCDs with a native pixel scale of 0.101 arcsec/pixel. Although the VIS PSF is slightly undersampled by the pixel grid, the Euclid observing strategy mitigates this by employing a dither pattern (typically four exposures per sky-pointing), allowing us to recover sub-pixel information during image coaddition. This process yields an effective angular resolution of 0.18 arcsec, with a system PSF FWHM (excluding pixelization effects) of 0.16 arcsec in coadded images.

The NISP instrument complements VIS by providing imaging in three NIR bands: Y E (949.6–1212.3 nm), J E (1167.6–1567.0 nm), and H E (1521.5–2021.4 nm). The NISP focal plane is composed of 4× 4 NIR HgCdTe detectors with a coarser pixel scale of 0.3 arcsec/pixel, which fully samples the PSF. For our simulation, we adopt the in-flight FWHM values measured during the performance verification phase: 0.35″ for the Y E band, 0.34″ for the J E band, and 0.35″ for the H E band.

The simulations are designed to replicate the data products from the Euclid Wide Survey, which is planned to cover more than 14,000 deg 2 of the extragalactic sky. The combination of VIS and NISP data poses a unique data-fusion challenge due to their heterogeneous nature. The VIS channel is expected to provide high-resolution (0.1″ pixels, 0.16″ FWHM) images, while the NISP channels are expected to provide lower-resolution (0.3″ pixels, ∼ 0.35” FWHM) images (see Table 2). However, the final Q1 data of Euclid have a median FWHM of 0.20″, 0.48″, 0.50″, and 0.54″ in I E, Y E, J E, and H E, respectively (Euclid Collaboration: Romelli et al., 2026; Marleau et al., 2026). These degenerations compared with simulations must be taken into account for future applications. The key simulation parameters for both instruments are listed in Table 2 (Euclid Collaboration: Cropper et al., 2025; Euclid Collaboration: Jahnke et al., 2025). To facilitate pixel-level operations in our tensor-based pipeline (e.g., gradient calculations), we simulate both VIS and NISP bands on a unified coordinate grid with a pixel scale of 0.1 arcsec/pixel (see Table 2), noting that NISP data can be downsampled to its native resolution if required.

4 Results

4.1 Comparison with lenstronomy

In this section, we compare our simulated images with those from the other lens simulation tool, lenstronomy (Birrer and Amara, 2018; Birrer et al., 2021). We use the lens parameters from our r-band image as input to lenstronomy to generate the corresponding lens image. The r-band images generated by our tools and lenstronomy are shown in Figure 3. In addition to the noise and amplitude differences (set according to the sky survey and SED), we can visually check that the simulated images are almost identical in geometry.

Compared with lenstronomy, we can simulate lens images more efficiently and incorporate color relations among multi-band data based on the types and redshifts of foreground and background galaxies. Based on the magnitude, zero points, and so on, those tuning parameters will make the lens image simulations and redshift measurements more friendly and directly correspond to surveys. In reality, a more sophisticated color gradient should be considered; thus, n and R h should be correlated to wavelength, type, and metallicity. For the modeling and ML sides, we choose to sample from the Gaussian distribution to cover the lens case as much as possible and avoid implicitly imposing a prior on the trained models. The tensor-based simulation can also let the user utilize the gradient information for lens modeling with GPU acceleration.

4.2 Application: deblending with GGSL-UNET

To demonstrate the utility of our simulation pipeline and assess the quality of the generated data, we employ the GGSL-UNet architecture (Zhong et al., 2025) to perform galaxy deblending and estimate the photometric magnitudes and redshifts of both foreground and background galaxies. The goal is to accurately disentangle foreground lens light from the background source in the noisy, PSF-convolved 2D multi-band observed images, which is a prerequisite for accurate photometry and mass estimation.

We generated a dataset of 40,000 multi-band images for each of the three surveys (KiDS, LSST, and Euclid) using the simulation pipeline described in Section 2. The dataset was split into a training set (90%) and a validation set (10%). As part of our simulations-modeling pipeline, we integrated and trained the GGSL-UNet using the L1 loss (mean absolute error, MAE) as the objective function. The optimization was performed with an initial learning rate of 1 0 − 3 and a decay rate of Γ = 0.5 every 10 epochs. Figure 4 illustrates the training and validation loss curves over 50 epochs. The consistent convergence of both training and validation losses indicates that the model effectively learns the morphological features of the lens systems without overfitting. The model with the minimum validation loss was selected for subsequent performance evaluation.

Given the input images, the GGSL-UNet outputs the foreground and background galaxy images, and contamination sources (not the lens system source) are also removed if present. In this work, we evaluated the trained models on a separate test set of 1,000 samples for each of the three surveys. Figure 5 shows the input noise images, simulated foreground/background components, deblended outputs reconstructed from GGSL-UNet with reduced chi-squared, and the residuals. Visual inspection and reduced chi-squared confirm that the network successfully separates the overlapping galaxy components and suppresses background noise and artifacts. To quantify the reconstruction quality, we compute the reduced chi-squared χ ν 2 statistics for the total image χ 1 2, the foreground component χ 2 2, and the background component χ 3 2. These are defined as Equations 23–25. where I l is the input mock observed image (simulated with noise), while I f and I b are predicted deblended foreground and background images. The term d o f represents the number of degrees of freedom, which is approximated by the number of pixels. The term σ denotes the noise map (standard deviation) of input lens images, and M f / b represents the noise-free mock models for foreground or background, which is expected for the total number of electrons (expected value) received on the CCD. It is worth noting that while the reduced χ 2 statistics are employed to evaluate the statistical pixel-level residuals, our GGSL-UNet architecture was natively trained under the L 1 loss (MAE) to preserve the global flux and sharp morphological boundaries. The combination of the L 1 optimization target and the post-training χ 2 test ensures both physical flux conservation and proper noise-level statistical consistency.

χ 1 2 = ∑ b a n d , p i x e l I l − I f − I b 2 σ 2 / d o f , (23)

χ 2 2 = ∑ b a n d , p i x e l M f − I f 2 σ 2 / d o f , (24)

χ 3 2 = ∑ b a n d , p i x e l M b − I b 2 σ 2 / d o f , (25)

Figure 6 shows the distributions of the chi-squared values for the test sets, with mean values reported in the top-right corner of panels. Overall, as shown in the figure, most chi-squared values for the deblended images are less than 1.0, indicating that the residuals are consistent with the noise level and that the deblending is statistically robust. For the impact of spatial resolution, the space-based telescope Euclid achieves the lowest average χ 2 and the smallest spread. This is attributed to its superior spatial resolution (e.g., VIS channel, 0.1″ pixels with PSF FWHM = 0.16″), which minimizes morphological degeneracy between the lens and source and allows for more precise deblending than ground-based surveys with relatively worse PSFs. One notable aspect is that while LSST images have a significantly higher SNR than KiDS due to greater depth, their average χ 2 values are comparable or slightly higher. This is an expected result, and there is a reason for it: in high-SNR regimes like LSST with small σ, the metric becomes highly sensitive to even minute systematic residuals (e.g., slight structural mismatches), whereas in lower-SNR regimes like KiDS, larger σ can mask these residuals, resulting in analytically lower χ 2 values. Nevertheless, visually and physically, the LSST reconstruction benefits from the higher depth (see Figure 5). For this reason, the Euclid Deep Survey, which observes a few fields with longer exposures, is expected to yield greatly enhanced-quality deblended images.

5 Conclusion

In this work, we have presented a GPU-accelerated, tensor-based framework for generating and deblending realistic strong-lensing images. This framework is designed to generate and deblend realistic, multi-band galaxy–galaxy strong lensing (GGSL) images across different survey configurations and to provide gradient information at the parameter level, enabling the user to constrain model parameters. The iterative backpropagation constraint or the Monte Carlo process can be accelerated on the GPU using this tensor-based framework. The primary motivation is to enable robust joint analysis of multi-survey data and to facilitate mass estimation in the absence of spectroscopic follow-up by implementing a self-consistent SED and redshift for the foreground and background galaxies, as well as the structure parameters in simulations. Another core innovation of our pipeline lies in its physical consistency: foreground and background sources are modeled using synthetic SEDs that incorporate star-forming and quiescent templates, ensuring that the simulated magnitudes and colors are self-consistent with the redshift evolution of the galaxies. The framework is highly flexible, allowing users to configure more user-friendly observational parameters for photometric redshift measurements corresponding to surveys, such as pixel scale, PSF, limiting magnitude, and zero points, to mimic arbitrary survey conditions.

Furthermore, we note a few limitations of the current framework that impact the generation of even more realistic lensing systems. First, the light profiles of both lens and source galaxies are modeled using parametric Sérsic profiles. While computationally efficient, this approach lacks the complex stellar substructures (e.g., spiral arms, tidal debris, and clumpy star-forming regions) present in real imaging data. Second, the mass distribution is simplified as a CIS profile. In reality, complex environmental effects (such as line-of-sight structures and external shear) can introduce subtle geometric distortions. Future developments will incorporate nonparametric pixelated source models or hydro-simulated galaxy templates to mitigate these limitations and further enhance the physical realism for next-generation automated lens modeling. Finally, the tensor coordinates are fixed and uniformly distributed over the field of view; denser grids should be assigned to higher-density regions, and adaptive grid accommodation at the lens center should be considered in future developments.

As a proof of concept, we applied this pipeline to simulate identical lens systems under the conditions of KiDS, LSST, and Euclid surveys, and then utilized a U-Net-based architecture (GGSL-UNet; Zhong et al., 2025) to perform image deblending. The outputs of the GGSL-UNet process exhibit a consistent loss across surveys and converge to stable values as the number of training iterations increases. The deblending chi-square calculation for KiDS and LSST yields an expected chi-square of approximately 1.0, and Euclid achieves the lowest average chi-square due to its better PSF, demonstrating that the different observational configuration simulations across surveys exhibit a consistent pattern. This shows the superior benefit of the follow-up foreground and background photometry redshift measurements by combining those deblended images. Although the images from the real survey are noisier, as contamination or artificial sources may be included in the lens system cut-out, the demonstration in GGSL-UNet shows that the neural network can successfully manage and remove those contaminants. Finally, the fully differentiable nature of our tensor-based simulation with GPU acceleration is a critical feature. It ensures that the pipeline is not only computationally efficient on GPUs but also ready for subsequent ML algorithms, such as integration into end-to-end, gradient-based inference frameworks (e.g., differentiable forward modeling) for future strong lensing science.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. Generative AI, Gemini, was used during the preparation of this manuscript for the purpose of language editing, proofreading, and refining sentence structure to enhance readability. No scientific content, data, or simulations were generated by AI. The authors reviewed, verified, and revised all AI-generated content and take full responsibility for the accuracy and integrity of the manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Jirong Mao, Yunnan Observatories, National Astronomical Observatories, Chinese Academy of Sciences (CAS), China

Reviewed by: Mario Gai, Osservatorio Astrofisico di Torino (INAF), Italy

Ruslan Gabbasov, Autonomous Metropolitan University, Mexico

Footnotes

Lenstronomy was run on a laptop AMD Ryzen 9 7940 CPU, while the GPU simulation was run on a laptop NVIDIA RTX 4070.