Paper deep dive
Neural posterior estimation of the neutrino direction in IceCube using transformer-encoded normalizing flows on the sphere
R. Abbasi, M. Ackermann, J. Adams, J. A. Aguilar, M. Ahlers, J. M. Alameddine, S. Ali, N. M. Amin, K. Andeen, C. Argüelles, Y. Ashida, S. Athanasiadou, S. N. Axani, R. Babu, X. Bai, A. Balagopal V., S. W. Barwick, V. Basu, R. Bay, J. J. Beatty, J. Becker Tjus, P. Behrens, J. Beise, C. Bellenghi, S. Benkel, S. BenZvi, D. Berley, E. Bernardini, D. Z. Besson, E. Blaufuss, L. Bloom, S. Blot, F. Bontempo, J. Y. Book Motzkin, C. Boscolo Meneguolo, S. Böser, O. Botner, J. Böttcher, J. Braun, B. Brinson, Z. Brisson-Tsavoussis, R. T. Burley, D. Butterfield, K. Carloni, J. Carpio, N. Chau, Z. Chen, D. Chirkin, S. Choi, A. Chubarov, B. A. Clark, G. H. Collin, D. A. Coloma Borja, A. Connolly, J. M. Conrad, D. F. Cowen, C. De Clercq, J. J. DeLaunay, D. Delgado, T. Delmeulle, S. Deng, P. Desiati, K. D. de Vries, G. de Wasseige, T. DeYoung, J. C. Díaz-Vélez, S. DiKerby, T. Ding, M. Dittmer, A. Domi, L. Draper, L. Dueser, D. Durnford, K. Dutta, M. A. DuVernois, T. Ehrhardt, L. Eidenschink, A. Eimer, C. Eldridge, P. Eller, E. Ellinger, D. Elsässer, R. Engel, H. Erpenbeck, W. Esmail, S. Eulig, J. Evans, P. A. Evenson, K. L. Fan, K. Fang, K. Farrag, A. R. Fazely, A. Fedynitch, N. Feigl, C. Finley, D. Fox, A. Franckowiak, S. Fukami, P. Fürst, J. Gallagher, E. Ganster, A. Garcia, M. Garcia, E. Genton, L. Gerhardt, A. Ghadimi, C. Glaser, T. Glüsenkamp, J. G. Gonzalez, S. Goswami, A. Granados, D. Grant, S. J. Gray, S. Griffin, K. M. Groth, D. Guevel, C. Günther, P. Gutjahr, C. Ha, A. Hallgren, L. Halve, F. Halzen, L. Hamacher, M. Handt, K. Hanson, J. Hardin, A. A. Harnisch, P. Hatch, A. Haungs, J. Häußler, K. Helbing, J. Hellrung, B. Henke, L. Hennig, F. Henningsen, L. Heuermann, R. Hewett, N. Heyer, S. Hickford, A. Hidvegi, C. Hill, G. C. Hill, R. Hmaid, K. D. Hoffman, A. Hollnagel, D. Hooper, S. Hori, K. Hoshina, M. Hostert, W. Hou, M. Hrywniak, T. Huber, K. Hultqvist, K. Hymon, A. Ishihara, W. Iwakiri, M. Jacquart, S. Jain, O. Janik, M. Jansson, M. Jin, N. Kamp, D. Kang, W. Kang, A. Kappes, L. Kardum, T. Karg, A. Karle, A. Katil, M. Kauer, J. L. Kelley, M. Khanal, A. Khatee Zathul, A. Kheirandish, T. Kim, H. Kimku, F. Kirchner, J. Kiryluk, C. Klein, S. R. Klein, Y. Kobayashi, S. Koch, A. Kochocki, R. Koirala, H. Kolanoski, T. Kontrimas, L. Köpke, C. Kopper, D. J. Koskinen, P. Koundal, M. Kowalski, T. Kozynets, A. Kravka, N. Krieger, T. Krishnan, K. Kruiswijk, E. Krupczak, A. Kumar, E. Kun, N. Kurahashi, C. Lagunas Gualda, L. Lallement Arnaud, M. J. Larson, F. Lauber, J. P. Lazar, K. Leonard DeHolton, A. Leszczyńska, C. Li, J. Liao, C. Lin, Q. R. Liu, Y. T. Liu, M. Liubarska, C. Love, L. Lu, F. Lucarelli, W. Luszczak, Y. Lyu, M. Macdonald, E. Magnus, Y. Makino, E. Manao, S. Mancina, A. Mand, I. C. Mariş, S. Marka, Z. Marka, L. Marten, I. Martinez-Soler, R. Maruyama, J. Mauro, F. Mayhew, F. McNally, K. Meagher, A. Medina, M. Meier, Y. Merckx, L. Merten, J. Mitchell, L. Molchany, S. Mondal, T. Montaruli, R. W. Moore, Y. Morii, A. Mosbrugger, D. Mousadi, E. Moyaux, T. Mukherjee, M. Nakos, U. Naumann, J. Necker, L. Neste, M. Neumann, H. Niederhausen, M. U. Nisa, K. Noda, A. Noell, A. Novikov, A. Obertacke, V. O'Dell, A. Olivas, R. Orsoe, J. Osborn, E. O'Sullivan, B. Owens, V. Palusova, H. Pandya, A. Parenti, N. Park, V. Parrish, E. N. Paudel, L. Paul, C. Pérez de los Heros, T. Pernice, T. C. Petersen, J. Peterson, S. Pick, M. Plum, A. Pontén, V. Poojyam, B. Pries, R. Procter-Murphy, G. T. Przybylski, L. Pyras, C. Raab, J. Rack-Helleis, N. Rad, M. Ravn, K. Rawlins, Z. Rechav, A. Rehman, I. Reistroffer, E. Resconi, S. Reusch, C. D. Rho, W. Rhode, L. Ricca, B. Riedel, A. Rifaie, E. J. Roberts, S. Rodan, M. Rongen, A. Rosted, C. Rott, T. Ruhe, L. Ruohan, D. Ryckbosch, J. Saffer, D. Salazar-Gallegos, P. Sampathkumar, A. Sandrock, G. Sanger-Johnson, M. Santander, S. Sarkar, M. Scarnera, M. Schaufel, H. Schieler, S. Schindler, L. Schlickmann, B. Schlüter, F. Schlüter, N. Schmeisser, T. Schmidt, A. Scholz, F. G. Schröder, S. Schwirn, S. Sclafani, D. Seckel, L. Seen, M. Seikh, S. Seunarine, P. A. Sevle Myhr, R. Shah, S. Shah, S. Shefali, N. Shimizu, B. Skrzypek, R. Snihur, J. Soedingrekso, D. Soldin, P. Soldin, G. Sommani, C. Spannfellner, G. M. Spiczak, C. Spiering, J. Stachurska, M. Stamatikos, T. Stanev, T. Stezelberger, T. Stürwald, T. Stuttard, G. W. Sullivan, I. Taboada, S. Ter-Antonyan, A. Terliuk, A. Thakuri, M. Thiesmeyer, W. G. Thompson, J. Thwaites, S. Tilav, K. Tollefson, J. A. Torres, S. Toscano, D. Tosi, K. Upshaw, A. Vaidyanathan, N. Valtonen-Mattila, J. Valverde, J. Vandenbroucke, T. Van Eeden, N. van Eijndhoven, L. Van Rootselaar, J. van Santen, J. Vara, F. Varsi, M. Venugopal, M. Vereecken, S. Vergara Carrasco, S. Verpoest, D. Veske, A. Vijai, J. Villarreal, C. Walck, A. Wang, E. H. S. Warrick, C. Weaver, P. Weigel, A. Weindl, J. Weldert, A. Y. Wen, C. Wendt, J. Werthebach, M. Weyrauch, N. Whitehorn, C. H. Wiebusch, D. R. Williams, L. Witthaus, G. Wrede, X. W. Xu, J. P. Yanez, Y. Yao, E. Yildizci, S. Yoshida, R. Young, F. Yu, S. Yu, T. Yuan, S. Yun-Cárcamo, A. Zander Jurowitzki, A. Zegarelli, S. Zhang, Z. Zhang, P. Zhelnin, P. Zilberman
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 96%
Last extracted: 4/26/2026, 4:41:04 PM
Summary
The paper presents a novel method for neutrino direction reconstruction in the IceCube detector using Amortized Neural Posterior Estimation (NPE). By utilizing a transformer encoder that maps to a specialized spherical normalizing flow on the 2-sphere, the method achieves state-of-the-art angular resolution for both 'tracks' and 'showers' event morphologies. The approach significantly outperforms traditional B-spline-based likelihood reconstructions in both speed and precision, enabling all-sky scans in seconds. The architecture incorporates C2-smooth rational-quadratic splines, scale transformations, and rotations to handle non-Gaussian uncertainties and complex event signatures.
Entities (7)
Relation Signals (5)
Normalizing Flow → definedon → 2-sphere
confidence 100% · maps to a normalizing flow on the 2-sphere.
IceCube → detects → Track
confidence 100% · the two main event morphologies in IceCube - tracks and showers
Abstract
Abstract:IceCube is a cubic-kilometer-scale neutrino detector located at the geographic South Pole. A precise directional reconstruction of IceCube neutrinos is vital for associations with astronomical objects. In this context, we discuss neural posterior estimation of the neutrino direction via a transformer encoder that maps to a normalizing flow on the 2-sphere. It achieves a new state-of-the-art angular resolution for the two main event morphologies in IceCube - tracks and showers - while being significantly faster than traditional B-spline-based likelihood reconstructions. All-sky scans can be performed within seconds rather than hours, and take constant computation time, regardless of whether the posterior extent is arc-minutes or spans the whole sky. We utilize a combination of $C^2$-smooth rational-quadratic splines, scale transformations and rotations to define a novel spherical normalizing-flow distribution whose parameters are predicted as a whole as the output of the transformer encoder. We test several structural choices diverting from the vanilla transformer architecture. In particular, we find dual residual streams, nonlinear QKV projection and a separate class token with its own cross-attention processing to boost test-time performance. The angular resolution for both showers and tracks improves substantially over the whole trained energy range from 100 GeV to 100 PeV. At 100 TeV deposited energy, for example, the median angular resolution improves by a factor of $1.3$ for throughgoing tracks, by a factor of $1.7$ for showers and by a factor of $2.5$ for starting tracks compared to state-of-the art likelihood reconstructions based on B-splines. While previous machine-learning (ML) efforts have managed to obtain competitive shower resolutions, this is the first time an ML-based method outperforms likelihood-based muon reconstructions above 100 GeV.
Tags
Links
- Source: https://arxiv.org/abs/2604.19846v1
- Canonical: https://arxiv.org/abs/2604.19846v1
Trouble viewing inline? Open PDF directly →
Full Text
140,239 characters extracted from source content.
Expand or collapse full text
Neural posterior estimation of the neutrino direction in IceCube using transformer-encoded normalizing flows on the sphere The IceCube Collaboration R. Abbasi1 M. Ackermann2 J. Adams3 J. A. Aguilar4 M. Ahlers5 J.M. Alameddine6 S. Ali7 N. M. Amin8 K. Andeen9 C. Argüelles10 Y. Ashida11 S. Athanasiadou2 S. N. Axani8 R. Babu12 X. Bai13 A. Balagopal V.8 S. W. Barwick14 V. Basu11 R. Bay15 J. J. Beatty16,17 J. Becker Tjus18 P. Behrens19 J. Beise20 C. Bellenghi21 S. Benkel2 S. BenZvi22 D. Berley23 E. Bernardini24 D. Z. Besson7 E. Blaufuss23 L. Bloom25 S. Blot2 F. Bontempo26 J. Y. Book Motzkin10 C. Boscolo Meneguolo24 S. Böser27 O. Botner20 J. Böttcher19 J. Braun28 B. Brinson29 Z. Brisson-Tsavoussis30 R. T. Burley31 D. Butterfield28 K. Carloni10 J. Carpio32,33 N. Chau4 Z. Chen34 D. Chirkin28 S. Choi11 A. Chubarov35 B. A. Clark23 G. H. Collin36 D. A. Coloma Borja24 A. Connolly16,17 J. M. Conrad36 D. F. Cowen37,38 C. De Clercq39 J. J. DeLaunay37 D. Delgado10 T. Delmeulle4 S. Deng19 P. Desiati28 K. D. de Vries39 G. de Wasseige40 T. DeYoung12 J. C. Díaz-Vélez28 S. DiKerby12 T. Ding32,33 M. Dittmer41 A. Domi35 L. Draper11 L. Dueser19 D. Durnford42 K. Dutta27 M. A. DuVernois28 T. Ehrhardt27 L. Eidenschink21 A. Eimer35 C. Eldridge43 P. Eller21 E. Ellinger44 D. Elsässer6 R. Engel26,45 H. Erpenbeck28 W. Esmail41 S. Eulig10 J. Evans23 P. A. Evenson8 K. L. Fan23 K. Fang28 K. Farrag46 A. R. Fazely47 A. Fedynitch48 N. Feigl49 C. Finley50 D. Fox37 A. Franckowiak18 S. Fukami2 P. Fürst19 J. Gallagher51 E. Ganster19 A. Garcia10 M. Garcia8 E. Genton10 L. Gerhardt52 A. Ghadimi25 C. Glaser6,20 T. Glüsenkamp50 J. G. Gonzalez8 S. Goswami32,33 A. Granados12 D. Grant53 S. J. Gray23 S. Griffin28 K. M. Groth5 D. Guevel28 C. Günther19 P. Gutjahr6 C. Ha54 A. Hallgren20 L. Halve19 F. Halzen28 L. Hamacher19 M. Handt19 K. Hanson28 J. Hardin36 A. A. Harnisch12 P. Hatch30 A. Haungs26 J. Häußler19 K. Helbing44 J. Hellrung18 B. Henke12 L. Hennig35 F. Henningsen35 L. Heuermann19 R. Hewett3 N. Heyer20 S. Hickford44 A. Hidvegi50 C. Hill21 G. C. Hill31 R. Hmaid46 K. D. Hoffman23 A. Hollnagel46 D. Hooper28 S. Hori28 K. Hoshina28 M. Hostert10 W. Hou26 M. Hrywniak50 T. Huber26 K. Hultqvist50 K. Hymon48 A. Ishihara46 W. Iwakiri46 M. Jacquart5 S. Jain28 O. Janik35 M. Jansson40 M. Jin10 N. Kamp10 D. Kang26 W. Kang55 A. Kappes41 L. Kardum6 T. Karg2 A. Karle28 A. Katil42 M. Kauer28 J. L. Kelley28 M. Khanal11 A. Khatee Zathul28 A. Kheirandish32,33 T. Kim56 H. Kimku54 F. Kirchner35 J. Kiryluk34 C. Klein2 S. R. Klein15,52 Y. Kobayashi46 S. Koch35 A. Kochocki12 R. Koirala8 H. Kolanoski49 T. Kontrimas21 L. Köpke27 C. Kopper35 D. J. Koskinen5 P. Koundal8 M. Kowalski2,49 T. Kozynets5 A. Kravka11 N. Krieger18 T. Krishnan10 K. Kruiswijk40 E. Krupczak12 A. Kumar2 E. Kun18 N. Kurahashi55 C. Lagunas Gualda21 L. Lallement Arnaud4 M. J. Larson23 F. Lauber44 J. P. Lazar40 K. Leonard DeHolton38 A. Leszczyńska8 C. Li28 J. Liao29 C. Lin8 Q. R. Liu53 Y. T. Liu38 M. Liubarska42 C. Love55 L. Lu28 F. Lucarelli57 W. Luszczak16,17 Y. Lyu15,52 M. Macdonald10 E. Magnus39 Y. Makino28 E. Manao21 S. Mancina24 A. Mand28 I. C. Mariş4 S. Marka58 Z. Marka58 L. Marten19 I. Martinez-Soler10 R. Maruyama59 J. Mauro40 F. Mayhew12 F. McNally60 K. Meagher28 A. Medina17 M. Meier46 Y. Merckx39 L. Merten18 J. Mitchell47 L. Molchany13 S. Mondal11 T. Montaruli57 R. W. Moore42 Y. Morii46 A. Mosbrugger35 D. Mousadi2 E. Moyaux40 T. Mukherjee26 M. Nakos28 U. Naumann44 J. Necker2 L. Neste50 M. Neumann41 H. Niederhausen12 M. U. Nisa12 K. Noda46 A. Noell19 A. Novikov8 A. Obertacke50 V. O’Dell28 A. Olivas23 R. Orsoe21 J. Osborn28 E. O’Sullivan20 B. Owens30 V. Palusova27 H. Pandya8 A. Parenti4 N. Park30 V. Parrish12 E. N. Paudel25 L. Paul13 C. Pérez de los Heros20 T. Pernice2 T. C. Petersen5 J. Peterson28 S. Pick2 M. Plum13 A. Pontén20 V. Poojyam25 B. Pries12 R. Procter-Murphy23 G. T. Przybylski52 L. Pyras11 C. Raab40 J. Rack-Helleis27 N. Rad2 M. Ravn20 K. Rawlins61 Z. Rechav28 A. Rehman8 I. Reistroffer13 E. Resconi21 S. Reusch2 C. D. Rho56 W. Rhode6 L. Ricca40 B. Riedel28 A. Rifaie44 E. J. Roberts31 S. Rodan62 M. Rongen35 A. Rosted46 C. Rott11 T. Ruhe6 L. Ruohan21 D. Ryckbosch43 J. Saffer45 D. Salazar-Gallegos12 P. Sampathkumar26 A. Sandrock44 G. Sanger-Johnson12 M. Santander25 S. Sarkar63 M. Scarnera40 M. Schaufel19 H. Schieler26 S. Schindler35 L. Schlickmann27 B. Schlüter41 F. Schlüter4 N. Schmeisser44 T. Schmidt23 A. Scholz21 F. G. Schröder8,26 S. Schwirn19 S. Sclafani23 D. Seckel8 L. Seen28 M. Seikh7 S. Seunarine62 P. A. Sevle Myhr40 R. Shah55 S. Shah22 S. Shefali45 N. Shimizu46 B. Skrzypek15 R. Snihur28 J. Soedingrekso6 D. Soldin11 P. Soldin19 G. Sommani18 C. Spannfellner21 G. M. Spiczak62 C. Spiering2 J. Stachurska43 M. Stamatikos17 T. Stanev8 T. Stezelberger52 T. Stürwald44 T. Stuttard5 G. W. Sullivan23 I. Taboada29 S. Ter-Antonyan47 A. Terliuk21 A. Thakuri13 M. Thiesmeyer28 W. G. Thompson10 J. Thwaites30 S. Tilav8 K. Tollefson12 J. A. Torres11 S. Toscano4 D. Tosi28 K. Upshaw47 A. Vaidyanathan9 N. Valtonen-Mattila18 J. Valverde9 J. Vandenbroucke28 T. Van Eeden2 N. van Eijndhoven39 L. Van Rootselaar6 J. van Santen2 J. Vara41 F. Varsi45 M. Venugopal26 M. Vereecken43 S. Vergara Carrasco3 S. Verpoest8 D. Veske58 A. Vijai23 J. Villarreal36 C. Walck50 A. Wang29 E. H. S. Warrick25 C. Weaver12 P. Weigel36 A. Weindl26 J. Weldert27 A. Y. Wen10 C. Wendt28 J. Werthebach6 M. Weyrauch26 N. Whitehorn12 C. H. Wiebusch19 D. R. Williams25 L. Witthaus6 G. Wrede35 X. W. Xu47 J. P. Yanez42 Y. Yao28 E. Yildizci28 S. Yoshida46 R. Young7 F. Yu10 S. Yu11 T. Yuan28 S. Yun-Cárcamo55 A. Zander Jurowitzki21 A. Zegarelli18 S. Zhang12 Z. Zhang34 P. Zhelnin10 and P. Zilberman28 1Department of Physics, Loyola University Chicago, Chicago, IL 60660, USA 2Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, D-15738 Zeuthen, Germany 3Dept. of Physics and Astronomy, University of Canterbury, Private Bag 4800, Christchurch, New Zealand 4Université Libre de Bruxelles, Science Faculty CP230, B-1050 Brussels, Belgium 5Niels Bohr Institute, University of Copenhagen, DK-2100 Copenhagen, Denmark 6Dept. of Physics, TU Dortmund University, D-44221 Dortmund, Germany 7Dept. of Physics and Astronomy, University of Kansas, Lawrence, KS 66045, USA 8Bartol Research Institute and Dept. of Physics and Astronomy, University of Delaware, Newark, DE 19716, USA 9Department of Physics, Marquette University, Milwaukee, WI 53201, USA 10Department of Physics and Laboratory for Particle Physics and Cosmology, Harvard University, Cambridge, MA 02138, USA 11Department of Physics and Astronomy, University of Utah, Salt Lake City, UT 84112, USA 12Dept. of Physics and Astronomy, Michigan State University, East Lansing, MI 48824, USA 13Physics Department, South Dakota School of Mines and Technology, Rapid City, SD 57701, USA 14Dept. of Physics and Astronomy, University of California, Irvine, CA 92697, USA 15Dept. of Physics, University of California, Berkeley, CA 94720, USA 16Dept. of Astronomy, Ohio State University, Columbus, OH 43210, USA 17Dept. of Physics and Center for Cosmology and Astro-Particle Physics, Ohio State University, Columbus, OH 43210, USA 18Fakultät für Physik & Astronomie, Ruhr-Universität Bochum, D-44780 Bochum, Germany 19I. Physikalisches Institut, RWTH Aachen University, D-52056 Aachen, Germany 20Dept. of Physics and Astronomy, Uppsala University, Box 516, SE-75120 Uppsala, Sweden 21Physik-department, Technische Universität München, D-85748 Garching, Germany 22Dept. of Physics and Astronomy, University of Rochester, Rochester, NY 14627, USA 23Dept. of Physics, University of Maryland, College Park, MD 20742, USA 24Dipartimento di Fisica e Astronomia Galileo Galilei, Università Degli Studi di Padova, I-35122 Padova PD, Italy 25Dept. of Physics and Astronomy, University of Alabama, Tuscaloosa, AL 35487, USA 26Karlsruhe Institute of Technology, Institute for Astroparticle Physics, D-76021 Karlsruhe, Germany 27Institute of Physics, University of Mainz, Staudinger Weg 7, D-55099 Mainz, Germany 28Dept. of Physics and Wisconsin IceCube Particle Astrophysics Center, University of Wisconsin—Madison, Madison, WI 53706, USA 29School of Physics and Center for Relativistic Astrophysics, Georgia Institute of Technology, Atlanta, GA 30332, USA 30Dept. of Physics, Engineering Physics, and Astronomy, Queen’s University, Kingston, ON K7L 3N6, Canada 31Department of Physics, University of Adelaide, Adelaide, 5005, Australia 32Department of Physics & Astronomy, University of Nevada, Las Vegas, NV 89154, USA 33Nevada Center for Astrophysics, University of Nevada, Las Vegas, NV 89154, USA 34Dept. of Physics and Astronomy, Stony Brook University, Stony Brook, NY 11794-3800, USA 35Erlangen Centre for Astroparticle Physics, Friedrich-Alexander-Universität Erlangen-Nürnberg, D-91058 Erlangen, Germany 36Dept. of Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA 37Dept. of Astronomy and Astrophysics, Pennsylvania State University, University Park, PA 16802, USA 38Dept. of Physics, Pennsylvania State University, University Park, PA 16802, USA 39Vrije Universiteit Brussel (VUB), Dienst ELEM, B-1050 Brussels, Belgium 40UCLouvain, Centre for Cosmology, Particle Physics and Phenomenology, CP3, Chemin du Cyclotron 2, 1348 Louvain-la-Neuve, Belgium 41Institut für Kernphysik, Universität Münster, D-48149 Münster, Germany 42Dept. of Physics, University of Alberta, Edmonton, Alberta, T6G 2E1, Canada 43Dept. of Physics and Astronomy, University of Gent, B-9000 Gent, Belgium 44Dept. of Physics, University of Wuppertal, D-42119 Wuppertal, Germany 45Karlsruhe Institute of Technology, Institute of Experimental Particle Physics, D-76021 Karlsruhe, Germany 46Dept. of Physics and The International Center for Hadron Astrophysics, Chiba University, Chiba 263-8522, Japan 47Dept. of Physics, Southern University, Baton Rouge, LA 70813, USA 48Institute of Physics, Academia Sinica, Taipei, 11529, Taiwan 49Institut für Physik, Humboldt-Universität zu Berlin, D-12489 Berlin, Germany 50Oskar Klein Centre and Dept. of Physics, Stockholm University, SE-10691 Stockholm, Sweden 51Dept. of Astronomy, University of Wisconsin—Madison, Madison, WI 53706, USA 52Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA 53Dept. of Physics, Simon Fraser University, Burnaby, BC V5A 1S6, Canada 54Dept. of Physics, Chung-Ang University, Seoul 06974, Republic of Korea 55Dept. of Physics, Drexel University, 3141 Chestnut Street, Philadelphia, PA 19104, USA 56Dept. of Physics, Sungkyunkwan University, Suwon 16419, Republic of Korea 57Département de physique nucléaire et corpusculaire, Université de Genève, CH-1211 Genève, Switzerland 58Columbia Astrophysics and Nevis Laboratories, Columbia University, New York, NY 10027, USA 59Dept. of Physics, Yale University, New Haven, CT 06520, USA 60Department of Physics, Mercer University, Macon, GA 31207-0001, USA 61Dept. of Physics and Astronomy, University of Alaska Anchorage, 3211 Providence Dr., Anchorage, AK 99508, USA 62Dept. of Physics, University of Wisconsin, River Falls, WI 54022, USA 63Dept. of Physics, University of Oxford, Parks Road, Oxford OX1 3PU, United Kingdom analysis@icecube.wisc.edu Abstract IceCube is a cubic-kilometer-scale neutrino detector located at the geographic South Pole. A precise directional reconstruction of IceCube neutrinos is vital for associations with astronomical objects. In this context, we discuss neural posterior estimation of the neutrino direction via a transformer encoder that maps to a normalizing flow on the 2-sphere. It achieves a new state-of-the-art angular resolution for the two main event morphologies in IceCube - tracks and showers - while being significantly faster than traditional B-spline-based likelihood reconstructions. All-sky scans can be performed within seconds rather than hours, and take constant computation time, regardless of whether the posterior extent is arc-minutes or spans the whole sky. We utilize a combination of C2C^2-smooth rational-quadratic splines, scale transformations and rotations to define a novel spherical normalizing-flow distribution whose parameters are predicted as a whole as the output of the transformer encoder. We test several structural choices diverting from the vanilla transformer architecture. In particular, we find dual residual streams, nonlinear QKV projection and a separate class token with its own cross-attention processing to boost test-time performance. The angular resolution for both showers and tracks improves substantially over the whole trained energy range from 100 GeV to 100 PeV. At 100 TeV deposited energy, for example, the median angular resolution improves by a factor of 1.31.3 for throughgoing tracks, by a factor of 1.71.7 for showers and by a factor of 2.52.5 for starting tracks compared to state-of-the art likelihood reconstructions based on B-splines. While previous machine-learning (ML) efforts have managed to obtain competitive shower resolutions, this is the first time an ML-based method outperforms likelihood-based muon reconstructions above 100 GeV. 1 Introduction The IceCube neutrino detector [36] is the world’s largest high-energy neutrino observatory and has made several breakthrough discoveries over the last decade. Among them, we find the association of high-energy neutrinos with specific astrophysical sources including the active galactic nucleus NGC 1068 [43] and the Galactic Plane [47]. These results were achieved through statistical analyses that rely fundamentally on two observables reconstructed for each individual neutrino event: the deposited energy and the arrival direction. The arrival direction is particularly important for point-source searches, where the flux sensitivity scales approximately with the angular resolution through improved background suppression. We focus on the direction reconstruction in this paper. 1.1 Motivation IceCube detects neutrinos via Cherenkov light [3] which is emitted from charged relativistic particles that originate in the neutrino interaction. The traditional approach to inferring the direction is a maximum-likelihood estimation based on the observed Cherenkov light distributions across the detector’s optical modules. Over the years, dedicated reconstructions have been developed for the distinct event morphologies that arise from different neutrino interaction channels. Three major morphology classes are useful to be differentiated in this context: shower-like events, throughgoing tracks, and starting tracks. Shower-like events, produced by charged-current νe _e interactions, neutral-current interactions of all flavors, and sub-PeV charged-current ντ _τ interactions, exhibit approximately spherical Cherenkov light patterns. The angular resolution of showers is typically on the order of 5 to 10 degrees and the associated uncertainty contours are often non-Gaussian. This class of events was instrumental in the discovery of the astrophysical diffuse neutrino flux [35] and neutrino emission from the Galactic Plane [47]. The Galactic Plane measurement relied on a neural-network-based likelihood approach for shower reconstruction [41], combined with CNN-based regression networks [39] in the event selection stage. Recently, a B-spline-based likelihood reconstruction achieved state-of-the-art angular resolution for showers [49] by modeling the shower extension and improving the treatment of systematic uncertainties. We refer to this method as Taupede2024 for the remainder of the paper. Throughgoing track events originate from charged-current νμ _μ interactions that occur outside the instrumented volume, producing muons that traverse the full detector. The resulting elongated light patterns typically yield angular resolutions below one degree, making this the most important event morphology for point-source studies. The primary reconstruction algorithm, referred to as SplineMPEMax in the following, models the muon track as an effective averaged energy-loss pattern [32] and has served as the backbone of all point-source analyses to date, including the identification of TXS 0506+056 [38] and NGC 1068 [43] as neutrino sources. At energies above a few TeV, and especially beyond 100 TeV, stochastic energy losses begin to dominate and the quasi-continuous modeling assumption becomes less accurate. An extension that explicitly models stochastic losses along the track was developed in [40]. However, its substantially longer runtime, numerical instabilities, and only marginal resolution improvements have so far precluded its adoption in practice. Starting track events, in which the νμ _μ interaction vertex lies within the instrumented volume, present a distinct reconstruction challenge. The Cherenkov light from the initial hadronic cascade overlaps with that of the typically short outgoing muon track, complicating directional inference. Nevertheless, this event class holds substantial promise for southern-sky neutrino astronomy in IceCube. Atmospheric muon veto techniques [35] can yield highly pure starting-track samples with a direct line of sight toward the Galactic Center, and the angular resolution is considerably better than that of pure showers. A recent southern-sky diffuse flux measurement incorporating a dedicated starting-track selection [48] has demonstrated their potential for southern-sky measurements. However, the full power of starting tracks to resolve Galactic structure has yet to be utilized. A unique challenge for reconstruction in IceCube, compared to artificially constructed detectors, is that the detection medium is a naturally formed glacier. While the ice at several kilometers depth is exceptionally transparent on average, inhomogeneous dust deposits introduce position-dependent variations in scattering and absorption [51] that must be accounted for in the reconstruction. Most of the aforementioned likelihood-based approaches rely on B-spline parameterizations of the expected Cherenkov photon arrival time distributions [55]. Incorporating the high-dimensional ice optical properties into these models is challenging. Beyond systematic uncertainties, complex event signatures such as high-energy tracks possess intrinsic degrees of freedom, notably the stochastic energy-loss profile, that are difficult to model in their own right [34][40]. A directional likelihood estimation must jointly profile over all these combined nuisance parameters, which is both numerically challenging and computationally expensive. An alternative to B-spline parameterizations that has emerged in recent years is to model the likelihood using neural networks [41]. This approach, combined with small feed-forward networks for the event selection stage, enabled the first measurement of neutrinos from the Galactic Plane [47]. Neural-network-based likelihoods typically offer an easier incorporation of systematic uncertainties into the modeling process. However, the challenge of parameterizing complex intrinsic degrees of freedom, such as the muon energy-loss profile, persists, and the computational cost of likelihood evaluation remains comparable to or greater than that of B-spline methods. In general, a low processing time is desirable for many applications. It is particularly important for real-time neutrino alerts [37] that are sent out to the astronomical community. Here, the goal is to obtain precise localization contours as fast as possible. Currently, expensive profile-likelihood scans can take up to several hours per event, which delays localization for astronomical follow-up observations by other experiments. One way to get a fast prediction is to switch from an explicit likelihood model to a neural-network regression model that operates in the inverse direction, i.e. predicts the direction from the data. However, such inverse models typically either perform point predictions or assume Gaussian uncertainties [39], which can lead to biases or a lack of precision. In this paper, we describe a method that also operates in the inverse direction but goes beyond point estimation and Gaussian posterior assumptions to predict the full directional posterior distribution. Rather than explicitly parameterizing nuisance degrees of freedom in the forward direction, the method implicitly learns to marginalize over them during training. At inference time, the complete posterior is obtained in a single neural-network forward pass, with the computational cost effectively offloaded to the training stage. The result is a reconstruction framework that combines the speed of direct inference with precise angular resolution, robust non-Gaussian uncertainty quantification, and the flexibility to incorporate complex systematic effects. Because it is both fast and precise, it is applicable to event selections, final-level analysis reconstructions, and real-time astronomical alerts. 1.2 Outline The proposed inference procedure is based on amortized neural posterior estimation (NPE) with conditional normalizing flows. In standard neural posterior estimation [29], one is typically interested in obtaining the posterior p(θ|xo)p(θ|x_o) for one particular event with observed data xox_o. The procedure involves iterative re-simulations of both a proposal prior and a posterior approximation, and is situated in the broader framework of simulation-based inference (SBI) [5]. Here, we are interested in all posteriors at the same time, and just apply one single approximate posterior training. Hence it is a form of amortized neural posterior estimation since the cost of estimation of a single posterior is amortized by a neural network. Conditional normalizing flows amortized with neural networks have previously shown promise for reconstructing the posterior distribution of the energy and direction of neutrinos with high fidelity [46], without being bound to Gaussian uncertainty assumptions. In that prior work we used a compression of the per-module photon data into a summary statistic that is fed into a graph neural network (GNN) encoding similar to [44]. It used a k-nearest neighbor approach to construct the graph followed by edge convolution [54]. The output of the encoding was then mapped to the normalizing-flow parameters. Here, we improve on that work for the directional reconstruction in several ways. First, we change the GNN to a transformer [53] encoding. Secondly, we develop a new efficient manifold normalizing-flow on the 2-sphere which is numerically stable in this conditional setting and is able to handle the extra requirement of modeling several orders of magnitude in scale simultaneously. To achieve this, we extend the recursive construction in [30] which relies on rational-quadratic splines [8]. We add a von-Mises-Fisher scaling function, remove the recursive conditioning scheme and make the rational-quadratic spline flow C2C^2-continuous. A schematic overview of the proposed algorithm is shown in Fig. 1. Depicted are the detected photon counts within the IceCube neutrino detector for a track-like example and a shower-like example event. The detector has a quasi-hexagonal shape and each detection unit or digital optical module (DOM) is indicated by a black dot. DOMs that are detecting photons are color coded with early-hit modules shown in red and late-hit modules shown in blue. The photon data is encoded via a transformer and mapped to the normalizing-flow parameters. More details of the detector and encoding pipeline are given in section 4. topside IceCube event (track)TransformerPosterior (localized) (a) track example topside IceCube event (shower)TransformerPosterior (all-sky) (b) shower example Figure 1: Schematic overview of transformer-based amortized neural posterior estimation for neutrino event directions in IceCube. Two example posteriors are depicted: a posterior (orthographic projection) localized on a scale of a few degrees for a track (a) and a more extended posterior (Mollweide projection) for a shower (b). The photon data is collected in optical modules that are hexagonally distributed in a cubic-kilometer volume. The color of the modules indicates time, where red is early and blue is late. The transformer maps the event-dependent photon data to the normalizing-flow parameters and thereby encodes all posteriors for events similar to the training dataset at once. The paper is structured as follows. In section 2 we introduce amortized neural posterior estimation. In section 3 we introduce normalizing flows on the 2-sphere, describe our new manifold normalizing flow and showcase the advantage of a manifold normalizing flow for astronomical skymap creation. In section 4 we describe the overall architecture, including the data preparation from raw IceCube data. In section 4.2, in particular, we motivate the switch from a GNN to a transformer encoding, introducing the concept of “Bayesian inductive bias”. In section 5 we describe the training procedure and introduce the concept of “probabilistic regularization”. Finally, in section 6 we show results and end with a discussion in section 7. 2 Amortized neural posterior estimation - a form of “ELBO-free” stochastic variational inference We can write Bayes’ theorem connecting observations x, for example data in the IceCube detector, and parameters θ as p(θ|x)=p(x|θ)⋅p(θ)p(x),p(θ|x)= p(x|θ)· p(θ)p(x), (1) where we call p(θ|x)p(θ|x) the posterior and p(x|θ)p(x|θ) the “data generating distribution” when seen as a probability distribution function (PDF)111We assume we have continuous data for simplicity and thus call it a “PDF”. The argument is of course also valid for discrete data which would technically be a probability mass function over x for fixed θ or “likelihood function” when seen as a function ℒ(θ)L(θ) of θ. For IceCube, we have independent photon observations per DOM and per photon, so we can write ℒ(θ)=∏j=1NDOM∏i=1Np,jpj(xi|θ),L(θ)= _j=1^N_DOM _i=1^N_p,jp_j(x_i|θ), (2) with number of DOMs NDOMN_DOM and number of photons for the jjth DOM Np,jN_p,j. The likelihood is a product of independent and identically distributed data (IID) factors per DOM. A standard approach to obtain the posterior from a measurement outcome is via Markov-Chain Monte Carlo (MCMC) sampling [31] of an objective function that is proportional to the likelihood function. However, this does not scale to high dimensions very well. In the 1990s, variational inference [20] emerged as an alternative to Markov-Chain Monte Carlo for Bayesian analysis, where the posterior solution is found via optimization instead of sampling. This “traditional” variational inference approach requires the Evidence Lower Bound (ELBO) to be formulated as a loss function. While this in principle allows the handling of higher dimensional scenarios, it is still inefficient in the sense that it requires a full ELBO optimization for a single measurement and corresponding posterior inference. In recent years, neural posterior estimation [29] has emerged as a way to directly approximate the posterior distribution of interest without an ELBO objective - hence one can also call it "ELBO-free" variational inference. It involves an iterative scheme of sampling from a “proposal prior” p~(θ) p(θ), running a simulator p~(x|θ) p(x|θ) to produce data “x” using the samples from the proposal prior as input parameters, and fitting an approximate posterior qψ(θ|x)q_ψ(θ|x) parametrized by a neural network via ψ. Here and in the following, x refers to all data in an event, not individual photons as in eq. 2. The procedure is repeated until the approximate posterior converges to the true posterior for a specific event datum xox_o. It is part of the broader scheme of simulation based inference [5]. One can also decide to just run a single training run and fit qψ(θ|x)q_ψ(θ|x) on samples drawn from the initial proposal prior, and thereby approximate all posteriors of the dataset at once. It assumes the proposal prior is sufficient or its influence becomes negligible once data is observed. This is a form of “amortized” neural posterior estimation, as all possible posterior predictions are amortized by a neural network with no extra refinement. We adapt this procedure here, since simulation data for IceCube is costly to produce and adaptive re-simulations typically not feasible. The loss function ℒNPE(ψ)L_NPE(ψ) in this scheme can be derived from the forward Kullback-Leibler (KL) divergence [22] between the joint distribution in the Monte Carlo simulation p~(x,θ)=p~(θ)⋅p~(x|θ) p(x,θ)= p(θ)· p(x|θ), from which we have samples in the simulation, and an approximation of the joint distribution that involves the posterior qψ(θ|x)q_ψ(θ|x), over the neural network parameters ψ. The KL divergence is a well-known quantity from information theory and is zero if the involved probability distributions are identical and positive otherwise. The implicit minimization of this quantity is performed via a Monte Carlo approximation of the integral and subsequent dropping of constant terms in the optimization parameters ψ as ψ⋆ ψ =argminψDKL(p~(x,θ)∥q(x)qψ(θ|x)) = _ψ\;D_KL\! ( p(x,θ)\; \|\;q(x)\,q_ψ(θ|x) ) (3) =argminψp~(x,θ)[logp~(x,θ)q(x)qψ(θ|x)] = _ψ\;E_ p(x,θ)\! [ p(x,θ)q(x)q_ψ(θ|x) ] (4) =argminψp~(x)[DKL(p~(θ|x)∥qψ(θ|x))]+const = _ψ\;E_ p(x) [D_KL\! ( p(θ|x)\; \|\;q_ψ(θ|x) ) ]+const (5) =argminψp~(x,θ)[−logqψ(θ|x)](drop terms constant in ψ) = _ψ\;E_ p(x,θ)\! [- q_ψ(θ|x) ] (drop terms constant in ψ) (6) ≈argminψ1N∑iN−logqψ(θi|xi)≡argminψℒNPE(ψ). ≈ _ψ\; 1N _i^N- q_ψ( _i|x_i)≡ _ψ\;L_NPE(ψ). (7) It can be seen that the loss function is a Monte Carlo approximation (eq. 7) of the conditional cross entropy (eq. 6) and shares the same minimum with the expected KL-divergence between the implicit simulation posterior p~(θ;x) p(θ;x) and the posterior approximation qψ(θ|x)q_ψ(θ|x) (eq. 5) because constant terms can be dropped. Since N is typically in the millions, in practice it is optimized stochastically over batches as in standard neural network training. The target parameters θ of the posterior can be chosen to be any subset of simulated parameters of interest. In our case we focus on the zenith and azimuth angles, while ignoring other parameters, such as the event position or energy, whose contributions appear as constant terms in eq. 5. Implicitly, these other parameters act as nuisance parameters and are marginalized in the result obtained by optimizing eq. 7. The choice for the posterior PDF is not predetermined and can be a mixture model as in the original NPE paper [29], for example. However, if a normalizing flow is chosen as the posterior approximation, it can be shown [14] that it not only allows efficient amortized posterior estimation to be performed, it also strictly generalizes supervised neural network training with a Mean-Squared-Error (MSE) loss. The latter can be viewed as a special case of NPE with a restricted affine normalizing flow. We therefore get a direct link between scalable Bayesian analysis and traditional MSE-based supervised regression for parameter prediction. An overview of these relationships is shown in Fig. 2. Neural Posterior Estimation Amortized neural posterior estimation with normalizing flows Supervised learning of parameters (MSE-loss) 1) PDF → normalizing flow 2) single fit, no iterative procedure affine normalizing flow with identity scaling Figure 2: Relation of amortized neural posterior estimation with normalizing flows (used here) to generic neural posterior estimation and standard supervised regression of parameters. 3 Conditional normalizing flows on the 2-sphere In order to perform neural posterior estimation for the direction of a neutrino, we need a conditional probability distribution on the 2-sphere. We base this distribution on a manifold normalizing flow. 3.1 Technical definitions A normalizing flow describes a PDF pu(θ)p_u(θ) using a diffeomorphic mapping θ=fu(zb)θ=f_u(z_b) with parameters u and so-called “base distribution” p0(zb)p_0(z_b) with a change-of-variable formula pu(θ)=p0(fu−1(θ))⋅|detJu−1(θ)|,p_u(θ)=p_0(f_u^-1(θ))·|detJ_u^-1(θ)|, (8) where JuJ_u is the Jacobian of the mapping Ju=dfu(zb)dzbJ_u= df_u(z_b)dz_b. Here, we call “zbz_b” the auxiliary base space and “θ” the target space since we are interested in describing a posterior over parameters θ. Eq. 8 then allows access to a valid probability distribution over θ in closed form. 3.2 Manifold distributions For manifold distributions it makes sense in this context to differentiate between intrinsic coordinates θint _int and embedding coordinates θemb _emb. For the 2-sphere, the intrinsic coordinates are the zenith and azimuth angles, i.e. θint=θzen,θazi _int=\ _zen, _azi\, while the embedding coordinates are the corresponding Euclidean embedding coordinates θemb=x,y,z _emb=\x,y,z\. Eq. 8 can then be extended to define a distribution over a non-Euclidean manifold [30], for example the 22-sphere, by a modification of the Jacobian determinant via |detJu−1(θ)| |detJ_u^-1(θ)| →det((Ju−1⋅P)T⋅Ju−1⋅P)[θemb] → det ((J_u^-1· P)^T· J_u^-1· P )[ _emb] (9) ≡detu−1[θemb]. ^-1_u[ _emb]. (10) Here, P(θemb)P( _emb) is a projection matrix that consists of the orthogonal vectors in the tangent plane at θemb _emb, which in the case of the sphere are two orthogonal vectors. The term Ju−1(θemb)J_u^-1( _emb) is the inverse Jacobian defined and interpreted in embedding coordinates. One can think of this modified structure as first projecting the Jacobian into the tangent space before taking the determinant. The final determinant then measures the local area change in the tangent space on the manifold. 3.3 Conditional PDFs We now use a neural network mapping u=gψ(X)u=g_ψ(X) to predict the flow function parameters u from data X and thereby describe a conditional PDF on the 2-sphere as pψ(θemb|X)=p0(fgψ(X)−1(θemb))⋅detgψ(X)−1[θemb].p_ψ( _emb|X)=p_0(f_g_ψ(X)^-1( _emb))·detJ^-1_g_ψ(X)[ _emb]. (11) This conditioning scheme requires an efficient normalizing flow that should not have more than ∼102–103 10^2--10^3 parameters since all of them are predicted by gψ(X)g_ψ(X). The function gψ(X)g_ψ(X) could be any neural network encoding architecture that is appropriate to handle the input data. As described in section 4.2 we use a transformer encoding for gψ(X)g_ψ(X), which among other things allows the handling of variable-sized input data and encodes the permutation symmetry of arguments to the data generation distribution, as discussed in section 4.2. For convenience, we will in some cases denote θemb _emb simply as θ. 3.4 A new normalizing flow The flow function fu(zb)f_u(z_b) we use is motivated by the conditional flow scheme described in [30] which describes a conditional flow on the 2-sphere by transforming the sphere into a cylinder. This construction has the advantage that the Jacobian determinant of the mapping from spherical coordinates to cylinder coordinates has determinant 1. It is then enough to just calculate contributions to the determinant from transformations on the cylinder. Let us call this transformation from spherical embedding coordinates to cylindrical embedding coordinates the “cylinder transformation” C(θemb)C( _emb). We can write it as C(θemb)=C(x,y,z)=(ρϕz)=(1atan(y/x)z).C( _emb)=C(x,y,z)= ( array[]cρ\\ φ\\ z array )= ( array[]c1\\ atan(y/x)\\ z array ). (12) The Jacobian determinant of this transformation is 1, even without invoking eq. 9. Next we apply smooth rational-quadratic splines on the z and ϕφ coordinates, respectively. Rational-quadratic splines (RQS) [8] are diffeomorphic functions defined on any interval of choice, and with appropriate boundary conditions can be turned into a diffeomorphism on the circle [30]. Therefore, they serve as flexible normalizing-flow generators on the cylinder height interval [−1,1][-1,1] and the cylinder angle [0,2π][0,2π]. A problem that arises with those functions is that they are not C2C^2 smooth, which usually leads to unphysical features in the resulting normalizing-flow PDF (see Fig. 15). We use constrained smooth splines that we derive in appendix B to get rid of such features. In this way we obtain a smooth RQS on the interval [−1,1][-1,1], rqsI(z)rqs_I(z), which we use for z. We also obtain a periodic smooth RQS, rqsA(ϕ|z)rqs_A(φ|z), on the angle [0,2π][0,2π] which we use for ϕφ. The angular RQS is conditioned on z. We condition it differently than proposed in [30] by using a fixed polynomial spline interpolator instead of a neural network (see Fig. 16). By an appropriate blending toward an identity mapping, this construction avoids singular features near the polar regions, an outcome that smoothing alone cannot accomplish. Further details are given in appendix B. Next we apply a scaling transformation that corresponds to the normalizing flow depiction of a von-Mises-Fisher (vMF) distribution. The vMF distribution [11] is a common symmetric distribution on the 2-sphere and can be written as a normalizing flow by starting with a flat base distribution, transforming to the cylinder, and applying the scaling function [19]222We use 1+z2 1+z2 which is equivalent to a uniform random variable as used in [19] since the flat distribution on the sphere corresponds to a uniform distribution on the height z of the cylinder. F(z)=1+(1/κ)⋅ln(1+z2+(1−1+z2)⋅e−2κ)F(z)=1+(1/κ)·ln ( 1+z2+ (1- 1+z2 )· e^-2κ ) (13) on the height z∈[−1,1]z∈[-1,1], followed by an inverse cylinder transformation and a rotation. The scaling transformation F(z)F(z), together with the cylindrical mappings, allows one to effectively “zoom in” on a particular region. For neutrinos, the expected directional posterior regions can span several orders of magnitude down to a fraction of a degree, which is challenging to model without such scaling functions. As final steps, we transform back from the cylinder to the sphere with C−1C^-1 and apply a rotation RR that we parametrize via householder reflections [17] (see eq. 37 in appendix C). The combination of all those flows written in a single block i yields fϕi(θemb)=[Ri∘C−1∘Fi∘rqsϕ,i∘rqsI,i⏟oncyl.coords.∘C](θemb),f_ _i( _emb)= [R_i\ \ C^-1\ \ F_i\ \ rqs_φ,i\ \ rqs_I,i_on\ cyl.\ coords.\ \ C ]( _emb)\ , (14) where ϕi _i denotes all parameters of the involved functions bundled together. It is to be read from right to left in order and can be interpreted as a flexible functional building block that maps coordinates on the sphere θemb=(x,y,z) _emb=(x,y,z) to new coordinates θemb′ _emb in an invertible manner. The Jacobian determinant of fi(θemb)f_i( _emb) is available in closed form and can be efficiently calculated since the determinant of the Jacobian of CC, C−1C^-1 and the rotation RiR_i is 1, while the Jacobian-determinants of the rational quadratic splines and the scaling function F(z)F(z) are absolute values of simple one-dimensional derivatives. For the final normalizing flow, we use 15 blocks of fϕi(θemb)f_ _i( _emb) (eq. 14) chained together, i.e. ftot(θemb)=[fϕ15∘…∘fϕ1](θemb),f_tot( _emb)= [f_ _15 … f_ _1 ]( _emb)\ , (15) which results in an expressive flow that can be multimodal, can express multiple orders of magnitude in scale due to the vMF scalings, and is efficient, i.e. it is described by a few 100 parameters packaged in (Φ1,…,Φ15)( _1,…, _15). We can therefore easily transform it into a conditional PDF by predicting those parameters with a neural network using eq. 11, which amortizes those flow parameters with the neural-network parameters ψ and yields the amortized flow function ftot,ψ(θemb)f_tot,ψ( _emb) (see Fig. 4(b) and Fig. 4(c) for an illustration). As depicted there, we also add a projection function fproj.(zb)f_proj.(z_b) in the very beginning in practice, which performs a non-learnable fixed stereographic projection from the 2-d plane to the sphere which has been motivated in [14] and effectively allows starting with a Gaussian base distribution. An alternative is to start with the flat distribution on the sphere as a base distribution and drop this projection. The flow has been implemented in the open-source github package jammy_flows [13]. 3.5 Sampling and skymaps Besides flexible PDF modelling, normalizing flows on the sphere offer a new capability of very efficient astronomical skymap creation. The bijective normalizing flow mapping ftot,ψf_tot,ψ transforms samples from the base space to the target space and thereby draws samples from the normalizing-flow distribution on the sphere. These samples can be used to define an adaptive grid on the 2-sphere via an adaptive HEALPIX binning. HEALPIX [15] is an equal-binning scheme on the 2-sphere that is widely used in astronomy. It can be adapted for an irregular binning scheme since pixels can be divided into equal-sized smaller pixels on demand, for example via the multi-order coverage [10] (MOC) format as implemented in the package mhealpy [26] which we utilize in the following. In the normalizing flow context, the irregular MOC binning can be created on the fly from the samples, which will automatically yield finer binning in regions where more detail is necessary. We impose the constraint that a maximum number of samples NmaxN_max is allowed in a given pixel, otherwise it is subdivided into 4 sub-pixels. On this irregular grid one can then employ the change-of-variable formula (eq. 8) to obtain the exact PDF value and draw smooth contour lines. The step-by-step procedure from samples to irregular grid and PDF evaluation is indicated in Fig. 3. (a) Full-sky posterior (b) Localized posterior Figure 3: Illustration of the constant-time skymap creation using normalizing flows irrespective of shape or size of the posterior. Drawing N samples from the posterior (left) defines adaptive multi-order coverage grid cells (center) via subdivision based on maximally allowing NmaxN_max samples per cell. We use N=10000N=10000 and Nmax=5N_max=5. At the cell locations the exact probability can be evaluated which produces smooth PDF maps and contours (right). In practice we find that ∼10000 10000 samples with NmaxN_max between 55 and 1010 work well to obtain a smooth PDF in constant time for the normalizing flow described in section 3.4, independent of the absolute size of the contour. While the binning is subject to variations due to the stochastic sampling process, numerical differences of the drawn contours are not noticable by eye with these settings. The scheme works without iterative brute-force scanning and is in particular useful for contour regions that are orders of magnitude smaller compared to the full sky as well as for irregular PDFs with disconnected regions. 4 Combined architecture for IceCube data Neutrino interactions lead to relativistic charged secondary particles like electrons and muons which emit Cherenkov photons [3] in transparent matter like the Antarctic ice. IceCube is a cubic-kilometer neutrino detector located at the South Pole. It consists of over 5000 digital optical modules (DOMs) embedded in the ice that detect such incoming Cherenkov photons with photomultiplier tubes (PMTs) [36]. A depiction of the hexagonal layout of these DOMs with two example events is shown in Fig. 1. 4.1 Data pre-processing The photons that reach a DOM have a specific time arrival distribution as depicted in Fig. 4(a). This distribution technically corresponds to reconstructed “effective photons” which come from a least-squares unfolding algorithm that involves the DOM detector responses [23]. In the following discussion we will not differentiate between reconstructed photons and true physical photons for simplicity. The shape of the distribution and the number of photons depend on event parameters like the flavor, energy or direction of the incoming neutrino. In a first step, noise is cleaned with an established cleaning algorithm. DOMs with detected photons that cannot causally be connected to a main “cluster” of DOMs are removed in this step. It is followed by a PMT afterpulse cleaning in which we only keep photons in a given DOM that arrive at most 5 microseconds after the first hit in the same DOM - this excludes late after pulses that typically arrive on a timescale of 6 microseconds or later after a physical photon. We found in various trainings consistently better results by applying such a combined cleaning first versus taking the raw effective photon data as input for the neural network. We then encode the surviving cleaned observed photon distribution in a fixed-length summary statistic following the ansatz originally described for a convolutional neural network (CNN) encoding [39]. On top of time quantiles and charge quantiles [39], we also add the absolute DOM position and the DOM position relative to the charge-weighted “center of gravity” position of the event. We furthermore add an identifier whether a PMT has normal or high quantum efficiency and differentiate between normal, unhit and saturated PMT DOMs via an extra 3-d one-hot encoding. In unhit or “empty” DOMs we set all the charge and timing information to zero. These are to be differentiated from the few percent of DOMs that are malfunctioning and even in principle would not be able to detect any photons, which are altogether excluded here. In saturated DOMs we remove all the charge and timing information since saturation is not well modeled in our simulation. We only keep the time of the first photon hit, which is reliable for saturated DOMs. Time entries are shifted to be aligned relative to the median time of the event and normalized by 1/10000.01/10000.0. Charge entries are calculated as Q~=ln(1+Q)⋅0.2 Q=ln(1+Q)· 0.2 which automatically maps zero to zero. Position entries and relative position entries are normalized by 1/500.01/500.0, e.g. x~=1500⋅x x= 1500· x. Examples of summary statistics are given in Fig. 4(a). In total, this per-DOM vector is 27-dimensional in all experiments and we denote it with TiT_i for DOM i. Photons hit DOM Reconstructed photon distributionNormal DOM →Ti=(x~,y~,z~,cx~,cy~,cz~,→~,T~1st,→~,PMTtype,1,0,0)\ \ \ → T_i=( x, y, z, c_x, c_y, c_z, Q, T_1st, T_o,PMT_type,1,0,0)Saturated DOM →Ti=(x~,y~,z~,cx~,cy~,cz~,→,T~1st,→,PMTtype,0,1,0)→ T_i=( x, y, z, c_x, c_y, c_z, 0, T_1st, 0,PMT_type,0,1,0)Empty DOM →Ti=(x~,y~,z~,cx~,cy~,cz~,→,0,→,PMTtype,0,0,1)\ \ \ \ → T_i=( x, y, z, c_x, c_y, c_z, 0,0, 0,PMT_type,0,0,1) (a) The top figure shows a photon arrival distribution with example time (“T”) and charge (“Q”) summary statistics. The bottom panel shows the different ways that this photon summary statistic is combined with other DOM-specific information into a vector representation TiT_i. The DOM position is denoted by x,y,zx,y,z, the difference of the DOM position to the center of gravity of the event by cx,cy,czc_x,c_y,c_z, charge quantities of the observed photon hits by → Q and temporal quantities by T1stT_1st and → T_o. PMTtypePMT_type is either 0 (normal PMT) or 1 (high-QE PMT). The final 3-dimensional one-hot part differentiates between “normal”, “saturated” and “empty” DOMs. The “∼ ” above a quantity represents variable-dependent rescaling. …T2T_2T1T_1TNT_NTransformer Encoder…T~2 T_2T~1 T_1T~N T_NAggregateMLPMLPFlow params: (ϕ1→,ϕ2→,…,ϕM→)[ψ]( _1, _2,…, _M)[ψ]Flow function: ftot,ψ(z)=[fϕM∘fϕM−1∘…∘fϕ1∘fproj.](z)f_tot,ψ(z)=[f_ _M f_ _M-1 … f_ _1 f_proj.](z) ψ (b) Encoding structure from DOM vector representations TiT_i into the flow function ftot,ψ(z)f_tot,ψ(z). See section 3.4 for more details on the flow definition. ftot,ψ(z)f_tot,ψ(z)ftot,ψ−1(θ)f_tot,ψ^-1(θ)p0(z)p_0(z)auxiliary base spacePosterior: pψ(θ|T1,T2,…,TN)p_ψ(θ|\T_1,T_2,…,T_N\)target space “θ” (c) The encoded flow function ftot,ψ(z)f_tot,ψ(z) is used to define the conditional posterior. It implicitly depends on the neural network parameters ψ and on the DOM summary statistics TiT_i. Figure 4: The data encoding pipeline from photon hits to posterior prediction. 4.2 Transformer encoding The transformer architecture [53] is a machine learning model that has gained traction over the past years. Originally used in neural language processing [53] for sequential language data in an encoder-decoder setting, its use case has expanded to all other data science due to its universality. Here, the transformer encoding is used as a compression algorithm for the IceCube data. The transformer encoding consists of a number of building blocks or layers NLN_L, where we use NL=20N_L=20 in all experiments. It operates on a set of vectors instead of a single vector, which in our case are the DOM embeddings TiT_i. By default, every layer uses the vanilla architecture with “pre”-layer normalization [57], followed by a self-attention block that exhanges information between tokens i and j and a final multi-layer perceptron (MLP) without dropout, that is applied per token i. For each self-attention block, and given the input vector XiX_i for each DOM i at a given layer L, we first define the linear query (Q), key (K) and value (V) embeddings for each token, sometimes also called “QKV embedding”: qi=Xi⋅WQ,ki=Xi⋅WK,vi=Xi⋅WV.q_i=X_i· W_Q, k_i=X_i· W_K, v_i=X_i· W_V. (16) The self-attention update between query i and key j is then given as αij=softmaxj(qikj⊤dk) _ij=softmax_j\! ( q_ik_j d_k ) (17) and the output is multiplied by the value embeddings as Xi′=[∑j=1nαijvj]⋅WOX_i = [ _j=1^n _ij\,v_j ]· W_O (18) with a final out-projection matrix WOW_O. The result is further processed per token via an MLP. Further details like the addition of an explicit positional encoding can be found in [53]. We also test several extensions to the vanilla transformer which are described in section 5.1. Initially, the per-DOM summary statistics TiT_i are fed into the transformer encoder as depicted in Fig. 4(b). Here, all the internal transformer layers are jointly denoted as “Transformer Encoder”, and transform tokens into transformed versions T~i T_i whose content is now an abstract encoding. In a final step tokens are aggregated by permutation invariant functions, which by default consist of forming the mean and standard deviation for every output feature over all the tokens and then mapping the result via an MLP to predict the normalizing-flow parameters. We also study alternative aggregation schemes with learnable aggregation tokens as described in section 5.1. Finally, the normalizing flow parameters define the target posterior pψ(θ|T1,…,TN)p_ψ(θ|T_1,…,T_N), where all parameters of the transformer and the final MLP are denoted as ψ. The posterior depends on the observed data summary statistics of each of the N DOMs in the event, T1,…,TNT_1,…,T_N (Fig. 4(c)), as well as the network parameters ψ and therefore structurally resembles an explicit posterior distribution with the mapping from data to parameters facilitated by the transformer. There are two probabilistic reasons in favor of transformers for variational inference using IceCube data and we call these “Bayesian inductive biases” in the following. 4.3 Bayesian inductive bias 1: Permutation invariance The first Bayesian inductive bias comes from the transformer permutation equivariance per token, which leads to invariance with an appropriate aggregation process. This is illustrated in Fig. 5(a). In general, Bayesian inference does not care in which order one looks at input parameters in the data distribution - they can be permuted, and the inference result would not change. This is even more explicit if the data is IID, which is the case for IceCube (see section 2). An encoding mechanism that respects this invariance to ordering, like a transformer with appropriate aggregation, therefore has some level of inbuilt advantages over other encoding mechanisms that do not respect it, like encoding schemes via a recurrent neural network (RNN) [27], for example. In early inference tests with normalizing flows we have used RNN encodings based on gated recurrent units [4], both using random DOM orderings at training and inference time or fixed orderings by mean time per DOM, and those consistently yielded worse results than GNNs or transformers. A GNN encoding based on k-nearest neighbors and edge convolutions [54] which has been used in IceCube before [44] can be made permutation invariant in the same manner as a transformer, but seems to lack generalizability to handle variable-sized inputs to the same extent. This will be discussed in the next section. 4.4 Bayesian inductive bias 2: dropping data factors Every neutrino interaction produces a different number of detected photons in the detector and the number of hit modules varies from event to event. The encoding scheme therefore has to handle variable-sized inputs. In the Bayesian picture (Fig. 5(b)), variable sized data inputs correspond to a varying number of data factors in the likelihood function. In IceCube, we sometimes have faulty DOMs so we additionally need to be able to remove inputs dynamically. A transformer can naturally handle variable-sized inputs, but so can a GNN or an RNN. However, when removing data factors in Bayes’ theorem we conjecture the transformer can handle those situations better than k-nearest-neighbor based GNNs because unlike the GNN no new spurious edge connections are created when removing data as illustrated in Fig. 5(b). This might partially explain the performance difference in section 6. We can additionally support prediction of the posterior for a varying number of inputs by probabilistically dropping inputs during training. We employ this as a form of data augmentation and probabilistic regularization during training (see section 5). Encoderp(θ|)=p(x1;θ)⋅p(x2;θ)⋅p(x3;θ)⋅…⋅p(θ)p()=p(x1;θ)⋅p(x3;θ)⋅p(x2;θ)⋅…⋅p(θ)p()p(θ|x)= [rgb]1,.5,0 [named]pgfstrokecolorrgb1,.5,0p(x_1;θ)· [rgb]0,0,1 [named]pgfstrokecolorrgb0,0,1p(x_2;θ)· [rgb]1,0,0 [named]pgfstrokecolorrgb1,0,0p(x_3;θ)·…· p(θ)p(x)= [rgb]1,.5,0 [named]pgfstrokecolorrgb1,.5,0p(x_1;θ)· [rgb]1,0,0 [named]pgfstrokecolorrgb1,0,0p(x_3;θ)· [rgb]0,0,1 [named]pgfstrokecolorrgb0,0,1p(x_2;θ)·…· p(θ)p(x) permute Bayes’ theorem with I.I.D. dataRNNData input:Posteriorp(θ|)p(θ|x)DOM 1→ 1x_1DOM 2→ 2x_2…Transformer (a) Illustration of permutation invariance of data arguments xix_i to the posterior. For independent data the data factors in Bayes’ theorem can be swapped explicitly. Permutation invariance is respected in a transformer encoding with appropriate aggregation, which gives inductive bias for posterior estimation. Encoding architectures that are not permutation invariant, like RNN encodings, do not have this inductive bias. GNNTransformer GNN TransformerData factors: x1x_1, x2x_2, x3x_3, x4x_4, …Data factors: x1x_1, x2x_2, x4x_4, …p(θ|)=p(x1;θ)⋅p(x2;θ)⋅p(x3;θ)⋅p(x4;θ)⋅…⋅p(θ)p()p(θ|x)= [rgb]1,.5,0 [named]pgfstrokecolorrgb1,.5,0p(x_1;θ)· [rgb]0,0,1 [named]pgfstrokecolorrgb0,0,1p(x_2;θ)· [rgb]1,0,0 [named]pgfstrokecolorrgb1,0,0p(x_3;θ)· p(x_4;θ)·…· p(θ)p(x)p(θ|)=p(x1;θ)⋅p(x2;θ)⋅p(x4;θ)⋅…⋅p(θ)p()p(θ|x)= [rgb]1,.5,0 [named]pgfstrokecolorrgb1,.5,0p(x_1;θ)· [rgb]0,0,1 [named]pgfstrokecolorrgb0,0,1p(x_2;θ)· p(x_4;θ)·…· p(θ)p(x) (b) Illustration of removing likelihood data factors in Bayes’ theorem. On the left side we use all data, including x3x_3, while on the right side we drop x3x_3 and predict a posterior without that data point. In contrast to the all-to-all connectivity of the transformer, the k-nearest neighbor edge-forming algorithm in the GNN creates new edges (pink edges) which more strongly affects the rest of the graph. Figure 5: Inductive biases of a transformer encoding for a) the data encoding in posterior estimation and b) the selective removal of data factors from the encoding. 5 Training details We focus on two types of events for training: charged-current νe _e-interactions (showers) and charged-current νμ _μ-interactions (tracks). An overview of the training datasets is given in table 1. The track dataset for training includes both throughgoing and starting tracks, so a single track model can handle both types. For showers, a νe _e charged-current model is also expected to work for neutral-current interactions or sub-PeV ντ _τ charged-current interactions. The lepton direction is taken as the target label, which above a few TeV is co-aligned with the neutrino direction due to kinematics. One can also train on the neutrino direction directly, but then one is more dependent on the spectral assumptions of the training dataset for low energies. The training datasets contain simulated neutrinos with energies ranging between 100100 GeV to 100100 PeV for both tracks and showers. For weighting purposes, we further calculate the deposited energy (see section 6 for details on the different weighting schemes, which are effectively treated as hyperparameters). For tracks, we define this as the sum of energy losses deposited within an “active” region which extends up to 350350m outside of the detector boundary. However, we only include tracks in training that cross within 5050m of the detector boundary at their point of closest approach. Similarly for showers, we require a shower to be within 5050m of the detector boundary. We found that not requiring such a cut would bias the training data too strongly towards events lying far outside the instrumented volume at high energies. Regular simulated neutrino events can contain coincident muons from air showers. Such events were removed to ensure pure signal events for training. The track training dataset contains all track event morphologies, including starting and throughgoing tracks, to obtain one unified track model after training. For testing, we use separate “throughgoing” and “starting” test datasets to later evaluate this trained model on those morphologies independently. In total, the shower training set comprises roughly 7 million events while the track training set comprises roughly 12 million events. For validation, we use 1000010000 events each with similar structure as the respective training datasets. All data is simulated using the newest ice model description FTP-v3 [51] and the default hole ice parametrization for FTP-v3 which corresponds to the central “Flasher unfolding” datapoint in the unified hole ice parametrization [45]. Training / validation Dataset Num. events Specifications showers ≈7≈ 7 million (train) 1000010000 (val) • only νe _e charged-current interactions • interaction max. 50m outside hexagon • no cosmic-ray air showers present • events pass filter for showers • 102GeV<Eν<108GeV10^2\ GeV<E_ν<10^8\ GeV tracks ≈12≈ 12 million (train) 1000010000 (val) • only νμ _μ charged-current interactions • all morphologies (starting/throughgoing) • dep. energy calculated within hexagon + 350m • no cosmic-ray air showers present • events pass filter for muons • 102GeV<Eν<108GeV10^2\ GeV<E_ν<10^8\ GeV Testing Dataset Num. events Specifications showers ≈10000≈ 10000 • well-contained showers, similar selection as [49] • contains both neutral and charged-current events • 104GeV<Eν<107GeV10^4\ GeV<E_ν<10^7\ GeV starting tracks ≈20000≈ 20000 • geometrically selected: starting within the hexagon • at least 300m of track length within the hexagon • similar to train dataset otherwise throughgoing tracks ≈20000≈ 20000 • geometrically selected as passing through the hexagon • at least 300m of track length within the hexagon • similar to train dataset otherwise Datasets for data and systematics checks (section 6.3) Dataset Specifications showers selection as developed for [47] with slightly more stringent final-level cuts tracks final-level selection as developed for [43] Table 1: Dataset specifications for the model training loop, split in training, validation and testing. Also shown are datasets used for data / Monte-Carlo comparisons and systematics checks. The “hexagon” refers to the quasi-hexagonal instrumented volume of IceCube. All Monte Carlo datasets use the FTP-v3 ice description [51]. A typical model with the described architecture from section 4 has a few million parameters, depending in detail on the chosen hyperparameters. Because a model training takes about one to two months on a RTX-3090 GPU machine in our setup, we decided to do manual hyperparameter variations based on outcomes of a set of parallel training runs. In total we trained about 60 models split over three consecutive training “periods”, where each period contains parallel runs. The hyperparameters from the best-performing run in each period were carried forward as “defaults” for the next period, with further variations introduced manually. For most trainings, we fix the normalizing flow to the one described in section 3.4 which has been found to be flexible enough to describe most relevant posterior shapes. The exception are runs that use a simple von-Mises Fisher distribution instead of a complex flow. We focus most of the hyperparameter variations on the exact architecture of the Transformer block as the encoding scheme has been found to be a major bottleneck in neural posterior estimation performance previously [14]. For run 1, we train some models on showers and some on tracks. For run 2, we focus solely on tracks. And for run 3, we again train on showers and tracks. The optimization process is performed via stochastic gradient descent (SGD) with ADAM [21] using default settings. We required flash attention v2 [6] to be able to train efficiently since we have a strongly varying token length between tens to thousands per batch item. For all training runs, we apply cosine-annealing [24] on the learning rate whose amplitude we reduce at fixed intervals of 200k steps. We found cosine-annealing to help with instabilities of previously used normalizing flows. The novel flow described in this paper likely does not require it anymore, but we retained it anyway. Additionally, some runs have a burn-in duration of 50k steps where the learning rate is linearly increased from zero to the initial learning rate, before starting with the scheduling. We anneal the learning rate towards a final learning rate ηf _f that is defined as ηf=2⋅BN _f=2· BN, with B the batch size and N the dataset size. This is motivated by the interpretation of SGD as stochastic variational inference over the weight posterior [25]. The batch size B is 80 through all trainings. We split each batch into 4 sections. In the first 20 items, we include all DOMs, in the second 20 items, we remove empty DOMs, in the third 20 items, we remove saturated DOMs, and in the last 20 items, we remove both empty and saturated DOMs. Additionally, with a 50%50\% chance we further remove a random amount of DOMs from the input for each batch item. This helps as a form of “probabilistic regularization” as the model learns to mimic the structure of Bayes’ theorem by training on all combinations of IID data factors, which is especially efficient for transformers via Bayesian inductive bias (see Fig. 5(b)). Dropping data randomly can also be seen as a form of data augmentation. Every time the same event is passed to the model during training, its input data to the model might look different. We perform stochastic weight averaging (SWA) [18] during the whole training run on the side to always have an averaged model available. For showers, we typically exit the run early due to overfitting at the highest energies (see a discussion in section 5.2 or Fig. 6), and the model used for testing is then the averaged model available at the best validation loss. For tracks, we also perform a few more 100k iterations of SWA at the end of training with the final learning rate ηf _f and the test model is then taken at the end of that. More details on how a model is selected for testing based on the validation loss are given in section 5.2. 5.1 Hyperparameters Our baseline transformer model in training period 1 uses 20 layers, pre-layer normalization, an internal embedding dimension of 9696 and a single head. The MLP part uses a hidden dimension of size 512. We aggregate at the end using the mean and diagonal variances of the final tokens and map the output to a “bottleneck” dimension of size 32. This representation is then mapped to the normalizing flow with a single MLP and hidden dimension 128. We vary this baseline multiple times which comprises training period 1. For the further periods we took the best models, varied them again, and in this way trained period 2 and 3. The hyperparameter extensions include a variation of the bottleneck dimension, an increase of the embedding dimension and the number of heads. We also checked standard sinusoidal positional encoding based on the (x,y,z)(x,y,z) DOM coordinates and relative value positional encoding [33] in several variations, the latter one having been used in the recent IceCube Kaggle challenge [2]. For relative positional encoding, we had to cap the maximum allowed number of input tokens to 200 due to memory constraints. We further explored using two simultaneous residual streams [56] and nonlinear QKV embedding, which potentially captures correlations between the Q,K and V tokens during the embedding stage. As an alternative to the standard aggregation scheme we test three different variants of class tokens, similar to their usage in vision transformers [7]. In order to test the influence of a fully-fledged normalizing flow compared to a simpler distribution, we also train standard vMF distributions with two different rotation parametrizations. Lastly, we compare against a non-transformer model based on GNNs with an architecture that has been described in [44]. More details on the hyperparameters are given in appendix D. 5.2 Validation curves (a) Best track model. (b) Best shower model Figure 6: Validation-loss curves for the best track and shower model. The x-axis shows optimization steps. Curves are shown for different energy ranges, along with the equally weighted average validation loss (dashed). The vertical dashed line marks the step at which the model is frozen for testing. The validation curves for tracks and showers look qualitatively slightly different as seen in Fig. 6 which shows the validation curves for the the best throughgoing track model and the best shower model. The validation is split up into curves for different deposited energy in log space. An average validation loss, which is calculated by weighting each energy-decade loss curve equally, is also shown. The model for further testing is selected based on the minimum of this average validation loss, as illustrated by the vertical dashed line. For tracks, we have more training data available at the highest energies, and there is minimal to no overfitting. The model for testing is typically obtained after SWA has run for an extended period late in the training run (see Fig. 6(a)). For showers, above a PeV there are only a few hundred thousand events in the training sample, while we have millions of events available around a TeV. During training, we train on all events equally weighted per deposited logarithmic energy, which means the high-energy events will be picked more often in batches and overfit earlier than at low-energy (Fig. 6(b)). For the same reason, the validation dataset also has less events at higher energies, and the loss for 100 TeV to 1 PeV (red curve) is actually below the loss for events between 1 TeV to 10 PeV (purple curve), while for enough statistics one should always expect an ordering based on energy since the associated posteriors become more compact. We accept this slight imprecision in validation and for the future plan to use a larger admixture at high energies. In general, there is always slight compromise for the tested shower models, where the low energy events could gain more from further training, while the high-energy events require an earlier cutoff to decrease the overfitting. This can be remedied for future studies with increasing the available training data at high energies and balancing out the statistics. 6 Results We test the shower models on a test dataset that has been originally introduced in [49] and consists of roughly 10000 well-contained showers between 10 TeV and 10 PeV. For tracks we have two different test datasets: roughly 20000 events sampled equally in log-energy that are starting in the detector and roughly 20000 tracks that are throughgoing. This separation was not made during training, where all track morphologies were used to train the models. All test datasets use the same ice description as in training, FTP-v3 [51]. More details are given in table 1. To give an impression of some of the predicted posteriors, example event views are shown in appendix A. 6.1 Test results (a) Shower test results (b) Throughgoing track test results Figure 7: Test results for different hyperparameters for showers and throughgoing tracks. Models are sorted by average total test loss (black). The best test loss per energy range is highlighted by a larger marker and a corresponding vertical cashed line to compare to other models. Options a), b) and c) (“Extra Info”) are described in the text in section 6. A detailed description of all parameters of each model is given in table 3 and 2. In Fig. 7 we summarize the test losses for showers and throughgoing tracks for all models. The results for starting tracks are summarized in Fig. 17 in appendix D. The models are sorted along the y-axis based on performance of the average total test loss. For starting tracks, that sorting uses the throughgoing track test results in order to see if the order changed dramatically, which it does not. Along the x-axis the test loss is split into different deposited energy regions. The best loss within a given energy loss is highlighted with a slightly larger marker. This is done to differentiate cases where the average model loss might be good for a specific model, but maybe another model has a better low energy performance. This is the case for model 9 and 11 for tracks, for example, as they were trained using a pure spectral weighting which over-emphasizes low energy events to the neglect of the PeV events and beyond. On the right of the plot we summarize relevant hyperparameter settings to give a quick overview in addition to the precise summary in tables 2 and 3. Important settings include the weighting scheme during training, the maximum number of DOMs given to the algorithm, and also the positional encoding for the transformer. One weighting scheme used is a weighting to a spectrum with spectral index of −1.8-1.8 in neutrino energy (we call this “si” weighting in Fig. 7), which was observed to approximately lead to a flat distribution in observed energy. In practice, however, it actually slightly overemphasized low-energy events in early training runs. Therefore, we introduced a second weighting that explicitly constructs weights based on a flat distribution in observed deposited energy (“flat”) - a scheme that turned out to produce an overemphasis on high-energy events while low-energy performance was lacking. Therefore, for training session 3, we also used weights that are a mixture of both weighting schemes (“si+flat”). This weighting turned out to give a good performance over all energies from 100100GeV to 100100 PeV. For relative value positional encoding , we had to limit the number of allowed input DOMs to be 200 - ordered by observed charge - due to memory constraints. The column “Extra Info” summarizes three specific hyperparameters that we found had the overall biggest influence on performance improvements. These are (a) nonlinear in-projection to the QKV tensor, (b) dual residual connections (“ResiDual”) [56] and (c) aggregating with an “improved class token”, a modified version of the class token in vision transformers [7] that interacts via cross-attention, as indicated in Fig. 7. A more detailed list of the used parameters and their description is found in table 2 and 3 in appendix D. For shower models, absolute position encoding (model 20) or no positional encoding (model 19) did not yield good results in training period 1 compared to relative value positional encoding with maximally 200 DOM input tokens. It seemed the information could not be adequately processed, although more training runs are necessary to settle this. However, switching to dual residual connections (model 5) changes this, even though a limitation of input tokens is then still favored (model 3). Adding nonlinear QKV projection in addition (model 1 and 2) gives another big boost which then allows the information of the full detector to be utilized without any input restriction. In the future it might be worthwhile to also try standard positional encodings again with these improvements. For track models, relative absolute or no positional encoding with the full DOM token input was always slightly better (model 18 or 20) compared to restriction in input tokens and relative value positional encoding (model 21). Again, using dual residual connections in combination with nonlinear QKV projection yields top performances (model 1 and 2). However, dual residual connections and a separate class aggregation gives nearly similar performance in this case (model 3). It might be interesting to combine all three of them in the future. In both cases of showers and tracks, a GNN encoding following the structure from [44] ranks last in model performance compared to all tested transformer architectures. 6.2 Angular resolution and coverage (a) Neutrino-induced showers (b) Throughgoing tracks (at least 300m distance within detector) (c) Starting tracks (at least 300m distance within detector) Figure 8: On the left, angular resolution (16%16\%, 50%50\%, 84%84\% quantiles) of the transformer normalizing flow (TNF) with and without saturated DOMs (TNF no sat.) compared to the state-of-the-art respective likelihood method based on B-splines, SplineMPEMax [32] for tracks and Taupede2024 [49] for showers. The improvement of median angular resolution of the new method versus the respective likelihood method is shown in the lower left figures. The right part shows expected versus observed coverage for all coverage probabilities. The corrected coverage for TNF assumes the von-mises-Fisher approximation for its calculation. In the following we calculate angular resolutions using the maximum of the posterior (MAP) estimate of the posterior estimate in order to compare to the respective maximum likelihood estimator. The transformer normalizing flow is abbreviated as “TNF” for simplicity. Figure 8(a) shows the median angular resolution for showers using TNF improves by up to a factor of two compared to the existing best B-spline based likelihood reconstruction Taupede2024 [49]. It should be noted that the Taupede2024 shower reconstruction by default excludes saturated DOMs and DOMs with a charge higher than 15 times the mean charge, while the TNF reconstruction has no exclusions. This has to be considered for the comparison mostly above 100 TeV, while below 100 TeV there are few saturated DOMs or high-charge outliers and the comparison can be interpreted as running on the exact same input. There is typically no uncertainty contour associated to the B-spline construction for showers, mostly because a profile-likelihood map is costly and can take several hours on a cluster (see Fig. 12) and an alternative via the Hessian has not been established in terms of numerical stability. In analyses applications [47], we typically use separate uncertainty predictors trained with high-level features for these reasons, but here we are interested in a self-consistent method comparison. That is why we do not include a coverage result for the B-spline applied to showers. The coverage for the transformer-based normalizing flow is calculated using the exact contours of the full-sky probability maps. It seems to perform well overall, but when split among energies (Fig. 9(a)) it becomes evident that it is good at low energies, but starts to undercover from a few 100 TeV upwards. This has to do with the lack of training statistics at high energies which shows up as earlier overfitting as discussed in section 5.2. Training on more simulations should alleviate this problem. For throughgoing tracks (Fig. 8(b)) the median angular resolution improves over the whole energy range compared to the standard method based on B-splines, SplineMPEMax [32], that has been in use for most track-based neutrino point source analyses in the past, e.g. in the detection of NGC 1068 [43]. An improvement or gain by a given factor means the angular resolution decreases by that factor and allows better localization of a neutrino source. In the relevant energy range from 1 to 100 TeV, the gain in angular resolution is about a factor of 1.11.1 to 1.31.3. This is the first time a machine-learning reconstruction yields consistently better results anywhere over the whole energy range since the introduction of SplineMPEMax in 2014. Additionally, at high energies there is no sign of a strong resolution floor. The standard B-spline reconstruction flattens out, since it does not fully take into account stochastic losses which dominate at PeV energies. The transformer-based normalizing flow, on the other hand, seems to perform well also in this regime. For starting tracks the performance gain in angular resolution is larger than for throughgoing tracks. The gain over the B-spline reconstruction reaches a factor of 2.52.5 at 100 TeV (Fig. 8(c)) and goes beyond a factor of 33 above a PeV. One reason is that the B-spline ansatz uses a chain of seeding reconstructions [40] which often gets stuck in local optima for starting tracks due to the extra hadronic shower from the neutrino interaction and shorter overall track lengths. The transformer-based normalizing flow seems to adapt to different morphologies automatically, as it was trained jointly on both throughgoing and starting tracks. In constrast to showers, for tracks we use a standard profile likelihood fit [28] to obtain an uncertainty estimator for SplineMPEMax. Typically, an energy-dependent correction is applied to widen the contours. Even after this correction, the coverage properties of the B-spline reconstruction are undercovering for the tails of the distribution, in particular for starting tracks. The TNF reconstruction for tracks is overcovering, which is conservative. We analyse the energy dependence of this behavior later in this section. It should be said that the track B-spline method [32] uses an outdated ice model for numerical reasons, and efforts are underway to improve it with a newer ice description in the near future. In all cases we observed that the inclusion of saturated DOMs as described in Fig. 4(b) helps at high energies for angular resolution. The inclusion of empty DOMs did not have a big influence, so we did not include it in the plots. For energy reconstruction, on the other hand, we expect the inclusion of empty DOMs to have an non-negligible contribution to the performance at all energies. Looking at the energy dependence of the track coverage for TNF (Fig. 9(b)), it has proper coverage for lower energies, but starts to have too large contours above a few 100 TeV. This effect is the opposite of the behavior for showers. The training statistics for tracks at PeV energies is a factor of a few larger than for showers, and this shows here. The over-coverage likely comes from the fact that the absolute resolution goes down to sub-degrees at high energies and several scales of resolution are combined during training. Events with small angular contours dominate the loss function, and start fluctuating more compared to lower energy events. We suspect further improvements with the exact parametrization of the flow, in particular the rotation parametrization, and a better tuned optimization routine might alleviate this issue in the future. It should be noted however, that over-coverage means the predicted uncertainty contours are too large, which can be acceptable and is in general conservative. (a) Neutrino-induced showers (b) Throughgoing tracks (at least 300m distance within detector) Figure 9: Coverage split up in different deposited energy bands. 6.3 Application on real data and systematics In order to test the robustness of the algorithm, we apply it on data and Monte Carlo datasets with various systematic uncertainties. For showers, we apply it on a small test dataset of 3 months livetime using a final level shower selection similar to the one used in [47] (see also table 1). Fig. 10 a) shows the shower data/MC comparison of the MAP of the zenith and azimuth posterior after matching the total Monte Carlo rate to the data rate. The Poisson uncertainty is indicated for visualization purposes by error bars on the data, while the uncertainty from finite weighted Monte Carlo statistics is indicated as colored bars for the respective Monte Carlo prediction using a Gamma approximation following [1]. The deviation in the lower part of the figure utilizes the combined uncertainty via an appropriate Gamma-Poisson Mixture. Cosmic-ray interactions in the atmosphere produce “atmospheric” neutrinos and muons which show up as irreducible background events in this selection. Muons are irreducible here because they enter the detector between strings or skim a corner of the detector so they appear as showers. The atmospheric neutrino fluxes assume a primary cosmic-ray flux model from [12] and atmospheric interactions following the conventional model in [9]. Also shown are “astrophysical” diffuse neutrino flux predictions assuming the measurement results from [52]. A few outliers are visible which is likely due to the low Monte Carlo statistics of the cosmic-ray muon background and unmodeled ice systematics. Overall the data/MC agreement is similar to expectations from comparisons with equivalent B-spline reconstructions. It is noteworthy that most of the remaining events from downgoing muons are correctly reconstructed as downgoing, even though the network was only trained on electron neutrino charged-current interactions and has never seen muons during training. Furthermore, Fig. 10 b)-e) show various systematics checks, where the angular resolution of the MAP is compared to the no-systematics baseline. The “ice systematics” curve corresponds to an average of different effects which comprise photon absorption, photon scattering, effective light yield and angular PMT acceptance related to the local ice structure. Overall the impact on the angular resolution by ice-related systematics is a little larger than detector related systematics, which is known behavior in existing shower reconstructions. The removal of the 5 or 10 highest charge DOMs per event affects mostly low-energy events, which can be explained by those events already having a low number of DOMs to begin with, and the removal of a fixed amount of high charge DOMs has a larger relative impact in this case. (a) (b) (c) (d) (e) Figure 10: Data / Monte Carlo comparison of the zenith and azimuth of the maximum of the posterior (MAP) (a) and effects of various systematics on angular resolution (b-e) for a slightly modified selection of the one used in [47]. The systematics checks are visualized as 16%16\% (dotted), 50%50\% (solid) and 84%84\% (dashed) quantiles to indicate the distributions. They include (b) ice model systematics, (c) the removal of “x” highest charge DOMs, (d) random variation (in terms of standard deviation) in (x,y,z)(x,y,z) coordinates of individual strings and (e) random variation (in terms of standard deviation) of absolute time offset in each DOM. For tracks, we use 2.5 months of livetime for data using the final sample selection developed in [43]. Fig. 11 shows similar data/MC and systematics comparisons as for showers, again after matching the overall Monte Carlo rate to the data rate. The Data/MC comparison looks good, similar to existing likelihood reconstructions. Overall, the relative ice-systematic effect is slightly smaller than for showers, while geometry effects play a bigger role. This makes sense, since track events on average contain more unscattered Cherenkov light which can reduce the impact of ice effects. On the other hand, the angular resolution of a few tenths of a degree is precise enough such that position and timing of the modules becomes more relevant compared to showers whose resolution is about an order of magnitude worse. (a) (b) (c) (d) (e) Figure 11: Data / Monte Carlo comparison of the zenith and azimuth of the maximum of the posterior (MAP) (a) and effects of various systematics on angular resolution (b-e) for the final track selection used in [43]. The systematics checks are visualized as 16%16\% (dotted), 50%50\% (solid) and 84%84\% (dashed) quantiles to indicate the distributions. They include (b) ice model systematics, (c) the removal of “x” highest charge DOMs, (d) random variation (in terms of standard deviation) in (x,y,z)(x,y,z) coordinates of individual strings and (e) random variation (in terms of standard deviation) of absolute time offset in each DOM. For specific analyses, it can make sense to include systematics during training. This can be achieved either by an explicit inclusion in the model, or by training on an ensemble of different systematics realizations and thereby effectively learn the marginalized posterior as demonstrated in [14]. We leave this for future work. 6.4 Timing Figure 12 shows the runtimes of TNF compared to the standard B-spline construction for both showers (Fig. 12(a)) and tracks (Fig. 12(b)). We split runtime calculations in three parts. The first part is moment prediction, which involves sampling from the PDF (by default we sample 10000 times) and calculating the first moments. We obtain the mean and the kappa parameter of the vMF approximation that closest matches the samples via a likelihood fit. The corresponding step for the B-spline is to obtain the full uncertainty contour. The second part is uncertainty evaluation, which just performs a single PDF evaluation — a quantity that can be important for re-evaluation of uncertainties during an analysis, for example a point source likelihood analysis. The third part is a skymap scan, which involves the skymap creation strategy described in section 3.4. It starts with 10000 samples, which then guides the creation of the HEALPIX evaluation grid which is on the order of 5000 pixels. For the B-spline skymap-scan, the whole sky is typically scanned with a refinement strategy and parallelized on a cluster. All the checks for TNF are performed for a pure CPU application with up to 4 cores and additionally with a Geforce GTX 3090 for the GPU check. Additionally we separate into total time and time without pre-processing. The preprocessing time is the time it takes to convert IceCube specific data to a format that feeds into the transformer or GNN encoder, and then also pass it through that encoding part, but before we activate the normalizing flow. For the respective B-spline reconstruction, pre-processing involves all reconstructions in the seeding chain that come before it. For moment prediction and uncertainty evaluation, we calculate results in batches and quote the time per batch item. For the skymap scan, we perform it for a single batch item at a time. For showers, the moment prediction of TNF is faster by at least a factor 10 than the likelihood counterpart. For uncertainty prediction the go-to algorithm for the B-spline is an all-sky profile-likelihood scan, so we can directly compare the skymap-scan running times, since they are roughly equal for the TNF in both cases as they are dominated by preprocessing time. Here we see TNF is faster by several orders of magnitude compared to the B-spline reconstruction, even if we fully utilize a cluster of several 100 compute nodes. For tracks the time spent in moment prediction is roughly on par between TNF and the B-spline based method. However, once we go to uncertainty evaluation, which in the B-spline case relies on a profile-likelihood scan around the best fit point [28], we have a speed up by more than an order of magnitude. For the skymap scan, the speed up is again several orders of magnitude as for showers. We can see that in all cases the proportion of time spent in preprocessing rises with energy. This could potentially be optimized in the future as it currently relies in parts on python implementations that could be ported to C/C++. (a) Timing comparison showers (b) Timing comparison tracks Figure 12: Runtimes of TNF for CPU and GPU in comparison to the B-spline method which runs on CPU. For TNF we allow for 4 CPU cores (CPU) with access to an additional RTX3090 (GPU). The total time is shown in solid, the time without preprocessing and encoding in dashed. 7 Discussion and Outlook In this paper, we combined a spherical normalizing flow with a transformer encoding to learn a model of the posterior distribution of the direction of leptons induced by high energy neutrinos in IceCube - an instance of amortized neural posterior estimation and ELBO-free variational inference. The posterior model takes as input the summary statistics of the raw photon data in each DOM, passes it as per-DOM data tokens through a transformer encoding, aggregates the result in a single vector, and in a final step maps that vector to a normalizing-flow posterior over the direction. To this end, we developed a novel spherical normalizing flow that combines smooth rational quadratic splines, scale transformations and rotations to learn a flexible conditional posterior that can span several orders of magnitude in extent for different events. For the rational quadratic splines in particular, the C2C^2 smoothness constraint together with a simplistic yet flexible spline-based conditioning were crucial for the stability of the learning process. We argue that the transformer encoding provides Bayesian inductive bias in two ways: a) it encodes the permutation invariance of the labeling of the input dimensions of the data PDF in Bayes’ theorem and b) it provides generalization capabilities when IID data factors (DOMs) are dropped and information is removed. After hyperparameter optimization of the transformer architecture in several iterative training sessions, we find the leading model for charged-current electron neutrinos (showers) and charged-current muon neutrinos (tracks) to be architecturally similar up to positional encoding. In both cases, the leading model mostly benefits from nonlinear QKV projection in each attention layer and dual residual pipelining (“ReSiDual”) compared to our vanilla transformer baseline. For tracks, a specific weighting scheme was additionally important to balance out low- and high energy tracks during the learning process. We compare the leading models to state-of-the-art B-spline-based likelihood reconstructions. The angular resolution is better than the B-spline-based likelihood reconstructions over the whole tested energy ranges. The gain is between a factor of 1.11.1 for 1 TeV throughgoing tracks up to a factor of 1.51.5 or 2.52.5 for 100 TeV showers or starting tracks, respectively. Processing times are also favorable, in particular for all-sky scans which are sent out in realtime alerts where the processing time can be pushed down from hours to seconds. This is achieved by synergies of the normalizing flow properties together with irregular MOC HEALPIX maps. The coverage for the trained posteriors is typically also better than the B-spline based counterparts, but not perfect. There is slight undercoverage at higher energies for showers, which stems from a lack of training data, and overcoverage at high energies for tracks, which stems from the interaction of the final large learning rate together with the specific flow parameterization. This prevents the neural-network prediction for high energy events with very small contour sizes to settle into local weight minima and effectively broadens their predicted contours compared to what they could be otherwise. In the near future, we expect these remaining issues to disappear with more training data and some finetuning of either the late-stage learning-rate scheduling or the precise flow parametrization. Furthermore, dedicated training on systematics datasets can incorporate complex systematic uncertainties and should make the application to real data robust. As the leading model for showers and tracks shares the same architecture - up to positional encoding - it is conceivable we will soon converge on hyperparameters that also agree in totality. It will then be possible to do a joint training of all neutrino interactions, including tau neutrinos which sit morphologically between showers and tracks. This would harness extra training statistics while being confident that neither individual class loses performance due to the model architecture. Looking further into the future, the algorithm is suited to face the challenges of heterogeneous detectors consisting of different types of optical modules like the IceCube Upgrade [50] and planned IceCube-Gen2 [42] with minimal changes. The IceCube collaboration acknowledges the significant contributions to this manuscript from Thorsten Glüsenkamp. The authors gratefully acknowledge the support from the following agencies and institutions: USA – U.S. National Science Foundation-Office of Polar Programs, U.S. National Science Foundation-Physics Division, U.S. National Science Foundation-EPSCoR, U.S. National Science Foundation-Office of Advanced Cyberinfrastructure, Wisconsin Alumni Research Foundation, Center for High Throughput Computing (CHTC) at the University of Wisconsin–Madison, Open Science Grid (OSG), Partnership to Advance Throughput Computing (PATh), Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS), Frontera and Ranch computing project at the Texas Advanced Computing Center, U.S. Department of Energy-National Energy Research Scientific Computing Center, Particle astrophysics research computing center at the University of Maryland, Institute for Cyber-Enabled Research at Michigan State University, Astroparticle physics computational facility at Marquette University, NVIDIA Corporation, and Google Cloud Platform; Belgium – Funds for Scientific Research (FRS-FNRS and FWO), FWO Odysseus and Big Science programmes, and Belgian Federal Science Policy Office (Belspo); Germany – Bundesministerium für Forschung, Technologie und Raumfahrt (BMFTR), Deutsche Forschungsgemeinschaft (DFG), Helmholtz Alliance for Astroparticle Physics (HAP), Initiative and Networking Fund of the Helmholtz Association, Deutsches Elektronen Synchrotron (DESY), and High Performance Computing cluster of the RWTH Aachen; Sweden – Swedish Research Council, Swedish Polar Research Secretariat, Swedish National Infrastructure for Computing (SNIC), and Knut and Alice Wallenberg Foundation; European Union – EGI Advanced Computing for research; Australia – Australian Research Council; Canada – Natural Sciences and Engineering Research Council of Canada, Calcul Québec, Compute Ontario, Canada Foundation for Innovation, WestGrid, and Digital Research Alliance of Canada; Denmark – Villum Fonden, Carlsberg Foundation, and European Commission; New Zealand – Marsden Fund; Japan – Japan Society for Promotion of Science (JSPS) and Institute for Global Prominent Research (IGPR) of Chiba University; Korea – National Research Foundation of Korea (NRF); Switzerland – Swiss National Science Foundation (SNSF). Appendix A Example events At low energies, showers are often single-string dominated, which introduces azimuthal degeneracy into the posterior (Fig. 13 a). At higher energies, this degeneracy is typically broken and posteriors are more Gaussian (Fig. 13 b)+c)). Most track events are visible in multiple strings that are not lying within a plane and have fairly Gaussian posteriors, as in the high-energy example shown in Fig. 14 a). However, due to the geometry there are outliers that are co-aligned with symmetry axes. Fig. 14 b) shows an event that directly passes along string lines that all lie within the same plane, which leads to a slight bimodal degerenacy in azimuth. Fig. 14 c) shows a low-energy contained muon that passes upwards close to a single string whose posterior is ring-like. (a) A low-energy shower that is mostly visible in a single string. (b) A shower which is visible in several neighboring strings. (c) A multi-PeV shower in the clear ice. Figure 13: Contours (left) and event view (right) for example shower events. The contours are indicated in solid for 68%68\% (black) and 95%95\% (white) contained probability mass. The corresponding contours of the von-Mises Fisher approximation, i.e. the second moment, are shown with dashed lines. The colors in the event view indicate time, with red being early and blue being late. (a) Bright horizontal track. (b) Muon track along string lines. (c) Low energy contained upgoing track along a single string line. Figure 14: Contours (left) and event view (right) for example track events. The contours are indicated in solid for 68%68\% (black) and 95%95\% (white) contained probability mass. The corresponding contours of the von-Mises Fisher approximation, i.e. the second moment, are shown with dashed lines. The colors in the event view indicate time, with red being early and blue being late. Appendix B Smooth rational-quadratic spline derivation (a) Unconstrained rational quadratic splines. (b) Smooth splines with two knot intervals on cylinder height. (c) Smoth splines with three knot intervals on cylinder height. (d) Smooth periodic splines with two knot intervals. One interval wraps the boundary. Figure 15: Depictions of rational quadratic splines with various constraints. Knots are shown as circles. In general the upper plot shows the inverse spline function and the corresponding normalizing flow PDF that follows from eq. 8 is shown in the lower part. The exception is a) which also depicts the forward transformation. All colored functions denote a real-1-parameter spline parametrization, rqsI(z)arqs_I(z)_a for the cylinder height z and rqsA(ϕ)arqs_A(φ)_a for the periodic angle ϕφ, where the parameter a is a left-over parametric degree of freedom that varies from negative values (i.e. blue) to positive ones (i.e. purple). Setting a=0a=0 equals an identity mapping (green). The forward function f(z)f(z) represented as a monotonic rational-quadratic spline is defined per knot segment k between knot points k and k+1k+1 and given by [8] f[ξ(z)]=yk+(yk+1−yk)[sk⋅ξ2+δk⋅ξ(1−ξ)]sk+[δk+1+δk−2sk]⋅ξ(1−ξ),f[ξ(z)]=y_k+ (y_k+1-y_k)[s_k·ξ^2+ _k·ξ(1-ξ)]s_k+[ _k+1+ _k-2s_k]·ξ(1-ξ), (19) where xkx_k and xk+1x_k+1 are the knot x-positions, yky_k and yk+1y_k+1 the knot y-positions, δk _k and δk+1 _k+1 the derivatives at the knots and ξ(z)=z−xkxk+1−xk∈[0,1]ξ(z)= z-x_kx_k+1-x_k∈[0,1] a relative coordinate between the knots. We also define the width wk=xk+1−xkw_k=x_k+1-x_k and height hk=yk+1−ykh_k=y_k+1-y_k of segment k and sk=hkwks_k= h_kw_k the linear slope between the knots in the segment. If the derivatives δk _k at the knot points are positive the resulting function is monotonically increasing [16]. The exact derivative at any point between the knots is given by [8] dfdz[ξ(z)]=sk2[δk+1ξ2+2sk⋅ξ(1−ξ)+δk(1−ξ)2][sk+[δk+1+δk−2sk]⋅ξ(1−ξ)]2 dfdz[ξ(z)]= s_k^2[ _k+1ξ^2+2s_k·ξ(1-ξ)+ _k(1-ξ)^2][s_k+[ _k+1+ _k-2s_k]·ξ(1-ξ)]^2 (20) and is by construction positive due to the underlying monotonicity. Additionally, the exact analytic inverses of both these functions exist [8] but we omit to write them down here for brevity. Using these definitions and multiple such segments we can define a flexible normalizing flow on an interval using eq. 8. In this default parametrization, the derivatives at the knot points δk _k are fitted and can be small or large independent of the knot positions, which in general results in the function only being C1C^1-smooth. The jumps in the second derivative lead to sharp features in the PDF via eq. 8 (see Fig. 15(a)). In the following, we impose that the second derivative at the knot points should be equal between two consecutive knot segments which leads to C2C^2 smoothness. Forming the second derivative of segment k and evaluating it at its upper knot point k+1k+1 yields dfkd2z[z=xk+1]=2[δk+1(δk+δk+1)−sk(δk+1+sk)]skwk. df_kd^2z[z=x_k+1]= 2 [ _k+1( _k+ _k+1)-s_k( _k+1+s_k) ]s_kw_k. (21) Forming the second derivative of segment k+1k+1 and evaluating it at its lower knot point k+1k+1 yields dfk+1d2z[z=xk+1]=2sk+1(δk+1+sk+1)sk+1wk+1−2δk+1(δk+1+δk+2)sk+1wk+1. df_k+1d^2z[z=x_k+1]= 2s_k+1( _k+1+s_k+1)s_k+1w_k+1- 2 _k+1( _k+1+ _k+2)s_k+1w_k+1. (22) The term wkw_k denotes the width of segment k, wk=xk+1−xkw_k=x_k+1-x_k. The two equations can be set equal to each other to impose continuity on the second derivative which constrains the inner first derivative parameter δk+1 _k+1, i.e. dfkd2z[z=xk+1]=dfk+1d2z[z=xk+1]. df_kd^2z[z=x_k+1]= df_k+1d^2z[z=x_k+1]. (23) For N segments, we have N−1N-1 such constraints which are quadratic in the respective inner derivatives, and the equations are coupled since the derivative δk+2 _k+2 or δk _k in the constraint 23 also appear in constraints involving neighboring segments k+2k+2 or k−1k-1. For N such segments chained after each other, we have in total 4N4N knot position values xk,xk+1,yk,yk+1x_k,x_k+1,y_k,y_k+1 and 2N2N derivatives at the segment boundaries. We also have 2(N−1)2(N-1) constraints for the segments to be continuous in function value (for both x and y coordinates to agree) and another N−1N-1 constraints to be continuous to 1st order. This leads to 6N−3(N−1)=3N+36N-3(N-1)=3N+3 free parameters in the default parametrization. If we fix the outer knot positions of the left-most and right-most segment, which on the cylinder and the circle has to be done, we can subtract 44 again to obtain 3N−13N-1 free parameters. We have not imposed the 2nd2nd-order constraint (eq. 23) yet. In the following we differentiate two cases: application on the interval z∈[−1,1]z∈[-1,1] and on the circle ϕφ. B.1 Interval on cylinder height z∈[−1,1]z∈[-1,1] For the normalizing flow we describe in section 3.4 the sphere is deformed to the cylinder, and for non-singular behavior the spline flow on the cylinder height z should be equal to an identity mapping at z=−1z=-1 and z=1z=1, which correspond to the poles in the spherical representation. This can be achieved by enforcing the two outer-most knot derivatives to be equal to 1. If we now additionally apply a 2nd-order smoothness constraint using eqs. 23 on all inner N−1N-1 inner knots connecting segments, we have overall pfree=3N−1−(N−1)−2=2N−2 p_free=3N-1-(N-1)-2=2N-2 (24) free parameters for the total spline flow. Since the constraints in eq. 23 lead to N−1N-1 coupled quadratic equations in N−1N-1 variables, they are not trivial to be solved analytically in all generality. In the following we solve them for N=2N=2 and for N=3N=3 with an extra symmetry constraint. Case N=2N=2: For 2 segments we have a single quadratic equation in the intermediate derivative δ1 _1. The single usable solution is given by δ1=−p2+(p24−q) _1=- p2+ ( p^24-q ) (25) with p=h0⋅(s1−1)+h1⋅(s0−1)h0+h1, p= h_0·(s_1-1)+h_1·(s_0-1)h_0+h_1, (26) q=−h0h1(h0(h0+h1)w12+h1(h0+h1)w02). q=-h_0h_1 ( h_0(h_0+h_1)w_1^2+ h_1(h_0+h_1)w_0^2 ). (27) Since q is negative, the second solution for δ1 _1 is always negative which violates the increasing monotonicity of the resulting spline function and it can therefore be omitted. We have 2 free parameters left to define the shape (eq. 24), which we choose to be the intermediate knot x and y values, x1x_1 and y1y_1. Since the lower and upper knot positions are fixed to be x0=−1x_0=-1 and x2=1x_2=1 on the cylinder height, they are contained in the width and height parameters. The width parameters, for example, are defined as w0=x1−x0=x1−(−1)=x1+1w_0=x_1-x_0=x_1-(-1)=x_1+1 and w1=x2−x1=1−x1w_1=x_2-x_1=1-x_1. For numerical stability, we found it useful to use a 1-parameter parametrization, where w0=w1=1w_0=w_1=1 and the height is parametrized by a real parameter a in log space such that h0=h1=1h_0=h_1=1 if a=0a=0. This allows the transformation to have a smooth 1-parameter set of curves that defines the normalizing flow (see Fig. 15(b)) and results in an identity for a=0a=0. Case N=3N=3, symmetric: For N=3N=3 segments there are four free parameters to fit according to eq. 24. In addition, from the two 2nd-order constraints we have two quadratic coupled equations in 2 variables, the two intermediate knot derivatives δ1 _1 and δ2 _2. In general, those two coupled equations are hard to solve analytically. However, we can impose two extra symmetry conditions that the first and third segment dimensions are equal, i.e. w0=w2w_0=w_2 and h0=h2h_0=h_2. In this case, there are four solutions for the two derivatives. Only two of those solutions are positive for both derivatives at the same time, and only one of those two solutions has the same result for the two derivatives, which is what we expect from symmetry considerations. This unique solution for the two coupled quadratic equations turns out to be δ1=δ2=−p2+(p24−q) _1= _2=- p2+ ( p^24-q ) (28) with p p = = h1(w0w1−h0(w0+w1))(2h0+h1)w0w1, h_1(w_0w_1-h_0(w_0+w_1))(2h_0+h_1)w_0w_1, (29) q q = = −h0h1h1w02+h0w12(2h0+h1)w02w12. -h_0h_1 h_1w_0^2+h_0w_1^2(2h_0+h_1)w_0^2w_1^2. (30) In addition, the symmetry constraint eliminates two from the four overall free parameters, and we end up with two free parameters to determine the shape of the flow function, which can for example be w0w_0 and h0h_0. Again, it is useful to fix the width w0w_0, this time to a third of the overall interval via w0=2.0/3.0w_0=2.0/3.0, and model the height h0h_0 as a relative real 1-parameter curve a in log space such that a=0a=0 leads to the identity mapping, as depicted in Fig. 15(c). In contrast to the 2-interval solution (Fig. 15(b)) such transformations lead to symmetric PDFs on z. B.2 Interval on cylinder angle ϕ∈[0,2π]φ∈[0,2π] The cylinder angle ϕφ is periodic. Instead of fixing the two boundary derivatives to 11 as was done for the height z, there is only the single constraint that the derivatives have to agree. This results in one more free parameter. On the other hand, also the second derivative must agree on the boundaries, a constraint that did not exist for the height, which results in one less free parameter. In total the number of free parameters therefore is similar to the situation for cylinder height z and equals also 2N−22N-2 with slightly different constraint equations. For N=2N=2, which corresponds to 3 knot positions, we can solve for the derivative at the inner knot δ1 _1 and one of the outer knots, i.e. δ0 _0. Using the first constraint the two outer knot derivatives are the same, so δ0=δ2 _0= _2. Using the two second derivative constraints at the inner and one outer knot we obtain two solutions, of which only one is always positive. The two derivative solutions for δ0 _0 and δ1 _1 are always equal for the unique positive solution, similar to the situation with three segments and 4 knots on the cylinder height, but with slightly different dependency on the widths and heights of the segments. Since δ0 _0 = δ2 _2, the result for all derivatives is δ0=δ1=δ2≡δ=−p2+(p24−q) _0= _1= _2≡δ=- p2+ ( p^24-q ) (31) with p=−h0h1(w0+w1)2(h0+h1)w0w1, p=- h_0h_1(w_0+w_1)2(h_0+h_1)w_0w_1, (32) q=−h0h1(h1w02+h0w12)2(h0+h1)w02w12. q=- h_0h_1(h_1w_0^2+h_0w_1^2)2(h_0+h_1)w_0^2w_1^2. (33) We fix w0=w1=πw_0=w_1=π and define the height again as a relative 1-parameter curve in log-space parametrized by a, where a=0a=0 leads to the identity. The solution turns out not to be symmetric over the interval [0,2π][0,2π]. We can however introduce an extra adaptive “shift” s to effectively turn the transformation into a symmetric one where the peak does not move. The total transformation is rqsA(ϕ)=fbase(ϕ)+smod2π,rqs_A(φ)=\f_base(φ)+s\ 2π, (34) with shift s, s=2π−h0−h1(12w0−π)(h0(12w0−π)−δw0(12w0+w1−π))(h0w12+2(12w0−π)(h0−δw0)(12w0+w1−π)), s=2π-h_0- h_1( 12w_0-π)(h_0( 12w_0-π)-δ w_0( 12w_0+w_1-π))(h_0w_1^2+2( 12w_0-π)(h_0-δ w_0)( 12w_0+w_1-π)), (35) where fbase(ϕ)f_base(φ) is the standard smooth rational quadratic spline without shift. Note that we have not added an extra knot, so there are in total still three knots or two intervals at use. One of the intervals however wraps around the 0/2π0/2π boundary with the included shift. We use this form in all experiments, as it is more stable for optimization due the independence of the peak position with the parametrized degree of freedom. The final transformation for the angle including this shift is illustrated in Fig. 15(d), which shows the 1-real-parameter parametrization rqsA(ϕ)arqs_A(φ)_a of this transformation in dependence of the free degree of freedom a. In [30] the azimuthal transformation is defined in dependence on the cylinder height z, by conditioning on z using an MLP. We found this can lead to numerical issues in the polar regions. In order to mitigate this, but still have some form of conditioning scheme implemented, we use a simple spline-based non-MLP conditioning that blends towards an identity mapping at the polar regions, as depicted in Fig. 16. Figure 16: Conditioning function a=fa∗(z)a=f_a^*(z) on the cylinder height z to produce a regularized parameter a to be used for the azimuthal angle flow function, rqsA(ϕ)arqs_A(φ)_a, depicted in Fig. 15(d). At the poles (z=−1z=-1 and z=1z=1) the output is always a=fa∗(z=−1/z=1)=0a=f_a^*(z=-1/z=1)=0 which enforces an identity mapping as rqsA(ϕ|z)a=0=ϕrqs_A(φ|z)_a=0=φ. The horizontal dashed lines correspond to specific values of a∗a^*. The function is just parametrized by single parameter a∗a^*, its overall normalization, and towards the poles, i.e. in the cylinder representation at z=−1z=-1 and z=1z=1, always tends towards 0 by construction. The parameter a∗a^* acts as a surrogate parameter for the actual parameter a=fa∗(z)a=f_a^*(z) which is then plugged into the parametrization in Fig. 15(d). The end result is the angular flow transformation rqsA(ϕ|z)a∗rqs_A(φ|z)_a^* that depends on the cylinder height z via the single parameter a∗a^*. At the poles the mapping a=fa∗(z=−1/z=1)=0a=f_a^*(z=-1/z=1)=0 and the corresponding polar angle transformation rqsA(ϕ)a=0=ϕrqs_A(φ)_a=0=φ is the identity (green curve in Fig. 15(d)). This scheme is quite effective in turning an extra ϕφ transformation on or off if desired and numerically stable since it smoothly blends to an identity at the poles irrespective of a∗a^*. The explicit form of fa∗(z)f_a^*(z) is a fifth-order polynomial spline in two intervals which was derived enforcing the necessary boundary conditions, fa∗(z)=a∗(6z5+15z4+10z3+1),z≤0,a∗(−6z5+15z4−10z3+1),z>0,f_a^*(z)= \ array[]la^*\,(6z^5+15z^4+10z^3+1),&z≤ 0,\\ a^*\,(-6z^5+15z^4-10z^3+1),&z>0~, array . i.e. it contains appropriate smoothness constraints at z=0z=0 and fulfills fa∗(z=−1/z=1)=0f_a^*(z=-1/z=1)=0. Appendix C Rotation Parametrizations The following are parametrizations of the rotation matrix RR discussed in section 3.4. The baseline rotation parametrization is based on householder reflections [17]. An individual householder reflection in 3-dimensional space can be defined by a 3-dimensional vector v as Rhh,i=−2v×vT/|v|2.R_h,i=1-2v× v^T/|v|^2. (36) For the full rotation we chain three such transformations as Rhh=Rhh,1∘Rhh,2∘Rhh,3,R_h=R_h,1 _h,2 _h,3, (37) which has in total 99 real parameters and technically ends up in an improper rotation matrix with determinant −1-1, since each householder reflection has negative parity. For pure von-Mises-Fisher distributions we use another representation which reads RF=(1−μ12/(1+μ3)−μ1μ2/(1+μ3)μ1−μ1μ2/(1+μ3)1−μ22/(1+μ3)μ2−μ1−μ2μ3)R_F= pmatrix1- _1^2/(1+ _3)&- _1 _2/(1+ _3)& _1\\ - _1 _2/(1+ _3)&1- _2^2/(1+ _3)& _2\\ - _1&- _2& _3 pmatrix (38) and is parametrized by a normalized mean vector μ→ μ. As briefly mentioned in section 3.4, we can describe a von-Mises-Fisher distribution via the change-of-variables formula which makes it a normalizing flow. In particular, the determinant of the Jacobian of the inverse of the block transformation in eq. 14 where the smooth spline transformations are identities, i.e. rqsϕ−1(x)=xrqs_φ^-1(x)=x and rqsI(x)=xrqs_I(x)=x, yields det(F−1) (J_F^-1 ) = = (dzF−1)[RF−1⋅x→]z ( ddzF^-1 )[R_F^-1· x]_z (39) = = κsinh(κ)eκ⋅[RF−1⋅x→]z κsinh(κ)e^κ·[R_F^-1· x]_z (40) = = κsinh(κ)eκ⋅μ→⋅x→. κsinh(κ)e^κ· μ· x. (41) since the Jacobian determinant of the rotation matrix and the cylinder transformations from eq. 14 are all equal to one, and the only term remaining after some algebra is dzF−1(z) ddzF^-1(z), evaluated at the z-component of RF−1⋅x→R_F^-1· x. Using eq. 8 we can combine this density update with a flat base distribution, p0=1/4πp_0=1/4π, to obtain the known von-Mises-Fisher density. We use this rotation parametrization only to describe the standard vMF density. In practice we found it to be unstable when used in longer normalizing-flow chains. Appendix D Additional test results and model hyperparameter settings This section lists the test results for starting tracks (Fig. 17). It also contains a summary of the hyperparameter settings used for the different models in track training (table 2) and shower training (table 3). The specific settings for each model are indexed by column, where the setting description to each column is given in table 4. Other abbreviations used in describing the parameter settings are given in table 5. Figure 17: Test results for different hyperparameters for starting tracks. Models are sorted by average total test loss from throughgoing tracks as shown in Fig. 7(b) which defines the “model ID”. The best test loss per energy range is highlighted by a larger marker and a corresponding vertical cashed line to compare to other models. Options a), b) and c) (“Extra Info”) are described in the text in section 6. A detailed description of all parameters of each model is given in table 2. Track training Numeric options (see tables 4 and 5 for definitions) Model ID 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 1 both all yes yes md no sin none no 96 1 no h no no no 2 both all yes yes md yes sin none no 96 1 no h no no no 3 both all no yes c3 no sin none no 96 1 no h no no no 4 both all no yes md no sin none no 96 1 no h no no no 5 both all no yes c1 no sin none no 96 1 no h no no no 6 both all no yes md no none none no 96 1 no h no no no 7 both all no yes c2 no sin none no 96 1 no h no no no 8 both all no yes md no sin none no 96 1 no h no no no 9 si all no yes md no sin none no 96 1 no h no no no 10 flat all no yes md no sin none yes 96 1 no h no no no 11 si all no no md no sin none yes 96 1 no h no no no 12 flat all no yes md no sin none yes 96 1 yes h no no no 13 flat all no no md no sin none yes 192 4 no h no no no 14 flat all yes no md no sin none yes 96 1 no h no no no 15 flat all no no md no sin none yes 96 1 no h no no no 16 flat all no no md no sin none no 96 1 no h no no no 17 flat all no no md no sin none yes 192 1 no h no no no 18 flat all no no md no sin none yes 96 1 no h no no no 19 flat all no no md no sin none yes 192 2 no h no no no 20 flat all no no md no none none yes 96 1 no h no no no 21 flat 200 no no md no none rel1rel_1 yes 96 1 no h no no no 22 flat all no no md no sin none yes 96 1 no h no no no 23 flat 200 no no md no none rel4rel_4 yes 96 1 no h no no no 24 flat 200 no no md no none none yes 96 1 no h no no no 25 flat all no no md no sin none yes 96 1 no h yes no no 26 flat 200 no no md no none rel6rel_6 yes 96 1 no h no no no 27 flat all no no md yes sin none yes 96 1 no h no no no 28 flat all no no md no sin none yes 96 1 no h no no no 29 flat all no no md no sin none yes 96 1 no h no yes no 30 flat all no no md no sin none yes 96 1 no h no no no 31 flat 200 no no md no none rel2rel_2 yes 96 1 no h no no no 32 flat all no no md no sin none yes 96 1 no h no no no 33 flat all no no md no sin none no 96 1 no h no no no 34 flat all no no md no sin none yes 96 1 no xyz no no yes 35 flat 200 no no md no none rel7rel_7 yes 96 1 no h no no no 36 flat 200 no no md no none rel5rel_5 yes 96 1 no h no no no 37 flat 200 no no md no none rel3rel_3 yes 96 1 no h no no no 38 flat all no no md no sin none yes 96 1 no h no no yes 39 flat all no no md no sin none yes 96 1 yes h no no no 40 both all - - - - - - no - - - h no no no Table 2: Options used for the different track training runs. The models are ordered according to their ID which is determined by average total test loss (see Fig. 7(b)). Explanations for the numeric option mapping are given in table 4 and for some abbreviations are given in table 5. Model 40 is a GNN and the related transformer-related options are marked with “-”. Shower training Numeric options (see tables 4 and 5 for definitions) Model ID 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 1 flat all yes yes md no none none no 96 1 no h no no no 2 flat all yes yes md yes none none no 96 1 no h no no no 3 flat 400 no yes md no none none no 96 1 no h no no no 4 flat 200 yes yes md no none rel1rel_1 no 96 1 no h no no no 5 flat all no yes md no none none no 96 1 no h no no no 6 flat 200 no no md no none rel1rel_1 yes 96 1 no h no no no 7 flat 200 no no md no none rel5rel_5 yes 96 1 no h no no no 8 flat 200 no yes c3 no none rel1rel_1 no 96 1 no h no no no 9 flat 200 no no md no none rel4rel_4 yes 96 1 no h no no no 10 flat 200 no yes c1 no none rel1rel_1 no 96 1 no h no no no 11 flat 200 no no md no none rel6rel_6 yes 96 1 no h no no no 12 flat 200 no yes md no none none no 96 1 no h no no no 13 flat 200 no yes c2 no none rel1rel_1 no 96 1 no h no no no 14 flat 200 no no md no none rel2rel_2 yes 96 1 no h no no no 15 flat 200 no yes md no none rel1rel_1 no 96 1 no h no no no 16 flat 200 no no md no none rel7rel_7 yes 96 1 no h no no no 17 flat 200 no no md no none none yes 96 1 no h no no no 18 flat 200 no no md no none rel3rel_3 yes 96 1 no h no no no 19 flat all no no md no none none yes 96 1 no h no no no 20 flat all no no md no sin none yes 96 1 no h no no no 21 flat 200 - - - - - - no - - - h no no no Table 3: Options used for the different shower training runs. The models are ordered according to their ID which is determined by average total test loss (see Fig. 7(a)). Explanations for the numeric option mapping are given in table 4 and for some abbreviations are given in table 5. Model 21 is a GNN and the related transformer-related options are marked with “-”. Numeric option abbreviations: 1: train data weighting 2: max DOMs 3: nonlinear in-projection 4: ReSi Dual 5: aggregation mode 6: add mean diff. for in-projection 7: abs. pos. encoding 8: rel. pos. encoding 9: bottleneck 10: transformer computing dim 11: numheads 12: extra layer norm 13: rotation type 14: randomize xyz of DOMs during training 15: ignore saturation for COG calc. 16: simple vMF flow Table 4: A list of numerical identifiers for different model hyperparameters. train data weighting flat Flat distribution in deposited energy. si Weighting data to a spectrum with spectral index -1.8. both Arithmetic mean of spectral index (-1.8) + equal weighting in deposited energy. rel. positional encoding rel1rel_1 Relative value position encoding in first layer, using standard relative value positional encoding. rel2rel_2 Relative value position encoding using just the first layer, also adding the absolute input token before projection to the positional encoding input. rel3rel_3 Relative value position encoding using just the first layer, but in parallel to another normal encoding. rel4rel_4 Relative value position encoding in first layer with both overall input and value token, and in layer 2) and 3) with just the previous value token. rel5rel_5 Relative value position encoding in layers 1)-5) using both the absolute input token before the transformer concatenated with the respective value token in each layer. rel6rel_6 Relative value position encoding in interleaved layers 1), 3), 6), 9) and 12) using both the absolute input token before the transformer concatenated with the respective value token in each layer. rel7rel_7 first layer with both overall input and value token, and in layer 5) and 10) with just the previous value token. aggregation mode md Mean + diagonal variance of all tokens (2 X vector size of just mean). c1 Learnable class token, added as a learnable bias in the first layer after QKV projection. c2 Learnable class token, added as an extra learnable token in the first layer, treated just as normal token. Like originally introduced in [7]. c3 Improved class token that is treated separately from all other tokens, interacts via cross attention and has its own separate MLP layers. In combination with “ReSi Dual”, it also uses two residual streams. rotation type h Three stacked householder reflectors for rotations. xyz Standard mean parametrization of rotation matrix. Table 5: Descriptions of some option abbreviations. References [1] C. A. Argüelles, A. Schneider, and T. Yuan (2019) A binned likelihood for stochastic models. 2019 (6), p. 30. External Links: ISSN 1029-8479, Document, Link Cited by: §6.3. [2] H. Bukhari, D. Chakraborty, P. Eller, et al. (2024) IceCube – Neutrinos in Deep Ice. 84 (6), p. 646. External Links: ISSN 1434-6052, Document, Link Cited by: §5.1. [3] P. A. Cherenkov (1934) Visible luminescence of pure liquids under the influence of γ-radiation. 2 (8), p. 451–454. External Links: Document Cited by: §1.1, §4. [4] K. Cho, B. van Merriënboer, C. Gulcehre, et al. (2014) Learning phrase representations using RNN encoder–decoder for statistical machine translation. In Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP), Doha, Qatar, p. 1724–1734. External Links: Link, Document Cited by: §4.3. [5] K. Cranmer, J. Brehmer, and G. Louppe (2020) The frontier of simulation-based inference. 117 (48), p. 30055–30062. External Links: Document, Link Cited by: §1.2, §2. [6] T. Dao (2024) FlashAttention-2: faster attention with better parallelism and work partitioning. In International Conference on Learning Representations (ICLR), Cited by: §5. [7] A. Dosovitskiy, L. Beyer, A. Kolesnikov, et al. (2021) An Image is Worth 16x16 Words: Transformers for Image Recognition at Scale. In International Conference on Learning Representations, External Links: Link Cited by: Table 5, §5.1, §6.1. [8] C. Durkan, A. Bekasov, I. Murray, et al. (2019) Neural Spline Flows. In Advances in Neural Information Processing Systems, Vol. 32, p. 7511 – 7522. External Links: Link Cited by: Appendix B, Appendix B, Appendix B, §1.2, §3.4. [9] A. Fedynitch, F. Riehn, R. Engel, T. K. Gaisser, and T. Stanev (2019) Hadronic interaction model sibyll 2.3c2.3c and inclusive lepton fluxes. 100, p. 103018. External Links: Document, Link Cited by: §6.3. [10] P. Fernique, T. Boch, T. Donaldson, et al. (2014) MOC - HEALPix Multi-Order Coverage map Version 1.0. External Links: Link, Document Cited by: §3.5. [11] R. Fisher (1953) Dispersion on a Sphere. 217 (1130), p. 295–305. External Links: Document Cited by: §3.4. [12] T. K. Gaisser (2012) Spectrum of cosmic-ray nucleons, kaon production, and the atmospheric muon charge ratio. 35 (12), p. 801–806. External Links: ISSN 0927-6505, Document, Link Cited by: §6.3. [13] T. Glüsenkamp Jammy Flows. Note: https://github.com/thoglu/jammy_flows Cited by: §3.4. [14] T. Glüsenkamp (2024) Unifying supervised learning and VAEs: coverage, systematics and goodness-of-fit in normalizing-flow based neural network models for astro-particle reconstructions. The European Physical Journal C 84, p. 163. External Links: Document, Link Cited by: §2, §3.4, §5, §6.3. [15] K. M. Górski, E. Hivon, A. J. Banday, et al. (2005) HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere. 622 (2), p. 759. External Links: Document, Link Cited by: §3.5. [16] J. A. Gregory and R. Delbourgo (1982) Piecewise Rational Quadratic Interpolation to Monotonic Data. 2 (2), p. 123–130. External Links: ISSN 0272-4979, Document, Link Cited by: Appendix B. [17] N. J. Higham (2002) Accuracy and Stability of Numerical Algorithms. Second edition, Society for Industrial and Applied Mathematics, . External Links: Document, Link Cited by: Appendix C, §3.4. [18] P. Izmailov, D. Podoprikhin, T. Garipov, et al. (2018) Averaging Weights Leads to Wider Optima and Better Generalization. In Proceedings of the Thirty-Fourth Conference on Uncertainty in Artificial Intelligence, UAI 2018, Monterey, California, USA, August 6-10, 2018, p. 876–885. External Links: Link Cited by: §5. [19] W. Jakob (2012) Numerically stable sampling of the von Mises Fisher distribution on S2S^2 (and other tricks). External Links: Link Cited by: §3.4, footnote 2. [20] M. I. Jordan, Z. Ghahramani, T. S. Jaakkola, et al. (1999) An Introduction to Variational Methods for Graphical Models. 37 (2), p. 183–233. External Links: ISSN 1573-0565, Document, Link Cited by: §2. [21] D. P. Kingma and J. Ba (2015) Adam: A Method for Stochastic Optimization. In Proceedings of the 3rd International Conference on Learning Representations, Cited by: §5. [22] S. Kullback and R. A. Leibler (1951) On information and sufficiency. 22 (1), p. 79–86. External Links: ISSN 00034851, Link Cited by: §2. [23] C. L. Lawson and R. J. Hanson (1974) Solving least squares problems. Englewood Cliffs, NJ: Prentice-Hall (English). External Links: ISBN 0-13-822585-0 Cited by: §4.1. [24] I. Loshchilov and F. Hutter (2017) SGDR: Stochastic Gradient Descent with Warm Restarts. In International Conference on Learning Representations, External Links: Link Cited by: §5. [25] S. Mandt, M. Hoffman, and D. Blei (2016) A Variational Analysis of Stochastic Gradient Algorithms. In Proceedings of The 33rd International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 48, New York, New York, USA, p. 354–363. External Links: Link Cited by: §5. [26] I. Martinez-Castellanos, L. P. Singer, E. Burns, et al. (2022) Multiresolution HEALPix Maps for Multiwavelength and Multimessenger Astronomy. 163 (6), p. 259. External Links: Document, Link Cited by: §3.5. [27] I. D. Mienye, T. G. Swart, and G. Obaido (2024) Recurrent Neural Networks: A Comprehensive Review of Architectures, Variants, and Applications. 15 (9). External Links: Link, ISSN 2078-2489, Document Cited by: §4.3. [28] T. Neunhöffer (2006) Estimating the angular resolution of tracks in neutrino telescopes based on a likelihood analysis. 25 (3), p. 220–225. External Links: ISSN 0927-6505, Document, Link Cited by: §6.2, §6.4. [29] G. Papamakarios and I. Murray (2016) Fast ε -free inference of simulation models with Bayesian conditional density estimation. In Proceedings of the 30th International Conference on Neural Information Processing Systems, NIPS’16, Red Hook, NY, USA, p. 1036–1044. External Links: ISBN 9781510838819 Cited by: §1.2, §2, §2. [30] D. J. Rezende, G. Papamakarios, S. Racaniere, et al. (2020) Normalizing Flows on Tori and Spheres. In Proceedings of the 37th International Conference on Machine Learning, Vol. 119, p. 8083–8092. External Links: Link Cited by: §B.2, §1.2, §3.2, §3.4, §3.4. [31] M. Richey (2010) The Evolution of Markov Chain Monte Carlo Methods. 117 (5), p. 383–413. External Links: Document, Link Cited by: §2. [32] K. Schatto (2014) Stacked searches for high-energy neutrinos from blazars with IceCube. Ph.D. Thesis, University of Mainz. Cited by: §1.1, Figure 8, Figure 8, §6.2, §6.2. [33] P. Shaw, J. Uszkoreit, and A. Vaswani (2018) Self-Attention with Relative Position Representations. In Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 2 (Short Papers), New Orleans, Louisiana, p. 464–468. External Links: Link, Document Cited by: §5.1. [34] The IceCube Collaboration (2014) Energy reconstruction methods in the IceCube neutrino telescope. 9 (03), p. P03009. External Links: Document, Link Cited by: §1.1. [35] The IceCube Collaboration (2014) Observation of High-Energy Astrophysical Neutrinos in Three Years of IceCube Data. 113, p. 101101. External Links: Document, Link Cited by: §1.1, §1.1. [36] The IceCube Collaboration (2017) The IceCube Neutrino Observatory: instrumentation and online systems. 12 (03), p. P03012. External Links: Document, Link Cited by: §1, §4. [37] The IceCube Collaboration (2017) The IceCube realtime alert system. 92, p. 30–41. External Links: ISSN 0927-6505, Document, Link Cited by: §1.1. [38] The IceCube Collaboration (2018) Neutrino emission from the direction of the blazar TXS 0506+056 prior to the IceCube-170922A alert. 361 (6398), p. 147–151. External Links: Document, Link Cited by: §1.1. [39] The IceCube Collaboration (2021) A convolutional neural network based cascade reconstruction for the IceCube Neutrino Observatory. 16 (07), p. P07041. External Links: Document, Link Cited by: §1.1, §1.1, §4.1. [40] The IceCube Collaboration (2021) A muon-track reconstruction exploiting stochastic losses for large-scale Cherenkov detectors. 16 (08), p. P08034. External Links: Document, Link Cited by: §1.1, §1.1, §6.2. [41] The IceCube Collaboration (2021) Combining Maximum-Likelihood with Deep Learning for Event Reconstruction in IceCube. In Proceedings of 37th International Cosmic Ray Conference — PoS(ICRC2021), Vol. 395, p. 1065. External Links: Document Cited by: §1.1, §1.1. [42] The IceCube Collaboration (2021) IceCube-Gen2: the window to the extreme Universe. Journal of Physics G: Nuclear and Particle PhysicsarXiv:2510.04762Astroparticle PhysicsJournal of InstrumentationACM Trans. Graph.Journal of High Energy PhysicsPhys. Rev. DAstroparticle PhysicsPhys. Rev. DAstroparticle PhysicsJournal of InstrumentationScienceJournal of InstrumentationThe Astronomical JournalInternational Virtual Observatory AlliancePhys. Rev. Lett.The American Mathematical MonthlyMachine LearningProceedings of the Royal Society of London Series AEPFL research reportIMA Journal of Numerical AnalysisThe Astrophysical JournalDokl. Akad. Nauk SSSRProceedings of the National Academy of SciencesJournal of InstrumentationJournal of InstrumentationComputer Physics CommunicationsPhys. Rev. Lett.Phys. Rev. DJournal of InstrumentationThe Annals of Mathematical StatisticsThe European Physical Journal CarXiv:2304.14802Information 48 (6), p. 060501. External Links: Document, Link Cited by: §7. [43] The IceCube Collaboration (2022) Evidence for neutrino emission from the nearby active galaxy NGC 1068. Science 378 (6619), p. 538–543. External Links: Document, Link Cited by: §1.1, §1, Table 1, Figure 11, Figure 11, §6.2, §6.3. [44] The IceCube Collaboration (2022) Graph Neural Networks for low-energy event classification & reconstruction in IceCube. 17 (11), p. P11003. External Links: Document, Link Cited by: §1.2, §4.3, §5.1, §6.1. [45] The IceCube Collaboration (2023) A model independent parametrization of the optical properties of the refrozen IceCube drill holes. In Proceedings of 38th International Cosmic Ray Conference — PoS(ICRC2023), Vol. 444, p. 1034. External Links: Document Cited by: §5. [46] The IceCube Collaboration (2023) Conditional normalizing flows for IceCube event reconstruction. In Proceedings of 38th International Cosmic Ray Conference — PoS(ICRC2023), Vol. 444, p. 1003. External Links: Document Cited by: §1.2. [47] The IceCube Collaboration (2023) Observation of high-energy neutrinos from the Galactic plane. Science 380 (6652), p. 1338–1343. External Links: Document, Link Cited by: §1.1, §1.1, §1, Table 1, Figure 10, Figure 10, §6.2, §6.3. [48] The IceCube Collaboration (2024) Characterization of the astrophysical diffuse neutrino flux using starting track events in IceCube. 110, p. 022001. External Links: Document, Link Cited by: §1.1. [49] The IceCube Collaboration (2024) Improved modeling of in-ice particle showers for IceCube event reconstruction. Journal of Instrumentation 19 (06), p. P06026. External Links: Document, Link Cited by: §1.1, Table 1, Figure 8, Figure 8, §6.2, §6. [50] The IceCube Collaboration (2024) The IceCube Upgrade: status and prospects for advances with GeV neutrinos. In Proceedings of 42nd International Conference on High Energy Physics — PoS(ICHEP2024), Vol. 476, p. 149. External Links: Document Cited by: §7. [51] The IceCube Collaboration (2025) State of the Ice Model in the IceCube Observatory. In Proceedings of 39th International Cosmic Ray Conference — PoS(ICRC2025), Vol. 501, p. 1013. Cited by: §1.1, Table 1, Table 1, §5, §6. [52] The IceCube Collaboration (2026) Improved measurements of the TeV-PeV extragalactic neutrino spectrum from joint analyses of IceCube tracks and cascades. 113, p. 062002. External Links: Document, Link Cited by: §6.3. [53] A. Vaswani, N. Shazeer, N. Parmar, et al. (2017) Attention is All you Need. In Advances in Neural Information Processing Systems, Vol. 30, p. . External Links: Link Cited by: §1.2, §4.2, §4.2. [54] Y. Wang, Y. Sun, Z. Liu, et al. (2019) Dynamic Graph CNN for Learning on Point Clouds. 38 (5). External Links: ISSN 0730-0301, Link, Document Cited by: §1.2, §4.3. [55] N. Whitehorn, J. van Santen, and S. Lafebre (2013) Penalized splines for smooth representation of high-dimensional monte carlo datasets. 184 (9), p. 2214–2220. External Links: ISSN 0010-4655, Document, Link Cited by: §1.1. [56] S. Xie, H. Zhang, J. Guo, et al. (2024) ResiDual: Transformer with Dual Residual Connections. External Links: Link Cited by: §5.1, §6.1. [57] R. Xiong, Y. Yang, D. He, et al. (2020) On layer normalization in the transformer architecture. In Proceedings of the 37th International Conference on Machine Learning, ICML’20. Cited by: §4.2.