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

Core-collapse supernovae (CCSN) are exciting astrophysical events linked to a broad range of physics, expanding from the explosions themselves to neutrinos and nuclear physics and the high density equation of state, including also gravitational waves and nucleosynthesis. Depending on the aspect in focus, different kinds of hydrodynamical simulations are required. Ideally, one would perform 3D, fully general relativistic  magneto-hydrodynamic simulations with accurate neutrino transport for many progenitor stars, rotation rates, and magnetic field configurations and strengths, and follow the evolution during several seconds after the explosion. This is clearly computationally impossible today, therefore some aspects may be sacrificed to gain insights \citep[see e.g.,][]{Janka:2012, Kotake:2012, Mueller:2019_2, Burrows:2021}. Here, we focus on the long-time evolution of several seconds after the explosion and investigate the dynamics, accretion, explosion energy, and approximate nucleosynthesis and their dependency on neutrino heating, rotation, and progenitor. Therefore, we perform 2D Newtonian simulations with a simple neutrino treatment, namely gray leakage. This allows us to describe and understand important trends occurring during the first seconds after the explosions. However, our simplifications prevent us to make solid quantitative conclusions.

The long-time evolution has been considered in many nucleosynthesis studies of supernova yields. Originally spherically symmetric explosions were artificially triggered by pistons \citep{Woosley:1995} or thermal-energy bombs \citep{Thielemann:1996, Nagataki:1998,Nomoto:2006}. In the last decades, there have been numerous studies based on spherically symmetric simulations with enhanced neutrino energy deposition
\cite[e.g.,][]{Scheck:2006, Arcones:2007,OConnor:2010,Suwa:2011, Ugliano:2012, Perego:2015, Pejcha:2015,Ertl:2016,Sukhbold:2016,Mueller:2016_1, Curtis:2019,Ebinger:2019,Ebinger:2020,Couch:2020}. However, there are important multidimensional effects such as convection that impact the nucleosynthesis, explosion energy, explosion morphology \citep[see e.g.,][ and references therein]{Janka:2012,Kotake:2012, Mueller:2019_2, Burrows:2021}, and neutrino-driven wind \citep{Arcones:2011}. Recently, it has been possible to perform 2D simulations with good neutrino transport, to follow the evolution for a few seconds after the explosion and to study the supernova nucleosynthesis \cite[see e.g., ][]{Harris:2017, Eichler:2018, Wanajo:2018, Sieverding:2020,Reichert:2021}. However, these studies cover a reduced number of models and often only for short times after explosion.  Models with longer times usually require a mapping to a larger grid with simpler input physics after shock revival \citep{Wongwathanarat:2015,Mueller:2018,Stockinger:2020,Bollig:2020}

Rotation is an additional and important ingredient that affects not only the explosion but also the long-time evolution. The impact of rotation on the explosion has been extensively studied based on 2D \citep{LeBlanc:1970,Fryer:2000,Kotake:2003} \citep[see. e.g, ][and references therein]{Marek:2009, Suwa:2010,Blondin:2017, Summa:2018,Vartanyan:2018} and 3D simulations \citep{Kuroda:2014,Nakamura:2014, Moesta:2014,Takiwaki:2016,Blondin:2017,Summa:2018}. In this paper, we also show how rotation impacts  the evolution during the first seconds after the explosion, by affecting the shock morphology and the mass accretion onto the proto-neutron star (PNS).

This paper is organised as follows. In Sect.~\ref{sec:method}, we present our simulation setup, initial conditions, and nucleosythesis network. An overview of the models and their evolution during the first second after bounce is discussed in Sect.~\ref{sec:exp}. We investigate the long-time evolution in Sect.~\ref{sec:longtime} including the evolution of the explosion energy (Sect.~\ref{sec:Eexp}), the impact of rotation and shock deformation (Sect.~\ref{sec:long-time_dshock}), and the evolution of accretion and neutrino-driven wind (Sect.~\ref{sec:long-time_observables}). An estimate of the nucleosynthesis is given in Sect.~\ref{sec:nucleosynthesis}. We conclude in Sect.~\ref{sec:conclusions}.

\section{Simulations and nucleosynthesis}
\label{sec:method}

\subsection{Hydrodynamics and neutrinos}
\label{sec:hydro_neut}

We employ the multiphysics FLASH code \citep{Fryxell:2000,Dubey:2009} to perform 2D (cylindrical geometry), CCSN simulations. The domain size for all simulations is $3.2 \cdot 10^{10}$~cm along the cylindrical axis and $1.6 \cdot 10^{10}$~cm perpendicular to it. Using the adaptive mesh refinement in FLASH \citep{MacNeice:2000}, we achieve a maximum resolution of $\approx 488$~m. Our setup is similar to previous core-collapse supernova studies with FLASH, see e.g., \cite{Couch:2014,Couch:2015_1}. Self-gravity is calculated with a multipole approximation \citep{Couch:2013_2} of Poisson's equation, without the effective general relativistic potential that was used in \citet{OConnor:2018_1}.

Before bounce, the deleptonization scheme of \cite{Liebendoerfer:2005_1} is used. After bounce, neutrinos are described by a gray, three-flavor neutrino leakage scheme with a ray-by-ray approximation, as in \cite{OConnor:2010, Couch:2014}. In order to facilitate  the explosions,  the neutrino heating  is artificially enhanced by a factor, $f_\mathrm{heat}$, that was introduced in previous studies \citep{OConnor:2010, Couch:2014, Couch:2015_1}. Here, we adjust slightly the implementation  to have the enhanced heating only in the gain layer and until the shock has reached a radius of $1000$~km.

Leakage schemes are not as accurate as the more sophisticated neutrino transport methods  \cite[see e.g., ][]{Just:2015, OConnor:2015, Kuroda:2016, Pan:2019,Glas:2019_1}. However, broad studies of multi-dimensional, long-time simulations up to $10$~s after bounce are yet not feasible with state-of-the-art transport methods, without sacrificing resolution. A consequence of the more approximate leakage scheme is an incorrect electron fraction, typically leaning towards more neutron-rich values \citep{Dessart:2009, Pan:2019}. We  discuss corrections to the electron fraction for nucleosynthesis calculations in Sec.~\ref{sec:nucleosynthesis}.

We use the LS220 equation of state (EOS) of \cite{Lattimer:1991} for the high density regions. In the outer layers of the progenitor and at late times, the density reaches very low values ($\rho \lesssim 10~\mathrm{g}\,\mathrm{cm}^{-3}$), therefore, we use a hybrid EOS approach with a transition to the Helmholtz EOS \citep{Timmes:1999, Timmes:2000, Fryxell:2000} at around  $5-1 \cdot 10^{5}~\mathrm{g}~\mathrm{cm}^{-3}$.

\subsection{Progenitor models and rotation}
\label{sec:pog_rot}

In this study, we use three different progenitors with ZAMS mass 11.2~$M_\odot$, 15~$M_\odot$, and 27~$M_\odot$  namely the s11.2, s15.0 and s27.0 models from \cite{Woosley:2002}. All three progenitor models are non-rotating. As we aim to explore the effect of rotation  on the explosion phase and on the long-time evolution, we superimpose a parametric rotational profile of the form
\begin{equation}
 \label{eq:rotation_parametric}
 \Omega(r) \, = \, \Omega_0 \cdot \frac{1}{1+(r/r_A)^2} \, .
\end{equation}
We set the characteristic radius $r_A$ to a fixed value of $3000$~km \citep[approximately the extent of the Fe and Si core in the s15 progenitor, see][]{Mueller:2004,Buras:2006} and use $\Omega_0$ as a free parameter to adjust the rotation strength. In addition to a non-rotating case ($\Omega_0 = 0$), we use six different rotation strengths:  $\Omega_0$~=~(0.01, 0.03, 0.06, 0.10, 0.20, 0.30)~$\cdot~(2\pi~\mathrm{rad \,s}^{-1})$. The three first values correspond to moderately rotating models \citep[compare to, e.g.,][]{Heger:2005,Ott:2006,Burrows:2007_2} and the three last ones to rapidly rotating.

We label our models with their corresponding progenitor, heating factor, and rotation strength. For example, the label s15\_F130\_R006 refers to a simulation of the s15.0 progenitor, with $f_\mathrm{heat} = 1.3$, and $\Omega_0 = 0.06 \cdot 2\pi$~rad/s.

\subsection{Tracer particles and nucleosynthesis network}
\label{sec:tracer_method}

In order to estimate the nucleosynthesis, we need the Lagrangian evolution of the matter ejected. Since the FLASH code is Eulerian, we use a tracer particle scheme  to track individual fluid elements (also called tracer particles, trajectories, mass elements). The tracers are not included in the simulation but calculated from the output (see e.g., \citet{Harris:2017} for detailed discussion). All tracer particles have equal mass and are distributed proportional to the density \cite[see][for tracers in FLASH code]{Dubey:2012}. For  the s15.0 progenitor, we initialize 21,281 particles at the beginning of the simulation. This corresponds to a mass of $1.76~\cdot~10^{-4}~M_\odot$ per tracer that is similar to the ``medium resolution'' case in \citet{Nishimura:2015}.  The amount of matter ejected by the tracers agrees within 1\% with the unbound material obtained directly from the hydrodynamics.

In addition to the evolution of density and temperature along the tracers, we also need the electron fraction. However, this is poorly determined in our simulations because we use a gray neutrino leakage scheme. Therefore, we correct the electron fraction to  estimate the nucleosynthesis within some uncertainties (see Sect.~\ref{sec:nucleosynthesis} for more details).

For the nucleosynthesis calculations, we use the nuclear reaction network code WinNet \citep{Winteler:2011,Winteler:2012}. The reaction rates include the JINA Reaclib compilation \citep{Cyburt:2010}, theoretical weak reactions from \citet{Langanke:2001}, and neutrino reactions from \citep[see also \citealt{Froehlich:2006b} for details about the neutrino reactions]{Langanke:2001b}. For all nucleosynthesis calculations, we evolve the electron fraction in nuclear statistical equilibrium (NSE) at a temperature of 20~GK and assume NSE down to 6.5~GK.

\section{Overview of models and explosions}
\label{sec:exp}
We present here all our models including 32 successful explosions (see Table~\ref{tab:model_overview} for an overview). We also include unsuccessful explosions in the table for completeness. These models were continued for 1-2~s after bounce but did not show signs of shock revival until then. First, we focus shortly on the explosion phase and discuss the long-time evolution in more detail later  (Sect.~\ref{sec:longtime}). We are aware that our results are only qualitative because the neutrino transport is approximated and we artificially increase the energy deposited by neutrinos to trigger explosions. Moreover, our simulations are 2D and this favours some artificial features and may amplify some of our conclusions, like the persistence of downflows. However, we are able to explore the long-time evolution for a big sample of models and we find trends and correlations that are robust and are present also in 3D simulations with better neutrino treatment \citep[e.g., see][]{Mueller:2015_2,Janka:2016,Summa:2018,Burrows:2019,Mueller:2019_1,Vartanyan:2019_1,Stockinger:2020,Bollig:2020}.

\begin{table*}
\caption{Overview of models}
\label{tab:model_overview}
\centering
\begin{threeparttable}
\begin{tabular}{cccccccccccc}
    name & prog. & $f_\mathrm{heat}$ & $\Omega_0$  & $t_\mathrm{end}$ & $t_\mathrm{exp}$\tnote{a} & $d_\mathrm{shock}$\tnote{b} & $E_\mathrm{exp}^\mathrm{100ms}$\tnote{  c} & $E_\mathrm{exp}^\mathrm{final}$ & $M_\mathrm{ej}$\tnote{d} & $M_\mathrm{PNS}$ & $t_\mathrm{wind}$\tnote{e} \\
    &&&[2$\pi$ rad/s]&[s] &[s] & & [B] & [B] & [$M_\odot$] & [$M_\odot$] & [s] \\
    \hline
    \hline
s15\_F100\_R000 & s15.0 & 1.00 & 0.00 & 1.55 & ... & ... & ... & ... & ... & 1.98 & ... \\
s15\_F105\_R000 & s15.0 & 1.05 & 0.00 & 1.54 & ... & ... & ... & ... & ... & 1.98 & ... \\
s15\_F105\_R003 & s15.0 & 1.05 & 0.03 & 1.55 & ... & ... & ... & ... & ... & 1.98 & ... \\
s15\_F105\_R020 & s15.0 & 1.05 & 0.20 & 1.45 & ... & ... & ... & ... & ... & 1.96 & ... \\
s15\_F108\_R000 & s15.0 & 1.08 & 0.00 & 8.96 & 1.07 & 1.26 & 0.04 & 1.16 & 0.76 & 2.16 & 0.00 \\
s15\_F110\_R000 & s15.0 & 1.10 & 0.00 & 8.14 & 0.93 & 1.06 & 0.04 & 1.41 & 0.90 & 2.10 & 0.00 \\
s15\_F110\_R001 & s15.0 & 1.10 & 0.01 & 1.63 & ... & ... & ... & ... & ... & 1.99 & ... \\
s15\_F110\_R003 & s15.0 & 1.10 & 0.03 & 1.54 & ... & ... & ... & ... & ... & 1.98 & ... \\
s15\_F110\_R020 & s15.0 & 1.10 & 0.20 & 1.46 & ... & ... & ... & ... & ... & 1.96 & ... \\
s15\_F115\_R000 & s15.0 & 1.15 & 0.00 & 9.25 & 0.92 & 1.01 & 0.04 & 1.34 & 0.93 & 2.20 & 0.00 \\
s15\_F115\_R003 & s15.0 & 1.15 & 0.03 & 1.61 & ... & ... & ... & ... & ... & 1.99 & ... \\
s15\_F120\_R000 & s15.0 & 1.20 & 0.00 & 7.89 & 0.28 & 0.94 & 0.14 & 1.21 & 0.86 & 2.03 & 0.00 \\
s15\_F120\_R001 & s15.0 & 1.20 & 0.01 & 6.74 & 0.30 & 0.87 & 0.19 & 1.20 & 0.76 & 1.95 & 0.00 \\
s15\_F120\_R003 & s15.0 & 1.20 & 0.03 & 9.31 & 0.36 & 0.29 & 0.15 & 0.70 & 1.28 & 1.80 & 2.44 \\
s15\_F120\_R006 & s15.0 & 1.20 & 0.06 & 6.03 & 0.38 & 0.77 & 0.11 & 0.63 & 0.88 & 1.86 & 0.00 \\
s15\_F120\_R010 & s15.0 & 1.20 & 0.10 & 1.61 & ... & ... & ... & ... & ... & 1.99 & ... \\
s15\_F120\_R020 & s15.0 & 1.20 & 0.20 & 1.46 & ... & ... & ... & ... & ... & 1.96 & ... \\
s15\_F130\_R000 & s15.0 & 1.30 & 0.00 & 9.26 & 0.23 & 0.59 & 0.15 & 0.79 & 1.42 & 1.75 & 1.73 \\
s15\_F130\_R001 & s15.0 & 1.30 & 0.01 & 4.60 & 0.22 & 0.37 & 0.27 & 1.00 & 0.82 & 1.85 & 0.00 \\
s15\_F130\_R003 & s15.0 & 1.30 & 0.03 & 8.63 & 0.23 & 0.37 & 0.22 & 0.75 & 1.52 & 1.68 & 3.48 \\
s15\_F130\_R006 & s15.0 & 1.30 & 0.06 & 4.17 & 0.24 & 0.25 & 0.28 & 0.73 & 0.59 & 1.87 & 0.00 \\
s15\_F130\_R010 & s15.0 & 1.30 & 0.10 & 9.08 & 0.28 & -0.00 & 0.24 & 0.64 & 1.40 & 1.72 & 5.10 \\
s15\_F130\_R020 & s15.0 & 1.30 & 0.20 & 9.25 & 0.32 & -0.13 & 0.19 & 0.77 & 1.34 & 1.75 & 2.46 \\
s15\_F130\_R030 & s15.0 & 1.30 & 0.30 & 9.30 & 0.36 & 0.20 & 0.15 & 0.67 & 1.26 & 1.76 & 3.60 \\
s15\_F150\_R000 & s15.0 & 1.50 & 0.00 & 7.96 & 0.17 & 0.64 & 0.30 & 1.21 & 1.13 & 1.87 & 0.00 \\
s15\_F150\_R001 & s15.0 & 1.50 & 0.01 & 7.69 & 0.17 & 0.68 & 0.36 & 1.12 & 1.31 & 1.71 & 0.16 \\
s15\_F150\_R003 & s15.0 & 1.50 & 0.03 & 9.43 & 0.21 & 0.24 & 0.40 & 1.09 & 1.59 & 1.69 & 3.37 \\
s15\_F150\_R006 & s15.0 & 1.50 & 0.06 & 9.03 & 0.17 & 0.19 & 0.36 & 0.91 & 1.56 & 1.69 & 5.34 \\
s15\_F150\_R010 & s15.0 & 1.50 & 0.10 & 8.00 & 0.21 & 0.22 & 0.31 & 0.73 & 1.55 & 1.69 & 4.81 \\
s15\_F150\_R020 & s15.0 & 1.50 & 0.20 & 8.94 & 0.22 & -0.21 & 0.37 & 0.83 & 1.57 & 1.67 & 4.81 \\
s15\_F150\_R030 & s15.0 & 1.50 & 0.30 & 9.40 & 0.21 & 0.18 & 0.16 & 0.81 & 1.42 & 1.69 & 3.80 \\
\hline
s11\_F100\_R000 & s11.2 & 1.00 & 0.00 & 4.36 & 0.16 & 0.36 & 0.04 & 0.36 & 0.41 & 1.33 & 0.00 \\
s11\_F100\_R003 & s11.2 & 1.00 & 0.03 & 4.43 & 0.25 & 0.19 & 0.04 & 0.44 & 0.33 & 1.39 & 0.00 \\
s11\_F100\_R020 & s11.2 & 1.00 & 0.20 & 4.63 & 0.25 & 0.61 & 0.02 & 0.64 & 0.31 & 1.40 & 0.00 \\
\hline
s27\_F105\_R000 & s27.0 & 1.05 & 0.00 & 3.79 & 1.42 & 0.93 & 0.05 & 0.67 & 0.27 & 2.22 & 0.00 \\
s27\_F105\_R003 & s27.0 & 1.05 & 0.03 & 1.82 & ... & ... & ... & ... & ... & 2.07 & ... \\
s27\_F105\_R020 & s27.0 & 1.05 & 0.20 & 1.73 & ... & ... & ... & ... & ... & 2.05 & ... \\
s27\_F110\_R000 & s27.0 & 1.10 & 0.00 & 3.06 & 0.31 & 0.63 & 0.07 & 0.33 & 0.17 & 1.93 & 0.00 \\
s27\_F110\_R003 & s27.0 & 1.10 & 0.03 & 1.30 & 0.54 & 0.39 & 0.05 & 0.23 & 0.09 & 1.87 & 0.00 \\
s27\_F110\_R020 & s27.0 & 1.10 & 0.20 & 2.06 & ... & ... & ... & ... & ... & 2.12 & ... \\
s27\_F120\_R000 & s27.0 & 1.20 & 0.00 & 2.89 & 0.29 & 0.72 & 0.07 & 0.65 & 0.39 & 2.02 & 0.00 \\
s27\_F120\_R003 & s27.0 & 1.20 & 0.03 & 1.79 & 0.34 & 0.59 & 0.10 & 0.55 & 0.19 & 1.91 & 0.00 \\
s27\_F120\_R020 & s27.0 & 1.20 & 0.20 & 2.11 & ... & ... & ... & ... & ... & 2.13 & ... \\
s27\_F130\_R000 & s27.0 & 1.30 & 0.00 & 3.47 & 0.21 & 0.73 & 0.07 & 1.24 & 0.55 & 2.09 & 0.00 \\
s27\_F130\_R003 & s27.0 & 1.30 & 0.03 & 2.85 & 0.20 & 0.55 & 0.10 & 0.90 & 0.38 & 2.00 & 0.00 \\
s27\_F130\_R020 & s27.0 & 1.30 & 0.20 & 3.82 & 0.29 & 0.16 & 0.07 & 0.46 & 0.23 & 1.92 & 0.00 \\
\hline
\end{tabular}
 \begin{tablenotes}
  \item[a] Post-bounce time when the maximum shock radius exceeds $600$~km
  \item[b] Explosion energy at $t_\mathrm{exp} + 100$~ms
  \item[c] Shock deformation parameter (Eq.~\ref{eq:dshock}) at $t_\mathrm{exp}$
  \item[d] Mass of ejected (unbound) matter at the end of the simulation
  \item[e] Total duration of neutrino-driven wind phases ($\dot{M}_\mathrm{acc,500km} = 0$)
 \end{tablenotes}
\end{threeparttable}
\end{table*}

The evolution of the shock radius is shown in Fig.~\ref{fig:rshk} for the three progenitors and for three rotation rates. The lightest progenitor, s11.2, explodes without additional neutrino heating within $0.25$~s for all variations of the rotation. Therefore, we only use $f_\mathrm{heat} = 1$ for this progenitor. For the heavier progenitors, the threshold for shock revival in non-rotating models is $f_\mathrm{heat} = 1.08$ and $1.05$ for s15.0 and s27.0, respectively. This is in agreement with estimates from progenitor compactness \citep{OConnor:2011}, which has values of $\xi_\mathrm{1.75} = 0.07$, $0.54$ and $0.53$ for the s11.2, s15.0 and s27.0 models, respectively \citep{Pan:2016,Summa:2016}. The explosion time varies depending on the progenitor and the heating factor and it is also affected by rotation.

\begin{figure}[ht]
 \centering
 \includegraphics[width=0.5\linewidth]{rshkMax_lin_s11_1.pdf}%
 \includegraphics[width=0.5\linewidth]{rshkMax_lin_s11_2.pdf}\\
  \includegraphics[width=0.5\linewidth]{rshkMax_lin_s15_1.pdf}%
 \includegraphics[width=0.5\linewidth]{rshkMax_lin_s15_2.pdf}\\
  \includegraphics[width=0.5\linewidth]{rshkMax_lin_s27_1.pdf}%
 \includegraphics[width=0.5\linewidth]{rshkMax_lin_s27_2.pdf}
 \caption{Early (left column) and long-time (right column) evolution of the maximum shock radii for  the three progenitors  (s11.2, s15.0, and s27.0 in top, middle, and bottom panels, respectively). Different rotation rates are indicated by different line style  an various heating factor by colours, as given in the figures. The dashed and dotted purple lines represent the models s15\_F120\_R001 and s15\_F150\_R030, respectively.}
 \label{fig:rshk}
\end{figure}

In the presence of rotation, we notice that shock revival generally occurs later, or fails overall \citep[see also, e.g., ][]{Fryer:2000,Thompson:2004,Summa:2018,Obergaulinger:2020}. This can be understood from a reduction of mass accretion through the stalled shock because of centrifugal forces. This leads to a diminished accretion luminosity and ultimately takes away support for the stalled shock.  Moreover, rotation contributes to shift matter from the poles towards the equatorial plane and this has an impact on the  neutrino luminosities in different directions. Figure~\ref{fig:rotation_angular_lumAll} indicates that, in the direction of the poles, the luminosities for rotating models are considerably smaller, in accordance with the lower density and mass accretion in this region. Notice that in 3D simulations and with magnetic fields, the smaller mass accretion rates and luminosities are often found in the equatorial plane, but the integrated luminosities are also smaller in the context of rapid rotation \citep{Summa:2018,Obergaulinger:2020}. There is no significant change in the equatorial plane. Both effects lead to a smaller integrated neutrino heating rate for rotating models. Anisotropic emission of neutrinos in rotating CCSN has already been reported in early multi-dimensional simulations \citep{Fryer:2000,Kotake:2003,Buras:2003,Thompson:2004,Marek:2009}. The described behaviour may be artificially enhanced in two-dimensional simulations, although, modern simulations in three dimensions also predict smaller luminosities in rapidly rotating cases \citep{Summa:2018}. However, non-axisymmetric spiral modes of the standing accretion shock instability (SASI) can increase the mass in the gain region and compensate the lack of specific heating, leading to a fundamentally different explosion behavior in 3D \citep{Nakamura:2014,Summa:2018,Kuroda:2020}.

\begin{figure*}[ht]
 \centering
 \includegraphics[width=0.9\linewidth]{rotation_angular_F150.pdf}
 \caption{Neutrino luminosities for non-rotating (black), moderately (blue), and rapidly rotating (red) for the s15.0 progenitor with a heating factor $f_\mathrm{heat} = 1.5$. The columns correspond to observer directions along the poles (left and right panels) and the equatorial plane (center panels). The different rows show the luminosity for  electron neutrinos (top), antineutrinos (middle), and heavy-lepton neutrinos (bottom).}
 \label{fig:rotation_angular_lumAll}
\end{figure*}

\section{Long-time evolution}
\label{sec:longtime}
The 32 exploding models are evolved for several seconds (see Table~\ref{tab:model_overview}) and we define the ``long-time'' phase starting one second after the explosion, i.e., $t > t_\mathrm{exp} + 1$~s. The long-time evolution is characterized by a continuous increase of the explosion energy due to long-lasting accretion onto the proto-neutron star. This accretion depends on the formation and evolution of downflows and it is correlated with rotation and to the morphology of the early explosion phase. In the following, we discuss the generation of explosion energy (Sect.~\ref{sec:Eexp}), how this is affected by rotation and by the explosion morphology (Sect.~\ref{sec:long-time_dshock}), and the evolution of key quantities during the long-time evolution phase (Sect.~\ref{sec:long-time_observables}).

\subsection{Explosion-energy generation}
\label{sec:Eexp}

The explosion energy ($E_\mathrm{exp}$) is calculated by adding the energy of all unbound matter\footnote{In the calculation of the explosion energy we do not consider the contribution of the outer layers  \cite[see e.g.,][]{Bruenn:2013,Bruenn:2016}. We estimate that this ``overburden'' contribution will be around $E_\mathrm{ov} \ll 0.1$~B in our models, which is not important for our overall qualitative description.}  and it is shown in Fig.~\ref{fig:eexp} for the three progenitors including different heating factors and rotation rates. We note that our definition corresponds to the ``diagnostic'' explosion energy, in the sense that it is not the final, saturated value. Shortly after the explosion, $E_\mathrm{exp}$ rapidly grows and in some cases it looks like it saturates and stays constant. However, for all models the explosion energy continues slowly increasing during seconds due to the long-lasting accretion \citep[see also][]{Bruenn:2013,Bruenn:2016, Nakamura:2015, Nakamura:2016, Mueller:2017, Summa:2016, Summa:2018, Harris:2017, OConnor:2018_1, Vartanyan:2018, Vartanyan:2019_1, Glas:2019_1, Mueller:2019_1, Stockinger:2020, Bollig:2020}. This is in contrast to one-dimensional models, that by definition cannot account for downflows and the explosion energy saturates promptly and stays constant \citep[e.g., see][]{Arcones:2007,Perego:2015,Couch:2020}. Even if long-lasting downflows in 2D simulations may become stable and artificially stay for too long time, there are also 3D simulations showing downflows and accretion at late times \citep{Burrows:2019,Vartanyan:2019_1,Stockinger:2020,Bollig:2020}. With our large set of models, we can investigate the  impact of different key aspects on the explosion energy.

\begin{figure}[ht]
 \centering
 \includegraphics[width=0.5\linewidth]{edia_few_s11_1.pdf}%
 \includegraphics[width=0.5\linewidth]{edia_few_s11_2.pdf}\\
  \includegraphics[width=0.5\linewidth]{edia_few_s15_1.pdf}%
 \includegraphics[width=0.5\linewidth]{edia_few_s15_2.pdf}\\
  \includegraphics[width=0.5\linewidth]{edia_few_s27_1.pdf}%
 \includegraphics[width=0.5\linewidth]{edia_few_s27_2.pdf}
 \caption{Same as in Fig.~\ref{fig:rshk}, but for the explosion energy. The dashed and dotted purple lines represent the models s15\_F120\_R001 and s15\_F150\_R030, respectively.}
 \label{fig:eexp}
\end{figure}

Downflows are a multi-dimensional, long-time, and angular dependent feature that is correlated to the explosion energy growth. We investigate the angular dependency of the downflows and explosion energy evolution defining the  explosion-energy growth rate as $\dot{E}_\mathrm{exp} \equiv \mathrm{d}E_\mathrm{exp}/\mathrm{d}t$ and comparing to the mass accretion at late times. Downflows are  derived from the mass accretion rate ($\dot{M}_\mathrm{acc}$) through a sphere with 500~km radius around the PNS. For model s15\_F120\_R001, Figure~\ref{fig:angular_mdot_dedia_F120_R001} shows the direction of downflows over time (upper panel) and how it is correlated to the explosion-energy generation (bottom panel). Before the explosion ($t \leq 0.3$~s post-bounce), there is only accretion. The explosion energy is first accumulated isotropically at shock revival. During the first seconds, the shock expands prominently towards the southern hemisphere leaving space for a long-lasting downflow from the northern hemisphere. Initially, the explosion energy grows also in directions of downflows, because we measure the accretion only at 500~km radius, but explosion energy is also accumulated below that radius. When the downflow is firmly established after $t \approx 1$~s, explosion energy is generated almost exclusively in the southern hemisphere.  The accreted matter acts as fuel for the supernova energy. The strong downflow from the northern hemisphere gradually changes its direction during $2~\mathrm{s} < t < 4~\mathrm{s}$ to the southern hemisphere and the explosion energy generation rate shifts its direction accordingly. The gravitational pull of the PNS decreases the velocity of some ejecta and leads to initially unbound matter to be bound again, which can be seen in slightly negative values of $\dot{E}_\mathrm{exp}$ in the southern hemisphere. Until the end of the simulation at $t = 6.74$~s, the downflow changes direction again and there are no signs of vanishing mass accretion yet.

\begin{figure}[ht]
 \centering
 \includegraphics[width=1\linewidth]{angular_mdot_dedia_F120_R001.pdf}
 \caption{Directionality of downflows and explosion energy generation in model s15\_F120\_R001. The upper panel shows mass downflows and outflows (relative strength, at $r = 500$~km) in a $t$-$\theta$-plane, i.e., direction-dependent, with equator being perpendicular to the cylindrical axis. The lower panel shows the explosion energy generation rate in the same plane. The direction in which explosion energy is generated depends on the
 direction of downflows. Time $t=0$~s corresponds to bounce.}
 \label{fig:angular_mdot_dedia_F120_R001}
\end{figure}

In some cases, we obtain explosions with less stable downflows where even a neutrino-driven wind (NDW) can form as shown in Fig.~\ref{fig:angular_mdot_dedia_F150_R030} for the model s15\_F150\_R030. The NDW phase ($3.5~\mathrm{s} < t < 6.5~\mathrm{s}$) is characterized by a vanishing mass accretion rate, $\dot{M}_\mathrm{acc} = 0$ and  matter is ejected in all directions (subsequently, we define NDW phases by having $\dot{M}_\mathrm{acc} = 0$ for at least $10$~ms duration). During this wind phase, considerably less explosion energy  is added. However, $\dot{E}_\mathrm{exp}$ does not vanish completely, since there is a continuous outflow of matter ejected by neutrinos when depositing energy in the layers around the PNS. At $t \approx 6.5$~s, the wind is terminated by a downflow, and consequently a new phase of explosion energy-generation sets in. This transition occurs almost instantaneously. For all models that develop a NDW (see Sect.~\ref{sec:long-time_observables}), we see typical integrated values of $\dot{E}_\mathrm{exp} \approx 0.04$~B~s$^{-1}$ during phases of no accretion.  Our results suggest that this ongoing interplay of accretion and ejection can continue for longer than previously thought. However, 3D simulations would be necessary to conclude the impact and duration of downflows and NDWs.

\begin{figure}[ht]
 \centering
 \includegraphics[width=1.\linewidth]{angular_mdot_dedia_F150_R030.pdf}
 \caption{Same as in Fig.~\ref{fig:angular_mdot_dedia_F120_R001}, but for model s15\_F150\_R030 which has a phase of no accretion corresponding to a neutrino-driven wind (NDW).}
 \label{fig:angular_mdot_dedia_F150_R030}
\end{figure}

\subsection{Impact of rotation and explosion morphology}
\label{sec:long-time_dshock}
The generation of explosion energy through mass accretion at late times is a robust mechanism in our models. At one second after the explosion (i.e., at the start of the long-time phase), the initial shock wave has reached a radius of $(1-2) \cdot 10^{4}$~km, has cooled down to $T\sim 0.5-1.5$~GK, and has reached densities of $\rho \sim 10^5~\mathrm{g~cm^{-3}}$. At this point, the explosion energy stored within the shock wave and the ejecta generated earlier is approximately constant and any additional $E_\mathrm{exp}$ generation originates from the vicinity of the PNS. Therefore, there is a correlation between $\dot{M}_\mathrm{acc}$ and  $\dot{E}_\mathrm{exp}$ for all exploding models as shown in Fig.~\ref{fig:scatter_dedia_vs_mdot}. In this figure both quantities are shown for the long-time evolution, i.e. $t > t_\mathrm{exp} + 1$~s. The gray dots correspond to individual times (each dot is a mean value over a 5~ms time interval) of all simulations and the colored symbols are obtained by averaging during the whole long-time phase. This introduces some bias for models with  short simulation times because the mass accretion is higher during the the first seconds. In any case,  there is a  correlation  of late-time accretion and explosion-energy growth rate.
In NDW phases of no accretion, the explosion energy growth rate adopts values of $\dot{E}_\mathrm{exp,NDW} = (0.035 \pm 0.007)~\mathrm{B~s^{-1}}$ (left panel of Fig.~\ref{fig:scatter_dedia_vs_mdot}). We note that this specific value might be influenced by our use of a gray neutrino leakage scheme and Newtonian gravity.

\begin{figure}[ht]
 \centering
 \includegraphics[width=1.0\linewidth]{scatter_dedia_vs_mdot_withNDW.pdf}
 \caption{Right panel: Growth-rate of the explosion energy versus mass accretion at late times (i.e.,$t>t_\mathrm{exp}+1$~s). The gray dots represent individual data points (shown as running averages of 5~ms length) from all simulations during the long-time evolution, while the colored symbols indicate the mean values from different models for that phase.  Left panel: Density distribution of $\dot{E}_\mathrm{exp}$ during NDW phases with $\dot{M}_\mathrm{acc} = 0$.}
 \label{fig:scatter_dedia_vs_mdot}
\end{figure}

The early shock morphology has a clear impact on the formation of stable downflows and thus on the late evolution of accretion. In order to quantify the shock morphology, we employ the shock deformation parameter introduced by \citet{Scheck:2006}. It uses the shock radius $R_\mathrm{s}$ as a function of the polar angle $\theta$ and is given by
\begin{equation}
 \label{eq:dshock}
 d_\mathrm{shock} \, = \, \frac{ \mathrm{max} [ R_\mathrm{s}(\theta) \, \mathrm{cos}(\theta) ] - \mathrm{min} [ R_\mathrm{s}(\theta) \, \mathrm{cos}(\theta) ] }
                    {2 \cdot \mathrm{max} [ R_\mathrm{s}(\theta) \, \mathrm{sin}(\theta) ]} - 1 \, .
\end{equation}
The parameter is equivalent to the ratio of the maximum shock diameters, parallel and perpendicular to the cylindrical axis. It can have positive and negative values, for a prolate and an oblate shock deformation, respectively. In the case of a spherical shock expansion, $d_\mathrm{shock}$ becomes zero. We find a  correlation between the shock deformation parameter at the time of shock revival and the rotation strength   (see Tab.~\ref{tab:model_overview}).  Non-rotating 2D simulations typically explode in a prolate morphology, which is a known feature of this geometry \citep[see e.g.,][]{Mueller:2015_2, Nakamura:2015, Bruenn:2016, Summa:2016, OConnor:2018_1, Vartanyan:2018}. However, the resulting accretion along preferred directions has also been observed in 3D simulations \citep{Burrows:2019,Vartanyan:2019_1}.

Rapidly rotating models tend to have less accretion due to  the increased centrifugal forces, which reduce the infall velocity of matter from the equatorial plane and therefore also the accretion rate. The impact of rotation is present in all phases of a CCSN. With increased rotation rate, the shock morphology becomes more spherical and eventually also slightly oblate, for rapidly rotating models. An exception here are simulations of the s11.2 progenitor, due to the smaller amount of angular momentum  (a factor of $10$ less than s15.0 models with the same rotation rate) and an earlier explosion time, which allows for less angular momentum to be accreted until shock revival.

We can investigate the relation between shock deformation (Eq.~\ref{eq:dshock}) and late mass accretion. Downflows usually originate from directions with moderate early shock expantion. In those directions,  part of the ejected matter  does not reach the escape velocity and eventually falls back onto the PNS. Figure~\ref{fig:scatter_mdot_vs_dshock} shows the relation of the early shock morphology $d_\mathrm{shock}(t_\mathrm{exp})$ and the mass accretion at late times. There is not a clear correlation of these two quantities but some trends. Simulations exploding with $d_\mathrm{shock}(t_\mathrm{exp}) \geq 0.5$ are typically moderately rotating or non-rotating (with the exception of the model s11\_F100\_R020) and consistently have an average mass accretion at late times of $\dot{M}_\mathrm{acc} > 0.01~M_\odot/\mathrm{s}$. Furthermore, we can separate all simulations in two groups depending on whether they have some phase of zero accretion. The models that develop a NDW are in the lower left corner region corresponding to lower late accretion and not extreme shock deformation.

\begin{figure}[ht]
 \centering
 \includegraphics[width=1.0\linewidth]{scatter_mdot_vs_dshock.pdf}
 \caption{Average late-time accretion rate versus the shock deformation parameter at explosion time for the three progenitors (different symbols) and various rotation strengths (different colours). Models within the shaded area have a phase of zero accretion or NDW during the longtime evolution.}
 \label{fig:scatter_mdot_vs_dshock}
\end{figure}

\subsection{Evolution in the long-time phase and the neutrino-driven wind}
\label{sec:long-time_observables}

As shown in previous section, there is a strong link between the shock deformation and the formation of stable downflows that critically affect the explosion energy. Here, we want to study the long-time evolution of the explosion energy and mass ejection. For the s15 progenitor, we compare  these quantities at two different times, $1$~s and $5$~s after shock revival.

\begin{figure}[ht]
 \centering
 \includegraphics[width=0.9\linewidth]{scatter_edia_change_vs_mdot.pdf}\\
 \includegraphics[width=0.9\linewidth]{scatter_mej_change_vs_mdot.pdf}
 \caption{Influence of the long-time evolution on the explosion energy (upper panel) and mass ejection (bottom panel) related to the mass accretion. Colors indicate the rotation strength. Open and close circles show the explosion energy at $1$~s and $5$~s after shock revival, respectively. The mass accretion is averaged   for the time interval $1 < t - t_\mathrm{exp} < 5$~s denoted as LT for long-time.}
 \label{fig:change_1_5s}
\end{figure}

The evolution of the explosion energy and its dependence on the mass accretion is shown in Fig.~\ref{fig:change_1_5s}, upper panel. At $1$~s after shock revival, the explosion energy of all models is  distributed at $(0.4 \pm 0.2)$~B, independent of the rotation strength and mass accretion. During the next $4$~s, rapidly rotating models increase their explosion energy by  $\sim 0.4$~B and end up with  values below $0.8$~B at $t_\mathrm{exp} + 5$~s. In contrast,  moderately rotating  and non-rotating models gain considerably more energy during this same time reaching around 0.9-1~B. The value for non-rotating models is also in good agreement with recent 3D studies \citep{Bollig:2020}. We note that for simulations that run for a longer time, the final explosion energy keeps increasing after $t_\mathrm{exp} + 5$~s and the trend with rotation is consistent also at later times. Therefore, we conclude that, for models with initial morphology favourable for large-scale downflows,  the saturation point of the explosion energy lies beyond $10$~s after the initial explosion.
\begin{figure}[ht]
 \centering
 \includegraphics[width=1.0\linewidth]{rotation_windTimes.pdf}
 \caption{Timelines of s15.0 progenitor simulations, where red areas mark the presence of an isotropic neutrino-driven wind phase. We omit the ``s15\_'' prefix in the model names.}
 \label{fig:wind_duration}
\end{figure}

The long-time evolution is also critical for nucleosynthesis (Sect.~\ref{sec:nucleosynthesis}), therefore here we investigate  the amount of mass ejected (Fig.~\ref{fig:change_1_5s}, bottom panel). At $1$~s after shock revival, models with lower accretion (corresponding to approximately spherical explosions, Fig.~\ref{fig:scatter_mdot_vs_dshock}) have ejected more mass than those with strong downflows and a prolate morphology. This is due to the more energetic shock expansion into the equatorial direction in spherical models compared to prolate explosions. This trend continues also in the next seconds. We calculate the average energy of ejected matter: For low accretion models (spherical explosions), we obtain values of $E_\mathrm{exp}/M_\mathrm{ej} \approx 0.5~\mathrm{B}/M_\odot$, while models with high accretion and prolate morphologies reach $E_\mathrm{exp}/M_\mathrm{ej} \approx 1.5~\mathrm{B}/M_\odot$ at $5$~s.

\begin{figure*}[!ht]
 \centering
 \includegraphics[width=0.8\linewidth]{wind_2D_twoPanel.pdf}
 \caption{Density (top) and radial velocity (bottom) of model s15\_F150\_R030 (left) at 6.3~s shortly before termination of the wind (compare to Fig.~\ref{fig:angular_mdot_dedia_F150_R030}~\&~\ref{fig:wind_duration}), and model s15\_F150\_R020 (right) at 8~s after bounce. Note the different axes scales for both models.}
 \label{fig:wind_2D_twoPanel}
\end{figure*}
Even if a considerable amount of mass is ejected in all exploding models, a neutrino-driven wind develops only in 12 of the 21 exploding models for the s15 progenitor. These isotropic NDW phases ($\dot{M}_\mathrm{acc} = 0$ for at least $10$~ms duration) are shown in red in Fig.~\ref{fig:wind_duration}. The first wind phases start to appear after $t \approx 3$~s. There are two kinds of NDW phases: short- and long-duration. The termination of a NDW is usually due to large accumulations of matter with negative radial velocity above the wind. We show examples of the two possible NDW phases in Fig.~\ref{fig:wind_2D_twoPanel}. Here, the symmetry axis is displayed horizontally and we show the density in the upper half of the domain, and the radial velocity in its lower half. The NDW regions are visible around the high-density region with the PNS in the center, and are characterized by high, supersonic velocities, a steep density gradient and the wind termination shock, where the density (velocity) increases (decreases) abruptly. In the example of model s15\_F150\_R030, one can see the matter accumulations with negative radial velocity and comparably high density, just above the wind termination shock. Long-duration winds can extend up to several $10^4$~km radius (see model s15\_F150\_R020 in Fig.~\ref{fig:wind_2D_twoPanel}). The duration of the wind depends on the heating factor and rotation rate. When increasing the heating factor and/or the rotation, there is lower mass accretion  and this allows a wind to form and last for several seconds. Other 2D and 3D simulations indicate also that there is not always a neutrino-driven wind after a successful explosion \citep[see, e.g., ][]{Bruenn:2016,OConnor:2018_1,Vartanyan:2019_1,Bollig:2020}.

\section{Nucleosynthesis}
\label{sec:nucleosynthesis}

Since our neutrino treatment is very simple, the following nucleosynthesis results are only approximate, even if we correct the neutrino properties affecting the evolution of the electron fraction. We focus only on the 21 exploding simulations of the 15~$M_\odot$ progenitor (see Table~\ref{tab:model_overview}). The amount of ejected mass  (Fig.~\ref{fig:change_1_5s}) corresponds to several thousands of ejected tracers per model. These ejected tracers can be divided into ``neutrino-processed'' and ``shock-processed''. Tracers are considered neutrino-processed, when their electron fraction changes by $|\Delta Y_\mathrm{e}| > 0.01$ compared to their starting value. Most of these tracers reach rather high temperature and enter NSE resetting the progenitor composition. Shock-processed particles retain their original electron fraction, their temperature is high enough to change slightly the progenitor composition but most of them stay too cold to enter NSE.

\begin{figure}[!ht]
 \centering
 \includegraphics[width=\linewidth]{F120_R006_allYeDist2.pdf}
 \caption{Modified  distributions for the initial $Y_e$ of model s15\_F120\_R006 when the trajectory cools down to $10$~GK, or at their peak temperature for cold tracers. The black filled bars correspond to the cold (shock-processed) tracers for all distributions. The blue bars show the distribution with the original neutrino energies. The red bars show the distributions with fixed neutrino energy difference and fixed number luminosity ratios, which are indicated above.}
 \label{fig:F120_R006_allYeDist}
\end{figure}

For the cold, shock-processed tracers that do not reach NSE, the evolution of $Y_e$ is not affected by neutrinos.  Therefore, for these tracers we use the electron fraction obtained in the simulations. In contrast, the hot tracers are sensitive to neutrino quantities that determine $Y_e$, namely neutrino and antineutrino energies and number luminosities \citep{Qian:1996}. Therefore, for those tracers we try different corrections to cover all possible conditions. We investigate two different corrections to the $Y_e$ based on one model and later extend our study to the remaining 20 explosions of the 15~$M_\odot$ progenitor. We select model s15\_F120\_R006 as reference here because it is the closest to the ``s15'' model in \citet{Wanajo:2018} when comparing $t_\mathrm{exp}$, $E_\mathrm{exp}$, $M_\mathrm{PNS}$ and $M_\mathrm{ejected}$. However, our neutrino energies and luminosities based on the neutrino leakage scheme lead to unrealistically neutron-rich conditions.

\begin{figure*}[!ht]
 \centering
  \includegraphics[width=0.9\linewidth]{Y_Z_ye.pdf}
 \caption{Abundances of all ejected neutrino-processed tracers of model s15\_F120\_R006 based on broad (left) and narrow with $RL_{\mathrm{n}} = L_{\mathrm{n},\bar{\nu}_e} / L_{\mathrm{n},\nu_e}=1$ (right) $Y_e$ distributions  shown in  Fig.~\ref{fig:F120_R006_allYeDist}. The colours of the lines indicate the electron fraction of the individual tracers at 5.8~GK.}
 \label{fig:ab_F120R006}
\end{figure*}

In the determination of the electron fraction, there are two critical quantities: the energy difference between antineutrinos and neutrinos ($\Delta \epsilon_\nu =  \langle\epsilon_{\bar{\nu}_e}\rangle - \langle\epsilon_{\nu_e}\rangle$) and the ratio of neutrino number luminosities ($RL_{\mathrm{n}} = L_{\mathrm{n},\bar{\nu}_e} / L_{\mathrm{n},\nu_e}$). Following previous studies based on simulations with accurate neutrino transport  \cite[see e.g.,][]{Liebendoerfer:2005_2,Bruenn:2016,Takiwaki:2016,Kotake:2018,Cabezon:2018,Just:2018,OConnor:2018_1,Summa:2018,Vartanyan:2018,Vartanyan:2019_2,Mueller:2019_2,Pan:2019,Powell:2020,Kuroda:2020}, one can find typical energy differences of $\Delta \epsilon_\nu =  \langle\epsilon_{\bar{\nu}_e}\rangle - \langle\epsilon_{\nu_e}\rangle \approx 2 - 2.5$~MeV. In our models we have on average $\Delta \epsilon_\nu = (2.3 \pm 0.6)$~MeV after the explosion time.

In our first approach, we use the original neutrino energies from our simulation and  correct the luminosities similarly to \citet{Sieverding:2020}, but adopting a constant luminosity ratio of $L_\mathrm{\mathrm{n},\bar{\nu}_e} /L_{\mathrm{n},\nu_e}= 1.25$ instead of $1$. This results in a $Y_e$ distribution of ejected matter  that is comparable to the studies of \citet{Wanajo:2018}  and \citet{Sieverding:2020}. The blue histogram in Fig.~\ref{fig:F120_R006_allYeDist} shows the corrected electron fraction distribution, with the cold tracers as black filled bars. The distributions for the hot tracers correspond to the initial electron fraction in the network calculations and corresponds to a temperature of $\sim 10$~GK. For the hot tracers, the  abundances obtained based on this distribution are shown in Fig.~\ref{fig:ab_F120R006} (left panel) where each line corresponds to an individual tracer and the colors indicate the $Y_e$ value at 5.8~GK, i.e. around the temperature when the approximation of NSE breaks down\footnote{Notice that between the initial $Y_e$ at  $\sim 10$~GK shown in Fig.~\ref{fig:F120_R006_allYeDist} and the $Y_e$ at 5.8~GK (Fig.~\ref{fig:ab_F120R006}), the network assumes NSE but the weak reactions are still evolved leading to an evolution of the electron fraction.}.  Iron-group nuclei dominate the final abundances and few tracers at the extremes of the $Y_e$ distribution reach the region of Sr, Y, and Zr.

In our second approach to correct the $Y_e$, we change our neutrino energies to match more closely the literature values. We subtract $4$~MeV from the $\bar{\nu}_e$ leakage energy and adopt a constant energy difference of $\Delta \epsilon_\nu =  \langle\epsilon_{\bar{\nu}_e}\rangle - \langle\epsilon_{\nu_e}\rangle = 2$~MeV. Moreover, we parametrize the ratio of number luminosities, $RL_\mathrm{n} = L_{\mathrm{n},\bar{\nu}_e} / L_{\mathrm{n},\nu_e} = 0.70, 0.85, 1.00, 1.15, 1.35$, where a smaller  ratio corresponds to more proton-rich conditions. This correction leads to very narrow $Y_e$ distributions as shown in Fig.~\ref{fig:F120_R006_allYeDist}. These distributions disagree with those from state-of-the-art simulations, however they are very useful to understand how the abundances depend on a given electron fraction in the supernova ejecta and to cover all possible conditions. Figure~\ref{fig:ab_F120R006} (right panel) shows the abundances for the hot tracers of model s15\_F120\_R006 based on the narrow distribution with $RL_\mathrm{n} = L_{\mathrm{n},\bar{\nu}_e} / L_{\mathrm{n},\nu_e} =1$. In that case as well, the majority of the produced nuclei lie in the iron peak, with a small portion of tracers producing nuclei with $Z>30$.

\begin{figure}[ht]
 \centering
 \includegraphics[width=1.0\linewidth]{Z_F120_R006_all_ye.pdf}
 \caption{ Integrated abundances of all trajectories of  model s15\_F120\_R006 and  different $Y_e$ distributions as given in Fig.~\ref{fig:F120_R006_allYeDist}.}
\label{fig:ab_overview_F120_R006}
\end{figure}

The abundances for cold  and hot tracers with different $Y_e$ treatment are shown in Fig.~\ref{fig:ab_overview_F120_R006} for model s15\_F120\_R006. In general, the abundances for cold and hot tracers show a clear iron peak. The cold component is characterized by the odd-even distribution from carbon to calcium, following the progenitor composition. The hot component reaches heavier elements than the cold one with a slight dependency of the abundances on the exact electron fraction.  The different assumptions for the initial $Y_e$ lead to some variability for the abundances beyond iron covering all expected possibilities for yields from neutrino-driven supernovae. The production of elements in the region of Sr, Y, Zr is more efficient for slightly neutron-rich conditions, corresponding to $RL_{\mathrm{n}} = 1.35$ \citep{Arcones:2014, Arcones.Montes:2011}. We observe a very similar behavior  and dependency of abundances on the $Y_e$ distribution for all  exploding models of the 15~$M_\odot$ progenitor. As a summary, we show in Fig.~\ref{fig:ab_overview}  the final abundances only for the narrow $Y_e$ distribution with $RL_\mathrm{n}=1$ for all models. Due to our necessary but artificial correction of the electron fraction, we do not find any clear correlation between the abundances and the explosion energy, accretion, and rotation. Simulations with detailed neutrino transport are necessary to narrow the uncertainties in the abundances and to link those to other astrophysical conditions. Even if our results are not completely conclusive for abundances beyond iron, they indicate that there is not a strong variability for iron group elements. This may be important to estimate uncertainties in the production of $^{56}$Ni and $^{44}$Ti as shown in Fig.~\ref{fig:all_Ni_Ti} where we present an overview of all models including cold and hot components and variations of $Y_e$.
\begin{figure}[ht]
 \centering
  \includegraphics[width=1.0\linewidth]{Z_all_ye_r100.pdf}
 \caption{ Integrated abundances of all trajectories for all exploding models of s15 and the narrow $Y_e$ distributions with $RL_\mathrm{n}=1$. The black lines represent the abundances from the shock--heated tracers from the different models.}
\label{fig:ab_overview}
\end{figure}

\begin{figure}[ht]
 \centering
 \includegraphics[width=1.0\linewidth]{yields_Ni_Ti_total.pdf}
 \caption{Yields of $^{56}$Ni, $^{44}$Ti obtained by post-processing the 21 exploding models. For each model, the yields are shown for the cold tracers (empty circles) and the neutrino-processed, hot tracers (filled symbols) for three different $Y_e$ narrow distributions corresponding to $RL_\mathrm{n}=$0.7, 0.85, and 1 (triangles, circles, and crosses, respectively). }
\label{fig:all_Ni_Ti}
\end{figure}

\section{Conclusions}
\label{sec:conclusions}

We have presented a broad study of long-time effects based on two-dimensional simulations of CCSN. We define long-time as starting one second after the explosion and following the evolution up to around 10~s post bounce. Our study is based on three  progenitors and variations of neutrino heating  and rotation rates. In total, we present 46 models with 32 of them exploding. In order to be able to run so many simulations for several seconds, we have used for the neutrinos a simple leakage scheme. We are aware that the 2D and neutrino treatment simplifications are critical to make any quantitative conclusion. However, we have found interesting features and trends that we expect to be present in 3D, detailed neutrino transport simulations. Our results indicate that the evolution during the first seconds after shock revival is not negligible, and that accretion and long-lasting downflows impact the growth of the explosion energy. The evolution during the long-time phase can be linked to initial conditions such as the rotation rate and the shock deformation. We find that neutrino-driven winds (NDWs) appear favorably in models with increased neutrino heating and rotation, and can be either long lasting (stable) or short lived (unstable). Improved simulations following the evolution  after the explosion are required to understand CCSN and the connection to observables.

During the explosion phase, both rotation and increased neutrino heating impact the shock acceleration and the initial explosion energy. Rotation weakens the explosion by decreasing the mass accretion and resulting accretion luminosity, which agrees with previous studies \citep{Fryer:2000, Kotake:2003, Buras:2003, Thompson:2004, Marek:2009,Summa:2018,Obergaulinger:2020}. On the other hand, an increased heating factor typically leads to more powerful explosions \citep{OConnor:2010,Couch:2014,Couch:2015_1}.

The long-time evolution is important to estimate  observables such as $E_\mathrm{exp}$, $M_\mathrm{ej}$, and $M_\mathrm{PNS}$. We have found that simultaneous accretion and ejection of matter can persist for several seconds after shock revival and delay the saturation of these observables. In particular, the late-time mass accretion onto the PNS is  correlated to the grow of the explosion energy, and governs the directionality of explosion-energy generation. Besides different progenitors, neutrino heating, and rotation strengths, we find that the shock deformation at the onset of the explosion affects the mass accretion of the following seconds. Rapidly rotating models typically explode with a less prolate deformation, suffer less from persistent downflows, and accumulate less explosion energy in the long-time phase. The amount of ejected mass during the long-time phase  is larger in fast rotating  than in  non-rotating models. We note that the prolate/oblate shock deformations in our models are heavily biased by our 2D geometry \citep{Mueller:2015_2,Nakamura:2015,Bruenn:2016,Summa:2016,OConnor:2018_1,Vartanyan:2018}.

The occurrence of neutrino-driven winds is limited in our models to simulations with the s15 progenitor  and to late times ($t > 3$~s). The total time spent in NDW phases is positively correlated with both the heating factor and rotation rate. However, we also see instances of unstable NDWs and continued accretion after such phases, typically leading to a further increase of the explosion energy. The collapse of a wind is typically induced by renewed accretion from the equatorial plane, which is a fundamentally multi-dimensional effect. Long-time simulations in 3D  are required to further investigate the NDW.

For an analysis of the nucleosynthesis in our models, we used a tracer particle scheme that we applied to all exploding simulations of the 15~$M_\odot$ progenitor. Due to our simplified neutrino treatment, we needed to correct the neutrino energies and luminosities to account for electron fractions that are consistent with modern transport schemes. We have presented the nucleosynthesis for cold trajectories that are not exposed to neutrinos as well as for hot trajectories assuming different electron fraction distributions. All trajectories produce predominately iron group elements with hot, neutrino-processed trajectories reaching elements beyond iron depending on the electron fraction. Our results are a good indication for the potential variability in the production of elements around Sr, Y, Zr and indicate that $^{56}$Ni and $^{44}$Ti are mainly produced by the hot component. In the future, longer simulations times  in three dimensions with detailed neutrino transport will be necessary to connect the nucleosynthesis with the astrophysical conditions (e.g., rotation, explosion energy, progenitor star).

\acknowledgements{
We thank M. Eichler, B. M\"uller, E. O'Connor, and A. Perego for helpful discussions. This research was supported by the ERC Starting Grant EUROPIUM-677912 and Deutsche Forschungsgemeinschaft through SFB 1245. This work has benefited from the COST Action “ChETEC” (CA16117) supported by COST (European Cooperation in Science and Technology) and from JINA Center for the Evolution of the Elements (National Science Foundation under Grant No. PHY-1430152). MO and MR acknowledge the support by the Spanish Ministry of Science, Education and Universities (PGC2018-095984-B-I00) and the Valencian Community (PROMETEU/2019/071). MO further acknowledges support from the Ministerio de Ciencia e Innovación via the Ramón y Cajal program (RYC2018-024938-I). S.M.C. is supported by the U.S. Department of Energy, Office of Science, Office  of Nuclear Physics, Early Career Research Program under Award Number DE-SC0015904. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research and Office of Nuclear Physics, Scientific Discovery through Advanced Computing (SciDAC) program under Award Number DE- SC0017955. This research was supported by the Exascale Computing Project (17-SC-20-SC), a collaborative effort of the U.S. Department of Energy Office of Science and the National Nuclear Security Administration. Calculations for this study were conducted on the Lichtenberg high performance computer of the TU Darmstadt (Project IDs 996 \& 1043).}


\section{Introduction}
\input{intro}
\section{Related Works}
\input{related}

\section{Preliminaries and Notations}
In this section, we introduce some important concepts, notations and background knowledge.
\paragraph{Setting of Two-player Turn-based Stochastic Games} In two-player turn-based games (2-TBSG), only one player takes his action at each step. Denote the two players as $P_1, P_2$, which are the max-player and the min-player respectively. We partition the whole action space as $\mathcal A=\mathcal A_1\cup\mathcal A_2$, where $\mathcal A_i$ is the state space of player $P_i$. Since the stochastic game is episodic under our setting, we denote $2H$ as the number of steps in one episode. At each step $h\in[2H]$, when $h$ is an odd number, it's the max-player $P_1$'s turn to observe the current state $s_h$ and take an action $a_h\in\mathcal A_1$, and then we receive the reward $r_h(s_h,a_h)$. Similarly, when $h$ is an even number, the min-player $P_2$ observes the current state $s_h$ and takes the action $a_h\in\mathcal A_2$, and then they receive the reward $r_h(s_h,a_h)$. After taking the action, the system makes transition to a new state $s_{h+1}\sim \mathbb{P}_h(\cdot|s_h, a_h)$. For each action $a\in\mathcal A$, denote $I(a)\in\{1,2\}$ as it indicates of the current player to play, so that $a\in \mathcal A_{I(a)}$.

\paragraph{General Notations} We denote the tabular episodic Markov Game as $\mathrm{MG}(2H,\mathcal S,\mathcal A_1,\mathcal A_2,\mathbb{P},r)$, where $2H$ is the number of steps in each episode, $\mathcal S$ is the set of states, and $(\mathcal A_1,\mathcal A_2)$ are the action sets of the max-player and the min-player respectively. Since we analyze the two-player turned base games in this paper, the odd number steps are the max-player's turn and the even number steps are the min-player's turn. $\mathbb{P}=\{\mathbb{P}_h\}_{h\in[2H]}$ is the collection of transition matrices and $\mathbb{P}_{h}(\cdot|s,a)$ outputs the probability distribution over states when $a$ is the action taken at step $h$ after state $s$. $r_h: \mathcal S\times\mathcal A\rightarrow [0,1]$ is a deterministic reward function at step $h$.

\paragraph{Markov policy, Q-function and Value function} We denote $(\pi,\mu)$ as the Markov policy of the max-player and the min-player respectively. Max-player's policy $\pi$ is a collection of $H$ mappings $\{\pi_h:\mathcal S\rightarrow \mathcal A_1\}_{h=1,3,\cdots,2H-1}$, which maps the current state to the action. Similarly, min-player's policy $\mu$ is a collection of $H$ mappings $\{\mu_h:\mathcal S\rightarrow \mathcal A_2\}_{h=2,4,\cdots,2H}$. We denote $\pi_h(s)$ and $\mu_h(s)$ to represent the following action to take under Markov policy $\pi,\mu$. Next, we define the well-known value functions and Q-functions. we denote $V_h^{\pi,\mu}:\mathcal S\rightarrow \mathbb{R}$ as the value function at step $h$ under Markov policy $\pi,\mu$ which calculates the expected cumulative rewards:
\[V_h^{\pi,\mu}(s):=\mathbb{E}_{\pi,\mu}\left[\sum_{h'=h}^{2H}r_{h'}(s_{h'},a_{h'})~:~s_h=s\right].\]
We also define $Q_h^{\pi,\mu}:\mathcal S\times\mathcal A$ as the Q-function at step $h$ so that $Q_{h}^{\pi,\mu}(s,a)$ calculates the cumulative rewards under policy $(\pi,\mu)$, starting from $(s,a)$ at step $h$:
\[Q_h^{\pi,\mu}(s,a):=\mathbb{E}_{\pi,\mu}\left[\sum_{h'=h}^{2H}r_{h'}(s_{h'},a_{h'})~:~s_h=s, a_h=a\right].\]
It's worth mentioned that the Markov policy of these two players are both deterministic, which means given a current state, their policies point to a specific action, instead of a probability distribution over all actions. That is because under the two-player turn based game setting, the Nash Equilibria is simply a pure strategy.

For simplicity, we introduce the commonly-used notation of operator $\mathbb{P}_h$ which is $[\mathbb{P}_{h}(V)](s,a):=\mathbb{E}_{s'\sim\mathbb{P}_h(\cdot|s,a)}V(s')$ for any value function $V$. By this definition, we have the Bellman equation as follows:
\[Q_h^{\pi,\mu}(s,a)=(r_h+\mathbb{P}_h V_{h+1}^{\pi,\mu})(s,a)\]
holds for all $(s,a,h)\in\mathcal S\times\mathcal A\times [2H]$. We define $V_{2H+1}^{\pi,\mu}(s)=0$ for all $s\in\mathcal S_{2H+1}$.

\paragraph{Best response and Nash equilibria}
Given the Markov policy $\pi$ of the max-player $P_1$, their exists a best response for the min-player, which we denote as $\mu^{\dagger}(\pi)$. It satisfies:
\[V_{h}^{\pi,\mu^{\dagger}(\pi)}(s)=\inf_{\mu}V_{h}^{\pi,\mu}(s)\]
holds for any $(s,h)\in\mathcal S\times [2H]$. For brevity, we denote $V_{h}^{\pi,\dagger}:=V_{h}^{\pi,\mu^{\dagger}(\pi)}$. Similarly, we can also denote $\pi^{\dagger}(\mu), V_{h}^{\dagger,\mu}$ as the best response for the max-player and the corresponding value function. Furthermore, it is well known that there exists Markov policies $\pi^*,\mu^*$ which are the best responses of each other, and they satisfy:
\[V_h^{\pi^*,\dagger}(s)=\sup_{\pi}V_{h}^{\pi,\dagger}(s),~~~~~V_{h}^{\dagger,\mu^*}(s)=\inf_{\mu}V_{h}^{\dagger,\mu}(s)\]
holds for $\forall (s,h)\in\mathcal S\times[2H]$.

\paragraph{Discounted 2-TBSG with Infinite Horizon} In such a game, the max-player $P_1$ and the min-player $P_2$ take turns to take action. However, different from standard 2-TBSG, there is only one episode in this game and there are infinite steps in this episode.

\paragraph{2-TBSG with Linear Function Expression}
A Markov Game $\mathrm{MG}(2H,\mathcal S,\mathcal A,\mathbb{P},r)$ is defined as linear when the probability transition kernels and the reward functions are linear with respect to a given feature map $\phi:\mathcal S\times\mathcal A\rightarrow \mathbb{R}^d$. Specifically, for each $h\in[2H]$, there exists an unknown vector $\mu_h\in\mathbb{R}^d$ and unknown measures $\theta_h=(\theta_h^{(1)},\theta_h^{(2)}, \ldots,\theta_h^{(d)})$ whose degree of freedom is $|S|\times d$, such that for $\forall (s,a)\in\mathcal S\times\mathcal A$:
\[\mathbb{P}_h(s'|s,a)=\langle\phi(s,a),\theta_h(s')\rangle~\text{and}~r(s,a)=\langle\phi(s,a), \mu_h\rangle.\]
For the complete version of the definition, we put it into Assumption \ref{assumption-linear} below.

\paragraph{Mathematical Notations} Let $f(n), g(n)$ be two positive series, we write $f(n)=\mathcal O(g(n))$ or $f(n)\lesssim g(n)$ if there exists a positive constant $C$ such that $f(n)\leqslant C\cdot g(n)$ for all $n$ larger than some $n_0\in\mathbb{N}$. Similarly, we write $f(n)=\Omega(g(n))$ or $f(n)\gtrsim g(n)$ if there exists a positive constant $C$ such that $f(n)\geqslant C\cdot g(n)$ for all $n$ larger than some $n_0\in\mathbb{N}$. If these two conditions hold simultaneously, we write $f(n) \asymp g(n)$ or $f(n)=\Theta(g(n))$. If logarithmic terms should be ignored, then we use the notations $\widetilde{\mathcal O}, \widetilde{\Omega}, \widetilde{\Theta}$.

\begin{algorithm}[!t]
\caption{Optimistic Nash Q-learning on Two-player Turn-based Stochastic Games}
\hspace*{0.02in}{\bf Initialize:}
Let $\overline{Q}_h(s,a)\leftarrow 2H,~ \underline{Q}_h(s,a)\leftarrow 0$, and $N_h(s,a)\leftarrow 0$ for all $(s,a)\in\mathcal S\times \mathcal A$.\\
\hspace*{0.02in}{\bf Define:}
$\alpha_t=\frac{2H+1}{2H+t},~\iota\leftarrow \log(SAT^2)$.\\
\begin{algorithmic}[1]
\FOR{episode $k\in[K]$}
\STATE{observe the initial state $s_1$}
\FOR{step $h\in[2H]$}
\STATE Take action $a_h\leftarrow\arg\max_{a'\in\mathcal A}\overline{Q}_{h}(s_h,a')$ if $h$ is an odd number, (i.e. $I(a_h)=1$),else take action $a_h\leftarrow\arg\min_{a'\in\mathcal A}\underline{Q}_{h}(s_h,a')$. After that, observe the next state $s_{h+1}$.
\STATE $t=N_h(s_h,a_h)\leftarrow N_h(s_h,a_h)+1$.
\STATE $\beta_t\leftarrow c\sqrt{(2H)^3\iota/t}$.
\STATE $\overline Q_h(s_h,a_h)\leftarrow (1-\alpha_t)\cdot\overline Q_h(s_h,a_h)+\alpha_t\cdot [r_h(s_h,a_h)+\overline{V}_{h+1}(s_{h+1})+\beta_t]$.
\STATE $\underline Q_h(s_h,a_h)\leftarrow (1-\alpha_t)\cdot\underline Q_h(s_h,a_h)+\alpha_t\cdot[r_h(s_h,a_h)+\underline{V}_{h+1}(s_{h+1})-\beta_t]$.
\STATE $\overline{V}_h(s_h)\leftarrow \overline{Q}_h(s_h,a_h),~\underline{V}_h(s_h)\leftarrow \underline Q_h(s_h,a_h)$.
\ENDFOR
\ENDFOR
\end{algorithmic}
\label{alg-1}
\end{algorithm}

\section{Regret Bound of Tabular 2-TBSG}
In 2-TBSG, we use a Nash Q-learning method: Algorithm \ref{alg-1} to obtain a policy tuple sequence $(\pi^k,\mu^k)$ to approximate a Nash equilibrium. In this algorithm, $\overline{Q}_h(s,a), \underline{Q}_h(s,a)$ are the upper and lower estimation of $Q_h^*(s,a)$ respectively. The max-player chooses action based on the upper estimation $\overline{Q}$ while the min-player chooses action based on the lower estimation $\underline{Q}$. By using the Upper Confidential Bound (UCB) technique, we can prove that
\[\underline{Q}_h(s,a)\leqslant Q_h^*(s,a)\leqslant \overline{Q}_h(s,a)\]
holds with high probability. The total regret of this algorithm is defined as:
\[\mathrm{Regret}(K)=\sum_{k=1}^{K}\left|\left(V_1^*-V_1^{\pi^k,\mu^k}\right)(s_1^k)\right|.\]
Here, $V_1^*$ is the abbreviation of $V_1^{\pi^*,\mu^*}$. In this section, we will focus on upper bounding the expected regret $\mathbb{E}[\mathrm{Regret}(K)]$. Notice that, unlike our related papers, there is an absolute value in our definition above, which is because under the two-player game setting, $\left(V_1^*-V_1^{\pi^k,\mu^k}\right)(s_1^k)$ can be either positive or negative.
\begin{Theorem}[Main Theorem 1: Logarithmic Regret Bound of Q-learning for Episodic 2-TBSG]
\label{thm-main-1}
After using Algorithm \ref{alg-1}, the expected total regret for episodic two-player turn-based stochastic game (2-TBSG) can be upper bounded by:
\[\mathbb{E}[\mathrm{Regret}(K)]\leqslant \mathcal O\left(\frac{H^6SA\log(SAT)}{\mathrm{gap}^+_{\min}}\right).\]
\end{Theorem}
Here, the definition of $\mathrm{gap}^+_{\min}$ will be left in the following section. Notice that the theorem is not a straight-forward extension of existing works since Algorithm \ref{alg-1} is originally proposed by us and it is the first gap dependent logarithmic regret bound for 2-TBSG. In \cite{xie2020learning}, the authors have proposed a similar algorithm, but their setting is 2-player zero-sum simultaneous-move Markov Games and prove a $\mathcal O(\sqrt{d^3H^3T\iota^2})$ regret bound guarantee. The policies of the two players are mixed strategies, which makes it impossible to define gaps. Next, we provide the proof sketch of Theorem \ref{thm-main-1} and leave all the proof details to appendix.

\subsection{First Step: Split the Total Regret into the Expected Sum of Gaps}
In the first step, we split the total regret defined above into several single-step sub-optimality gaps.
\begin{Lemma}\label{thm-split}
\begin{equation}
\left(V_1^*-  V_1^{\pi^k,\mu^k}\right)(s_1^k)=\mathbb{E}\left[\sum_{h=1}^{2H}\mathrm{gap}_h(s_h^k,a_h^k)~|~\pi^k,\mu^k\right].
\end{equation}
Here, $\mathrm{gap}_h(s_h^k,a_h^k):=V_h^*(s_h^k)-Q_h^*(s_h^k,a_h^k)$.
\end{Lemma}
It is obvious to find that: when $a_h^k\in\mathcal A_1$, which means $a_h^k$ is the max-player's action, then $\mathrm{gap}_h(s_h^k,a_h^k)\geqslant 0$. In contrast, when $a_h^k\in\mathcal A_2$, which means $a_h^k$ is the min-player's action, then $\mathrm{gap}_h(s_h^k,a_h^k)\leqslant 0$. This is because the optimal strategy group $(\pi^*,\mu^*)$ stands for both the optimal point for max-player and min-player, though at different directions. Therefore, we introduce a new notation:
\begin{equation*}
\begin{aligned}
&\mathrm{gap}^{+}_h(s_h,a_h):=\left|\mathrm{gap}_h(s_h,a_h)\right|\\
&~~~=\begin{cases}V_h^*(s_h)-Q_h^*(s_h,a_h)&\mbox{if $a_h\in\mathcal A_1$ (or $h$ is odd)}\\ Q_h^*(s_h,a_h)-V_h^*(s_h)&\mbox{if $a_h\in\mathcal A_2$ (or $h$ is even)}\end{cases}.
\end{aligned}
\end{equation*}
Combining the two equations above, we obtain that:
\begin{align}\label{eqn-split2}
\mathbb{E}[\mathrm{Regret}(K)]&=\sum_{k=1}^{K}\left|\mathbb{E}\left[\sum_{h=1}^{2H}\mathrm{gap}_h(s_h^k,a_h^k)~|~\pi^k,\mu^k\right]\right|\notag\\
&\leqslant \sum_{k=1}^{K}\mathbb{E}\left[\sum_{h=1}^{2H}\mathrm{gap}^{+}_h(s_h^k,a_h^k)~|~\pi^k,\mu^k\right].
\end{align}

\subsection{Second Step: Concentration Inequalities and the Extension}
According to an existing lemma in a related work \cite{bai2020provable}, we have the following conclusion as its special case in the 2-TBSG setting:
\begin{Lemma}[Lemma 11 of \cite{bai2020provable}]
For any $p\in(0,1]$, we let $\iota=\log(SAT/p)$. Then, with probability at least $1-p$, Algorithm 1 has the following property:
\begin{equation}\label{eqn-gap}
\begin{aligned}
&\overline{Q}_h^k(s,a)\geqslant Q_h^*(s,a)\geqslant \underline{Q}_h^k(s,a), \\
&\overline{V}_h^k(s)\geqslant V_h^*(s)\geqslant \underline{V}_h^k(s)
\end{aligned}
\end{equation}
holds for $\forall s\in\mathcal S, h\in[2H], a\in\mathcal A, k\in[K]$. We call the event above $\varepsilon_{\mathrm{conc}}$ and then $P(\varepsilon_{\mathrm{conc}})\geqslant 1-p$.
\end{Lemma}

This lemma indicates that $\overline{V},\overline{Q}$ and $\underline{V},\underline{Q}$ proposed by Algorithm 1 are exactly the upper bound and lower bound of the optimal point (Nash equilibrium) $V^*, Q^*$, with high probability. Similar to the episodic setting, we can make the following conclusion:
\begin{equation}\label{eqn-5}
\mathrm{gap}^{+}_h(s_h^k,a_h^k)\leqslant\overline{Q}_h^k(s_h^k,a_h^k)-\underline{Q}_h^k(s_h^k,a_h^k)
\end{equation}
holds for $\forall h\in[2H], k\in[K]$ as long as $\varepsilon_{\mathrm{conc}}$ holds.
Then, Equation (\ref{eqn-5}) is transformed to:
\begin{equation}\label{eqn-6}
\mathrm{gap}^{+}_h(s_h^k,a_h^k)=\mathrm{clip}\left[\mathrm{gap}^+_h(s_h^k,a_h^k)~|~\mathrm{gap}^+_{\min}\right]\leqslant \mathrm{clip}\left[\overline{Q}_h^k(s_h^k,a_h^k)-\underline{Q}_h^k(s_h^k,a_h^k)~|~\mathrm{gap}^+_{\min}\right].
\end{equation}

\subsection{Third Step: Peeling}
Like \cite{yang2020q}, we separate all the gaps $\overline{Q}_h^k(s_h^k,a_h^k)-\underline{Q}_h^k(s_h^k,a_h^k)$ into different intervals and count them individually. Note that when the gap  $\overline{Q}_h^k(s_h^k,a_h^k)-\underline{Q}_h^k(s_h^k,a_h^k)$ belongs to $[0, \mathrm{gap}_{\min})$, then it will be clipped to 0 by Equation (\ref{eqn-6}). For the other gaps, we divide them into $N$ different intervals $[\mathrm{gap}_{\min}, 2\mathrm{gap}_{\min}), \cdots, [2^{N-1}\mathrm{gap}_{\min}, 2^{N}\mathrm{gap}_{\min})$. Here, $N=\lceil\log_{2}(2H/\mathrm{gap}_{\min})\rceil$. The following lemma tells us the upper bound of the counting number in each interval.
\begin{Lemma}[Bounded Counting Number of Each Interval]\label{lemma-1} Under the concentration event $\varepsilon_{\mathrm{conc}}$, for each $n\in[N]$, we denote:
\[C^{(n)}:=\left|\left\{(k,h):\overline{Q}_h^k(s_h^k,a_h^k)-\underline{Q}_h^k(s_h^k,a_h^k)\in \Lambda_n\right\}\right|.\]
where $\Lambda_n=[2^{n-1}\mathrm{gap}_{\min}, 2^n\mathrm{gap}_{\min})$. Then, we have the following upper bound:
\[C^{(n)}\leqslant \mathcal O\left(\frac{H^6SA\iota}{4^n\mathrm{gap}_{\min}^2}\right)\]
where $\iota=\log(SAT/p)$ is the logarithmic term.
\end{Lemma}
Once we have this lemma proved, we can easily estimate the upper bound of expected cumulative regret.
\begin{align}
&\mathbb{E}[\mathrm{Regret}(K)]\leqslant \mathbb{E}\left[\sum_{k=1}^{K}\sum_{h=1}^{2H}\mathrm{gap}_h^{+}\left(s_h^k,a_h^k\right)\right]\notag\\
\leqslant& \sum_{\varepsilon_{\mathrm{conc}}}\mathbb{P}(\mathrm{traj})\cdot\sum_{k=1}^{K}\sum_{h=1}^{2H}\mathrm{clip}\left[(\overline{Q}_h^k-\underline{Q}_h^k)(s_h^k,a_h^k)~|~\mathrm{gap}_{\min}\right]+\sum_{\mathrm{traj}\notin\varepsilon_{\mathrm{conc}}}\mathbb{P}(\mathrm{traj})\cdot 2TH\notag\\
\leqslant& \sum_{n=1}^{N}2^n\mathrm{gap}_{\min}C^{(n)}+p\cdot 2TH\notag\\
\leqslant& \sum_{n=1}^{N}\mathcal O\left(\frac{H^6SA\iota}{2^n\mathrm{gap}_{\min}}\right)+p\cdot 2TH \leqslant \mathcal O\left(\frac{H^6SA\log(SAT)}{\mathrm{gap}_{\min}}\right).\notag
\end{align}
In the last step, we let $p=\frac{1}{T}$, and then $\iota=\log(SAT^2)=\mathcal O(\log(SAT))$. This leads to our main theorem. However, the proof of Lemma \ref{lemma-1} is difficult. We rely on a general lemma about the upper bound of the weighted sum of $(\overline{Q}_h^k-\underline{Q}_h^k)(s_h^k,a_h^k)$.

\begin{Lemma}[Peeling Argument]\label{lemma-2} Under the event $\varepsilon_{\mathrm{conc}}$, the following holds for $\forall h\in[2H]$ and a weight sequence $\{w_{k,h}\}_{k\in[K]}$ which satisfies: $0\leqslant w_{k,h}\leqslant w,~\sum_{k=1}^{K}w_{k,h}\leqslant C$, it holds that:
\begin{equation}\label{eqn-lemma2}
\sum_{k=1}^{K}w_{k,h}\left(\overline{Q}_h^k-\underline{Q}_h^k\right)(s_h^k,a_h^k)\leqslant 4ewSAH^2+60c\sqrt{SACewH^5\iota}.
\end{equation}
\end{Lemma}

After proving this lemma. We can make \[w_{k,h}^{n}=\mathbb{I}\left[\left(\overline{Q}_h^k-\underline{Q}_h^k\right)(s_h^k,a_h^k)\in\Lambda_n\right]\] and
\[C_h^{(n)}=\sum_{k=1}^{K}\mathbb{I}\left[\left(\overline{Q}_h^k-\underline{Q}_h^k\right)(s_h^k,a_h^k)\in\Lambda_n\right],\]
where $\Lambda_n = [2^{n-1}\mathrm{gap}_{\min},2^{n}\mathrm{gap}_{\min})$. Then, the sequence $\{w_{k,h}^{n}\}_{k\in[K]}$ is a $\left(1, C_h^{(n)}\right)$-sequence. According to Lemma \ref{lemma-2},  we know that:
\begin{equation*}
2^{n-1}\mathrm{gap}_{\min}C_h^{(n)}\leqslant\sum_{k=1}^{K}w_{k,h}^{n}\left(\overline{Q}_h^k-\underline{Q}_h^k\right)(s_h^k,a_h^k)\leqslant 4eSAH^2+60c\sqrt{SAC_h^{(n)}eH^5\iota},
\end{equation*}
which leads to the fact that:
\begin{equation}
C_h^{(n)}\leqslant \mathcal O\left(\frac{H^5SA\iota}{4^n\mathrm{gap}_{\min}^2}\right).
\end{equation}
Finally, after summing them up:
\[C^{(n)}=\sum_{h=1}^{2H}C_h^{(n)}=\mathcal O\left(\frac{H^6SA\iota}{4^n(\mathrm{gap}^+_{\min})^2}\right),\]
which comes to our conclusion of Lemma \ref{lemma-1}.

\paragraph{Connection Between the Vanilla Regret and Duality Regret.} Here we remark on the connection between the vanilla regret and duality regret.
\begin{itemize}
\item In this paper, we analyze the vanilla regret $|V_1^*-V_1^{\pi^k, \mu^k}|$. Since the game setting we are analyzing is 2-player turn-based stochastic game (2-TBSG), and $(\pi^*,\mu^*)$ is a Nash Equilibrium rather than the optimal policy. The vanilla regret actually measures the distance between the two value functions proposed by $(\pi^*, \mu^*)$ and $(\pi^k, \mu^k)$.
\item Another type of regret is the duality regret $V_1^{\dagger, \mu^k}-V_1^{\pi^k, \dagger}$. It can measure how close the policy $(\pi^k,\mu^k)$ is to a Nash Equilibrium.
\end{itemize}

Therefore, it's an important question whether a small vanilla regret implies a small duality regret. (Notice that the inverse would not be true since there might be more than one Nash Equilibria.) To answer this question, we propose the following proposition.

\begin{Proposition}
\label{prop-1}
If for $\forall s\in\mathcal S, h\in [2H]$, the vanilla regret
\[\left|V_h^*(s)-V_h^{\pi,\mu}(s)\right|<\frac{1}{2}\cdot\mathrm{gap}_{\min},\]
we can conclude that $(\pi,\mu)$ is a Nash Equilibrium, which means its duality regret is 0.
\end{Proposition}

\begin{proof}[Proof of Proposition \ref{prop-1}] We are going to prove $(\pi_h, \mu_h)=(\pi^*_h,\mu^*_h)~\forall h\in [2H]$ by using the method of induction. When $h=2H$, notice that for $\forall s\in\mathcal S$:
\[V_{2H}^{\pi,\mu}(s)= r(s, \mu_{2H}(s))=Q_{2H}^*(s, \mu_{2H}(s)).\]
Since $\frac{1}{2}\cdot\mathrm{gap}_{\min}>\left|V_{2H}^*(s)-V_{2H}^{\pi,\mu}(s)\right|=\left|V_{2H}^*(s)-Q_{2H}^*(s,\mu_{2H}(s))\right|$, therefore we have $\mu_{2H}^*(s)=\mu_{2H}(s)$ according to the definition of $\mathrm{gap}_{\min}^+$. It means that $(\pi,\mu)$ and $(\pi^*,\mu^*)$ make the same decisions at the $2H$'s step. On the other hand, assume that $(\pi_t,\mu_t)=(\pi_t^*,\mu_t^*)$ holds for $t = h+1, h+2, \ldots, 2H$, then for $t=h$: if $h$ is an odd number, then it's max-player's turn to take action, and it holds that for $\forall s\in\mathcal S$:
\[V_{h}^{\pi,\mu}(s)=r(s,\pi(s))+\mathbb{E}_{s'|s,\pi(s)} V_{h+1}^{\pi,\mu}(s') \overset{(a)}{=} r(s,\pi(s))+\mathbb{E}_{s'|s,\pi(s)} V_{h+1}^*(s') = Q_h^*(s,\pi_h(s)).\]
Since $\frac{1}{2}\cdot\mathrm{gap}_{\min}>\left|V_h^*(s)-V_h^{\pi,\mu}(s)\right|=\left|V_h^*(s)-Q_h^*(s,\pi_h(s))\right|$, we can conclude that $\pi_h(s)=\pi^*_h(s)$ according to the definition of $\mathrm{gap}_{\min}^+$. Therefore, $(\pi,\mu)$ and $(\pi^*,\mu^*)$ also make the same decisions at the $h$'s step, which finishes our induction. When $h$ is an even number, we can finish our induction in the same way, and that comes to the conclusion of Proposition \ref{prop-1} since $(\pi,\mu)=(\pi^*,\mu^*)$ and then the policy pair $(\pi,\mu)$ is a Nash Equilibrium. It leads to fact that:
\[V_1^{\dagger, \mu}=V_1^{\pi,\mu}=V_1^{\pi,\dagger},\]
which means policy pair $(\pi,\mu)$ has zero duality regret.

\end{proof}

The conclusion shows that when vanilla regret is sufficiently small for all $h\in [2H], s\in\mathcal S$, then the policy pair is a Nash Equilibrium. It bridges the gap between the vanilla regret and the duality regret, and makes it more reasonable and convincing for us to work on the upper bound for the expected sum of vanilla regrets.

\section{Regret Bound of Discounted 2-TBSG}
In this section, we study the discounted 2-TBSG with infinite horizon. In order to obtain a policy tuple sequence $(pi^k,\mu^k)$ to approximate a Nash equilibrium, we propose the following Algorithm \ref{alg-2}, which is similar to Algorithm \ref{alg-1}.
We also use the UCB technique to establish a upper estimation $\hat{Q}(s,a)$ and a lower estimation $\breve{Q}(s,a)$ of the optimal Q-function $Q^*(s,a)$.

\begin{algorithm}[!t]
\caption{Optimistic Nash Q-learning on Discounted 2-TBSG with Infinite Horizon}
\hspace*{0.02in}{\bf Initialize:}
Let $\overline{Q}(s,a), \hat{Q}(s,a)\leftarrow 1/(1-\gamma)$ and  $\underline{Q}(s,a), \breve{Q}(s,a)\leftarrow 0$ for all $(s,a)\in\mathcal S\times \mathcal A$. Also $N(s,a)\leftarrow 0$.\\
\hspace*{0.02in}{\bf Define:}
$\iota(k)=\log(SAT(k+1)(k+2)), \alpha_k=\frac{H+1}{H+k}$ where $H=\frac{\log(2/(1-\gamma)\mathrm{gap}^+_{\min})}{\log(1/\gamma)}$.\\
\begin{algorithmic}[1]
\STATE{Observe the initial state $s_1$.}
\FOR{episode $t\in[T]$}
\STATE Take action $a_t\leftarrow\arg\max_{a'\in\mathcal A}\overline{Q}(s_t,a')$ if $t$ is an odd number, (i.e. $I(a_t)=1$), else take action $a_t\leftarrow\arg\min_{a'\in\mathcal A}\underline{Q}(s_t,a')$. After that, observe the reward $r(s_t,a_t)$ as well as the next state $s_{t+1}$.
\STATE $k=N(s_t,a_t)\leftarrow N(s_t,a_t)+1$.
\STATE $b_k\leftarrow \frac{c_2}{1-\gamma}\sqrt{H\iota(k)/k}$,~~~~~~~Here $c_2$ is a constant that can be set to $4\sqrt{2}$.
\STATE $\overline Q(s_t,a_t)\leftarrow (1-\alpha_k)\cdot\overline Q(s_t,a_t)+\alpha_k\cdot[r(s_t,a_t)+\gamma\hat{V}(s_{t+1})+b_k]$.
\STATE $\underline Q(s_t,a_t)\leftarrow (1-\alpha_k)\cdot\underline Q(s_t,a_t)+\alpha_k\cdot[r(s_t,a_t)+\gamma\breve{V}(s_{t+1})-b_k]$.
\STATE $\hat{Q}(s_t,a_t)\leftarrow \min(\hat{Q}(s_t,a_t), \overline{Q}(s_t,a_t)), ~\breve{Q}(s_t,a_t)\leftarrow \max(\breve{Q}(s_t,a_t), \underline{Q}(s_t,a_t)) $.
\STATE $\hat{V}(s_t)\leftarrow \hat{Q}(s_t,a_t), \breve{V}(s_t)\leftarrow \breve{Q}(s_t,a_t)$.
\ENDFOR
\end{algorithmic}
\label{alg-2}
\end{algorithm}

Since the V-function denotes the expected discounted sum of rewards given the initial state $s$ and the Q function denotes the expected discounted sum of rewards given the initial state $s$ and the initial action $a$, so they can be described as:
\begin{equation*}
\begin{aligned}
V^{\pi}(s) &:= \mathbb{E}\left[\sum_{i=1}^{\infty}\gamma^{i-1}\cdot r(s_i, a_i)~:~s_1=s\right],\\
Q^{\pi}(s,a) &:=\mathbb{E}\left[\sum_{i=1}^{\infty}\gamma^{i-1}\cdot r(s_i, s_i)~:~s_1=s, a_1 = a\right].
\end{aligned}
\end{equation*}
where $\gamma$ is the discounted ratio. In the following statements, $\overline{Q}_t, \overline{V}_t$ stand for the $\overline{Q}, \overline{V}$ functions in the $t$-th iteration, and so does the other subscripts.

\subsection{Sub-optimality Gap and the Splitting of Total Regret}
Similar to the episodic setting, given $(s,a)\in\mathcal S\times\mathcal A$, define $\mathrm{gap}(s,a)$ as:
\[\mathrm{gap}(s,a):=V^*(s)-Q^*(s,a).\]
Notice that, when $a\in\mathcal A_1$, which means the action $a$ is taken by max-player $P_1$, then $V^*(s)\geqslant Q^*(s,a)\Rightarrow \mathrm{gap}(s,a)\geqslant 0$. In contrast, when $a\in\mathcal A_2$, which means the action $a$ is taken by min-player $P_2$, then $V^*(s)\leqslant Q^*(s,a)\Rightarrow\mathrm{gap}(s,a)\leqslant 0$. Here, we can introduce a notation $\mathrm{gap}^+(s,a):=|\mathrm{gap}(s,a)|$ which stands for the absolute value of $\mathrm{gap}(s,a)$. Also, we denote $\mathrm{gap}^+_{\min}$ as the minimum non-zero absolute gap:
\[\mathrm{gap}^+_{\min}:=\min_{s,a}\{\mathrm{gap}^+(s,a)~:~\mathrm{gap}^+(s,a)\neq 0\} > 0.\]
In the main theorem proposed in the following section, we will estimate an upper bound of the expected total regret for the first $T$ steps
\[\mathrm{Regret}(T):=\sum_{t=1}^{T}\left|(V^*-V^{\pi_t, \mu_t})(s_t)\right|.\]

\subsection{Main Theorem}
In this section, we propose our main theorem in the infinite-horizon setting. Unlike dual gap regret \cite{xie2020learning}, the total regret above proposed by \cite{liu2020regret} follows the sample complexity definition in \cite{sm2003on} and directly compares the actual value function and the value function from the first $t$ iterations. Similar to the episodic setting, we can obtain a gap-dependent logarithmic upper bound for the expected total regret.
\begin{Theorem}[Main Theorem 2: Logarithmic Regret Bound of Q-learning for Infinite-horizon Discounted 2-TBSG] After using Algorithm \ref{alg-2}, the expected total regret for infinite-horizon two-player turn-based stochastic game can be upper bounded by:
\begin{equation*}
\mathbb{E}[\mathrm{Regret}(T)]\leqslant \mathcal O\left(\frac{SA}{\mathrm{gap}^+_{\min}(1-\gamma)^5\log(1/\gamma)}\cdot\log\frac{SAT}{\mathrm{gap}^+_{\min}(1-\gamma)}\right).
\end{equation*}
\label{thm-main-2}
\end{Theorem}
Next, we introduce the proof sketch of Theorem \ref{thm-main-2}.

\subsection{First step: Splitting the Regret into Expected Sum of Gaps} \label{sec:split_regret}
The splitting of the regret is just the same as that in the finite horizon setting:
\begin{align}\label{eqn-7}
\mathbb{E}[\mathrm{Regret}(T)]&=\mathbb{E}\left[\sum_{t=1}^{T}\left|\sum_{h=0}^{\infty}\gamma^h\mathrm{gap}(s_{t+h},a_{t+h})\right|\right]\notag\\
&\leqslant\mathbb{E}\left[\sum_{t=1}^{T}\sum_{h'=t}^{+\infty}\gamma^{h'-t}\mathrm{gap}^+(s_{h'},a_{h'})\right].
\end{align}
Notice that when $t$ is an odd number, it's the max-player's turn to take action, so $\mathrm{gap}(s_t,a_t)\geqslant 0$ and then:
\begin{equation*}
\begin{aligned}
&~~\mathrm{gap}^+(s_t,a_t)=\mathrm{gap}(s_t,a_t)=V^*(s_t)-Q^*(s_t,a_t)\\
&=Q^*(s_t,a^*)-Q^*(s_t,a_t)\leqslant \hat{Q}_t(s_t,a^*)-\breve{Q}_t(s_t,a_t)\\
&\leqslant \left(\hat{Q}_t-\breve{Q}_t\right)(s_t,a_t).
\end{aligned}
\end{equation*}
Similarly, when $t$ is an even number, it's the min-player's turn to take action, so $\mathrm{gap}(s_t,a_t)\leqslant 0$, and then:
\begin{equation*}
\begin{aligned}
&~~\mathrm{gap}^+(s_t,a_t)=-\mathrm{gap}(s_t,a_t)=Q^*(s_t,a_t)-V^*(s_t)\\
&=Q^*(s_t,a_t)-Q^*(s_t,a^*)\leqslant \hat{Q}_t(s_t,a_t)-\breve{Q}_t(s_t,a^*)\\
&\leqslant \left(\hat{Q}_t-\breve{Q}_t\right)(s_t,a_t).
\end{aligned}
\end{equation*}
Therefore, we can conclude that $\mathrm{gap}^+(s_t,a_t)\leqslant \left(\hat{Q}_t-\breve{Q}_t\right)(s_t,a_t)$. By the definition of $\mathrm{gap}^+_{\min}$, we have:
\begin{equation*}
\mathrm{gap}^+(s_t,a_t)=\mathrm{clip}[\mathrm{gap}^+(s_t,a_t)|\mathrm{gap}^+_{\min}]\leqslant \mathrm{clip}\left[\left(\hat{Q}_t-\breve{Q}_t\right)(s_t,a_t)|\mathrm{gap}^+_{\min}\right].
\end{equation*}
Combine it with Equation (\ref{eqn-7}), we obtain that:
\begin{equation}
\label{eqn-8}
\mathbb{E}[\mathrm{Regret}(T)]\leqslant \mathbb{E}\left[\sum_{t=1}^{T}\sum_{h'=t}^{+\infty}\gamma^{h'-t}\cdot\mathrm{clip}\left[\left(\hat{Q}_{h'}-\breve{Q}_{h'}\right)(s_{h'},a_{h'})|\mathrm{gap}^+_{\min}\right]\right].
\end{equation}

\subsection{Second step: Concentration Properties}
Extended from \cite{dong2019q}, we can obtain the following lemma which shows that Algorithm \ref{alg-2} satisfies bounded learning error with high probability.
\begin{Lemma}[Concentration Property]
\label{lemma-3}
When applying Algorithm \ref{alg-2}, event $\mathcal E_{\mathrm{conc}}$ occurs with probability at least $1-\frac{1}{T}$. Here, $\mathcal E_{\mathrm{conc}}$ occurs if  $\forall (s,a,t)\in\mathcal S\times\mathcal A\times\mathbb{N}$:
\begin{equation*}
\begin{aligned}
0 &\leqslant \left(\hat{Q}_t-Q^*\right)(s,a)\leqslant \left(\overline{Q}_t-Q^*\right)(s,a)\leqslant \frac{\alpha_{n^t}^0}{1-\gamma}+\sum_{i=1}^{n^t}\gamma\alpha_{n^t}^i\left(\hat{V}_{\tau(s,a,i)}-V^*\right)(s_{\tau(s,a,i)})+\beta_{n^t},\\
0 &\leqslant \left(Q^*-\breve{Q}_t\right)(s,a)\leqslant \left(Q^*-\underline{Q}_t\right)(s,a)\leqslant \frac{\alpha_{n^t}^0}{1-\gamma}+\sum_{i=1}^{n^t}\gamma\alpha_{n^t}^i\left(V^*-\breve{V}_{\tau(s,a,i)}\right)(s_{\tau(s,a,i)})+\beta_{n^t}.
\end{aligned}
\end{equation*}
Here, $\iota(k)=\log(SAT(k+1)(k+2))$ and $\beta_k=\frac{c_3}{1-\gamma}\sqrt{\frac{H\iota(k)}{k}}$.
\end{Lemma}

Then, under $\mathcal E_{\mathrm{conc}}$, we know that:
\begin{equation}\label{eqn-inf-conc}
0\leqslant \left(\hat{Q}_t-\breve{Q}_t\right)(s,a)\leqslant \frac{2\alpha_{n^t}^0}{1-\gamma}+2\beta_{n^t}+\sum_{i=1}^{n^t}\gamma\alpha_{n^t}^i\left(\hat{V}_{\tau(s,a,i)}-\breve{V}_{\tau(s,a,i)}\right)(s_{\tau(s,a,i)}).
\end{equation}

\subsection{Third step: Peeling} \label{sec:peeling}
Similar to Lemma \ref{lemma-2}, we upper bound the weighted sum of upper-lower gaps and then bound the counting number of gaps in different intervals, just like we did in the episodic setting.
\begin{Lemma}[Peeling Argument]
\label{lemma-4}
Under the event $\mathcal E_{\mathrm{conc}}$, the following holds for any weighted sequence $\{\omega_t\}$ which satisfies: $0\leqslant \omega_t\leqslant \omega, \sum_{t=1}^{+\infty}\omega_t\leqslant C$.
\begin{equation*}
\sum_{t=1}^{+\infty}\omega_t\left(\hat{Q}_t-\breve{Q}_t\right)(s_t,a_t)\leqslant \frac{\gamma^H C}{1-\gamma}+\mathcal O\left(\frac{\sqrt{\omega SAHC\iota(C)}+\omega SA}{(1-\gamma)^2}\right).
\end{equation*}
\end{Lemma}

After that, we classify the positive gaps into different intervals. Since all the positive gaps $\mathrm{gap}^+(s,a)\in[\mathrm{gap}^+_{\min}, 1/(1-\gamma))$ and this interval can be separated into $N$ intervals \[\Lambda_n=[2^{n-1}\mathrm{gap}^+_{\min}, 2^n\mathrm{gap}^+_{\min})\]
where $n=1,2,\ldots,N$. Here, $N=\left\lceil\log_2\left(\frac{1}{\mathrm{gap}^+_{\min}(1-\gamma)}\right)\right\rceil$. Under the event $\mathcal E_{\mathrm{conc}}$, for $n\in[N]$, we define:
\[C^{(n)}:=\left|\left\{t\in\mathbb{N}_{+}~: \left(\hat{Q}_t-\breve{Q}_t\right)(s_t,a_t)\in\Lambda_n \right\}\right|.\]
By using the sequence
\[\omega_t^{(n)}:=\mathbb I\left[\left(\hat{Q}_t-\breve{Q}_t\right)(s_t,a_t)\in \Lambda_n\right],\]
we can upper bound the $C^{(n)}$ by using Lemma \ref{lemma-4}.
\begin{Lemma}
\label{lemma-5}
We can upper bound the $C^{(n)}$ by:
\[\mathcal O\left(\frac{SA}{4^n(\mathrm{gap}^+_{\min})^2(1-\gamma)^4\log(1/\gamma)}\log\left(\frac{SAT}{(1-\gamma)\mathrm{gap}^+_{\min}}\right)\right).\]
\end{Lemma}

Finally, we come to our main theorem. According to Equation (\ref{eqn-8}), we know that: if the trajectory satisfies the $\mathcal E_{\mathrm{conc}}$ condition, then:
\begin{align}
\label{eqn-9}
&~~\mathrm{Regret}(T)\leqslant \frac{1}{1-\gamma}\sum_{t=1}^{+\infty}\mathrm{clip}\left[\left(\hat{Q}_t-\breve{Q}_t\right)(s_t,a_t)|\mathrm{gap}^+_{\min}\right]\leqslant \frac{1}{1-\gamma}\sum_{n=1}^{N}2^n\mathrm{gap}^+_{\min}C^{(n)}\notag\\
&\leqslant \sum_{n=1}^{N} \mathcal O\left(\frac{SA}{2^n \mathrm{gap}^+_{\min}(1-\gamma)^5\log(1/\gamma)}\cdot\iota\right)= \mathcal O\left(\frac{SA}{\mathrm{gap}^+_{\min}(1-\gamma)^5\log(1/\gamma)}\cdot \iota\right).
\end{align}
Here, $\iota=\log\left(\frac{SAT}{\mathrm{gap}^+_{\min}(1-\gamma)}\right)$ is the logarithmic term on $T$. For the other trajectories outside $\mathcal E_{\mathrm{conc}}$, we have a trivial upper bound:
\begin{equation}
\label{eqn-10}
\mathrm{Regret}(T)\leqslant \sum_{t=1}^{T}\sum_{h'=t}^{+\infty}\gamma^{h'-t}\cdot\left(\hat{Q}_{h'}-\breve{Q}_{h'}\right)(s_{h'},a_{h'})\leqslant \sum_{t=1}^{T}\sum_{h'=t}^{+\infty}\frac{\gamma^{h'-t}}{1-\gamma}\leqslant \frac{T}{(1-\gamma)^2}.
\end{equation}
Now we combine Equation (\ref{eqn-9}) and Equation (\ref{eqn-10}), we obtain that:
\begin{align}
\label{eqn-11}
&\mathbb{E}[\mathrm{Regret}(T)]=\mathbb{E}\left[\sum_{t=1}^{T}\sum_{h'=t}^{+\infty}\gamma^{h'-t}\mathrm{gap}(s_{h'},a_{h'})\right]\notag\\
&\leqslant\mathbb{P}(\overline{\mathcal E_{\mathrm{conc}}})\cdot\frac{T}{(1-\gamma)^2}+\mathbb{P}(\mathcal E_{\mathrm{conc}})\cdot \mathcal O\left(\frac{SA}{\mathrm{gap}^+_{\min}(1-\gamma)^5\log(1/\gamma)}\cdot\iota\right)\notag\\
&\leqslant \mathcal O\left(\frac{SA}{\mathrm{gap}^+_{\min}(1-\gamma)^5\log(1/\gamma)}\cdot\iota\right).
\end{align}
which comes from $\mathbb{P}(\overline{\mathcal E_{\mathrm{conc}}})\leqslant 1/T$. Theorem \ref{thm-main-2} is proved and it provides us an upper bound which is logarithmically dependent on $T$.

\section{Episodic 2-TBSG with Linear Function Approximation}
In this section, we analyze the gap dependent total regret bound of episodic 2-TBSG under the linear function expression assumption. Here, we have two different settings, centralized version and independent version. As follows, we are going to introduce these two settings and our corresponding algorithms one by one. After that, gap dependent regret bound will be provided in both settings.

\subsection{Centralized Setting and Algorithm}
\begin{Assumption}[Assumption 4.1 from \cite{he2020logarithmic}]
\label{assumption-linear}
A Markov Game $\mathrm{MG}(2H,\mathcal S,\mathcal A,\mathbb{P},r)$ is defined as linear when the probability transition kernels and the reward functions are linear with respect to a given feature map $\phi:\mathcal S\times\mathcal A\rightarrow \mathbb{R}^d$ where $d$ is the feature dimension. Specifically, for each $h\in[2H]$, there exists an unknown vector $\mu_h\in\mathbb{R}^d$ and unknown measures $\theta_h=(\theta_h^{(1)},\theta_h^{(2)}, \ldots,\theta_h^{(d)})$ whose degree of freedom is $|S|\times d$, such that for $\forall (s,a)\in\mathcal S\times\mathcal A$:
\[\mathbb{P}_h(s'|s,a)=\langle\phi(s,a),\theta_h(s')\rangle~\text{and}~r(s,a)=\langle\phi(s,a), \mu_h\rangle.\]
For simplicity, we assume that $\|\phi(s,a)\|_2\leqslant 1, \|\mu_h\|_2\leqslant \sqrt{d}$ and $\|\theta_h(\mathcal S)\|\leqslant \sqrt{d}$.
\end{Assumption}

In the centralized setting, a central controller controls both players, and this central controller's goal is to learn a Nash Equilibrium. In our algorithm, both the max-player and the min-player update their policies $(\pi^k,\mu^k)$ according to the history information. Under Assumption \ref{assumption-linear}, we know that for any policy $\pi$, the action-value function $Q_h^{\pi}(s,a)$ is a linear function $\langle\phi(s,a), \theta_h^{\pi}\rangle$ with respect to the feature $\phi(s,a)$. In order to estimate the optimal action-value function $Q_h^*(s,a):=\langle\phi(s,a), \theta_h^*\rangle$, we only have to estimate the parameters $\theta_h^*$. Since Assumption \ref{assumption-linear} only gives a condition on the linear structure of the stochastic game, it is still a two-player turn-based general sum stochastic game with finite horizon. Therefore, we have exactly the same definition on the cumulative regret (or total regret) as Section 4:
\[\mathrm{Regret}(K)=\sum_{k=1}^{K}\left|\left(V_1^*-V_1^{\pi^k,\mu^k}\right)(s_1^k)\right|.\]

In the following Least Square Value Iteration on 2-TBSG (LSVI-2TBSG) algorithm, we introduce two new variables $\overline{w}_h^k, \underline{w}_h^k$, which are the upper and lower estimations of $\theta_h^*$ in the $k$-th episode. They are computed by solving the following regularized least-square problems:
\begin{equation}
\begin{aligned}
\overline{w}_h^k &\leftarrow \arg\min_{w\in\mathbb{R}^d} \lambda\|w\|^2 + F_1(w),\\
\underline{w}_h^k &\leftarrow \arg\min_{w\in\mathbb{R}^d} \lambda\|w\|^2 + F_2(w),
\end{aligned}
\end{equation}
where
\[F_1(w)=\sum_{i=1}^{k-1}\left[\phi(s_h^i, a_h^i)^{\top}w-r_h(s_h^i,a_h^i)-\overline{V}_{h+1}^k(s_{h+1}^i)\right]^2,\]
\[F_2(w)=\sum_{i=1}^{k-1}\left[\phi(s_h^i, a_h^i)^{\top}w-r_h(s_h^i,a_h^i)-\underline{V}_{h+1}^k(s_{h+1}^i)\right]^2.\]
Actually, these two least-square problems can be easily solved as:
\begin{equation*}
\begin{aligned}
\overline{w}_h^k&=(\Lambda_h^k)^{-1}\sum_{i=1}^{k-1}\phi(s_h^i,a_h^i)\left[r_h(s_h^i,a_h^i)+\overline{V}_{h+1}^k(s_{h+1}^i)\right],\\
\underline{w}_h^k&=(\Lambda_h^k)^{-1}\sum_{i=1}^{k-1}\phi(s_h^i,a_h^i)\left[r_h(s_h^i,a_h^i)+\underline{V}_{h+1}^k(s_{h+1}^i)\right],
\end{aligned}
\end{equation*}
where $\Lambda_h^k=\sum_{i=1}^{k-1}\phi(s_h^i,a_h^i)\phi(s_h^i,a_h^i)^{\top}+\lambda I$. Then, we update the estimated Q-values by:
\begin{equation*}
\begin{aligned}
\overline Q_h^k(s,a)&= \min(2H, \phi(s,a)^{\top}\overline{w}_h^k + \beta\cdot T(s,a)),\\
\underline Q_h^k(s,a)&= \max(0, \phi(s,a)^{\top}\underline{w}_h^k - \beta\cdot T(s,a)),
\end{aligned}
\end{equation*}
where $T(s,a)=\sqrt{\phi(s,a)^{\top}(\Lambda_h^k)^{-1}\phi(s,a)}$ can be regarded as a UCB term. For the pseudo-code of LSVI-2TBSG algorithm, we leave it to the appendix.

\subsection{Main Theorem for the Centralized Setting}
In this section, we propose the main theorem under the linear function expression assumption in the centralized setting. When using LSVI-2TBSG algorithm, the expected total regret can be upper bounded by the following theorem:
\begin{Theorem}[Main Theorem 3: Logarithmic Regret Bound of LSVI-2TBSG (Centralized)]
Under Assumption \ref{assumption-linear}, after using LSVI-2TBSG algorithm (centralized)
\[\mathbb{E}[\mathrm{Regret}(K)]\leqslant 1+\frac{Cd^3H^5\log(16dK^2(K+1)H^3)}{\mathrm{gap}^+_{\min}}\iota,\]
where $\iota=\log\left(\frac{Cd^3H^4\log(4dKH)}{(\mathrm{gap}^+_{\min})^2}\right)$ is a logarithmic term.
\label{thm-main-3}
\end{Theorem}
By using similar techniques, we can prove the theorem and propose the first gap-dependent logarithmic regret bound under the linear MDP assumption. We leave the technical proof of this theorem to the appendix.

\subsection{Independent Setting and Algorithm}
In the independent setting, we do not have a central controller who controls both players. We can only control the max-player and play against the min-player whose policies are arbitrary but potentially adversarial. Since we only control the max-player, our goal is not to learn a Nash Equilibrium, but to maximize the reward of the max-player. Because of the differences between the centralized setting and the independent setting, we are going to redefine the gaps and regret functions in this section.

Since we can not get access to the min-player's policies and the Markov model of the game a priori, we are interested in the exploitability of max-player:
\[\mathrm{Explicit}(\pi^k, \mu^k):=V_1^{\dagger, \mu^k}(s_1^k)-V_1^{\pi^k, \mu^k}(s_1^k),\]
which measures how much better the max-player can perform. Then, the cumulative regret can be defined as the sum of exploitability in different episodes:
\[\mathrm{Regret}_{\mu}(K):=\sum_{k=1}^{K}\mathrm{Explicit}(\pi^k, \mu^k)=\sum_{k=1}^{K}\left(V_1^{\dagger, \mu^k}(s_1^k)-V_1^{\pi^k, \mu^k}(s_1^k)\right).\]

Also, we need to redefine the gap. Previously, the gap is defined as $\mathrm{gap}_h(s,a)=|V_h^*(s)-Q^*_h(s,a)|$ and $\mathrm{gap}_{\min}^+:=\min_{h,s,a}\{\mathrm{gap}_h(s,a)>0\}$. However, in the independent setting, we can not control the min-player so it is not suitable to only consider $(\pi,\mu)=(\pi^*,\mu^*)$ since the Nash Equilibrium point is our final target. In order to measure the gap caused by max-player, we define:
\[\mathrm{gap}^{\mu}_h(s,a)=|V_h^{\dagger, \mu}(s)-Q_h^{\dagger, \mu}(s,a)|.\]
Here, we still need the absolute value since $V_h^{\dagger, \mu}(s)-Q_h^{\dagger, \mu}(s,a)$ is non-negative when $h$ is an odd number (and it's the max-player's turn to take action) while non-positive when $h$ is an even number. The minimal gap can be obtained after taking minimum over all the $(h,s,a)$ tuples and all possible pure strategy $\mu$:
\[\mathrm{gap}_{\min}^+:=\min_{\mu,h,s,a}\{\mathrm{gap}^{\mu}_h(s,a)>0\}.\]
Notice that the total number of pure strategies is finite, so the minimal gap above is positive and well-defined. Similar as the LSVI-2TBSG algorithm, we introduce a new variable $w_h^k$, which is the upper estimation of $\theta_h^*$ in the $k$-th episode. Here, we do not need the lower estimation $\underline{w}_h^k$ since the min-player's policy is beyond our control. In each episode, the $w^k_h$ is computed by solving the following regularized least-square problem:
\[w_h^k\leftarrow \arg\min_{w\in\mathbb{R}^d} \lambda\|w\|^2 + \sum_{i=1}^{k-1}\left[\phi(s_h^i,a_h^i)^{\top}w-r_h(s_h^i,a_h^i)-\overline{V}_{h+1}^k(s_{h+1}^i)\right]^2. \]
Actually, it can be solved as:
\[w_h^k=(\Lambda_h^k)^{-1}\sum_{i=1}^{k-1}\phi(s_h^i,a_h^i)\left[r_h(s_h^i,a_h^i)+\overline{V}_{h+1}^k(s_{h+1}^i)\right], \]
where $\Lambda_h^k=\sum_{i=1}^{k-1}\phi(s_h^i,a_h^i)\phi(s_h^i,a_h^i)^{\top}+\lambda I$. Then, we update the estimated Q-values by:
\[\overline Q_h^k(s,a)= \min(2H, \phi(s,a)^{\top}\overline{w}_h^k + \beta\cdot T(s,a)),\]
where $T(s,a)=\sqrt{\phi(s,a)^{\top}(\Lambda_h^k)^{-1}\phi(s,a)}$ is the UCB term. We leave the pseudo-code of the algorithm to the appendix.

\subsection{Main Theorem for the Independent Setting}
In this section, we propose the main theorem under the linear function expression assumption in the independent setting. The expected cumulative regret can be upper bounded by a logarithmic term:
\begin{Theorem}[Main Theorem 4: Logarithmic Regret Bound of LSVI-2TBSG (Independent)]
Under Assumption \ref{assumption-linear}, after using LSVI-2TBSG algorithm (independent)
\[\mathbb{E}[\mathrm{Regret}(K)]\leqslant 1+\frac{Cd^3H^5\log(16dK^2(K+1)H^3)}{\mathrm{gap}^+_{\min}}\iota,\]
where $\iota=\log\left(\frac{Cd^3H^4\log(4dKH)}{(\mathrm{gap}^+_{\min})^2}\right)$ is a logarithmic term.
\label{thm-main-4}
\end{Theorem}

We provide a logarithmic regret upper bound for 2-TBSG under the linear MDP assumption, in both centralized and independent settings. To the best of our knowledge, this is the very first gap dependent upper bound for 2-TBSG, which makes our results novel and complete.

\section{Conclusion}
We gave the first set of gap-dependent logarithmic regret bounds for two-player turn-based stochastic Markov games in both tabular and episodic cases, and in both centralized setting and independent setting.
One fruitful future direction is to extend our analysis to more general settings~\cite{wang2020reinforcement,jiang2017contextual,du2021bilinear,jin2021bellman}.