11institutetext: 1 INAF-Osservatorio Astronomico di Roma, Via Frascati 33,
I-00078 Monteporzio Catone, Italy
2 Université Paris-Saclay, CNRS, Institut d’astrophysique spatiale, 91405,
Orsay, France
3 Université Paris Cité, CNRS, Astroparticule et Cosmologie, F-75013 Paris,
France
4 Instituto de Astrofísica de Canarias (IAC); Departamento de Astrofísica,
Universidad de La Laguna (ULL), E-38200, La Laguna, Tenerife, Spain
5 Instituto de Astrofísica de Canarias, Calle Vía Láctea s/n, E-38204, San
Cristóbal de La Laguna, Tenerife, Spain
6 Université Paris-Cité, 5 Rue Thomas Mann, 75013, Paris, France
7 Université PSL, Observatoire de Paris, Sorbonne Université, CNRS, LERMA,
F-75014, Paris, France
8 School of Physics and Astronomy, University of Nottingham, University Park,
Nottingham NG7 2RD, UK
9 Instituto de Física de Cantabria, Edificio Juan Jordá, Avenida de los
Castros, E-39005 Santander, Spain
10 Departamento de Física Teórica, Atómica y Óptica, Universidad de
Valladolid, 47011 Valladolid, Spain
11 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências,
Universidade de Lisboa, Tapada da Ajuda, PT-1349-018 Lisboa, Portugal
12 Department of Physics and Astronomy, Purdue University, 525 Northwestern
Ave., West Lafayette, IN 47907, USA
13 Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
University of Manchester, Oxford Road, Manchester M13 9PL, UK
14 Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-
Universität München, Scheinerstrasse 1, 81679 München, Germany
15 European Southern Observatory, Alonso de Cordova 3107, Casilla 19001,
Santiago, Chile
16 Department of Astronomy, University of Geneva, ch. dÉcogia 16, CH-1290
Versoix, Switzerland
17 Institut d’Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université,
98 bis boulevard Arago, 75014 Paris, France
18 Canada-France-Hawaii Telescope, 65-1238 Mamalahoa Hwy, Kamuela, HI 96743,
USA
19 Instituto de Matemática Estatística e Física, Universidade Federal do Rio
Grande, 96203-900, Rio Grande, RS, Brazil
20 University of Nottingham, University Park, Nottingham NG7 2RD, UK
21 Aix-Marseille Univ, CNRS, CNES, LAM, Marseille, France
22 ICRAR, M468, University of Western Australia, Crawley, WA 6009, Australia
23 INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, I-80131
Napoli, Italy
24 Institute of Cosmology and Gravitation, University of Portsmouth,
Portsmouth PO1 3FX, UK
25 Institut für Theoretische Physik, University of Heidelberg, Philosophenweg
16, 69120 Heidelberg, Germany
26 INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via
Piero Gobetti 93/3, I-40129 Bologna, Italy
27 Dipartimento di Fisica e Astronomia ”Augusto Righi” - Alma Mater Studiorum
Università di Bologna, via Piero Gobetti 93/2, I-40129 Bologna, Italy
28 INFN-Sezione di Bologna, Viale Berti Pichat 6/2, I-40127 Bologna, Italy
29 Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1,
D-85748 Garching, Germany
30 Dipartimento di Fisica, Università di Genova, Via Dodecaneso 33, I-16146,
Genova, Italy
31 INFN-Sezione di Roma Tre, Via della Vasca Navale 84, I-00146, Roma, Italy
32 Department of Physics ”E. Pancini”, University Federico II, Via Cinthia 6,
I-80126, Napoli, Italy
33 Dipartimento di Fisica, Universitá degli Studi di Torino, Via P. Giuria 1,
I-10125 Torino, Italy
34 INFN-Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy
35 INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, I-10025 Pino
Torinese (TO), Italy
36 INAF-IASF Milano, Via Alfonso Corti 12, I-20133 Milano, Italy
37 Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of
Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain
38 Port d’Informació Científica, Campus UAB, C. Albareda s/n, 08193 Bellaterra
(Barcelona), Spain
39 Institut d’Estudis Espacials de Catalunya (IEEC), Carrer Gran Capitá 2-4,
08034 Barcelona, Spain
40 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans,
s/n, 08193 Barcelona, Spain
41 INFN section of Naples, Via Cinthia 6, I-80126, Napoli, Italy
42 Dipartimento di Fisica e Astronomia ”Augusto Righi” - Alma Mater Studiorum
Universitá di Bologna, Viale Berti Pichat 6/2, I-40127 Bologna, Italy
43 INAF-Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, I-50125,
Firenze, Italy
44 Centre National d’Etudes Spatiales, Toulouse, France
45 Institut national de physique nucléaire et de physique des particules, 3
rue Michel-Ange, 75794 Paris Cédex 16, France
46 Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, UK
47 ESAC/ESA, Camino Bajo del Castillo, s/n., Urb. Villafranca del Castillo,
28692 Villanueva de la Cañada, Madrid, Spain
48 European Space Agency/ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma,
Italy
49 Univ Lyon, Univ Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822,
F-69622, Villeurbanne, France
50 Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique
Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix,
Switzerland
51 Mullard Space Science Laboratory, University College London, Holmbury St
Mary, Dorking, Surrey RH5 6NT, UK
52 Departamento de Física, Faculdade de Ciências, Universidade de Lisboa,
Edifício C8, Campo Grande, PT1749-016 Lisboa, Portugal
53 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências,
Universidade de Lisboa, Campo Grande, PT-1749-016 Lisboa, Portugal
54 Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK
55 INFN-Padova, Via Marzolo 8, I-35131 Padova, Italy
56 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, Astrophysique,
Instrumentation et Modélisation Paris-Saclay, 91191 Gif-sur-Yvette, France
57 INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, I-34143
Trieste, Italy
58 Istituto Nazionale di Astrofisica (INAF) - Osservatorio di Astrofisica e
Scienza dello Spazio (OAS), Via Gobetti 93/3, I-40127 Bologna, Italy
59 Istituto Nazionale di Fisica Nucleare, Sezione di Bologna, Via Irnerio 46,
I-40126 Bologna, Italy
60 INAF-Osservatorio Astronomico di Padova, Via dell’Osservatorio 5, I-35122
Padova, Italy
61 Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029
Blindern, N-0315 Oslo, Norway
62 Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The
Netherlands
63 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak
Grove Drive, Pasadena, CA, 91109, USA
64 von Hoerner & Sulger GmbH, SchloßPlatz 8, D-68723 Schwetzingen, Germany
65 Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark
66 Institut d’Astrophysique de Paris, 98bis Boulevard Arago, F-75014, Paris,
France
67 Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg,
Germany
68 Aix-Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France
69 Université de Genève, Département de Physique Théorique and Centre for
Astroparticle Physics, 24 quai Ernest-Ansermet, CH-1211 Genève 4, Switzerland
70 Department of Physics and Helsinki Institute of Physics, Gustaf Hällströmin
katu 2, 00014 University of Helsinki, Finland
71 NOVA optical infrared instrumentation group at ASTRON, Oude Hoogeveensedijk
4, 7991PD, Dwingeloo, The Netherlands
72 Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71,
53121 Bonn, Germany
73 Department of Physics, Institute for Computational Cosmology, Durham
University, South Road, DH1 3LE, UK
74 University of Applied Sciences and Arts of Northwestern Switzerland, School
of Engineering, 5210 Windisch, Switzerland
75 European Space Agency/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The
Netherlands
76 Department of Physics and Astronomy, University of Aarhus, Ny Munkegade
120, DK-8000 Aarhus C, Denmark
77 Centre for Astrophysics, University of Waterloo, Waterloo, Ontario N2L 3G1,
Canada
78 Department of Physics and Astronomy, University of Waterloo, Waterloo,
Ontario N2L 3G1, Canada
79 Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5,
Canada
80 Space Science Data Center, Italian Space Agency, via del Politecnico snc,
00133 Roma, Italy
81 Institute of Space Science, Bucharest, Ro-077125, Romania
82 Departamento de Astrofísica, Universidad de La Laguna, E-38206, La Laguna,
Tenerife, Spain
83 Dipartimento di Fisica e Astronomia ”G.Galilei”, Universitá di Padova, Via
Marzolo 8, I-35131 Padova, Italy
84 Departamento de Física, FCFM, Universidad de Chile, Blanco Encalada 2008,
Santiago, Chile
85 Centre for Electronic Imaging, Open University, Walton Hall, Milton Keynes,
MK7 6AA, UK
86 AIM, CEA, CNRS, Université Paris-Saclay, Université de Paris, F-91191 Gif-
sur-Yvette, France
87 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas
(CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain
88 Universidad Politécnica de Cartagena, Departamento de Electrónica y
Tecnología de Computadoras, 30202 Cartagena, Spain
89 Infrared Processing and Analysis Center, California Institute of
Technology, Pasadena, CA 91125, USA
90 INAF-Osservatorio Astronomico di Brera, Via Brera 28, I-20122 Milano, Italy
91 Junia, EPA department, F 59000 Lille, France
92 SISSA, International School for Advanced Studies, Via Bonomea 265, I-34136
Trieste TS, Italy
93 IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2,
34151 Trieste, Italy
94 INFN, Sezione di Trieste, Via Valerio 2, I-34127 Trieste TS, Italy
95 Dipartimento di Fisica e Scienze della Terra, Universitá degli Studi di
Ferrara, Via Giuseppe Saragat 1, I-44122 Ferrara, Italy
96 Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Via Giuseppe
Saragat 1, I-44122 Ferrara, Italy
97 Institut de Physique Théorique, CEA, CNRS, Université Paris-Saclay F-91191
Gif-sur-Yvette Cedex, France
98 Dipartimento di Fisica - Sezione di Astronomia, Universitá di Trieste, Via
Tiepolo 11, I-34131 Trieste, Italy
99 NASA Ames Research Center, Moffett Field, CA 94035, USA
100 INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, I-40129 Bologna,
Italy
101 INFN-Bologna, Via Irnerio 46, I-40126 Bologna, Italy
102 Institut de Recherche en Astrophysique et Planétologie (IRAP), Université
de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, F-31400 Toulouse, France
103 Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire
Lagrange, Bd de l’Observatoire, CS 34229, 06304 Nice cedex 4, France
104 Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH
Aachen University, D-52056 Aachen, Germany
105 Department of Physics & Astronomy, University of California Irvine, Irvine
CA 92697, USA
106 University of Lyon, UCB Lyon 1, CNRS/IN2P3, IUF, IP2I Lyon, France
107 INFN-Sezione di Genova, Via Dodecaneso 33, I-16146, Genova, Italy
108 INAF-Istituto di Astrofisica e Planetologia Spaziali, via del Fosso del
Cavaliere, 100, I-00100 Roma, Italy
109 Instituto de Física Teórica UAM-CSIC, Campus de Cantoblanco, E-28049
Madrid, Spain
110 Department of Physics, P.O. Box 64, 00014 University of Helsinki, Finland
111 Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical
Institute (AIRUB), German Centre for Cosmological Lensing (GCCL), 44780
Bochum, Germany
112 Department of Physics, Lancaster University, Lancaster, LA1 4YB, UK
113 Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France
114 Department of Physics and Astronomy, University College London, Gower
Street, London WC1E 6BT, UK
115 Astrophysics Group, Blackett Laboratory, Imperial College London, London
SW7 2AZ, UK
116 Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, Avenue des
Martyrs, 38000, Grenoble, France
117 Centre de Calcul de l’IN2P3, 21 avenue Pierre de Coubertin F-69627
Villeurbanne Cedex, France
118 Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 2,
I-00185 Roma, Italy
119 Zentrum für Astronomie, Universität Heidelberg, Philosophenweg 12, D-
69120 Heidelberg, Germany
120 Department of Mathematics and Physics E. De Giorgi, University of Salento,
Via per Arnesano, CP-I93, I-73100, Lecce, Italy
121 INFN, Sezione di Lecce, Via per Arnesano, CP-193, I-73100, Lecce, Italy
122 INAF-Sezione di Lecce, c/o Dipartimento Matematica e Fisica, Via per
Arnesano, I-73100, Lecce, Italy
123 Institute for Computational Science, University of Zurich,
Winterthurerstrasse 190, 8057 Zurich, Switzerland
124 Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The
University of Edinburgh, Edinburgh EH9 3FD, UK
125 Université St Joseph; Faculty of Sciences, Beirut, Lebanon
126 Department of Astrophysical Sciences, Peyton Hall, Princeton University,
Princeton, NJ 08544, USA
127 Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of
Helsinki, Helsinki, Finland
128 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700
AV Groningen, The Netherlands
129 Department of Mathematics and Physics, Roma Tre University, Via della
Vasca Navale 84, I-00146 Rome, Italy
130 Cosmic Dawn Center (DAWN)
131 Niels Bohr Institute, University of Copenhagen, Jagtvej 128, 2200
Copenhagen, Denmark
132 Departement of Physics and Astronomy, University of British Columbia,
Vancouver, BC V6T 1Z1, Canada

#  Euclid preparation. XXV. The Euclid Morphology Challenge – Towards model-
fitting photometry for billions of galaxies

Euclid Collaboration: E. Merlin0000-0001-6870-89001<EMAIL_ADDRESS>M.
Castellano0000-0001-9875-82631 H. Bretonnière0000-0001-9935-91092,3 M.
Huertas-Company0000-0002-1416-84834,5,6,7 U. Kuchner0000-0002-0035-52028 D.
Tuccillo9 F. Buitrago0000-0002-2861-981210,11 J. R.
Peterson0000-0001-5471-960912 C.J. Conselice13 F. Caro1 P.
Dimauro0000-0001-7399-28541 L. Nemani1 A. Fontana0000-0003-3820-28231 M.
Kümmel14 B. Häußler0000-0002-1857-208815 W. G. Hartley16 A. Alvarez
Ayllon0000-0002-1353-792916 E. Bertin0000-0002-3602-366417,18 P. Dubath16 F.
Ferrari0000-0002-0056-197019 L. Ferreira0000-0002-8919-079X20 R.
Gavazzi0000-0002-5540-693521,17 D. Hernández-Lang14 G.
Lucatelli0000-0002-2410-177613 A. S. G. Robotham0000-0003-0429-357922 M.
Schefer16 C. Tortora0000-0001-7958-653123 N. Aghanim2 A. Amara24 L. Amendola25
N. Auricchio0000-0003-4444-865126 M. Baldi0000-0003-4145-194327,26,28 R.
Bender0000-0001-7179-062629,14 C. Bodendorf29 E.
Branchini0000-0002-0808-690830,31 M. Brescia0000-0001-9506-568032,23 S.
Camera0000-0003-3399-357433,34,35 V. Capobianco0000-0002-3309-769235 C.
Carbone36 J. Carretero0000-0002-3130-020437,38 F. J.
Castander0000-0001-7316-457339,40 S. Cavuoti0000-0002-3787-419623,41,32 A.
Cimatti42,43 R. Cledassou0000-0002-8313-223044,45 G.
Congedo0000-0003-2508-004646 L. Conversi0000-0002-6710-847647,48 Y.
Copin0000-0002-5317-751849 L. Corcione0000-0002-6497-588135 F.
Courbin0000-0003-0758-651050 M. Cropper0000-0003-4571-946851 A. Da
Silva0000-0002-6385-160952,53 H. Degaudenzi0000-0002-5887-679916 J. Dinis52,53
M. Douspis2 F. Dubath16 C.A.J. Duncan54,13 X. Dupac47 S.
Dusini0000-0002-1128-066455 S. Farrens0000-0002-9594-938756 S. Ferriol49 M.
Frailis0000-0002-7400-213557 E. Franceschi0000-0002-0585-659126 P. Franzetti36
S. Galeotta0000-0002-3748-511557 B. Garilli0000-0001-7455-875036 B.
Gillis0000-0002-4478-127046 C. Giocoli0000-0002-9590-796158,59 A.
Grazian0000-0002-5688-066360 F. Grupp29,14 S.V.H. Haugan0000-0001-9648-726061
H. Hoekstra62 W. Holmes63 F. Hormuth64 A. Hornstrup0000-0002-3363-093665 P.
Hudelot66 K. Jahnke0000-0003-3804-213767 S. Kermiche0000-0002-0302-573568 A.
Kiessling0000-0002-2590-127363 T. Kitching0000-0002-4061-459851 R. Kohley47 M.
Kunz0000-0002-3052-739469 H. Kurki-Suonio0000-0002-4618-306370 S.
Ligori0000-0003-4172-460635 P. B. Lilje0000-0003-4324-779461 I. Lloro71 O.
Mansutti0000-0001-5758-465857 O. Marggraf0000-0001-7242-385272 K. Markovic63
F. Marulli0000-0002-8850-030327,26,28 R. Massey0000-0002-6085-378073 H.J
McCracken0000-0002-9489-776517 E. Medinaceli0000-0002-4040-778326 M.
Melchior74 M. Meneghetti0000-0003-1225-708426,28 G. Meylan50 M.
Moresco0000-0002-7616-713627,26 L. Moscardini0000-0002-3473-671627,26,28 E.
Munari0000-0002-1751-594657 S.M. Niemi75 C. Padilla0000-0001-7951-016637 S.
Paltani16 F. Pasian57 K. Pedersen76 W.J. Percival0000-0002-0644-572777,78,79
G. Polenta0000-0003-4067-919680 M. Poncet44 L. Popa81 L.
Pozzetti0000-0001-7085-041226 F. Raison29 R. Rebolo5,82 A.
Renzi0000-0001-9856-197083,55 J. Rhodes63 G. Riccio23 E.
Romelli0000-0003-3069-922257 E. Rossetti27 R. Saglia0000-0003-0378-703214,29
D. Sapone84 B. Sartoris14,57 P. Schneider72 A. Secroun0000-0003-0505-371068 G.
Seidel0000-0003-2907-353X67 C. Sirignano0000-0002-0995-714683,55 G.
Sirri0000-0003-2626-285328 J. Skottfelt0000-0003-1310-828385 J.-L.
Starck0000-0003-2177-779486 P. Tallada-Crespí87,38 A.N. Taylor46 I.
Tereno52,11 R. Toledo-Moreo0000-0002-2997-485988 I.
Tutusaus0000-0002-3199-039969 L. Valenziano0000-0002-1170-010426,28 T.
Vassallo0000-0001-6512-635857 Y. Wang0000-0002-4749-298489 J.
Weller0000-0002-8282-201014,29 A. Zacchei0000-0003-0396-119257 G.
Zamorani0000-0002-2318-301X26 J. Zoubian68 S. Andreon0000-0002-2041-878490 S.
Bardelli0000-0002-8900-029826 A. Boucaud0000-0001-7387-26333 C. Colodro-Conde4
D. Di Ferdinando28 J. Graciá-Carpio29 V. Lindholm70 N.
Mauri0000-0001-8196-154842,28 S. Mei0000-0002-2849-559X3 C. Neissner37 V.
Scottez66,91 A. Tramacere0000-0002-8186-379316 E. Zucca0000-0002-5845-813226
C. Baccigalupi0000-0002-8211-163092,93,57,94 A. Balaguera-Antolínez4,82 M.
Ballardini95,96,26 F. Bernardeau97 A. Biviano93,57 S.
Borgani0000-0001-6151-643998,93,57,94 A.S. Borlaff0000-0003-3249-443199 C.
Burigana95,100,101 R. Cabanac0000-0001-6679-2600102 A. Cappi26,103 C.S.
Carvalho11 S. Casas0000-0002-4751-5138104 G.
Castignani0000-0001-6831-068727,26 A.R. Cooray105 J. Coupon16 H.M.
Courtois0000-0003-0509-1776106 O. Cucciati0000-0002-9336-755126 S. Davini107
G. De Lucia0000-0002-6220-910457 G. Desprez16 J.A. Escartin29 S.
Escoffier0000-0002-2847-749868 M. Farina108 K. Ganga0000-0001-8159-82083 J.
Garcia-Bellido0000-0002-9370-8360109 K. George0000-0002-1734-845514 G.
Gozaliasl0000-0002-0236-919X110 H. Hildebrandt0000-0002-9814-3338111 I.
Hook0000-0002-2960-978X112 O. Ilbert21 S. Ilić113,44,102 B. Joachimi114 V.
Kansal86 E. Keihanen0000-0003-1804-771570 C.C. Kirkpatrick70 A.
Loureiro0000-0002-4371-0876115,114,46 J. Macias-Perez0000-0002-5385-2763116 M.
Magliocchetti108 G. Mainetti117 R. Maoli118,1 S. Marcin74 M.
Martinelli0000-0002-6943-77321 N. Martinet0000-0003-2786-779021 S. Matthew46
M. Maturi25,119 R.B. Metcalf0000-0003-3167-257427,26 P.
Monaco0000-0003-2083-756498,93,57,94 G. Morgante26 S.
Nadathur0000-0001-9070-310224 A.A. Nucita120,121,122 L. Patrizii28 V. Popa81
C. Porciani0000-0002-7797-250872 D. Potter0000-0002-0757-5195123 A.
Pourtsidou0000-0001-9110-555046,124 M. Pöntinen0000-0001-5442-2530110 P.
Reimberg66 A.G. Sánchez0000-0003-1198-831X29 Z.
Sakr0000-0002-4823-3757102,25,125 M. Schirmer0000-0003-2568-999467 M.
Sereno0000-0003-0302-032526,28 J. Stadel0000-0001-7565-8622123 R. Teyssier126
C. Valieri28 J. Valiviita0000-0001-6225-3693127 S.E. van
Mierlo0000-0001-8289-2863128 A. Veropalumbo0000-0003-2387-1194129 M.
Viel0000-0002-2642-570793,57,94,92 J. R. Weaver0000-0003-1614-196X130,131 D.
Scott0000-0002-6878-9840132

(September 2022)

The European Space Agency’s Euclid mission will provide high-quality imaging
for about $1.5$ billion galaxies. A software pipeline to automatically process
and analyse such a huge amount of data in real time is being developed by the
Science Ground Segment of the Euclid Consortium; this pipeline will include a
model-fitting algorithm, which will provide photometric and morphological
estimates of paramount importance for the core science goals of the mission
and for legacy science. The Euclid Morphology Challenge is a comparative
investigation of the performance of five model-fitting software packages on
simulated Euclid data, aimed at providing the baseline to identify the best
suited algorithm to be implemented in the pipeline. In this paper we describe
the simulated data set, and we discuss the photometry results. A companion
paper (Euclid Collaboration: Bretonnière et al. 2022) is focused on the
structural and morphological estimates. We created mock Euclid images
simulating five fields of view of 0.48 deg2 each in the
$I_{\scriptscriptstyle\rm E}$ band of the VIS instrument, containing a total
of about one and a half million galaxies (of which 350 000 have nominal
signal-to-noise ratio above $5$), each with three realisations of galaxy
profiles (single and double Sérsic, and ‘realistic’ profiles obtained with a
neural network); for one of the fields in the double Sérsic realisation, we
also simulated images for the three near-infrared $Y_{\scriptscriptstyle\rm
E}$, $J_{\scriptscriptstyle\rm E}$ and $H_{\scriptscriptstyle\rm E}$ bands of
the NISP-P instrument, and five Rubin/LSST optical complementary bands ($u$,
$g$, $r$, $i$, and $z$), which together form a typical data set for a Euclid
observation. The images were simulated at the expected Euclid Wide Survey
depths. To analyse the results we created diagnostic plots and defined metrics
to take into account the completeness of the provided catalogues, and the
median biases, dispersions, and outlier fractions of their measured flux
distributions. Five model-fitting software packages (DeepLeGATo, Galapagos-2,
Morfometryka, ProFit, and SourceXtractor++) were compared, all typically
providing good results. Of the differences among them, some were at least
partly due to the distinct strategies adopted to perform the measurements. In
the best case scenario, the median bias of the measured fluxes in the
analytical profile realisations is below 1% at signal-to-noise ratio above 5
in $I_{\scriptscriptstyle\rm E}$, and above 10 in all the other bands; the
dispersion of the distribution is typically comparable to the theoretically
expected one, with a small fraction of catastrophic outliers. However, we can
expect that real observations will prove to be more demanding, since the
results were found to be less accurate on the most realistic realisation. We
conclude that existing model-fitting software can provide accurate photometric
measurements on Euclid data sets. The results of the challenge are fully
available and reproducible through an online plotting tool.

###### Key Words.:

Galaxies: structure – Galaxies: evolution – Cosmology: observations

##  1 Introduction

The European Space Agency’s Euclid mission (Laureijs et al., 2011, mission
RedBook), due to start operations in 2023, is designed to provide accurate
photometric, spectroscopic and morphological data (in particular cosmic shear
and clustering distributions) for billions of galaxies across 15 000 deg2 of
sky, using them as tracers to study the properties of the dark components of
the Universe.

To this end, a processing pipeline is being assembled by the Science Ground
Segment, a team that is in charge of releasing the data to the community. This
pipeline is ready to ingest, process and analyse the raw imaging data from the
satellite on a daily basis; optical data from external ground-based
instruments (Rubin/LSST, DECAM, CHFT, Pan-STARRS, OmegaCAM, Subaru; see Euclid
Collaboration: Scaramella et al., 2022) will also be used to complement the
optical and near-infrared images obtained by the two satellite photometers VIS
(observing in $I_{\scriptscriptstyle\rm E}$, a broad optical band; see Cropper
et al., 2016) and NISP-P (observing in the three near-infrared – NIR – bands
$Y_{\scriptscriptstyle\rm E}$, $J_{\scriptscriptstyle\rm E}$, and
$H_{\scriptscriptstyle\rm E}$; Maciaszek et al., 2016; Euclid Collaboration:
Schirmer et al., 2022), allowing for high-quality photometric redshift
estimates. The final step of the image analysis pipeline will produce a global
catalogue containing all the astrometric, photometric and morphological
information about each source detected in the $I_{\scriptscriptstyle\rm E}$
images (plus an additional sample of NIR-detected sources). This catalogue
will then be exploited for scientific use by the Euclid Collaboration, and it
will also be released to the community for legacy use.

The pipeline currently implements two photometric techniques (aperture and
template-fitting, performed with Euclid-specific versions of two public
software tools, a-phot and t-phot respectively; see Merlin et al., 2015, 2016,
2019), and a module to estimate so-called CAS morphological parameters
(Concentration/Asymmetry/Smoothness: non-parametric morphological features
that can be used to distinguish between discs, ellipticals, compact, diffuse,
symmetric/asymmetric or clumpy objects by means of a dimensional reduction,
see Conselice, 2003). However, the pipeline is foreseen to also include a
profile model-fitting algorithm. The Euclid Morphology Challenge (EMC) was
organized with the aim of analysing and comparing the performance of various
model-fitting software tools on Euclid data, in order to establish the
foundations for choosing the tools that will be integrated into the official
processing pipeline. The final choice will be driven by many factors,
including computational performance, robustness of the algorithm, and
compatibility with the current version of the pipeline; however, the accuracy
of the parameter estimates will of course be the main driver. Therefore,
assessing the performance of the different software packages on simulated
data, for which the ground truth is known, is a necessary and fundamental step
for a sound selection. Eight development teams of model-fitting software
packages were invited to participate to the challenge, and five provided at
least partial results.

In this paper we present the data set created for the EMC, and we discuss the
results concerning photometry. In fact, albeit not being the central focus of
the challenge, flux measurements obtained via model-fitting techniques will
have great relevance, providing a crucial complement to the more
straightforward methods already included in the pipeline. A companion paper is
dedicated to the analysis of such morphological estimates (Euclid
Collaboration: Bretonnière et al. 2022, EMC2022b hereafter).

This paper is structured as follows. Section 2 describes the technique used to
create the simulated data set for the Challenge, with some technical details
given in Appendix A. In Sect. 3 we briefly present the software tools taking
part in the Challenge, and in Sect. 4 we describe the methods used to analyse
and rank the data provided by the participants. The results are then presented
in Sect. 5, where we investigate the general accuracy of the photometric
measurements, and the reliability of the estimated uncertainty budgets, with a
further focus on each software package’s performance given in Appendix B.
Finally, Sect. 6 presents a summary of the work and provides conclusions.

All magnitudes are given in the AB system.

##  2 Simulating the Euclid universe

Simulated data sets are being produced and used by the Euclid Science Ground
Segment to test the full processing pipeline from image reduction to data
analysis. These simulations consist of raw single exposures, including
observational features and defects, and they must be processed and stacked to
reach the nominal depth and be ready for scientific analysis, with background
light and defects removed. To simplify this complex procedure, and to have all
details under control, for the EMC we decided to produce a tailored data set,
directly simulating background-subtracted images at the expected nominal
depths of the final stacked mosaics in all bands. With this approach, we were
also free to try different options, producing simulations with single and
double Sérsic analytical profiles, and also with realistic morphologies. In
this section we explain the procedure we followed to obtain all these
simulated data sets.

###  2.1 Catalogues and images creation

We started by creating mock cosmological catalogues with the code Egg (v1.3.1,
Schreiber et al., 2017). Egg uses the statistical distributions of real
galaxies as detected and classified in the five CANDELS fields (Grogin et al.,
2011; Koekemoer et al., 2011) to build a simulated catalogue of a patch of the
sky, complete with the properties of the objects as observed by a chosen set
of pass-band filters, with a chosen pixel resolution, and to a chosen limiting
magnitude. We refer the reader to the paper describing the code for a detailed
description of its workflow; here we provide a short summary. Egg draws
redshifts and stellar masses from observed galaxy stellar mass ($M_{*}$)
functions, and subsequently attributes a star-formation rate (SFR) to each
galaxy from the observed SFR–$M_{*}$ main sequence; dust attenuation, optical
colours and simple disc-plus-bulge morphologies are obtained from empirical
relations established from the high-quality Hubble and Herschel observations
of the CANDELS fields. Random scatter is introduced in each step to reproduce
the observed distributions of each parameter. Finally, based on these
observables, a suitable panchromatic spectral energy distribution (SED) is
selected for each galaxy and synthetic photometry is produced by integrating
the redshifted SED over the chosen broad-band filters. The galaxies are
created as two-component objects, with a bulge and a disc both described by a
Sérsic (1968) profile,

$I(r)\propto\exp{\left[-b_{n}(r/r_{\rm{e}})^{1/n}\right]},\ $ (1)

where Sérsic indices for the bulge and disc components are $n_{\rm bulge}=4$
and $n_{\rm disc}=1$. The output catalogue contains the physical and observed
properties of the galaxies within a field of view (FoV) corresponding to the
chosen area; the objects are placed at random positions with a fixed angular
two-point correlation function, neglecting large-scale clustering beyond
$None$ (i.e., beyond $\sim$1 Mpc at $z>0.5$).

Figure 1: Filter transmission curves used to simulate the images (normalised
to arbitrary units) in this work. The four filters $I_{\scriptscriptstyle\rm
E}$, $Y_{\scriptscriptstyle\rm E}$, $J_{\scriptscriptstyle\rm E}$, and
$H_{\scriptscriptstyle\rm E}$ form the Euclid set, while $u$, $g$, $r$, $i$,
and $z$ are external complementary filters (from LSST/VRO). The NIR curves are
early estimates of the actual transmission functions, which are described in
their latest updated version in Euclid Collaboration: Schirmer et al. (2022);
however, since we only use these data to compute the input fluxes, the
differences are not relevant for the present work.

We created five catalogues, each one with a size of $0.482$ deg2 ($None$ per
side, which is comparable to the typical area on which each single photometric
catalogue will be extracted from real data), with limiting magnitude
$I_{\scriptscriptstyle\rm E}=27.1$ (the nominal 1$\sigma$ limit as given in
the mission RedBook). The total number of simulated galaxies was about $1.5$
million. We then used an Egg built-in script to obtain all the observational
properties of the sources. In particular, for each galaxy the following
parameters are given: position of the centroid in pixels; total flux in the
simulated band; bulge-to-total flux ratio (${\rm b}/{\rm t}$); scale length of
the bulge and of the disc (defined as the radius at which the component is a
factor of $e$ less bright than it is at its center); axis ratio for both
components; and position angle for both components. For the first of the five
fields (F0), we produced nine lists, to include a full multi-wavelength
realisation of a Euclidian sky patch: one for each of the four Euclid bands
$I_{\scriptscriptstyle\rm E}$, $Y_{\scriptscriptstyle\rm E}$,
$J_{\scriptscriptstyle\rm E}$ and $H_{\scriptscriptstyle\rm E}$, plus five for
the Rubin/LSST bands $u$, $g$, $r$, $i$ and $z$. The filter transmission
curves are shown in Fig. 1. For the other four fields (F1–4) we only produced
the $I_{\scriptscriptstyle\rm E}$ list, since the main purpose of these
simulations is the morphological analysis, which with real data will mostly be
performed on the $I_{\scriptscriptstyle\rm E}$ images, given that it will be
the band with the highest resolution and depth. We point out that in the
multi-band realisation the morphological parameters do not change across the
spectrum, while total fluxes and ${\rm b}/{\rm t}$ do; this information was
not explicitly shared with the participants.

We fed these catalogues to GalSim (Rowe et al., 2015, v2.2.1), a Python
package that produces simulated astronomical images (it is also used for the
official Euclid simulations). We created the images at their expected native
pixel scale: $None$ for $I_{\scriptscriptstyle\rm E}$ ($25\,000\times 25\,000$
pixels); $None$ for NIR bands; and $None$ for LSST bands. After the procedure
described in the following paragraphs, we resampled all the images to the
$I_{\scriptscriptstyle\rm E}$ pixel scale, since this is the procedure that
will be followed in the real pipeline. To simulate the effects of point spread
functions (PSFs), we provided GalSim with the Euclid Mission Database models
for the $I_{\scriptscriptstyle\rm E}$ and NIR bands (as provided by the
corresponding Euclid working groups; both of them are over-sampled by a factor
of 6), while for LSST we provided custom simulated PSFs created using PhoSim
(Peterson et al., 2015), at the expected observed pixel scale of $None$ (with
no over-sampling). The approximate FWHMs of these PSFs are $None$
($I_{\scriptscriptstyle\rm E}$), $None$ (NIR), and $None$–$None$ (LSST).

We produced noiseless galactic profiles with GalSim, with pixel values in
$\mathrm{\SIUnitSymbolMicro}$Jy/pixel. We simulated three different sets of
images, all having identical sets of coordinates and total fluxes of the
galaxies, which we describe in the following; see also EMC2022b for further
details on this.

  * •

Double Sérsic profiles (DS), directly using the standard output of Egg, which
consists of a catalogue formatted to be used with SkyMaker (Bertin, 2009). In
particular, this means that the dimensions of the objects are given as scale
lengths. On the contrary, GalSim requires half-light radii; the two values
coincide for the bulges, while the conversion factor is 1.678 for the discs
(see e.g. Graham & Driver, 2005), so we applied this correction before
simulating the images. Also, GalSim requires that the fluxes of the two
components are given separately, while Egg outputs a total magnitude and a
${\rm b}/{\rm t}$; therefore, to assign a flux to each component we simply
used the relations $f_{\rm bulge}={\rm b}/{\rm t}\;f_{\rm tot}$ and $f_{\rm
disc}=(1-{\rm b}/{\rm t})\;f_{\rm tot}$, where $f_{\rm tot}=10^{-0.4(m-{\rm
ZP})}$ (where $m$ is the magnitude of the sources as given in the Egg
catalogue and ZP is the zero-point of the image, see Sect. 2.1.1).

  * •

Single Sérsic profiles (SS), in which galaxies are modeled with a single
Sérsic index, defined using the ${\rm b}/{\rm t}$ values from the Egg
catalogue as $n_{\rm total}=(1-{\rm b}/{\rm t})\,n_{\rm disc}+{\rm b}/{\rm
t}\;n_{\rm bulge}=3\;{\rm b}/{\rm t}+1$. To compute the single effective
radius $r_{\rm e,tot}$ from the two values given in the 2-component catalogue,
we used the following formula, calibrated empirically to obtain a good visual
match between the two realisations: $r_{\rm tot}=[{\rm b}/{\rm
t}\;r_{\mathrm{e,b}}]^{\alpha}+[(1-{\rm b}/{\rm
t})\,r_{\mathrm{e,d}}]^{\alpha}$, where $\alpha=0.8$ if
$r_{\mathrm{e,b}}<r_{\mathrm{e,d}}$ (98% of the cases), and $\alpha=2.0$
otherwise. Finally, position angles and axis ratios already had the same
values for bulges and discs, so we simply kept them unchanged.

  * •

Realistic morphologies (RM), in which galaxy stamps are created by means of a
neural network using a variational auto-encoder trained on observed COSMOS
galaxies, as described in full detail in Lanusse et al. (2021) and Euclid
Collaboration: Bretonnière et al. (2022); each simulated galaxy mimics the
properties of its corresponding analytical realisation. In this data set, the
biggest and brightest objects ($r_{\rm e}>$$$, $I_{\scriptscriptstyle\rm
E}<20.5$) are not simulated due to technical limitations; the list of excluded
sources was provided to the participants, and accounts for approximately $1\%$
of the total simulated galaxies. Also, the position angles of the galaxies are
not constrained to be close to those of the Egg catalogue (and therefore they
were not considered in the final analysis of the results). We point out that
this is the first time that such a demanding test has been performed: the
codes must provide an analytical fit on non-analytical shapes for which a
ground-truth value is known. This is inherently a very challenging task.
Moreover, the method used to create the images is not perfect. The
conditioning of the latent space with galaxy morphology is not always exact,
and this can introduce a systematic bias with respect to the input values (see
the discussion in Euclid Collaboration: Bretonnière et al., 2022), although
the consistency is fully guaranteed in a statistical sense; some level of
scatter remains on a object-by-object basis, meaning that the comparison with
the input catalogue must be taken with caution. For more details, see
EMC2022b.

Figure 2: Distribution of redshifts as a function of $I_{\scriptscriptstyle\rm
E}$ magnitude (input values) in F0. The orange line is the running mean, the
vertical dashed line the 5$\sigma$ limit. Figure 3: Distribution of bulge-to-
total ratios as a function of $I_{\scriptscriptstyle\rm E}$ magnitude in the
input catalogue of F0; the colours encode the global scale length, defined
here as $r={\rm b}/{\rm t}\;r_{\mathrm{e,bulge}}+(1-{\rm b}/{\rm
t})\,r_{\mathrm{e,disc}}$. Figure 4: Distributions of various input parameters
for the $I_{\scriptscriptstyle\rm E}$ band, in all five simulated fields, for
all the separate components and realisations in the Challenge, as described in
the legend and in the labels of the panels. For comparison, we also show the
corresponding parameter distribution in COSMOS.

Figure 2 shows the redshift distribution as a function of
$I_{\scriptscriptstyle\rm E}$ input magnitude. Egg outputs the redshift of
each simulated galaxy, and although this information is not explicitly used in
the present work, it might nevertheless be useful to have an idea of the
global distribution; indeed, most of the analysis and figures will be
presented as a function of the input $I_{\scriptscriptstyle\rm E}$ magnitude,
which correlates with redshift. For example, looking at the plot one can see
that galaxies at $z=3$ typically begin to be detectable at
$I_{\scriptscriptstyle\rm E}=23$. Figure 3 shows the distribution of simulated
galaxies in the magnitude-size-${\rm b}/{\rm t}$ space for the same field
(note how there is a non-negligible fraction of bulge-only objects with ${\rm
b}/{\rm t}=1$). Finally, in Fig. 4 we show the distributions of various input
parameters for the $I_{\scriptscriptstyle\rm E}$ band components and
realisations (magnitude, effective radius, axis ratio, Sérsic index, bulge-to-
total ratio and redshift), for the full data set (the five fields), down to
the nominal 5$\sigma$ limit $I_{\scriptscriptstyle\rm E}=25.35$ (Laureijs et
al., 2011); we also show the COSMOS distributions (Mandelbaum et al., 2012)
for reference. As expected, simulations and observations agree remarkably
well, with the exception of the ${\rm b}/{\rm t}$ distribution, which is more
skewed towards disc galaxies in the Egg catalogue. This might have some impact
on the analysis of the results (see Sect. 5).

By construction, we could not simulate irregulars, which are estimated to
constitute less than $10\%$ of the galaxies at $z<1$, but up to $70\%$ at
$z=3$ (e.g. Huertas-Company et al., 2015, , although recent preliminary
results from the James Webb Space Telescope seem to indicate a lower number).
This is an obvious but unavoidable limitation of this work. The RM realisation
can provide a hint about how model-fitting codes can deal with non-analytical
shapes.

We then added a field of stars, to include the effects of their presence as
contaminants in the fitting procedures. To obtain a realistic distribution, we
took their celestial coordinates (converted to pixel positions) from one FoV
of the official Euclid simulations used for Scientific
Challenges,111Scientific Challenges are official Euclidean benchmark tests
performed to check and validate the progress of the work in preparation for
the launch of the satellite and simply placed PSF stamps at the positions of
the sources, scaling their flux to match the catalogue magnitudes. We excluded
very bright stars ($I_{\scriptscriptstyle\rm E}<15$), in order to avoid that
large regions of the simulations were affected by their presence, and also
because – given the limited extension of the PSF stamps – they would saturate
creating artificial defects on the images. The fraction of pixels
significantly contaminated by stellar light (that is, where the surface
brightness from stellar light is more than the 1$\sigma$ surface brightness
per pixel) is approximately $1\%$ in the $I_{\scriptscriptstyle\rm E}$ images,
$5\%$ in the NIR bands, and from $1\%$ to $7\%$ in LSST bands ($u$ and $z$
respectively).

####  2.1.1 Observational noise

Band | $m_{\rm lim}$ | SBbkg |  $t_{\rm exp}$ [s]
---|---|---|---
$I_{\scriptscriptstyle\rm E}$ | 24.6 | 22.33 | 4 $\times$ 590
$Y_{\scriptscriptstyle\rm E}$ | 23.0 | 22.10 | 4 $\times$ 88
$J_{\scriptscriptstyle\rm E}$ | 23.0 | 22.11 | 4 $\times$ 90
$H_{\scriptscriptstyle\rm E}$ | 23.0 | 22.28 | 4 $\times$ 54
$u$ | 23.6 | 22.70 | 150
$g$ | 24.5 | 22.00 | 150
$r$ | 23.9 | 20.80 | 150
$i$ | 23.6 | 20.30 | 150
$z$ | 23.4 | 19.40 | 150
Table 1: Parameters used to simulate the images. $\mbox{m}_{\rm lim}$ is the
10$\sigma$ limiting magnitude within a 2$\mathrm{\SIUnitSymbolArcsecond}$
aperture, SBbkg is the background surface brightness, and $t_{\rm exp}$ is the
total exposure time of the final mosaic (these are not updated to the latest
estimates of the actual in-flight values). See text for more details.

Once the seed images containing the sources were produced, we proceeded to add
simulated observational noise. First of all we replaced the smooth analytical
profiles with stochastic realisations from a Poissonian distribution, to
simulate the effects of photon shot noise. We used the Python module
scipy.stats.poisson.rvs for this purpose. We paid particular attention to keep
the units of the images always consistent during the whole process: we first
converted the noiseless images from $\mathrm{\SIUnitSymbolMicro}$Jy/pixel to
observational units, using the correct image observational ZP at 1 second, to
obtain consistent Poissonian realisations, which depend on the total exposure
time. Only after this step did we convert the images back to
$\mathrm{\SIUnitSymbolMicro}$Jy. To calculate the ZPs we followed the
procedure described in Euclid Collaboration: Martinet et al. (2019), which we
summarize in Appendix A.1. The same method was also used to produce an empty
sky map containing only Gaussian noise, simulating the observational
background at the desired depth. Since the images were simulated with zero
background light, the Gaussian noise must have zero mean, and the standard
deviation of the pixel values defines the depth of the final simulated image.

The values used in this procedure are summarized in Table 1. The exposure
times used in our simulations the for $I_{\scriptscriptstyle\rm E}$ and for
the NIR bands were taken from Laureijs et al. (2011), and it is worth pointing
out that the actual in-flight values will be slightly different ($4\times 560$
seconds for $I_{\scriptscriptstyle\rm E}$ and $4\times 88$ seconds for all NIR
bands). The exposure times for the LSST bands were estimated from simulated
data created for one of the internal validation Scientific Challenges, and are
representative of an early release (LSST final data after ten years of
observations will of course be much deeper). The limiting magnitudes and
background surface brightness values, which are needed in the computations,
are set to be consistent with the expected values for the Euclid Wide Survey,
and they were taken from a dedicated study by J. C. Cuillandre (priv. comm.);
these too are slightly different from the current best estimated values, which
can be found in Euclid Collaboration: Scaramella et al. (2022). These small
inconsistencies are due to the fact that the present work began before the
latest estimates had been made available; however, they have negligible impact
on the scientific results of the EMC. It is worth pointing out that the images
were produced with homogeneous noise levels, i.e we did not simulate regions
of different depths. Finally, we summed each image containing the Poissonian
realisations of the galaxies and stars with the corresponding ‘empty’ Gaussian
sky noise.

####  2.1.2 Rebinning

As mentioned, we simulated all the images at their native pixel scales ($None$
for $I_{\scriptscriptstyle\rm E}$, $None$ for NIR, and $None$ for LSST), and
we then rebinned F0 NIR and LSST images to the $I_{\scriptscriptstyle\rm E}$
pixel scale using Swarp (Bertin et al., 2002). This is consistent with the
real pipeline workflow, with some differences: in the pipeline, a newer
version of the software named Swarp++ is used, and single-exposure images are
combined to create the mosaics, allowing information to be gained in the
process. We also rebinned the PSFs accordingly. All the rebinning processes
were performed using the BILINEAR interpolation mode. We note that this
resampling procedure introduces artifacts in the noise map (in particular,
pixel correlations) that alter the apparent signal-to-noise ratio (S/N) of the
map, so the actual uncertainties of the measurements must be computed using a
dedicated RMS map, which we discuss next.

###  2.2 RMS maps

It is common practice to assign uncertainties to the measurements performed on
scientific images by means of a weight or an RMS map, which is often obtained
from first principles during the data reduction chain. When this is not
possible, it can be easily determined, at least to a first approximation, by
measuring the RMS of the pixel values in ‘empty’ regions of the (non-rebinned)
science frame – although such a measurement only provides information on the
noise due to the unresolved sky background. For the sake of the EMC goals, we
wanted to factorise out any possible source of complication, and therefore
ready-to-use RMS maps were provided to the participants, along with the
scientific images. This is also again consistent with the pipeline
architecture. The procedure to build the RMS map is descibed in Appendix A.2.

The RMS maps were produced at the native pixel scales of the scientific
images, and we checked that the pre-resampling S/Ns are consistent with the
expected values. This is shown in Fig. 5, where we plot the S/N estimated for
each simulated source and check that it is equal to the expected value at the
limiting magnitude (star symbols); for this test we used a-phot, forcing the
measurements within 2$\mathrm{\SIUnitSymbolArcsecond}$ apertures (as for the
definition of S/N adopted in this work) at the true input positions of the
sources. Note that some distributions of points are overlapping (the three NIR
bands have the same expected depth, and so do two of the LSST bands). The
overall agreement with the expected values is very accurate. Finally, we
proceeded to resample the maps along with the scientific images, again using
Swarp, checking that the S/N values of the resampled images are correct when
the RMS maps are used to estimate uncertanties.

Figure 5: Signal-to-noise ratios of the simulated images at the native pixel
scales, measured with forced photometry in 2$\mathrm{\SIUnitSymbolArcsecond}$
diameter apertures at the input positions of the sources. Star symbols show
the expected values, which are also reported in the legend. Figure 6: A small
region of the DS F0 realisation in the nine bands; all are shown with the same
colour scale. Figure 7: A small crop of the three realisations of
$I_{\scriptscriptstyle\rm E}$ images. Left to right: DS, SS, RM. Notice that
in RM the orientations of the galaxies are not the same as in the other two
realisations.

###  2.3 Challenge set-up

The scientific and RMS images were finally uploaded to a private online
repository for the participants to download, together with lists containing
the IDs and input positions of the centroids of the simulated objects down to
various nominal S/N levels ($100$, $10$, $5$, $1\sigma$). This was done to
factor out possible inaccuracies in object detection and deblending, so that
the challenge could be actually focused on the accuracy of fitting photometry
and morphology without adding any further possible source of error. We point
out that to obtain the lists we simply applied different cuts to the input
$I_{\scriptscriptstyle\rm E}$ magnitudes, so they must be considered as coarse
reference levels rather than accurate estimations of detection significance.

As a visual example, Fig. 6 shows a small crop of the DS F0 images in all nine
bands, while Fig. 7 shows a small crop of the DS, SS and RM F0
$I_{\scriptscriptstyle\rm E}$ images. A sample of the ground-truth input
values of the simulated galaxies was also provided to the participants (in
particular, a small portion of the multi-band F0 data set and the whole F4
$I_{\scriptscriptstyle\rm E}$ data set), allowing for a check that their
procedures were reasonably correct without evident errors.

The output requested from each participant consisted of the estimates of (i)
flux, (ii) Sérsic index (in the SS and RM realisations), (iii) half-light
semi-major axis, (iv) axis ratio, and (v) position angle, each with a
corresponding 1$\sigma$ uncertainty, for each component of the simulated
objects. In particular, while the SS and RM realisations only required a
single fit with free $n$, for the DS realisation we asked for two estimates,
namely: one fit with fixed indices ($n_{\rm disc}=1$ and $n_{\rm bulge}=4$,
consistent with the way the images were simulated); and another with $n_{\rm
disc}=1$ and $n_{\rm bulge}$ left free to vary.

As mentioned, the output estimates were required for the objects belonging to
the list including only ${\rm S/N}>5$ (in $I_{\scriptscriptstyle\rm E}$)
sources. Analysis of objects at lower S/N were explicitly mentioned as an
optional output that would not influence the final comparison among the
software packages.

##  3 Model-fitting software packages

Eight development teams of model-fitting software packages were invited to
participate to the Challenge; of them, five (DeepLeGATo, Galapagos-2,
Morfometryka, ProFit, and SourceXtractor++) provided at least partial results.
All but one are based on parametric methods, using functional forms to fit the
observed light distributions, the exception being DeepLeGATo, which exploits a
convolutional neural network (CNN).

Here we briefly summarize the basic properties and features of these five
software packages, and point out a few important details about the procedures
each one followed. It is instructive and important to notice how various
subtleties in the interpretation of the requests, usage of the provided input
data sets, and processing methods used by each participant led to differences
in the format and accuracy of the provided outputs.

###  3.1 DeepLeGATo

DeepLeGATo (Tuccillo et al., 2018) is a software package for estimating galaxy
structure based on a supervised deep-learning approach. The code uses CNNs to
perform a simple regression between an image centered on a given galaxy and
its structural parameters, providing results in very short times (see Appendix
B.1). In the version used in this work, the training was performed with images
of fixed size (128$\times$128 pixels) independently of the galaxy effective
radius; this likely caused sub-optimal performance on the largest and
brightest galaxies. The images used for training were not the ones provided
with the EMC data-set; instead, they were idealised analytic 1- or 2-component
Sérsic profiles, convolved with the PSF of the corresponding band, to which
noise having similar properties to that in the EMC data was added. A different
training was performed for each structural parameter with slightly different
architectures, which are variations of standard CNNs (see Tuccillo et al.,
2018, for more details). The fits were performed at the positions of the
sources as given in the input files.

Because the loss function used for training is a standard mean square error,
only point estimates were provided; therefore in the implementation used for
this work the uncertainty budget (i.e. the uncertainties on the estimates) was
not computed. Noticeably, the 5$\sigma$ input source lists were subdivided
into S/N bins using the other input catalogs (10 and 100$\sigma$), and
different parameters were used for the fits in different bins; this can be
seen in the distributions of the points in the plots (see Appendix B). While
this makes the version of the software used in the Challenge not directly
suitable for an implementation in the Euclid pipeline, because the S/N of real
sources cannot be known a priori, it is worth pointing out that the code
remains under development, and more recent releases work without the need for
this fine-tuning of the parameters for different input data. Only
$I_{\scriptscriptstyle\rm E}$ SS and DS fits were provided.

###  3.2 Galapagos-2

Galapagos-2 (Häußler et al., 2022)222https://github.com/MegaMorph/galapagos is
an updated and enhanced version of Galapagos (Barden et al., 2012). It
provides a wrapper around either Galfit (Peng et al., 2002, 2010) for single-
band fits, or GalfitM (Häußler et al., 2013) for either single- or multi-band
fits, and it is specifically designed to carry out fully automated fitting on
all objects in a large survey. Starting from the input images and a simple
setup file, it employs SExtractor (Bertin & Arnouts, 1996) for object
detection, and then uses this information to automatically set up the fits.
The postage stamp size used for each fit/object depends on the estimated size
of the object, set up by the user. Using some limited input from the setup,
e.g. enlargement factors to conservatively increase the size and shape of the
object estimated by SExtractor, Galapagos-2 takes care of neighbouring objects
(deblending and fitting of bright nearby objects, masking of fainter and more
distant objects), estimates the sky background level with a sophisticated and
robust scheme (see Barden et al., 2012, for details), and sets up the fit
using SExtractor values as initial guesses to run the fitting algorithm
GalfitM for all objects, starting from the brightest objects and using the PSF
provided. Once an object has been fit, it is merely subtracted from the fits
of nearby objects. This significantly speeds up the process overall, as these
bright, large objects take the longest to fit, but this only needs to be done
once. In a fully automated pipeline, it then reads out the result and provides
one final catalogue, which contains fitting information for all bands. GalfitM
itself uses a Levenberg–Marquardt minimisation to derive the best-fit
parameters and uncertainties. In the multi-band realisations of the EMC, all
bands are connected via physically reasonable polynomials and fitted
simultaneously, to reduce the degrees of freedom of the fit and make full use
of the multi-wavelength information. The software requires one additional
input image compared to what was provided, namely a weight image to flag bad
pixels. Since no bad pixels were in the data, this was trivially created as a
uniform image of the correct size.

All the requested outputs were provided; however, for DS F0 only a
simultaneous fit was run, therefore including $I_{\scriptscriptstyle\rm E}$ in
the multi-band fitting process – in other words, there is no isolated fit for
$I_{\scriptscriptstyle\rm E}$ DS F0. This causes the $I_{\scriptscriptstyle\rm
E}$ fit for F0 to be substantially different from that of other fields; for
this reason, it was decided not to include F0 in the analysis of
$I_{\scriptscriptstyle\rm E}$-only results for all codes (see below).

###  3.3 Morfometryka

Morfometryka (Ferrari et al., 2015), written in Python, was primarily designed
to measure non-parametric morphological quantities, but as a bonus it performs
single Sérsic model fitting. The software takes as input a galaxy stamp (plus
the PSF model), estimates the background with an iterative algorithm, segments
the sources and defines the target. Then, it filters out external sources
using the code GalClean (de Albernaz Ferreira & Ferrari, 2018). From the
segmented region it calculates basic geometrical parameters (e.g. centroid,
position angle, axis ratio) using light-profile moments. Then it performs
photometry, measuring fluxes within ellipses with the aforementioned
parameters (contextually masking out point sources over the ellipse annulus
with a sigma clipping criterion). From the luminosity growth curve it
establishes the Petrosian radius, inside which all the measurements are made.
The Sérsic fit is performed on the 1-D luminosity profile; for robustness, the
1-D outputs are used as inputs for a 2-D Sérsic fit of the galaxy pixels.
Finally, it measures several morphometric parameters, e.g. concentrations,
asymmetries, Gini and M20 (the former is a coefficient quantifying the
inequality among values of a frequency distribution, in this case of pixel
values; the latter is the second order moment, i.e. the flux values weighted
by the their square distance to the center, of the 20% brightest pixels; see
Lotz et al., 2004), entropy, spirality, curvature among others). In a
forthcoming version, the luminosity profile curvature (Lucatelli & Ferrari,
2019) will be used to provide a more robust input to a parametric model-based
fit of the light profile, eventually replacing the 1-D Sérsic fit as a metric,
mainly to mitigate a long lasting problem of Sérsic index determination (see
discussion in EMC2022b). Only the $I_{\scriptscriptstyle\rm E}$ SS fit was
provided.

###  3.4 ProFit

Profit (Robotham et al., 2017) is a software package designed to perform
Bayesian two-dimensional photometric galaxy profile modelling. It consists of
a low-level C++ library accessible via a command-line interface and documented
API, along with high-level R (v.3.6.1) and Python interfaces. The fitting
process for each object starts running the source finder ProFound (Robotham et
al., 2018)333https://github.com/asgr/ProFound on a $500\times 500$ pixel
cutout centred on each target, to create a segmentation map and find nearby
sources requiring simultaneous modelling; the output also provides some
reasonable initial guess for the profile solution. The actual fitting is then
performed by the Highlander core
software444https://github.com/asgr/Highlander, which combines a genetic
algorithm step with a CHARM (Turchin, 1971; Smith, 1984) Markov chain Monte-
Carlo (MCMC) process, repeated twice; each one is run for 100 steps (where
model realisations are modified by the number of free parameters also). The
CHARM algorithm is particularly useful on highly covariant parameter search,
but it is computationally expensive, because a single iteration requires
sampling all parameters. Since the kind of fitting used for ProFit is
relatively low in the number of parameters, but sometimes quite highly
covariant in the posterior, CHARM has proven to be a powerful exploration
tool. The provided solution is the combination of parameters that generate the
maximum likelihood given the per-pixel Data - Model residual. The parameter
priors are implicitly assumed to be uniform. Errors are estimated from the
final MCMC run, with full covariance matrix information available.

All requested data was provided. Partially building upon the effort put in the
EMC, the whole ProFit pipeline has recently been developed into a new package,
ProFuse (Robotham et al., 2022).

###  3.5 SourceXtractor++

SourceXtractor++ (Bertin et al., 2020; Kümmel et al.,
2020)555https://github.com/astrorama/SourceXtractorPlusPlus is a ground-up re-
write of the widely used SExtractor2 software (Bertin & Arnouts, 1996),
written in C++ with a strong focus on extensibility and model-fitting
photometry; the software is under active development and the results submitted
to the challenge represent a snapshot in this process (the version of
SourceXtractor++ used in the EMC is 0.12).

Each SourceXtractor++ run includes two stages, detection and measurement. The
detection stage follows the same procedure as used in SExtractor2. Detection
parameters need to be optimised for a compromise between the completeness of
the true object list and the number of spurious objects extracted or
deblended; over-extraction of sources impacts the performance of the run-time
required, and may also reduce the accuracy of morphological measurements if
objects are over-deblended. The parameters for the EMC were tuned aiming for
good overall performance, and therefore not for reaching $100\%$ completeness
of the input $5\sigma$ source list. The SS and DS simulations produce slightly
different distributions in apparent extent and surface brightness for the
galaxy images, and so the parameters governing detection were optimised
separately for the different simulations.

Measurement in one or several bands is controlled via a Python configuration
file with flexible model fitting at its core, which allows for the
simultaneous source analysis over a large number of FITS files with different
pixel grids. Various components (point source, exponential disc, free Sérsic ,
etc.) can be used individually or in combination; reasonable priors must
therefore be provided to the fitting engine in order to cover the range of
parameter values and provide sensible fits. The chosen priors for the EMC are
described in Appendix B.5.2.

Noticeably, the SourceXtractor++ pipeline used for the EMC includes a pre-
processing of the images, namely the extraction and usage of PSFs from images,
performed using the PSFEx software (Bertin, 2011), and a re-downsampling to
the original pixel scales of the NIR and LSST images ($None$ and $None$,
respectively). This procedure was allowed by the guidelines of the EMC, given
that no additional input data were used. However, given that no other
participants did anything similar, we decided to also check the performance of
the package on standard, non-pre-processed data, finding overall good
agreement with a few differences that should be taken into account when
considering the overall processing cost of the pipeline. We discuss this in
Appendix B.5 (see also EMC2022b). Additionally, SourceXtractor++ priors were
obtained by comparing the output distribution of a given morphological
parameter with the equivalent distribution for the provided samples of the
input true catalogues; the priors on parameters were iteratively adjusted
under the constraint that a simple analytical transfer function is required to
map each distribution to a Gaussian. Each parameter was calibrated
independently, without including covariances; only the statistical
distributions were used (i.e. there was no object-by-object comparison in the
process). A detailed description of the calibration of the priors is given in
Appendix B.5.2. All requested data were provided, except for the multi-band DS
fit with free $n_{\rm bulge}$.

##  4 Diagnostic metric

Given the high dimensionality of the output data, a straightforward comparison
of the results was not feasible. In order to obtain a reasonably comprehensive
overview of the quality of the performance, we defined an ad-hoc metric. The
participants provided catalogues that were matched to the input ones by means
of the unique ID of each source. Then for each run we proceeded to estimate
the difference between the input and the measured fluxes of each object, and
computed averaged statistical diagnostics.

Importantly, to compute such statistics we used a subset of sources from the
$5\sigma$-limit list, including only those for which all software packages
provided a meaningful fit. Recall that the default request was to provide a
fit for all the sources in the nominal $5\sigma$ list; with some minor caveats
and exceptions mentioned in Sect. 3, all participants obliged to this, and
some also provided results for lower S/N sources. However, not all the fits
were successful (i.e. some were given as NaN or default values in the output
catalogues), and some were flagged as ‘bad’ or ‘unreliable’ in one or more
codes. In general, these sources were not included in the lists used to
evaluate the accuracy of the results, which therefore only included the
objects for which all the software packages had provided a reliable fit. An
important exception is the case of Galapagos-2, which outputs several quality
flags describing which component of a galaxy can be considered as reliably
fit. In particular, for the DS runs a total of five flags were provided:
USE_FLAG_SS indicates that the single Sérsic fit is usable, as it did not run
into fitting constraints; USE_FLAG_BULGE_CONSTR and USE_FLAG_DISK_CONSTR serve
the same purpose for the bulge and disc components, respectively; in addition,
USE_FLAG_BULGE_BRIGHT and USE_FLAG_DISK_BRIGHT indicate whether the bulge and
the disc are relatively bright enough (${\rm b}/{\rm t}>0.2$ and ${\rm b}/{\rm
t}<0.8$, respectively), that their fit could in general be trusted, with the
additional difficulty that ${\rm b}/{\rm t}$ itself is defined via such a fit.
However, all of these flags were ignored in our analysis, to avoid excessive
complications in the definition of the common set of fitted sources within the
submissions. While this choice certainly impacts the statistics of the
results, because galaxies are taken into account that are known to violate
fitting constraints (they are 4% and 13% of the total number of objects, for
single and double component fits, respectively), we found that the effects on
the overall analysis was marginal. We stress that a general user of
Galapagos-2 shall consider these flags, according to their purposes; a
thorough description of the flags is provided in Häußler et al. (2022). See
the Appendix of EMC2022b for a comprehensive discussion on this topic.

To estimate the impact for each software tool, an additional term evaluating
the completeness fraction of the output catalogues with respect to the full
5$\sigma$ list was included in the global metric, as described below.

To build the metric, we started considering the relative flux difference of
each object with respect to the input true flux, i.e. $\delta f=(f_{\rm
meas}-f_{\rm true})/f_{\rm true}$. We then used $\delta f$ to evaluate three
diagnostics: the bias $\mathcal{B}$; the dispersion $\mathcal{D}$; and the
outlier fractions $\mathcal{O}$. In summary, the three diagnostics were first
averaged over the sources belonging to bins of input magnitude (we used 15
bins to divide the full interval of simulated magnitudes, from 14 to 28), to
quantify the impact of S/N; then these averages (normalised with weighting
factors) were summed, and further combined with the completeness $\mathcal{C}$
diagnostic, finally yielding a global score $\mathcal{S}$ for each field and
realisation. For $I_{\scriptscriptstyle\rm E}$, the values computed for each
simulated field were finally averaged to obtain a single figure; while this is
not strictly correct from a statistical point of view, given that the results
in the different fields were very similar we assume that the outcome is
sufficiently accurate. In more detail, the four quantities were defined as
follows.

  * •

The bias $\mathcal{B_{\rm bin}}$ is the median of $\delta f$ in each bin of
input magnitude, computed considering only the objects having $|\delta f|\leq
5\sigma_{\rm true}$, where $\sigma_{\rm true}$ is the standard deviation of
$\delta f$ for an ideal distribution of fluxes, that we obtained by perturbing
the input true values with a random realisation of observational Gaussian
noise consistent with the expected depth of each image (we imposed a minimum
value corresponding to $\delta f=0.02$ (2%), to avoid unrealistically small
values of $\sigma_{\rm true}$ at the bright end). An unbiased measurement
would yield $\mathcal{B_{\rm bin}}=0$. We then define the average value of the
bias as the weighted mean of its values across the magnitude bins,
$\mathcal{B}=\sum_{\rm bins}w_{\rm bin}\mathcal{|B_{\rm bin}|}$, where $w_{\rm
bin}$ is a weighting factor given by the fraction of objects in each bin of
true magnitude (to give more weight to highly populated bins) multiplied by
the logarithm of the median S/N in that bin (to give more weight to the fit of
bright objects). Note that while the values of $\mathcal{B_{\rm bin}}$ can be
positive or negative, $\mathcal{B}$ is defined to be positive.

  * •

The dispersion $\mathcal{D_{\rm bin}}$ is the ratio between $\sigma_{\rm
meas}$, i.e. the standard deviation of the distribution of $\delta f$ (again
only including objects within 5$\sigma_{\rm true}$) and $\sigma_{\rm true}$,
in each magnitude bin. The average dispersion is defined as
$\mathcal{D}=\sum_{\rm bins}w_{\rm bin}\mathcal{D_{\rm bin}}$.

  * •

The outlier fraction $\mathcal{O_{\rm bin}}$ is the number of objects having
$|\delta f|>5\sigma_{\rm true}$ divided by the total number of fitted objects
in each bin. These objects fall outside the expected distribution, and we
assume that their large bias is due to some systematic error in their
measurement (e.g. strong contamination from neighbours, or catastrophic
failure of the fit). Therefore they were not included in the statistics of
“well-behaved” sources, and were instead isolated into a separate diagnostic.
The average outlier fraction is defined as $\mathcal{O}=\sum_{\rm bins}w_{\rm
bin}\mathcal{O_{\rm bin}}$.

  * •

The completeness $\mathcal{C}$ is simply the number of objects for which a
successful fit was provided, divided by the total number of objects in the
input list of ${\rm S/N}>5$ sources (we do not weight this quantity by the
magnitude bins). A galaxy is considered not to be fit if there is no entry in
the provided output catalogue, or if the challenge participant flagged that
galaxy as a ‘bad fit’ (see discussions above). Each software package has
different ways of identifying unreliable fits, and we refer the reader to the
publications describing each code for additional information. Here, we simply
trusted the participants’ verdict on the reliability of their fits.

We point out that the definitions used in EMC2022b are very similar, but since
it is difficult to construct meaningful expectations for ideal perturbed
distributions of morphological parameters (corresponding to the $\sigma_{\rm
true}$ we use here for the fluxes), some differences were introduced. The
interested reader should therefore pay attention to these details.

Finally, the global diagnostic $\mathcal{S}$ for each run is defined as

$\mathcal{S}=(1-\mathcal{C})+k_{\mathcal{B}}\mathcal{B}+k_{\mathcal{D}}(\mathcal{D}-1)+k_{\mathcal{O}}\mathcal{O}$
(2)

where we subtract 1 from $\mathcal{D}$ when computing the final global
statistics, because when the dispersion is ‘ideal’ the ratio with $\sigma_{\rm
true}$ is 1, and we want the value of all diagnostics to be close to zero for
ideal fits. In this expression, the $k$ factors are multiplicative constants
assigned to each of the three diagnostics in an attempt to reasonably weight
their relative contributions.

We chose $k_{\mathcal{B}}=20.0$, $k_{\mathcal{D}}=0.6667$ and
$k_{\mathcal{O}}=5.0$. While these choices are to some extent arbitrary, it is
worth pointing out that the large differences in these values do not reflect
the actual weight given to each diagnostic; on the contrary, they were chosen
exactly to try and reach a reasonable balance between the three weighted
quantities. We argue that a fit might be defined as ‘optimal’ if it has e.g.
$\mathcal{C}=1.0$ (100% completeness above 5$\sigma$), $\mathcal{B}<0.015$
(1.5% median bias), $\mathcal{D}<1.33$ (dispersion no larger than 4/3 of that
from the perturbed true fluxes), and $\mathcal{O}<0.1$ (10% of outliers) in
all bins of magnitude; in the case of these exact values, applying the chosen
weights one gets $\mathcal{S}\simeq 0.3+0.22+0.5=1.02$ (the bin weights
$w_{\rm bin}$ are not relevant here). So, we see that if $\mathcal{S}\leq 1$
the fit can be considered as optimal; $\mathcal{S}\leq 1.33$ is very good;
$1.33<\mathcal{S}\leq 1.67$ is good; and $1.667<\mathcal{S}\leq 2.0$ is
acceptable. Finally, with this metric values of $\mathcal{S}$ much larger than
2 indicate a bad overall fit. Note that when marginalizing the contributions
for a 100% complete fit, $\mathcal{S}=2.0$ can be due to: a 10% overall
offset; a 3$\sigma_{\rm true}$ standard deviation; or a 40% outlier fraction.

The diagnostics were evaluated automatically by means of Python scripts, but
the results were also visualized graphically, to allow for sanity checks and
for a quick grasp of any particular features. Figure 8 shows an example of a
diagnostic plot that we used to analyse one of the provided output catalogues;
similar plots were created for the outputs for each field and realisation that
each participant provided. Each dot is a single fitted galaxy, and its $\delta
f$ is plotted against its true input magnitude in the considered band. The
dots are colour-coded by the true bulge-to-total ratio (which for the SS and
RM realisations is a proxy for the Sérsic index, $n=3\;{\rm b}/{\rm t}+1$).
For each bin of magnitude, the median, standard deviation, and outlier
fraction of the distribution were computed, and the values were then used to
compute the diagnostics described above; the dotted lines show the
$1\sigma_{\rm true}$ and $5\sigma_{\rm true}$ levels. Specific examples are
given in Appendix B.

Figure 8: Example of a diagnostic plot used to analyse each output catalogue
provided by the participants to the EMC. Each dot shows the $\delta f$ value
of a fitted galaxy, colour-coded by the ${\rm b}/{\rm t}$ value from the
original Egg catalogue (colour-bar on the right), as a function of the input
true magnitude in the considered band. The black solid line is the running
median of the distribution, which should be identically zero for a perfect
fit. The dotted lines show the (positive and negative) 1 and 5 $\sigma_{\rm
true}$ levels, used to compute the diagnostics described in the text and
reported in the top left corner of the plot (see Sect. 4). Figure 9: Visual
summary of the $\mathcal{S}$-values listed in Table 2.

Figure 10: Visual summary of the diagnostic quantities: (left to right)
absolute value of the median bias ($\mathcal{B}$), average ratio of
$\sigma_{\rm meas}$ and $\sigma_{\rm true}$ ($\mathcal{D}$), and outlier
fraction ($\mathcal{O}$). In the bias panel, the points corresponding to
ProFit DS multi-band runs are not shown, being off the scale (their values are
0.49 and 0.44, respectively). See text for details.

Because of a technical problem not related to the performance of the code and
gone unnoticed during the run, the processing of F1 by DeepLeGATo was
interrupted before the end of the input list, and the corresponding catalogue
was therefore incomplete. To ease the comparisons, and considering that the
problem was not caused by a bug in the code (since the processing of all the
other fields ended smoothly), we decided to remove F1 entirely from the
analysis process, after checking that this would not favour one of the codes
with respect to the others.

Finally, as already mentioned the Galapagos-2 runs on DS F0 were performed in
a simultaneous multi-band fit, including $I_{\scriptscriptstyle\rm E}$ with
the other bands; this caused the results to be significantly different from
those in the other four $I_{\scriptscriptstyle\rm E}$ fields. To avoid any
impact on the evaluation, and considering the many different approaches of the
other participants (DeepLeGATo only provided the $I_{\scriptscriptstyle\rm E}$
fits, Morfometryka did not provide the fit, ProFit only provided non-
simultaneous fits, and SourceXtractor++ provided both a simultaneous and a
separated fit), we decided to remove F0 from the analysis of
$I_{\scriptscriptstyle\rm E}$ DS.

In summary, $I_{\scriptscriptstyle\rm E}$ SS and RM were analysed by averaging
F0, F2, F3, and F4 (resulting in 212 000 objects for SS, and 204 229 for RM,
where the bright galaxies were not simulated); $I_{\scriptscriptstyle\rm E}$
DS by averaging F2, F3 and F4 (207 064 objects); and the other bands on DS F0
alone (because it was the only field simulated with a multi-band data set; it
contains 70 700 objects).

##  5 Results

Software package | SS | RM | DSb4 | DSbf | DSMb4 | DSMbf
---|---|---|---|---|---|---
DeepLeGATo | 1.44 (0.97) | – | 1.57 (0.96) | – | – | –
Galapagos-2 | 1.33 (0.90) | 2.37 (0.75) | 1.91/2.03 (0.95) | 1.91/2.51 (0.95) | 1.65/1.90 (0.54) | 1.90/2.00 (0.95)
Morfometryka | 2.35 (0.84) | – | – | – | – | –
ProFit | 1.09 (1.00) | 2.42 (0.94) | 3.04 (1.00) | 2.53 (1.00) | 4.02 (0.73) | 3.78 (0.73)
SourceXtractor++ | 0.81 (0.96) | 2.11 (0.87) | 1.19 (0.97) | 0.95 (0.97) | – | 0.99 (0.97)
Table 2: Software package, runs and $\mathcal{S}$-values considering the
common list of sources (see text for details). Lower $\mathcal{S}$ means
better performance. The acronyms in the columns refer to the different
realisations described in Sect. 2.1, with the following additional
specifications for the DS runs: b4 = fixed Sérsic index for bulge ($n=4$); bf
= free Sérsic index for bulge; M = multi-band (i.e. NIR and LSST bands,
excluding $I_{\scriptscriptstyle\rm E}$; no M is for $I_{\scriptscriptstyle\rm
E}$ band only). For the Galapagos-2 DS runs, we give the numbers of the best
and worst performance (first and second number, respectively), corresponding
to either the double Sérsic fit or a single Sérsic fit, which the software
produces in all cases (see Appendix B.2). The numbers in parenthesis are the
completeness $\mathcal{C}$, i.e. the fraction of sources from the input
5$\sigma$ list having a successful measurement in the output catalogue.

In this section we discuss the results obtained with the different software
packages. First of all, it is worth pointing out again that the complexity of
the challenge caused a significant scatter in the interpretation of its goals
and spirit by the participants. This caused substantial differences in the
adopted approaches, level of processing and output formats between them.
Together with the high dimensionality of the data set, this makes a direct and
comprehensive comparison of the results very challenging. In other words, the
strategies and techniques adopted by the participants influenced the overall
accuracy of the provided output, and this must be taken into account in the
analysis, to ensure a fair overview of each code’s capabilities and
limitations. Nevertheless, we believe it is possible to draw some interesting
general conclusions from the comparison. We will describe the overall outcomes
in the following, with some particular cases discussed in more detail, when
necessary.

Individual diagnostic plots for all the different runs are available in an on-
line interactive
tool.666https://share.streamlit.io/hbretonniere/euclid_morphology_challenge
Further discussions on the results provided by each participating team are
provided in Appendix B, together with a summary of the computational times and
memory workload required by each software package.

###  5.1 Global outcome

In the following we separately analyse the three realisations SS, DS, and RM.
We separate the multi-band data set from the $I_{\scriptscriptstyle\rm
E}$-only DS fits, since the results are significantly different (we identify
the multi-band case with the addition of the letter ‘M’ to the acronym DS
whenever necessary). The values of the global metric $\mathcal{S}$ for each
realisation are listed in Table 2 (where the values of the completeness factor
$\mathcal{C}$ are also reported, to allow for a better evaluation of its
impact), and shown in Fig. 9. A visual summary of the values of each
diagnostic used to compute $\mathcal{S}$ (i.e. $\mathcal{B}$, $\mathcal{D}$
and $\mathcal{O}$), averaged over the magnitude bins and fields, is given in
Fig. 10. For the multi-band case, these values are the average of the ones
obtained in each of the nine individual bands. Note that in these plots the
values are shown before the weighting by the $k$ factors in Eq. 2, so the
relative difference between the values obtained by any two codes for the three
diagnostics does not straightforwardly reflect their final difference in
$\mathcal{S}$.

Because the global diagnostic quantities only provide a crude overview of the
results, being averages over the input parameter space obtained with arbitrary
weights, in Sect. 5.2 we also present a collection of summary plots (Figs. 11
to 14), showing the trends of the diagnostics as a function of the input
magnitudes in all the cases of interest.

We want to begin emphasizing that each code proved to have points of strength
and of weakness, so the comparison of the global score $\mathcal{S}$ and of
its factors is only intended as a quick overview, and should by no means be
taken as a rigorous evaluation and ranking of the software packages. We fully
acknowledge that this is a simplified view and therefore alternative metrics,
tailored to specific science cases or assigning different weights to the
considered diagnostics, could result in different conclusions. That said, we
can claim that with our metric all software packages provided acceptable to
good results in at least some of the realisations. Some differences are
present in a few cases (e.g. a particular realisation or band, the faint end
of the simulated distribution of galaxies, etc.), but the outputs provided are
typically fairly accurate. Unsurprisingly, in the DS realisation the best
results were obtained in the $I_{\scriptscriptstyle\rm E}$ runs, for most of
the software packages, given its high resolution and depth. The multi-band
data set proved to be more demanding, because of the lower S/N and resolution
of the images, and also of the resampling procedure which introduced noise
correlations. All the participants reached at least 95% global completeness in
the $I_{\scriptscriptstyle\rm E}$ DS data set; in the other realisations, some
lower scores were obtained (see the numbers in parenthesis in Table 2).

In the $I_{\scriptscriptstyle\rm E}$-only runs, the best global results by all
codes were obtained on the SS realisation, with all software packages reaching
values of $\mathcal{S}$ between 1.0 and 1.7, with the exception of
Morphometryka ($\mathcal{S}=2.54$), penalised by a significant bias
$\mathcal{B}$ caused by systematic underestimation of fluxes of bright bulge-
dominated galaxies, and of all objects at faint magnitudes (see discussions in
Sect. 5.2 and in Appendix B.3). On the DS realisation of the
$I_{\scriptscriptstyle\rm E}$ band, SourceXtractor++ and DeepLeGATo reached
$\mathcal{S}<2.0$, although it must be recalled that the scheme adopted by the
latter (dividing the sources according to their nominal S/N bin) introduces
peculiar features in the distribution of $\delta f$ (see Appendix B.1);
Galapagos-2 reached $\mathcal{S}\simeq 2.0$, while ProFit was penalised by a
large $\mathcal{B}$ caused by a strong overestimation of fluxes at faint
magnitudes (see Appendix B.4). Interestingly, the fits with free $n_{\rm
bulge}$ were in general slightly better than those with fixed $n_{\rm
bulge}=4$ (except for Galapagos-2), despite the fact that the simulated
galaxies had indeed $n_{\rm bulge}=4$.

Finally, the results for the RM realisation were less accurate than the ones
for the other $I_{\scriptscriptstyle\rm E}$ images. This should be expected,
given the inherently more difficult task of fitting analytical profiles on
complex realistic morphologies. Interestingly, here the three codes that
provided results have very similar trends; this seems to imply that when
dealing with realistic galaxy shapes, the impact of prior calibrations, pre-
processing of the images, and robustness of the algorithm become of secondary
importance with respect to the inherent difficulty of the task.

The results for the multi-band data set show the evident impact of the
different strategies followed by the participants. SourceXtractor++ (which
included an image and PSF pre-processing pipeline, calibrated its priors on
the provided example ground-truth data, and fitted all bands simultaneously)
obtained optimal results, comparable to those for the
$I_{\scriptscriptstyle\rm E}$ band fitted alone. Galapagos-2 did perform a
simultaneous fit, but without the image and PSF pre-processing obtained sub-
optimal results. ProFit, which performed a separate fit on each band, had
weaker performance due to not using information from the
$I_{\scriptscriptstyle\rm E}$ image in the multiband fits: faint galaxies
detected in $I_{\scriptscriptstyle\rm E}$ are likely to have very low S/N or
could even be undetected in most NIR and LSST bands, and without the
$I_{\scriptscriptstyle\rm E}$ parameter to constrain the fit, it is very
difficult to properly model their light profiles in these bands. This outcome
is important for highlighting how the synergy with Euclid can significantly
improve the accuracy of Rubin/LSST measurements, as already pointed out by
several studies (see e.g. Rhodes et al., 2017; Capak et al., 2019).

It is worth mentioning here that any statistical result showing a dependence
on the bulge fraction of the sources is biased by the low fraction of
simulated bulge-dominated galaxies with respect to the real Universe
distribution (see Sect. 2.1), so the overall performance on real data might be
worse.

Figure 11: Trends of the completeness $\mathcal{C}$ with input
$I_{\scriptscriptstyle\rm E}$ magnitude of the output provided by the five
participants to the challenge, for (top to bottom) DS, SS, and RM realisations
(in the latter, points start at $I_{\scriptscriptstyle\rm E}=20.5$ because
bright objects were not included in the simulated images).

###  5.2 Trends of the diagnostics with input magnitudes

In Fig. 11 the trends of completeness $\mathcal{C}$ are shown for each
software package as a function of $I_{\scriptscriptstyle\rm E}$ in bins of
magnitude, for the three realisations DS, SS and RM (we remind the reader that
for the latter only objects with $I_{\scriptscriptstyle\rm E}>20.5$ were
simulated).

As one can expect, the fraction of successfully fitted objects decreases with
increasing magnitude for all codes, since fainter objects are generally harder
to detect and fit; the only exception is DeepLeGATo, for which bright objects
are more often prone to failure. This is a consequence of the fixed stamp
sizes that are used in the current version of the software, so that very
bright (and large) objects may be larger than the stamp size.

Overall, down to $I_{\scriptscriptstyle\rm E}\simeq 23.5$ all codes
successfully fit more than $95\%$ of the galaxies; the completeness then
typically decreases, with the noticeable exception of ProFit which always
stays close to $100\%$ for the SS and DS realisations. This is likely due to
the different thresholds used to perform detection in the software packages
(we recall that Galapagos-2 uses SExtractor and ProFit uses ProFound), leading
to different efficiencies in the actual detection of sources.

In each panel of Figs. 12 to 14, the trends of one of the diagnostics
$\mathcal{B_{\rm bin}}$, $\mathcal{D_{\rm bin}}$, and $\mathcal{O}_{\rm bin}$
are shown as a function of magnitude in the considered band, for all the
software packages that have provided a fit in the relevant realisations. Some
individual cases of particularly notable behaviour are described in more
detail in Appendix B.

####  5.2.1 Bias in $I_{\scriptscriptstyle\rm E}$ runs

In the $I_{\scriptscriptstyle\rm E}$ DSb4 runs, for Galapagos-2 and ProFit,
the typical absolute bias in the measured flux $\mathcal{B_{\rm bin}}$ is
below $1\%$ for bright sources ($I_{\scriptscriptstyle\rm E}<23$), increasing
to $10$–$15$% at the faint end. This is likely due to contamination from
nearby brighter sources, and/or to the inherent difficulty in fitting low S/N
objects with analytical profiles. DeepLeGATo has a less stable trend, with
slightly larger values of bias for intermediate and bright sources, but a
lower bias at the faint end. Finally, SourceXtractor++ shows the most stable
trend, without a strong overestimation at the faint end, despite a more
pronounced average bias of $2$–$3\%$.

In the SS runs, in most cases we see more stable trends, with monotonic trends
towards overestimation of fluxes with decreasing brightness reaching about 5%
bias at the faint end. Noticeably, ProFit behaves differently, slightly
underestimating intermediate magnitude sources before starting to overestimate
at the faint end. An even more striking exception is given by Morphometryka
(SS was the only provided output data set), which has a clearly declining
monotonic trend, reaching $\delta f\simeq-15\%$ at $I_{\scriptscriptstyle\rm
E}\simeq 25.3$.

The situation is completely different in the RM realisation, where all three
of the codes that provided results have similar declining trends in the bias
$\mathcal{B}$, with faint sources typically having fluxes underestimated by
about $10$–$15$% at $I_{\scriptscriptstyle\rm E}=24$–$25$. This is at odds
with the results from the other realisations. We checked that this is not an
issue in the simulations: while minor inconsistencies between the input true
fluxes in the catalogue and the actual realisations of the sources in the
images can be present because of the simulation method (see Sect. 2.1), we
found that the impact of this is negligible, with a typical mean offset of
about 0.05% and some scatter in the values that is not sufficient to explain
the global trends of the measurements. We postpone further investigation on
this topic to future work, given that it does not strongly impact the present
analysis; the trend is almost identical for the three considered software
packages, leaving the comparison among them essentially unchanged.

####  5.2.2 Dispersion in $I_{\scriptscriptstyle\rm E}$ runs

The dispersion $\mathcal{D_{\rm bin}}$ for all codes is typically comparable
to $\sigma_{\rm true}$ at the faint end in all realisations. In DSb4, at
$I_{\scriptscriptstyle\rm E}<23$ there is a hierarchy of performance with
DeepLeGATo reaching values around 3.0 (again probably because of the limited
dimensions of the stamps), Profit reaching around 2.0, SourceXtractor++ 1.5
and Galapagos-2 going from 0.5 to 1.5. In SS this hierarchy is less pronounced
with all codes, including Morfometryka, staying below $2.0$ at all magnitudes;
the evident exception is again DeepLeGATo. In the RM case again we see similar
(and sub-optimal) trends for the three codes that provided results.

####  5.2.3 Outliers in $I_{\scriptscriptstyle\rm E}$ runs

The outlier fraction $\mathcal{O_{\rm bin}}$ in DSb4 stays below 10% at all
magnitudes for SourceXtractor, and goes from very low values to 20% at faint
magnitudes for Galapagos-2. Profit reaches 30%, while again DeepLeGATo suffers
from the limited dimensions of the stamps, reaching 100% outliers at the
bright end, while remaining close to zero at $I_{\scriptscriptstyle\rm E}>23$.
In SS, SourceXtractor and Galapagos-2 stay below $10\%$ at all magnitudes,
while Morfometryka has large values (around 50%) at the bright end, likely
because of a sub-optimal estimate of bulge-dominated sources (see Appendix
B.3). Again in the RM case there are no major differences between the quality
of the performance, except for SourceXtractor performing better on
$I_{\scriptscriptstyle\rm E}<22$ objects.

####  5.2.4 The free Sérsic bulge case

There is no substantial difference between the fixed and free $n_{\rm bulge}$
cases in the DS realisation. We do not show the trends for the free $n_{\rm
bulge}$ case; as mentioned, the latter yields slightly better results than the
fixed $n_{\rm bulge}=4$ case for SourceXtractor++, reducing the bias at all
magnitudes, and for ProFit, which has better trends for all diagnostics; for
Galapagos-2, the bias is almost identical in the two runs, while the
dispersion and the outlier fraction have opposite trends (the free case having
many more bright outliers but, simultaneously, a lower dispersion ratio for
the few “well-behaved” sources), resulting in a similar final score (see Fig.
10).

####  5.2.5 The multi-band data set

In the multi-band data set (for which only Galapagos-2, ProFit, and
SourceXtractor++ provided results; see Figs. 13 and 14) the trends are
generally similar to those of the $I_{\scriptscriptstyle\rm E}$ case, with
larger values of bias at the faint end for all codes in the NIR and LSST
bands. In particular, Profit reaches $\mathcal{B_{\rm bin}}\simeq 1.5$ in the
NIR bands, and approximately 4 in the LSST $u$ and $g$ bands, likely because
of the independent fits performed on low S/N sources (we do not show the
corresponding points in the plots, for readability); there is also a
particular trend in the LSST bands, with sources being underestimated at
intermediate magnitudes before turning to a strong overestimation at the faint
end (we were not able to find an easy explanation for this). Galapagos-2 and
SourceXtractor obtain better (and similar) results in NIR thanks to the
simultanous fit, reaching $\mathcal{B_{\rm bin}}\simeq 0.1$; in LSST they
behave similarly as well, with SourceXtractor++ performing clearly better only
in the $u$ band. Here the free $n_{\rm bulge}$ fit generally yields slightly
worse results in dispersion and outlier fraction for Galapagos-2, and slightly
better results for ProFit, than the fixed $n_{\rm bulge}$ fit (not shown in
these plots; see again Fig. 10).

Figure 12: DSb4 (top), SS (center), and RM (bottom) summary plots for
$I_{\scriptscriptstyle\rm E}$. Left to right: median bias $\mathcal{B}$;
dispersion (standard deviation) $\mathcal{D}$; and fraction of outliers
$\mathcal{O}$. Each line and colour correspond to a different code as
indicated in the legend; note the different $y$-axis scales. Figure 13: Multi-
band summary plots for $I_{\scriptscriptstyle\rm E}$ and the NIR bands. Left
to right: median bias $\mathcal{B}$; dispersion (standard deviation)
$\mathcal{D}$; and fraction of outliers $\mathcal{O}$. Top to bottom:
$I_{\scriptscriptstyle\rm E}$, $Y_{\scriptscriptstyle\rm E}$,
$J_{\scriptscriptstyle\rm E}$, and $H_{\scriptscriptstyle\rm E}$. Each line
and colour correspond to a different code as indicated in the legends; note
the different $y$-axis scales. Figure 14: Same as Fig. 13, for LSST bands. Top
to bottom: $u$, $g$, $r$, $i$, $z$. Note the different $y$– axis scales.

###  5.3 Separated bulge and disc estimates

Figure 15: Trends for bulges only, in the multi-band DSb4 fit, for
$I_{\scriptscriptstyle\rm E}$ and the NIR bands. Left to right: median bias
$\mathcal{B}$; dispersion (standard deviation) $\mathcal{D}$; and fraction of
outliers $\mathcal{O}$. Top to bottom: $I_{\scriptscriptstyle\rm E}$,
$Y_{\scriptscriptstyle\rm E}$, $J_{\scriptscriptstyle\rm E}$, and
$H_{\scriptscriptstyle\rm E}$. Each line and colour correspond to a different
code as indicated in the legends; note the different $y$-axis scales. Figure
16: Same as Fig. 15, but fot discs only. Figure 17: Trends for bulges only, in
the multi-band DSb4 fit (as in Fig. 15), for LSST bands. Top to bottom: $u$,
$g$, $r$, $i$, and $z$. Figure 18: Same as Fig. 17, but for discs only.

So far we have considered total fluxes, but it is also instructive to
investigate how the software packages performed in the separate flux
measurements of the two components of each galaxy (bulge and disc) in the DS
realisation. It is worth stressing that both estimates can be individually
worse than the total flux one, if their sum is close to the true total value,
but the partition among bulge and disc is not well recovered; this is linked
to the accuracy of the morphological parameter estimates, discussed in
EMC2022b.

In Figs. 15 to 18 we show the summary plots for bulges and discs separately.
In $I_{\scriptscriptstyle\rm E}$, evident features are the bias trends for
Galapagos-2 and ProFit, both overestimating the flux of bulges at the faint
end and underestimating that of discs at all magnitudes, by approximately 5%
down to $I_{\scriptscriptstyle\rm E}\simeq 22$ and then getting worse;
DeepLeGATo and SourceXtractor++ show a strong underestimating trend for faint
bulges, while showing a reasonable accuracy for discs at all magnitudes
(DeepLeGATo overestimating their flux by a few percent). The dispersion and
the outlier fraction are generally larger for bulges (often reaching very high
values) than for discs, in particular at intermediate and faint magnitudes,
where discs show nearly optimal values. For bulges, the dispersion has very
similar trends for all codes, while for discs SourceXtractor++ performs
better, followed by Galapagos-2, while ProFit and DeepLeGATo suffer at the
bright end. SourceXtractor++ also yields the best performance concerning the
outlier fraction, both for bulges and for discs.

Similar observations can be made for the other bands. In the NIR bands, the
bias of bulges for Galapagos-2 and SourceXtractor++ show a similar
underestimation of 3–5% at intermediate magnitudes, and all codes show a
strong overestimation at the faint end, particularly dramatic for ProFit. For
discs, SourceXtractor++ stays close to zero quite firmly, Galapagos-2 shows an
underestimation at the faint end balancing out the bulge overestimation, while
ProFit again shows an overestimation. In LSST bands, Galapagos-2 shows the
same opposite trend for bulges and discs at faint magnitudes; interestingly,
SourceXtractor++ tends towards the reverse (underestimating bulges and
counterbalancing with an overestimation of discs). ProFit still shows the
strong rising trend in all cases.

The trend of the dispersion for bulges is very similar for the three codes,
all having $\mathcal{D}\simeq 2.0$–$3.0$ for the bulk of the bins in all
bands; for discs, better values are found towards the faint end, again with
SourceXtractor++ performing slightly better than the other two codes. Finally,
bulges typically have very high outlier fractions, typically above $50\%$ in
all bins in the NIR bands, while having lower values at the bright end in the
LSST bands, but always peaking at values above $75\%$ at intermediate
magnitudes; on the other hand, discs have a lower fraction of outliers at the
faint end.

###  5.4 Colours in the multi-band data set

Figure 19: Colours estimates for Galapagos-2. Each small sub-panel on the left
part of the plot shows the measured versus true colour (magnitude difference)
for a pair of bands: the colour is always $x$-axis label $-$ $y$-axis label.
The numbers in each panel are the NMAD and $\eta$, as defined in the text. The
upper larger panel shows the values of NMAD and $\eta$ for the colours of all
bands with respect to $I_{\scriptscriptstyle\rm E}$, as reported in the
legend. The rightmost larger panel shows the $I_{\scriptscriptstyle\rm E}-g$
versus $I_{\scriptscriptstyle\rm E}-H_{\scriptscriptstyle\rm E}$ colour-colour
diagram; points are the measured colours for individual galaxies, while the
black lines are density contour levels of the corresponding true colours
distribution. Figure 20: Same as Fig. 19, for ProFit. Figure 21: Same as Fig.
19, for SourceXtractor++.

Thus far, we have focused on the accuracy of the flux estimates. However, for
the multi-band data set it is also interesting to examine how accurately the
three software packages that provided results were able to recover flux ratios
among the bands, i.e. colours. Here too, the multi-dimensionality of the data
requires some effort to summarize the results in a few informative plots. In
each small sub-panel of Figs. 19 (Galapagos-2), 20 (ProFit), and 21
(SourceXtractor++), we show the comparison between a measured colour and its
corresponding true value. For each fitted galaxy (colour-coded by its true
$I_{\scriptscriptstyle\rm E}$ magnitude) the measured colour is the $y$-axis
coordinate, and the true colour is the $x$-axis coordinate (so a perfect
colour estimate would result in a diagonal line from the bottom left to the
top right corners; for reference, we plot this line in each panel). The
figures show all possible combinations between the nine bands, giving an
overview of the results. In each small sub-panel, the considered colour is
given by the $x$-axis band label and the $y$-axis band labels, so for example
the first sub-panel in the upper left corner is the $I_{\scriptscriptstyle\rm
E}-Y_{\scriptscriptstyle\rm E}$ colour, the one below it is the
$I_{\scriptscriptstyle\rm E}-J_{\scriptscriptstyle\rm E}$ colour, and so on. A
distribution of points with a narrow vertical strip, which is often seen (e.g.
in colours including $I_{\scriptscriptstyle\rm E}$), imply that while the true
colours are close to zero, the measured ones have a large dispersion; if the
distribution of points follows the diagonal the colour is estimated with good
accuracy.

In each panel the values of two statistical diagnostics are also shown.
Defining $\delta c=|{\rm colour}_{\rm meas}-{\rm colour}_{\rm true}|/(1+{\rm
colour}_{\rm true})$, the first number is the normalised median absolute
deviation, NMAD $=1.48\mathcal{M}(\delta c)$ considering the sources with
$\delta c\leq 0.2$ ($\mathcal{M}$ is the median of the distribution); and the
second one is $\eta$, the percentage of outliers having $\delta c>0.2$. The
upper bigger panel in the figures shows these two values for the colours
including $I_{\scriptscriptstyle\rm E}$, as reported in the legend. As a
further example, the lower bigger panels show a colour-colour diagram,
$I_{\scriptscriptstyle\rm E}-g$ versus $I_{\scriptscriptstyle\rm
E}-H_{\scriptscriptstyle\rm E}$; again, plotted are all measured colours of
individual galaxies, colour-coded by their input $I_{\scriptscriptstyle\rm E}$
magnitudes; the density contours shown as solid lines represent the true
colour distribution.

Obviously, fainter sources have more scattered distributions. In all cases,
the colours including the $u$ band are clearly the least accurate, with large
NMAD and $\eta$; this is easily understandable given the fact that deep-sky
galaxies are typically fainter in blue bands, because of intrinsic properties,
dust absorption, and cosmological redshifting; also, the band itself is quite
shallow ($m_{\rm lim}=23.6$ at 10$\sigma$).

These diagrams again show how the impact of the different approaches used by
the three participants affect the results. Profit, using independent fits in
each band, obtains the least accurate results both in terms of NMAD and
$\eta$. Galapagos-2 and SourceXtractor++ perform significantly better,
probably because of the multi-band simultaneous fit. Between the two,
SourceXtractor++ yields slightly better results, most likely due to the
careful prior calibration used (see Appendix B.5.2).

###  5.5 Uncertainty budgets

Figure 22: Summary plots showing the fraction of sources satisfying the
relation $|f_{\rm meas}-f_{\rm true}|<\sigma_{\rm meas}$ in each magnitude
bin, for the $I_{\scriptscriptstyle\rm E}$ band runs. Expected fractions
should be equal to 0.683 (horizontal dashed line). Left to right: DSb4, SS,
RM. Each symbol and colour correspond to a different code, as indicated in the
legends. Figure 23: Same as Fig. 22, for the DSb4 multi-band run. Each panel
refers to a band, as indicated in the titles (note that
$I_{\scriptscriptstyle\rm E}$ is the last panel, for readability of the plot).

We computed a final quantity of interest, the fraction of sources that have a
measurement of flux consistent with the true value within their nominal
uncertainty budget, as provided in the output catalogues: i.e., $|f_{\rm
meas}-f_{\rm true}|<\sigma_{\rm meas}$ (e.g. Häußler et al., 2007). This
fraction should in principle be 68.3% (assuming a Gaussian distribution of the
measurements, a reasonable assumption given that the noise map is
predominantly Gaussian at least when the photon noise is not relevant, i.e.
for intermediate and faint galaxies).

For the codes that provided magnitudes, we computed the errors on fluxes as
$\sigma_{f}=0.921f\,\sigma_{m}$. For the double component runs, we estimated a
global uncertainty with the usual approach, i.e. $\sigma_{\rm tot}={\rm
b}/{\rm t}\;\sigma_{\rm bulge}+(1-{\rm b}/{\rm t})\,\sigma_{\rm disc}$, since
we cannot consider the errors on the two components as independent. We show
the results in Figs. 22 ($I_{\scriptscriptstyle\rm E}$) and 23 (multi-band
DSb4). It is evident that uncertainties are typically underestimated for
bright sources despite their per pixel RMS being enhanced by the contribution
of photon noise, while for faint sources they are mostly correct, in a
statistical sense.

A possible explanation for these trends is that bright sources are fit with
less uncertainty because of their high S/N, but the uncertainty budget (which
is based on the information in the RMS map) does not account for systematic
errors that cause the actual offsets in the measurements. Note that in the
evaluation of the bias $\mathcal{B}$ we considered the relative offset $\delta
f$, which is typically small for bright sources even when the absolute offset
$|f_{\rm meas}-f_{\rm true}|$ is large.

At the faint end, ProFit yields substantially less accurate results than
Galapagos-2 and SourceXtractor++, in all bands. It must be pointed out that
the version of DeepLeGATo that was used for the EMC did not include a
consistent algorithm for the uncertainty budget estimation. Therefore, we have
not included DeepLeGATo in this analysis. A new version including error
estimation is under development. The uncertainty estimation in Morfometryka
takes into account the RMS map, however it only represents the fit parameters
uncertainty, since this is extracted from the covariance matrix of the fit.
For fluxes and magnitudes, the quoted uncertainty is only due to the fitted
radius and underlying region used to integrate the flux.

##  6 Conclusions

This is the first of two papers presenting and discussing the Euclid
Morphology Challenge, a project aimed at comparing the results of model-
fitting software packages on a set of simulated Euclid observations, with the
aim of providing a first quantitative study to define the best-suited
algorithm to be included in the Euclid Science Ground Segment processing
pipeline. Five model-fitting tools (DeepLeGATo, Galapagos-2, Morfometryka,
ProFit, and SourceXtractor++) were tested in the challenge, providing partial
or complete results. It is worth stressing that the results obtained for the
EMC and discussed in this paper and in EMC2022b were used by the developers to
improve and in some cases significantly upgrade the software packages, which
is per se a relevant outcome of this work.

While the companion paper (Euclid Collaboration: Bretonnière et al. 2022)
focuses on results concerning the morphological parameters, in this work we
have presented the simulated data set, and focused on the results concerning
photometric estimates only.

We used the code Egg (Schreiber et al., 2017) to produce mock galaxy
catalogues, which were then exploited to create synthetic images in three
different realisations: two using GalSim (Rowe et al., 2015), to generate
single and double component Sérsic analytical profiles (SS and DS
respectively, the latter with $n_{\rm bulge}$=4 and $n_{\rm disc}$=1 and
varying bulge-to-total luminosity ratios); and one with the method described
in Euclid Collaboration: Bretonnière et al. (2022) and Lanusse et al. (2021)
to obtain realistic morphologies by means of a CNN technique (RM). We did not
simulate irregular galaxies, which constitute a substantial fraction of the
high-redshift population, and we did not include AGNs. Furthermore, the
distribution of the bulge-to-total ratios of the simulated galaxies is more
skewed towards disc-dominated objects than in the real Universe (see Fig. 4);
this implies that any result showing a dependence on this feature must be
considered as statistically optimistic.

We created five fields of view in the Euclid $I_{\scriptscriptstyle\rm E}$
band with an area of 0.482 deg2 each, and for one of them we also simulated
eight additional bands in the double Sérsic realisation (Euclid NIR
$Y_{\scriptscriptstyle\rm E}$, $J_{\scriptscriptstyle\rm E}$, and
$H_{\scriptscriptstyle\rm E}$, plus five Rubin/LSST bands $u$, $g$, $r$, $i$
and $z$). The images were produced with zero background flux (i.e., already
‘backgroud subtracted’). Gaussian noise realisations mimicking the expected
depths of the Euclid Wide Survey were added to the noiseless images. The
analytical profiles of the galaxies were finally randomized, following a
consistent Poissonian distribution. RMS maps including photon noise were also
simulated, and provided to the participants along with the PSFs of each field,
and lists of objects having nominal ${\rm S/N}>5$ in $I_{\scriptscriptstyle\rm
E}$, for which the actual centroid positions in pixel coordinates were given.
The requested output consisted of the photometric (and morphological)
estimates for all the objects in these lists.

To allow for a quantitative analysis of the results and for a comparison of
the accuracy of the codes, we defined a metric as described in Sect. 4,
computing three quantities related to, respectively, the median bias in the
flux measurements ($\mathcal{B}$), the mean dispersion of the measured flux
distribution ($\mathcal{D}$), and the outlier fraction ($\mathcal{O}$); these
quantities were evaluated on a subset of the data set including only the
sources for which all participants provided a fit, and we also computed a
further value quantifying the fraction of sources of the input list for which
a software package succeeded to provide a successful fit ($\mathcal{C}$).
Finally, we combined them in a single figure of merit $\mathcal{S}$; the lower
the value of $\mathcal{S}$, the better the global fit, with $\mathcal{S}\leq
1$ being considered optimal performance. The global results are summarised in
Table 2 and Figs. 9 to 14. We also analysed how accurately colours were
retrieved, and the reliability of the uncertainty estimates. The computational
times and some in-depth analysis for each individual software package are
discussed in the Appendix.

In general, all the participants provided acceptable to good results in at
least some of the tests, with a few differences from case to case. A thorough
comparison could only be performed with the results of ProFit, Galapagos-2,
and SourceXtractor++ (the first two provided all the requested outputs, and
the third only lacked one part of a set, namely the free $n_{\rm bulge}$ run
in the double Sérsic multi-band realisation). DeepLeGATo provided fits only
for the $I_{\scriptscriptstyle\rm E}$ analytical realisations, yielding good
results with the caveat of a S/N-dependent calibration of the models, while
Morfometryka provided data only for the single Sérsic set; therefore, their
contribution to the challenge can only be considered as indicative of their
potential.

Considering the $\mathcal{S}$ values, among the three codes that provided the
complete output the best performance was obtained by SourceXtractor++, in all
realisations. It is worth stressing again that the pipeline it adopted for the
EMC was the most tuned on the data set, with the inclusion of a substantial
pre-processing of the data, and the priors being modeled using the samples of
true values that were provided along with the simulated images. In the
analytical realisations (SS and DS), SourceXtractor++ and Galapagos-2 reached
typical values in the bias of the measured fluxes with respect to the true
values below 1% at S/N above 10 in all bands, and above 5 in
$I_{\scriptscriptstyle\rm E}$ (in at least one of the possible configurations
of the codes for which results were provided). For these two codes, the
dispersion is typically slightly larger than, but comparable to, the one
expected from a theoretical distribution obtained by perturbing the input true
fluxes with simulated observational noise consistent with the nominal depths
of the images; a small fraction of outliers (objects having a relative error
in the fitted flux larger than 5 times the expected theoretical distribution)
is always present. ProFit also yielded good results in SS, but was less
accurate on the DS realisation (Galapagos-2 also being sub-optimal at the
faint end in $I_{\scriptscriptstyle\rm E}$ for this realisation). The three
software packages performed similarly on the RM data set, and in all cases
results were poorer with respect to the analytical realisations; this seems to
imply that the overall quality of the fits might be sub-optimal when dealing
with real data. However, further investigation is needed to assess to what
extent the quality of the simulation might have impacted the results.

All codes tend to underestimate the uncertainties at the bright end, but
statistically provide a reasonable estimate for intermediate and faint
sources. Finally, a relevant outcome is the importance of the simultaneous
fitting of the multi-band data set: using the information provided by the deep
and high resolution $I_{\scriptscriptstyle\rm E}$ band helps to obtain more
accurate measurements on the shallower NIR bands, and noticeably also on the
ground-based LSST bands, both in terms of absolute fluxes and colours. Once
again this highlights how the synergy between the two surveys will be of
paramount important for the success of both.

We did not attempt to investigate the performance of the model-fitting
algorithms in the context of the Euclid Deep Survey (see Laureijs et al.,
2011). It is reasonable to expect that the results should be similar to those
obtained for the Wide Survey at equivalent S/N, although the increased level
of contamination and blending might worsen the performance. Dedicated
simulations would be required to quantitatively confirm this hypothesis, but
this is beyond the scope of the present study.

The results of the EMC set the baseline to decide which algorithm is to be
implemented in the Euclid pipeline. One possible way to proceed might be to
use the computational power of DeepLeGATo to provide fast and reliable priors
to SourceXtractor++; we will investigate this option in the next future.
However, given the many different approaches and techniques adopted by the
participants (different use of priors, different pre-processing strategies,
different approaches for multi-band fitting, etc.) and the other parameters
not included in the metric that should be taken into account (such as the
computing time, required resources, compatibility with the current pipeline,
accuracy of uncertainty budget estimates etc.), it is important to stress that
the results presented in this work and in EMC2022b must be interpreted with
caution. We followed one possible approach to analyse in a compact way a very
complex and multi-layered data set; and we emphasise that some of the software
packages might be better suited than the ones obtaining the best score here,
for other specific science cases. We invite the interested reader to check the
full set of results using the online tool.

Future work will include testing at least some of the software packages in
more realistic environments using Euclid official simulations, and finally the
implementation in the Euclid pipeline of the chosen algorithms.

###### Acknowledgements.

The Euclid Consortium acknowledges the European Space Agency and a number of
agencies and institutes that have supported the development of Euclid, in
particular the Academy of Finland, the Agenzia Spaziale Italiana, the Belgian
Science Policy, the Canadian Euclid Consortium, the French Centre National
d’Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt, the Danish
Space Research Institute, the Fundação para a Ciência e a Tecnologia, the
Ministerio de Ciencia e Innovación, the National Aeronautics and Space
Administration, the National Astronomical Observatory of Japan, the
Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the
Romanian Space Agency, the State Secretariat for Education, Research and
Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom
Space Agency. A complete and detailed list is available on the Euclid web site
(http://www.euclid-ec.org). H. Hildebrandt is supported by a Heisenberg grant
of the Deutsche Forschungsgemeinschaft (Hi 1495/5-1) as well as an ERC
Consolidator Grant (No. 770935).

## References

  * Barden et al. (2012) Barden, M., Häußler, B., Peng, C. Y., McIntosh, D. H., & Guo, Y. 2012, MNRAS, 422, 449
  * Bertin (2009) Bertin, E. 2009, Mem. Soc. Astron. Italiana, 80, 422
  * Bertin (2011) Bertin, E. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 442, Astronomical Data Analysis Software and Systems XX, ed. I. N. Evans, A. Accomazzi, D. J. Mink, & A. H. Rots, 435
  * Bertin & Arnouts (1996) Bertin, E. & Arnouts, S. 1996, A&AS, 117, 393
  * Bertin et al. (2002) Bertin, E., Mellier, Y., Radovich, M., et al. 2002, in Astronomical Society of the Pacific Conference Series, Vol. 281, Astronomical Data Analysis Software and Systems XI, ed. D. A. Bohlender, D. Durand, & T. H. Handley, 228
  * Bertin et al. (2020) Bertin, E., Schefer, M., Apostolakos, N., et al. 2020, in ADASS XXIX, ASP Conf. Ser. (San Francisco: ASP)
  * Capak et al. (2019) Capak, P., Cuillandre, J.-C., Bernardeau, F., et al. 2019, arXiv:1904.10439
  * Conselice (2003) Conselice, C. J. 2003, ApJS, 147, 1
  * Cropper et al. (2016) Cropper, M., Pottinger, S., Niemi, S., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9904, Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave, ed. H. A. MacEwen, G. G. Fazio, M. Lystrup, N. Batalha, N. Siegler, & E. C. Tong, 99040Q
  * de Albernaz Ferreira & Ferrari (2018) de Albernaz Ferreira, L. & Ferrari, F. 2018, MNRAS, 473, 2701
  * Euclid Collaboration: Bretonnière et al. (2022) Euclid Collaboration: Bretonnière, H., Huertas-Company, M., Boucaud, A., et al. 2022, A&A, 657, A90
  * Euclid Collaboration: Martinet et al. (2019) Euclid Collaboration: Martinet, N., Schrabback, T., Hoekstra, H., et al. 2019, A&A, 627, A59
  * Euclid Collaboration: Scaramella et al. (2022) Euclid Collaboration: Scaramella, R., Amiaux, J., Mellier, Y., et al. 2022, A&A, 662, A112
  * Euclid Collaboration: Schirmer et al. (2022) Euclid Collaboration: Schirmer, M., Jahnke, K., Seidel, G., et al. 2022, arXiv e-prints, arXiv:2203.01650
  * Ferrari et al. (2015) Ferrari, F., de Carvalho, R. R., & Trevisan, M. 2015, ApJ, 814, 55
  * Graham & Driver (2005) Graham, A. W. & Driver, S. P. 2005, PASA, 22, 118
  * Grogin et al. (2011) Grogin, N. A., Kocevski, D. D., Faber, S. M., et al. 2011, 197, 35
  * Häußler et al. (2013) Häußler, B., Bamford, S. P., Vika, M., et al. 2013, MNRAS, 430, 330
  * Häußler et al. (2007) Häußler, B., McIntosh, D. H., Barden, M., et al. 2007, ApJS, 172, 615
  * Häußler et al. (2022) Häußler, B., Vika, M., Bamford, S. P., et al. 2022, arXiv:2204.05907
  * Huertas-Company et al. (2015) Huertas-Company, M., Pérez-González, P. G., Mei, S., et al. 2015, ApJ, 809, 95
  * Koekemoer et al. (2011) Koekemoer, A. M., Faber, S. M., Ferguson, H. C., et al. 2011, ApJS, 197, 36
  * Kümmel et al. (2020) Kümmel, M., Schefer, M., Apostolakos, N., et al. 2020, in ADASS XXIX, ASP Conf. Ser. (San Francisco: ASP)
  * Lanusse et al. (2021) Lanusse, F., Mandelbaum, R., Ravanbakhsh, S., et al. 2021, MNRAS, 504, 5543
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv:1110.3193
  * Lotz et al. (2004) Lotz, J. M., Primack, J., & Madau, P. 2004, AJ, 128, 163
  * Lucatelli & Ferrari (2019) Lucatelli, G. & Ferrari, F. 2019, MNRAS, 489, 1161
  * Maciaszek et al. (2016) Maciaszek, T., Ealet, A., Jahnke, K., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9904, Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave, ed. H. A. MacEwen, G. G. Fazio, M. Lystrup, N. Batalha, N. Siegler, & E. C. Tong, 99040T
  * Mandelbaum et al. (2012) Mandelbaum, R., Hirata, C. M., Leauthaud, A., Massey, R. J., & Rhodes, J. 2012, MNRAS, 420, 1518
  * Merlin et al. (2016) Merlin, E., Bourne, N., Castellano, M., et al. 2016, A&A, 595, A97
  * Merlin et al. (2015) Merlin, E., Fontana, A., Ferguson, H. C., et al. 2015, A&A, 582, A15
  * Merlin et al. (2019) Merlin, E., Fortuni, F., Torelli, M., et al. 2019, MNRAS, 490, 3309
  * Peng et al. (2002) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2002, AJ, 124, 266
  * Peng et al. (2010) Peng, C. Y., Ho, L. C., Impey, C. D., & Rix, H.-W. 2010, AJ, 139, 2097
  * Peterson et al. (2015) Peterson, J. R., Jernigan, J. G., Kahn, S. M., et al. 2015, ApJS, 218, 14
  * Rhodes et al. (2017) Rhodes, J., Nichol, R. C., Aubourg, É., et al. 2017, ApJS, 233, 21
  * Robotham et al. (2022) Robotham, A. S. G., Bellstedt, S., & Driver, S. P. 2022, ProFuse: Galaxies and components modeler, Astrophysics Source Code Library, record ascl:2204.018
  * Robotham et al. (2018) Robotham, A. S. G., Davies, L. J. M., Driver, S. P., et al. 2018, MNRAS, 476, 3137
  * Robotham et al. (2017) Robotham, A. S. G., Taranu, D. S., Tobar, R., Moffett, A., & Driver, S. P. 2017, MNRAS, 466, 1513
  * Rowe et al. (2015) Rowe, B. T. P., Jarvis, M., Mandelbaum, R., et al. 2015, Astronomy and Computing, 10, 121
  * Schreiber et al. (2017) Schreiber, C., Pannella, M., Leiton, R., et al. 2017, A&A, 599, A134
  * Sérsic (1968) Sérsic, J. L. 1968, Atlas de galaxias australes (Cordoba, Argentina: Observatorio Astronomico, 1968)
  * Smith (1984) Smith, R. L. 1984, Operations Research, 32, 1296
  * Tuccillo et al. (2018) Tuccillo, D., Huertas-Company, M., Decencière, E., et al. 2018, MNRAS, 475, 894
  * Turchin (1971) Turchin, V. F. 1971, Theory of Probability & Its Applications, 16, 720

##  Appendix A Technical information

Here we provide detailed information about the technical realisation of the
noise background and of the RMS maps for the data set.

###  A.1 Noise maps and zero-points

We created Gaussian background noise maps using SkyMaker (Bertin 2009). The
software requires the value of the expected background surface brightness
(${\rm SB}_{\rm bkg}$, in magnitude arcsec-2) and the observational ZP at 1
second of exposure as inputs, to produce a map at the desired depth. The
conceptual steps to compute the ZP, following Euclid Collaboration: Martinet
et al. (2019), are the following. Given a desired limiting magnitude $m_{\rm
lim}$ within a given area $A$ (e.g. in square arcseconds), the corresponding
limiting flux for the total exposure time is

$S=f_{{\rm lim},A}=\frac{t_{\rm exp}}{g}10^{({\rm ZP}-m_{\rm lim})/2.5}\;,$
(3)

where $g$ is the gain of the detector. The uncertainty per pixel from the
background for the same source is

$\sigma_{\rm bkg,pixel}=\sqrt{f_{\rm bkg,pixel}+\sigma_{\rm
RON}^{2}}\simeq\sqrt{l^{2}\frac{t_{\rm exp}}{g}10^{({\rm ZP}-{\rm SB}_{\rm
bkg})/2.5}}\;,$ (4)

where $l$ is the pixel-scale, and for the second step we assume that the
contribution of the read-out noise $\sigma_{\rm RON}$ is negligible. Assuming
that $\sigma_{\rm bkg,pixel}$ is constant, the total noise within the aperture
is

$N=\sqrt{\sum_{{\rm pixels},A}\sigma_{\rm bkg,pixel}^{2}}=\sqrt{A}\sigma_{\rm
bkg,pixel}\;;$ (5)

so if $r$ is the radius of the aperture relative to $A$, $r=\sqrt{A/\pi}$, the
background S/N reads

${\rm S/N}=\frac{f_{{\rm lim},A}}{\sqrt{A}\sigma_{\rm bkg,pixel}}=\frac{t_{\rm
exp}10^{({\rm ZP}-m_{\rm lim})/2.5}}{g\sqrt{\frac{\pi
r^{2}}{l^{2}}}\sqrt{l^{2}\frac{t_{\rm exp}}{g}10^{({\rm ZP}-{\rm SB}_{\rm
bkg})/2.5}}}\;.$ (6)

Finally, the formula can be inverted to make the expression for ZP explicit:

$ZP=2.5\operatorname{log_{10}}\left[\frac{(S/N)^{2}\pi rg}{t_{\rm
exp}}\right]-2m_{\rm lim}-SB_{\rm bkg}\;,$ (7)

where it can be noticed that the pixel-scale terms $l^{2}$ have canceled out.

We computed the ZP of each image putting $\mbox{S/N}=10$ and considering the
corresponding $m_{\rm lim}$ computed in 2$\mathrm{\SIUnitSymbolArcsecond}$
apertures (so $r=1$$\mathrm{\SIUnitSymbolArcsecond}$ in the above formulas),
SBbkg, and exposure times as given in Table 1. We point out again that this
method is used to compute the ZP for 1 second of exposure (to obtain the ZP at
the given exposure time, a term $2.5\operatorname{log_{10}}t_{\rm exp}$ should
be added).

###  A.2 RMS maps and photon noise

To obtain the RMS map for each scientific image, first of all we produced a
map with a constant value equal to the standard deviation of the original
Gaussian background noise image, which only provides information on the noise
due to the unresolved sky background and undetected faint sources.

The flux of the photons arriving from galaxies is a Poissonian distribution
and contributes to the uncertainty on the actual flux estimate. This
contribution (photon noise) was added to the background noise map with the
following procedure. Given an arbitrary exposure time $t$, the pixel values
$S_{t}=C_{t}+{\rm sky}_{t}$ are due to the sum of the counts coming from
sources ($C_{t}$) and from the background (${\rm sky}_{t}$). The noise per
pixel $N_{t}$ is found by summing in quadrature the sky contribution ${\rm
RMS}_{{\rm sky},t}$ and the photon noise contribution ${\rm RMS}_{{\rm
source},t}=\sqrt{C_{t}})$:

$N_{t}=\sqrt{{\rm RMS}^{2}_{{\rm
sky},t}+\left(\sqrt{C_{t}}\right)^{2}}=\sqrt{{\rm RMS}^{2}_{{\rm
sky},t}+C_{t}}\;.$ (8)

Considering a science image normalised to $t=1$ second, in which pixels have
values $S_{1}=S_{t}/t$, with $N_{1}$ being the corresponding RMS map pixel
values, the S/N of each pixel must be conserved, i.e.

$\frac{S_{t}}{N_{t}}=\frac{S_{1}}{N_{1}}\;.$ (9)

This relation leads to the procedure to build the noise map $N_{1}$ of the
normalised image as

$N_{1}=\frac{N_{t}S_{1}}{S_{t}}=\frac{N_{t}}{t}\;.$ (10)

By combining Eqs. 8 and 10, and considering that $C_{t}=C_{1}t$ we finally
obtain

$N_{1}=\sqrt{\frac{{\rm RMS}^{2}_{{\rm sky},t}+C_{t}}{t^{2}}}=\sqrt{{\rm
RMS}^{2}_{{\rm sky},1}+\frac{C_{1}}{t}}\;.$ (11)

##  Appendix B Computational times and separate analysis of the performance of
each software package

In this section we provide information about the computational times and
memory workload required by each software package to complete the runs, and
analyse in more detail some particular cases in which the results were
especially interesting, adding some final remarks where appropriate.

###  B.1 DeepLeGATo

Figure 24: Examples of diagnostic plots for DeepLEGaTo. Top,
$I_{\scriptscriptstyle\rm E}$ SS F4 run; bottom, $I_{\scriptscriptstyle\rm E}$
DS F4 run. See Fig. 8 for a description of this and the subsequent similar
plots. Note that the colour coding is given by the value of the bulge fraction
${\rm b}/{\rm t}$ in the original Egg double-Sersic catalogue; it can be