Flow chart of the fully differentiable method used in this work and implemented in HERCULENS, which merges analytical profiles with pixelated components. We indicate the typical number of parameters N for the analytical potential η_ψ and the pixelated potential ψ_pix. Partial derivatives are computed and propagated for all parameters. These derivatives are then used to update parameter values in the direction indicated by the gradient of the loss function L. We note that the blurring operator B (see Eq. (4)) is also used to generate the model m but is not shown in the diagram to avoid clutter. This work focuses on modeling the lensing potential ψ, and we refer to Joseph et al. (2019) and Galan et al. (2021) for the modeling of surface brightness distributions.