Unscrambling the lensed galaxies in JWST images behind SMACS0723

The first deep field images from the James Webb Space Telescope (JWST) of the galaxy cluster SMACS~J0723.3-7327 reveal a wealth of new lensed images at uncharted infrared wavelengths, with unprecedented depth and resolution. Here we securely identify 14 new sets of multiply imaged galaxies totalling 42 images, adding to the five sets of bright and multiply-imaged galaxies already known from Hubble Space Telescope data. We find examples of arcs crossing critical curves, allowing detailed community follow-up, such as JWST spectroscopy for precise redshift determinations, and measurements of the chemical abundances and of the detailed internal gas dynamics of very distant, young galaxies. One such arc contains a pair of compact knots that are magnified by a factor of hundreds, and features a microlensed transient. We also detect an Einstein cross candidate only visible thanks to JWST's superb resolution. Our parametric lens model is available through the following link: https://www.dropbox.com/sh/gwup2lvks0jsqe5/AAC2RRSKce0aX-lIFCc9vhBXa?dl=0, and will be regularly updated using additional spectroscopic redshifts. The model is constrained by 16 of these sets of multiply imaged galaxies, three of which have spectroscopic redshifts, and reproduces the multiple images to better than an rms of $0.5^{\prime \prime}$, allowing for accurate magnification estimates of high-redshift galaxies. The intracluster light extends beyond the cluster members, exhibiting large-scale features that suggest a significant past dynamical disturbance. This work represents a first taste of the enhanced power JWST will have for lensing-related science.


INTRODUCTION
At long last, the first deep images from the James Webb Space Telescope (JWST) were delivered on July 11, 2022. The target is the central region of SMACS J0723.3-732 (SMACS0723 hereafter; z d = 0.39), which is part of the southern extension of the MACS sample (Ebeling et al. 2010;Repp & Ebeling 2018). SMACS0723 is a massive strong-lensing (SL) galaxy cluster which is rich in galaxy images distorted by the gravitational lensing effect, as was seen in recent Hubble Space Telescope (HST) imaging taken as part of the Reionization Lensing Cluster Survey (RELICS; Coe et al. 2019). The P lanck-derived mass is high 8.39 × 10 14 M , and there were nine galaxies at z > 5.5 and two candidates at z = 7 based on their photometric redshift estimation (Salmon et al. 2020;Strait et al. 2021). Three far-infrared sources have also been detected in the field using the Herschel Space Observatory .
A lens model for this cluster based on these previous HST images was published recently (Golubchik et al. 2022). This model uses the Light-Traces-Mass (LTM) method, which assumes that the galaxy light traces the underlying stellar and dark matter (DM) but does not assume a parametrized mass function for the DM (Broadhurst et al. 2005;Zitrin et al. 2009Zitrin et al. , 2015. In Golubchik et al. (2022), the authors identify five strongly lensed galaxies and derive the spectroscopic redshifts for three of those systems using publicly-available ESO Very Large Telescope (VLT) Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) data. Two additional models are also available, one based on Lenstool (Jullo et al. 2007;Sharon et al. 2022) and the other on Glafic (Oguri et al. 2012). These models, however, do not use the spectroscopic information from MUSE. This first set of images of SMACS0723 from JWST represents a milestone for not only astronomy, but science in general. As with its predecessor, HST, this first image from JWST reveals the distant universe with incredible detail. Some of the most prominent arcs detected with the HST RELICS program now show multiple unresolved substructures in JWST previously unseen (see for instance system 4 in Figure 1). These substructures facilitate the identification of systems of strongly lensed images, where an image system consists of all the gravitationally-lensed counter-images from the same background galaxy. The greater sensitivity of JWST to small flux fluctuations and improved spatial resolution enable searches for the fainter caustic or micro-caustic crossing events of stars at z < 2 such as Icarus (Kelly et al. 2018), or Warhol (Chen et al. 2019;Kaurov et al. 2019), stars between 2 < z < 5 such as Godzilla , stars at z ≈ 6 such as Earendel (Welch et al. 2022) or even further up to the first stars as proposed in Windhorst et al. (2018).
JWST is already providing breakthrough results in the study of the most distant galaxies. Pointing JWST towards gravitational lenses results in a more powerful combined telescope with an effective diameter a factor |µ| times larger, where |µ| is the absolute magnification of the lensed background object. For galaxies, typical magnifications can reach factors of a few tens, making the combination JWST+SL cluster similar to a space telescope analogue to JWST but with with 20-30 m diameter. In the context of high-z object detection, the lensed source counts are initially lower than a blank field due to the reduction in search volume from magnification. At the same time, the blank field luminosity function flattens with higher redshift, allowing the number of lensed galaxies to overtake the number of blank field galaxies (e. g., Mahler et al. 2019;Salmon et al. 2020;Pascale et al. 2022).
For much smaller background objects like magnified stars near caustics or in general for sub-pc structures near caustics, the underlying magnification factor can be of order ∼ 1000. This translates into effective apertures for JWST of ≈ 200 meters! Even without the ultra-high magnification boosts near the critical lines, SL has allowed us to probe several magnitudes deeper than blank fields, i. e., down to rest-frame UV luminosities M UV −13 magnitudes (Bouwens et al. 2017(Bouwens et al. , 2022aLivermore et al. 2017;Atek et al. 2018;Ishigaki et al. 2018;Vanzella et al. 2021) and stellar masses M 10 6 M (Bhatawdekar et al. 2019;Kikuchihara et al. 2020;Furtak et al. 2021;Strait et al. 2021). With the redder wavelength range, greater sensitivity, and resolving power of the JWST, we can expect many more galaxies at even higher redshifts to be detected in the coming months. In order to fully characterize them and study their physics, it is beneficial to have different models which are constructed in independent works. After this paper was submitted to the archive, we have become aware of two new lens models based on Lenstool which include both MUSE spectroscopic redshift constraints Figure 1. Image stamps depicting the 19 image multiplicities in the SMACS0723 field, as labeled. The publicly-available color image of SMACS0723 is depicted, which is valuable for the identification of image family members by their similar colors and morphological components. We also require that the model be able to reproduce the positions of the members of each image family. A 5 bar is shown for reference. and JWST imaging (Mahler et al. 2022;Caminha et al. 2022).
In this paper we identify 14 new image families using the JWST data and present a new parametric SL model. The paper is organized as follows: In §2 we detail the NIRCam imaging and photometry needed to provide constraints for the lens model, which in turn is discussed in §3. In §4 we present and discuss the results. The conclusion appears in §5. Throughout this work we use a ΛCDM cosmology with Ω M = 0.3, Ω Λ = 0.7, and H 0 = 70 km s −1 Mpc −1 . Unless otherwise stated, we use AB magnitudes (Oke & Gunn 1983), and errors correspond to 1σ.

OBSERVATIONS AND DATA
Observations with the Near Infrared Camera (NIR-Cam; e. g., Rieke et al. 2005) aboard the JWST were executed as part of the Director's Discretionary Time on 2022 June 06 (PI: Pontoppidan; Program ID 2736). Exposures were taken in F090W, F150W, and F200W in the short-wavelength (SW), and F356W, F277W, and F444W in the long-wavelength (LW) channels, totalling 12.5 h of integration time. Near Infrared Spectrograph (NIRSpec; Jakobsen et al. 2022) observations were made in the same field as a part of the same program. Two exposures were taken in each of the F170LP and F290LP filters, totalling 11787 sec of integration time in each filter. SMACS0723 was also observed using the Mid Infrared Instrument (MIRI; Bouchet et al. 2015;Rieke et al. 2015) and the Near Infrared Imager and Slitless Spectrograph (NIRISS; Doyon et al. 2012). The JWST data analyzed in this work can be found on MAST: 10.17909/rwdx-k029.
The raw data from all four instruments were processed using the JWST science calibration pipeline which performed the background subtraction, flat-fielding, correc-tion for cosmic ray hits, correction for image distortions, re-pixelization by the drizzle approach onto a common astrometric reference frame, co-addition to make the mosaic, and combination into a color image. For the NIRCam data, we reduce the raw data products independently from uncal files to i2d mosaics using a custom pipeline based on JWST pipeline version 1.6.2 alongside with the Calibration Reference Data System (CRDS) version 0942. The pipeline includes updates based on on-flight calibrations that account for zero point offsets between the filters and the NIRCam modules (Rigby et al. 2022), critical for accurate photometry and derived properties. Furthermore, this pipeline includes additional steps over the default JWST pipeline to improve on its astrometry and background subtraction. For a detailed description on the additional steps see Adams et al. (2022).
Using the NIRCam images based on our independent reduction with improved calibration, we extract photometry for our targets of interest using forced photometry of small, 0.32 arcsecond diameter circular apertures centered on the peak flux of the target in each band. The photometry is aperture corrected for a point spread function using the simulated point spread functions (PSFs) from WebbPSF (Perrin et al. 2014). This information is then fed into the photometric redshift code LePhare (Arnouts et al. 1999;Ilbert et al. 2006). LePhare operates by fitting a large grid of galaxy spectral energy distributions to the photometry provided. We run LePhare with a suite of galaxy templates from (Bruzual & Charlot 2003), allowing for E(B-V) values between 0 and 1.5 (Calzetti et al. 2000), redshifts between 0 and 12 and applying the (Madau 1995) treatment for absorption from the IGM. The redshifts provided in Table 1 are the redshift of peak probability and errors generated from a χ 2 grid produced by LePhare based on the templates, redshifts and dust attenuation combinations available. Several entries are missing values as a result of various sources of contamination or by a lack of detection. We refer to the Appendix for more details. We note that these NIRCam data can be compared with bluer data using HST taken as a part of the RELICS program (PI: D. Coe), and the Spitzer Space Telescope (PI: M. Bradac).
We also made use of the available NIRSpec data set. In particular, we analyzed those spectra which were acquired at the positions of our image systems to confirm our image system designations. We refer to §3.1 for more details. We also consulted the MIRI images to test for achromaticity within each of our image systems, for those systems that were detected. Close-up view of system 5 with the z = 1.43 critical curve marked in the dashed line on this F090W+F150W+F200W composite color image. Two compact sources in the inset, knots A and B, have a separation of ∼ 0.16 , and have sizes that are consistent with the width of the PSF. The knots bracket the critical curve that must pass in between them at an angular separation of ∼ 0.08 . We estimate high magnification factors as a result of this close proximity to the critical curve of |µ| ∼750 for each source. A red knot is also detected at a position slightly offset from knot A, Tr, which may be a microlensed transient. We refer to §4 for more information regarding this potential microlensed transient.
Overall, the arcs in this field make up a rich tapestry of distorted images typical of lensed sources and at a level of detail that has never been seen before. The increased spatial resolution reveals lensing constraints even down to the small substructures within the arcs, since many of these small scale sources are equally multiply-imaged, such as the double arc in system 5 (Figs. 1 and 2). Meanwhile on the scale of the cluster, thanks to the remarkably dark infrared sky background of JWST images, we can get a high signal-to-noise view of the intra-cluster light (ICL). This cluster has an ICL that is elongated and nonuniform. We refer to §4 for a discussion of how the ICL's relatively unique features suggest that this cluster is not dynamically-relaxed.

SL MODELING OF SMACS0723
We generate a new SL model exploiting the JWST data in which we identify the new sets of multiple images, as is described in §3.1. Cluster member galaxies used in the modeling are described briefly in §3.2. The modeling method is described in §3.3.

New multiple images in SMACS0723
We identify new lensing constraints taking advantage of the superior quality of the JWST images. We start by adopting the first five multiple image systems published in Golubchik et al. (2022), following their numbering scheme. These systems serve also as a model- Figure 3. Critical curve at z = 1.45 from the SMACS 0723 lens model is overlaid as a color composite image generated from the six existing JWST NIRCam filters depicting the central 1.5 by 1 arcminutes of the cluster. The model incorporates strong-lensing constraints from the 19 image families identified in total, 14 of which are discovered by this work using the JWST data). The two "+" symbols mark the positions of the two dark matter halos used in this model, and the 4 magenta rings mark the galaxies whose weight is left free in the model. Image 2.2 is rare for containing ≥ 10 bright knots located at large projected galactocentric radii (inset). independent guide to search for new systems. For example, galaxies at close projected separations and at similar redshifts typically maintain similar spatial distributions when distorted by a gravitational lens. Moreover, amongst the galaxy images, counter images with negative and positive parities are in general symmetric to one another with respect to the critical curve. After visual inspection of the new JWST/NIRCam image, and corroboration by available photometric and/or spectroscopy redshift constraints and model-predicted locations, we identify 14 new sets of multiply lensed galaxies which are reported in Table 1.
Fortunately, there are spectroscopic redshifts available for some of the image systems. The redshifts of image systems 1, 2, and 5 are z = 1.450 ± 0.001, z = 1.378 ± 0.001 and z = 1.425 ± 0.001, respectively and were measured with MUSE (Golubchik et al. 2022). Note that while JWST/NIRSpec spectroscopic redshifts of 35 galaxies have been reported in this field with redshifts as high as 8.3 as reported by the JWST press release, an investigation of the NIRSpec spectroscopy un-covered spectra for only one of our 42 images from our 19 image families: image 4.1. However, an inspection of the spectrum in each of the spectral band-passes did not produce any significant features that were present in multiple visits and therefore did not yield a redshift. While valuable, further spectroscopic redshifts are however not crucially necessary for securing our image system identifications which relies on a large set of image multiplicities vetted in several ways. We note that NIRISS spectroscopy of SMACS0723 is also available but not used in this study. Its redshift constraints will inform future lens models.

Cluster member galaxies
For our model we use the same the cluster member galaxy selection as in Golubchik et al. (2022), which was based on the red sequence shown therein (their Fig. 4). In addition, we add here 3-4 galaxies not used in the initial model by Golubchik et al. (2022) which lie near the bright star 20 north of the BCG (which hindered somewhat their previous detection), and are included as they may affect the reproduction of nearby counter images. A revision of this selection using spectroscopic data (such as e.g., the NIRISS data), is deferred to future work.

SL modeling method
The LTM model presented in Golubchik et al. (2022) required a very strong external shear, indicating that the cluster's central mass distribution is highly elongated. Since the LTM approach is limited in the intrinsic matter ellipticity it can accommodate, and for comparison, we make use here of a new, fast parametric method that we recently constructed (A. Zitrin; in preparation; the method can also be referred to as analytic, i. e., it is not limited to a grid's resolution). The method is similar in nature to other parametric lens modeling techniques such as Lenstool (Jullo et al. 2007), Glafic (Oguri et al. 2010) or GLEE (Halkola et al. 2006;Grillo et al. 2015). In such parametric methods, dark matter halos can be elliptical which more easily accounts for the required elongation (which is also evident from the distribution of arcs).
The models in our new parametric method have two main components. First, the cluster members galaxies are modeled with double pseudo elliptical mass distributions (dPIE; Elíasdóttir et al. 2007) which are assumed spherically symmetric with the exception of the BCG, and are defined following the prescriptions of Jullo et al. (2007) and Zitrin et al. (2013): Where the σ is the velocity dispersion, r cut is the cutoff radius, r core is the core radius, and L * is the typical luminosity of a galaxy at the cluster redshift. We fix r * core = 0.2 kpc for all dPIEs, use an L * equivalent to a galaxy of m ref F 814W = 22.13, and assume a constant mass-to-light ratio (α = 0.5), while σ * and r * cut are left free to be optimized by the model.
The second modeling component is the set of cluster DM halos. These can in principle be modeled either as elliptical Navarro-Frenk-White (NFW) profiles (Navarro et al. 1996), or as pseudo elliptical mass distributions (PIEMD; e. g., Monna et al. 2015); or dPIEs. The third and last component which can be added if necessary is a two-component external shear.
The method is similar to our previous parametric implementation outlined in Zitrin et al. (2015) which has been well vetted, but is not limited to the assigned grid resolution. The new improved version used here has been already applied to various clusters and has also been tested on simulated clusters, accommodating both image-and source-plane minimization. The minimization of the model is done via a Monte-Carlo Markov Chain (MCMC) with a Metropolis-Hastings algorithm (e. g., Hastings 1970). We include annealing in this procedure, and the chain typically runs for several dozen thousand steps after the burn-in stage. Errors are calculated from the same MCMC chain.

SMACS0723 SL model
For modeling SMACS0723 by this method, we constrain the model using 16 of the 19 identified image systems and a total of 48 images (2 sets of multiple images are used from system 5; see Table 1), including spectroscopic redshifts for 3 systems and leaving all other system redshifts to be fit by the model. We represent cluster members as dPIE profiles, as explained above, and the cluster DM halos as PIEMDs (i. e., the main free parameters for each halo aside from the ellipticity and its position, are a core radius and a velocity dispersion). We use two DM halos, centered on the two central galaxies respectively (marked by "o" symbols in Fig. 3). We do not incorporate an external shear. Since galaxies can deviate from the assumed scaling relations, we leave the four brightest galaxies to be freely weighted (marked by rings in Fig 3), which means each of their velocity dispersions σ is scaled by a weighting parameter optimized by the model (see Table 2; G1-G4). The BCG's ellipticity and position angle are also left free to be optimized by the model. In total, there are 29 free parameters in the model: the normalization for the galaxy velocity dispersion σ * , the normalization for the galaxy cut-off radius r * cut , 4 free parameters for each of the 2 dark matter halo PIEMDs, 4 galaxy weights, the BCG ellipticity and position angle, and the 13 redshifts for systems without spectroscopic redshifts. For the lensed sources, each multiply-imaged system should converge to a single position in the source plane. To enforce this constraint, we minimize the differences between the image positions in the source plane following the prescription of Keeton (2010), which usually converges to a solution relatively quickly and without loss of quality, as shown by Keeton (2010) in their successful Hubble Frontier Fields (HFF; Lotz et al. 2017) lens models. Given the 48 lensed images, this constraint amounts to 62 SL constraints on the model (two for each lensed image with the exception of one image per multiplicity), which gives a total of 33 degrees of freedom for the model.

RESULTS AND DISCUSSION
The best-fit SL model, which is the one for which χ 2 is minimized, is presented in Fig. 3. Critical lines are overlaid for z s = 1.45 (i. e., the redshift of image system 1). The figure illustrates the elongated configuration of SMACS0723. The model was minimized using a positional uncertainty of σ pos = 0.5 , resulting in a χ 2 of 159.4 (χ 2 ν = 4.8) in the image plane, and an rms of 0.48 in reproducing the positions of the multiple images; the optimized model parameters can be found in Table 2.
The effective Einstein radius at a given redshift is computed as θ E = A/π, with A being the area enclosed within the critical curves. At z = 2, θ E = 18.4 ± 1.8 , and the mass contained inside the critical curve is (5.91 ± 0.83) ×10 13 M , which is overall consistent with the values of θ E = 16.9 ± 2 , and M = (4.15 ± 0.58)× 10 13 M obtained in the HST LTM model (Golubchik et al. 2022). The uncertainties on the Einstein radius and mass are similar for each of the models and are limited by the systematics.
The JWST/NIRCam images uncover fine details of the faint substructures within the arcs. In particular, two compact sources in system 5, knots A and B (5.A and 5.B in Table 1), are spatially-resolved on opposite sides of the critical curve, constraining the position of the critical curve in our lens model (Fig. 2). The image pair has a close angular separation to the critical curve of 0.08 , implying extremely high magnifications for each image and making it an ideal system in the search for caustic crossing events. Near to the critical curve, magnification follows an inverse relationship with distance µ = µ 0 /D crit , where µ 0 is dependent on the slope of the lensing potential and can be measured from the lens model (Welch et al. 2022;Diego 2019). Uncertainty in µ 0 is accounted for by sampling the model MCMC chain, while the uncertainty in D crit is assumed to be the uncertainty in the centroid position of the images. Our parametric model implies a magnification of µ = 737 +1553 −454 for each image, whose positions are reproduced to within 0.02 by the model.
Magnification near the critical curve is highly sensitive to systematics, and may vary significantly between modeling approaches ). Meneghetti et al. (2017) demonstrated that the uncertainty in magnification can reach greater than 30% at µ > 10 across all modeling approaches. Hence statistical errors in magnification from the model are most likely underestimate the true error which is dominated by systematics. Assuming these images are unresolved, this magnification could imply parsec or even subparsec sizes for this object. Given the arc's close proximity to the cluster center, however, microlensing from intracluster stars could work to smooth out the critical curve, placing an upper limit on the persistent magnification possible (Venumadhav et al. 2017;Dai 2021).
Consider again the image pair consisting of knots A and B in Fig. 2. There is a red knot offset from knot A that is detected only in some of the NIRCam bands, called Tr (see appendix section B). We present the argument that Tr is neither knot A nor knot B. Assuming knots A and B are counter images, the pair will always be detected together in any given band. Since the separation of these mirrored knots is about 0.16 , which is larger than the FWHM of any JWST filter in this dataset (F444W, FWHM = 0.145 , is the largest), JWST should be able to spatially resolve this pair. We detect knots A and B in all the SW filters, but detect only a single peak in each of the LW filters near to the position of knot A (see appendix, section B and Fig. 6). Since the image pair is not detected together through the LW filters, this single peak must refer to Tr. Hence, Tr is distinct from knots A and B. The lack of a counterimage, and its placement in close proximity to the critical curve, suggest Tr is a microlensing event. However we cannot rule out that it is a red foreground galaxy. Given the very red colors and small angular sizes in- Figure 5. Projected mass density of SMACS0723. We show κ, the surface mass density distribution in units of the critical density for strong lensing, for the source redshift of system 1, at z = 1.45, plotted on a grid of constant RA and Dec (dashed lines). The resulting mass distribution exhibits a clear ellipticity along the same axis as the ICL as shown in Fig. 4 . volved, confirmation of Tr as a transient would require JWST follow-up observations. We note that there are no other similar dual-color detections in any other objects in system 5 that would suggest an image misalignment.
The giant arcs also present some rather unusual features, such as system 2 with a spectroscopic redshift of z = 1.38 first measured by Golubchik et al. (2022). Image 2.2 has >10 prominent clumps all of a similar red color and extending to high galactocentric radii (see Fig. 3, inset). While significant sub-kpc structure of star formation regions is not unusual for star-forming galaxies at z = 1 − 3 (or even sub-100pc, e. g., Johnson et al. 2017), the very large separations cast some doubt that the extended knots are identical in nature to those in what might be the plane of a galactic disk. These clumps also appear to be different (redder) than those typically detected in clumpy galaxies (e. g., Shibuya et al. 2016). The physical characteristics are more consistent with numerical predictions for globular clusters. For example, Sameie et al. (2022) found in hydrodynamic simulations that the longest lived globular clusters likely formed at high redshifts beyond the half-light radius (up to hundreds of parsecs away from their host galaxy). Pozzetti et al. (2019) go on to predict the color evolution for high redshift globular clusters, but also point out that NIRCam colors are most likely insensitive to the forma-tion redshift. Vanzella et al. (2017) have reported detections of candidate young (<10 Myr) globular clusters which are blue and at high redshifts of z = 3-6. Based on this information, it is tempting to speculate that the clumps seen in system 2 may be redder and older (∼Gyr) counterparts of those young candidate globular clusters detected at the much higher redshifts. Based on SED-fitting from 0.4 − 4.4µm, Mowla et al. (2022) found evidence that the clumps were consistent with evolved globular clusters with ages of 3.9 − 4.1 Gyr, implying formation a mere ∼0.5 Gyr after the Big Bang. We note that another possible explanation, is that the red clumps fitting the description of those seen in system 2 may arise as a result of ram pressure stripping of dwarf galaxies in the dense environments near in projection to its disk (e. g., Mayer et al. 2006;Boselli et al. 2022). Further follow-up, such as detailed source plane reconstruction and photometric SED modeling, is needed to uncover the nature of these sources.
A striking feature of this new JWST/NIRCam image is the morphology of the baryonic component responsible for the intracluster light (ICL), which exhibits a loop-like feature in the northwest component and a large lobe in the southeast component. We perform our own reduction of the raw images to mitigate the variations in sky background across chips as described in §2, although some background noise can also be seen in each of the four corners of the image due to unoptimized sky subtraction in the JWST reduction pipeline (Fig. 4). The ICL is due to stars stripped away from their galaxies but still gravitationally bound to the cluster. Earlier work has suggested that the ICL is a good tracer of the DM distribution (or vice-versa) since both stars in the ICL and DM are expected to behave as collisionless particles and thus respond only to gravity (Montes & Trujillo 2018). Our initial SL model exhibits a similar macro-structure to the ICL as seen in Fig. 4. The ICL also extends beyond the current set of SL constraints, but in combination with weak lensing derived from the JWST images it should be possible to produce a lens model that covers the entire range of the ICL in a way that gives insights as to the intriguing apparent correlation between the ICL and the DM distribution. In particular, the presence of these large-scale features in the ICL may motivate detailed simulations to examine the merging scenarios that constrain the time since the major merger and estimate the relative velocities involved.
We note that after this paper was submitted to the archive, and prior to submission to ApJ, we have become aware of two other papers on SL modeling in this field (Mahler et al. 2022;Caminha et al. 2022). There are some differences in the results presented in these other papers, but we have checked and found that including their results into our lens model does not alter the results of this study.

SUMMARY
We presented a new JWST parametric model for the massive cluster SMACS0723, incorporating new multiple image family constraints identified in the JWST press-release color image and in the NIRCam data set. Our model builds upon the existing lens models (Golubchik et al. 2022) by increasing the SL constraints with 14 newly identified image families, which represents a factor of three improvement over the HST model. This parametric model has a large impact because it provides a starting point for modeling SL clusters in the era of JWST. The mass map reveals a somewhat complex and extended mass profile and, in the detection of the image pair knots A and B opposite the critical curve in system 5, also opens a route to the study of caustic transients in the discovery of an apparent microlensed transient adjacent to knot A. The ICL has a shape that is broadly-similar to that of the galaxies and the model, but with rare larger-scale features that appear to retain the memory of its dynamical history, a question which can be investigated further with numerical simulations. In the future, spectroscopic and photometric redshifts will bolster the reliability of our lens model, allowing also a more precise placement of the critical curve needed to engage in caustic transient studies.

ACKNOWLEDGEMENTS
We wish to thank the anonymous referee, whose suggestions have improved this manuscript. We would also like to thank the RELICS and PEARLS collaborations for data products or discussions that have stimulated this work, and of course the people behind the JWST, its instruments and public release of this incredibly beautiful first data set.
This work is based on observations made with the NASA/ESA Hubble Space Telescope (HST) and NASA/ESA/CSA James Webb Space Telescope (JWST) obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST, and NAS 5-26555 for HST. These observations are associated with program 14096 for HST, and 2736 for JWST. This work is also based on observations made with ESO Telescopes at the La Silla Paranal Observatory obtained from the ESO Science Archive Facility. The authors thank the Centre for Astronomy and Particle Theory (CAPT) of University of Nottingham for providing all computational infrastructure to run the JWST Reduction pipeline, and Philip Parry for technical support. LF acknowledges financial support from CAPES -Brazil. We also thank Anthony Holloway and Sotirios Sanidas at JBCA for critical help with computer infrastructure that made this work possible.
This research made use of Astropy, 1 a communitydeveloped core Python package for Astronomy (Astropy Collaboration et al. 2013;Price-Whelan et al. 2018) as well as the packages NumPy (van der Walt et al. 2011), SciPy (Virtanen et al. 2020, matplotlib (Hunter 2007), and some of the astronomy MATLAB packages (Ofek 2014).

A. ARC SYSTEMS IN SMACS0723
We present the list of all multiply-imaged systems below. The image families are vetted in multiple ways: by their similar morphological components, similar colors, similar redshifts (photometric and/or spectroscopic), and/or consistency with the lens model. The IDs are given in the first column, where the "?" indicates an image member candidate, by which we mean that it fails one or more of the above criteria, or there is more than one candidate that fits these criteria. The photometric redshift estimates, z N IRCam , are computed using the SED fitting code LePhare, as discussed in the main text, with matched-aperture photometry of the six bands of NIRCam imaging as described in Adams et al. (2022). An entry has no z N IRCam value if the lensed source was significantly contaminated by a stellar diffraction spike, a bright cluster member galaxy or other projection effects, or if it was too faint to be detected in two or more of the NIRCam bands. The model-predicted redshifts, z model and z N IRCam are in relative agreement with the exception of five systems. The discrepancies in systems 3, 4, 12, and 15 can be attributed to a degeneracy in the SED fitting for z=6-7 and z=1.5-2.5, which comes as a result of the Lyman Break and Balmer Break both falling between F090W and F150W at these redshifts respectively. The inclusion of F115W could break this degeneracy in most cases, and may be necessary for any future surveys of the redshift distribution of background galaxies. The other outlier is image 16.1, for which z N IRCam − z model > 3σ. This system is faint and elongated, and is most likely not well represented by the 0.32 circular aperture imposed on these data, in addition to being contaminated by the ICL. We also note that that the lens models are limited by the arc systematics. Lens predictions may change as image system constraints become improve, for example by the measurement of new spectroscopic redshifts.

B. IMAGE PAIR IN SYSTEM 5
Knots A and B of system 5 are counterimages which are situated opposite the critical curve and thus likely highly magnified. As shown in Fig. 6, the image pair appears in each of the SW filters and is spatially resolved, with an angular separation of ∼ 0.16 . In SW filter F200W, the microlensed transient, Tr, is detected as well; it is bright and offset from the position of knot A, resulting in the red knot seen in Fig. 2. In the LW filters, the knots A and B drop out, leaving only the redder Tr visible as a single peak. Note-Column 1: ID; Columns 2 & 3: Right Ascension and Declination; Column 4: Redshift. For systems 1, 2 and 5 we quote the spectroscopic redshift from MUSE (Golubchik et al. 2022), and note that the NIRISS spectrum for system 1 confirms the stated redshift. Column 5: The redshift of the system as predicted from the SL model. Candidate images whose identification is not secure are marked with "?." Note that these were therefore not used in the minimization. a The first five image systems are drawn from Golubchik et al. (2022).
b These are the photometric redshift estimates using the six bands of NIRCam imaging. c These redshifts are spectroscopic values from Golubchik et al. (2022).
d These systems were not used to constrain the parametric SL model, but were confirmed in the WSLAP+ model. e Mahler et al. (2022) find a spectroscopic redshift of 1.99 for this system using MUSE, but this is not included in the model. Figure 6. The image pair A and B of system 5 across all six available JWST filters. In each filter we center on the positions of A and B with rings of radius 0.16 , which is roughly equal to their angular separation. The F090W, F150W, F200W filters exhibit two separate knots consistent with the marked positions of A and B. F200W, however, also shows the bright source Tr that is slightly offset from the position of knot A. Only a single peak is detected in each of the LW filters, which is also slightly offset from the position of knot A. The insets depict the F200W and F444W images with the diffuse arc background subtracted off and the position of knot B indicated by a "+" sign. Since knot A and knot B are counterimages of one another and must appear together, the Tr object seen in F200W and the LW filters is instead understood to be a candidate microlensing event.