\section{Introduction}

	The explosive growth of distributed \jpang{generation} in distribution networks together with the advancement of communication and control technology at the consumer level have gradually transformed \jpang{the traditionally} passive consumers into ``prosumers", which can both produce and consume energy \cite{chen2019energy}.
	Then, energy trading with the main grid or between prosumers is inevitable \jpang{since energy surplus and shortage are bound to exist} \cite{Wang2021Distributed}. In this situation, the peer-to-peer (P2P) market, \jpang{which operates} in a distributed manner, is more popular due to the ever-increasing number of prosumers, \jpang{in which} the various prosumers can be self-organized to operate economically and reliably under a given market mechanism \cite{Morstyn2018peer-to-peer}. In addition, the \jpang{increasing penetration and an aggravating volatility of} renewable generation calls for online market clearing methods. In this paper, we intend to investigate the distributed online energy trading market \jpang{for prosumers}.

	For such P2P energy trading markets, they are usually formulated as generalized Nash games, where each prosumer maximizes its profit with coupling constraints, \jpang{e.g.,} global power balance \cite{chen2019energy,Wang2021Distributed,Cui2020FairPeer-to-Peer,Nespoli2018rational,Chen2021SameFlexibility,Belgioioso2022Operationally}. Then, clearing the resulting P2P market corresponds to finding the generalized Nash equilibrium (GNE) of the energy trading game. For example, in \cite{chen2019energy}, the energy sharing game among prosumers is formulated with full information, and \cite{Wang2021Distributed} further designs a fully distributed algorithm based on Nesterov's methods to seek the GNE with only partial-decision information. In \cite{Cui2020FairPeer-to-Peer}, a P2P energy market is formulated as a generalized Nash game, where the \jpang{prosumers who share payments} are mutually coupled and influenced. \jpang{Following this,} \cite{Nespoli2018rational} and \cite{Chen2021SameFlexibility} further consider system-level grid constraints. \jpang{Lastly, in} \cite{Belgioioso2022Operationally}, a P2P energy market of prosumers is formulated as a generalized aggregative game with global coupling constraints.
	The aforementioned works have made great progress in the distributed GNE seeking for the P2P energy trading market. However, they usually focus on \jpang{only} one time section and provide offline solutions to solve the game. Due to the volatility of renewable generations and the complexity of load profiles, \jpang{both current and future operation status changes much more over time, requiring much faster algorithms}, i.e., online GNE tracking.

	In this paper, we formulate a P2P energy trading market among prosumers in the distribution network and propose a distributed online algorithm to track the GNE of the market. The major contributions are as follows.
	\begin{itemize}
		\item A P2P energy trading market is modeled as a generalized Nash game with both individual and \jpang{coupled} time-varying constraints. Moreover, we prove the uniqueness of the GNE of this market at any time section.
		\item A novel distributed online algorithm is proposed to track the GNE, where each prosumer can make decisions only using local variables and neighboring information. This reduces the communication burden and makes it easier to implement in practice.
		\item \jpang{We prove a sublinear regret bound, i.e.,} that the regret of the online algorithm can be bounded by a sublinear function of learning time, indicating that the online algorithm suffers minimal ``loss in hindsight''.
	\end{itemize}

	The rest of this paper is organized as follows. In Section \uppercase\expandafter{\romannumeral2}, the P2P energy trading game is formulated. \jpang{Section \uppercase\expandafter{\romannumeral3} introduces and analyzes the performance} of a distributed online algorithm to track the GNE of the game in a time-varying environment.

	Numerical results are presented in Section \uppercase\expandafter{\romannumeral4} to verify the effectiveness of our algorithm. Finally, Section \uppercase\expandafter{\romannumeral5} concludes the paper.

	\textit{Notations}:
	In this paper, $\mathbb{R}_{+}^n$ is the $n$-dimensional (nonpositive) Euclidean space. For a column vector $x\in \mathbb{R}^n$ (matrix $A_{m\times n}\in \mathbb{R}^{m\times n}$), its transpose is denoted by $x^\mathsf{T}$($A^\mathsf{T}$). For a matrix $A$, $ [A]_{i,j}$ stands for the entry in the $i$-th row and $j$-th column of A.
	For vectors $x,y\in \mathbb{R}^n$, $x^\mathsf{T}y=\left\langle x,y \right\rangle$ is the inner product of $x,y$, while $\otimes$ represents the Kronecker product.
	$\left\|x \right\|=\sqrt{x^\mathsf{T}x}$ is the Euclidean norm.
	The identity matrix with dimension $n$ is denoted by $I_n$. Sometimes, we also omit $n$ to represent the identity matrix with the proper dimension. $\mathbf{0}_{n}, \mathbf{1}_{n}$ are all zero and all one vectors with dimension $n$, respectively. The Cartesian product of the sets $\Omega_i, i=1, \cdots, n$ is denoted by $\prod_{i=1}^n\Omega_i$.
	Given a collection of $y_i$ for $i$ in a certain set $Y$, the vector composed of $y_i$ is defined as ${y}={col}(y_{i}):=(y_1, y_2, \cdots, y_n)^\mathsf{T}$. The projection of $x$ onto a set $\Omega$ is defined as $  	\mathcal{P}_{\Omega}(x):=\arg \min_{y\in \Omega}\left\|x-y \right\|$.

	\section{Problem Formulation}

	\subsection{Network model}

	We consider a distribution network with a group of prosumers, denoted by the set $\mathcal{N}=\{1,2,...,N\}$. For each prosumer, its load demand can be satisfied by its own generation and trading with the main grid or its neighboring prosumers. The trading edge is denoted by $\mathcal{E}\subseteq \mathcal{N}\times \mathcal{N}$. For a prosumer $ i $, the set of its neighbors is denoted by $ \mathcal{N}_i =\{\mathcal{N}_1^1,\dots,\mathcal{N}_1^{N_i}\}$ with $\left|\mathcal{N}_i\right|=N_i$. If $ j\in \mathcal{N}_i $, prosumers $ i $ and $ j $ can trade and communicate directly. Otherwise, direct trading and communication are not allowed. Then, the trading network is modeled as an undirected graph $\mathcal{G}=(\mathcal{N},\mathcal{E})$. The adjacency matrix of $\mathcal{G}$ is denoted by $W$ with elements $w_{i,j}$. If $j\in \mathcal{N}_i$, the weight $w_{i,j}$ satisfies $w_{i,j}=w_{j,i}>0$. Otherwise, $w_{i,j}=w_{j,i}=0$.  The Laplacian matrix of the communication graph is denoted by $L$ and we have $\textbf{1}^{\rm T}L=0$, where $ \textbf{1} $ is an all-one vector.
	Moreover, the graph $\mathcal{G}$ is assumed to be connected.

	For the weights, we have the following assumption, which implies that every row sum of $W$ is identical.
	\begin{assumption}\label{assumption_weight}
		The weight $w_{i,i}>0$ and $\sum_{j\in \mathcal{N}}w_{i,j}=w_0>0$ for all $i\in \mathcal{N}$.
	\end{assumption}

	\subsection{Prosumer model }
	The scenario is that each prosumer is equipped with dispatchable generation, a non-dispatchable load, and an energy storage system (ESS). To maintain power balance, it can generate electricity, \jpang{charge or discharge from} the ESS, \jpang{and/or} trade with the main grid or neighboring prosumers. In this paper, we focus on the time horizon $\mathcal{T}=\{1,2,...,T\}$. Here, we will introduce them in detail.

	\emph{Dispatchable generation:} The power generated by dispatchable generation units of prosumer $i$ at time $t$, denoted by $p_i^g(t)$, is limited by
	\begin{align}\label{generation_limit}
		p_i^{g,min}\leq p_i^g(t)\leq p_i^{g,max},\quad \forall i\in \mathcal{N},t \in \mathcal{T}
	\end{align}
	where $p_i^{g,min}$ and $p_i^{g,max}$ are minimum and maximum local generation, respectively. Its generation cost is as follows.
	\begin{align}\label{generation_cost}
		f_i^g\left(p_i^g(t)\right)=a_i^g \left(p_i^g(t)\right)^2+b_i^gp_i^g(t)
	\end{align}
	where $a_i^g>0$ and $b_i^g$ are constants.

	\emph{Energy Storage Systems (ESS):} The ESS profile is constrained by the following dynamics.
	\begin{align}
		&0\leq p_i^c(t)\leq p_i^{c,max},\quad \forall i\in \mathcal{N},t \in \mathcal{T}\label{charge_limit}\\
		&0\leq p_i^d(t)\leq p_i^{d,max},\quad \forall i\in \mathcal{N},t \in \mathcal{T}\label{discharge_limit}\\
		&s_i(t+1)=s_i(t)+\frac{\Delta t}{e_i^{cap}}\left(\eta_i^c p_i^c(t)-\frac{1}{\eta_i^d}p_i^d(t)\right)\label{soc}\\
		&s_i^{min}\leq s_i(t)\leq s_i^{max},\quad \forall t\in \mathcal{T}\label{soc_limit}
	\end{align}
	where $p_i^c(t)$, $p_i^d(t)$, and $s_i(t)$ are the charging, discharging power, and state of charge (SoC) of the ESS $i$, respectively. $\Delta t$, $e_i^{cap}$, $\eta_i^c$ and $\eta_i^d$ are sampling time, ESS maximum storage capacity, and (dis)charging efficiencies, respectively. Moreover, $s_i^{min}$ and $s_i^{max}$, with $0<s_i^{min}<s_i^{max}<1$, denote the minimum and maximum SoC, while $p_i^{c,max}$ and $p_i^{d,max}$ denote the maximum (dis)charging power.

	Each prosumer might also minimize the usage of its ESS to reduce its degradation. Depending on the efficiency of a storage unit, there are losses based on usage that usually grow quadratically in power. \jpang{For simplicity, we disregard the effects due to SoC levels.} As defined in \cite{Hans2019Hierarchical}, the corresponding cost function is
	\begin{align}\label{es_cost}
		f_i^{es}\left(p_i^c(t),p_i^d(t)\right)=a_i^c \left(p_i^c(t)\right)^2+a_i^d \left(p_i^d(t)\right)^2
	\end{align}
	where $a_i^c$ and $a_i^d$ are both positive constants.

	\emph{Trading with the main grid:}  Let $p_i^{mg}(t)$ be the power purchased from the main grid at time $t$ and $p^{mg}(t)=col\left\{p_i^{mg}(t)\right\}_{i\in\mathcal{N}}$. Similar to \cite{Atzeni2013Demand-Side}, \jpang{we set grid cost as}
	\begin{align}\label{grid_cost}
		C_t(p^{mg}(t))=c_t^{mg}\left(\sum\nolimits_{i\in \mathcal{N}}p_i^{mg}(t)\right)^2
	\end{align}
	where $c_t^{mg}$ is a \jpang{time-varying cost coefficient}, since the energy production varies along the time period according to the energy demand and the availability of distributed energy sources. Then the cost assigned to prosumer $i$ is
	\begin{align}\label{mg_cost}
		f_i^{mg}(p^{mg}(t))&=\frac{p_i^{mg}(t)}{\sum_{i\in \mathcal{N}}p_i^{mg}(t)}C_t(p^{mg}(t))\nonumber\\
		&=c_t^{mg}p_i^{mg}(t)\sum\nolimits_{i\in \mathcal{N}}p_i^{mg}(t)
	\end{align}
	Moreover, the total power exchanged with the main grid is limited, i.e.,
	\begin{align}\label{mg_coupling}
		p^{mg,min}\leq \sum\nolimits_{i\in \mathcal{N}}p_i^{mg}(t)\leq p^{mg,max},\quad \forall t \in \mathcal{T}
	\end{align}

	\emph{Trading with neighbors:} The trading cost with neighboring prosumers \jpang{of prosumer $i$} is
	\begin{align}\label{trade_cost}
		f_i^{tr}\left(\left\{p_{i,j}^{tr}(t)\right\}\right)=\sum_{j\in \mathcal{N}_i}\left[a^{tr}\left(p_{i,j}^{tr}(t)\right)^2+d_{i,j}^{tr}p_{i,j}^{tr}(t)\right]
	\end{align}
	where $p_{i,j}^{tr}(t)$ is the power purchased from prosumer $j$ at time $t$, $d_{i,j}^{tr}=d_{j,i}^{tr}>0$ is the price and $a^{tr}$ is a small positive constant, which represents the tax cost incurred by using the energy sharing platform.

	Disregarding loss on the power lines, the sum of the trading power of prosumer $i$ and $j$ at time $t$ should be 0.
	\begin{align}\label{trade_coupling}
		p_{i,j}^{tr}(t)+p_{j,i}^{tr}(t)=0,\quad \forall (i,j)\in \mathcal{E},t \in \mathcal{T}
	\end{align}
	Furthermore, trade between prosumers is limited by
	\begin{align}\label{trade_limit}
		p_{i,j}^{tr,min}\leq p_{i,j}^{tr}(t)\leq p_{i,j}^{tr,max},\quad \forall (i,j)\in \mathcal{E},t \in \mathcal{T}
	\end{align}
	where $p_{i,j}^{tr,min}\leq 0$ and $p_{i,j}^{tr,max}\geq 0$ are the minimum and maximum \jpang{tradeable power between prosumers $i$ and $j$}.

	Denote by $p_i^{l}(t)$ the undispatchable load demand, and the local power balance for each prosumer $i$ is
	\begin{align}\label{power_balance}
		&p_i^g(t)-p_i^c(t)+p_i^d(t)+p_i^{mg}(t)+\sum\nolimits_{j\in \mathcal{N}_i}p_{i,j}^{tr}(t)=p_i^{l}(t)
	\end{align}
	\subsection{Energy trading game}
	Before giving the game model, we first simplify notations. The decision variable of prosumer $i$ is denoted by
	\begin{align*}
		x_i(t):=col\left(p_i^g(t),p_i^c(t),p_i^d(t),p_i^{mg}(t),{p_{i,j}^{tr}(t)}\right)\in\mathbb{R}^{n_i}
	\end{align*}
	where $n_i=4+N_i$ and $\sum\nolimits_{i\in\mathcal{N}}n_i=n$.

	Define a sparse matrix $E_i$, where the rows of $E_i$ correspond to every trading edge in $\mathcal{E}$ one by one. Let the $k$-th row of $E_i$ corresponds to $\left(I_k,J_k\right)$ in $\mathcal{E}$, then the elements $[E_i]_{k,l}$ of $E_i$ are assigned as follows.
	\begin{align*}
		[E_i]_{k,l}=\left\{
		\begin{matrix}
			1,&\text{If } \left\{I_k,J_k\right\}=\left\{i,\mathcal{N}_i^l\right\} \text{ and } I_k<J_k. \\
			-1,&\text{If } \left\{I_k,J_k\right\}=\left\{i,\mathcal{N}_i^l\right\} \text{ and } I_k>J_k.\\
			0,&\text{Otherwise.}
		\end{matrix}
		\right.
	\end{align*}

	Let
	\begin{align*}
		A_i&=\begin{bmatrix}
			\mathbf{0}_{m\times 3}&
			\begin{matrix}
				-1&\mathbf{0}^\mathsf{T}_{N_i}\\
				1&\mathbf{0}^\mathsf{T}_{N_i}\\
				\mathbf{0}_{\sum_{i\in \mathcal{N}}N_i}&E_i
			\end{matrix}
		\end{bmatrix}\\
		b_i&=\begin{bmatrix}
			-\frac{p^{mg,min}}{N}&
			\frac{p^{mg,min}}{N}&
			\mathbf{0}^\mathsf{T}_{\sum_{i\in \mathcal{N}}N_i}
		\end{bmatrix}^\mathsf{T}\\
		g_i(x_i(t))&=A_ix_i(t)-b_i,\quad \forall i\in \mathcal{N},t \in \mathcal{T}
	\end{align*}
	where $m=2+\sum_{i\in \mathcal{N}}N_i$.
	Then coupling constraints \eqref{mg_coupling} and \eqref{trade_coupling} can be reformulated as
	\begin{align}\label{coupling_constraint}
		\sum\nolimits_{j\in \mathcal{N}}g_i(x_i(t))\leq 0
	\end{align}

	In \eqref{coupling_constraint}, the sparse matrix $E_i$ is designed to address the equality constraint \eqref{trade_coupling} by transforming it into two equivalent inequality constraints.

	Similarly, let
	\begin{align*}
		G_i=\begin{bmatrix}
			1&-1&1&1&\mathbf{1}^\mathsf{T}_{N_i}
		\end{bmatrix}
	\end{align*}
	Then equality constraint $\eqref{power_balance}$ can be reformulated as
	\begin{align}\label{equality_constraint}
		G_ix_i(t)-p_i^{l}(t)=0,\quad \forall i\in \mathcal{N},t \in \mathcal{T}
	\end{align}

	Moreover, the domain set of $x_i(t)$ is denoted by
	\begin{align*}
		&\chi_i=\left\{x_i(t)|x_i(t)\ \text{satisfies}\  \eqref{generation_limit},\eqref{charge_limit},\eqref{discharge_limit},\eqref{soc},\eqref{soc_limit},\eqref{trade_limit}\right\}\\
		&\Omega_i=\chi_i\cap \left\{x_i(t)|x_i(t)\text{ satisfies }\eqref{equality_constraint}\right\}\\
		&X=\prod\nolimits_{i\in\mathcal{N}}\Omega_i\cap \left\{x(t)|x(t)\text{ satisfies }\eqref{coupling_constraint}\right\}
	\end{align*}
	where $\chi_i$ is the set of all of the local inequality constraints, $\Omega_i$ considers the \jpang{reformulated power balance constraint of \eqref{power_balance}}, while $X$ includes coupling constraints.

	In the energy trading game, each prosumer intends to minimize its cost while maintaining the global power balance. The optimization problem of each prosumer is
	\begin{subequations}\label{Energy_trading_game}
		\begin{align}\label{Energy_trading_game_obj}
			\min_{x_i(t)} \ &J_{i,t}\left(x_i(t),x_{-i}(t)\right)=f_i^g\left(p_i^g(t)\right)+f_i^{es}\left(p_i^c(t),p_i^d(t)\right)\nonumber\\
			&\qquad+f_i^{mg}\left(p^{mg}(t)\right)+f_i^{tr}\left(\left\{p_{i,j}^{tr}(t)\right\}\right)\\
			\text{ s.t.}\ & x_i(t)\in X_i(x_{-i}(t))
		\end{align}
	\end{subequations}
	where $x_{-i}(t):=\operatorname{col}\left(x_1(t),...,x_{i-1}(t),x_{i+1}(t),...,x_N(t)\right)$ and $X_i(x_{-i}(t)):=\left\{x_i(t)\in \Omega_i| (x_i(t),x_{-i}(t))\in X\right\}$.

	In summary, the energy trading game is represented as
	\begin{itemize}
		\item Player: all prosumers, denoted by $\mathcal{N}=\{1,2,...,N\}$.
		\item Strategy: decision variable $x_i$.
		\item Payoff: the disutility function $J_{i,t}\left(x_i(t),x_{-i}(t)\right)$.
	\end{itemize}
	Due to the global coupling constraints \eqref{coupling_constraint}, it is a generalized Nash game. The corresponding GNE is defined as
	\begin{align*}
		x_i^*(t)&\in \arg \min\; J_{i,t}\left(x_i(t),x_{-i}^*(t)\right),  \text{s.t. } x_i(t)\in X_i\left(x_{-i}^*(t)\right)
	\end{align*}

	Regarding the game \eqref{Energy_trading_game}, we have following assumptions.
	\begin{assumption}\label{assumption_domain_set}
		$\chi_i$ is a  non-empty, compact and convex set.
	\end{assumption}
	\begin{assumption}\label{assumption_feasibility}
		Given any $x_{-i}(t)$, problem \eqref{Energy_trading_game} is feasible.
	\end{assumption}
	Since the constraints of problem \eqref{Energy_trading_game} are all affine, the commonly assumed Slater's condition is simplified to only require feasibility, and therefore Assumption \ref{assumption_feasibility} suffice.

	\subsection{Uniqueness of the GNE}
	The pseudo-gradient of $\{J_{i,t}\}_{i\in\mathcal{N}}$ is defined as
	\begin{align}\label{pseudo_gradient_definition}
		&F_t(x)=col\left(\left\{\nabla_{x_i}J_{i,t}(x_i,x_{-i})\right\}\right)\nonumber\\
		&=col\left(\left\{\begin{matrix}
			2a_i^gp_i^g+b_i^g\\
			2a_i^cp_i^c\\
			2a_i^dp_i^d\\
			c_t^{mg}\left(2p_i^{mg}+\sum_{j\in \mathcal{N},j\neq i}p_j^{mg}\right)\\
			col\left(\left\{2a^{tr}p_{i,j}^{tr}+d_{i,j}^{tr}\right\}_{j\in \mathcal{N}_i}\right)
		\end{matrix}\right\}\right)
	\end{align}
	Define $\underline {a}^g=\underset{i\in\mathcal{N}}{\min}\,a_i^g$, $\overline {a}^g=\underset{i\in\mathcal{N}}{\max}\,a_i^g$, \jpang{with $\underline {a}^c$, $\overline {a}^c$, $\underline {a}^d$, $\overline {a}^d$ defined similarly}. In addition, $\underline {c}^{mg}=\underset{t\in \mathcal{T}}{\min}\, c_t^{mg}$, $\overline {c}^{mg}=\underset{t\in \mathcal{T}}{\max}\, c_t^{mg}$. Then, we can prove that the pseudo-gradient is strongly monotone and Lipschitz continuous.
	\begin{lemma}\label{lemma_monotonicity_continuity}
		For $\forall t\in \mathcal{T}$, the pseudo-gradient $F_t(x)$ has following properties\jpang{:}
		\begin{enumerate}[]
			\item $F_t(x)$ is strongly monotone with parameter $0<\eta\leq \min\left\{2\underline {a}^g,2\underline {a}^c,2\underline {a}^d, \underline {c}^{mg},2a^{tr}\right\}$, i.e., $\langle F_t(x)-F_t(y),x-y\rangle\geq \eta||x-y||_2^2$;
			\item  $F_t(x)$ is $\theta-$Lipschitz continuous, i.e., $\left\|F_t(x)-F_t(y)\right\|_2\leq \theta ||x-y||_2$ with $\theta\geq \max\left\{2\overline {a}^g,2\overline {a}^c,2\overline {a}^d,N\overline {c}^{mg},2a^{tr}\right\}$.
		\end{enumerate}
	\end{lemma}

	\begin{proof}
		For 1), taking any two variables $x^1, x^2$, then
		\begin{align}
			&\langle F_t(x^1)-F_t(x^2),x^1-x^2\rangle\nonumber\\
			&=\sum\nolimits_{i \in \mathcal{N}}2a_i^g\left\|p_i^{g,1}-p_i^{g,2}\right\|_2^2+\sum\nolimits_{i \in \mathcal{N}}2a_i^c\left\|p_i^{c,1}-p_i^{c,2}\right\|_2^2\nonumber\\
			&\quad+\sum\nolimits_{i \in \mathcal{N}}2a_i^d\left\|p_i^{d,1}-p_i^{d,2}\right\|_2^2+c_t^{mg}h^THh\nonumber\\
			&\quad+\sum\nolimits_{i \in \mathcal{N}}\sum\nolimits_{j \in \mathcal{N}_i}2a^{tr}\left\|p_{i,j}^{tr,1}-p_{i,j}^{tr,2}\right\|_2^2
		\end{align}
		where $h=col\left(p_i^{mg,1}-p_i^{mg,2}\right)$ and $H=I_N+\textbf{1}_{N\times N}$. Since the minimal eigenvalue of $H$ is $1$, we have $h^THh\geq \sum_{i \in \mathcal{N}}\left\|p_i^{mg,1}-p_i^{mg,2}\right\|_2^2$. Then $\langle F_t(x^1)-F_t(x^2),x^1-x^2\rangle\geq\eta\left\|x^1-x^2\right\|_2^2$ with $0<\eta\leq \min\left\{2\underline {a}^g,2\underline {a}^c,2\underline {a}^d, \underline {c}^{mg},2a^{tr}\right\}$. Therefore, $F_t(x)$ is strongly monotone with parameter $\eta$.

		Similarly, the second assertion could be obtained, which is omitted due to the space limit.
	\end{proof}
	A GNE with the same Lagrangian multipliers for all the agents is called variational GNE (v-GNE) \cite{Facchinei2010Generalized}, which has the economic interpretation of no price discrimination \cite{Kulkarni2012variational}. For $\forall t\in\mathcal{T}$,  every solution $x^*(t)\in X$ to the following variational inequality is a v-GNE of game $\eqref{Energy_trading_game}$.
	\begin{align}\label{variational _inequality}
		\langle F_t(x^*(t)),x-x^*(t)\rangle\leq 0,\qquad \forall x\in X
	\end{align}

	\begin{theorem}\label{lemma_equilibrium_existence}
		For the generalized Nash game $\eqref{Energy_trading_game}$, there exists a unique v-GNE.
	\end{theorem}
	\begin{proof}
		Following \cite{Facc2010Nashequilibria}, since $J_{i,t}(x_i(t),x_{-i}(t))$ is differentiable and convex with respect to $x_i(t)$ for any $x_{-i}(t)$, if {Assumption} {\ref{assumption_domain_set}} and {\ref{assumption_feasibility}} hold, the v-GNE of $\eqref{Energy_trading_game}$ exists. Moreover, by the strong monotonicity of $F_t(x)$, the uniqueness of v-GNE is guaranteed.
	\end{proof}

	\section{Online Algorithm}

	In this section, we first propose an online distributed algorithm based on a consensus algorithm and primal-dual strategy to solve the problem \eqref{Energy_trading_game}. Then, we prove that the regret of the proposed algorithm is bounded by a sublinear function of the learning time.
	\subsection{Algorithm design}
	Recalling the objective function \eqref{Energy_trading_game_obj}, $f_i^{mg}\left(p^{mg}(t)\right)$ is associated with decisions of all of the other prosumers. To solve game \eqref{Energy_trading_game}, full information is needed, i.e., one prosumer needs to communicate with all of the other prosumers. However, this is difficult to realize in practice due to communication limits. This section designs a distributed algorithm with only partial-decision information, where the prosumers only need to exchange information with its neighbors. To this end, we endow each prosumer with an auxiliary variable $\bar{x}^i(t)$ that provides an estimate of the decisions of other prosumers at time $t$. Moreover, \jpang{$\bar{x}_j^i(t)$ represents prosumer $i$'s estimate of $j$'s decision} and $\bar{x}^i_{-i}(t)=col\left(\left\{\bar{x}_j^i(t)\right\}_{j\in\mathcal{N}_i}\right)$.  Clearly, we have $\bar{x}_i^i(t)=x_i(t)$.

	\jpang{Firstly, note that the} Lagrangian of problem \eqref{Energy_trading_game} is
	\begin{align}\label{Lagrangian_function}
		&L_{i,t}\left(x_i(t),\lambda_i(t),\mu_i(t);x_{-i}(t)\right)=J_{i,t}\left(x_i(t),x_{-i}(t)\right)\nonumber\\
		&\ +\lambda_i^T(t)\sum\nolimits_{i \in \mathcal{N}}g_i(x_i(t))+\mu_i(t)(G_ix_i(t)-p_i^l(t))
	\end{align}
	where $\lambda_i(t)$ and $\mu_i(t)$ are Lagrange multipliers.

	The iterative process of $x_i(t)$, $\bar{x}_{-i}^i(t)$ and dual variables $\left(\lambda_i(t),\mu_i(t)\right)$ is shown in Algorithm \ref{alg1}, where $0<\rho(t)<1$ is a stepsize or the so-called learning rate, which decreases over $t$, and  $c=\frac{1}{w_0}$.
	\begin{algorithm}[t]
		\caption{Online Algorithm for P2P Energy Trading}
		\label{alg1}
		\begin{algorithmic}
			\State \textbf{Initialization:} $x_i(0)\in \chi_i$, $\bar{x}_{-i}^i(0)\in\mathbb{R}^{n-n_i}$, $\lambda_i(0)=0$,
			$\mu_i(0)=0$
			\For{$t=1$ to $T$}

			\begin{subequations}
				\begin{align}
					&x_i(t+1)=(1-\rho(t))x_i(t)\nonumber\\
					&\qquad+\rho(t)\mathcal{P}_{\chi_i}\left\{x_i(t)
					-\rho(t)[\nabla_{x_i}J_{i,t}(x_i(t),\bar{x}_{-i}^i(t))\right.\nonumber\\
					&\qquad+ \rho(t)(A_i^T\lambda_i(t)
					+G_i^T\mu_i(t))\nonumber\\
					&\qquad\left.+c\sum\nolimits_{j\in \mathcal{N}_i}(x_i(t)-\bar{x}_{i}^j(t))]\right\}\label{algorithm_xi}\\
					&\bar{x}_{-i}^i(t+1)=\bar{x}_{-i}^i(t)\nonumber\\
					&\qquad-c \rho(t) \sum\nolimits_{j \in \mathcal{N}_i}w_{i,j}\left(\bar{x}_{-i}^i(t)-\bar{x}_{-i}^j(t)\right)\label{algorithm_x-i}\\
					&\lambda_i(t+1)=\mathcal{P}_{\mathbb{R}_+^m}\bigg\{\left(1-\rho(t)\right) \sum\nolimits_{j \in \mathcal{N}_i\cup \{i\}} \frac{w_{i j}}{w_0} \lambda_j(t)\nonumber\\
					&\qquad+\rho(t)\left[A_i\left(2 x_i(t+1)-x_i(t)\right)-b_i\right] \bigg\}\label{algorithm_lambda}\\
					&\mu_i(t+1)=\left(1-\rho(t)\right) \mu_i(t)\nonumber\\
					&\qquad+\rho(t)\left[G_i\left(2 x_i(t+1)-x_i(t)\right)-p_i^l(t)\right]\label{algorithm_mu}
				\end{align}
			\end{subequations}
			\EndFor
		\end{algorithmic}
	\end{algorithm}

	The update for $x_i(t)$ in Algorithm \ref{alg1} employs the projected primal-dual gradient decomposition method combined with the consensus approach \cite{wang2019distributed2,wang2019distributed,wang2019distributed_variation}. The update of $\bar{x}_{-i}^i(t)$ can be regarded as a discrete-time integration for the consensus error of the local estimation \cite{pavel2019distributed}. At each sampling time $t$, prosumer $i$ gets $\nabla_{x_i}J_{i,t}\left(x_i(t),x_{-i}(t)\right)$ by using $x_i(t), \bar{x}_{-i}^i(t)$ and updates $x_i(t+1)$ with \eqref{algorithm_xi}. Meanwhile, the estimation $\bar{x}_{-i}^i(t+1)$ is updated by communication with neighboring prosumers by \eqref{algorithm_x-i}. Then, dual variables $\left(\lambda_i(t),\mu_i(t)\right)$ are updated using the latest updated local information $x_i(t+1)$ with \eqref{algorithm_lambda} and \eqref{algorithm_mu}, respectively.

	\begin{remark}
		Algorithm \ref{alg1} is fully distributed with only partial-decision information. Each prosumer makes a decision only based on local information and communication with its neighbors, which is easy to implement in practice.
		Compared with the existing work \cite{Kaihong2021Environments} only considering the time-varying objective function, we further include the time-varying constraints. Moreover, compared with the algorithm in \cite{2021arXiv210506200M}, the update is much simpler without the need of solving an optimization problem at each iteration, which reduces the computation cost.
	\end{remark}

	\subsection{Regret analysis}
	In this subsection, we will prove that the regret of Algorithm \ref{alg1} is bounded by a sublinear function of the learning time. First, we give some notations.
	Under {Assumption} {\ref{assumption_domain_set}},  $\left\|x_i\right\|$, $\left\|g_i(x_i)\right\|$, $\left\|\nabla_{x_i}J_{i,t}(x_i,x_{-i})\right\|$, $\left\|\nabla_{x_i}g_i(x_i)\right\|$ and  $\left\|G_ix_i-p_i^l(t)\right\|$ are bounded. Then, for $\forall i\in \mathcal{N}$ and $\forall t\in \mathcal{T}$, we define
	\begin{align*}
		\kappa_1&=\sup\nolimits_{x_i\in\chi_i}\left\|x_i\right\|,\ 			\kappa_2=\sup\nolimits_{x_i\in\chi_i}\left\|g_i(x_i)\right\|\\
		\kappa_3&=\sup\limits_{x_i\in\chi_i}\left\|\nabla_{x_i}J_{i,t}(x_i,x_{-i})\right\|,\
		\kappa_4=\sup\limits_{x_i\in\chi_i}\left\|\nabla_{x_i}g_i(x_i)\right\|\\
		\kappa_5&=\sup\limits_{x_i\in\chi_i}\left\|G_ix_i-p_i^l(t)\right\|,\
		\kappa_6=\max\limits_{i\in\mathcal{N}}\sup_{\left\|x_i\right\|=1}\left\|G_ix_i\right\|
	\end{align*}

	We use the dynamic regret to evaluate the performance of the online algorithm, which is defined as follows.
	\begin{align}\label{regret_definition}
		R_i(T)=\sum\nolimits_{t=1}^T\left[J_{i,t}\left(x_i(t), x_{-i}^*(t)\right)-J_{i,t}\left(x^*(t)\right)\right]
	\end{align}
	where $x^*(t)=col(x_i^*(t))$ is the v-GNE of \eqref{Energy_trading_game} at time $t$.

	An online algorithm is \jpang{generally considered to perform well} if the regret increases sublinearly \cite{2021arXiv210506200M}, \cite{Shahrampour2018MirrorDescent}, i.e.,
	\begin{align}\label{regret_limit}
		\lim _{T \rightarrow \infty} \frac{R_i(T)}{T}=0,\quad \forall i \in \mathcal{N}
	\end{align}

	However, it would be impossible to keep the regret \eqref{regret_definition} increasing sublinearly if the v-GNE sequence $\left\{x^*(0),x^*(1),\dots,x^*(N)\right\}$ of \eqref{Energy_trading_game} fluctuates drastically. Therefore, motivated by \cite{Hall2015Online} and \cite{Kaihong2020Pseudoconvex}, we adopt the following accumulation to describe the difficulty of tracking the v-GNE sequence:
	\begin{align}\label{varphi_T_definition}
		\varPhi_T=\sum\nolimits_{t=0}^T\left\|x^*(t+1)-x^*(t)\right\|
	\end{align}
	Sublinear regret is only possible when $\varPhi_T$ is small.

	The next result shows that by implementing decreasing $\rho(t)$, the regret of Algorithm \ref{alg1} is bounded by $\varPhi_T$ and $\sum_{t=1}^T\sqrt{\rho(t)}$.
	\begin{theorem}\label{Theorem_1}
		Suppose that {Assumption} {\ref{assumption_weight}} and {\ref{assumption_domain_set}} hold, for $\forall i\in \mathcal{N}$, the regret of Algorithm \ref{alg1} is bounded by
		\begin{align}\label{theorem_1}
			R_i(T) \leq \mathcal{O}\left(\sqrt{T\left(\frac{\varPhi_T+1}{\rho^2(T)}+\sum\nolimits_{t=1}^T\sqrt{\rho(t)}\right)}\right)
		\end{align}
	\end{theorem}
	The proof of {Theorem \ref{Theorem_1}} is given in the appendix.
	\begin{remark}
		From {Theorem \ref{Theorem_1}}, we can get the sublinearly bounded dynamic regrets of Algorithm \ref{alg1}. To this end, take $\rho(t)=\frac{K}{(at+b)^\alpha}$ with $a,b>0$, $0<K\leq 1$ and $0<\alpha<\frac{1}{2}$. If $\varPhi_T$ is sublinear with $T^{1-\frac{\alpha}{2}}$, then $R_i(T)$ is sublinear with $T$. Note that the sublinearity of $\varPhi_T$ is a general assumption, which is widely adopted in online optimization and online game \cite{Kaihong2021Environments,2021arXiv210506200M,Kaihong2020Pseudoconvex}. In the P2P energy market, it mainly limits the fluctuation of the undispatchable load demand, which results in a GNE sequence with limited deviation.

	\end{remark}

	\section{Case Studies}

	In this section, numerical simulation on a case with six prosumers (microgrids) is introduced to verify the effectiveness of the proposed algorithm.

	\begin{figure}[t]
		\centering
		\includegraphics[width=0.18\textwidth]{topology}
		\caption{Communication and trading graph $\mathcal{G}$ for the case study.}
		\label{fig_topology}
	\end{figure}

	Each prosumer communicates and trades with its neighboring prosumers via a connected graph $\mathcal{G}$, which is shown in Fig. \ref{fig_topology}. The adjacency matrix $W$ is set to be $w_{i,j}=w_{j,i}=\frac{1}{{N}_i+1}$ and $\sum_{j\in \mathcal{N}}w_{i,j}=1$ for $\forall i\in \mathcal{N}$.
	The time interval is set to be $\Delta t=1\ \text{min}$ and $T=1440$, which indicates that the optimization period is 24 hours. The learning rate is set to be $0.8\sqrt[3]{\frac{5}{0.1t+5}}$.

	We assume that prosumer 1, 3, 4, and 6 are equipped with PVs and $p_i^l(t)$ of these prosumers is equal to the difference between their loads and PV generations, while $p_i^l(t)$ of prosumer 2 and 5 only have a load demand. The minute-sampled profile of PV generation is obtained from \cite{website_PV}, which is collected in Utah, from the U.S.. We use the data on 16th September 2013. Since each prosumer is located in the same area, we assume that the photovoltaic generation curves of different prosumers are the same, but with different amplitudes. The profiles of the daily power consumption of each prosumer are from \cite{website_load}. Fig. $\ref{fig_pv_load}$ depicts the load and PV generation of prosumer 1 from 6:00 am to 6:00 pm.

	\begin{figure}[t]
		\centering
		\includegraphics[width=0.3\textwidth]{P1_6to18.jpg}
		\caption{Load and PV generation of prosumer 1.}
		\label{fig_pv_load}
	\end{figure}

	\begin{figure}[t]
		\centering
		\includegraphics[width=0.45\textwidth]{regret_6.jpg}
		\caption{The trajectories of the average regrets of each prosumer. Note log-scale on y-axis for the inlay figure.}
		\label{fig_reg}
	\end{figure}

	Fig. \ref{fig_reg} shows the trajectories of the average regrets of each prosumer from 6:00 am to 6:00 pm. We set $t=0$ at 6:00 am here and therefore the length of the learning period is 720. As shown in Fig. $\ref{fig_reg}$, $R_i(t)/t$, $i\in\mathcal{N}$ approximately decays to zero after $t=120$ (i.e., after two hours). The downtrend of the logarithm of $R_i(t)/t$, $i\in\mathcal{N}$ also validates this property, which is consistent with the result in {Theorem \ref{Theorem_1}}.

	\section{Conclusion}
	In this paper, we propose a distributed online algorithm to track the GNE of the P2P energy trading game in a time-varying environment. We prove that by an appropriate choice of the decreasing learning rate, the regret of the proposed algorithm is bounded by a sublinear function of learning time. Simulation results with six prosumers verify the effectiveness of the proposed algorithm. In future work, we may focus on the effect of communication delay on the performance of online algorithms.


\section*{}}

\begin{document}
\normalem

\title{Dichotomy of heavy and light pairs of holes in the $t-J$ model}

\author{A. Bohrdt}
\email[Corresponding author email: ]{bohrdt@fas.harvard.edu}
\address{ITAMP, Harvard-Smithsonian Center for Astrophysics, Cambridge, MA 02138, USA}
\affiliation{Department of Physics, Harvard University, Cambridge, Massachusetts 02138, USA}

\author{E. Demler}
\affiliation{Institut f\"{u}r Theoretische Physik, ETH Zurich, 8093 Zurich, Switzerland}

\author{F. Grusdt}
\affiliation{Department of Physics and Arnold Sommerfeld Center for Theoretical Physics (ASC), Ludwig-Maximilians-Universit\"at M\"unchen, Theresienstr. 37, M\"unchen D-80333, Germany}
\affiliation{Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4, D-80799 M\"unchen, Germany}

\date{\today}

\begin{abstract}

 A key step in unraveling the mysteries of materials exhibiting unconventional superconductivity is to understand the underlying pairing mechanism. While it is widely agreed upon that the pairing glue in many of these systems originates from  antiferromagnetic spin correlations \cite{Scalapino1999,Moriya2003,OMahony2022}, a microscopic description of pairs of charge carriers remains lacking.
 Here we use state-of-the art numerical methods to probe the internal structure and dynamical properties of pairs of charge carriers in  quantum antiferromagnets in four-legged cylinders. Exploiting the full momentum resolution in our simulations, we are able to distinguish two qualitatively different types of bound states: a highly mobile, meta-stable pair, which has a dispersion proportional to the hole hopping $t$, and a heavy pair, which can only move due to spin exchange processes and turns into a flat band in the Ising limit of the model.
Understanding the pairing mechanism can on the one hand pave the way to boosting binding energies in related models \cite{Bohrdt2022_bilayer}, and on the other hand enable insights into the intricate competition of various phases of matter in strongly correlated electron systems \cite{Scalapino2012,Qin2019}.

\end{abstract}

\maketitle

\textbf{Introduction.--}
Following the discovery of high temperature superconductivity in the cuprates, understanding the mechanism by which pairs of charge carriers can form in a system with repulsive interactions has been a key question in the field.
Motivated by experimental results on the cuprate materials, a lot of theoretical and numerical work has focused on identifying the potential pairing symmetry \cite{Scalapino1995,Tsuei2000} as well as the binding energies in these microscopic models \cite{Chernyshev1998,Leung2002}.
Despite a vast research effort over several decades, the existence of a superconducting phase in the simplest model describing interacting electrons, the Fermi-Hubbard model, remains debated \cite{Qin2019}.
Competing orders, such as charge density waves and stripes, contribute to the difficulty in realizing as well as understanding superconductivity \cite{Scalapino2012}.
In order to unravel the competition between different orders, and thus the conditions for the existence of a superconducting phase, it is essential to gain a deeper understanding of the nature of pairs of charge carriers.
The existence of pairs close to half-filling does not imply that for a finite density of holes, the system necessarily realizes a $d$-wave paired state. Instead, a finite number of charge carriers can for example self-organize into a charge or pair density wave state \cite{Fradkin2015}. However, understanding whether and how pairs form in the two-hole problem is crucial to the subsequent understanding of self-organization of many holes.

Here we approach the question of the underlying binding mechanism from a new perspective: through novel spectroscopic tools, we search for bound states of charge carriers in a quantum antiferromagnet and directly probe their internal structure.
In particular, we numerically simulate rotational two-hole spectra, where different angular momenta can be imparted on the system, using time-dependent matrix product states. Crucially, these rotational spectra go beyond the standard pairing correlations through the momentum resolution they provide. The momentum dependence of the peaks in the spectral function enables direct insights into the effective mass of the pairs, which is an essential property for understanding their ability to condense at finite doping and temperature.

We study pairing between holes doped into the two-dimensional $t-J$ model, which corresponds to the enigmatic Fermi-Hubbard model to second order in $t/U$ (up to next-nearest neighbor hopping terms, where $U$ is the on-site interaction) and describes electrons in cuprates \cite{Zhang1988}:
    \begin{multline}
\H_{t-J} = - t ~ \hat{\mathcal{P}} \sum_{\ij} \sum_\sigma \l  \hat{c}^\dagger_{\vec{i},\sigma} \c_{\vec{j},\sigma} + \text{h.c.} \r \hat{\mathcal{P}} + \\
+ J \sum_{\ij} \hat{\vec{S}}_{\vec{i}} \cdot \hat{\vec{S}}_{\vec{j}} - \frac{J}{4} \sum_{\ij} \hat{n}_{\vec{i}} \hat{n}_{\vec{j}},
\label{eqtJModel}
\end{multline}
where $\hat{\mathcal{P}}$ projects to the subspace with maximum single occupancy per site; $\hat{\vec{S}}_{\vec{j}}$ and $\hat{n}_{\vec{j}}$ denote the on-site spin and density operators, respectively. In our numerical simulations, we consider a 40 site long, four-legged cylinder.

In order to probe a possible bound state of two charges, we consider an extension of conventional angle-resolved photoemission spectroscopy (ARPES). In particular, we excite the system by creating not one, but two charges while simultaneously imparting angular momentum on the system. The resulting spectra thus directly contain information about the existence of possible bound states, their ground state energy, as well as their dispersion relation.
In our numerical matrix product state calculations, we find well-defined peaks in the rotational spectral function for all angular momenta, for spin singlet as well as triplet pairs, and throughout an extended frequency range.

In order to gain a deeper understanding of the rotational two-hole spectra, we also consider the conceptually simpler $t-J_z$ model, where the $SU(2)$ invariant spin interactions are replaced by Ising type interactions. Without additional spin dynamics, a direct comparison of our numerical results to an effective theory describing pairs of charge carriers bound by strings is possible, yielding excellent agreement in terms of the existence as well as the dispersion of the various bound states we observe.
In particular, we discover a strongly dispersive low-energy peak, with a dispersion scaling with the hole hopping $t$, as well as completely flat bands at competitive energies. We attribute the flat bands to destructive interference of pairs with $d$-wave symmetry \cite{supp}.

Upon introducing spin dynamics, the flat bands develop into weakly dispersive bands, whereas the $t$-dependent feature remains largely unchanged. We thus discover two qualitatively different kinds of bound states: highly dispersive peaks, including a high energy feature with strong spectral weight in the $s$-wave spectra; and a weakly dispersive band, which has a high amount of spectral weight in the $d$-wave spectra. The dispersion of the latter is determined by the spin coupling $J$. The emergence of a slow time-scale set by $J$ is intuitive in the case of a single hole, which forms a spinon-holon bound state and can thus only move as fast as the spin excitation \cite{Bohrdt2020_dynamics}. In contrast, it is surprising to find a bound state of two holes, which in principle can individually move with a hopping amplitude $t$, but when bound into a pair have a dispersion given by the spin exchange.

The remainder of this paper is organized as follows. We start by introducing the rotational two-hole spectra. We then discuss results for the $t-J_z$ model, where the SU(2) invariant spin interactions are replaced by Ising-type interactions. We discuss the features found in the numerically obtained spectra in detail and compare to a semi-analytical theoretical description of pairs of charge carriers \cite{Grusdt2022_arxiv}. Finally, we consider the full $t-J$ model.
\\

\begin{figure}
\centering
\includegraphics[width=0.99\linewidth]{spectrum_tJ_m0.pdf}
\caption{
\textbf{Rotational spectroscopy of two holes} in a singlet state in the $t-J$ model with $t/J=3$, on a $40\times 4$ cylinder, based on a time evolution up to $T_{max}/J=3$ and bond dimension $\chi=1200$. (a) Sketch of the response probed by the rotational spectrum. (b) The upper, middle, and lower plots correspond to $k_y=0,\pi/2,\pi$ and data is shown as a function of momentum $k_x$ and frequency $\omega/J$ for $m_4=0$; $m_4=2$ exhibits similar features, see \cite{supp}. Gray dashed lines correspond to a cosine dispersion $-2 J \alpha \cos (k_x) + b_J$, black line corresponds to a cosine dispersion $-2 t \alpha \cos(k_x) + b_t $, where $\alpha = 0.33 $ in both cases, $b_J=-11$, and $b_t= -9$.
}
\label{fig:fig1}
\end{figure}

\textbf{Rotational Spectra.--}
In order to probe the internal structure of pairs of charge carriers, we study rotational spectra.
We define an operator $\hat{\Delta}_{m_4}(\vec{j},\sigma,\sigma')$ that creates a pair of holes on the bonds adjacent to site $\vec{j}$ with discrete $C_4$ angular momentum $m_4=0,1,2,3$ as
\begin{equation}
    \hat{\Delta}_{m_4}(\vec{j},\sigma,\sigma') = \sum_{\vec{i}: \ij} e^{i \frac{\pi}{2} m_4 \varphi_{\vec{i} - \vec{j}}}    \c_{\vec{i},\sigma'} \c_{\vec{j},\sigma},
\end{equation}
with $\varphi_{\vec{r}} = {\rm arg}(\vec{r})$ the polar angle of $\vec{r}$.
In order to annihilate a spin-singlet, we define the singlet pair operator (and similar for triplets) as
\begin{equation}
    \hat{\Delta}_{m_4}^{(s)}(\vec{j}) =  \hat{\Delta}_{m_4}(\vec{j},\uparrow,\downarrow) - \hat{\Delta}_{m_4}(\vec{j},\downarrow,\uparrow).
    \label{eq:singlet}
\end{equation}

The simplest term creating a spin-singlet excitation with discrete angular momentum $m_4$, charge two, and total momentum $\vec{k}$ is directly given by the spatial Fourier transform of the singlet pair operator as
\begin{equation}
\hat{\Delta}_{m_4}^{(s)}(\vec{k}) =  \sum_{\vec{j}} \frac{e^{- i \vec{k} \cdot \vec{j}}}{\sqrt{V}}\hat{\Delta}_{m_4}^{(s)}(\vec{j})
\label{eqDefRmkSinglet}
\end{equation}
with volume $V$. The discrete angular momentum $m_4$ is a good quantum number at $C_4$ invariant momenta $\vec{k} = (0,0), (\pi,\pi)$ only.
Based on this operator, we now consider the rotational Green's function
\begin{equation}
\mathcal{G}_{\rm rot}^{(m_4)}(\vec{k},t) = \theta(t)  \bra{\Psi_0} \hat{\Delta}_{m_4}^{(s) \dagger}(\vec{k},t) \hat{\Delta}_{m_4}^{(s)}(\vec{k},0) \ket{\Psi_0},
\label{eq:rotationalGF}
\end{equation}
which we calculate using time-dependent matrix product states \cite{Kjall2013,Zaletel2015,Paeckel2019}. The corresponding two-hole rotational spectrum, $- \pi^{-1} {\rm Im} \mathcal{G}_{\rm rot}^{(m_4)}(\vec{k},\omega)$, in Lehmann representation is
\begin{equation}
 A_{\rm rot}^{(m_4)}(\vec{k},\omega) =  \sum_{n} \delta \l \omega-E_n+E_0^0 \r  | \bra{\Psi_n} \hat{\Delta}_{m_4}^{(s)}(\vec{k})  \ket{\Psi_0^0} |^2,
 \label{eqDefArot}
\end{equation}
where $\ket{\Psi_0^0}$ ($E_0^0$) is the ground state (energy) of the un-doped system and $\ket{\Psi_n}$ ($E_n$)  are the eigenstates (eigenenergies) with two holes.

The two-hole rotational spectral function defined above is closely related to the dynamical pairing correlations frequently considered in the literature \cite{Dagotto1990pair,Poilblanc1994,Qin2019},
\begin{equation}
    P(\omega) = \int dt ~e^{i\omega t} \left< \hat{\Delta}_{m_4}^{(s) \dagger}(t) \hat{\Delta}_{m_4}^{(s)}(0) \right>,
\end{equation}
where $    \hat{\Delta}_{m_4}^{(s)} = \sum_{\vec{j}} \hat{\Delta}_{m_4}^{(s)}(\vec{j})$. We, however, consider the full momentum dependence of the pairing correlations, which enables direct insights into the center-of-mass dispersion of pairs of charge carriers.

The resulting rotational spectra thus directly probe the existence of bound states and their internal structure: If a bound state of two holes with long-lived rotational excitations exists, the rotational spectra should exhibit well-defined coherent peaks. If on the other hand such bound states do not exist, the excitation with the rotational operator $\hat{\Delta}_{m_4}(\vec{k})$ will lead to a broad continuum in the corresponding spectral function.

In Fig.~\ref{fig:fig1}, we show the two-hole spectral function without additional angular momentum, i.e. $m_4=0$ ($s$-wave), for the $t-J$ model for momenta $k_y=0,\pi/2$, and $\pi$ and $0\leq k_x \leq \pi$. We find a well-defined coherent peak at low energies for all momenta, indicating the existence of a bound state. The spectrum furthermore reveals a plethora of different features, including a highly dispersive band (black line) as well as bands with a dispersion proportional to the spin-exchange $J$ (gray dashed lines). At momentum $\mathbf{k}=(\pi,\pi)$, the spectral weight vanishes for all energies since $\hat{\Delta}_{0}^{(s)}(\vec{k}=(\pi,\pi))=0$.

In order to gain a deeper understanding of these intriguing results, we take a step back and analyze the conceptually simpler $t-J_z$ model in the following section.
\\

\begin{figure}
\centering
\includegraphics[width=0.99\linewidth]{spectrum_tJz_m0.pdf}
\caption{
\textbf{Two hole rotational spectra in the $t-\rm{XXZ}$ model}
for $t/J_z = 3$ and $J_\perp/J_z=0.1$ and $m_4=0$ ($s$-wave) on a $40 \times 4$ cylinder, based on time evolution up to $T_{max}/J_z = 10$ and bond dimension $\chi=600$. The colormap corresponds to numerical matrix product state simulations of the singlet two-hole rotational spectrum, blue lines are geometric string theory predictions for the position of states (all shifted by $-0.35J$), and the black line is a cosine fit. The upper, middle, and lower plots correspond to $k_y=0,\pi/2,\pi$ and data is shown as a function of momentum $k_x$ and frequency $\omega/J_z$.
The overall ground state energy for two holes is marked by orange circles for $k_y=0$, and the green dashed line in the top panel corresponds to twice the energy of a single hole.
}
\label{fig:tJz}
\end{figure}
\textbf{The $t-\rm{XXZ}$ model.--}
We now consider a modification of the $t-J$ model, where the SU(2) invariant spin interactions are replaced by in-plane and Ising-type spin interactions with coupling constants $J_\perp$ and $J_z$:
\begin{multline}
\H_{t-\rm{XXZ}} = \sum_{\ij} \left( J_\perp \left(\hat{S}^x_{\vec{i}} \hat{S}^x_{\vec{j}}+\hat{S}^y_{\vec{i}} \hat{S}^y_{\vec{j}}\right)  +J_z\hat{S}^z_{\vec{i}} \hat{S}^z_{\vec{j}} \right) -
\\
- t ~ \hat{\mathcal{P}} \sum_{\ij} \sum_\sigma \l  \hat{c}^\dagger_{\vec{i},\sigma} \c_{\vec{j},\sigma} + \text{h.c.} \r  \hat{\mathcal{P}} - \frac{J_z}{4} \sum_{\ij} \hat{n}_{\vec{i}} \hat{n}_{\vec{j}}  .
\label{eqtJzModel}
\end{multline}
In the limit of $J_\perp \ll J_z$, also called the $t-J_z$ model, the lack of spin dynamics facilitates our theoretical understanding. Experimentally, the anisotropic interactions can for example be realized by employing Rydberg interactions \cite{GuardadoSanchez2020Ry}.

\begin{figure}
\centering
\includegraphics[width=0.99\linewidth]{peakPos_full.pdf}
\caption{
\textbf{Strongly dispersive pair state} in the $t-\rm{XXZ}$ model
for $t/J_z=1,2,3$ and $J_\perp/J_z=0.1,1.0$ and $m_4=0$ ($s$-wave). The symbols correspond to the position of the lowest energy peak at $k_y=\pi$ extracted from numerical matrix product state simulations of the singlet two-hole rotational spectrum. Yellow crosses correspond to the isotropic case, $J_z=J_\perp=J$ with $t/J=3$. All data points are shifted vertically to collapse at $k_x=0$.
Blue dashed lines are geometric string theory predictions for the position of the lowest states. For $t/J_z=2,3$, the lowest two states predicted by the geometric string theory are shown.
Black line is a cosine fit, $0.62 \cos(k_x) + 0.72$, to the extracted peak positions for $t/J_z=3$, $J_\perp/J_z=0.1$.
}
\label{fig:peak_pos}
\end{figure}

Remarkably, the two-hole spectral function, Fig.~\ref{fig:tJz}, exhibits a highly dispersive peak with a mass proportional to $1/t$; I.e., we find a long-lived, tightly bound state of two holes, which can move as fast as the hole hopping $t$. This is in stark contrast to the case of a single hole in the same model, which has a very high effective mass $\gg 1/t$ and thus an almost flat dispersion \cite{Grusdt2018tJz}, since it can only move due to Trugman loops \cite{Trugman1988}, which are higher order processes.

In Fig.~\ref{fig:peak_pos}, we further analyze the scaling of the mass of the pair by analyzing the position $\omega_{\rm max}$ of the lowest energy peak at $k_y=\pi$ as a function of $k_x$ for different values of $t/J_z=1,2,3$ and $J_\perp/J_z=0.1$. Note that here, the frequency $\omega_{max}$ is shown in units of the hole hopping $t$. We added a constant to the curves for different values of $t/J_z$ to show  the remarkable agreement in the overall shape of the dispersive peak for different values of $t/J_z$, thus highlighting the scaling with the hole hopping $t$.

The lowest lying peak for $k_y=0$ and $k_y=\pi/2$ is -- within our numerical resolution -- completely flat. Note that the situation of two unbound, and thus approximately immobile, holes should have a very small matrix element in the spectral function considered here, and therefore cannot account for the pronounced flat band peaks we find in the two-hole spectra. This is further corroborated by the direct comparison of the energy for two holes ($E_{2h}-E_{0h} = -12.08 J_z$, orange circles in top panel Fig.~\ref{fig:tJz}) and twice the energy of a single hole ($2 \cdot (E_{1h}-E_{0h}) = -11.01 J_z$, green dashed line in top panel Fig.~\ref{fig:tJz}).

In a companion paper \cite{Grusdt2022_arxiv}, we extend the geometric string theory developed for a single hole \cite{Grusdt2018tJz,Grusdt2019} to the case of two charge carriers. In particular, this geometric string theory approach describes the properties of two holes bound together by a string of displaced spins, and thus provides estimates of the energy and dispersion relation of such pair states, see blue lines in Figs.~\ref{fig:tJz} and \ref{fig:peak_pos}. Note that the existence of a state at a given energy does not imply that said state is visible in the spectral function, since the spectral weight, i.e. the overlap with the excitation we consider, can still be zero.

The geometric string theory correctly predicts the highly dispersive peak, as well as the existence of completely flat bands. The latter have also been predicted by a related theoretical study \cite{Shraiman1988}.
\\

\textbf{Results for the $t-J$ model.--}
In the isotropic case, $J_\perp = J_z$, the strongly dispersive peak remains visible, see Fig.~\ref{fig:fig1} bottom. In Fig.~\ref{fig:peak_pos}, we compare the peak position for the isotropic $t-J$ model with the $t-\rm{XXZ}$ model with $J_\perp/J_z=0.1$ and find that the momentum dependence along the $x$-direction of the lowest lying peak for $k_y=\pi$ is qualitatively very similar between the two cases. This indicates in particular that also in the $t-J$ model, a highly mobile, tightly bound pair state exists.

\begin{figure}
\centering
\includegraphics[width=0.99\linewidth]{peakPos_JperpDep.pdf}
\caption{
\textbf{Weakly dispersive pair} in the $t-\rm{XXZ}$ model
for $t/J_z=3$, $J_\perp/J_z=0.1,0.5,1.0$, and $m_4=2$ ($d$-wave). The symbols corresponds to the position of the lowest energy peak at $k_y=0$ extracted from matrix product state simulations of the singlet two-hole rotational spectrum.
Dashed lines are a cosine fit to the extracted peak positions with prefactor as indicated in the legend, and an additional offset (not indicated).
}
\label{fig:peak_pos2}
\end{figure}

The flat bands, particularly visible in the $J_\perp/J_z=0.1$ case at momenta $k_y = 0,\pi/2$, acquire a dispersion approximately proportional to $J_\perp^2/J_z$, as can be seen in the top and middle plot in Fig.~\ref{fig:fig1}. We analyze this behavior in more detail by explicitly extracting the peak position $\omega_{\rm max}$ at $k_y=0$ for $t/J_z=3$ and different values of $J_\perp/J_z$ in Fig.~\ref{fig:peak_pos2}.

Again, the geometric string theory \cite{Grusdt2022_arxiv} correctly captures the highly dispersive peak with mass proportional to $1/t$ (black lines) discussed in Fig.~\ref{fig:peak_pos}. However, since this theoretical description does not account for spin dynamics, it does not predict the dispersion proportional to $J_\perp^2/J_z$ (gray dashed lines). The corresponding bands in the $t-J_z$ model are flat, as predicted by the geometric string theory.

\begin{figure}
\centering
\includegraphics[width=0.99\linewidth]{spectrum_ky0_m0vs2.pdf}
\caption{
\textbf{Angular momentum dependence} of the rotational two-hole spectra in the $t-J$ model with $t/J=3$ for $m_4=0,2$ (left, right column) and $k_y=0$, obtained from time-dependent matrix product state simulations of the singlet two-hole rotational spectrum.
}
\label{fig:m_dep_tJ}
\end{figure}

\textbf{Different angular momenta.--}
Thus far, we considered the case of $s$-wave pairs, i.e. $m_4=0$ and found a well-defined coherent peak at low energies.
In Fig.~\ref{fig:m_dep_tJ}, we show the spectral function for the $t-J$ model at $t/J=3$ at momentum $k_y=0$ for $m_4=0,2$, corresponding to $s$- and $d$-wave excitations.
Note, however, that the $40 \times 4$ cylinder geometry used in our numerical simulations breaks the $C_4$ symmetry, and therefore $m_4=0$ and $m_4=2$ excitations can hybridize.

The spectral weight exhibits a strong dependence on the excitation, see Fig.~\ref{fig:m_dep_tJ}. Before, we have identified two qualitatively different types of bound states: a strongly dispersive band, see Fig.~\ref{fig:peak_pos}, and a weakly dispersive band, see Fig.~\ref{fig:peak_pos2}. By considering different rotational excitations, we find that the weakly dispersive pair, realizing the overall ground state, has most spectral weight for the $d$-wave ($m_4=2$) excitation, whereas for the $s$-wave ($m_4=0$) excitation, a large fraction of the spectral weight appears in a strongly dispersive high energy feature. This high-energy feature is also predicted by the geometric string theory \cite{Grusdt2022_arxiv}.

The hybridization of different $m_4$ excitations in the cylinder geometry leads to the appearance of the weakly dispersive feature, found with strong spectral weight for $m_4=2$, also at $m_4=0$, albeit with a smaller spectral weight.
The overall lowest energy peak appears at momentum $\mathbf{k}=(0,0)$ and $m_4=0,2$.

We note that in early exact diagonalization studies on small systems \cite{Dagotto1990pair,Poilblanc1994}, the integration over all momenta, or the consideration of only momentum zero, leads to a much sharper low energy peak in the $d$-wave than in the $s$-wave case. We attribute this to the widely different dispersions of the $s$- and $d$-wave peaks revealed here, and the exact $C_4$ symmetry in small systems.

The geometric string theory correctly predicts the $d$-wave character of the weakly dispersive band, as well as the accumulation of spectral weight at $\vec{k}=(0,0)$ \cite{Grusdt2022_arxiv}. A more detailed comparison of spectral weights for different values of $m_4$ shows good qualitative agreement at low energies (see supplementary material).

So far, we considered a singlet excitation. In order to annihilate a triplet instead, we define the triplet pair operator as
\begin{equation}
    \hat{\Delta}_{m_4}^{(t)}(\vec{j}) =  \hat{\Delta}_{m_4}(\vec{j},\uparrow,\downarrow) + \hat{\Delta}_{m_4}(\vec{j},\downarrow,\uparrow).
    \label{eq:triplet}
\end{equation}
Upon considering the corresponding triplet spectral function, we  find that the lowest lying peaks are at higher energies than in the case of the singlet spectral function. This finding furthermore suggests that the lowest energy peak in the singlet spectral function cannot be attributed to unbound states of two holes.
\\

\textbf{Summary and Outlook.--}
In this work, we extensively studied the properties of pairs of charge carriers in the $t-J$ and $t-\rm{XXZ}$ model through rotational spectra. We find well defined coherent peaks at low energies for all momenta and angular momenta $m_4$. Our work provides the first extensive numerical study of the mass of pairs of charge carriers in extended systems. We have revealed two qualitatively different types of bound states: First, a weakly dispersive peak with a dispersion approximately proportional to $J_\perp^2/J_z$, which has most spectral weight for a $d$-wave excitation. And second, a highly dispersive peak, which corresponds to tightly bound pairs with a mass proportional to $1/t$.
We find the same signatures of these light pairs of charge carriers in the $t-J_z$ and the SU(2) invariant $t-J$ models. The band corresponding to these light pairs, as well as the bands corresponding to heavy pairs, are qualitatively captured by a semi-analytic geometric string theory approach \cite{Grusdt2022_arxiv}.

An intriguing direction for future research is the direct experimental probe of the two-hole rotational states. Understanding the pairing mechanism in the Fermi-Hubbard and related models has been one of the key motivations in the development of quantum simulators, and in particular cold atoms in optical lattices \cite{Hofstetter2002}. In the past two decades, remarkable progress has been made in the field \cite{Tarruell2018,Bohrdt2021_review}, and several proposals to probe the pairing symmetry have been put forward \cite{Rey2009,Kitagawa2011}. More recently, the single-hole spectral function has been measured experimentally \cite{Brown2019,Bohrdt2018}. Using additional lattice modulations, the two-hole rotational spectral function considered here could be accessed experimentally. In solid state experiments, the $s$-wave two-hole spectral function can be accessed through coincidence angle-resolved photo-emission spectroscopy \cite{Su2020}, which relies on simultaneous measurements of two photo-electrons and provides direct insights into the pair Green’s function.
A different approach is based on Anderson-Goldman pair tunneling in a tunnel junction setup \cite{Anderson1967a,Anderson1972,Scalapino1970}: to study the structure of individual pairs in a strongly underdoped quasi-2D material as considered here, we propose to tunnel-couple the latter to a probe-superconductor along $z$-direction. Momentum resolution can in this case be obtained through an in-plane magnetic field.

The observation of light as well as heavy pairs in the spectra shown here furthermore suggests a real-space and -time experiment.
Upon slowly releasing two holes next to each other in a cold atom experiment, a low energy state of the pair can be prepared. In the ensuing time evolution, we predict the pairs to spread through the system in two distinct wave-fronts, corresponding to the light and the heavy pair, respectively.

Finally, the existence of flat or weakly dispersive bands opens a new avenue to understand the many competing orders found experimentally in cuprate materials as well as numerically in Fermi-Hubbard and $t-J$ models at finite doping \cite{Corboz2014,Qin2019,Arovas2022}. In a next step, we will investigate the dichotomy between two types of light and heavy pairs in the Fermi-Hubbard model.

\textbf{Acknowledgements.--}
We thank Immanuel Bloch, Antoine Georges, Markus Greiner, Mohammad Hafezi, Lukas Homeier, and Ulrich Schollw\"ock for fruitful discussions. This research was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -- EXC-2111 -- 390814868, by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programm (Grant Agreement no 948141) — ERC Starting Grant SimUcQuam, by the ARO (grant number W911NF-20-1-0163), and by the NSF through a grant for the Institute for Theoretical Atomic, Molecular, and Optical Physics at Harvard University and the Smithsonian Astrophysical Observatory.

\section*{References}