\documentclass[prd, twocolumn, nofootinbib, floatfix]{revtex4-1}
%\usepackage{silence}
%\usepackage{draftwatermark}
\usepackage{dcolumn}    % For aligning numbers at the decimal point

% Define a new column type "d" for numeric columns.
% The syntax is: D{input sep}{output sep}{number of digits}
% For example, d{6} expects up to 6 digits (you can adjust this as needed).
\newcolumntype{d}[1]{D{.}{.}{#1}}
\usepackage{lmodern} % math, rm, ss, tt
\usepackage[T1]{fontenc}
\usepackage{textcomp}
\usepackage[scr=boondoxo,frak=boondox,bb=boondox]{mathalfa}
\pdfoutput=1
\usepackage{makecell, multirow}
\setcellgapes{3pt}
\usepackage{hyperref}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{tabularx}
\usepackage{tikz}
\usepackage[vlines]{tabularht}
\usepackage{adjustbox}
\usepackage{tikz}
\usepackage[]{draftfigure}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{float} 
\usepackage{color}
\usepackage{color,soul}
\usepackage{ulem}
\usepackage{multirow}
\usepackage{epsfig}
\usepackage{epstopdf}
\usepackage{rotating}
\usepackage{afterpage}
\usepackage{soul}
\usepackage{fancyvrb}
\usepackage{lineno}
\usepackage{amsmath}
\usepackage{xspace}
\usepackage{hyperref}
\usepackage{url}
\usepackage{amssymb,latexsym,mathrsfs}
\usepackage{multirow,bigdelim, booktabs}
\usepackage{bigints}
\usepackage{soul}
\usepackage{graphicx}
\usepackage{subcaption}
%\usepackage{subfig}
\usepackage{dcolumn}
\usepackage{bm}
\usepackage{epsfig}
\usepackage{array}
\usepackage{caption}
\usepackage{color}
\usepackage{multirow}
\usepackage{makecell}
\usepackage{textcomp}
\usepackage[title]{appendix}
\usepackage{tabularx}
\usepackage{natbib}
\usepackage{graphicx}
\setcitestyle{authoryear,open={(},close={)}} 

\captionsetup{justification=raggedright,singlelinecheck=false}
\hypersetup{
    colorlinks=true,
    linkcolor=blue,
    citecolor=blue,
} 


\definecolor{mediumpurple}{rgb}{0.58, 0.44, 0.86}





\begin{document} 
%\linenumbers
\title{  A Fully Photometric Approach to  Type Ia Supernova Cosmology in the LSST Era: Host Galaxy  Redshifts and  Supernova  Classification
}

\input{PLASTICC_sims_newcommand.tex}



\author{
Ayan Mitra${}^{1,2}$,
Richard Kessler${}^{3,4}$,
Rebecca~C.~Chen${}^{5,6,7}$,
Alex~Gagliano${}^{8,9,10}$,
Matthew~Grayling${}^{11}$,
Surhud~More${}^{12}$,
Gautham~Narayan${}^{2}$,
Helen~Qu${}^{13}$,
Srinivasan~Raghunathan${}^{1,14}$,
Alex~I.~Malz${}^{15}$,
Michelle~Lochner${}^{16}$,
\\
The LSST Dark Energy Science Collaboration \\
%
${}^{1}$ Center for AstroPhysical Surveys, National Center for Supercomputing Applications, Urbana, IL 61801, USA \\
${}^{2}$ Department of Astronomy, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA \\
${}^{3}$ Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA \\
${}^{4}$ Department of Astronomy and Astrophysics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA \\
${}^{5}$ Kavli Institute for Particle Astrophysics \& Cosmology, P.O. Box 2450, Stanford University, Stanford, CA 94305, USA \\ 
${}^{6}$ SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA \\
${}^{7}$ Brinson Prize Fellow \\
${}^{8}$ The NSF AI Institute for Artificial Intelligence and Fundamental Interactions \\
${}^{9}$ Center for Astrophysics \textbar{} Harvard \& Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA \\
${}^{10}$ Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA \\ 
${}^{11}$ Institute of Astronomy and Kavli Institute for Cosmology, Madingley Road, Cambridge CB3 0HA, UK \\
${}^{12}$ The Inter-University Centre for Astronomy and Astrophysics (IUCAA), Post Bag 4, Ganeshkhind, Pune 411007, India \\
${}^{13}$ Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY 10010, USA \\
${}^{14}$ Department of Physics \& Astronomy, University of California, One Shields Avenue, Davis, CA 95616, USA\\
${}^{15}$ Space Telescope Science Institute, Baltimore, MD, USA \\
${}^{16}$ Department of Physics and Astronomy, University of the Western Cape, Bellville, Cape Town 7535, South Africa
}


    

\date{\today} 

\begin{abstract}
The upcoming Vera C. Rubin Observatory's Legacy Survey of Space and Time (LSST) is expected to discover 
 {nearly a} million  Type Ia supernovae (SNe~Ia), offering an unprecedented opportunity to constrain dark energy. 
The vast majority of these events will lack spectroscopic classification and redshifts, 
necessitating a fully photometric approach  {to maximize cosmology constraining power.}
We present detailed simulations based on the 
Extended LSST Astronomical Time Series Classification Challenge (\acro), 
and a cosmological analysis using photometrically classified SNe~Ia with host galaxy photometric redshifts.
This dataset features realistic multi-band light curves,  {non-SNIa contamination,} 
host mis-associations, and transient-host correlations across the 
high-redshift Deep Drilling Fields (DDF) ($\sim 50$ deg$^2$).
We  {also} include a  {spectroscopically confirmed} 
low-redshift sample based on the Wide Fast Deep (WFD) fields.
%
We employ a 
joint SN+host photometric redshift fit,  
a neural network based photometric classifier (SCONE),
and BEAMS with Bias Corrections (BBC) methodology to construct a bias-corrected Hubble diagram. 
We produce statistical + systematic covariance matrices, 
and perform cosmology fitting with a prior   {using}    Cosmic Microwave Background    {constraints}.
We fit and present results for the \wCDM\ dark energy model, 
and the more general Chevallier-Polarski-Linder (CPL) \wwCDM\ model.  
%
With a simulated sample of ${\sim}6000$  {events}, we achieve  {a} 
Figure of Merit (FoM)  {value} of about   {$150$}, 
 {which is significantly larger than the \DESVYR\ FoM of $54$}.
%
 {Averaging analysis results over 25 independent samples, we find small but significant biases
indicating a need for further analysis testing and development.}
%
\end{abstract}

\maketitle
\tableofcontents
\section{Introduction }
\label{sec:intro}

The discovery of dark energy and the accelerating expansion of the Universe, 
first identified through the observation of a few dozen Type Ia supernovae (SNe Ia) \cite{adam,perl}, marked a pivotal moment in modern cosmology. Subsequent SNIa surveys yielded constraints on the dark energy equation of state that  agree well with a cosmological constant (\(w=-1\)) \citep{SNLS:2005qlf,SupernovaSearchTeam:2004lze, 2018ApJ...859..101S, ESSENCE:2007acn,  conley11, kessler3, 2012ApJ...746...85S, Betoule2014}, 
where $w$ is the ratio of pressure to energy density for dark energy 
(for reviews see e.g. \citet{Frieman:2008sn, 2013PhR...530...87W}). 

Few-$\sigma$ hints of evolving dark energy 
have emerged from analyses with higher statistics and better systematics handling from 
Pantheon+ \citep{pantheon_new}
the Dark Energy Survey (DES)  \citep[\DESVYR]{DES:2024jxu},
the Dark Energy Spectroscopic Instrument (DESI) \citep[DESI]{DESI:2024mwx},
and the  UNION {3} sample Bayesian analysis \citep{Rubin:2023ovl}.
While these deviations from a cosmological constant have not reached discovery significance, they have spurred considerable interest in the possibility that dark energy may be dynamic rather than static. 
  {Current sample statistics include $\sim 10^3$ SNe, and }
upcoming projects such as  the 
Rubin Legacy Survey of Space and Time \citep[LSST]{Cahn2009,ivezic} and the  
Roman Space Telescope \citep{roman1,w2} will 
 significantly  {improve} the statistical   {precision on dark energy parameters,} and   {a corresponding reduction in systematics is needed to make full use of these samples.} 



Historically, SN~Ia cosmology surveys 
relied on intensive spectroscopic follow-up to confirm each transient’s type and to secure a host-galaxy spectroscopic redshift  {(spec-$z$)}. 
These spectroscopically confirmed samples 
 {include  ${\sim}10^3$  SN~Ia, and}  {contribute to}  
composite data sets like the Pantheon compilation \citep{Scolnic2018} and   {the} Pantheon+ extension \citep{pantheon+,pantheon_new}. 

Previous photometric surveys%
\footnote{
CfA: Center for Astrophysics Supernova Program;
CSP: Carnegie Supernova Project;
SNLS: Supernova Legacy Survey;
ESSENCE: Equation of State: SupErNovae trace Cosmic Expansion;
SDSS-II: Sloan Digital Sky Survey-II Supernova Survey;
Pan-STARRS: Panoramic Survey Telescope and Rapid Response System;
DES-SN: Dark Energy Survey Supernova Program.
}
(CfA, CSP, SNLS, ESSENCE, SDSS-II, Pan-STARRS, DES-SN)
acquired spectroscopic information for only a small fraction of their samples, and
this challenge will be even more pronounced in the LSST era, which will have a significantly higher detection rate.  {The 4MOST TiDES \citep{4most_tides} program anticipates $18000$ spectroscopic SNe observations, 
 {which} 
is a very small fraction of the total number of transients.} 
Therefore, for the vast majority of newly discovered SNe  it has become necessary to develop  photometric classification and a formalism to analyse a contaminated Hubble diagram. 




 Pan-STARRS1 \citep{Jones2018,FSS:2018cey} used the PSNID classifier \citep{sako}  and the  Bayesian Estimation Applied to Multiple Species (BEAMS) \citep{kunz}  formalism to determine cosmology constraints from a photometrically identified sample. For \DESVYR\ they used more advanced photometric classifiers:  SuperNNova (SNN) \footnote{\url{https://github.com/supernnova/SuperNNova}} \citep{snn} and SCONE \citep{scone} (see Sec.~\ref{subsec:scone}). Both  have exhibited very high levels of purity and accuracy $> 98\%$.
 They also used an extension of the BEAMS formalism, ``BEAMS with Bias Correction {s}'' (BBC), that is described in  \citet{bbc}[hereafter KS17]. 
 
 




While photometric classification has been successfully used to deliver cosmological results, 
photometric redshifts  {(photo-$z$)} have not. 
A photo-$z$ method using the SN+host galaxy has been proposed 
\citep{Kessler2010,Palanque2010,zbeams}, 
and \citet[M23]{Mitra:2022ykq} performed the first rigorous SNIa-cosmology analysis with photo-$z$'s and systematics, but did not include non-SNIa contamination.  {\citet{Chen25} examined the photo-$z$ bias in the \DESVYR\ sample by comparing $w$ results using spec-$z$'s vs photo-$z$. 
\citet{Ruhlmann-Kleider:2022vum} measured $\Omega_m$ using photo-$z$'s in SNLS data, and 
more recently, \citet{2025arXiv250815899K}
analysed simulated data  with photo-$z$ using a simulation based  inference technique.


As part of the DESC\footnote{\href{https://lsstdesc.org/}{Dark Energy Science Collaboration}} effort to prepare  {for} cosmology analyses with the LSST data,
here we simulate LSST data that includes both host galaxy photo-$z$'s and non-SNIa contamination,
and we incorporate SN+host photo-$z$ and BEAMS into the analysis. This analysis supersedes M23, and is applicable to real data. M23 presented a thorough light curves to cosmology analysis with photometric redshift ($\zphot$) using  LSST simulated data. However,  M23 was focused on 
evaluating systematics related to photo-z, and therefore the analysis is not practical on real data because of the simplifying assumption that
every SNIa is perfectly classified (i.e. spectroscopic sample).
%%%%%

A  full-scale analysis involving both photometric classification and 
$\zphot$ with extensive systematics has not yet been done.  In this paper, 
 {we simulate SNIa with host galaxy photometric redshifts and non-SNIa contamination,
and perform an analysis using SN+host photo-$z$s and photometric classification.}
%
The historical trajectory and evolving methodologies of SNIa cosmology research are summarised in Table~\ref{tab:sn_history}. 
 {While we anticipate spectroscopic resources for a small fraction of the events (eg. 4MOST TiDES \citep{4most_tides}), 
here we exclude all spectroscopic information in the high-$z$ simulation so that the analysis is
more sensitive to potential biases from a photometric analysis.}

\begin{table*}[ht]
\centering
\begin{tabular}{|p{3.5cm}|p{4cm}|p{7cm}|}
\hline
\textbf{Classification} & \textbf{Redshift} & \textbf{Survey/Project} \\ \hline \hline
Spectroscopic  & Spec Only (SN or Host) & Two decades of results \newline 
Hi-z \citep{adam}, SCP \citep{perl}
\newline SNLS \citep{SNLS:2005qlf}, ESSENCE \citep{ESSENCE:2007acn},  SDSS-II \citep{2009ApJS..185...32K}, PS1 \citep{panstarrs},   {RAISIN \citep{Jones:2022mvo}}, \newline Pantheon \citep{2018ApJ...859..101S, Scolnic2018}, JLA \citep{Betoule2014}, UNION3 \citep{Rubin:2023ovl}, DES-3YR \citep{DES3YR} etc.\\ \hline
Photometric  & Spec Only (SN or Host) & PS1 \citep{Jones2018}, DES-SN5YR \citep{DES:2024hip}, Amalgame \citep{Popovic24} \\ \hline
Spectroscopic  & Spec or Photo (Host Only) & DES-5YR + redMaGiC \citep{Chen2022}, LSST  {forecast} [M23] \citep{Mitra:2022ykq}. \\ \hline
Photometric  & Spec or Photo (Host Only) &   { Roman forecast \citep{Kessler:2025eib} } and \textbf{ {this work}}.\\ 
&&\\ \hline 
\end{tabular}
\caption{ {History of SNIa-based $w$ measurements (with systematics)  from previous surveys and from forecasts for future surveys}.  {History starts with spectroscopic redshift and classifications and evolves to photometric classifications and redshift.}}
\label{tab:sn_history}
\end{table*}

This paper is organised as follows, in sec.~\ref{sec:sim} {,} we describe the simulation process
 to generate the data and  bias corrections. 
In sec.~\ref{sec:anaysis} {,} we describe the main analysis of this paper 
including light curve fitting,  selection cuts, systematic uncertainties, photometric classification  using SCONE, 
 {and BBC to produce a bias-corrected Hubble diagram (HD).}
 {In Sec.~\ref{sec:cosmology}, the HD is}  used for  
dark energy  {parameter} estimation. 
%
Finally in sec.~\ref{sec:conclude} {,} we present the conclusion of this analysis. 





% =========================================================================
%\section{Overview of LSST and Dark Energy Science Collaboration } 
%\label{sec:lsst}



\section{Simulated Data}
\label{sec:sim}

\subsection{Overview}
Ideally, we would work with simulated images from the Data Challenge 2 (DC2) \citep{dc2} \citep{2020MNRAS.497..210S} or 
OpenUniverse2024 \citep{LSSTDarkEnergyScience:2025lah},
and obtain light curves from
the LSST difference imaging analysis (DIA)\footnote{\URLDIA} 
based on \citet{DIA1}.
While such an analysis is most realistic, it is very CPU intensive and repeating the analysis is not practical.
Instead,  {we generate catalogue level simulations of lightcurves  corresponding to  }  the output of DIA   {using a two-detection trigger}.

We simulate and analyse \NSAMPLE\ independent samples, where each sample includes two subsets:
\textbf{(1)} a complete low-$z$ $(z\le 0.08)$ sample of {\spec}ally confirmed SNe~Ia with accurate \spec\ redshifts, and \textbf{(2)} a high-$z$  $(z< 1.55)$ SN~Ia sample with non-SNIa contamination from SNII/Ib/Ic and peculiar SNIax and 91bg-like.  The high-$z$ subset includes host-galaxy photo-$z$'s, but does not have \spec\ redshifts nor identification. 

The low-$z$ cadence  is  from the the Wide Fast Deep (WFD) observing strategy from DC2. 
The high-$z$ cadence is from ${\sim5}0$~deg$^2$ Deep Drilling Fields (DDF)\footnote{\URLDDF} 
 {used in the PLaSTiCC data challenge \citep{plasticc_H2020, plasticc_K2019},\footnote{More recent LSST cadences are available \href{https://www.lsst.org/content/charge-survey-cadence-optimization-committee-scoc}{here}}
 {and corresponds to 3 years, between $59570$-$60675$ MJD days.}} 
 {While we anticipate that the number of supernova in WFD will significantly exceed that in DDF, 
here we focus on DDF to better study  {photometric methods}, 
 {and to forecast constraints for a potential early science analysis using a few years of DDF where we
expect templates to be available well before WFD.}
}  

%%%


To model host properties and incorrect host associations, we use host galaxy libraries ({\HOSTLIB}s) from the \acro\ (Extended LSST Astronomical Time Series Classification Challenge) \citep{2023AAS...24111701N}. The simulation for SNIa and detection efficiency   is the same in both subsets. 
 

\subsubsection*{ELaSTiCC}
\label{sec:elasticc} 
To help prepare \spec\ follow-up for transients and hosts found by Rubin,
\acro\ was developed to enable early and accurate photometric classification of transients using both light curves and contextual information \citep{Foley:2013jrk,Baldeschi:2020whr, Gagliano:2020ucg}. 
\acro\ includes  $\sim$5 million detected transients and $\sim$50 million alerts across 30+ classes, 
embedding realistic transient–host correlations from the CosmoDC2 catalog \citep{korytov}. 
 {While we don't directly use the final \acro\ simulations here, we leverage the host galaxy library  developed  for \acro\ which includes host-SN correlations, and we include a few 
 modifications
 
to model host mis-associations.}
%

%--------------


\subsection{Implementation}
\label{sec:snana}
The SNANA  {(Supernova Analysis package)} \footnote{\url{https://github.com/RickKessler/SNANA}} simulation includes  three main steps
as described in \citet{kessler2019}:
\textbf{(a)} estimate the true source magnitude from a rest-frame SED model that incorporates
redshifting and cosmic expansion ($\OM=0.315$, $w=-1$, flatness),
lensing magnification, peculiar velocity, Milky way extinction, and 
host galaxy extinction; 
\textbf{(b)}  model the flux uncertainty and random fluctuations using
sky noise,  zero-point, and PSF 
for each observation;
\textbf{(c)} model DIA detection efficiency vs. SNR,  {and trigger.}
%


We simulate SNe~Ia  using the SALT3 SED model \citep{salt,salt3}
with a NIR extension from \citet{Pierel:2022pqc} 
 {that covers wavelength range 2,000-20,000\AA\ in the rest-frame. 
To avoid classification artifacts from $u$-band dropouts at high-$z$,
the SED model is linearly extrapolated down to zero flux at 500\AA.}


 {To model instrinsic scatter about the mean SALT3 model, } we use a dust based  model from  \citet{BS20}, along with the stretch and colour populations from  \citet{Popovic:2021yuo}. 


The stretch-luminosity parameter is $\alphaTrue=0.15$, and the intrinsic color-luminosity parameter 
$\betaIntTrue=2.12$.
The combination of intrinsic color variation and dust results in an effective fitted $\beta\sim 3$  {in the analysis}  
(sec 3.4  in \citet{DES:2024hip}),
but there is no true effective $\beta$ value to compare with.


Correlations between SNe and host-galaxy mass are modelled using \acro\ {\HOSTLIB}s and probability distributions \citep{lokken}.

The average $5\sigma$ limiting magnitudes per co-added night,  and their average time between nights for each of the \bands\ passbands, are listed in Table~\ref{t_lowz} 
for both the low-$z$ (WFD) and the high-$z$ (DDF) subsets. On average the DDF is  $\sim 1.0$ mag deeper  than the WFD, with $g$ and $r$ being the deepest bands. Compared to DDF, the WFD cadence has $30\%$ fewer co-added observations on average.   For  host galaxies  we use a $1$ year co-added depth which corresponds to a depth of $i \sim 25.4$.   

 {
A simulated event is recorded if it passes the 2-detection trigger,
where any two observations, separated by at least 30 minutes, are detected
based on the DIA efficiency vs. detection in Fig. 17 of \citet{Sanchez}.
The magnitudes where the DIA efficiency is 50\% are $23.66, 24.69, 24.06, 23.45, 22.54, 21.62$
for $u,g,r,i,z,Y$ bands, respectively.}


\subsubsection*{Rates}

For each SN type in the high-$z$ sample, we use the redshift range, survey duration, sky area, and
volumetric rate to compute the physical rate defined as the true number of events without instrumental selection.

For SN~Ia, we adopt a redshift-dependent volumetric rate model,
\begin{eqnarray}
     R(z)  & = & 2.5  \times 10^{-5} (1+z)^{1.5}~{\RateUnit} ~~(z<1) \label{eq:rate_lowz} \\
     R(z) & = & 9.7  \times 10^{-5} (1+z)^{-0.5}~{\RateUnit} ~(z>1)~ \label{eq:rate_highz}
\end{eqnarray}
based on \citet{Dilday2008} at $z<1$ and \citet{Hounsell2018} at $z>1$. 
 {The $z>1$ rate is a continuous extension that decreases with increasing redshift and is consistent with the large rate uncertainties in \citet{2014AJ....148...13R}  and  \citet{Strolger:2015kra}}

For the low-$z$ sample, we do not use measured rates; instead, we arbitrarily generate $4200$ SNe~Ia, with a redshift distribution following Eq.~\ref{eq:rate_lowz}, such that the number of selected events is roughly 10\% of the high-$z$ Ia sample.

Peculiar SNe are simulated following \citet{plasticc_K2019} and \citet{Vincenzi2021}. The rate of 91bg-like SNe is set to $\sim$12\% of the Ia rate with the same redshift dependence as Eq.~\ref{eq:rate_lowz}. The Iax rate follows the cosmic star-formation history \citep{Madau:2014bja}, normalized to $6\times 10^{-6}$~yr$^{-1}$~Mpc$^{-3}$ at $z=0$.

Core-collapse SNe rates are based on \citet{Strolger:2015kra} and \citet{Shivvers:2016bnc} (see also \citet{DES:2024hip, Vincenzi2021}). We define $\rcc$ as the green curve in Fig.~6 of \citet{Strolger:2015kra}, with Type II and Type Ibc rates taken as $0.7\times \rcc$ and $0.3\times \rcc$, respectively.


\begin{table}
\begin{center}
\caption{Average depth and time between observations.}
\begin{tabular}{ | c | c  c | c  c |} 
 \hline
%          & Depth         & Gap   &   Depth     & Gap\\ 
%          & $(5\sigma)$   & time  & $(5\sigma)$ & time\\ 
          & \multicolumn{2}{c|}{WFD} & \multicolumn{2}{c|}{DDF} \\
 Filter   & depth\tablenote{$5\sigma$ limiting magnitude.}   
          & gap\tablenote{Average time (days) between visits, excluding seasonal gaps.}         & depth   & gap \\   %[0.5ex] 
 \hline
    $u$ & $23.84$ & $10.5$  &  $25.05$ & $5.3$\\
 \hline
   $g$  & $24.80$ & $11.9$&  $25.52$ & $7.3$\\
 \hline
    $r$ &$24.21$ &$8.2$ & $25.60$ &$7.3$ \\
 \hline
   $i$ & $23.57$ & $8.6$&  $25.19$ & $7.3$\\
 \hline
    $z$&$22.65$  &$9.0$      &$24.79$  &$7.3$ \\
 \hline
    $Y$&$21.79$  &$11.2$  &$23.83$  &$7.4$ \\
%    \hline
%    Mean&$23.5$  &$8.4$  &$25.0$  &$7.0$ \\
 \hline
\end{tabular}
\label{t_lowz}
\end{center}
\end{table}    
    

\subsubsection*{Host Galaxies }
\label{sec:wgtmap}
 {We use the host galaxy libraries from \citet{lokken} for our study. Starting from the   {transient-}host  correlations  measured across ~$16,000$ galaxies in Pan-STARRS DR1 \citep{Gagliano:2020ucg, lokken}   {they} employed an approximate nearest neighbour algorithm to construct   {a} synthetic library containing   {a few} million CosmoDC2 galaxies that preserve these observed correlations \citep{korytov}.}
   {This procedure was applied to construct four separate host libraries ({\HOSTLIB}) for: 
SNIa,  II, Ibc and    {peculiar SNIa (Iax and 91bg)}.}
  To select the true host and possible host mismatch, all potential galaxies within a $10^{\prime\prime}$  radius are included for the $d_{DLR}$ analysis described in ~\ref{sec:host_gala_misassociation}. Roughly   half the galaxies have no neighbours within $10^{\prime\prime}$,  $\sim 1/3$  have  a single neighbour and  the remainder have 2 or more neighbours. 



  To describe the brightness profile of each galaxy, a single Sersic profile was estimated from \citet{Gagliano:2020ucg}.  The \HOSTLIB\ Sersic  parameters are assumed to represent the true galaxy profile, without instrumental smearing. A true SN location  near a host galaxy is overlaid using this profile. 
To incorporate PSF effects we smear with an average 
 FWHM  $\sim 1.0^{\prime\prime}$.
Thus there are two Sersic profiles: one for selecting SN location and the other for computation of $d_{DLR}$ in the analysis. 


 Since each  \HOSTLIB\ is weighted based on the transient type {,} the galaxy density is not accurate and therefore the neighbour fractions are not accurate. To have a host mismatch fraction that is similar to   \desfv\  \citep{DES:2023tfm}, we increase the Sersic size by  a factor of three. To avoid this artifact in future simulations,   an unweighted \HOSTLIB\ should be used eg. the diffsky\footnote{\URLDIFFSKY} catalogue from {\tt OpenUniverse2024}.  
 


% added by Alex M:
  {The photo-$z$ PDFs are true posteriors, derived jointly from the underlying redshift-photometry relation of the transient-specific host catalogs, rather than being the output of any specific photo-z estimators that come bundled with their own assumptions of priors/training sets that may be unrealistically well-matched to the test set relative to what we will have with real, non-simulated data.
This underlying redshift-photometry relation was modeled using a normalizing flow with the \texttt{PZFlow} package \citep{2024AJ....168...80C}.}  
  {For compact storage in the {\tt HOSTLIB},}
the photo-$z$   {PDF}   {is converted to}
%information is provided to the {\tt HOSTLIB} by defining 
$11$ quantiles 
%redshift columns 
corresponding to an integrated cumulative density function's probabilities of $0,\ 10\% ...100\%$.



  We define three metrics using the quantity $\Dz = |\zphot-\ztrue|/(1+\zphot)$: 
  \begin{enumerate}
      \item $\Dzbias=$ mean bias on $\Dz$,
      \item $\sigIQR=$ RMS  of the inter quartile distribution of $\Dz$ divided by 1.349,
      \item $\fout=$ the outlier fraction    of events satisfying $|\Dz| > 0.1$.
  \end{enumerate}
  {Our metric values are } $\Dzbias=0.001$,  $\sigIQR=0.127$ and $\fout=0.46$   {(see Fig.~\ref{fig:photoz_res}(a))}.

\subsubsection*{Weight Maps}
Unlike the {\tt PLaSTiCC} {\HOSTLIB} used in \citet{Mitra:2022ykq}, the use of \acro\ {\HOSTLIB}s enables modelling of the   host galaxy correlations. Here we weight each transient host as a  function of log of stellar mass ({\tt LOGMASS}) and the star formation rate ({\tt SFR}). The corresponding probability distribution  is called a  {\tt WGTMAP}. 
 {A {\tt WGTMAP}  for each transient class is provided by \citet{Gagliano:2020ucg} who used catalogue data from Pan-STARRS  \citep{ 2016arXiv161205560C}.} 


For each transient class, Fig.~\ref{Fig:host_sfr} shows the generated number of events as a function of  
  {observed}\footnote{  {Ideally, the true {\tt SFR} and {\tt LOGMASS} would be used in the simulation,
and observed values used in the analysis (Sec.~\ref{sec:anaysis}). Our simplification of using observed distributions in 
the {\tt WGTMAP} is expected to have a very small effect.} }  
{\tt SFR} and {\tt LOGMASS}  for the host galaxies selected by {\tt SNANA}. The green points  show a representative sample   from the  \acro\ \HOSTLIB, and the yellow  points show the host galaxies  selected by the {\tt WGTMAP}. For all transient classes, the {\tt WGTMAPS} preferentially select higher {\tt SFR}, and for the $91$bg {\tt WGTMAP} preferentially selects higher {\tt LOGMASS}. 



\begin{figure*}[htp]
    \centering
    \begin{minipage}[b]{0.29\textwidth}
        \centering
        \includegraphics[width=\linewidth]{Plot/scatter_density_curve_SNIa.png}
        \subcaption{SNIa}
    \end{minipage}
    \hspace{0.02\textwidth} % Adjust horizontal spacing
    \begin{minipage}[b]{0.29\textwidth}
        \centering
        \includegraphics[width=\linewidth]{Plot/scatter_density_curve_SNII.png}
        \subcaption{SNII}
    \end{minipage}
    \hspace{0.02\textwidth}
    \begin{minipage}[b]{0.335\textwidth}
        \centering
        \includegraphics[width=\linewidth]{Plot/scatter_density_curve_SNIbc.png}
        \subcaption{SNIbc}
    \end{minipage}

    \vspace{0.02\textwidth} % Adjust vertical spacing

    \begin{minipage}[b]{0.31\textwidth}
        \centering
        \includegraphics[width=\linewidth]{Plot/scatter_density_curve_SNIax.png}
        \subcaption{SNIax}
    \end{minipage}
    \hspace{0.02\textwidth}
    \begin{minipage}[b]{0.31\textwidth}
        \centering
        \includegraphics[width=\linewidth]{Plot/scatter_density_curve_SNIa91bg.png}
        \subcaption{SNIa91bg}
    \end{minipage}
\caption{Density of simulated galaxies as a function  {of} {\tt SFR}  {and} {\tt LOGMASS}. 
Each panel shows a different transient class. 
The yellow (green)  points correspond to \acro\ {\tt HOSTLIB} with (without) the {\tt WGTMAP}. 
To avoid saturation we show a random $10\%$ of the \HOSTLIB.  }
\label{Fig:host_sfr}
\end{figure*}



\begin{table}[ht!]
\centering
\caption{Summary of simulation statistics vs. sample and transient type.}
\renewcommand{\arraystretch}{1.2}
%\resizebox{\columnwidth}{!}{%
\begin{tabular}{|c|c|c|c|c|}
\hline
\multirow{3}{*}{\textbf{Sample}} & \multirow{3}{*}{\textbf{Transient}} & \multicolumn{3}{|c|}{\textbf{Number of events}}\\\cline{3-5}
 &  & \textbf{Total}  & \textbf{Passed} & \textbf{Passed} \\
 & & \textbf{generated}\footnotemark[1] & \textbf{trigger}\footnotemark[2] & \textbf{cuts}\footnotemark[3]\\
\hline\hline
Low-$z$  &  \multirow{2}{*}{Ia}  & \multirow{2}{*}{4200}  & \multirow{2}{*}{1133}   & \multirow{2}{*}{691 (0.165)}\\
 $<0.08$      &   &    &   & \\\hline\hline
\multirow{3}{*}{High-$z$} & All   & 233703 & 23479  & 5034 (0.022)\\
   & Ia & 35719  & 11393 & 4389 (0.13)\\
\multirow{3}{*}{ < 1.55} & Iax & 18909  & 579   & 141 (0.008)\\
   & 91bg & 4798   & 399   & $2\ (\sim 0)$\\
   & II & 121994 & 8643  & 313 (0.003)\\
   & Ibc & 52283  & 2465  & 189 (0.004)\\
\hline\hline
 \textbf{Final} & \multicolumn{4}{c|}{\textbf{5080 (Ia) + 645 (CC) = 5725}} \\
\hline\hline
\end{tabular}
%}
\label{T4}
\footnotetext[1]{ Total number of generated SNe}
\footnotetext[2]{Two or more detections separated by more than $30$ minutes. }
\footnotetext[3]{After selection cuts; numbers in `()' are fractions with respect to  ``Total generated'' column. }
\end{table}

\bigskip



\begin{figure*}[tb]
    \centering
    \makebox[1.\textwidth]{\includegraphics[width=1.0\textwidth]{Plot/SN_piechart_final.pdf}}
    \caption{ Pie chart showing the percentage of true SN  types.  
    \textbf{(a)} Generated  {physical rate}, 
    \textbf{(b)} Passes trigger and 
    \textbf{(c)} After analysis cuts ( last three columns of  Table~\ref{T4}). 
    The evolution of Ia fraction (orange) from left to right illustrates the increasing purity 
    with detection and analysis. }
    \label{fig:sn_pie}
\end{figure*}




% - - - - - - BiasCor SIM - - - - - -
\subsubsection*{Sample Sizes }
\label{subsec:ana_simbias}

A summary of average simulation statistics is shown in Table \ref{T4}.  For low-$z$ we generate $4200$ events, and  $\sim 1/4$ satisfy the $2$-detection trigger. For the high-$z$ sample we generate nearly $36,000$ SNIa, and  $\sim 1/3$ satisfy the trigger. The high-$z$ sample also includes  nearly $200,000$ non-SNIa, and  $\sim 6\%$ satisfy the trigger. Compared to the non-SNIa, the SNIa trigger efficiency is much higher because these events are significantly brighter. 

To visualize the relative number of events in the high-$z$ sample Fig.~\ref{fig:sn_pie} shows several pie-charts corresponding to  Table~\ref{T4}. 
Comparing physical rates, SNe~Ia corresponds to $15\%$ of all SNe (Fig.~\ref{fig:sn_pie}a), and after the trigger requirement, nearly $50\%$ of the SNe are type Ia (Fig.~\ref{fig:sn_pie}b). Examples of simulated light curves at different redshifts 
 are shown  in Fig.$1$ of \citet{Mitra:2022ykq}.   


To train our photometric classifier (Sec.~\ref{subsec:scone}) and to implement distance bias corrections in BBC   {to account for selection effects} (Sec.~\ref{subsec:ana_bbc}),  {we perform a separate biasCor simulation to generate a large sample of millions of  events.} The simulation procedure is 
 identical to that used for the simulated data. 
When simulated data are replaced with real data, the biasCor simulation is still needed.
 
% ==============================================================
% ==============================================================
\section{Analysis}
\label{sec:anaysis}

The SNIa-cosmology analysis steps are shown in Fig.~\ref{fig:analysis},
and described below. This analysis is similar to  M23, with additional complexities, from using  photometric classifier {\tt SCONE}  and the {\tt BEAMS} formalism.

While M23  performed  separate analyses for $\zphot$ and $\zspec$ samples, here we focus our study on  a single $\zphot$ analysis with no spectroscopic redshift or spectroscopic classification
 {in the high-$z$ subset}. 
While we anticipate some spectroscopic resources for LSST-DDF, here we make the most conservative assumption of no resources in order to be more sensitive to potential photo-$z$ related problems in the analysis.



\begin{figure*}[tb]
    \centering
    \makebox[1.\textwidth]
    {\includegraphics[width=1.0\textwidth]{Plot/DIA.pdf}}
    \caption{ Flowchart showing the cosmology analysis steps. }
    \label{fig:analysis}
\end{figure*}





% - - - - - LC fit - - - - -
\subsection{Lightcurve Fitting}
\label{subsec:ana_lcfitz} 

To standardize the brightness of SNIa, each light curve is fit
   {with} the same model used in the simulations: 
the SALT3 model in \citet{Pierel:2022pqc} that was trained using the formalism in \citet{salt3}.
This fitting process extracts the following parameters for each event:   time of peak brightness ($t_0$),  amplitude ($x_0$),  stretch ($x_1$),   color ($c$). For the low-$z$ subset we use a fixed redshift $\zspec$. For the high-$z$ subset we float  $\zphot$ as a fifth parameter using the host galaxy $\zphot$ pdf as a prior, and using the fitting formalism in \citet{Kessler2010}    {to measure the SN+Host photo-$z$.} 
SALT3 light curve fits for several events are shown by the smooth curves in Fig. $1$ of M23.

  
Following Eq.~1 in \citet{Kessler2010}, the 5-parameter SALT3 fit 
uses {\sc minuit} \citep{wfit2} to minimize the following $\chi^2$,
%

%


\begin{equation}
    \chi^2 =  \left( \sum_i \chisqflux \right ) + \chisqhost + \chisqsyst ~.
\label{eq:chi2tot}    
\end{equation}

$\chisqflux$  is the flux chi-squared for the $i^{th}$ observation, 
\begin{equation}
\chisqflux =  
       \frac{ \left [\FiData - \FiModel(\xvecSALT) \right ]^2}{\sigFi^2} + \logsigmodel  
    \label{eq:chiflux}
\end{equation}
 where 
 $\FiData$ is the flux for the $i^{th}$ observation,
 $\FiModel$ is the SALT3  model flux computed from the fitted parameters $\xvecSALT=\{t_0,x_0,x_1,c,\zphot\}$,
  {and $\sigFi$ is the quadrature sum of the flux and model uncertainties. 
 Since $\sigFi$ depends on fitted parameters, the 2nd term in eq.~\ref{eq:chiflux} is needed,}
 where $\sigFiTilde = \sigFi$ after  the first fit iteration so that
the $\logsigmodel$ term is close to zero in the second and third fit iteration. 


The contribution from the  host galaxy prior is 
\begin{equation}
     \chisqhost = -2\log[\probzhost]
\end{equation}
where $\probzhost$ is the host galaxy probability for the redshift value to be $\zphot$,
  {as determined from the {\tt PZFlow} PDF unpacked from the quantiles.}


The terms $\chisqflux$ and $\chisqhost$ are the conventional terms defined in \citet{Kessler2010} and used in many $\zphot$ studies including M23. Here we add a new term, $\chisqsyst$ to address a subtle artifact 
from fitting with systematics in Sec.~\ref{subsec:ana_cuts_syst}. 
%
The underlying artifact is that a small systematic perturbation (e.g., changing filter zeropoint) 
can lead to an unphysical large change in SALT3-fitted values that result in either rejecting the
event, or assigning such a large systematic uncertainty that the event is effectively ignored in a cosmology fit.
While this artifact occurs rarely for any one systematic,
the cumulative impact from many systematics can lead to a   significant loss   {(see sec.~\ref{sec:cov})}, 
and this loss is not currently modelled in our BBC method. 
To reduce this loss, we use the reference fitted parameters without a systematic shift as a prior 
in the SALT3 fit for each systematic shift,
\begin{equation}
    \chisqsyst = \left[ (\pk- \pkref ) / \sigkref \right]^2
    \label{eq:chisqsyst}
\end{equation}
where $\pkref$ are the fitted parameters ($k=1,5$) for the reference,
$\pk$ are for fits with a systematic perturbation, and $\sigkref$ is the fitted 
uncertainty from the reference fit.


After each {\sc minuit} fit iteration, the wavelength range for each LSST passband is 
transformed to the rest-frame using fitted $\zphot$.
If the central rest-frame wavelength $\meanlamb$ is outside of the 
SALT3 model range (2800-17000\AA),
the passband is dropped in the next fit iteration;
a previously dropped passband can be added if it is within the model wavelength range.
If any passband is dropped or added, the fit iteration is repeated to ensure that a consistent set of passbands are included in the fit.  

  The fitted photo-$z$ uncertainty is not explicitly used in our analysis. Compared to a spec-$z$  fit, the impact of fitting the photo-$z$  is  to increase the uncertainty on the SALT3 fitted parameters.  


% - - - - - cuts & systematics - - - - -
\subsection{Selection Requirements \& Systematic Uncertainties}
\label{subsec:ana_cuts_syst}


We apply the following selection requirements (cuts) based on previous analyses:
\begin{enumerate}
  \setlength\itemsep{0.1em}
   \item at least three bands with maximum signal to noise ratio SNR$>4$ 
   \item $|x_1| < 3.0$
   \item $|c| < 0.3$
   \item stretch uncertainty $\sigma_{x1} < 1.0$
   \item time of peak brightness uncertainty $\sigma_{t0} < 2.0$~days
   \item $\Pfit > 0.05$\footnote{$\Pfit$ is the SALT3 fit probability 
          computed from $\chi^2$ and number of degrees of freedom.}. No explicit cuts are applied on {\tt SCONE} classification. 
   \item $0.01 < \zphot < 1.4$ \label{cut_commonlast}
   \item Fitted bands satisfy $2800 <\meanlamb / (1+\zphot) < 17000 $ the valid range of the rest frame wavelengths for the SALT3 model; same cut is applied for  $\zphot\pm \sigma ({\zphot})$ \label{cut_fitband}
   \item valid bias correction for all systematic variations (see  Section~\ref{subsec:ana_bbc}). \label{cut_biascor}
\end{enumerate}

 





  

 The statistics after cuts are presented in the last column of Table~\ref{T4}, and the transient type fractions  (for high-$z$) are shown in  Fig.~\ref{fig:sn_pie}. The high-$z$ true-Ia fraction is $15\%$ for  physical rates (Fig.~\ref{fig:sn_pie}a),  increases to $\sim 50\%$  after two detection triggers (Fig.~\ref{fig:sn_pie}b), and increases again to $\sim 90\%$ after selection cuts (Fig.~\ref{fig:sn_pie}c). 

 
The  $\zphot$ residual vs. $\ztrue$ is shown in Fig.~\ref{fig:photoz_res}a
for all galaxies in the catalogue, 
and in Fig.~\ref{fig:photoz_res}b for host galaxies
after SN~Ia trigger and selection cuts. 
 {Fig.~\ref{fig:photoz_res}c and Fig.~\ref{fig:photoz_res}d  shows the performance for SN + Host combined fit, broken up into SNIa and non-Ia samples, that are used in the Hubble diagram.}  { 
Both the core resolution and the outlier fraction for Ia  are much  better ($\sim 5$x)  than
for the galaxies in Fig.~\ref{fig:photoz_res}a.} 


Figure~\ref{fig:z_dist} shows the redshift distribution in the final cosmology sample for  1 of the \NSAMPLE\ samples, used for the Hubble diagram.  
\begin{figure*}[!ht]
  \begin{minipage}{\textwidth}
    \centering
    \includegraphics[width=.4\textwidth]{Plot/photoz1_all.pdf}\quad
    \includegraphics[width=.4\textwidth]{Plot/photoz2.pdf}\\
    \includegraphics[width=.4\textwidth]{Plot/photoz_SNplusHOST_IaOnly.pdf}\quad
    \includegraphics[width=.4\textwidth]{Plot/photoz_SNplusHOST_non-Ia.pdf}
    \caption{Photo-$z$ residual ($\zphot-\ztrue$) vs. $\ztrue$ for
  (a) full host galaxy catalogue, 
  (b) subset of host galaxy catalogue after trigger and cuts,  
  (c) combined SALT3 fitted SN + HOST photo-$z$ for SNIa only and (d) same but for non-Ia. 
   The $\sigIQR$ and $\fout$ numbers on each panel are computed for $0.4<\ztrue < 1.4$.
  The source of photo-$z$ is indicated on each  {panel.}  } 
    \label{fig:photoz_res}
  \end{minipage}\\[1em]
\end{figure*}

\begin{figure}[h]
    \centering
    \includegraphics[width=\columnwidth]{Plot/z_dist.pdf}
    \caption{
     {True redshift distributions for 1 of the 25 samples:
    low-$z$ in dashed blue, high-$z$ true SNe~Ia in black, and high-$z$ true non-Ia in green.
    } 
    }
    \label{fig:z_dist}
\end{figure} 





The systematics considered in this analysis are   {the same as those described in Sec.~IV B of} \citet{Mitra:2022ykq}. Table~\ref{tb:syst_sources} summarizes the systematics used and their corresponding values.  The systematics include variations in galactic extinction, calibration  and host photo-$z$. 
 {The photo-$z$ systematic is a naive coherent shift of $0.01$ based on previous DES weak lensing cosmology  analyses \citep{PZ2021_DES3YR_WL},
but this systematic evaluation needs more effort based on how photo-$z$'s are trained for SN hosts, 
particularly at higher redshifts where the training set is limited. 


  {Missing systematics include  
  \begin{enumerate}
\item   astrophysical model of intrinsic scatter, which dominates the systematic error budget in \desfv ; 
\item   photometric classification (e.g, different classifiers and training sets),
  which contributes $\sim 1/3$  to  the total \DESVYR\ systematic 
  (see Table 7 in \citet{DES:2024hip}); 
\item   {photo-$z$ pdf variations that are more complex than a fixed shift; }
\item   {simulated redshift dependence of volumetric rate;  }
\item   {cosmology parameters used in simulated bias corrections.} 
  \end{enumerate}
}   

The first \ref{cut_commonlast} cuts  are applied only to the nominal analysis without systematics; for each systematic variation we do not apply these \ref{cut_commonlast} cuts but instead process  the list of events passing cuts from the nominal analysis. For example, if an event has fitted SALT3
color parameter $c=0.299$, and migrates to $c=0.3001$ for a calibration systematic, this event is preserved  without applying cuts that require $|c|<0.3$.  

Requiring valid fit bands (cut \ref{cut_fitband}) for each systematic results in a $ 3.3 \%$ loss from requiring  the same events passing all cuts.
Requiring a  valid biasCor (cut \ref{cut_biascor} and Sec.~\ref{sec:bbc}) for each systematic results in a $0.51\%$ loss.



\begin{table}
\caption{Source\rk{s} of Systematic Uncertainty } \vspace{-0.5cm}
\begin{center}
\begin{tabular}{ | l | l | l | c | } 
 \hline
Row  & Label & Decription & Value \footnote{Shift (or scale) applied to simulated data before each re-analysis} \\
 \hline\
 1 & StatOnly & no systematic shifts &  --- \\
 \hline\
 2 &  MWEBV & shift $E(B-V)$ & 5\%  \\ 
 \hline\
 3 &  CAL\_HST  & HST  calibration offset & $0.007\times\lambda$\\ \
  4 & CAL\_ZP & LSST zero point shift & $5$ mmag \\ \
5 & CAL\_WAVE & LSST Filter shift & $5$ \AA\\
  \hline\
6 & zSPEC  & shift $\zspec$ redshifts  &   $5\times 10^{-5}$\\ \
 7 & zPHOT  & shift $\zphot$ redshifts  &   $0.01$ \\ 
 %8 & zPHOTERR & scale host $\zphot$ uncertainty &  $1.2$\\ 
 \hline
\end{tabular}
\label{tb:syst_sources}
\end{center}
\end{table}
\subsubsection{Host Galaxy Association}
\label{sec:host_gala_misassociation}
Host galaxy matching is done via the Directional Light Radius (DLR) method 
  { \citep{Sullivan_2010, DES:2016aon}}. 
A normalized, dimensionless parameter \(d_{DLR}\) is defined as the ratio of the angular separation \(\Delta\theta\) between the supernova (SN) and DLR defined as the galaxy's centroid to the galaxy's effective radius in the direction of the SN. DLR is derived from the smeared Sersic profile.
Mathematically, this is expressed as:
\[
d_{DLR} = \frac{\Delta\theta}{\text{DLR}}.
\]
% 
 {The two galaxies with the smallest $d_{DLR}$ are tagged as  $\tt HOST1$ and $\tt HOST2$, where $\tt HOST1$ has the smaller $d_{DLR}$ if $\tt HOST2$ exists.}
A host  match is defined if  {$\tt HOST1$   has $d_{DLR}<4$}.  

 {Defining hostless events for a photo-$z$ analysis is ill-defined because
we cannot perform a photo-$z$ fit on hostless events. Nonetheless, we make a rough
estimate here by applying the SNR and sampling cuts in Sec.~\ref{sec:anaysis}, but no lightcurve fit.
The hostless fraction is   {$2\%$}, which is consistent with the 2\% hostless fraction for \DESVYR\ \citep{DES:2023tfm}.  
If we restrict our
true redshift to the \DESVYR\ range ($z<0.9$), our hostless fraction
stays at 2\% and is still consistent with \DESVYR. For events passing cuts, Table~\ref{tab:host_mis} shows the breakdown of events
with 1 or 2 matches with $\ddlr < 4$ and compares the host galaxy mismatch statistics of 
Host 1 with the numbers reported in \citet{DES:2023tfm}. }
  {Figure~\ref{fig:host_misac} 
illustrates the distributions of DLR and $r$-band magnitude
while showing the Table~\ref{tab:host_mis} nominal (second column) cuts.} 

\begin{table}[htbp]
\centering
\begin{tabular}{|c|c|c|c|}
\hline
 & LSST &LSST&    \\
Condition & SIM &SIM II\footnote{ $z<0.9$,  {to compare with DES} range. 
}
  & DES\footnote{see sec. 5.1 in \citet{Qu:2024fja}.}    \\
\hline
 {$>=1$} Host  {match}\footnote{ {  {atleast one match with DDLR$<4$; for the denominator, } zero matches include events with DDLR$>4$, but does not include
  host-less events}} & 98 &  {98} & 98   \\
$N_{\rm match}{=}2$ all\footnote{ {2 host matches with}  $\ddlr$<4} & 12 & 12  &9  \\ % both host 1 and host 2 DDLR < 4
$N_{\rm match}{=}2$ wrong\footnote{  Two host matches with $\ddlr <4$, and the galaxy with the smaller $\ddlr$ is the incorrect host. This value is not available for DES. } & 1.0 & 1.1 & ---  \\ 
Host depth  {($i$-band)} & 25.4 && 25.5   \\
\hline
No Host & 2 \% &2 \% & 2 \%   \\ 
\hline
\end{tabular}
\caption{ Statistics of host galaxy (mis)match  (in $\%$) in this analysis (second column).  In the third column we compared it with DES analysis (fourth column) from \citet{DES:2023tfm}.  
}
\label{tab:host_mis}
\end{table}

\begin{figure*}[ht!]
    \centering
    \begin{tikzpicture}
        % Place the figure inside TikZ
        \node [anchor=south west, inner sep=0] (image) at (0,0) {\includegraphics[width=\textwidth]{Plot/Host-Mis.pdf}};      
        % Overlay the watermark at the top right of the figure
        %\node at (image.north east) [opacity=0.2, anchor=north east, xshift=-10, yshift=-10] {\includegraphics[width=0.2\textwidth]{Plot/incomplete.jpg}};
    \end{tikzpicture}
    \caption{  Distributions illustrating
    host galaxy (mis)-association. 
 We find  
 $\sim 98\%$ host galaxy match and $\sim 2\%$ of no host matches.    $12\%$ `double match all',  is the aggregate for all matches for $\rm HOST\_\ddlr < 4$ and $11\%$ for the sum of all correct host matches, satisfying $\rm HOST\_\ddlr < 4$  condition. The grey vertical dashed line in the first figure marks $\ddlr = 4$. 
    }
    \label{fig:host_misac}
\end{figure*}


% - - - - - - Classifier - - - - - -
\subsection{Photometric Classification: SCONE}
\label{subsec:scone}
Deep learning based photometric classifiers have outperformed   {their} predecessors such as  light curve template matching algorithms \citep{p2} or feature extraction with either sequential cuts \citep{campbell} or machine learning methods \citep{Moller2016, Lochner2016, Dai2018}.   {\DESVYR\ used neural network based photometric classifiers, SCONE\footnote{Supernova Classification with a Convolutional Neural Network} \citep{scone} and SuperNNova \cite{Moller2020}}


We use SCONE  as our photometric classifier  {to determine the probability ($\Pia$) that each event is a SN~Ia.} SCONE   requires photometric data only, without the need for  accurate redshift  {estimation},    {and it } has relatively low computational and dataset size requirements for achieving high accuracy.    % These cuts are explained in details in \citet{Qu:2024fja}. 




For training the network with  $150$ iterations, we generate a simulation consisting of ${\sim}50,000$ Ia and ${\sim}50,000$ non-Ia after selection requirements. Even for real data, the procedure of training on simulated data is the same to ensure that 
the training data has the same redshift range and selection effects as the data.
The training process  took ${\sim}3$ hours. The network attained $98\%$ accuracy for both training and validation  {and a LOSS of 0.07}. Fig.~\ref{fig:scone_pred} shows $\Pia$\   distributions and they match our qualitative expectations; 
for true Ia's the probability is peaked at $p\sim1$ and for true non-Ia the peak is at $p\sim0$.


\begin{figure}[h]
    \centering
    \includegraphics[width=\columnwidth]{Plot/SCONE_PRED.pdf}
    \caption{Distribution of SCONE  {$\Pia$}
    for true SNIa (black) and for true non-Ia (dashed blue histogram). 
    }
    \label{fig:scone_pred}
\end{figure}

%  - - - - - - BBC - - - - - - - - - -
\subsection{BEAMS with Bias Corrections (BBC)}
\label{subsec:ana_bbc}
 
 BBC  {incorporates} BEAMS  { methodology and} produces a Hubble diagram that is corrected for selection effects and also corrected for non-Ia contamination. 

\subsubsection{ BEAMS}
\label{sec:beams}
The BEAMS framework was introduced by \citet[][]{kunz, 2012MNRAS.421..913N}, \citet{hlozek12, Knights:2012if}. 
 In the BEAMS framework, the total likelihood function includes two terms: one for the SN Ia population, $\mathcal{L}_\mathrm{Ia}$, and another for the contaminant population, $\mathcal{L}_\mathrm{CC}$ \footnote{The CC notation includes all non-Ia, but here we follow the KS17 notation.}, expressed as:
\begin{equation}
    \prod_{i=1}^{N_{\mathrm{SNe}}} \big( \mathcal{L}^i_\text{Ia}+\mathcal{L}^i_\text{CC} \big) .
    \label{eq:SUM}
\end{equation}
The likelihood components for each event $i$  are:
\begin{equation}
  \begin{aligned}
    \mathcal{L}^i_\text{Ia} &=  {\Pia}^i\times D_\text{Ia}(z_i, \mu_{\text{obs},i}, \mu_{\text{ref}}),\\
    \mathcal{L}^i_\text{CC} &= (1-P^i_{\text{Ia}}) \times D_{\text{CC}}(z_i, \mu_{\text{obs},i}, \mu_{\text{ref}}),
  \end{aligned}
  \label{eq:likelihoods}
\end{equation}
where $D_{\text{Ia}}$ and $D_{CC}$  are the SNIa and non-Ia Hubble residual likelihoods.   {A subtle caveat in Eq 8 is that the same light curve data are used for $\Pia$ and $D_{Ia}$, $D_{CC}$ but potential covariances are not accounted for.} The quantities in  Eq.~\ref{eq:likelihoods} can  be   used to determine a  BEAMS probability,
\begin{equation}
   \mathcal{P}_{\tt BEAMS} =  \frac{ {\Pia}^iD^i_{\text{Ia}}}{  P^i_{\text{Ia}}D^i_{\text{Ia}} + (1- P^i_{\text{Ia}})D^i_{\text{CC} }} \label{eq:PBIa},
\end{equation}
which quantifies the probability of a given SN belonging to the Ia population versus the contaminant population, taking into account both the classification probability ($P^i_{\text{Ia}}$) and the Hubble residual information. 

%----- BBC-------
\subsubsection{BBC}
\label{sec:bbc}
  BBC is a fitting framework that  accounts for selection effects using a simulation, incorporates BEAMS to account for non-Ia contamination, and implements the \citet{2011ApJ...740...72M} method of measuring distances independently of cosmological parameters. 
BBC reads the SALT3 fitted parameters (high-$z$ and low-$z$) 
from the data and biasCor simulation (sec.~\ref{subsec:ana_simbias}),
and produces a bias-corrected Hubble diagram, both unbinned and in redshift bins.
For each event, the measured distance modulus is based on the \citet{Trip1998} equation,
\begin{equation}
   \mu = m_B + \alpha x_1 - \beta c + \gamma G_{\rm host} + M_0 + \mubias~,   \label{eq:mu}
\end{equation}
where $m_B \equiv -2.5\log_{10}(x_0)$, 
$\{\alpha,\beta,M_0\}$ are global nuisance parameters, and
$\mubias = \mu - \mutrue$ is determined from the biasCor simulation
in a 4-dimensional space of $\{z,x_1,c,x_M\}$  {\citep{Popovic:2021cwq}},   {where $x_M$ is the host galaxy stellar mass}. A valid bias correction is required for each event, resulting in a  ${<}1 \%$ loss, usually for extreme values of $x_1,c,z$ where the biasCor statistics are very low. The term $\gamma G_{\rm host}$ captures residual dependency between SNIa luminosity and host galaxy magnitude. The distance uncertainty ($\sigmu$) is computed from Eq.~3 of KS17. 

 


There are two subtle issues concerning the use of $\zphot$ and its   {measured} uncertainty $\sigz$.
First, 
the calculated distance error from $\sigma_z$ ($\sigmu^z$ in Eq.~3 of KS17) 
is an overestimate because it does not account for the correlated color error that 
reduces the distance error.
By floating $\zphot$ in the SALT3 fit, redshift correlations propagate to the other SALT3 parameter uncertainties,
and therefore we set $\sigmu^z=0$,   {however the peculiar velocity uncertainty is included}.
%
The second issue concerns the $\mu_\mathrm{{bias}}$ computation, where
$\mutrue$ is computed at SALT3-fitted $\zphot$ rather than the true redshift. 

To avoid a dependence on cosmological parameters, the BBC fit is performed in \NZBIN\ logarithmically-spaced redshift bins.
The fitted parameters include the global nuisance parameters 
($\alpha,\beta,M_0$) 
and bias-corrected distances
in \NZBIN\ redshift bins. The unbinned Hubble diagram is obtained from Eq.~\ref{eq:mu} using
the fitted parameters, and the distance uncertainty is computed from Eq.~(10) in  \citet{Kessler:2023toi}. 


%----- Cov Mat-------
\subsubsection{Covariance Matrix}
\label{sec:cov}
The uncertainties in SN distance measurements are encapsulated in a covariance matrix, which quantifies both statistical and systematic uncertainties. This covariance matrix ($\mathbf{C}$) is 
\begin{equation}
    \mathbf{C} = \mathbf{C}_{\text{stat}}+\mathbf{C}_{\text{syst}}
    \label{eq:cov}
\end{equation}
where $\mathbf{C}_{\text{stat}}$  is the diagonal covariance representing measurement errors.
$\mathbf{C}_{\text{sys}}$  is  the systematic covariance matrix 
\citep{conley11},
\begin{equation}
    C^{ij}_{\text{syst}} = \sum_{S}\left(\Delta \mu_{\text{obs},S}^{i}\right)\left(\Delta \mu_{\text{obs},S}^{j}\right) ,
\end{equation}
where $i$  and $j$ are SN indices,  and  $\Delta\mu_{\text{obs},S}$ is the change in  SN distance moduli resulting from systematic variation $S$. 

Because of cuts~\ref{cut_fitband} and \ref{cut_biascor} in Sec.~\ref{subsec:ana_cuts_syst}, the SN sample for each systematic is slightly different and therefore $\mathbf{C}_{\text{sys}}$ cannot be computed. Following previous analyses we select the common events in all systematics and use this sample for the Hubble diagram and  $\mathbf{C}$.  For redshift systematics that result in migration to another redshift bin, the original (no syst) redshift bin is preserved for the BBC fit.  


 {Table~\ref{tb:bbc} summarizes the loss from the common-event cut for two different priors 
($1\sigma$ and $2\sigma$) in Eq.~\ref{eq:chisqsyst}. The  $1\sigma$ prior is used in our analysis, 
while the other prior is a diagnostic to see the impact of a looser prior. 
The final column is a comparison for a flat prior case, which is the standard used for previous analyses 
e.g., \DESVYR. Table~\ref{tb:bbc} breaks down this loss into three sources,}
 {and the dominant source of loss is from analysis cuts.}
This implicit common-event cut is not modelled in the bias correction simulation because we do not perform the CPU-intensive task of evaluating systematics for the large biasCor sample. 
This implicit cut rejects ${\sim}4\%$ of the sample ($2^{nd}$ column of  Table~\ref{tb:bbc}) and  could be a  potential source of unmodelled bias. 
 {In comparison a flat prior results ${\sim}12\%$ loss, illustrating that using a $1\sigma$ prior
significantly reduces this loss.}


\begin{table}[h]
\begin{center}
\caption{
Analysis loss from requiring valid biasCor and requiring that the same events pass cuts for all systematic variations.}
\begin{tabular}{  | c  c | c  c |c |  c| } 
 \hline
    
         Data cuts not applied &       & $\%$ Loss 
          &     & $\%$ Loss& $\%$ Loss\\  
         to biasCor sample &       &  $1\sigma$
          &     &$2\sigma$ &Flat\\   
          &       &   prior\footnote{ {This prior is used in our analysis;  see} $\chisqsyst$ in Eq.~\ref{eq:chisqsyst}.  }
          &     & prior\footnote{ {Additional diagnostic that is not used in the analysis.}}& prior\footnote{All previous analyses using {\tt SNANA}, including M23 and \DESVYR\ , used a flat prior.}\\          
 \hline
      Analysis cuts on systematics &  & $ 3.3\%$ & &$ 3.9\%$& $ 7.5\%$\\
      Valid biasCor in Nominal &  &   $0.44\%$ & &$0.89\%$& $0.86\%$\\
      Valid biasCor in systematics\footnote{for stat-only analysis without systematics, this loss is zero by definition.} &  &   $0.07\%$ & &$0.37\%$& $4.00\%$\\
      \hline
      &  &  &&& \\
     Total &  & $\sim$  $4\%$ & &$\sim 5.2\%$& $\sim 12.4\%$\\
 \hline
% \hline
\end{tabular}
\label{tb:bbc}
\end{center}
\end{table}


\subsubsection{Cosmology Fitting and Figure of Merit}
\label{subsec:ana_wfit}
In cosmology fitting, 
we minimize the following $\chi^{2}$:
\begin{equation}
    \chi^2 = \Delta\mu^{T}\mathbf{C}^{-1}\Delta\mu,
    \label{eq:chi}
\end{equation}
where $C$ is the covariance matrix (Eq.~\ref{eq:cov}), and $\Delta\mu$ is the vector of differences between observed and theoretical distance moduli.
\newcommand{\urlwfit}{\url{https://github.com/RickKessler/SNANA/blob/master/src/wfit.c}}
We use the same fast minimization code\footnote{\urlwfit} that was used in  \citet{Mitra:2022ykq}. We  approximate a
Cosmic Microwave Background (CMB) prior using the $R$-shift parameter (e.g., see Eq. 69
in \citet{komatsu}) computed from the same cosmological 
parameters that were used to generate the SNe~Ia. The $R$-uncertainty is
$\sigR=0.006$, tuned to approximate the constraining power of \citet{Planck2018}.
We fit with  \wCDM\ and \wwCDM\ (CPL \citep{cpl2}) models, where $w=[w_0+w_a(1-a)]$ and $a = 1/(1+z)$.


  {
Prior to 2021, most cosmology results were based on a redshift binned HD
to reduce computation time compared to an unbinned HD. After \citet{Binning_is_sinning}
showed that fitting an unbinned Hubble diagram gives optimal constraints,
cosmology results have been based on unbinned HDs (e.g.,  {Pantheon+ \citep{pantheon+}} and {\desfv} \citep{DES:2024hip}).
In anticipation of much larger HDs in the future that would require enormous 
memory and CPU resources for the optimal unbinned approach\footnote{For example, a future WFD+DDF HD with $\sim 10^5$ events results in a covariance matrix $\mathrm{C}$ with $\sim 10 $ billion elements, which makes the $\chi^2$ computation (Eq.~\ref{eq:chi}) very CPU intensive.}, 
\citet{Kessler:2023toi} proposed a rebinned approach to optimise constraints with a much smaller Hubble diagram, and hence use much less CPU resources. In our implementation, the rebinned analysis requires approximately four times less wall  time than the corresponding unbinned analysis.
}
%
We include rebinned results here with \nrebinx\ stretch bins, and \nrebinc\ colour bins. 
The maximum Hubble diagram size is $14 \times \nrebinx 
 \times \nrebinc = \HDSIZEREBIN$; after excluding empty bins the average fitted Hubble diagram size is  $\hdsizerebin$.  

 
 

For the \wwCDM\ model, the FoM is computed based on the dark energy task force 
(DETF) definition in  {\citet{Albrecht2006_FoM,Wang2008_FoM}},
\begin{equation}
   \text{FoM} \simeq \frac{1}{ \sigma (w_a) \sigma(w_0) \sqrt{1-\rho^2} }~,
    \label{eq:FoM}
\end{equation} 
%
where $\rho$ is the reduced covariance between $w_0$ and $w_a$. 

We fit the Hubble diagram and evaluate FoM for the three binning methods: binned,  unbinned,  and rebinned. M23 reported that unbinned results were pathological for unknown reasons and therefore not shown. Here we  {no longer } find  {pathological unbinned results, but we cannot trace what caused the problem in M23.} 
 

%  - - - - - - wfit - - - - -
 


%%%%%%%%%%%%%%%%%%%%%% 
\section{Cosmology Results  }
\label{sec:cosmology}



For one of the \NSAMPLE\ statistically independent samples, we show the Hubble diagram produced by the BBC fit,
both binned and unbinned, in Fig.~\ref{fig:HD}. The Hubble residuals with respect to the
true cosmology, $\Delta\mu = \mu - \mu_{\rm true}$, are consistent with zero.  
 {There is evidence of a redshift dependence as indicated by the $2.2\sigma$ significance in the 
slope   {and its uncertainty from }  a linear fit of Hubble residual vs redshift. 
  {Fitting a  slope to the combined }  $25$  data sets 
  {results in} a strong $6.5\sigma$  {significance in the slope, which leads to biased} 
cosmological parameters as shown below.
}



\begin{figure*}[!ht]
    \centering
    \includegraphics[width=\textwidth]{Plot/HD_phot.pdf}
    \caption{Redshift binned and unbinned Hubble diagram for one of the \NSAMPLE\ statistically independent samples.
   The lower panel shows Hubble residual $\Delta\mu$ with respect  {to} the true cosmology,  {and}
   the error bar shows the  {mean uncertainty}    in each BBC redshift bin (same redshift bins as in Fig.~\ref{fig:musig_compare}). 
   The  red-coloured events show the  {non-Ia} contamination. 
   The corresponding colour bar on the right shows the probability derived from BEAMS. 
   The number of events with $\mathcal{P}$(CC) $>0.75$ is $652$, 
   which is $\sim 11.4\%$ of the  {entire}  $\zphot$ data sample of $5725$ events.  
    }
    \label{fig:HD}
\end{figure*} 


\begin{figure}[h]
    \centering
    \includegraphics[width=\columnwidth]{Plot/sig_Mu_Compare_cc.pdf}
    \caption{
    Number of events (top) and BBC-fitted distance uncertainty (bottom) per 
     {BBC} redshift bin. 
    }
    \label{fig:musig_compare}
\end{figure} 


For BBC-fitted nuisance parameters we find $\alpha-\alphaTrueSym$ = $-0.00095\pm 0.00057$ (averaged over $25$ samples), we cannot however   {evaluate} the bias for $\beta$ since $\beta_{true}$ is 
 {not defined}
as described previously in Sec.~\ref{sec:snana}. We find $\sigint = 
 0.112 $  but there is no true $\sigint$ for comparison. Next, we show the BBC fitted distance uncertainties ($\sigmu$) in redshift bins 
(Fig.~\ref{fig:musig_compare}{b}). 
Owing to the larger statistics in high-$z$, the  {weighted} distance uncertainty 
 {per $0.1$ wide redshift bin} is in  {the} $0.01-0.02$ range at $z>0.6$. 

 

For the subsections below, we define $w$-bias to be   $w - w_{\rm true}$ where $w$ is from the \wCDM\  cosmology fit. 
A similar definition is used for $w_0$ and $w_a$ for the \wwCDM\ model.
{The bias uncertainty is the standard deviation of the  $w$-bias values divided  by $\sqrt{\NSAMPLE}$.}


% - - - - - - - - - - - - - - - - - - -     
\subsection{ \wCDM\  Results}
\label{sec:wcdm}
 {For the \wCDM\ cosmology fits, Table~\ref{tb:wcdm_results} summarizes the 
$w$-bias and associated uncertainties for the binned, unbinned, and rebinned analyses, 
each averaged over 25 realizations that include a CMB prior.} 
%
We find  {for all binning methods that there is a significant stat only bias of $\AVGwbias \sim 0.01$-$0.03$ with a significance of $\sim$ $3$-$8\sigma$ .}  {With systematics the biases are slightly smaller and the significance is in the $\sim$ $2$-$6\sigma$ range.} 


 {Overall, the mean $w$-uncertainty across all binning methods is  $\AVGwsig{\sim}0.020$ for the stat only case.
%
Including systematics the binned method results in $ \AVGwsig{\sim}0.025$, and the unbinned and rebinned methods result in a slightly smaller uncertainty of  $\AVGwsig\sim0.022$.  The larger binned uncertainty is expected.}
 Fig.~\ref{fig:wbias} shows the $w$ bias for all the three binning schemes.  
\begin{figure}[h]
    \centering
    \includegraphics[width=\columnwidth]{Plot/w_bias.pdf}
    \caption{ Bias in $w$CDM cosmology  {for each HD binning method}. 
    The grey dashed vertical line  {shows}  the  {true} reference $w=-1$.   }
    \label{fig:wbias}
\end{figure} 

% - - - - - - - - - - - - - - - - - - -     


\begin{table}[ht!]
\begin{center}
\caption{Summary of $w$CDM Cosmology bias for Binned, Unbinned, and Rebinned  {methods.}  }
\begin{tabular}{|c|cc|cc|}
\hline
\multirow{2}{*}{\textbf{Case}} 
 & \multicolumn{2}{c|}{\textbf{Stat only}} 
 & \multicolumn{2}{c|}{\textbf{Stat+Syst}} \\
\cline{2-5}
 & $\AVGwbias$\tablenote{Average bias among $\NSAMPLE$ samples with uncertainty of std/$\sqrt{\NSAMPLE}$} & $\AVGwsig$\tablenote{Average fitted uncertainty among $\NSAMPLE$ samples.}
 & $\AVGwbias$ & $\AVGwsig$ \\
\hline
Binned  & $0.0225 \pm 0.0037$ & $0.0200$  & $0.0200 \pm 0.0054$  & $0.0248$ \\
Unbinned  & $0.0126 \pm 0.0038$ & $0.0205$  & $0.0078 \pm 0.0039$  & $0.0222$ \\
Rebinned  & $0.0256 \pm 0.0032$ & $0.0198$  & $0.0199 \pm 0.0034$  & $0.0216$ \\
\hline
\end{tabular}
\label{tb:wcdm_results}
\end{center}
\end{table}




\subsection{ \wwCDM\ Results}
\label{sec:w0wa}

 {For the \wwCDM\ cosmology fits, Table~\ref{tb:w0wa_results} summarizes the average $w_0$ and $w_a$ biases and associated uncertainties for the binned, unbinned, and rebinned analyses, each averaged over 25 realizations that include a CMB prior. %The uncertainty on $\langle w_0\text{-bias}\rangle$ and $\langle w_a\text{-bias}\rangle$ is the rms of the 25 measurements divided by $\sqrt{25}$. 
The $\sigma_{w_0}$ and $\sigma_{w_a}$  {rows show}
the average fitted uncertainties across the 25 measurements. 
We find for all three binning methods a mild but measurable stat-only bias in $w_0$, ranging from $0.02$ to $0.03$ ($\sim$ $1$-$2$$\sigma$), and in $w_a$, ranging from $0.12$ to $0.24$ ($\sim$2-4$\sigma$). These biases are consistent in sign across the binned, unbinned, and rebinned analyses. 
When including systematics, the $w_0$-bias  slightly increases, varying between $0.02$ and $0.06$ ($\sim$1.5-3.5$\sigma$), while the $w_a$-bias increases modestly to $0.18$--$0.27$ ($\sim$3-4$\sigma$).}

 {For the binned analysis, stat+syst $\AVGFoM=$ \AVGFOMzphotsyst , is a $30\%$ degradation compared  to the stat-only  $\AVGFoM$  \AVGFOMzphotstat. The equivalent rebinned analysis goes  from $196$ to \AVGFOMzphotsystR , making $\sim19\%$ degradation.}  For the unbinned analysis  $\AVGFoM=$ degrades from $196$ to \AVGFOMzphotsystU, giving a degradation of $\sim14\%$. We suspect that the marginal improvement in the FoM for the rebinned analysis over the unbinned analysis arises from Python’s numerical precision handling and the treatment of the covariance matrix during inversion. 
 Therefore on average FoM values are $\sim\!190$ for the stat-only fits and $\sim\!150$ for the stat+syst case, indicating a typical degradation of about 25--30\% 
  {with} systematics. 
 The  {contours and bias}  for a single sample  {are}  shown in Fig.~\ref{fig:de_compare_stat} (stat only) and Fig.~\ref{fig:de_compare} (stat+syst). 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{table*}[ht!]
\begin{center}
\caption{Summary of  {bias, uncertainty and FoM for} $w_0w_a$CDM cosmology fits  {using}
binned, unbinned, and rebinned  {methods.}   }
\begin{tabular}{|l|c|c|c|c|c|c|}
\hline
& \multicolumn{3}{c|}{Stat only} & \multicolumn{3}{c|}{Stat+Syst} \\
\cline{2-7}
& Binned & Unbinned & Rebinned & Binned & Unbinned & Rebinned \\
\hline
{$\AVGwwbias$\tablenote{Average bias among $\NSAMPLE$ samples with uncertainty of std/$\sqrt{\NSAMPLE}$}}   & $-0.0221 \pm 0.0152$ & $-0.0175 \pm 0.0153$ & $-0.0314 \pm 0.0152$ & $-0.0242 \pm 0.0154$ & $-0.0573 \pm 0.0161$ & $-0.0307 \pm 0.0150$ \\
$\AVGwabias$   & $0.1874 \pm 0.0623$ & $0.1239 \pm 0.0631$ & $0.2424 \pm 0.0598$ & $0.1787 \pm 0.0646$ & $0.2724 \pm 0.0650$ & $0.2093 \pm 0.0579$ \\
{$\AVGwwsig$\tablenote{Average fitted uncertainty among $\NSAMPLE$ samples.}}   & $0.0632$ & $0.0636$ & $0.0641$ & $0.0767$ & $0.0732$ & $0.0724$ \\
$\AVGwasig$   & $0.2665$ & $0.2739$ & $0.2670$ & $0.3167$ & $0.3008$ & $0.2994$ \\
$\AVGFoM$   & $193$ & $182$ & $196$ & $131$ & $156$ & $159$ \\
\hline
\end{tabular}
\label{tb:w0wa_results}
\end{center}
\end{table*}



\begin{figure}[h]
    \centering
    \includegraphics[width=\columnwidth]{Plot/w0wa_bias_single.pdf}
    \caption{Similar to Fig.~\ref{fig:wbias} but for bias in $w_0 w_a$CDM cosmology. 
    Shown only for stat+systematic case. 
    }
    \label{fig:wabias}
\end{figure} 




\begin{figure*}[ht]
    \centering
    \begin{subfigure}[t]{0.48\textwidth}
        \centering
        \includegraphics[width=\columnwidth]{Plot/DE_Ellipse-stat_only.pdf}
        \caption{
             {For stat-only analysis,} \ww\  {68\% (dark blue) and 95\% (light blue) confidence }
            contours for a single SNIa data sample combined with CMB prior. 
            The crosses show the  bias  {and $\pm 1\sigma$ uncertainty} 
            from averaging results over the \NSAMPLE\ simulated data samples. 
        }
        \label{fig:de_compare_stat}
    \end{subfigure}
    \hfill
    \begin{subfigure}[t]{0.48\textwidth}
        \centering
        \includegraphics[width=\columnwidth]{Plot/DE_Ellipse-syst_unbin.pdf}
        \caption{
            Same as in left  {panel},  but for stat+syst case 
             {and unbinned HD.}
        }
        \label{fig:de_compare}
    \end{subfigure}
    \caption{Comparison of $2\sigma$ $(95\%$ confidence$)$ contours for two cases: (a) stat only, and (b) stat+syst.}
    \label{fig:de_compare_sidebyside}
\end{figure*}




Since there are many systematics contributing to the decrease in $\AVGFoM$,
we quantify the impact of each systematic ``$i$'' by 
recomputing the covariance matrix separately for each systematic  ($\mathbf{C}_{\text{syst, i}}$), 
and repeating the cosmology fit for each $\mathbf{C}_{\text{syst, i}}$.
We finally compute the FoM ratios
\begin{equation}
    \RatioFoM = \FoMSysti/\FoMStat~,
\end{equation}
where $\FoMSysti$ is the FoM from including only systematic $i$,
and $\FoMStat$ is the FoM without systematic uncertainties.
Note that $\RatioFoM \leq 1$. 

\begin{table}
\caption{FoM-Ratio $\RatioFoM$ for Each Systematic (\wwCDM\ model) 
 {and for each HD binning method.}
}
\begin{center}
\begin{tabular}{| l | c | c | c |}
 \hline
 Systematic(s) & \multicolumn{3}{c|}{$\RatioFoM$} \\
 \cline{2-4}
               & Binned & Unbinned & Rebinned \\ [0.5ex]
 \hline
None (stat only) & $1.00$ & $1.00$ & $1.00$ \\
zSHIFT & $0.95$ & $0.98$ & $0.98$ \\
%ZERRSCALE & $0.95$ & $0.99$ & $0.99$ \\
Photo\_shift & $0.95$ & $0.99$ & $0.99$ \\
MWEBV & $0.96$ & $0.99$ & $0.99$ \\
CAL\_WAVE & $0.90$ & $0.99$ & $0.95$ \\
Cal\_ZP & $0.80$ & $0.90$ & $0.89$ \\
CAL & $0.78$ & $0.88$ & $0.85$ \\
\hline\hline
Stat + All Syst  & $0.72$ & $0.85$ & $0.84$ \\
 \hline
\end{tabular}
\label{tb:syst_breakdown}
\end{center}
\end{table}

Table~\ref{tb:syst_breakdown}  {shows $\RatioFoM$ for all three binning methods.}
The FoM degradation is  dominated by the calibration systematics. 
Photo-$z$ systematics are not the dominant source of FoM degradation,  {in part because the}
anti-correlation between photo-$z$ and the SALT3 color parameter acts to self-correct distance errors. 
Similar findings were reported by \citet{Chen2022}. 
Thus, the overall impact of photo-$z$ systematics remains modest due to these combined effects.


\subsection{Comparison with \desfv \ and improvements with DESI }
\label{sec:comparison_with_desfv}


In Fig.~\ref{fig_w0wa_desi_des_planck_lsst_data_mock}, we provide constraints on  $w_{0}$ and $w_{a}$ from 
 {\desfv\ (yellow) and 3-year LSST-DDF (blue),}
and also in combination with CMB and baryonic acoustic oscillations (BAO) measurements. 
%
%
 {Compared to} the \desfv{} SNe only sample,  the LSST-DDF sample results in
 $\times16$ 
 improvement in the FoM for $w_{0}, w_{a}$. 
We also show the constraints when the SNe are combined with the CMB power spectra measurements from 
{\it Planck} \citep{rosenberg22} and the BAO measurements from the recent DESI-DR2 release \citep{desi_dr2_bao}. 
When combined with CMB and BAO, we find a $\times1.4$ improvement in the FoM 
when 
 {replacing}
DES  {with} 
the LSST sample \citep[Also see][]{raghunathan25}. 
The improvement in this case is less dramatic compared to the SN-only case because the 
 {constraining power is primarily}
from the different directions of degeneracies probed by CMB, BAO, and SN measurements. 



To combine the LSST-DDF SN sample with DESI-DR2 measurements, 
we impose a consistent cosmology by replacing the BAO data vectors with those
calculated using the \wwCDM\ cosmology with values set to {\it Planck} 2018 cosmology and $(w_{0}, w_{a}) = (-1, 0)$.
 {The original DESI DR2 covariances are used, 
and we do not  tweak the CMB measurements from {\it Planck}.}
%
\iffalse
\rkout{This is to avoid any differences in the actual cosmological parameters inferred from DESI-DR2, which are in mild tension with the standard $\Lambda {\rm CDM}$ \citep{desi_dr2_bao}, and the fiducial values used for the LSST SNe simulations. 
We do not make any such tweak for the CMB measurements from {\it Planck}. 
To this end, we replace the BAO data vectors from the actual DESI-DR2 distance measurements with the ones calculated using the $w_{0} w_{a} \Lambda {\rm CDM}$ cosmology with values set to {\it Planck} 2018 cosmology and $(w_{0}, w_{a}) = (-1, 0)$. 
The measurement covariance is same as the ones released by DESI DR2.}
\fi 


\begin{figure*}[h]
    \centering
    {\includegraphics[width=.8\textwidth, keepaspectratio]{Plot/w0wa_desi_des_planck_lsst_data_mock.pdf}}
    \caption{Comparison of constraints on the dark energy EoS parameters 
$w_{0}$ and $w_{a}$ 
 {using SNe~Ia from a simulated 3-year LSST-DDF sample from this work (light blue)}
and  {from} \DESVYR\ 
 {from \citet{DES:2024hip} (yellow)}.
The joint constraint combining DESI~DR2~BAO and {\it Planck}~CMB with 
\DESVYR\ SNe is shown in maroon,  {and the green lines show the best fit values.
The navy blue contour shows the corresponding joint constraint with LSST-DDF SNe replacing \desfv,
and  BOA  replaced with those calculated using $w_0,w_a = (-1,0).$}
%  % end rkout
    }
    \label{fig_w0wa_desi_des_planck_lsst_data_mock}
\end{figure*}







%clearpage 
\subsection{Bias Tests }
\label{sec:biases}
%\clearpage
%\subsection{Robustness Tests}
%\label{sec:robust_test}

The cosmology bias led us to implement 3 code improvements and perform 1 test.
Unfortunately this effort did not reduce the bias, but we nonetheless 
describe this effort below to inform future analyses.
%

 {\begin{itemize}
    \item In previous cosmology analyses using SNANA and pippin, a flat prior on each fitted parameter was used to run SALT3 light curve fitting on both nominal and  systematic variants. As described in sec.~\ref{subsec:ana_lcfitz}, requiring the same events in all systematics results in a $12\%$ loss. 
    Here we update the light curve fitting code to use a prior on systematic variants  
    (see $\chisqsyst$ term in eq.~\ref{eq:chi2tot} and ~\ref{eq:chisqsyst}), 
    which reduces this common-event loss to $4\%$. 
    %%%%%%%%%%%%%%%%%%%
    \item We found and fixed a subtle bug in evaluating the  $D_{CC}$ term in eq. \ref{eq:likelihoods}; 
    the distances were not bias corrected. This fix results in  {few milli-mag}
     changes in the distances. The small impact is not surprising because  the rms of the distance bias correction  is $\stdmubiasccprior$ mag, which is much smaller than the non-Ia Hubble residual rms of $\stddmuccprior$ mag. 
    % 
    \item The $D_{CC}$ term in eq.~\ref{eq:likelihoods} has  {previously}     been evaluated in redshift bins. However, we noticed that this term has significant dependence on colour and stretch and therefore we updated the BBC code to evaluate $D_{CC}$ on a $3D$ grid of $z,c,x_1$. 
        %%%%%
    \item As part of our bias testing, we found that the bias disappears if we 
      exclude all true non-SNIa events,   {which is equivalent to a perfect classifier}.
    Based on this clue, we ran a test requiring  $\Pia>0.01$ (sec.~\ref{subsec:scone}),
     {which} rejects $75\%$ of the non-SNIa. 
\end{itemize}}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%












%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}
\label{sec:conclude} 

 {}



This work represents one of the first comprehensive LSST-era cosmology 
 {analyses on a simulated high-$z$ sample}
based solely on photometric supernova data.
 {The analysis} combines photometric Type Ia classification with host-galaxy photometric redshifts. 

However, residual biases in recovered cosmological parameters indicate that further 
 {analysis development is needed.}
%
  {In our analysis,} systematics from   {astrophysical and } SN model training, especially limitations of the SALT3 model   {(sec.~\ref{subsec:ana_cuts_syst})},   {are not  included. These effects must be incorporated in future studies. Systematics arising from the photometric classifier (also not studied in this analysis) have been shown to contribute roughly one third of the \DESVYR\ systematic error budget.}
Statistical uncertainties therefore dominate over systematic contributions (sec.~\ref{sec:cosmology})   {in our analysis}.
%
%
Finally, the demonstrated efficiency of our rebinned HD approach, with  significant computational savings, highlights a scalable path for analysing the much larger LSST SN sample to come.


% {Discuss about missing systematics.}

 


\acknowledgments 
Author contributions are listed below. \\

A.~Mitra: co-lead project, SNANA simulations and analysis, writing. 
R.~Kessler: co-lead project, software, analysis, writing. 
R.~Chen: Internal review, and writing assistance.
A.~Gagliano: Analysis of host-galaxy catalogs, manuscript revisions.
M.~Grayling: Internal review, and writing assistance.
S.~More: Conceptualization.
G.~Narayan: Conceptualization and writing assistance.
H.~Qu: SCONE software support. 
S.~Raghunathan: Joint constraints using DESI-DR2. 
A.~I.~Malz: DESC Builder producing photo-$z$ prior information in mock data (data curation, writing).
M.~Lochner: DESC Builder for early contributions to supernova and significant optimization of observing strategy, manuscript revision.
\medskip


%
RK acknowledges pipeline scientist support from the LSST Dark Energy Science Collaboration. 
ML acknowledges support from the South African Radio Astronomy Observatory and the National Research Foundation (NRF) towards this research. Opinions expressed and conclusions arrived at, are those of the authors and are not necessarily to be attributed to the NRF.
This work was completed in part with resources provided by the University of Chicago’s Research Computing Center. This work is also supported by the National Science Foundation under Cooperative Agreement PHY-2019786 (The NSF AI Institute for Artificial Intelligence and Fundamental Interactions, http://iaifi.org/).

This paper has passed an internal review by the DESC  and we thank the DESC internal reviewers: Matthew Grayling and Rebecca Chen. 

The DESC acknowledges ongoing support from the Institut National de 
Physique Nucl\'eaire et de Physique des Particules in France; the 
Science \& Technology Facilities Council in the United Kingdom; and the
Department of Energy, the National Science Foundation, and the LSST 
Corporation in the United States.  DESC uses resources of the IN2P3 
Computing Center (CC-IN2P3--Lyon/Villeurbanne - France) funded by the 
Centre National de la Recherche Scientifique; the National Energy 
Research Scientific Computing Center, a DOE Office of Science User 
Facility supported by the Office of Science of the U.S.\ Department of
Energy under Contract No.\ DE-AC02-05CH11231; STFC DiRAC HPC Facilities, 
funded by UK BIS National E-infrastructure capital grants; and the UK 
particle physics grid, supported by the GridPP Collaboration.  This 
work was performed in part under DOE Contract DE-AC02-76SF00515.

\appendix
%\section{Comparison with DES}

\bibliography{main}
\bibliographystyle{mnras}%{apsrev}{aastex701}{plainnat}
\end{document}