11institutetext: Institute for Astronomy, University of Edinburgh, Royal
Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK 22institutetext: Department
of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK 33institutetext:
Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy,
University of Manchester, Oxford Road, Manchester M13 9PL, UK 44institutetext:
Mullard Space Science Laboratory, University College London, Holmbury St Mary,
Dorking, Surrey RH5 6NT, UK 55institutetext: Aix-Marseille Université, CNRS,
CNES, LAM, Marseille, France 66institutetext: Universität Innsbruck, Institut
für Astro- und Teilchenphysik, Technikerstr. 25/8, 6020 Innsbruck, Austria
77institutetext: Universität Bonn, Argelander-Institut für Astronomie, Auf dem
Hügel 71, 53121 Bonn, Germany 88institutetext: Université Paris-Saclay, CNRS,
Institut d’astrophysique spatiale, 91405, Orsay, France 99institutetext:
Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth
PO1 3FX, UK 1010institutetext: INAF-Osservatorio Astronomico di Brera, Via
Brera 28, 20122 Milano, Italy 1111institutetext: INAF-Osservatorio di
Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, 40129
Bologna, Italy 1212institutetext: Dipartimento di Fisica e Astronomia,
Università di Bologna, Via Gobetti 93/2, 40129 Bologna, Italy
1313institutetext: INFN-Sezione di Bologna, Viale Berti Pichat 6/2, 40127
Bologna, Italy 1414institutetext: Max Planck Institute for Extraterrestrial
Physics, Giessenbachstr. 1, 85748 Garching, Germany 1515institutetext:
Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-
Universität München, Scheinerstrasse 1, 81679 München, Germany
1616institutetext: INAF-Osservatorio Astrofisico di Torino, Via Osservatorio
20, 10025 Pino Torinese (TO), Italy 1717institutetext: Dipartimento di Fisica,
Università di Genova, Via Dodecaneso 33, 16146, Genova, Italy
1818institutetext: INFN-Sezione di Genova, Via Dodecaneso 33, 16146, Genova,
Italy 1919institutetext: Department of Physics ”E. Pancini”, University
Federico II, Via Cinthia 6, 80126, Napoli, Italy 2020institutetext: INAF-
Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli,
Italy 2121institutetext: INFN section of Naples, Via Cinthia 6, 80126, Napoli,
Italy 2222institutetext: Instituto de Astrofísica e Ciências do Espaço,
Universidade do Porto, CAUP, Rua das Estrelas, PT4150-762 Porto, Portugal
2323institutetext: Dipartimento di Fisica, Università degli Studi di Torino,
Via P. Giuria 1, 10125 Torino, Italy 2424institutetext: INFN-Sezione di
Torino, Via P. Giuria 1, 10125 Torino, Italy 2525institutetext: INAF-IASF
Milano, Via Alfonso Corti 12, 20133 Milano, Italy 2626institutetext: INAF-
Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone,
Italy 2727institutetext: INFN-Sezione di Roma, Piazzale Aldo Moro, 2 - c/o
Dipartimento di Fisica, Edificio G. Marconi, 00185 Roma, Italy
2828institutetext: Institut de Física d’Altes Energies (IFAE), The Barcelona
Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona),
Spain 2929institutetext: Port d’Informació Científica, Campus UAB, C. Albareda
s/n, 08193 Bellaterra (Barcelona), Spain 3030institutetext: Institute for
Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University,
52056 Aachen, Germany 3131institutetext: Institute of Space Sciences (ICE,
CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona, Spain
3232institutetext: Institut d’Estudis Espacials de Catalunya (IEEC), Edifici
RDIT, Campus UPC, 08860 Castelldefels, Barcelona, Spain 3333institutetext:
Dipartimento di Fisica e Astronomia ”Augusto Righi” - Alma Mater Studiorum
Università di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy
3434institutetext: European Space Agency/ESRIN, Largo Galileo Galilei 1, 00044
Frascati, Roma, Italy 3535institutetext: ESAC/ESA, Camino Bajo del Castillo,
s/n., Urb. Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid,
Spain 3636institutetext: Université Claude Bernard Lyon 1, CNRS/IN2P3, IP2I
Lyon, UMR 5822, Villeurbanne, F-69100, France 3737institutetext: Institute of
Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne
(EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland 3838institutetext:
UCB Lyon 1, CNRS/IN2P3, IUF, IP2I Lyon, 4 rue Enrico Fermi, 69622
Villeurbanne, France 3939institutetext: Departamento de Física, Faculdade de
Ciências, Universidade de Lisboa, Edifício C8, Campo Grande, PT1749-016
Lisboa, Portugal 4040institutetext: Instituto de Astrofísica e Ciências do
Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, 1749-016
Lisboa, Portugal 4141institutetext: Department of Astronomy, University of
Geneva, ch. d’Ecogia 16, 1290 Versoix, Switzerland 4242institutetext: INAF-
Istituto di Astrofisica e Planetologia Spaziali, via del Fosso del Cavaliere,
100, 00100 Roma, Italy 4343institutetext: Université Paris-Saclay, Université
Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France 4444institutetext:
Institut de Ciencies de l’Espai (IEEC-CSIC), Campus UAB, Carrer de Can
Magrans, s/n Cerdanyola del Vallés, 08193 Barcelona, Spain 4545institutetext:
INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34143 Trieste,
Italy 4646institutetext: Istituto Nazionale di Fisica Nucleare, Sezione di
Bologna, Via Irnerio 46, 40126 Bologna, Italy 4747institutetext: INAF-
Osservatorio Astronomico di Padova, Via dell’Osservatorio 5, 35122 Padova,
Italy 4848institutetext: Institute of Theoretical Astrophysics, University of
Oslo, P.O. Box 1029 Blindern, 0315 Oslo, Norway 4949institutetext: Higgs
Centre for Theoretical Physics, School of Physics and Astronomy, The
University of Edinburgh, Edinburgh EH9 3FD, UK 5050institutetext: Jet
Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove
Drive, Pasadena, CA, 91109, USA 5151institutetext: von Hoerner & Sulger GmbH,
Schlossplatz 8, 68723 Schwetzingen, Germany 5252institutetext: Technical
University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark
5353institutetext: Cosmic Dawn Center (DAWN), Denmark 5454institutetext:
Institut d’Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98
bis boulevard Arago, 75014 Paris, France 5555institutetext: Max-Planck-
Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
5656institutetext: Department of Physics and Helsinki Institute of Physics,
Gustaf Hällströmin katu 2, 00014 University of Helsinki, Finland
5757institutetext: Aix-Marseille Université, CNRS/IN2P3, CPPM, Marseille,
France 5858institutetext: AIM, CEA, CNRS, Université Paris-Saclay, Université
de Paris, 91191 Gif-sur-Yvette, France 5959institutetext: Leiden Observatory,
Leiden University, Einsteinweg 55, 2333 CC Leiden, The Netherlands
6060institutetext: 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 6161institutetext: Department of Physics, P.O. Box 64, 00014
University of Helsinki, Finland 6262institutetext: Helsinki Institute of
Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland
6363institutetext: NOVA optical infrared instrumentation group at ASTRON, Oude
Hoogeveensedijk 4, 7991PD, Dwingeloo, The Netherlands 6464institutetext:
Dipartimento di Fisica ”Aldo Pontremoli”, Università degli Studi di Milano,
Via Celoria 16, 20133 Milano, Italy 6565institutetext: INFN-Sezione di Milano,
Via Celoria 16, 20133 Milano, Italy 6666institutetext: Dipartimento di Fisica
e Astronomia ”Augusto Righi” - Alma Mater Studiorum Università di Bologna, via
Piero Gobetti 93/2, 40129 Bologna, Italy 6767institutetext: Department of
Physics, Institute for Computational Cosmology, Durham University, South Road,
DH1 3LE, UK 6868institutetext: 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 6969institutetext: Université Paris Cité, CNRS, Astroparticule
et Cosmologie, 75013 Paris, France 7070institutetext: University of Applied
Sciences and Arts of Northwestern Switzerland, School of Engineering, 5210
Windisch, Switzerland 7171institutetext: European Space Agency/ESTEC,
Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands 7272institutetext: School of
Mathematics, Statistics and Physics, Newcastle University, Herschel Building,
Newcastle-upon-Tyne, NE1 7RU, UK 7373institutetext: Department of Physics and
Astronomy, University of Aarhus, Ny Munkegade 120, DK-8000 Aarhus C, Denmark
7474institutetext: Waterloo Centre for Astrophysics, University of Waterloo,
Waterloo, Ontario N2L 3G1, Canada 7575institutetext: Department of Physics and
Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada
7676institutetext: Perimeter Institute for Theoretical Physics, Waterloo,
Ontario N2L 2Y5, Canada 7777institutetext: Université Paris-Saclay, Université
Paris Cité, CEA, CNRS, Astrophysique, Instrumentation et Modélisation Paris-
Saclay, 91191 Gif-sur-Yvette, France 7878institutetext: Space Science Data
Center, Italian Space Agency, via del Politecnico snc, 00133 Roma, Italy
7979institutetext: Centre National d’Etudes Spatiales – Centre spatial de
Toulouse, 18 avenue Edouard Belin, 31401 Toulouse Cedex 9, France
8080institutetext: Institute of Space Science, Str. Atomistilor, nr. 409
Măgurele, Ilfov, 077125, Romania 8181institutetext: Instituto de Astrofísica
de Canarias, Calle Vía Láctea s/n, 38204, San Cristóbal de La Laguna,
Tenerife, Spain 8282institutetext: Departamento de Astrofísica, Universidad de
La Laguna, 38206, La Laguna, Tenerife, Spain 8383institutetext: Dipartimento
di Fisica e Astronomia ”G. Galilei”, Università di Padova, Via Marzolo 8,
35131 Padova, Italy 8484institutetext: INFN-Padova, Via Marzolo 8, 35131
Padova, Italy 8585institutetext: Departamento de Física, FCFM, Universidad de
Chile, Blanco Encalada 2008, Santiago, Chile 8686institutetext: Satlantis,
University Science Park, Sede Bld 48940, Leioa-Bilbao, Spain
8787institutetext: Centro de Investigaciones Energéticas, Medioambientales y
Tecnológicas (CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain
8888institutetext: Instituto de Astrofísica e Ciências do Espaço, Faculdade de
Ciências, Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal
8989institutetext: Universidad Politécnica de Cartagena, Departamento de
Electrónica y Tecnología de Computadoras, Plaza del Hospital 1, 30202
Cartagena, Spain 9090institutetext: Institut de Recherche en Astrophysique et
Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard
Belin, 31400 Toulouse, France 9191institutetext: Kapteyn Astronomical
Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The
Netherlands 9292institutetext: INFN-Bologna, Via Irnerio 46, 40126 Bologna,
Italy 9393institutetext: Infrared Processing and Analysis Center, California
Institute of Technology, Pasadena, CA 91125, USA 9494institutetext: IFPU,
Institute for Fundamental Physics of the Universe, via Beirut 2, 34151
Trieste, Italy 9595institutetext: INAF, Istituto di Radioastronomia, Via Piero
Gobetti 101, 40129 Bologna, Italy 9696institutetext: Department of Mathematics
and Physics E. De Giorgi, University of Salento, Via per Arnesano, CP-I93,
73100, Lecce, Italy 9797institutetext: INAF-Sezione di Lecce, c/o Dipartimento
Matematica e Fisica, Via per Arnesano, 73100, Lecce, Italy 9898institutetext:
INFN, Sezione di Lecce, Via per Arnesano, CP-193, 73100, Lecce, Italy
9999institutetext: Institut für Theoretische Physik, University of Heidelberg,
Philosophenweg 16, 69120 Heidelberg, Germany 100100institutetext: Université
St Joseph; Faculty of Sciences, Beirut, Lebanon 101101institutetext: Institut
d’Astrophysique de Paris, 98bis Boulevard Arago, 75014, Paris, France
102102institutetext: Junia, EPA department, 41 Bd Vauban, 59800 Lille, France
103103institutetext: SISSA, International School for Advanced Studies, Via
Bonomea 265, 34136 Trieste TS, Italy 104104institutetext: INFN, Sezione di
Trieste, Via Valerio 2, 34127 Trieste TS, Italy 105105institutetext: ICSC -
Centro Nazionale di Ricerca in High Performance Computing, Big Data e Quantum
Computing, Via Magnanelli 2, Bologna, Italy 106106institutetext: Instituto de
Física Teórica UAM-CSIC, Campus de Cantoblanco, 28049 Madrid, Spain
107107institutetext: CERCA/ISO, Department of Physics, Case Western Reserve
University, 10900 Euclid Avenue, Cleveland, OH 44106, USA 108108institutetext:
Laboratoire Univers et Théorie, Observatoire de Paris, Université PSL,
Université Paris Cité, CNRS, 92190 Meudon, France 109109institutetext:
Dipartimento di Fisica e Scienze della Terra, Università degli Studi di
Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy 110110institutetext:
Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Via Giuseppe
Saragat 1, 44122 Ferrara, Italy 111111institutetext: Dipartimento di Fisica -
Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, 34131 Trieste,
Italy 112112institutetext: NASA Ames Research Center, Moffett Field, CA 94035,
USA 113113institutetext: Kavli Institute for Particle Astrophysics & Cosmology
(KIPAC), Stanford University, Stanford, CA 94305, USA 114114institutetext: Bay
Area Environmental Research Institute, Moffett Field, California 94035, USA
115115institutetext: Minnesota Institute for Astrophysics, University of
Minnesota, 116 Church St SE, Minneapolis, MN 55455, USA 116116institutetext:
Institute Lorentz, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The
Netherlands 117117institutetext: Institute for Astronomy, University of
Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA 118118institutetext:
Department of Physics & Astronomy, University of California Irvine, Irvine CA
92697, USA 119119institutetext: Department of Astronomy & Physics and
Institute for Computational Astrophysics, Saint Mary’s University, 923 Robie
Street, Halifax, Nova Scotia, B3H 3C3, Canada 120120institutetext:
Departamento Física Aplicada, Universidad Politécnica de Cartagena, Campus
Muralla del Mar, 30202 Cartagena, Murcia, Spain 121121institutetext:
Department of Computer Science, Aalto University, PO Box 15400, Espoo, FI-00
076, Finland 122122institutetext: Ruhr University Bochum, Faculty of Physics
and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological
Lensing (GCCL), 44780 Bochum, Germany 123123institutetext: Université Paris-
Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France 124124institutetext: Univ.
Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, Avenue des Martyrs, 38000,
Grenoble, France 125125institutetext: Department of Physics and Astronomy,
Vesilinnantie 5, 20014 University of Turku, Finland 126126institutetext: Serco
for European Space Agency (ESA), Camino bajo del Castillo, s/n, Urbanizacion
Villafranca del Castillo, Villanueva de la Cañada, 28692 Madrid, Spain
127127institutetext: ARC Centre of Excellence for Dark Matter Particle
Physics, Melbourne, Australia 128128institutetext: Centre for Astrophysics &
Supercomputing, Swinburne University of Technology, Victoria 3122, Australia
129129institutetext: W.M. Keck Observatory, 65-1120 Mamalahoa Hwy, Kamuela,
HI, USA 130130institutetext: Department of Physics and Astronomy, University
of the Western Cape, Bellville, Cape Town, 7535, South Africa
131131institutetext: Dipartimento di Fisica, Sapienza Università di Roma,
Piazzale Aldo Moro 2, 00185 Roma, Italy 132132institutetext: Centro de
Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto,
Portugal 133133institutetext: Zentrum für Astronomie, Universität Heidelberg,
Philosophenweg 12, 69120 Heidelberg, Germany 134134institutetext: Dipartimento
di Fisica, Università di Roma Tor Vergata, Via della Ricerca Scientifica 1,
Roma, Italy 135135institutetext: INFN, Sezione di Roma 2, Via della Ricerca
Scientifica 1, Roma, Italy 136136institutetext: Department of Astrophysics,
University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland
137137institutetext: Department of Astrophysical Sciences, Peyton Hall,
Princeton University, Princeton, NJ 08544, USA 138138institutetext: Niels Bohr
Institute, University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark
139139institutetext: Cosmic Dawn Center (DAWN)

#  Euclid preparation

LensMC, weak lensing cosmic shear measurement with forward modelling and
Markov Chain Monte Carlo sampling

Euclid Collaboration: G. Congedo0000-0003-2508-0046<EMAIL_ADDRESS>Euclid preparationEuclid preparation L. Miller0000-0002-3376-6200 Euclid
preparationEuclid preparation A. N. Taylor Euclid preparationEuclid
preparation N. Cross Euclid preparationEuclid preparation C. A. J. Duncan
Euclid preparationEuclid preparationEuclid preparationEuclid preparation T.
Kitching0000-0002-4061-4598 Euclid preparationEuclid preparation N.
Martinet0000-0003-2786-7790 Euclid preparationEuclid preparation S. Matthew
Euclid preparationEuclid preparation T. Schrabback0000-0002-6987-7834 Euclid
preparationEuclid preparation M. Tewes0000-0002-1155-8689 Euclid
preparationEuclid preparation N. Welikala Euclid preparationEuclid preparation
N. Aghanim Euclid preparationEuclid preparation A. Amara Euclid
preparationEuclid preparation S. Andreon0000-0002-2041-8784 Euclid
preparationEuclid preparation N. Auricchio Euclid preparationEuclid
preparation M. Baldi0000-0003-4145-1943 Euclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparation S. Bardelli0000-0002-8900-0298 Euclid preparationEuclid
preparation R. Bender0000-0001-7179-0626 Euclid preparationEuclid
preparationEuclid preparationEuclid preparation C. Bodendorf Euclid
preparationEuclid preparation D. Bonino Euclid preparationEuclid preparation
E. Branchini0000-0002-0808-6908 Euclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparation M.
Brescia0000-0001-9506-5680 Euclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparation J.
Brinchmann0000-0003-4359-8797 Euclid preparationEuclid preparation S.
Camera0000-0003-3399-3574 Euclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparation V.
Capobianco0000-0002-3309-7692 Euclid preparationEuclid preparation C.
Carbone0000-0003-0125-3563 Euclid preparationEuclid preparation V. F. Cardone
Euclid preparationEuclid preparationEuclid preparationEuclid preparation J.
Carretero0000-0002-3130-0204 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation S. Casas0000-0002-4751-5138 Euclid
preparationEuclid preparation F. J. Castander0000-0001-7316-4573 Euclid
preparationEuclid preparationEuclid preparationEuclid preparation M.
Castellano0000-0001-9875-8263 Euclid preparationEuclid preparation S.
Cavuoti0000-0002-3787-4196 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation A. Cimatti Euclid preparationEuclid preparation
C. J. Conselice Euclid preparationEuclid preparation L.
Conversi0000-0002-6710-8476 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation Y. Copin0000-0002-5317-7518 Euclid
preparationEuclid preparation F. Courbin0000-0003-0758-6510 Euclid
preparationEuclid preparation H. M. Courtois0000-0003-0509-1776 Euclid
preparationEuclid preparation M. Cropper Euclid preparationEuclid preparation
A. Da Silva0000-0002-6385-1609 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation H. Degaudenzi0000-0002-5887-6799 Euclid
preparationEuclid preparation A. M. Di Giorgio0000-0002-4767-2360 Euclid
preparationEuclid preparation J. Dinis Euclid preparationEuclid
preparationEuclid preparationEuclid preparation F. Dubath0000-0002-6533-2810
Euclid preparationEuclid preparation X. Dupac Euclid preparationEuclid
preparation M. Farina Euclid preparationEuclid preparation S.
Farrens0000-0002-9594-9387 Euclid preparationEuclid preparation S. Ferriol
Euclid preparationEuclid preparation P. Fosalba0000-0002-1510-5214 Euclid
preparationEuclid preparationEuclid preparationEuclid preparation M.
Frailis0000-0002-7400-2135 Euclid preparationEuclid preparation E.
Franceschi0000-0002-0585-6591 Euclid preparationEuclid preparation S.
Galeotta0000-0002-3748-5115 Euclid preparationEuclid preparation B.
Garilli0000-0001-7455-8750 Euclid preparationEuclid preparation B.
Gillis0000-0002-4478-1270 Euclid preparationEuclid preparation C.
Giocoli0000-0002-9590-7961 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation A. Grazian0000-0002-5688-0663 Euclid
preparationEuclid preparation F. Grupp Euclid preparationEuclid
preparationEuclid preparationEuclid preparation S. V. H.
Haugan0000-0001-9648-7260 Euclid preparationEuclid preparation M. S. Holliman
Euclid preparationEuclid preparation W. Holmes Euclid preparationEuclid
preparation F. Hormuth Euclid preparationEuclid preparation A.
Hornstrup0000-0002-3363-0936 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation P. Hudelot Euclid preparationEuclid preparation
K. Jahnke0000-0003-3804-2137 Euclid preparationEuclid preparation E.
Keihänen0000-0003-1804-7715 Euclid preparationEuclid preparation S.
Kermiche0000-0002-0302-5735 Euclid preparationEuclid preparation A.
Kiessling0000-0002-2590-1273 Euclid preparationEuclid preparation M.
Kilbinger0000-0001-9513-7138 Euclid preparationEuclid preparation B. Kubik
Euclid preparationEuclid preparation K. Kuijken0000-0002-3827-0175 Euclid
preparationEuclid preparation M. Kümmel0000-0003-2791-2117 Euclid
preparationEuclid preparation M. Kunz0000-0002-3052-7394 Euclid
preparationEuclid preparation H. Kurki-Suonio0000-0002-4618-3063 Euclid
preparationEuclid preparationEuclid preparationEuclid preparation S.
Ligori0000-0003-4172-4606 Euclid preparationEuclid preparation P. B.
Lilje0000-0003-4324-7794 Euclid preparationEuclid preparation V.
Lindholm0000-0003-2317-5471 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation I. Lloro Euclid preparationEuclid preparation D.
Maino Euclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparation E. Maiorano0000-0003-2593-4355
Euclid preparationEuclid preparation O. Mansutti0000-0001-5758-4658 Euclid
preparationEuclid preparation O. Marggraf0000-0001-7242-3852 Euclid
preparationEuclid preparation K. Markovic0000-0001-6764-073X Euclid
preparationEuclid preparation F. Marulli0000-0002-8850-0303 Euclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparation R. Massey0000-0002-6085-3780 Euclid
preparationEuclid preparation S. Maurogordato Euclid preparationEuclid
preparation H. J. McCracken0000-0002-9489-7765 Euclid preparationEuclid
preparation E. Medinaceli0000-0002-4040-7783 Euclid preparationEuclid
preparation S. Mei0000-0002-2849-559X Euclid preparationEuclid preparation M.
Melchior Euclid preparationEuclid preparation M. Meneghetti0000-0003-1225-7084
Euclid preparationEuclid preparationEuclid preparationEuclid preparation E.
Merlin0000-0001-6870-8900 Euclid preparationEuclid preparation G. Meylan
Euclid preparationEuclid preparation M. Moresco0000-0002-7616-7136 Euclid
preparationEuclid preparationEuclid preparationEuclid preparation B. Morin
Euclid preparationEuclid preparation L. Moscardini0000-0002-3473-6716 Euclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparation E. Munari0000-0002-1751-5946 Euclid
preparationEuclid preparation S.-M. Niemi Euclid preparationEuclid preparation
J. W. Nightingale0000-0002-8987-7401 Euclid preparationEuclid
preparationEuclid preparationEuclid preparation C. Padilla0000-0001-7951-0166
Euclid preparationEuclid preparation S. Paltani Euclid preparationEuclid
preparation F. Pasian Euclid preparationEuclid preparation K. Pedersen Euclid
preparationEuclid preparation W. J. Percival0000-0002-0644-5727 Euclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparation V. Pettorino Euclid preparationEuclid
preparation S. Pires0000-0002-0249-2104 Euclid preparationEuclid preparation
G. Polenta0000-0003-4067-9196 Euclid preparationEuclid preparation M. Poncet
Euclid preparationEuclid preparation L. A. Popa Euclid preparationEuclid
preparation L. Pozzetti0000-0001-7085-0412 Euclid preparationEuclid
preparation F. Raison0000-0002-7819-6918 Euclid preparationEuclid preparation
R. Rebolo Euclid preparationEuclid preparationEuclid preparationEuclid
preparation A. Renzi0000-0001-9856-1970 Euclid preparationEuclid
preparationEuclid preparationEuclid preparation J. Rhodes Euclid
preparationEuclid preparation G. Riccio Euclid preparationEuclid preparation
E. Romelli0000-0003-3069-9222 Euclid preparationEuclid preparation M.
Roncarelli0000-0001-9587-7822 Euclid preparationEuclid preparation E. Rossetti
Euclid preparationEuclid preparation R. Saglia0000-0003-0378-7032 Euclid
preparationEuclid preparationEuclid preparationEuclid preparation D.
Sapone0000-0001-7089-4503 Euclid preparationEuclid preparation B. Sartoris
Euclid preparationEuclid preparationEuclid preparationEuclid preparation P.
Schneider0000-0001-8561-2679 Euclid preparationEuclid preparation A.
Secroun0000-0003-0505-3710 Euclid preparationEuclid preparation G.
Seidel0000-0003-2907-353X Euclid preparationEuclid preparation S.
Serrano0000-0002-0211-2861 Euclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparation C.
Sirignano0000-0002-0995-7146 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation G. Sirri0000-0003-2626-2853 Euclid
preparationEuclid preparation L. Stanco0000-0002-9706-5104 Euclid
preparationEuclid preparation P. Tallada-Crespí0000-0002-1336-8328 Euclid
preparationEuclid preparationEuclid preparationEuclid preparation D.
Tavagnacco0000-0001-7475-9894 Euclid preparationEuclid preparation I. Tereno
Euclid preparationEuclid preparationEuclid preparationEuclid preparation R.
Toledo-Moreo0000-0002-2997-4859 Euclid preparationEuclid preparation F.
Torradeflot0000-0003-1160-1517 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation I. Tutusaus0000-0002-3199-0399 Euclid
preparationEuclid preparation E. A. Valentijn Euclid preparationEuclid
preparation L. Valenziano0000-0002-1170-0104 Euclid preparationEuclid
preparationEuclid preparationEuclid preparation T. Vassallo0000-0001-6512-6358
Euclid preparationEuclid preparationEuclid preparationEuclid preparation A.
Veropalumbo0000-0003-2387-1194 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation Y. Wang0000-0002-4749-2984 Euclid
preparationEuclid preparation J. Weller0000-0002-8282-2010 Euclid
preparationEuclid preparationEuclid preparationEuclid preparation G.
Zamorani0000-0002-2318-301X Euclid preparationEuclid preparation J. Zoubian
Euclid preparationEuclid preparation E. Zucca0000-0002-5845-8132 Euclid
preparationEuclid preparation A. Biviano0000-0002-0857-0732 Euclid
preparationEuclid preparationEuclid preparationEuclid preparation M.
Bolzonella0000-0003-3278-4607 Euclid preparationEuclid preparation A.
Boucaud0000-0001-7387-2633 Euclid preparationEuclid preparation E.
Bozzo0000-0002-8201-1525 Euclid preparationEuclid preparation C.
Burigana0000-0002-3005-5796 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation C. Colodro-Conde Euclid preparationEuclid
preparation D. Di Ferdinando Euclid preparationEuclid preparation J. Graciá-
Carpio Euclid preparationEuclid preparation N. Mauri0000-0001-8196-1548 Euclid
preparationEuclid preparationEuclid preparationEuclid preparation C. Neissner
Euclid preparationEuclid preparationEuclid preparationEuclid preparation A. A.
Nucita Euclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparation Z. Sakr0000-0002-4823-3757
Euclid preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparation V. Scottez Euclid preparationEuclid
preparationEuclid preparationEuclid preparation M. Tenti0000-0002-4254-5901
Euclid preparationEuclid preparation M. Viel0000-0002-2642-5707 Euclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparation M. Wiesmann Euclid preparationEuclid preparation
Y. Akrami0000-0002-2407-7956 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation V. Allevato Euclid preparationEuclid preparation
S. Anselmi0000-0002-3579-9583 Euclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparation C.
Baccigalupi0000-0002-8211-1630 Euclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparation M. Ballardini Euclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparation S. Borgani0000-0001-6151-6439 Euclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparation A. S.
Borlaff0000-0003-3249-4431 Euclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparation S.
Bruton0000-0002-6503-5218 Euclid preparationEuclid preparation R.
Cabanac0000-0001-6679-2600 Euclid preparationEuclid preparation A. Cappi
Euclid preparationEuclid preparationEuclid preparationEuclid preparation C. S.
Carvalho Euclid preparationEuclid preparation G. Castignani0000-0001-6831-0687
Euclid preparationEuclid preparationEuclid preparationEuclid preparation T.
Castro0000-0002-6292-3228 Euclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparation G. Cañas-Herrera0000-0003-2796-2149 Euclid
preparationEuclid preparationEuclid preparationEuclid preparation K. C.
Chambers0000-0001-6965-7789 Euclid preparationEuclid preparation A. R.
Cooray0000-0002-3892-0190 Euclid preparationEuclid preparation J. Coupon
Euclid preparationEuclid preparation S. Davini Euclid preparationEuclid
preparation G. De Lucia0000-0002-6220-9104 Euclid preparationEuclid
preparation G. Desprez Euclid preparationEuclid preparation S. Di
Domizio0000-0003-2863-5895 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation H. Dole0000-0002-9767-3839 Euclid
preparationEuclid preparation A. Díaz-Sánchez0000-0003-0748-4768 Euclid
preparationEuclid preparation J. A. Escartin Vigo Euclid preparationEuclid
preparation S. Escoffier0000-0002-2847-7498 Euclid preparationEuclid
preparation I. Ferrero0000-0002-1295-1132 Euclid preparationEuclid preparation
F. Finelli0000-0002-6694-3269 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation L. Gabarra Euclid preparationEuclid
preparationEuclid preparationEuclid preparation J. García-
Bellido0000-0002-9370-8360 Euclid preparationEuclid preparation E.
Gaztanaga0000-0001-9632-0815 Euclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparation F.
Giacomini0000-0002-3129-2814 Euclid preparationEuclid preparation G.
Gozaliasl0000-0002-0236-919X Euclid preparationEuclid preparationEuclid
preparationEuclid preparation D. Guinet0000-0002-8132-6509 Euclid
preparationEuclid preparation A. Hall0000-0002-3139-8651 Euclid
preparationEuclid preparation H. Hildebrandt0000-0002-9814-3338 Euclid
preparationEuclid preparation S. Ilić0000-0003-4285-9086 Euclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparation A. Jimenez Muñoz0009-0004-5252-185X Euclid
preparationEuclid preparation S. Joudaki Euclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparation J. J. E. Kajava0000-0002-3010-8333 Euclid preparationEuclid
preparationEuclid preparationEuclid preparation V. Kansal Euclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparation D. Karagiannis Euclid preparationEuclid
preparation C. C. Kirkpatrick Euclid preparationEuclid preparation L.
Legrand0000-0003-0610-5252 Euclid preparationEuclid preparation J. Macias-
Perez Euclid preparationEuclid preparation G. Maggio0000-0003-4020-4836 Euclid
preparationEuclid preparation M. Magliocchetti0000-0001-9158-4838 Euclid
preparationEuclid preparation R. Maoli0000-0002-6065-3025 Euclid
preparationEuclid preparationEuclid preparationEuclid preparation M.
Martinelli0000-0002-6943-7732 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation C. J. A. P. Martins0000-0002-4886-9261 Euclid
preparationEuclid preparationEuclid preparationEuclid preparation M.
Maturi0000-0002-3517-2422 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation L. Maurin0000-0002-8406-0857 Euclid
preparationEuclid preparation R. B. Metcalf0000-0003-3167-2574 Euclid
preparationEuclid preparationEuclid preparationEuclid preparation M.
Migliaccio Euclid preparationEuclid preparationEuclid preparationEuclid
preparation P. Monaco Euclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparationEuclid preparation G. Morgante Euclid preparationEuclid preparation
S. Nadathur0000-0001-9070-3102 Euclid preparationEuclid preparation L.
Patrizii Euclid preparationEuclid preparation A. Peel0000-0003-0488-8978
Euclid preparationEuclid preparation A. Pezzotta Euclid preparationEuclid
preparation V. Popa Euclid preparationEuclid preparation C.
Porciani0000-0002-7797-2508 Euclid preparationEuclid preparation D.
Potter0000-0002-0757-5195 Euclid preparationEuclid preparation M.
Pöntinen0000-0001-5442-2530 Euclid preparationEuclid preparation P.
Reimberg0000-0003-3410-0280 Euclid preparationEuclid preparation P.-F. Rocci
Euclid preparationEuclid preparation A. G. Sánchez0000-0003-1198-831X Euclid
preparationEuclid preparation J. A. Schewtschenko Euclid preparationEuclid
preparation A. Schneider0000-0001-7055-8104 Euclid preparationEuclid
preparation E. Sefusatti0000-0003-0473-1567 Euclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparationEuclid
preparation M. Sereno0000-0003-0302-0325 Euclid preparationEuclid
preparationEuclid preparationEuclid preparation P. Simon Euclid
preparationEuclid preparation A. Spurio Mancini0000-0001-5698-0990 Euclid
preparationEuclid preparation J. Stadel0000-0001-7565-8622 Euclid
preparationEuclid preparation J. Steinwagner Euclid preparationEuclid
preparation G. Testera Euclid preparationEuclid preparation R.
Teyssier0000-0001-7689-0933 Euclid preparationEuclid preparation S.
Toft0000-0003-3631-7176 Euclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparation S.
Tosi0000-0002-7275-9193 Euclid preparationEuclid preparationEuclid
preparationEuclid preparationEuclid preparationEuclid preparation A.
Troja0000-0003-0239-4595 Euclid preparationEuclid preparationEuclid
preparationEuclid preparation M. Tucci Euclid preparationEuclid preparation C.
Valieri Euclid preparationEuclid preparation J. Valiviita0000-0001-6225-3693
Euclid preparationEuclid preparationEuclid preparationEuclid preparation D.
Vergani0000-0003-0898-2216 Euclid preparationEuclid preparation

LensMC is a weak lensing shear measurement method developed for Euclid and
Stage-IV surveys. It is based on forward modelling to deal with convolution by
a point spread function with comparable size to many galaxies; sampling the
posterior distribution of galaxy parameters via Markov Chain Monte Carlo; and
marginalisation over nuisance parameters for each of the 1.5 billion galaxies
observed by Euclid. The scientific performance is quantified through high-
fidelity images based on the Euclid Flagship simulations and emulation of the
Euclid VIS images; realistic clustering with a mean surface number density of
$\qty{250}{arcmin^{-2}}$ ($I_{\scriptscriptstyle\rm E}<29.5$) for galaxies,
and $\qty{6}{arcmin^{-2}}$ ($I_{\scriptscriptstyle\rm E}<26$) for stars; and a
diffraction-limited chromatic point spread function with a full width at half
maximum of $None$ and spatial variation across the field of view. Objects are
measured with a density of $\qty{90}{arcmin^{-2}}$ ($I_{\scriptscriptstyle\rm
E}<26.5$) in $4500\,\unit{deg^{2}}$. The total shear bias is broken down into
measurement (our main focus here) and selection effects (which will be
addressed elsewhere). We find: measurement multiplicative and additive biases
of $m_{1}=(-3.6\pm 0.2)\times 10^{-3}$, $m_{2}=(-4.3\pm 0.2)\times 10^{-3}$,
$c_{1}=(-1.78\pm 0.03)\times 10^{-4}$, $c_{2}=(0.09\pm 0.03)\times 10^{-4}$; a
large detection bias with a multiplicative component of $1.2\times 10^{-2}$
and an additive component of $-3\times 10^{-4}$; and a measurement PSF leakage
of $\alpha_{1}=(-9\pm 3)\times 10^{-4}$ and $\alpha_{2}=(2\pm 3)\times
10^{-4}$. When model bias is suppressed, the obtained measurement biases are
close to Euclid requirement and largely dominated by undetected faint galaxies
($-5\times 10^{-3}$). Although significant, model bias will be straightforward
to calibrate given the weak sensitivity.

###### Key Words.:

Gravitational lensing: weak; Cosmology: observations; Methods: data analysis

##  1 Introduction

Weak gravitational lensing by large-scale structure is a mature cosmological
tool to measure its distribution of dark matter and study dark energy through
the evolution with redshift (Schneider, 2006; Kilbinger, 2015; Mandelbaum,
2018). Weak lensing is particularly sensitive to modifications of the theory
of gravity or the emergence of physics beyond the concordance $\Lambda$-cold
dark matter model, which affect the clustering of dark matter (Amendola et
al., 2018).

Galaxies surveys from the ground, such as the Dark Energy Survey (DES, Abbott
et al., 2018b), the Kilo Degree Survey (KiDS, Kuijken et al., 2015), and the
Hyper Suprime-Cam survey (HSC, Aihara et al., 2017), are now achieving
constraints on the dark matter sector (primarily in the
$\Omega_{\text{m}}$–$\sigma_{8}$ parameter space and their combination
$S_{8}$) to a few percent (Abbott et al., 2018a; Hildebrandt et al., 2017;
Hikage et al., 2019). After extensive consistency checks and sensitivity
studies, recent lensing measurements from galaxy surveys have been shown to be
broadly in agreement with each other, but in mild tension with the Planck
satellite at the $2\,\sigma$ level or greater (Aghanim et al., 2020; Joudaki
et al., 2020; Amon et al., 2022; Heymans et al., 2021; Asgari et al., 2021;
Loureiro et al., 2022), including the latest joint analysis of DES and KiDS
(Abbott et al., 2023).

In the coming years, galaxy surveys will enter a new regime of area, depth,
and image quality. The space-based Euclid telescope (survey area of
$14\,000\,\unit{\deg^{2}}$, full width at half maximum resolution of $None$,
depth of 24.5, $I_{\scriptscriptstyle\rm E}$+$Y_{\scriptscriptstyle\rm
E}J_{\scriptscriptstyle\rm E}H_{\scriptscriptstyle\rm E}$ filters, see
Laureijs et al., 2011; Cropper et al., 2016; Euclid Collaboration: Scaramella
et al., 2022); the planned space-based Roman telescope
($1700\,\unit{\deg^{2}}$, $None$, 26.5, YJH+F184, see Spergel et al., 2015);
and the ground-based Rubin Observatory Legacy Survey of Space and Time (LSST,
$18\,000\,\unit{\deg^{2}}$, $None$, 27.5, $ugrizy$, see Ivezić et al., 2019)
will substantially increase the number of observable galaxies compared to
current surveys. The systematic observation of a billion galaxies or more
across one third of the visible sky will then be possible for the first time.
The combined effect of improved survey area and angular resolution will be an
enhanced ability to probe both the large and small scales via weak lensing and
galaxy clustering, allowing us to constrain cosmological models and dark
energy to percent level precision (Mandelbaum et al., 2018; Euclid
Collaboration: Blanchard et al., 2020), or even an order of magnitude better
when combined with data from Planck (Euclid Collaboration: Ilić et al., 2022).

With such a dramatic improvement in precision that will be achieved in the
coming years, experiments are now focussing on understanding the accuracy of
their analyses (see, e.g., Euclid Collaboration: Paykari et al., 2020). Along
with theory uncertainties, the cosmic shear measurement and redshift
estimation are the most challenging aspects of any large-scale weak lensing
surveys. The concern of this paper is on cosmic shear measurement, which
provides the necessary data for the weak lensing cosmological analysis
(Kilbinger, 2015). In order to achieve of order one percent precision on the
dark energy equation of state, a billion galaxies or more with median redshift
around $1$ need to be observed. This observation has to be carried out
consistently so the same shape measurement procedure is applied to all
objects. This measurement has to be conducted with outstanding accuracy to
satisfy the stringent requirement of $2\times 10^{-3}$ and $3\times 10^{-4}$
on the measured multiplicative and additive shear biases that were set in the
early development phase of weak lensing space telescopes (Massey et al., 2012;
Cropper et al., 2013).

Throughout the years, a number of shear measurement methods have been
developed, tested on data challenges, and applied to real data. These can be
categorized in two main classes: non-parametric and parametric. Among the non-
parametric is Kaiser-Squires-Broadhurst (KSB, Kaiser et al., 1995) based on
weighted moments of image data. Because of its simplicity, these estimators
were used for the very first attempts at measuring cosmic shear in the early
2000s. These methods are fast, so can be quickly calibrated, but sensitive to
effects which need to be characterised. Later on, with better precision, it
became clear that more effects needed to be factored in, particularly a
realistic point spread function (PSF) and the sensitivity of bias on PSF
ellipticity, known as leakage. Parametric methods, based on forward modelling
and model fitting, soon appeared to be better suited to accurately incorporate
such real data features building on solid statistical grounds. In recent
years, a more systematic use of model fitting techniques has been observed
across all major lensing surveys. The Bayesian-inspired shape method lensfit
(Miller et al., 2007; Kitching et al., 2008; Miller et al., 2013) has been
extremely successful, first in the Canada-France Hawaii Telescope Lensing
Survey (CFHTLenS, Heymans et al., 2012), and then more recently also KiDS.
This is based on forward modelling and marginalisation over galaxy nuisance
parameters. A similar method is IM3SHAPE (Zuntz et al., 2013), a maximum
likelihood estimator based again on analytic forward models that has been
applied to DES. While lots of real-data effects, including the PSF, are
accounted for and can directly be built in parametric methods, any in-built
correction clearly introduces extra computational overhead as the parameter
probability distribution needs to be sampled accurately.

While lensing measurements have become more precise over time, accuracy has
also needed to be examined more carefully. Methods have been compared in data
challenges and have run on common simulations with increasing realism, such as
in the Shear TEsting Programme (STEP, Heymans et al., 2006; Massey et al.,
2007) and the Gravitational LEnsing Accuracy Testing (GREAT, Bridle et al.,
2010; Kitching et al., 2012; Mandelbaum et al., 2015). With no method
outshining in absolute terms, and methods being better at some aspects of the
measurement but worse at others, it has become evident that some form of
calibration is still necessary. The field has become reliant on galaxy
simulations more than ever. Sophisticated high-fidelity simulations now need
to reproduce the realism of actual observations as close as possible, so all
biases from detection, measurement, and selection can be fully captured
(Fenech Conti et al., 2017; Kannawadi et al., 2019; Euclid Collaboration:
Martinet et al., 2019; MacCrann et al., 2021; Li et al., 2023a, b).
Calibration naturally raises the question about how sensitive our results are
to the assumptions we make in our simulations (Hoekstra et al., 2017), or how
large these simulations need to be to meet the desired precision (Pujol et
al., 2019; Jansen et al., 2024). Other methods now rely on some form of
calibration that is directly built in the measurement process. Galaxy images
used in calibration are either simulated internally, inferred from real data,
or a combination of the two. Metacalibration (Sheldon & Huff, 2017; Huff &
Mandelbaum, 2017) derives internal estimates of the sensitivity of the
ellipticity estimator to input shears and has been extremely successful on DES
Year 3 (Gatti et al., 2021); Bayesian Fourier Domain (BFD, Bernstein &
Armstrong, 2014; Bernstein et al., 2016) estimates the Taylor coefficients of
the galaxy likelihood expanded over shear with information about moments
measured from calibration fields; a similar implementation to BFD uses forward
modelling (Sheldon, 2014); the KiDS self-calibration (Fenech Conti et al.,
2017) derives internal estimates of the ellipticity bias from noise-free
galaxy images; MomentsML relies on simulated images to train shear-predicting
artificial neural networks (Tewes et al., 2019). Because many selection biases
happen before the shear measurement introduces its own bias, the field has
gradually become more aware that those will probably be the limiting factor in
future lensing surveys. Further work around Metacalibration has led to
Metadetection to address the issue (Sheldon et al., 2020). An application of
the method to Euclid-like simulations has shown that while selection biases
may exceed requirements, the outlook is still positive with demonstrated
success at handling detection and blending biases (Hoekstra, 2021; Hoekstra et
al., 2021; Melchior et al., 2021).

The impact of neighbours to lensing measurements has also become one of the
most important issues that current and future surveys will need to address. In
space, the large number density of detected galaxies of about
$\qty{30}{arcmin^{-2}}$ ($I_{\scriptscriptstyle\rm E}<24.5$) is compensated by
a good image resolution, so the impact of neighbours may not be as bad as on
ground. In fact, due to the worse resolution on the ground, the impact of
neighbours is serious, affecting $60\%$ of the sample (Bosch et al., 2017;
Arcelin et al., 2020). DES Year 1 results (Samuroff et al., 2017) showed that
the total neighbour bias can be as large as $9\%$ (reaching $80\%$ at a close
distance to the neighbour, if uncorrected). Cuts to the catalogue to remove
those objects can reduce the total bias to $1\%$ (reaching $30\%$ at a close
distance to the neighbour), however at the cost of reducing the effective
number density by $30\%$ and leaving a residual bias on $S_{8}$ of
$2\,\sigma$. While Metacalibration has been extremely successful in a few
idealised cases (e.g. isolated galaxies) and Metadetection in the handling of
blending and detection bias, a suite of advanced simulations have been
required for the tomographic calibration of DES Year 3 (MacCrann et al.,
2021). These have shown that an external calibration is still required, as the
unresolved neighbour introduces correlation between the two galaxies at
different redshift, plus the Metacalibration shear responses will be biased by
the presence of the neighbour itself. Therefore calibration simulations have
been made necessary to correctly capture neighbour bias and the interplay
between shear and redshift. However, these simulations assume uniform random
distributions of galaxies that have been reweighted to mimic clustering, which
raises the question that the inferred bias may likely be underestimated. The
most recent simulations by KiDS have realistic clustering of $N$-body
simulations, mimic a number of measurement effects, and address the shear-
redshift interplay (Li et al., 2023a). With a larger number density it is
expected that the situation may be more serious in future surveys.

In this paper we present our advanced shear measurement method LensMC,
specifically developed for Euclid, that builds upon the knowledge and success
of ground-based measurement at handling real data effects. Similarly to
lensfit, it adopts a mean estimator. Contrarily to lensfit, it does not
marginalise over nuisance parameters with numerical approximations. With full
flexibility in the choice of the prior, all the marginalisation in LensMC is
performed by Markov Chain Monte Carlo (MCMC) sampling, for individual galaxies
or jointly across groups of neighbouring galaxies. While IM3SHAPE returns the
maximum of the likelihood estimated via the Levenberg-Marquardt algorithm,
LensMC employs a combination of large-scale and small-scale algorithms,
conjugate gradient and simplex methods, to estimate a suitable initial guess
for the MCMC sampling, thus requiring no information about the model
derivatives and dramatically reducing the sensitivity on the initial guess
(which is assumed always the same for all galaxies). The galaxy models in
LensMC are rendered directly in the Fourier space, hence only a single Fourier
transform is required. A recent profile-fitting method, The Farmer, has also
drawn attention recently (Weaver et al., 2023). This method is a maximum
likelihood estimator, whose position and shape initial guesses are provided by
the detection method. It includes a decision tree based on $\chi^{2}$ values
to classify objects on their likely type and provides error estimates via
Cramer-Rao bounds. Preliminary results are encouraging, however the method has
not been tested to full space-based cosmic shear accuracy yet.

Accurate cosmic shear measurement requires controlling the bias from a number
of sources. Key examples are: source clustering, faint objects, neighbours,
PSF leakage, image artefacts, and cosmic rays. Additionally, any forward-
modelling methods is plagued by potential model bias. One of the main sources
of model bias is addressed here. In summary, LensMC employs: (i) forward
modelling to deal with Euclid image undersampling and convolution by a PSF
with comparable size to many galaxies; (ii) joint measurement of object groups
to correctly handle bias due to neighbours; (iii) masking out objects
belonging to different groups; (iv) MCMC sampling to sample the posterior in a
multi-dimensional parameter space, calculate weights, and correctly
marginalise ellipticity over nuisance parameters and other objects in the same
group. We particularly focus on the realism of our simulations, including
clustering, stars, object detection, the handling of neighbours due to the
high number density, and the use of realistic undersampled chromatic PSF
images with spatial variation across the field of view. We do not include
further real data effects such as non-linearities or cosmic rays as these will
be addressed separately. Also we assume the same broadband PSF as obtained
from a spectral energy distribution (SED) of an SBc-type galaxy at redshift of
1 in both simulations and measurements.

Sect. 2 introduces our method and the practical advantages in addressing real
data problems. Our forward-modelling method is sufficiently fast to analyse
the typical data volume of Stage-IV surveys and can be applied to the
complexity of Euclid measurements, including undersampled data and a complex
PSF, while accounting for the full degrees of freedom in the galaxy modelling.
Additionally, it allows the proper handling of resolved neighbours by joint
measurement and masking of more distant galaxies, stars, and artefacts in the
images. Sect. 3 describes the simulations used for our intensive testing of
the method. The images are fiducial realisations of the Euclid VIS detector
and galaxies are rendered based on $N$-body simulations with full variability
of the morphological properties. All galaxies are convolved with a realistic
pre-flight PSF model with full spatial variation, but ignoring the chromatic
variability. Sect. 4 presents the main results of this testing. When model
bias, chromaticity, and selection biases are suppressed, the obtained biases
are close to Euclid requirement. This measurement bias is largely dominated by
undetected faint galaxies in the images. The bias is also found to be stable
and mostly insensitive to the many effects in the simulations. As the Euclid
analysis will also need to correct for other artefacts in the images, the
residual bias will be straightforward to calibrate through image simulations.
Once we include the model bias by allowing the full variability in the galaxy
models, the overall bias becomes significant. However, since the sensitivity
is weak (the derivative of the bias with respect to the assumptions made in
the simulations appears negligible), it will then be straightforward to also
calibrate the model bias through image simulations. Sect. 5 discusses the main
findings and draws the conclusions of our work.

##  2 Method

The main physical quantity of interest in weak lensing is the reduced cosmic
shear (Kilbinger, 2015),

$g=\frac{\gamma}{1-\kappa}\leavevmode\nobreak\ ,$ (1)

where $\kappa$ and $\gamma$ are convergence and shear (both related to the
gravitational potential), and $g\approx\gamma$ in the weak lensing regime. The
related observable in weak lensing is the ellipticity of a galaxy,

$\epsilon=\frac{a-b}{a+b}\mathrm{e}^{2\mathrm{i}\varphi}\leavevmode\nobreak\
,$ (2)

where $a$ and $b$ are, respectively, the semi-major and semi-minor
axes,111This holds true only for elliptical isophotos, but the ellipticity
remains well-defined if one specifies how it is measured, i.e., it becomes
method dependent. $\varphi\in\left[0,\pi\right)$ is the orientation angle of
the galaxy, and $|\epsilon|\leq 1$. The effect of weak lensing is to distort
the ellipticity of a source galaxy, $\epsilon_{\text{s}}$, by the canonical
transformation (Seitz & Schneider, 1997),

$\epsilon=\frac{\epsilon_{\text{s}}+g}{1+\epsilon_{\text{s}}\,g^{*}}\leavevmode\nobreak\
,$ (3)

where all spin-2 quantities are expressed in complex notation, e.g.,
$\epsilon=\epsilon_{1}+\mathrm{i}\,\epsilon_{2}$, where $\epsilon_{1}$
quantifies the distortion along $x$ and $y$, and $\epsilon_{2}$ along the
coordinate axes rotated by $\pi/4$. As it is customary in weak lensing, we
will refer to $\epsilon_{\text{s}}$ as the intrinsic ellipticity of the source
galaxy, and $\epsilon$ as the lensed or observed ellipticity. The ellipticity
in Eq. (3) is a point estimate for shear in that information on the underlying
cosmic shear can be derived in a statistical sense as a sample average,
$g=\langle\epsilon\rangle$, which holds whenever the distribution of
orientation angles is uniform, e.g., when there are no astrophysical intrinsic
alignments (Joachimi et al., 2015) or shear dependent selection effects.

In weak lensing measurements we infer the reduced shear through sample
averages. In this work, we use the ellipticity as a point estimator for shear
and the problem of measuring ellipticity can be formulated fully in Bayesian
terms. Suppose we have a pixel data vector, $\boldsymbol{D}$, and a model for
the galaxy brightness distribution,
$\boldsymbol{I}=\boldsymbol{I}(\epsilon,\theta,\phi)$, as a function of
ellipticity, $\epsilon$, intrinsic nuisance parameters, $\theta$, and linear
nuisance parameters, $\phi$.222The parameter vectors $\theta$ and $\phi$ are
not to be confused with angular coordinates. Here $\theta$ represents non-
linear parameters (such as object size and position offsets) and $\phi$
represents linear parameters (such as component fluxes) that can be
analytically marginalised out. Thanks to Bayes’ theorem, we can define a joint
posterior as follows,

$p(\epsilon,\theta,\phi|\boldsymbol{D})=\frac{p(\boldsymbol{D}|\epsilon,\theta,\phi)p(\epsilon,\theta,\phi)}{p(\boldsymbol{D})}\leavevmode\nobreak\
,$ (4)

where $p(\boldsymbol{D}|\epsilon,\theta,\phi)$ is the likelihood,
$p(\epsilon,\theta,\phi)$ is the prior on ellipticity and nuisance parameters,
and $p(\boldsymbol{D})$ is the evidence or marginal likelihood,

$p(\boldsymbol{D})=\int
p(\boldsymbol{D}|\epsilon,\theta,\phi)p(\epsilon,\theta,\phi)\,\text{d}\epsilon\,\text{d}\theta\,\text{d}\phi\leavevmode\nobreak\
.$ (5)

We can then construct the ellipticity marginal posterior

$p(\epsilon|\boldsymbol{D})=\frac{1}{p(\boldsymbol{D})}\int
p(\boldsymbol{D}|\epsilon,\theta,\phi)p(\epsilon,\theta,\phi)\,\text{d}\theta\,\text{d}\phi\leavevmode\nobreak\
,$ (6)

marginalising out the nuisance parameters. Common choices of estimators are
the maximum likelihood or maximum posterior probability, but these are usually
biased if the distribution is not Gaussian. However, the bias can be predicted
in simple cases of low dimensionality or when the probability function is
fully analytic (Cox & Snell, 1968; Hall & Taylor, 2017). Another option that
has been successful in ground-based surveys (Miller et al., 2013) is to set
our estimator to the mean of the posterior distribution,

$\hat{\epsilon}=\int\epsilon\,p(\epsilon|\boldsymbol{D})\,\text{d}\epsilon\leavevmode\nobreak\
.$ (7)

We will adopt this definition as it has some useful properties:

  1. 1.

as the nuisance parameters are marginalised out, their impact on the
ellipticity estimator via their correlation is mitigated;

  2. 2.

overfitting333This is the tendency of some estimators, in particular the
maximum of the probability distribution function, to interpret random
fluctuations in the noise as actual signal in the data. is inherently reduced
as we pick an average representative of all possible realisations that are
statistically equivalent;

  3. 3.

following on from the previous point, we expect the mean estimator to be, in
general, less biased than the maximum estimators;

  4. 4.

the mean of the distribution can be estimated through MCMC sampling techniques
(see Sect. 2.3); such estimators satisfy the central limit theorem and
therefore converge to the true mean.

We will discuss more about those points later in this section. Whatever choice
is made, any estimator can be seen as a non-linear mapping between
$\boldsymbol{D}$ and $\epsilon$. Therefore even if $\boldsymbol{D}$ were to be
Gaussian distributed, the estimator will not, hence leading to a fundamental
bias in the measurement, which has been commonly referred to as noise bias
(Melchior & Viola, 2012; Refregier et al., 2012; Viola et al., 2014).
Moreover, as the shear is estimated through an sample average over a
population of galaxies with varying morphological properties and complex
selections, the properties of the shear bias will be different from that of
galaxy ellipticity. Assuming shear is small, it is customary in the field to
model the shear bias on each component with a linear model (Guzik & Bernstein,
2005; Huterer et al., 2006; Heymans et al., 2006),

$\hat{g}_{i}=(1+m_{i})\,g_{i}+c_{i}+n_{i}\leavevmode\nobreak\ ,$ (8)

where $m_{i}$ and $c_{i}$ are the multiplicative and additive biases for the
$i$-th component, $g_{i}$ is the true shear, $\hat{g}_{i}$ is an estimate of
it, and $n_{i}$ is the corresponding statistical noise. The transformation in
Eq. (8) should in principle have $m_{i}$ replaced by a matrix $m_{ij}$ to
model any potential cross-talk between shear components. Alternatively, it
could be rewritten as a spin-2 equation (Kitching & Deshpande, 2022),

$\hat{g}=(1+m_{0})\,g+m_{4}\,g^{*}+c+n\leavevmode\nobreak\ .$ (9)

However, generalising upon previous work, $m_{0}$ and $m_{4}$ are now spin 0
and 4 complex operators, where $m_{0}=|m_{0}|\exp[2\,i\,\sigma_{0}]$,
$m_{4}=|m_{4}|\exp[4\,i\,\sigma_{4}]$, and $\sigma_{0}$ and $\sigma_{4}$ are
angles. Physically, this added flexibility allows for complete mode-mixing:
$m_{0}$ models a dilation and rotation of the true shear, whereas $m_{4}$
models a reflection around the axis determined by its phase. We defer the
application of this approach to future work. Multiplicative terms can be
induced by, e.g., non-Gaussianities in the posterior (skewness at first
order), caused by, e.g., pixel noise and a small galaxy size relative to the
PSF. Additive terms are due to anisotropies induced by, e.g., the PSF and its
spatial variability across the field of view. This effect is referred to as
PSF leakage and is defined by the dependence of $c_{i}$ on the PSF ellipticity
$\epsilon_{\text{PSF},i}$ (see, e.g., Gatti et al., 2021, and references
therein),

$\alpha_{i}=\frac{\,\text{d}c_{i}}{\,\text{d}\epsilon_{\text{PSF},i}}\leavevmode\nobreak\
.$ (10)

We focus primarily on the $c$ dependence on $\epsilon_{\text{PSF}}$ as instead
the $m$ dependence is negligible as long as the PSF is stable and its
variation in size is within a percent levels. Earlier studies (Massey et al.,
2012; Cropper et al., 2013) have set out requirements for space-based missions
on $m$ and $c$ based on a top-down error breakdown from cosmology to two-point
statistics. For Euclid, the requirement is on the statistical error on bias,
$\sigma_{m}<2\times 10^{-3}$ and $\sigma_{c}<3\times 10^{-4}$. That is roughly
equivalent to saying that a shear of $1\%$ should be measured with a
fractional accuracy and precision of $0.2\%$. Note that the requirement is an
order of magnitude more stringent than current ground-based experiments
(Hildebrandt et al., 2017). The detailed breakdown of the total budget on $m$
and $c$ into various error terms (Cropper et al., 2013) suggests we can set
the required statistical error on the bias due to the measurement alone to
$\sigma_{m}<5\times 10^{-4}$ and $\sigma_{c}<5\times 10^{-5}$. Therefore, in
order to measure $|g|\approx 0.03$ with a residual post-calibration
multiplicative bias smaller than $\sigma_{m}$, one will need at least
$n=\sigma_{\epsilon}^{2}\,|g|^{-2}\sigma_{m}^{-2}\approx 4\times 10^{8}$
galaxies,444This assumes we measure shear with accuracy given by
$\hat{g}=(1+m)\,g$. The standard deviation of the measured shear scales as
$\sigma_{\epsilon}/\sqrt{n}$ and therefore we require that
$\sigma_{\epsilon}/\sqrt{n}\lesssim\sigma_{m}|g|$. where
$\sigma_{\epsilon}\approx 0.3$ is the ‘shape noise’, i.e., the standard
deviation of the per-component intrinsic ellipticity distribution. Obviously
measurement noise and intrinsic scatter in the morphological properties will
also need to be factored in. A ballpark estimate for the Euclid requirement on
PSF leakage that we will be using as benchmark in our analysis is
$\sigma_{\alpha}\lesssim\sigma_{c}/\left|\delta\epsilon_{\text{PSF}}\right|$,
where $|\delta\epsilon_{\text{PSF}}|\approx 0.1$ is the order of magnitude
(absolute) variation in PSF ellipticity across the field of view, which yields
$\sigma_{\alpha}<1\times 10^{-3}$ if we assume a budget of $\sigma_{c}<1\times
10^{-4}$. This derivation may be too conservative, as a full propagation of
the errors and biases through to cosmological parameters has been demonstrated
to be able to capture the spatial pattern imprinted by the PSF and other
effects (Euclid Collaboration: Paykari et al., 2020). Other surveys have
implemented other solutions such a first-order expansion on PSF ellipticity
and PSF model residuals in KiDS (Heymans et al., 2006; Giblin et al., 2021) or
the angular correlations between PSF ellipticity and size implemented in the
rho statistics in DES (Jarvis et al., 2020).

In the next subsections we address the key elements of the LensMC measurement
method: galaxy modelling, PSF convolution, likelihood, sampling, and a further
discussion about handling real data effects. We emphasise the role of joint
measurement of objects to address neighbour bias, which is a concern for
current and upcoming surveys, and also our MCMC strategy to sample a multi-
dimensional parameter space and marginalise each lensing target over all
nuisance parameters and other objects.

###  2.1 PSF-convolved galaxy models

We assume 2D-projected galaxy models as a mixture of two circular Sérsic
profiles (Sérsic, 1963), for disc,

$I_{\text{d}}(r)\propto\exp\left(-\frac{r}{r_{\text{e}}}\right)\leavevmode\nobreak\
,$ (11)

and bulge,

$I_{\text{b}}(r)\propto\exp\left[-a_{\text{b}}\left(\frac{r}{r_{\text{h}}}\right)^{\frac{1}{n_{\text{b}}}}\right]\leavevmode\nobreak\
,$ (12)

where $r$ is the distance from the centre, $r_{\text{e}}$ is the exponential
scale length of the disc, $r_{\text{h}}$ is the bulge half-light radius,
$n_{\text{b}}=1$, and $a_{\text{b}}\approx 2\,n_{\text{b}}-0.331$ (Peng et
al., 2002). The bulge Sérsic index is fixed to $1$ based on recent multi-
wavelength observations of the Hubble CANDELS/GOODS-South field (Welikala et
al., in prep.)555Their work highlights that both the dust in the disc and the
3D modelling of the galaxy influence the inferred bulge structural
parameters., while bulge profiles with $n_{\text{b}}=4$ (de Vaucouleurs) are
instead more typical for early-type galaxies at low redshift. Both profiles
are normalised so that their integral is 1. In the measurement, $r_{\text{e}}$
plays the role of object size parameter, and we fix the bulge-to-disc scale
length ratio to $r_{\text{h}}/r_{\text{e}}=0.15$ based on the same Hubble
Space Telescope measurements (Welikala et al., in prep.). Finally, we impose a
hard cut-off on the surface brightness profile at $r_{\max}/r_{\text{e}}=4.5$
since observations indicate that galaxies have truncated surface brightness
distributions (Van der Kruit & Searle, 1981, 1982). The parameters
$n_{\text{b}}$, $r_{\max}/r_{\text{e}}$, and $r_{\text{h}}/r_{\text{e}}$ are
assumptions made in the modelling that can potentially lead to large biases in
presence of a mismatch in the assumed Sérsic index compared to simulations
(Simon & Schneider, 2017). We stress that our choice of fixed values are based
on recent observations, and the model bias due to incorrect assumptions are
often intertwined with the particular simulation setup and its complexity. A
detailed investigation of sensitivity of the calibration to bulge parameters
is presented later in this paper.

The Sérsic model introduced above is an isotropic profile with zero
ellipticity. To make it anisotropic, i.e., elliptical with ellipticity
$\epsilon=\epsilon_{1}+\mathrm{i}\,\epsilon_{2}$, we use the following
distortion matrix

$\mathsf{S}=\frac{\bar{r}_{0}}{q_{\epsilon}\,r_{0}}\begin{pmatrix}1-\epsilon_{1}&-\epsilon_{2}\\\
-\epsilon_{2}&1+\epsilon_{1}\end{pmatrix}\leavevmode\nobreak\ ,$ (13)

where $r_{0}/\bar{r}_{0}$ is the scale factor necessary to scale a model of
size $\bar{r}_{0}$ to the desired size $r_{0}$. Because observed galaxy shapes
are a 2D visual projection of an intrinsically 3D distribution, we introduce
the additional scale factor, $q_{\epsilon}=1-|\epsilon|$, to make the profile
semi-major axis invariant under ellipticity transformation.666Instead,
$q_{\epsilon}=1$ would make the profile area invariant under an ellipticity
transformation. The choice of $q_{\epsilon}$ leads to different shear bias
properties that can have significant impact on the final calibration (Fenech
Conti et al., 2017). Discs and bulges typically show different intrinsic
ellipticity. As discs will be observed more elliptical if edge-on, their
ellipticity is primarily driven by inclination angle. In contrast, bulges are
spheroids that are almost invariant under inclination, so they will appear
more circular. In the measurement, we still apply the same ellipticity to both
components as part of our 2D modelling, but we are aware that a positive
ellipticity gradient from the intrinsically 3D distribution would induce a
bias if not fully captured (Bernstein, 2010). Our ellipticity estimate will
therefore be a proxy of the inclination angle, especially for disc-dominated
galaxies. Any residual ellipticity gradient, if significant, will have to be
addressed separately as part of the calibration.

The Euclid telescope, optical elements, and detector introduce distortions of
the input galaxy brightness distribution, which must be corrected. The effect
is mostly convolutive, which tends to blur the galaxy image further. An
example of a typical PSF image for a space-based telescope like Euclid is
given in Fig. 1. This is: (i) close to being diffraction limited; (ii)
undersampled due to the half width being comparable with the pixel size; (iii)
chromatic due to the integration over a wide range of wavelengths in the VIS
filter; (iv) SED dependent due to integration being weighted by a combination
of galaxy bulge and disc SEDs;777With galaxy bulges being, on average, redder
than discs. However, PSF images will be generated from the total galaxy SED.
(v) spatially variant across the field of view due to optical distortions at
the exit pupil; and (vi) epoch variant due to varying Solar aspect angle
throughout the mission inducing thermal distortion on the optics. A
comprehensive study of the modelling will be presented elsewhere (Duncan et
al., in prep.). A smaller contribution also comes from the CCD pixel response
function, which models the response of the detector pixel as an integrated
measure of the incoming flux illuminating individual pixels. In forward
modelling, we include all the convolutive effects due to the telescope PSF and
CCD, individually for bulge, disc, and for each image exposure. Colour
gradients originate from incorrectly using the total galaxy SED when
generating the PSF image, while the bulge and disc will have naturally
different SEDs (Semboloni et al., 2013; Er et al., 2018). Using individually
PSF-convolved model components may help to control colour gradients, if
individual SEDs were available. However, the impact of colour dependence and
gradients on our analysis is not addressed here since we assume the same
broadband PSF as obtained from an SED of an SBc-type galaxy at redshift 1 in
both simulations and measurements. Further non-linear distortions, such as in
the case of charge transfer inefficiency (CTI, Rhodes et al., 2010) and the
brighter-fatter effect (BFE, Antilogus et al., 2014), are typically corrected
for at the image pre-processing stage, but residuals could still affect the
shear measurement (Massey et al., 2014; Israel et al., 2015), which we do not
include here.

Figure 1: Example of a Euclid chromatic PSF for an assumed SED of a SBc-type
galaxy at redshift 1. The flux in the image is rescaled to its maximum value,
and oversampled by a factor 3 with respect to the native VIS pixel size of
$None$. Diffraction spikes are clearly visible at a significant distance from
the centre, despite the blurring due to the chromaticity. The full width at
half maximum (FWHM), including its variation across the field of view, is
$$$^{+0.0040}_{-0.0019}$, comparable with the Euclid pixel size, so images
will be undersampled at the Nyquist limit. The ellipticity is
$\epsilon_{1,\text{PSF}}=0.017^{+0.038}_{-0.024}$ and
$\epsilon_{2,\text{PSF}}=0.001^{+0.042}_{-0.020}$, with the superscript and
subscript denoting absolute ranges.

Galaxy modelling for large-volume surveys like Euclid requires fast and
efficient rendering of the images. All operations described so far are best
suited to work in the Fourier space. We adopt a similar approach to galsim
(Rowe et al., 2015). Consider the generic profile $I(r)$, which could be
either Eq. (11) or (12). Because of its isotropy, the 2D Fourier transform is
the 1D Hankel transform,

$\tilde{I}(k)=2\pi\int_{0}^{\infty}I(r)\,\text{J}_{0}(kr)\,r\,\text{d}r\leavevmode\nobreak\
,$ (14)

where $k$ is the Fourier frequency (sampled on an oversampled grid) and
$\text{J}_{0}$ is the Bessel function of the first kind. We call
$\tilde{I}(k)$ the template model which any profile with arbitrary choice of
parameter values (ellipticity, size, and position offset) can be derived
from.888For computational purposes, the transform is calculated once and saved
in a cached 1D array to minimise computing time. To render a galaxy profile
with parameters $\epsilon_{1}$, $\epsilon_{2}$, and $r_{\text{e}}$, we apply
the inverse distortion matrix $\mathsf{S}^{-1}$ to coordinates in Fourier
space, so anisotropic coordinates are now defined by
$\boldsymbol{k}^{\prime}=\mathsf{S}^{-1}\boldsymbol{k}$. A position shift by
$\boldsymbol{\delta r}=(\delta x,\delta y)$ from the centre999As it will be
explained later in the section, in reality the modelling actually includes
position shifts from right ascension, $\alpha$, and declination, $\delta$.
Therefore the position parameters will be $\Delta\alpha$ and $\Delta\delta$.
is equivalent to the phase $\boldsymbol{k}^{\prime}\cdot\boldsymbol{\delta
r}$. The sheared-stretched-shifted model becomes

$\tilde{I}^{\prime}(\boldsymbol{k^{\prime}})=\frac{1}{\det\mathsf{S}}\,\mathrm{e}^{-2\pi\mathrm{i}\,\boldsymbol{k}^{\prime}\cdot\,\boldsymbol{\delta
r}}\,\tilde{I}(k^{\prime})\leavevmode\nobreak\ ,$ (15)

where $\tilde{I}(k^{\prime})$ is the template calculated at the sheared-
stretched Fourier mode $\boldsymbol{k^{\prime}}$. Since isotropy is lost
through the operation above, $\tilde{I}^{\prime}(\boldsymbol{k^{\prime}})$ is
no longer a Hankel transform but a full Fourier transform, which is a function
of the vector $\boldsymbol{k^{\prime}}$. Note that if $\tilde{I}(k)$ is
precomputed and stored in a cached array, then this step will also involve
interpolating at the new coordinates. The interpolation error in the convolved
model can be seen at a large radius from the centre, but is around the same
level of the precision used to store the model itself. We alleviate the
problem of undersampling by calculating the PSF and galaxy model on
coordinates with a common oversampling factor of 3. It is worth noting that,
contrarily to the analytic approach adopted by galsim, our template models are
numerical arrays obtained by a Fourier transform of the isotropic model
arrays. The main reason for this approach is that the theoretical definition
of Hankel transform of Eq. (14) is invalidated by the finite limit of
integration and the oversampling of the model, which make our models a bit
more realistic. Finally, the oversampled models for PSF and galaxy are
multiplied together, the convolved model is downsampled by the same factor to
the actual pixel scale,101010Downsampling in real space corresponds to
aliasing in Fourier space, i.e., $n$-folding the transform and summing up. and
the downsampled convolved model is inverse fast Fourier transformed to real
space.

We further optimise the calculation of the models for speed by drawing on
square images of size proportional to the galaxy size being rendered, so
larger galaxies require larger arrays. Given a cut-off at
$r_{\max}=4.5\,r_{\text{e}}$ from the centre, a minimum image half-size of
$2\,r_{\max}$ is required to avoid aliasing from the Fourier transform.
Therefore the minimum image size will always be larger than
$18\,r_{\text{e}}$. As galaxies are expected to have mean size of $None$, but
spanning the whole range from just below resolution, $None$, to the largest
(although rare) galaxies, $None$, we define a template bank of pre-calculated
Hankel transforms of circular objects of different scale lengths
$\bar{r}_{0}=\\{0.2,0.4,0.8,1.6,3.2\\}\,\unit{\arcsec}$ and model arrays of
different sizes $\\{3.2,4.8,6.4,9.6,12.8,19.2,25.6,38.4\\}\,\unit{\arcsec}$.
To avoid aliasing from the interpolation, the template scale length
$\bar{r}_{0}$ is required to be slightly larger than the galaxy size being
rendered; we require $\bar{r}_{0}>\sqrt{2}\,r_{\text{e}}$. We apply the
operations of shear and stretch to this template bank to get any object of
arbitrary ellipticity and size before the convolution with the PSF takes place
as previously discussed.

The final galaxy model is a linear mixture of PSF-convolved co-centred
components. Let us label the profiles with a subscript ‘d’ for disc and ‘b’
bulge,

$I(r)=F_{\text{d}}\,I_{\text{d}}(r)+F_{\text{b}}\,I_{\text{b}}(r)\leavevmode\nobreak\
,$ (16)

where $F_{\text{d}}$ and $F_{\text{b}}$ are disc and bulge fluxes, and
$F=F_{\text{d}}+F_{\text{b}}$ is the total flux. If B/T is the bulge fraction,
then the fluxes are also defined by $F_{\text{d}}=F(1-\text{B/T})$ and
$F_{\text{b}}=F\;\text{B/T}$.

To summarise, given pre-computed template models for disc and bulge, we can
generate a galaxy model with a desired ellipticity, size, and position by
carrying out all the operations with simple algebra in Fourier space on
oversampled coordinates, and then take one Fourier transform each time at the
end. This is sufficiently fast for an intensive, repeated calculation of the
same model with varying realisations of galaxy parameters (ellipticity, size,
position offset, and fluxes). However, we keep $n_{\text{b}}$,
$r_{\max}/r_{\text{e}}$, and $r_{\text{h}}/r_{\text{e}}$ fixed in our
modelling as allowing too much freedom would induce strong degeneracies
between parameters and complicate the measurement substantially. We address
the model bias sensitivity in Sect. 4.4.

###  2.2 Likelihood

Suppose we have multi-exposure image data vectors
$\boldsymbol{D}=\\{\boldsymbol{D}_{\text{exp}}\\}$.111111On average 4
exposures for the Euclid Wide Survey and 64 for the Deep Survey fields. We
wish to estimate the model,
$\boldsymbol{I}=\\{\boldsymbol{I}_{\text{exp}}\\}$, that best represents the
available exposures. The model
$\boldsymbol{I}=\boldsymbol{I}(\epsilon,\theta,\phi)$ is a function of
ellipticity $\epsilon=(\epsilon_{1},\epsilon_{2})$, nuisance parameters
$\theta=(r_{\text{e}},\delta x,\delta y)$,121212As explained, though, we model
positions in world coordinates. and linear flux nuisance parameters
$\phi=(F_{\text{d}},F_{\text{b}})$.

Assuming Gaussian data,131313The Gaussian approximation holds true in the
limit of large counts in the image. the likelihood can be written as

$\ln
p(\boldsymbol{D}|\epsilon,\theta,\phi)\propto-\frac{1}{2}\Big{[}\boldsymbol{D}-\boldsymbol{I}(\epsilon,\theta,\phi)\Big{]}^{\top}\mathsf{C}^{-1}\Big{[}\boldsymbol{D}-\boldsymbol{I}(\epsilon,\theta,\phi)\Big{]}\leavevmode\nobreak\
,$ (17)

where $\mathsf{C}$ is the noise covariance matrix usually estimated from the
data as a block diagonal matrix, and the normalisation constant,
$1/2\ln[(2\pi)^{\lambda}\det\mathsf{C}]$ ($\lambda$ is the dimensionality),
has been ignored. The noise is intrinsically non-stationary since various
noise sources (such as the Poisson noise141414Poisson noise is a significant
noise source especially for bright objects. This term is included in the
simulations, but not in the measurement as it would require prior knowledge of
the distribution profile that is being measured. from the objects in the
image) vary spatially. Because the model is linear in the component fluxes,
$\boldsymbol{I}=F_{\text{d}}\,\boldsymbol{I_{\text{d}}}+F_{\text{b}}\,\boldsymbol{I_{\text{b}}}$,
it is straightforward to integrate over the fluxes,
$\phi=(F_{\text{d}},F_{\text{b}})$, and we have the following marginalised
likelihood,

$\ln
p(\boldsymbol{D}|\epsilon,\theta)\propto\frac{1}{2}\,\mathcal{F}^{-1}_{ij}(\epsilon,\theta)\,\rho_{i}(\epsilon,\theta)\,\rho_{j}(\epsilon,\theta)\leavevmode\nobreak\
,$ (18)

where $i$ indexes the model component,
$\rho_{i}(\epsilon,\theta)=\boldsymbol{D}^{\top}\mathsf{C}^{-1}\boldsymbol{I}_{i}(\epsilon,\theta)$
is a $2\times 1$ vector, $\boldsymbol{I}_{i}=\partial\boldsymbol{I}/\partial
F_{i}$ ($i=\text{d},\text{b})$, and $\mathcal{F}_{ij}$ is the $2\times 2$
Fisher matrix,

$\mathcal{F}_{ij}(\epsilon,\theta)=\boldsymbol{I}_{i}(\epsilon,\theta)^{\top}\mathsf{C}^{-1}\boldsymbol{I}_{j}(\epsilon,\theta)\leavevmode\nobreak\
.$ (19)

Note that because the right-hand side of Eq. (18) is quadratic in $\rho_{i}$
and $\mathcal{F}_{ij}$ is positive definite, we find $\ln
p(\boldsymbol{D}|\epsilon,\theta)>0$. A full derivation of the marginalised
likelihood, including edge cases and implementation, can be found in Appendix
A. The dimensionality of the problem has now been reduced from 7 free
parameters to 5: ellipticity, size, and position offsets.151515Having assumed
that the two components are co-centred and the bulge size is locked to the
disc size by a fixed rescaling.

Forward modelling provides solid grounds for a further generalisation to
measuring multiple objects jointly, especially if they are observed within a
short angular separation such as for neighbours. We label each likelihood with
the index $\omega$ running through the objects being jointly measured, $\ln
p_{\omega}(\boldsymbol{D}|\epsilon,\theta,\phi)$. The joint likelihood is then

$\ln p(\boldsymbol{D}|\\{\epsilon,\theta,\phi\\}_{\omega})=\sum_{\omega}\ln
p_{\omega}(\boldsymbol{D}|\epsilon,\theta,\phi)\leavevmode\nobreak\ ,$ (20)

where $\\{\epsilon,\theta,\phi\\}_{\omega}$ is the set of all parameters for
all the objects being measured. In the above equation we assume the
independence of the individual likelihoods. This is a fair assumption since
close neighbours will very often be so due to random visual alignment.
Consequently, those galaxies will be at different redshift and have different
shear. A much smaller fraction will include tidal interaction. In this case,
the galaxies will be at the same redshift and have the same shear. It it then
expected to have some degree of correlation between the individual
likelihoods. In more extreme but much rarer cases, the galaxies will be
tidally interacting and therefore our Sèrsic modelling would break down
entirely as we do not include any extra correlation term. Despite affecting a
very small fraction of objects, dedicated simulations would be required to
assess the impact on shear bias. Also, it is worth noting that we need to be
careful with the marginalisation of the individual likelihoods. The main issue
lies in the marginalised likelihood of Eq. (18). This relies on calculating
$\rho_{i}(\epsilon,\theta)$ for the various model components. However, when
multiple objects are present in the same neighbourhood, this quantity will
effectively introduce correlation between the likelihoods of the two objects.
Therefore the statistical independence required to multiply likelihoods
together will not be ensured. We have verified in testing that not
marginalising individual likelihoods is indeed the correct approach to the
problem. The joint likelihood is defined in a $7\times N$-dimensional
parameter space, where $N$ is the number of objects being measured jointly,
with $N=2$ being a typical number found in testing. For increased stability,
we first optimise the likelihood for fluxes and positions offsets, and then
also for ellipticity and sizes. This proves to be very robust as opposed to
iterating over individual objects after masking neighbours to achieve a
reliable initial guess (Drlica-Wagner et al., 2018). One key benefit of MCMC
is that it marginalises the ellipticity of an object over all remaining
nuisance parameters, which include object nuisance parameters as well as the
parameters of the other objects included in the joint sampling (see Sect.
2.3).

Our prior is based on enforcing hard bounds on all parameters:
$0\leq|\epsilon|<1$ given that the modulus of ellipticity in Eq. (2) cannot
exceed 1; $0\leq r_{\text{e}}\leq$$$, where the upper bound is based on
observations made in the Hubble Deep Fields; position offsets are restricted
to $\pm\,$$$ since the accuracy to which detections are made is typically sub-
pixel; and fluxes are positive. A more informative prior could be derived from
real observations in the future. A summary of all parameters, being free,
fixed or derived, is presented in Table 1.

Table 1: Summary of free, fixed and derived parameters in the modelling. The value should be interpreted as an initial guesses or constant depending on whether the parameter is allowed to vary or remain constant. | Free | Fixed | Derived | Initial value | Bounds | Unit |  Description
---|---|---|---|---|---|---|---
$\epsilon_{1}$ | ✓ |  |  | 0 | $|\epsilon|<1$ |  |  first component of ellipticity in the tangent plane to $(-\alpha,\delta)$ as defined in Eq. (13)
$\epsilon_{2}$ | ✓ |  |  | 0 | $|\epsilon|<1$ |  |  second component of ellipticity in the tangent plane to $(-\alpha,\delta)$ as defined in Eq. (13)
$r_{\text{e}}$ | ✓ |  |  | 0.3 | $[0,2]$ | arcsec |  effective radius: disc scale length as defined in Eq. (13)
$\Delta\alpha$ | ✓ |  |  | 0 | $[-0.3,0.3]$ | arcsec |  offset in right ascension as a phase prefactor in Eq. (15)
$\Delta\delta$ | ✓ |  |  | 0 | $[-0.3,0.3]$ | arcsec |  offset in declination as a phase prefactor in Eq. (15)
$\alpha$ |  |  | ✓ |  |  | degree |  right ascension
$\delta$ |  |  | ✓ |  |  | degree |  declination
$F_{\text{d}}$ | ✓ |  |  |  | $F_{\text{d}}\geq 0$ |  |  disc flux: marginalised over or free as defined in Eq. (16)
$F_{\text{b}}$ | ✓ |  |  |  | $F_{\text{b}}\geq 0$ |  |  bulge flux: marginalised over or free as defined in Eq. (16)
$F$ |  |  | ✓ |  | $F\geq 0$ |  |  total flux
B/T |  |  | ✓ |  | $[0,1]$ |  |  bulge fraction as defined in text after Eq. (16)
S/N |  |  | ✓ |  | $\text{S/N}\geq 0$ |  |  signal-to-noise ratio
$I_{\scriptscriptstyle\rm E}$ |  |  | ✓ |  |  |  |  magnitude: depends on the assumed zero-point, exposure gain and integration time
$n_{\text{b}}$ |  | ✓ |  | 1 |  |  |  Sérsic index of the bulge as defined in Eq. (12)
$a_{\text{b}}$ |  | ✓ |  | depends on $n_{\text{b}}$ |  |  |  Sérsic coefficient of the bulge as defined in Eq. (12), see Peng et al. (2002)
$r_{\text{h}}/r_{\text{e}}$ |  | ✓ |  | 0.15 |  |  |  bulge half-light radius to effective radius ratio as defined in text after Eq. (12), see Welikala et al., in prep.
$r_{\text{max}}/r_{\text{e}}$ |  | ✓ |  | 4.5 |  |  |  truncation radius to effective radius ratio as defined in text after Eq. (12)

###  2.3 Massive Markov Chain Monte Carlo sampling

Shear measurement poses serious difficulties in identifying the best strategy
to sample the posterior probability distribution of Eq. (4), assuming the
likelihood of Eq. (18) or (20),

  1. 1.

since the lensing sample is very broad in morphological properties, it will
contain both low and high S/N objects, whose posterior probability
distribution can be either very broad or very narrow; hence any sampling
strategy must be robust to this variability;

  2. 2.

the dimensionality of the problem is $7\times N$, where $N$ is the number of
objects being measured jointly to mitigate neighbour bias; sampling must then
be resilient to the large dimensionality, and provide marginalisation and
error estimations with minimum overheads;

  3. 3.

the shape of the distribution is a strong function of object parameters, such
as ellipticity and size, and therefore it varies significantly across the
sample; without prior knowledge of the physical properties of each galaxy, any
sampling method must run in a consistent, robust way;

  4. 4.

given the large sample size of order $10^{9}$ galaxies, sampling the posterior
is a computational challenge, so a trade-off between method complexity,
runtime, and access to computing resources must be identified;

  5. 5.

the sampling must be completely automated, without human supervision, and no
fine tuning of sampling parameters and options is allowed.

Considering all these challenges, our priority must be the average convergence
property of the sampler. The best strategy identified is MCMC, which allows us
to sample the posterior generated from the marginalised likelihood of Eq. (18)
for an individual object, or the joint likelihood of Eq. (20) for a group of
objects in an efficient and consistent way, for all $10^{9}$ objects in the
sample (hence the adjective ‘massive’). More importantly, MCMC seems to be the
best choice to tackle neighbours, particularly as an estimate can be found in
a high dimensional parameter space. It is worth noting that another key
benefit of MCMC sampling is that it is both a maximisation and sampling
procedure. The maximisation happens during the burn-in phase where the sampler
tries to reach the region of higher probability. The actual sampling happens
in the later stage of the chain after the burn-in phase. The marginalisation
over nuisance and error estimation are then natural by-products with no extra
overhead. This implies that not only can ellipticity estimates be marginalised
over object nuisance parameters, but also over other object parameters in the
joint group, hence minimising the impact of neighbours in the final
ellipticity estimate.

When searching for such an algorithm that could potentially suit our needs, we
have considered a number of potential candidates that are widely used in
cosmology and other fields (MacKay, 2002). The development of various sampling
methods has primarily been driven by the quest to achieve lower auto-
correlation and higher acceptance rate (Hastings, 1970; Swendsen & Wang, 1986;
Skilling, 2006; Goodman & Weare, 2010; Foreman-Mackey et al., 2013;
Betancourt, 2017; Karamanis & Beutler, 2021; Lemos et al., 2022). Although
appealing, all these methods do suffer from increased complexity, which is the
limiting factor in large-scale applications, where ‘large’ in this context
implies runs repeated over a billion times. Even for the most sophisticated
methods, it is often realised that a good initial guess is the key for good
sampling of the posterior.

For shear measurement on the scale of large galaxy surveys, there is not much
room for sampling complexity. The method has to be light enough and yet robust
to all the posteriors that need to be sampled. Furthermore, the likelihood
runtime limits the maximum number of samples that can be drawn for each galaxy
without having an overall impact on the survey analysis runtime. The
likelihood runtime is mostly dominated by the model component generation. The
measurement is dominated by sampling the posterior, plus some additional pre-
processing, therefore to limit the galaxy runtime to around few seconds per
galaxy, the MCMC sampling must be sparse. This is considered acceptable as we
are interested in accurate shear estimates, which are found by averaging over
ellipticity measurements. We choose an improved version of the Metropolis-
Hastings algorithm, which has been modified in two ways: (i) it incorporates
some of the ideas of parallel tempering (Swendsen & Wang, 1986; Sambridge,
2013), so the parameter space can be sampled on a larger scale initially, and
then on the smaller scale; (ii) it updates the proposal distribution during
the burn-in phase, automatically tuning it to find a good match with the
target distribution. Let $\vartheta$ be the generic vector of all parameters.
As explained in Sect. 2.2, this is $\vartheta=(\epsilon,\theta)$ for
individual galaxy measurement or $\vartheta=(\epsilon,\theta,\phi)$ when
jointly measuring groups of galaxies. Taking the logarithm of Bayes’ theorem
we have

$\ln p(\vartheta|\boldsymbol{D})\propto\ln p(\boldsymbol{D}|\vartheta)+\ln
p(\vartheta)\leavevmode\nobreak\ ,$ (21)

where the evidence has been ignored since our method is invariant to it. The
parameter vector at the current iteration step is denoted with
$\vartheta_{t}$, where $t=0,\ldots$ is the MCMC sample index. We also define a
tempering function, $T_{t}$, as a function of $t$. This acts as a Boltzmann
temperature and its expression will be defined later in the text. When the
temperature is high, $T_{t}\gg 1$, sampling from the posterior is equivalent
to sampling globally from the prior. When the temperature is gradually
reduced, as in annealing, $T_{t}\rightarrow 1$, we begin sampling directly
from the posterior. We define such a tempering function for application during
the burn-in phase only, and make sure $T_{t}=1$ for the final part of the
chain where we will take sample from. The method goes as follows:

  1. 1.

at $t$, draw a new sample $\vartheta^{\prime}_{t}$ from the proposal
distribution $q(\vartheta^{\prime}_{t}|\vartheta_{t})$; here we assume a
symmetric Gaussian proposal with mean $\vartheta_{t}$ and a pre-defined
diagonal covariance of $10^{-4}$ on all parameters (in units of arcsec for
size and position offsets);

  2. 2.

calculate the logarithm of the acceptance ratio

$\ln A=\frac{\ln p(\boldsymbol{D}|\vartheta^{\prime}_{t})-\ln
p(\boldsymbol{D}|\vartheta_{t})}{T_{t}}+\ln p(\vartheta^{\prime}_{t})-\ln
p(\vartheta_{t})\leavevmode\nobreak\ ;$ (22)

  3. 3.

accept or reject $\vartheta^{\prime}_{t}$ with probability $A$, i.e., draw $u$
from the uniform distribution on $[0,1]$ and accept $\vartheta^{\prime}_{t}$
if $u<\min(1,A)$; to speed things up and avoid calculating the likelihood
outside the prior, we immediately reject $\vartheta^{\prime}_{n}$ if
$p(\vartheta^{\prime}_{n})=0$.

For consistency, all posteriors are sampled from an initial guess that is a
circular galaxy of mean size, $r_{\text{e}}=$$$, and zero offset from the
nominal position in the sky. If we were to run any MCMC method from this point
onward we would end up with varying autocorrelation time depending on how far
the actual galaxy is from the initial guess, therefore we would need to wait
longer for very elliptical, small, or large galaxies. To improve the
convergence of the chains within a smaller number of iterations, we get a
better initial guess by running an initial maximisation of the posterior
before the actual MCMC. We run the conjugate-gradient search method (Powell,
1964) restricted to only 100 function evaluations, and then a downhill simplex
search (Nelder & Mead, 1965). The burn-in phase of the MCMC starts right
afterwards. During this phase the temperature is gradually lowered to 1. We
adopt the following power law cooling scheme (Cornish & Porter, 2007),

$T_{t}=\begin{cases}10^{\,f_{\text{heat}}(1-t/t_{\text{cool}})},&\text{if}\leavevmode\nobreak\
t<t_{\text{cool}}\\\ 1,&\text{otherwise}\end{cases}\leavevmode\nobreak\ ,$
(23)

where $f_{\text{heat}}=10$ is the heat factor and $t_{\text{cool}}=100$ is
cooling-down sample index. The parameter $t_{\text{cool}}$ represents the
number of samples it takes for the tempering function to become 1. Note that
$T_{0}=10^{\,f_{\text{heat}}}$ and $T_{t_{\text{cool}}}=1$. We begin with a
diagonal Gaussian proposal of width 0.01 on all parameters, which is then
recalculated from the chains every 100 samples and rescaled by the factor
$2.4\,{\lambda}^{-1/2}$, with $\lambda$ being the number of parameters
(Dunkley et al., 2005). The burn-in phase lasts for a total of 500 samples,
which is long enough for the tempering function to become 1, the proposal
covariance to be recalculated 5 times, and the chain to stabilise and reach
the high probability region (well before we start accumulating the final chain
samples). The final phase lasts for an additional $N_{\scriptscriptstyle\rm
MC}=200$ samples. Again, this is plenty to estimate both the mean and
covariance of the chains with enough precision, but we recognise that sampling
noise may still be non-negligible.

A good quantitative way to test the convergence of the chains is to
investigate their auto-correlation. We do so for a variety of galaxies and
results are shown in Appendix B. A less quantitative way is to verify that the
acceptance rate is within the expected range. We have also increased the final
200 samples up by a factor 5, without noticing any significant difference in
the shear results. For further verification, we have compared the method with
our implementation of affine invariant (Goodman & Weare, 2010) and parallel
tempering (Swendsen & Wang, 1986; Sambridge, 2013). While these methods
produce better ellipticity chains, they have not shown any significant
advantage in recovering shear, but increased complexity and therefore runtime
overhead.

Once the samples are drawn from the distribution function, calculating the
mean and width of the marginalised distribution becomes straightforward. Our
point estimate for ellipticity component marginalised over nuisance is the
mean of the chain after the burn-in phase,

$\hat{\epsilon}_{i}=\frac{1}{N_{\scriptscriptstyle\rm
MC}}\sum_{t}\epsilon_{i,t}\leavevmode\nobreak\ ,$ (24)

where $i=1,2$, and $\epsilon_{1,t}$ and $\epsilon_{2,t}$ are the two
ellipticity chains. The marginalised ellipticity covariance matrix is

$C_{ij}=\frac{1}{N_{\scriptscriptstyle\rm
MC}-1}\sum_{t}(\epsilon_{i,t}-\hat{\epsilon}_{i})\,(\epsilon_{j,t}-\hat{\epsilon}_{j})\leavevmode\nobreak\
.$ (25)

We calculate the averaged per-component variance (ignoring negligible
covariance between components),

$C_{\epsilon}=\frac{1}{2}(C_{11}+C_{22})\leavevmode\nobreak\ ,$ (26)

and choose to define the galaxy shear weight by

$w=\frac{1}{C_{\epsilon}+\sigma_{\epsilon}^{2}}\leavevmode\nobreak\ ,$ (27)

where $\sigma_{\epsilon}$ is the assumed shape noise, i.e., the width of the
1D intrinsic ellipticity distribution, ideally estimated in tomographic bins
from deeper measurements. We note the negligible sensitivity to the choice of
the $1/2$ factor in $C_{\epsilon}$.161616In fact, ignoring the $1/2$ factor
would lead to a redefinition of the weight,
$w=1/[(C_{11}+C_{22})/2+\sigma_{\epsilon}^{2}]\propto
1/[(C_{11}+C_{22})+\sigma_{\epsilon}^{\prime 2}]$ with
$\sigma_{\epsilon}^{\prime}=\sqrt{2}\,\sigma_{\epsilon}^{2}$, but results show
weak sensitivity to the value assumed for $\sigma_{\epsilon}$, as it will be
demonstrated at the end of Sect. 4.2. The MCMC provides a convenient and
efficient way to calculate both the mean and width of the ellipticity
posterior at no extra computational cost. The weights can then be used to
define sample averages, such as in 1-point estimates:

$\hat{g}_{i}=\frac{1}{\sum_{k}w_{k}}\sum_{k}w_{k}\,\hat{e}_{i,k}\leavevmode\nobreak\
,$ (28)

where $k$ indexes the galaxies in the lensing catalogue. The generalisation to
2-point estimates is straightforward. Please note that any choice of weight
leads to shear bias due to correlation with the measured ellipticity, and this
is tested later in the paper.

We have also implemented the self-calibration of ellipticity proposed by
Fenech Conti et al. (2017)171717The same correction can also be proved to map,
within some approximations, to other studies (Cox & Snell, 1968; Refregier et
al., 2012; Hall & Taylor, 2017) via importance sampling and likelihood ratio
while checking the quality of the sampling weights (Wraith et al., 2009),
without finding strong evidence of improvement. As results will be dominated
by other larger effects, we leave out further discussion from this paper.

###  2.4 Handling real data

Handling real data requires being careful with additional aspects of the
measurement. For instance, throughout our discussion we have proposed that our
sampling strategy is best suited to handle the presence of neighbours, i.e.,
resolved objects181818In this context ‘resolved’ implies that the object has
been detected and at least partially deblended so that our measurement can be
applied to all reported object positions. close to the lensing targets.
However, the situation is complicated by the fact that there is more variety
in real data as the brightness distribution of an object can be contaminated
in different ways depending on the nature of the nearby objects:

  1. 1.

neighbours (resolved galaxies or stars);

  2. 2.

blends (unresolved galaxies or stars);

  3. 3.

any other contamination (cosmic rays, transients, or ghosts).

Each case leads to a particular type of bias, and therefore we deal with close
objects in two ways. First, neighbours are grouped with a friend-of-friend
algorithm to a maximum angular separation of $r_{\text{friend}}=$$$. If the
separation is too small, the objects will be mostly measured in isolation,
therefore they will still be affected by the neighbours due to improper
masking. If the separation is too large, the groups will begin to be
unmanageable in size, but the benefit in controlling the neighbour bias will
be negligible. We find that $None$ is a good trade-off between measuring $N$
close neighbours jointly within a default postage stamp size of $None$, and
the non-negligible overhead in sampling a $7\times N$ dimensional posterior.
The joint analysis of object groups also gives us a good control of neighbour
bias, leading to a correction of $m\approx-7\times 10^{-4}$ as it will be
shown later in the paper.191919We do not attempt to optimise our choice of
$r_{\text{friend}}$ due to a number of other effects being more substantial
than this. Second, the segmentation maps and masks that are usually provided
with the data (Bertin et al., 2020; Kümmel et al., 2020) are combined in a
binary map and passed on to the likelihood to mask out objects in different
groups. Detector artefacts or cosmic rays are also masked out in a similar
way. In this case, to be even more conservative, we further dilate the masks
by one pixel so most of contamination bias should be taken care of. But
masking also helps partially with blends because objects that are false
negatives by the detection strategy may sometimes be true positives by the
masking procedure and therefore be masked out. Blending with faint galaxies
has been demonstrated to be relevant when trying to calibrate methods that are
particularly sensitive to the problem (Euclid Collaboration: Martinet et al.,
2019). We demonstrate that, to some extent, this is also the case in our
simulations where we measure objects deeper than the Euclid nominal survey
depth, as we will show in the next section.

Real images have a background level that needs to be subtracted. LensMC uses
the background estimates and noise maps that the Euclid processing provides,
but residual local background variations are subtracted at the individual
object group level. This is implemented by post-masking median subtraction.
Likewise, the standard deviation of the background noise is estimated after
masking. We measure galaxies and stars with the same model described earlier
in this section. We find that good star-galaxy separation is based on
selecting objects whose measured size is greater than the PSF size. This
method fits well with joint measurement of groups, however at the price of
rejecting faint galaxies that would anyway have negligible weight or be hard
to distinguish from faint stars. More details will be given in Sect. 4.

A measurement is made in sky coordinates using the supplied world coordinate
system (WCS) solution, which includes both linear and non-linear distortions
(Greisen & Calabretta, 2002). We assume the default coordinate system
$(-\alpha,\delta)$, where $\alpha$ is the right ascension and $\delta$ is the
declination. We measure position offsets from the provided nominal position in
arcsec. The resulting $\alpha$ and $\delta$ are reported in degrees, and
$r_{\text{e}}$ in arcsec. Likewise, the measured ellipticity is defined in the
tangent plane to the $(-\alpha,\delta)$ coordinate system centred at the
object position. We use the WCS to estimate a local linear approximation of
the mapping from sky coordinates to tangent plane coordinates at the nominal
position. We define 9 points in a square grid of size $None$ in pixel
coordinates centred at the nominal position, map them to sky coordinates, and
finally map the sky coordinates back to the undistorted tangent plane. The
Jacobian matrix, which models the local linear approximation of the mapping,
is the least square solution to the mapping from sky coordinates to tangent
coordinates. As part of this procedure, we also calculate the astrometric
offsets due to the exposures being dithered differently. The brightness model
is then correctly generated taking into account both the local distortion and
astrometric offsets so all the observables are measured uniquely in tangent
plane to sky coordinates.

When reporting our measurement we always compute $\chi^{2}=-2\ln p/\nu$, where
$p$ is the likelihood of Eq. (18) or (20) calculated at the mean estimate and
$\nu$ is the number of degrees of freedom. The $\chi^{2}$ will not in general
follow the theoretical distribution for a number of reasons. The noise is only
approximately Gaussian and non-Gaussianities will always be present in the
data. For instance, key examples are the Poisson noise from the background and
the object, digitalisation noise, non-linear artefacts, modelling mismatches,
or failures in the sampling. Nonetheless, the $\chi^{2}$ metric defined in
this way is still a good statistical measure of the quality of the
measurement. We also compare the $\chi^{2}$ calculated above with the same
quantity, which we call $\chi^{2}_{\text{bkg}}$, after having masked out all
the objects, which is expected to be just noise. Objects will be flagged up if
the $\chi^{2}$ is not consistent with the background. Following an F-test
procedure, we calculate the test statistic
$(\chi^{2}/\nu)\,(\chi^{2}_{\text{bkg}}/\nu_{bkg})^{-1}$ and reject the null
hypothesis (the measured $\chi^{2}$ is consistent with the background) if the
p-value is less than 0.01. Nonetheless we find that the impact of flagged
objects is negligible, so we usually include them in our results. However,
that may not be true for real data where the contamination from data artefacts
will be more important.

The measurement includes estimation of the object magnitude based on the
supplied zero-point, gain, and exposure time. Each exposure may come with its
own values as these varies both spatially and temporally, therefore it is
important to normalise the data to common units. As the data is measured in
analogue-to-digital units (ADU), we multiply each exposure by
$G\,10^{-(I_{\scriptscriptstyle\rm E,\,0}-\bar{I}_{\scriptscriptstyle\rm
E,\,0})/2.5}/\tau$, where $G$ is the gain in $\unit{e^{-}/ADU}$,
$I_{\scriptscriptstyle\rm E,\,0}$ is the magnitude zero-point,
$\bar{I}_{\scriptscriptstyle\rm E,\,0}$ is the average magnitude zero-point
across the exposures, $\tau$ is the exposure time, and the data is now in
normalised photoelectron count rate of $\unit{e^{-}/s}$. The flux, F, is then
measured in the same units, and we can estimate the magnitude as follows

$I_{\scriptscriptstyle\rm
E}=-2.5\operatorname{log_{10}}\frac{F}{\unit{e^{-}/s}}+\bar{I}_{\scriptscriptstyle\rm
E,\,0}\leavevmode\nobreak\ .$ (29)

The specific values for zero-point, gain and exposure time assumed in our
simulations will be provided in Sect. 3.

Analysing a volume of about 1.5 billion galaxies for Euclid will be a massive
computational challenge, especially if employing MCMC to sample the posterior.
Our measurement on highly realistic images runs at about 5 seconds per galaxy
per exposure per computing core, including joint measurement of
groups.202020The overhead of the joint measurement is about half of the quoted
total. We have discussed the benefits of using a fast, efficient
implementation of MCMC in the previous section. Here we want to highlight the
fact that all the pre- and post-processing described above add very little
overhead to the measurement. We find that the maximum posterior does suffer
from a large bias of $m\approx-1\times 10^{-2}$, which is about twice the bias
obtained when using the mean of the MCMC samples. Since the bias tends to
increase with magnitude, we interpret it as the maximum posterior estimate of
the ellipticity being more prone to noise bias. This is further evidence that
the MCMC can mitigate noise bias by consistent sampling and marginalisation of
a multi-dimensional posterior, in particular when jointly measuring groups of
objects, with the full complexity of real data and at the modest cost of
slightly more overhead.212121Compared to the maximum estimate, the MCMC adds
only $40\%$ to the total runtime.

##  3 Simulations

In order to validate our measurement method in a realistic setup, we design a
suite of simulations that incorporate most of the real data effects that
future lensing surveys like Euclid will need to account for. It is essential
then to bring in as much realism as possible. One problem that all shear
methods have to deal with is clustering that leads to close neighbours, which
is a concern for Euclid, Rubin, and present surveys as well. Because the
inferred bias depends on the details about the realism of clustering of faint
galaxies, this has to be incorporated in simulations particularly for
calibration purposes (Kannawadi et al., 2019; Euclid Collaboration: Martinet
et al., 2019). To make our custom simulations realistic and bring in all those
effects we are most concerned about, we utilise the exquisite, state-of-the-
art Flagship simulation mock galaxy catalogue (Potter et al., 2017, Euclid
Collaboration: Castander et al, in prep.),222222We ingest the catalogue
version 2.1.10 retrieved from the official website cosmohub.pic.es (Carretero
et al., 2017; Tallada et al., 2020). developed specifically for Euclid. The
same Flagship simulation is also used for the Euclid Science Ground Segment
simulations (Euclid Collaboration: Serrano et al., in prep.). Flagship
provides, in particular, a realistic distribution of galaxy morphologies, and
clustering of galaxies obtained through a full $N$-body dark matter
simulation.

The morphological parameters and spatial distributions are provided over an
octant of the full sky, which is just less than 40% of the Euclid Wide Survey.
Here we use values for the provided disk ellipticity and orientation angle,
disc scale length, VIS flux, bulge fraction, and position over a region
defined by $150^{\circ}<\alpha<230^{\circ}$ and
$15^{\circ}<\delta<85^{\circ}$. We also select all galaxies that are
classified in the catalogue as being either central or satellite in the halo,
and exclude quasars or high-redshift galaxies. Figure 2 shows the joint size-
magnitude distribution of galaxies. A significant fraction of the galaxies
have intrinsic effective radii similar to the PSF, which has a FWHM of
$$$^{+0.0040}_{-0.0019}$, and therefore appear only marginally resolved in the
PSF-convolved images. It is worth highlighting that because of the very faint
magnitude limit ($I_{\scriptscriptstyle\rm E}<29.5$, but complete to
$I_{\scriptscriptstyle\rm E}<27$) a significant fraction of the objects will
be too faint to be detected, but these will be still included in the
background noise. In addition to galaxies provided by Flagship, we also
simulate a uniform spatial distribution of stars up to
$I_{\scriptscriptstyle\rm E}<26$. Figure 3 shows the number count of galaxies
and stars. The galaxy count is obtained from Flagship and compared against a
polynomial model to VIS-corrected magnitudes in the GOODS-South field up to 26
(Giavalisco et al., 2004) and Ultra Deep Field beyond 26 (Beckwith et al.,
2006). Stars are drawn from a polynomial model of $i$ magnitudes generated
with the Besançon model (Czekaj et al., 2014) in an area of
$\qty{10}{deg^{2}}$ around the North Ecliptic Pole. Overall, we obtain a
number density of $\qty{250}{arcmin^{-2}}$ ($I_{\scriptscriptstyle\rm
E}<29.5$) for galaxies and $\qty{6}{arcmin^{-2}}$ ($I_{\scriptscriptstyle\rm
E}<26$) for stars.

Figure 2: Input magnitude-size distribution of galaxies. The data points are
the mean $r_{\text{e}}$ as a function of $I_{\scriptscriptstyle\rm E}$. Also
shown are the standard deviation of $r_{\text{e}}$ and
$I_{\scriptscriptstyle\rm E}$ in each bin. The horizontal band denotes the PSF
FWHM and its variation across the field of view. A significant fraction of the
galaxies have intrinsic effective radii similar to the PSF, especially at the
Euclid wide field (WF) depth. Figure 3: Input differential number count for
galaxies and stars. Solid is the measured galaxy count from Flagship. Dashed
is a polynomial model of VIS-corrected magnitudes in the GOODS South
($I_{\scriptscriptstyle\rm E}<26$) and Ultra Deep Field
($I_{\scriptscriptstyle\rm E}>26$). Stars are drawn from a polynomial model of
$i$ magnitudes generated with the Besançon model in the North Ecliptic Pole.
The cumulative number counts are $\qty{250}{arcmin^{-2}}$
($I_{\scriptscriptstyle\rm E}<29.5$) for galaxies and $\qty{6}{arcmin^{-2}}$
($I_{\scriptscriptstyle\rm E}<26$) for stars.

We define sky patches of size $None$ (about $\qty{40}{arcmin^{2}}$), broadly
corresponding to the size of a single Euclid CCD, but also include an adjacent
area of extra 10% (called buffer/guard region) all around the patch to draw
objects in the image that will not be part of the measurement. When selecting
the morphological properties from the Flagship catalogue (ellipticity, disc
scale length, bulge fraction, position, and flux), their AB flux is first
converted to AB magnitude, then further converted to VIS photoelectron count
rate via the magnitude-flux relation of Eq. (29)). We assume a constant
magnitude zero-point of $I_{\scriptscriptstyle\rm E,\,0}=25.719$, which has
been calculated using Euclid as-designed system throughput data. Star
positions are drawn from the count model uniformly in each patch.

All galaxies in each patch have ellipticity assigned by Flagship. In principle
we could use the cosmic shear from Flagship, however in this work we apply the
same constant shear to all galaxies in a patch, with the shear varying from
one patch to another according to a uniform distribution on a circle of radius
$|g|=0.02$ and random orientation. This choice mimics the typical shear
expected for a real survey and also minimises the error on multiplicative
bias. We assume an (infinitely thin) annular distribution as opposed to a disc
distribution or an even more realistic log-normal distribution because we want
to minimise the statistical error on multiplicative bias, and obviously be as
cosmology agnostic as possible. On the other hand, a variable shear field
might in principle introduce an additional correlation with neighbour bias,
particularly if neighbours at different redshifts are considered (MacCrann et
al., 2021). However, capturing realistic clustering is the most important
aspect of the simulations, which is what we focus on in this work. Similarly
to previous work (Bridle et al., 2010; Kannawadi et al., 2019), we apply shape
noise mitigation by making, in total, 4 clones of each patch with all
ellipticities rotated by $45^{\circ}$ while maintaining the same overall
constant shear, which gives us significant leverage on the required simulation
volume. It is worth noting, though, that a varying shear could also be dealt
with in a shear response approach, leading to an increased effective sample
size and reduce simulation volume in calibration (Pujol et al., 2019; Jansen
et al., 2024).

We set up a suite of simulations for each of 9 realisations of the PSF image
drawn at different positions in the field of view. While varying the PSF
image, we keep the objects at the same positions. We assume a Euclid PSF model
for a fixed SBc-type galaxy SED at $z=1$, the median of the distribution. The
mean ellipticity and its variation across the field of view is:
$\epsilon_{1,\text{PSF}}=0.017^{+0.038}_{-0.024}$ and
$\epsilon_{2,\text{PSF}}=0.001^{+0.042}_{-0.020}$, with the superscript and
subscript denoting absolute ranges. We note that this variation, if not
included in the modelling, would be responsible for an error in the shear
measurement that would far exceed science requirements. We do not include PSF
mismodelling in our simulations as the current Euclid requirement on PSF
ellipticity error is already quite stringent, but will be addressed elsewhere.
The Euclid Wide Survey is designed to take 4 dithered exposures (pointings),
plus two extra short exposures, of the same sky area. Most often these will be
taken in the same observation. Hence the PSF model is not expected to vary too
much across the exposures, but the images will be different because taken at
different positions in the field of view.

We generate Euclid detector images containing galaxies rendered with the
brightness model of Sect. 2.1 with varying ellipticity, $r_{\text{e}}$,
position, and fluxes. In our initial tests, we make our results insensitive to
model bias by construction, and therefore we fix $n_{\text{b}}$,
$r_{\text{h}}/r_{\text{e}}$ and $r_{\max}/r_{\text{e}}$. Later on, we address
model bias sensitivity by allowing $n_{\text{b}}$, and
$r_{\text{h}}/r_{\text{e}}$ to vary. For stars, we use a restricted model with
zero ellipticity and $r_{\text{e}}$, so we effectively render point-like PSF
images. For the measurement, we will be using the same galaxy profile (with
fixed bulge parameters) for all detected objects.

The pixel photoelectron noise is given by

$\sigma_{\text{px}}^{2}(x,y)=\left(R_{\text{bkg}}+R_{\text{dark}}\right)\tau+\lambda_{\text{obj}}(x,y)+\sigma_{\text{read}}^{2}\leavevmode\nobreak\
.$ (30)

The first term is Poisson noise from a constant zodiacal light background and
dark current, with rates $R_{\text{bkg}}=\qty{0.225}{e^{-}/s}$ and
$R_{\text{dark}}=\qty{0.001}{e^{-}/s}$, and exposure time $\tau=\qty{565}{s}$.
The second term, $\lambda_{\text{obj}}(x,y)$, is spatially varying Poisson
noise from all the objects in the image, which is non-negligible in the Euclid
VIS images. The third term is Gaussian noise from the CCD readout with a
constant $\sigma_{\text{read}}=\qty{4.5}{e^{-}}$. We assume that all noise
sources are uncorrelated.232323For instance, the readout noise could
potentially be correlated with non-negligible pixel covariance, which we
ignore in this work. In generating the images, we also apply a bias of
$\qty{2000}{e^{-}}$ (about as expected for Euclid), and finally digitise the
data. Digitalisation corresponds to dividing the image by a gain of
$\qty{3.1}{e^{-}/ADU}$242424The in-orbit detectors will have slightly larger
gain, likely around $\qty{3.4}{e^{-}/ADU}$, but this is not expected to change
any of our results. and floor truncating it to nearest integer, which itself
adds uniform noise of variance $1/12\,\unit{ADU}$. We set a tangent projection
as our WCS at the centre of the patch, draw 4 undithered exposures and stack
them up by taking their average. These images will be used by the object
detection for the main results presented here, but we will also include a
discussion about the dithering. An example of stacked CCD image is shown in
Fig. 4.

Figure 4: Example of LensMC-Flagship image. The input galaxy distribution is
provided by Flagship and stars are drawn from a model. We emulate the VIS
detector by including realistic image properties and noise, but we do not
include non-linearities, CTI, BFE, or cosmic rays. To aid the visualisation of
the faint objects, the image has been clipped between the 10th and 90th
percentiles, illustrating the sheer number of objects and their clustering.
The image size is 4096 pixels.

##  4 Results

To quantify the performance of LensMC on our realistic LensMC-Flagship
simulations, we run the measurement on about $45\,000$ random patches, which
is equivalent to an area of about $\qty{500}{deg^{2}}$, with mean number
density, according to Fig. 3, of $\qty{250}{arcmin^{-2}}$
($I_{\scriptscriptstyle\rm E}<29.5$) for galaxies and $\qty{6}{arcmin^{-2}}$
($I_{\scriptscriptstyle\rm E}<26$) for stars. We measure the same area (with
the objects at the same positions) 9 times for varying noise realisations and
PSF across the field of view, totalling an equivalent, effective area of
$4500\,\unit{deg^{2}}$. We run all our simulations across the GridPP UK
network (Faulkner et al., 2005; Britton et al., 2009).252525Testing has taken
6 months, with our final run averaging 15 000 cores/day for two weeks. A
qualitative test of the measurement performance is shown in Fig. 5. After the
galaxy models have been subtracted from the image data, the residuals look
consistent with noise, for galaxies measured individually or jointly in
groups, despite the presence of neighbours. More quantitative tests will be
discussed as part of the validation presented in Appendix C.

| Data | Model(s) | Residuals
---|---|---|---
One galaxy |  |  |
One galaxy |  |  |
Two galaxies |  |  |
Two galaxies |  |  |
Three galaxies |  |  |
Figure 5: Examples of measurement performance. The target galaxies are denoted
with a cross. The image residuals look consistent with noise, for galaxies
measured individually or jointly in groups, despite the presence of
neighbours. All images have a size of 128 pixels.

###  4.1 Selection

For our baseline test, we run SourceXtractor++ (Bertin et al., 2020; Kümmel et
al., 2020)262626Version 0.19.2 with default settings as used in Euclid. We do
not test the sensitivity to changes in these settings and that will be the
focus of future work. to detect galaxies and stars in each of the undithered
stacked images. The code attempts to deblend detected objects and produces a
detection catalogue with a total number density of $\qty{88}{arcmin^{-2}}$
($I_{\scriptscriptstyle\rm E}<26.5$), and $\qty{34}{arcmin^{-2}}$
($I_{\scriptscriptstyle\rm E}<24.5$). Figure 6 contains the galaxy magnitude
selection applied by SourceXtractor++. This shows the number count of the
objects in the simulation and after the detection by SourceXtractor++. The
detection catalogue is complete to the magnitudes of interest, apart from
false positives of about $\qty{6}{arcmin^{-2}}$ ($I_{\scriptscriptstyle\rm
E}<24.5$), probably due to a combination of sub-optimal detection and
mismatching with the true input catalogue in presence of neighbours at those
magnitudes.

Figure 6: Differential number count for all objects in the simulation and
after the detection by SourceXtractor++. The detection catalogue is mostly
complete to $I_{\scriptscriptstyle\rm E}<24.5$, apart from false positives of
about $\qty{6}{arcmin^{-2}}$ ($I_{\scriptscriptstyle\rm E}<24.5$). The
distribution starts rolling off at 26, so a large fraction of faint objects is
not detected.

The detections are grouped (with $r_{\text{friend}}=$$$) according to their
reported SourceXtractor++ positions. LensMC goes through each object group and
measures the object parameters starting from an initial guess at the provided
SourceXtractor++ positions. If the size of the group is 1, LensMC will measure
the object in isolation and masks out neighbours through the supplied
segmentation maps. Instead, if it is greater than 1, LensMC will measure the
objects jointly, while masking out neighbours belonging to other groups. We
match the input catalogue with the measurement catalogue and within a maximum
angular distance of $None$ from the measured object (which also corresponds to
the LensMC maximum search region around the detected object). The few measured
objects that do not get a useful match to within that distance are then
flagged up and removed from the analysis. We test the sensitivity to the
maximum match distance without noticing any appreciable change to the bias.

A key selection applied to the measurement catalogue is the star-galaxy
separation. As found in applications to real data, the object size is an
excellent statistic to discriminate between galaxies and stars (Sevilla-Noarbe
et al., 2018). Therefore we classify objects according to measured
$r_{\text{e}}>r_{\text{s/g}}$ where $r_{\text{s/g}}=$$$, which is slightly
larger than the pixel size and image resolution. Note that we apply our star-
galaxy separation to broadband data simulated with a fixed choice of SED
representative of a typical galaxy at redshift 1. However, this does not test
how well the star-galaxy separation works with a broad range of galaxy SEDs,
and also with a clear distinction between galaxy and star SEDs. We quantify
the performance of our separation by calculating: i) $N_{\text{g}}$, the
number of true positives, i.e., galaxies correctly identified as such; ii)
$N_{\text{s}}$, the number of true negatives, i.e., stars correctly identified
as such; iii) $N_{\neg{\text{g}}}$, the number of false positives, i.e., stars
wrongly identified as galaxies; iv) $N_{\neg{\text{s}}}$, the number of false
negatives, i.e., galaxies wrongly identified as stars. The above numbers are
always defined in the measurement catalogue. The true positive rate (TPR) and
false positive rate (FPR) are

TPR
$\displaystyle=\frac{N_{\text{g}}}{N_{\text{g}}+N_{\neg{\text{s}}}}\leavevmode\nobreak\
,$ (31) FPR
$\displaystyle=\frac{N_{\neg{\text{g}}}}{N_{\neg{\text{g}}}+N_{\text{s}}}\leavevmode\nobreak\
.$ (32)

Realistic values of $\text{FPR}>0$ and $\text{TPR}<1$ are always linked to
type I and II errors in the shear analysis. Type I is the inclusion of stars
in the lensing sample, hence leading to potentially large negative
multiplicative bias. Type II is the omission of galaxies (with potentially
large shear signal) from the lensing sample which introduces selection bias
and also a dilation in statistical error.

For the sample of detected objects to the detection limit
($I_{\scriptscriptstyle\rm E}<26.5$) we find $\text{TPR}=93.3\%$,
$\text{FPR}=4.6\%$, purity of 99.8%,272727Astronomical completeness coincides
with TPR, but $\text{purity}=N_{\text{g}}/(N_{\text{g}}+N_{\neg{\text{g}}})$.
and a star fraction of 6.6%.282828$\text{Star
fraction}=N_{\text{s}}/(N_{\text{g}}+N_{\text{s}})$. The TPR gives us the
frequentist probability of a positive being a galaxy, so
$\text{TPR}=p(+|\text{g})$. Similarly, $\text{FPR}=p(+|\text{s})$. Bayesian
posterior probabilities provide a more meaningful interpretation of those
numbers. The prior probability of an object being a galaxy is $p(\text{g})$
and a star is $p(\text{s})=1-p(\text{g})$ (i.e., the star fraction). Applying
Bayes’ theorem we get the probability of a galaxy given a positive detection,

$p(\text{g}|+)=\frac{p(+|\text{g})\,p(\text{g})}{p(+|\text{g})\,p(\text{g})+p(+|\text{s})\,p(\text{s})}\leavevmode\nobreak\
,$ (33)

and similarly for $p(\text{s}|+)$. With the numbers above we obtain
$p(\text{g}|+)=99.7\%$ and $p(\text{s}|+)=0.3\%$ for all objects in the
detection catalogue. A more relaxed FPR of about $20\%$ would still give us
$p(\text{g}|+)=99\%$ and $p(\text{s}|+)=1\%$, given the strong imbalance
between the galaxy and star samples. These numbers give us reassurance that
once an object has been classified as a galaxy, there is an average
$3\,\sigma$ confidence that it will indeed be a galaxy for the entire sample
up to the detection limit ($I_{\scriptscriptstyle\rm E}<26.5$).

Figure 7 shows the size distribution of true galaxies and true stars and the
operating curve (true positive rate vs false positive rate for varying
threshold) of our classification with either a horizontal line or a cross to
denote the default threshold, $r_{\text{s/g}}=$$$. Both plots provide solid
justification for our choice of $r_{\text{s/g}}$, but confusion is evident
around $I_{\scriptscriptstyle\rm E}=24$, which might explain most of the false
positives. The area under the operating curve at the bottom of Fig. 7 is
large, and the curve itself is reasonably flat for a wide range of false
positive rate suggesting excellent discrimination and weak sensitivity on the
threshold (in that the shear bias does not appreciably change for a wide range
of threshold values around the nominal one). However, in real measurements, an
optimal value could be inferred from external data or simulations, hence
allowing for a dramatic reduction of the false positives at the expense of a
modest reduction in the true positives. Our TPR, FPR, and operating curve of
Fig. 7 are consistent or better than the best estimators presented in Sevilla-
Noarbe et al. (2018), although a key caveat in our work is likely to be that
we have not included a full colour variation of galaxy and star SEDs, which
would lead to variable PSFs and potentially harder separation. Moreover, we do
not investigate any effect due to star density variation, which might well
change by a factor 2 or 3 going from the high to the low latitudes.

Figure 7: Star-galaxy separation of detected objects. (Top) Observed size-
magnitude distribution of true galaxies and stars. The data points are the