\section{Introduction}
	As a long-lasting problem, phase optimization has been studied in recent decades for distribution networks~\cite{TPUVO, OPBDS, TPOPA, feron_2021}.
	A more balanced power supply can improve power delivery efficiency, reduce power line hazards during transmission, and increase decision flexibility for a variety of ancillary services. A balanced three-phase power flow can be produced by dynamically controlling the power supply for load units from each phase line. Such a control mechanism is inherently compatible with an EV charging network. Indeed, phase optimization for EV charging can produce substantial economic benefits, especially in a large-scale and high-density (LSHD), power-delivery adaptive, and information-sensitive network. A three-phase adaptive charging network (ACN)~\cite{ACN2018} includes all those features that strongly motivate the necessity of phase optimization. A real-world ACN often has a three-phase configuration, equipped with software-defined charging control. Our work relies on the Caltech ACN, a multi-story charging facility at Caltech campus~\cite{li2019acndata,lee2021acnsim}.

	\begin{figure}
	\centering
		\includegraphics[width=0.4\textwidth]{python_images/parking_system.png}
		\caption{The ACN phase optimization process. \textit{Step 1.} EV $i$ sends information \textit{a priori}: e.g., arrival time $a_i$, parking duration $d_i$, and energy demand $e_i$. \textit{Step 2.} the controller computes the optimal phase selection $X_i\in\{ab,bc,ca\}$ and charging strategy $A_i$. \textit{Step 3.} the controller returns a set $S_i$ of all empty slots whose corresponding phases match $X_i$. \textit{Step 4.} EV $i$ chooses a slot from $S_i$ and starts charging according to $A_i$.}
		\label{fig: PC}
	\end{figure}

    In theory, ACN phase optimization is a problem that can be solved by either the demand or the supply side. Existing charging infrastructures, such as~\cite{ACN2018},  assume fixed phase selections determined by demand providers such as EV drivers, and thereby a globally optimal solution can be exactly computed by the supply side. In contrast to this, when the supplier of charging facilities is able to control the phase selections, i.e., the energy demand allocated to phases becomes a decision variable, the charging system can be operated in a more energy and cost-efficient way by maintaining a balanced or near-balanced network. The control of phase selections can be achieved by either deploying PV inverters with a balancing inverter~\cite{LBEVC} or developing smart mobile applications to inform a driver of his desirable parking/charging locations. In fact, balancing phases via either PV inverters or providing phase suggestions to drivers require to solve a phase optimization problem. In this work,  we will formally give a mathematical formulation for the phase optimization problem and provide algorithms for proactive phase allocation.

	In practice, we consider an ACN phase optimization process, as illustrated in Figure~\ref{fig: PC}. Based on charging inputs provided by an EV driver upon arrival, an ACN controller computes optimal phase(s) within seconds and shares the suggested location(s) to the driver through a mobile app or payment machines at the entrance of the charging station. Those available locations have chargers mapping to one of the phases, i.e., $\{ab, bc, ca\}$ to which the vehicle will connect. It is worth emphasizing that such a phase-optimization and the demand-allocation process can significantly improve charging efficiency and reduce electricity cost if the following assumptions are satisfied: (1). Each phase is connected with a sufficient amount of ACN chargers. (2). Charger-to-phase connection is immutable. (3) EV drivers' charging inputs are accurate, i.e., the charging duration and energy demand. (4) EV drivers follow suggested charging/parking instructions. Note that (1) and (2) are reasonable assumptions for an LSHD ACN; (3) and (4) are practical assumptions by penalizing the violation of ACN instructions, e.g., by disabling the charging service at places corresponding to the phases that are not allocated to a driver.

	\textit{Related Work.}
    There is vast literature on the EV charging problem for ACNs. The authors of \cite{li2019acndata} proposed model predictive control-based methods for EV charging. Optimal EV charging as second-order-cone programming (SOCP) was proposed and solved in~\cite{7905971} and \cite{8274166}. The authors of~\cite{lee2020adaptive} illustrated detailed ACN architectures and proposed a novel model-predictive-control (MPC) algorithm that handles unbalanced three-phase infrastructure. With historical ACN charging records, reinforcement learning was used to generate penalty-based terminal functions that facilitate energy-cost efficiency~\cite{li2021learning}.
    Despite that many ACN charging algorithms have been studied and implemented in practice, balancing three-phase EV charging demands for ACNs has not yet attracted enough attention, except for a few works that tackled imbalanced EV connections. For instance,~\cite{EMEV, BCDS} used optimization-based methods to balance renewable penetrations, e.g., to minimize the peak-valley gap. Similarly,~\cite{PVUI, Rezai2016GridablePE} handled imbalanced phase connections by analyzing voltage sensitivity and demand response. However, none of them have studied the phase optimization problem by finding an optimal dispatch of committed EVs for real-world ACNs.

	\textit{Main Contributions.}
    The main contributions of this paper are three-fold. First, to find optimal phase allocations and charging strategies that minimize a given cost function, we formulate a mathematical model for the ACN phase optimization problem in Section~\ref{pf}. To the best of our knowledge, our work is among the first that considers such an ACN phase optimization problem. Second, we propose an algorithm that finds a near-optimal solution by solving a bi-level optimization and show that it is optimal in the zero laxity case. Moreover, we design an online approach for optimizing and allocating phases in a real-world ACN and validate the efficacy of our algorithms using realistic ACN data.

The rest of this paper is organized as follows. In \autoref{pf}, we describe the model and formulate our phase optimization problem. Next, we demonstrate the intractability of this problem in \autoref{ppop} and propose our method in \autoref{opscad}. Finally, we provide our experimental results in \autoref{sec:exp}.

	\section{Problem Formulation} \label{pf}
	\subsection*{Notation and conventions}
	\makebox[3cm]{$A_{(i)}, A^{(j)}$}  $i^{th}$ row, $j^{th}$ column of $A$ \par
	\makebox[3cm]{$A \leq B$}  $A_{(i,j)} \leq B_{(i,j)}, \forall i,j$ \par
	\makebox[3cm]{ $\mathbb{R}_{+}, \mathbb{Z}_{+}$}  Non-negative real, integer fields \par

	\subsection{EV charging model} \label{evcs}

	Suppose there is a set of $N$ electric vehicles (EVs), denoted by $\mathcal{N} \vcentcolon= \{1, \cdots, N\}$ that arrive at an adaptive charging network (ACN) during a discrete time horizon $\mathcal{T} \vcentcolon= \{1, \cdots, T\}$.
	Each EV $i$ has a charging-profile $\mathbf{ev}_i \in \mathbb{R}_{+}^3$ defined as a triple $[a_i, d_i, e_i]$ where  $a_i$ (hrs) denotes the EV's arrival time; $d_i$  (hrs) is the charging duration and $e_i$ (kWh) is its initial energy demand. We define the matrix of all $\mathbf{ev}$ as $D \in \mathbb{R}_{+}^{N \times 3}$, where $D_{(i)} \vcentcolon= \mathbf{ev}_i$.
Based on $D$, an $N \times T$ binary matrix $E$ with
	\begin{equation}
		\label{eqn:E}
		E_{(i, j)} =
		\begin{cases}
			1      & \quad \text{if } a_i \leq j \leq a_i + d_i\\
			0  & \quad \text{otherwise}
		\end{cases}
	\end{equation} has its
entry $E_{(i,j)}$ indicates whether the EV $i$ is presented at the ACN at time $j$.
	Moreover, we define an \textit{EV charging matrix} $A \in \mathbb{R}_{+}^{N \times T}$, whose $(i, j)$-th entry is the power (kW) charged for the EV $i$ at time $j$. Finally, we use a $3 \times N$ binary \textit{phase selection matrix} $X$ to denote the assignment of phases for EVs. The set consisting of all such $X$ is
	\begin{equation}
		\mathcal{X} \vcentcolon= \Big\{X \in \{0,1\}^{3 \times N} | \displaystyle \sum ^{3}_{i = 1} X_{(i, j)} = 1, \forall i \in \mathcal{N}\Big\}.
	\end{equation}

In the sequel, with the matrices $D,E,A$ and $X$ and the set $\mathcal{X}$ defined above, we describe the adaptive charging network (ACN) model used in this paper.

	\begin{figure}
		\centering
		\includegraphics[width=0.45\textwidth]{python_images/acn_circuit.png}
		\caption{The ACN circuit model involving infrastructure constraints of the Delta-to-Wye step-down transformer and ACN chargers.}
		\label{fig: caltechAcn}
	\end{figure}

	\subsection{The ACN model} \label{pfacnf}

	In the following, we describe the constraints of the ACN circuit model depicted in Figure~\ref{fig: caltechAcn}, which has been widely used to model real-world charging facilities such as the California Parking Garage~\cite{7905971,lee2020adaptive}. All charging variables and network constraints are in the unit of power (kW). At each time step $t$, we have the phase power vector  shown in the figure $[P^0_{ab}, P^0_{bc},P^0_{ca}] = (XA)^{(t)}$. Our model considers the charging constraints from the EVSE capacity and EV sessions, along with the network constraint induced by the transformer's capacity. They are the main constraints of the real-world Caltech ACN~\cite{ACN2018}.

	\subsubsection{Charging constraints}
First, the charging matrix $A$ must satisfy the charger constraint $\mathcal{C}_{\mathrm{r}}$ such that
	\begin{equation}
		\label{eqn:C_r}
		\mathcal{C}_{\mathrm{r}} \vcentcolon= \big\{A \in \mathbb{R}_{+}^{N \times T} | \mathbf{0}_{N \times M} \leq A \leq r_{\max} E\big\}.
	\end{equation}
We assume all chargers have the same maximum power rate $r_{\max} \in \mathbb{R}_{+}$ without loss of generality.
Moreover, $A$ needs to satisfy the demand constraints $\mathcal{C}_{\mathrm{d}}$ such that
	\begin{equation}
		\label{eqn:C_d}
		\mathcal{C}_{\mathrm{d}} \vcentcolon= \Big\{A \in \mathbb{R}_{+}^{N \times T} | \Big(\displaystyle\sum_{j = 1}^{T} A^{(j)}\Big)\Delta_T \leq D^{(3)}\Big\}.
	\end{equation}

\subsubsection{Network constraints}\label{nc}

The step-down transformer capacity induces line power limits on the primary side $c_2$ and the secondary side $c_1$ in Figure~\ref{fig: caltechAcn}. We denote the set of all $(X,A)\in \mathcal{X} \times \mathbb{R}_{+}^{N \times T}$ satisfying this network constraint by
	\begin{equation}
		\label{eqn:C_soc}
		\begin{aligned}
			\mathcal{C}_{\mathrm{soc}} & \vcentcolon=
			& \left\{\begin{array}{c | c}
				(X, A)  &
				\left|\left[
				\begin{array}{c}
					\Phi_1 \\
					\hline
					\Phi_2
				\end{array}
				\right] X A\right| \leq C_{\max}
			\end{array}\right\}.
		\end{aligned}
	\end{equation}
In~\eqref{eqn:C_soc}, $|\cdot|$ is the element-wise magnitude operator on the complex matrix, $C_{\max} \in \mathbb{R}^{6 \times T}$ is the line-power limit matrix, and the fraction symbol between $\Phi_1$ and $\Phi_2$ denotes a vertical concatenation of block matrices. The elements of each row depend on the transformer's capacity and e thline power supply. Without loss of generality, mimicking the model used in~\cite{ACN2018,lee2020adaptive}, we consider the effects of the transformer's capacity, i.e. for any $j \in \mathcal{T}$, $C_{\max(i,j)} = c_1$ if $i \in \{1, 2, 3\}$ and $c_2$ if $i \in \{4, 5, 6\}$. The matrix $\Phi_1$ transforms $P_i^0$ to $P_j^1$ and $\Phi_2$ transforms $P_i^0$ to $P_j^2$ ($i \in \{ab, bc, ca\}, j \in \{a, b, c\}$). Given a turning-ratio $n_r$, $\Phi_1$ and $\Phi_2$ can be explicitly defined as
	\begin{equation}
	\label{eqn:Phi1}
	\Phi_1 \vcentcolon=
	\begin{bmatrix}
		1\angle 30 \degree & 0\angle -90 \degree & -1\angle 150 \degree \\
		-1\angle 30 \degree & 1\angle -90 \degree & 0\angle 150 \degree \\
		0\angle 30 \degree & -1\angle -90 \degree & 1\angle 150 \degree
	\end{bmatrix},
    \end{equation}
    \begin{equation}
	\label{eqn:Phi2}
	\Phi_2 \vcentcolon=
	\begin{bmatrix}
		\frac{1}{n_r}\angle 30 \degree & \frac{1}{n_r}\angle -90 \degree & -\frac{2}{n_r}\angle 150 \degree \\
		-\frac{2}{n_r}\angle 30 \degree & \frac{1}{n_r}\angle -90 \degree & \frac{1}{n_r}\angle 150 \degree \\
		\frac{1}{n_r}\angle 30 \degree & -\frac{2}{n_r}\angle -90 \degree & \frac{1}{n_r}\angle 150 \degree
	\end{bmatrix}.
    \end{equation}

For notational simplicity, the parameters in our ACN model are summarized as a $4$-tuple: $\mathsf{ACN} \vcentcolon= \{D, \mathcal{C}_{\mathrm{r}}, \mathcal{C}_{\mathrm{d}}, \mathcal{C}_{\mathrm{soc}}\}$.

	\subsection{Phase optimization for optimal charging}\label{pfocf}
	\label{sec:pop}

	Subject to the constraints described above, the goal of this work is
	to maximize the aggregated energy delivery in $\mathcal{T}$. We define a class of linearly parameterized objective functions $f_{\mathbf{w}}: \mathbb{R}^{N \times T} \to \mathbb{R}$ of the EV charging matrix $A$, where
$
	 	f_{\mathbf{w}}(A) \vcentcolon= -\sum^{N}_{i = 1} (A\mathbf{w})_i
$
	 and $\mathbf{w} \in \mathbb{R}_+^T$ is  non-negative.
\begin{remark}
The objective   $f_{\mathbf{w}}$ captures some important properties in the EV charging problem and covers various engineering objectives:
	\begin{enumerate}
		\item if $w_i = 1, \forall i\in \mathcal{T}$, $f_{\mathbf{w}}$ is the power aggregation function. For convenience, we denote such an $f_{\mathbf{w}}$ as $f_{\mathbf{1}}$;
		\item if $w_i = T - i + 1, \forall i \in \mathcal{T}$, $f_{\mathbf{w}}$ is the quick-charge function defined in Eq. 2 of \cite{ACN2018};
		\item if $w_i = p_i \geq 0, \forall i \in \mathcal{T}$, $f_{\mathbf{w}} \Delta_T$ can be used to represent the total energy price paid for the electricity suppliers.
	\end{enumerate}
\end{remark}

The phase optimization problem considered in this paper is presented below. For the ease of our presentation, we fix $\mathbf{w}=\mathbf{1}$ in the remainder of this paper, and our results can be straightforwardly extended to the general case.
	\begin{mini!}|s|[2]
		{\begin{aligned}
				& X \in \mathbb{Z}_{+}^{3 \times N}\\
				& A \in \mathbb{R}_{+}^{N \times T}\\
		\end{aligned}}
		{f_{\mathbf{1}}(A) \label{P1:0}}
		{\label{P:1}}
		{}
		\addConstraint{X}{\in \mathcal{X}, \label{P1:1}}
		\addConstraint{A}{\in \mathcal{C}_{\mathrm{r}} \cap \mathcal{C}_{\mathrm{d}}, \label{P1:2}}
		\addConstraint{(X, A)}{\in \mathcal{C}_{\mathrm{soc}}. \label{P1:3}}
	\end{mini!}
where~\eqref{P1:2} and~\eqref{P1:3} encode the charging and network constraints defined in Section~\ref{pfacnf} respectively.
For notational convenience,	we define the feasible set as $\mathcal{S}$, the set of all optimal solutions $(X^*, A^*)$ as $\mathcal{S}^*$, the set of all $X^*$ as $\mathcal{X}^*$, and the set containing all $A^*$ as $\mathcal{A}^*$.

\section{Preliminaries} \label{ppop}
In this section, we discuss existing approaches for solving the phase optimization formulation~\eqref{P1:0}-\eqref{P1:3}.

\subsection{Existence, non-uniqueness and non-convexity}

	\begin{theorem}[Existence and non-uniqueness]\label{thm:1} The set of optimal solutions $\mathcal{S}^*$ for the optimization~\eqref{P1:0}-\eqref{P1:3} is non-empty and it is a non-singleton set.
	\end{theorem}
	\begin{proof}
		(1) Existence: Given any $X \in \mathcal{X}$, let $\mathcal{A}(X)$ be a set that contains all $A$ such that $A \in \mathcal{C}_{\mathrm{r}} \cap \mathcal{C}_{\mathrm{d}}$ and $(X, A) \in \mathcal{C}_{\mathrm{soc}}$. Since $f_{\mathbf{1}}$ is continuous in $\mathbb{R}^{N \times T}$, it is continuous in $\mathcal{A}(X)$. Moreover, $\mathcal{A}(X)$ is a compact set in $\mathbb{R}^{N \times T}$. According to the Weierstrass extreme value theorem, the infimum of $f_{\mathbf{w}}$ must be attained in $\mathcal{A}(X)$, denoted as $A({X}^*)$. Since this is true for any $X$ and $\mathcal{X}$ is countably finite, $\exists X^* \in \mathcal{X} \text{ s.t. } f_{\mathbf{1}}(A^*({X^*})) \geq f_{\mathbf{1}}(A({X}^*)), \forall X \in \mathcal{X}$. Therefore, $(X^*, A^*({X^*}))$ attains the optimal value. (2) Non-uniqueness: Given any optimal $(X^*, A^*)$, one can verify that any permutation of rows of $X^*\in\mathcal{X}$ is also optimal.
	\end{proof}
    Despite the existence of an optimal solution, finding such a pair $(X^*, A^*)$ is an intractable task because (1) the problem~\eqref{P1:0}-\eqref{P1:3} is an MIP since $X \in \mathbb{Z}_{+}^{3 \times N}$; (2) the set $\mathcal{C}_{\mathrm{soc}}$ is on $(X, A)$. Next, we discuss several existing approaches to solve the optimization~\eqref{P1:0}-\eqref{P1:3} and argue they are not efficient enough for our phase optimization problem.

\subsection{Existing approaches}
\subsubsection{Brute-force interior-point method}\label{socpx}
We first consider a direct approach for solving~\eqref{P1:0}-\eqref{P1:3}.
	Given a fixed phase selection matrix $X \in \mathcal{X}$, the optimization~\eqref{P1:0}-\eqref{P1:3} is a second-order-cone program (SOCP) formulated in \cite{lee2020adaptive}. Solving the SOCP gives us an exact optimal solution in polynomial time by the interior-point method, according to \cite{Kuo2004}. We denote this SOCP algorithm as $\mathsf{SOCP}({X})$ given $X\in\mathcal{X}$. Since $\mathcal{X}$ is countably finite, a global optimal solution $(X^*, A^*)$ can be exactly found by iterating over all $X \in \mathcal{X}$ in a  brute-force way. We denote this brute-force algorithm as $\mathsf{BFSOCP}$. Due to the brute-force approach,  $\mathsf{BFSOCP}$ is not an efficient method, especially when $N$ and $T$ are large.
	\begin{remark} [EXPTIME]\label{rk:exp}
If the time complexity of $\mathsf{SOCP}$ is $O(p(NT))$, where $p(\cdot)$ is a polynomial function, then the time complexity of $\mathsf{BFSOCP}$ is  $O(3^{p(NT)})$.
	\end{remark}

\subsubsection{Approximation by simulated annealing}\label{sec:absa}
Another approach to solving the optimization~\eqref{P1:0}-\eqref{P1:3} is to approximate $(X^*, A^*)$ by simulated annealing, which has been applied to combinatorial problems such as the $q$-coloring problem and the traveling salesman problem. Although this heuristic approach is more efficient than the $\mathsf{BFSOCP}$, it still requires a long time to converge or escape from the local minima. Therefore it becomes unsuitable for the online implementation of our phase optimization problem. In Section~\ref{sec:exp}, we will demonstrate the performance of simulated annealing against other methods.

\section{Phase Optimization Algorithms}
	\label{opscad}

In this section, we present our algorithms for the phase optimization problem defined in \eqref{P1:2}-\eqref{P1:3}.

\subsection{PXA algorithm for phase optimization}
\label{sofpo}

As discussed in Section~\ref{socpx}, if we have an approach to obtain some optimal phase selection matrix $X^*$ effectively, then we obtain an $(X^*, A^*) \in \mathcal{S}^*$ by solving the $\mathsf{SOCP}({X^*})$ defined in Section~\ref{socpx}. Motivated by this, we propose the following bi-level optimization to generate our estimate of $X^*$, which can be regarded as a charging constraint relaxation of the non-convex optimization~\eqref{P1:0}-\eqref{P1:3}:
	\begin{mini!}|s|[2]
		{\begin{aligned}
				& X \in \mathbb{Z}_{+}^{3 \times N}\\
				& P \in \mathbb{R}_{+}^{3 \times T}\\
		\end{aligned}}
		{g_{\mathbf{1}}(P) \label{P2:0}}
		{\label{P:2}}
		{}
		\addConstraint{X}{\in \mathcal{X}, \label{P2:1}}
		\addConstraint{PM}{\leq XW, \label{P2:2}}
		\addConstraint{\left|\left[
			\begin{array}{c}
				\Phi_1 \\
				\hline
				\Phi_2
			\end{array}
			\right]P \right|}{\leq C_{\max}. \label{P2:3}}
	\end{mini!}
In~\eqref{P2:0}-\eqref{P2:3}, the variable $P \in \mathbb{R}_{+}^{3 \times T}$ is  an  \textit{aggregated phase power matrix}, whose row vector is the power charged at $t \in \mathcal{T}$ at each phase. To minimize charging loss as we minimize $f_{\mathbf{1}}$ in~\eqref{P1:0}-\eqref{P1:3}, the aggregated power function $g_{\mathbf{1}}: \mathbb{R}^{3 \times T} \to \mathbb{R}$ is defined as
$
		g_{\mathbf{1}}(P) \vcentcolon= -\sum^{N}_{i = 1} \sum^{T}_{j = 1} P_{(i,j)}.
$ Two matrices
	$M \in \mathbb{R}_{+}^{T \times K}$ and $W \in \mathbb{R}_{+}^{N \times K}$ in~\eqref{P2:2} form $K$ linear constraints over $P$ and $X$ as a relaxation of the charging constraint in~\eqref{P1:2}. We will discuss detailed constructions of $M$ and $W$ later in Section~\ref{sec:analysis}. The second-order cone constraint regarding $P$ can be written as~\eqref{P2:3}. We denote the optimal value of the optimization~\eqref{P2:0}-\eqref{P2:3} as $p^*_2$, and the sets containing all feasible $P$ as $\mathcal{P}$ and optimal $P^*$ as $\mathcal{P}^*$.
	\begin{remark}[\textsf{MISOCP}]
The optimization~\eqref{P2:0}-\eqref{P2:3} is a mixed-integer second-order cone programming (MISOCP) since $g_{\mathbf{1}}(P)$ is linear on $P$, ~\eqref{P2:1} is an integral polytope,~\eqref{P2:2} is a mixed-integer polytope on $P$ and $X$ and~\eqref{P2:3} is a second-order cone on $P$. Note that a MISOCP can be efficiently solved by existing optimization solvers. Moreover, we can show that the branch-and-bound method guarantees the convergence to a global optimum with mild conditions.
\end{remark}

By minimizing the variables $P,X$ and $A$ of the bi-level optimizations~\eqref{P2:0}-\eqref{P2:3} and~\eqref{P1:0}-\eqref{P1:3}, we present our \textit{PXA algorithm} in Algorithm~\ref{alg:PXA}. As a MISOCP solver, the computational complexity of PXA is exponential in $N$. However, in practice, the PXA computes much faster than both $\mathsf{BFSOCP}$ and Simulated annealing, especially in the online setting.

\begin{algorithm}[h]
\DontPrintSemicolon
\SetKwBlock{Begin}{function}{end function}
\KwIn{ACN parameters $\mathsf{ACN}=\{D, \mathcal{C}_{\mathrm{r}}, \mathcal{C}_{\mathrm{d}}, \mathcal{C}_{\mathrm{soc}}\}$}
\KwOut{Phase selection $X$ and charging decision $A$}
Initialize $M$ and $W$\\
\Begin($\text{PXA} {(}\mathsf{ACN}|M, W {)}$){
\textit{Step 1}: Solve~\eqref{P2:0}-\eqref{P2:3} with $M, W$ and obtain $X$\\
\textit{Step 2}:
Solve $\mathsf{SOCP}({X})$ and obtain $A$}
\caption{{PXA} for Phase Optimization}
\label{alg:PXA}
\end{algorithm}

	\subsection{MPC-based PXA Algorithm for online phase optimization }\label{sopo}
In real-world EV charging applications, we need to implement the PXA algorithm shown in Algorithm~\ref{alg:PXA} in an online setting, as presented below.

	 Given any $t \in \mathcal{T}$, we observe a fleet of EVs which has already been connected to the ACN (old) and a group of incoming EVs (new). Denote the number of those old and new EVs by $N_{\mathsf{old}}, N_{\mathsf{new}}$. The number of all observable EVs at time $t$ satisfies $N_{\mathsf{old}} + N_{\mathsf{new}} \leq N$. We denote the fixed phase selection matrix for those connected EVs by $X_{\mathsf{old}} \in \mathbb{R}^{3 \times N_{\mathsf{old}}}$. Then $X_{\mathsf{new}} \in \mathbb{R}^{3 \times N_{\mathsf{new}}}$ for those new EVs as our decision variables.
	 We also define the \textit{EV-profile matrix} at time $t$ as $D(t)$ such that $D(t)^{(1)} = \max(t\mathbf{1}_{N\times 1}, D^{(1)})$, $D(t)^{(2)} = \min(D^{(1)} + D^{(2)} - t, D^{(2)})$ and $D(t)^{(3)} = D^{(3)} - \Delta\mathbf{e}(t)$, where $\Delta\mathbf{e}(t)$ is a vector whose entry $i$ is the energy sum of EV $i$ that has been charged in $\{1,\cdots, t - 1\}$. This online problem is essentially the same as the optimization~\eqref{P2:0}-\eqref{P2:3} except that a subset of columns of $X$, i.e., $X_{\mathsf{old}}$ has been fixed.
	 Using the idea of MPC, the online MPC-based PXA algorithm ({PXA-MPC}) described above is presented in Algorithm~\ref{alg:1}.
	 \begin{algorithm}
	 	\caption{{{PXA-MPC}} (Online Phase Optimization)}\label{alg:1}
	 	\KwIn{$\mathsf{ACN}= \{D, \mathcal{C}_{\mathrm{r}}, \mathcal{C}_{\mathrm{d}}, \mathcal{C}_{\mathrm{soc}}\}$}
	 	\KwOut{Phase selections $\{X(t)|t\in\mathcal{T}\}$ and charging decisions $\{A(t)|t\in\mathcal{T}\}$}
	 	Initialize
	 	$A^*(0) \gets \mathbf{0}_{N \times T}$,
	 	$X(0) \gets \mathbf{0}_{3 \times N}$, $M$ and $W$\\
	 	\For{$t \in \mathcal{T}$}{

	 		\textit{Step 1}: Solve~\eqref{P2:0}-\eqref{P2:3} with  $M$, $W$, $D(t)$, $X_{\mathsf{old}}$, $\mathsf{ACN}$ and obtain  $X_{\mathsf{new}}$

	 		\textit{Step 2}: Update $X(t) = [X_{\mathsf{old}} |  X_{\mathsf{new}}] $

	 		\textit{Step 3}: Solve $\mathsf{SOCP}({X})$ with $D(t), X(t), \mathsf{ACN}$. Obtain a solution $A(t)$

	 		\textit{Update system states}:

	 		$D^{(2)} \gets \min(D^{(1)} + D^{(2)} - t - 1, D^{(2)})$

	 		$D^{(1)} \gets \max((t + 1)\mathbf{1}_{N\times 1}, D^{(1)})$

	 		$D^{(3)} \gets D^{(3)} - {A(t)}^{(1)} \Delta_T$
	 	}
	 \end{algorithm}

	\section{Theoretical Results}
	\label{sec:analysis}
In this section, we provide an analysis of the PXA algorithm, design of the matrices $M$ and $W$, and show that the solution given by the PXA algorithm is optimal under certain laxity conditions.

	\subsection{Constructions of matrices $M$ and $W$}\label{nzlcs}
We propose the following constructions of $M$ and $W$ as a charging constraint relaxation to approximate an optimal $X^*$.

Define $M(m,T) \in \{0, 1\}^{n \times {T \choose m}}$ a  binary \textit{$(m,T)$-selection matrix} if all column vectors of $M(m,T)$ form all possible combinations of selecting $m$ out of $T$ entries. If the $i$-th entry is selected by the $j$-th column, $M_{(i,j)}(m,T)=1$ and $M_{(i,j)}(m,T)=0$ otherwise.
A binary $(m,T)$-selection matrix is used to represent the aggregated phase power for time steps in $\mathcal{T}$. For instance, as shown in Figure~\ref{fig:pxa}(a) with $T = 6$, a charging constraint of $P$ can be described with a column vector $[1, 0, 1, 0, 1, 0]^{\top}$ of $M(3,6)$. Similarly, for all $m\in \mathcal{T}$, we will use every column of $M(m,T)$ to describe a charging constraint in the optimization~\eqref{P2:0}-\eqref{P2:3}. Denote the concatenation of $M(m,T)$ by $M_{\mathcal{T}}$:
	\begin{equation}
\nonumber
		M_{\mathcal{T}} \vcentcolon= \left[	M(1,T) | M(2,T) | \cdots | M(T,T) \right] \in \{0, 1\}^{T \times 2^{T}-1}.
	\end{equation}
Moreover, define $2^{T}-1$ horizontal concatenations of power demand vector $D^{(3)}(\Delta_T^{-1})$ as $V \in \mathbb{R}^{N \times 2^{T}-1}$ such that $V\vcentcolon= [D^{(3)}(\Delta_T^{-1}) | \cdots | D^{(3)}(\Delta_T^{-1})]$. Denote $C = \min\{A, B\}$ as the minimum-element operation on two matrices $A$ and $B$, e.g., for any $i,j, C_{(i,j)} = \min\{A_{(i,j)}, B_{(i,j)}\}$. Now we can define a set of $P$ as
	\begin{equation}
	\nonumber
		\begin{aligned}
			&\mathcal{P}_3 \vcentcolon= \{P \in \mathbb{R}_{+}^{3 \times T} | PM_{\mathcal{T}} \leq XW_3\}, \text{ where }\\
			&   W_3 \vcentcolon= \min\{EM_{\mathcal{T}}r_{\max}, V\}
		\end{aligned}
	\end{equation}
so that the aggregated power must be bounded by both the sum of chargers' maximum rates and the total energy demand of EVs that exist in the selected time steps.
\begin{remark}
We note the following relations:
	\begin{enumerate}
		\item The maximum-rate constraint of $P$ from ACN chargers, i.e. $PI_T \leq XEr_{\max}$, is a superset of $\mathcal{P}_3$, since every column of $I_T$ is a column of $M_{\mathcal{T}}$;
		\item For every feasible point $(X, A)$ in the optimization~\eqref{P1:0}-\eqref{P1:3}, $(X, P)$ is feasible in $\mathcal{P}_3$, where $P = XA$. $(X, P)$ also satisfies the second-order cone constraint, thus being feasible in the optimization~\eqref{P2:0}-\eqref{P2:3}.
	\end{enumerate}
\end{remark}

Despite that characterizing the constraint~\eqref{P2:2} with $P \in \mathcal{P}_3$ forms an acceptable relaxed constraint as it considers both the charger's maximum rate and the EVs' energy demands, $M_{\mathcal{T}}$ has exponentially many ($2^{T}-1$) columns so that $P\in \mathcal{C}_p$ contains more than $2^{T}-1$ linear constraints. To reduce the number of constraints, we can carefully select $m \ll 2^{T}-1$ columns of $M_{\mathcal{T}}$ to form a relaxation of Constraint~\ref{P2:2}. In our experiments,  we set $\tilde{M} \vcentcolon= \left[I | E^T | \mathbf{1}_{T \times 1}\right]$ ($m = T + N + 1$) and its corresponding $W_2$ as a sub-matrix of $W_3$ by selecting the corresponding columns in  $W_3$. Figure~\ref{fig:pxa} (b) shows the relation of the feasible sets with different choices of $M$ and $W$. Although we cannot guarantee  $X \in \mathcal{X}^*$ by using $\tilde M$, we demonstrate in Section~\ref{sec:exp} the result given by the PXA algorithm using $\tilde M$ generates a near-optimal solution. Moreover, in the next section, we prove that in the zero laxity case, considering the relaxation $\mathcal{P}_1$ shown in Figure~\ref{fig:pxa}(b), solving the optimization~\eqref{P2:0}-\eqref{P2:3} by setting $M = I_T$ and $W = Er_{\max}$  gives an optimal $X^*\in \mathcal{X}^*$.

	\begin{figure}
	\centering
		\includegraphics[width=0.43\textwidth]{python_images/constraint_sets_example.png}
		\caption{(a) shows an example of charging constraint formed by summing up the phase power at time step $t=1,3$ and $5$. Every two EV sessions are allocated at one phase and their energy demand equals to a proportion of $r_{\max}$. (b) shows the relations of the feasible sets $\mathcal{P}_1,\mathcal{P}_2$ and $\mathcal{P}_3$ formed by different constructions of $M$.}
		\label{fig:pxa}
	\end{figure}

	\subsection{Zero laxity analysis}
	\label{sec:zero-laxity}
Consider the case when all EV sessions have zero laxity, i.e., $\Delta_T r_{\max} d_i = D_{(i)}$, for any  $i \in \mathcal{N}$. Note that the zero-laxity case frequently occurs in specific environments, e.g., fast-charging service centers in roadside rest areas or shopping plazas, wherein most EV drivers leave upon or before their EVs complete charging, thus making EV charging sessions almost zero-laxity. Denote by $I_T$ a $T\times T$ identity matrix. The following theorem holds.

		\begin{theorem}[Global optimality]\label{thm:l0o}
Suppose $M = I_T$, $W = Er_{\max}$ and $\Delta_T r_{\max} d_i = D_{(i)}$ for any $i \in \mathcal{N}$. The optimal solution $X$ of the optimization~\eqref{P2:0}-\eqref{P2:3} satisfies that $X \in \mathcal{X}^*$.
		\end{theorem}

		\begin{proof}
			It is equivalent to prove that, given any $P^* \in \mathcal{P}^*$, there exists a pair $(X, A) \in \mathcal{S}^*$ such that $P^* = XA$. For any $P^* \in \mathcal{P}^*$, there must exist an $X \in \mathcal{X}$ such that the constraint~\eqref{P2:2} is satisfied. We construct a matrix $A \in \mathbb{R}^{N \times T}$ with the following steps: (1) denote the column sum of $X$ as $\mathbf{n} \vcentcolon= [n_1, n_2, n_3]^{\top} = \sum^{N}_{i = 1} X^{(i)}$; (2) for all $i \in \mathcal{T}, j \in \mathcal{N}$,
			\begin{equation}
				\label{eqn:A_c}
				A_{(i, j)} =
				\begin{cases}
					\frac{P^*_{(m, j)}}{n_m}      & \ \text{if } X_{(m, i)} = 1 \text{, where } m \in \{1,2,3\},\\
					0  & \ \text{otherwise}.
				\end{cases}
			\end{equation}
			Notice that $P^* = XA$ and $P^*I \leq XEr_{\max}$. For all $m \in \{1,2,3\}$ and $j \in \mathcal{T}$, we can verify that $P^*_{(m, j)}\leq n_m r_{\max}$, implying that $A_{(i, j)} \leq r_{\max}$. Thus, $A \in \mathcal{C}_{\mathrm{r}}$. Since $\mathcal{C}_{\mathrm{r}} = \mathcal{C}_{\mathrm{r}} \cap \mathcal{C}_{\mathrm{d}}$ and the constraint~\eqref{P2:3} implies $(X, A) \in \mathcal{C}_{\mathrm{soc}}$, we conclude that $(X, A)$ is feasible in the optimization~\eqref{P1:0}-\eqref{P1:3}. We still need to show $p^*_1 = f_{\mathbf{1}}(A)$ to complete the proof. We show this by deriving a contradiction and suppose that there exists some $\tilde A$ that is feasible in the optimization~\eqref{P1:0}-\eqref{P1:3} such that $f_{\mathbf{1}}(\tilde A) < f_{\mathbf{1}}(A)$, thus $p^*_1 \neq f_{\mathbf{1}}(A)$. Since $\tilde A$ is feasible, then there exists an $\tilde X \in \mathcal{X}$ such that $(\tilde X, \tilde A) \in \mathcal{C}_{\mathrm{soc}}$. If we denote $\tilde P \coloneqq \tilde X \tilde A$, it implies $\tilde P$ satisfies the constraint~\eqref{P2:3}. In addition, $\tilde A \in \mathcal{C}_{\mathrm{r}}$ leads to $\tilde X \tilde A \leq \tilde X Er_{\max}$, and hence $\tilde P I_T \leq \tilde X Er_{\max}$. Combining both, $(\tilde X, \tilde P)$ is feasible in the optimization~\eqref{P2:0}-\eqref{P2:3} and $g_{\mathbf{1}}(\tilde P) = f_{\mathbf{1}}(\tilde A) < f_{\mathbf{1}}(A) = g_{\mathbf{1}}(P^*)$. However, $P^*$ is an optimal solution that attains $p^*_2$, implying $g_{\mathbf{1}}(\tilde P) < p^*_2$. Therefore, the existence of $\tilde A$ leads to a contradiction. As a result, $f_{\mathbf{1}}(A) = p^*_1$ implies  that $(X, A) \in \mathcal{S}^*$.
		\end{proof}
	Theorem~\ref{thm:l0o} shows that solving the optimization~\eqref{P2:0}-\eqref{P2:3} gives an exact $X\in\mathcal{X}^*$ in the zero-laxity charging scenarios, thus solving the optimization~\eqref{P1:0}-\eqref{P1:3}. The following theorem shows that solving the optimization~\eqref{P2:0}-\eqref{P2:3} by the branch-and-bound method guarantees the convergence to a global optimum with a mild condition  $p^*_1 \Delta_T = -\sum^N_{i=0}D_{(i, 3)}$ that ensures the feasibility of EVs' energy demands.
	\begin{theorem}\label{thm:eos}
	Suppose $\Delta_T r_{\max} d_i = D_{(i)}$ for any $i \in \mathcal{N}$.	If $p^*_1 \Delta_T = -\sum^N_{i=0}D_{(i,3)}$, solving the optimization~\eqref{P2:0}-\eqref{P2:3} with  $M = I_T$, $W = Er_{\max}$ by branch-and-bound attains $p^*_1$.
	\end{theorem}
	\begin{proof}
		Since Theorem~\ref{thm:l0o} proves that the optimal solution $X$ of the optimization~\eqref{P2:0}-\eqref{P2:3} belongs to $\mathcal{X}^*$, it is enough to show that the branch-and-bound always obtains such an $X$. Since $p^*_1 \Delta_T = -\sum^N_{i=0}D_{(3,i)}$, there exists a $(X^*, A^*)$ feasible in the optimization~\eqref{P1:0}-\eqref{P1:3} such that $f_{\mathbf{1}}(A^*) = p^*_1$. From the proof of Theorem~\ref{thm:l0o}, $(X^*, X^*A^*)$ is feasible in the optimization~\eqref{P2:0}-\eqref{P2:3} and $p^*_2 = g_{\mathbf{1}}(X^*A^*)$. Therefore, we have an integral optimal solution when solving the SOCP relaxation at every branching process. Furthermore, the branch-and-bound process on $X$ is equivalent to a ternary-tree search due to $X \in \{0, 1\}^{3 \times N}$ and $\sum^3_{i = 1} X_{(i, j)} = 1$, for any $j \in \mathcal{N}$. Specifically, at each level $j$ of this ternary-tree, the algorithm computes an optimal solution of the SOCP relaxation with fixes the $j$-th EV's phase selection as $ab, bc, ca$ respectively. We can verify that (i) the branching trajectories of such a ternary-tree are towards all $X \in \mathcal{X}$, thus including the trajectory towards $X^*$; (ii) since $p^*_2$ equals to the total power demand, the bounding criteria is always no less than $p^*_2$, which guarantees that the branching trajectory towards $X^*$ will not be pruned. Combining (i) and (ii), we conclude that the branch-and-bound converges to $p^*_2$, which equals to $p^*_1$.
	\end{proof}

  	\section{Experimental Results} \label{sec:exp}
	\subsection{Experiment setups}
	 We use CVXPY~\cite{diamond2016cvxpy} for modelling the optimization. GUROBI~\cite{gurobi} is used as the branch-and-bound solver for the optimization~\eqref{P2:0}-\eqref{P2:3}, and ECOS~\cite{domahidi2013ecos} for solving the SOCP. We use the ACN-data~\cite{li2019acndata}, which is a dataset collected from a real-world ACN at Caltech. The values of ACN parameters in $\mathsf{ACN}=\{D, \mathcal{C}_{\mathrm{r}}, \mathcal{C}_{\mathrm{d}}, \mathcal{C}_{\mathrm{soc}}\}$ are summarized in Table~\ref{tab:1}. We use  $M = \tilde M$ and $W = W_2$ defined in Section~\ref{nzlcs} in our experiments. We split EV sessions overnight and treat all EV sessions in one day as a charging episode. The charging data $D$ is selected in the pre-COVID period, from date 2018 May. 1st to 2018 Jul. 31st ($91$ episodes in total).

	\subsection{Benchmarks}
	 We consider four benchmarks. (1) {(EV RDM.)} The phases decided by the incoming EV drivers as in the ACN database. (2) {(UNI. RDM.)} We draw each phase connection uniformly at random from $\{ab, bc, ca\}$. (3) {(RRB. RDM.)} We apply the round-robin method based on the arrival time. (4) {(WST.)} The most unbalanced phase selection by allocating all sessions to one phase line. (5) {(M.C.S.A.)} The Markov chain simulated annealing method introduced in Section~\ref{sec:absa} with $10000$ iterations per day.
	 \begin{figure}
	 	\centering
	 	\begin{subfigure}[b]{0.47\textwidth}
	 		\centering
	 		\includegraphics[width=\textwidth]{python_images/test_cases.png}
	 		\caption{}
	 	\end{subfigure}
	 	\hfill
	 	\begin{subfigure}[b]{0.47\textwidth}
	 		\centering
	 		\includegraphics[width=\textwidth]{python_images/maximum_energy_charge.png}
	 		\caption{}
	 	\end{subfigure}
	 	\hfill
	 	\caption{(a) shows $9$ EV charging instances for various phase optimization algorithms; (b) highlights the performance of the {PXA-MPC} (Algorithm~\ref{alg:1}), compared with the charging outcomes for various phase optimization algorithms.}
	 	\label{fig:onoff}
	 \end{figure}

    \begin{table}[!ht]

    \renewcommand{\arraystretch}{1.3}
    \centering
    \caption{Parameters of the ACN used in our experiments.}
    \label{tab:1}
    \begin{tabular}{|c|c|c|c|c|c|c|}
    \hline
        Episode &  $r_{\max}$& $\mathcal{C}_{\max}$& $V$ & $n_r$& $T$& $\Delta_T$ \\
    \hline
        Day & $3\mathrm{kW}$ & $20\mathrm{kW}$ & $220\mathrm{V}$ & $4$ & $24\mathrm{hrs}$ & $0.2\mathrm{hr}$ \\
    \hline
    \end{tabular}
    \end{table}

\subsection{Near-optimality of the PXA algorithm}
We implement the PXA algorithm in Algorithm~\ref{alg:PXA} for $9$ toy charging scenarios with $T=12$.
The first $4$ scenarios are zero-laxity charging (see Section~\ref{sec:zero-laxity}), and the rest are non-zero laxity charging scenarios. Knowing the ACN parameters $\mathsf{ACN}=\{D, \mathcal{C}_{\mathrm{r}}, \mathcal{C}_{\mathrm{d}}, \mathcal{C}_{\mathrm{soc}}\}$, the PXA algorithm always satisfies $100\%$ of the energy demands, as shown in Figure~\ref{fig:onoff} (a). The uncharged energy demands for {WST.}, {UNI. RDM.}, {RRB. RDM.} are $50\%$, $19.1\%$ and $9.2\%$ respectively.

\subsection{Energy efficiency of the {PXA-MPC} algorithm}
We implement the {PXA-MPC} (Algorithm~\ref{alg:1}) with the quick-charge function defined in Section~\ref{sec:pop} and assume the controller only gets to know charging sessions when their arrival times are no $\Delta_T$ larger than the current time $t$. In Figure~\ref{fig:onoff} (b), regarding the online performance (solid), the {PXA-MPC} algorithm charges an additional $7.52\%$ of weekly total energy demand compared against the {EV RDM.}; $2.86\%$ compared against the {UNI.RDM} and $1.34\%$ compared against the {RRB. RDM.}. Offline performance (dashed) is also considered. In particular, we observe that the {PXA-MPC} algorithm is near-optimal. Moreover, the {PXA-MPC} takes on average $343.52$ ms to compute at each time step with a standard deviation of $120.42$ ms. In Figure~\ref{fig:balance} (a), we see the {EV RDM.} (top) induces an unbalanced three-phase energy dispatch. While the {PXA-MPC} algorithm (bottom) allocates the energy equally to each phase, especially when the EV demands are high from $10$ am to $8$ pm. We further demonstrate that the {PXA-MPC} algorithm is capable of satisfying more energy demands compared with other phase optimization methods. In Figure~\ref{fig:balance}(b), we show the trade-off between the line capacity and demand satisfaction. We observe that with the same line capacity, the {PXA-MPC} algorithm achieves the highest energy efficiency. In practice, the phases of an ACN can be highly unbalanced, and the {WST.} corresponds to the most unbalanced case when all energy demands are allocated to a single phase.

	\begin{figure}
	 	\centering
	 	\begin{subfigure}[b]{0.45\textwidth}
	 		\centering
	 		\includegraphics[width=\textwidth]{python_images/phase_energy.png}
	 		\caption{}
	 	\end{subfigure}
	 	\hfill
	 	\begin{subfigure}[b]{0.45\textwidth}
	 		\centering
	 		\includegraphics[width=\textwidth]{python_images/capacity_variation.png}
	 		\caption{}
	 	\end{subfigure}
	 	\hfill
	 	\caption{(a) shows the average energy charged over $91$ episodes at phases $\{ab,bc,ca\}$. The top figure shows the daily trajectory from EV drivers' phase selections and the bottom figure shows the distribution of energy among phases using the {PXA-MPC}; (b) shows demand satisfaction rates with different line-power limit  $C_{\max}$.}
	 	\label{fig:balance}
	 \end{figure}

	\subsection{Electricity cost minimization}
    Since future information on the energy price is usually unknown, we assume no knowledge of the price vector and use the quick charge function to compute $X^*$ instead. We demonstrate how the {{PXA-MPC}} algorithm can be used to minimize electricity costs. We fix the line capacity as $\mathcal{C}_{\max} = 30$. It can be seen in Figure~\ref{fig:price}, the {{PXA-MPC}} satisfies $99.86\%$ of the energy demand, while the {EV RDM.}, {UNI. RDM.} and the {RRB. RDM.} meet $96.29\%$, $98.85\%$ and $99.51\%$ respectively. Furthermore, the average electricity price of the {{PXA-MPC}} algorithm is $33.85$ cents per kWh, which is constantly lower than that of the {EV RDM.} ($35$ cents/kWh), the {UNI. RDM.} ($34.3$ cents/kWh), and the {RRB. RDM.} ($34.09$ cents/kWh). Compared with the {M.C.S.A}, the {{PXA-MPC}} tends to satisfy higher charging demand in lieu of saving more costs. The phase selection matrix $X^*$ computed from the quick charge function helps the {{PXA-MPC}} satisfy a higher charging demand than that of the {M.C.S.A}. In conclusion, the {{PXA-MPC}} serves more but costs less compared with other online methods.
	\begin{figure}
		\centering
		\includegraphics[width=0.45\textwidth]{python_images/minimum_price_charge.png}
		\caption{Comparison of  weekly average prices (solid line, left $y$-axis) for {EV RDM.} (red), {UNI. RDM.} (green) and {{PXA-MPC}} (blue), together with their energy delivery in percentage (dashed line, right $y$-axis).}
		\label{fig:price}
	\end{figure}
	\section{Conclusion and Future Directions} \label{c}
		In this paper, we formulate a phase optimization problem for EV charging. We provide a phase optimization algorithm, the PXA algorithm (Algorithm~\ref{alg:PXA}) and an online implementation, the {PXA-MPC} algorithm (Algorithm~\ref{alg:1}). We prove the optimality of the PXA algorithm for the zero-laxity case and demonstrate its efficacy on real-world EV charging data.
Our research opens three interesting directions. Note that our {PXA-MPC} algorithm selects $T + N + 1$ out of $2^T - 1$ linear constraints to solve a linear relaxation of the charging constraint. First, it would be interesting to verify if an optimal solution can be obtained with all linear constraints, i.e., if $\mathcal{P}_3=\overline{\mathcal{P}}$. Second, in addition to the sub-matrix $\tilde M$ defined in Section~\ref{nzlcs}, finding other ways to select those linear constraints and figuring out their optimality would be useful. Finally, it is worth investigating more applications of p

{\footnotesize


\section{Introduction}
The study of Noether symmetries and corresponding first integrals of underlying dynamical system has important mathematical significance as well as it has profound physical background. Discussion of Lie point transformations which leave the action integral of the damped harmonic oscillator (DHO) invariant has widened in recent years. Among the enormous and useful results concerning DHO, the discussion of Noether symmetries and the relevant conserved quantities have attracted attention of many researchers. If a Lie point transformation leaves the action integral invariant, the classical Noether's theorem \cite{Noe} provides a corresponding first integral. Another significance of these conserved quantities is that they help us in establishing solutions of underlying systems. Moreover, Lagrangian and Hamiltonian functions play essential role in deriving equations of motion. Therefore, in order to study relation of symmetries and conservation laws using Noether's theorem, one needs the corresponding equation of motion that arise from the action principle \cite{Gold}. A number of studies have been reported to construct integral of motions of DHO \cite{Cho,And,Cer}. In the current study, we intend to investigate the symmetry properties of DHO equation that is derivable from a variational principle. We shall see that the full eight parameter group that leaves the equation of motion invariant, only five of them leave the action integral invariant and become Noether symmetries. Consequently, we obtain only five conserved quantities for DHO. It is commonly known that maximum number of conserved quantities of a dynamical system often leads to find the solution of underlying equation. In this paper, we find many results on Noether symmetries and conservation laws of DHO and point out some interesting results related to Noether symmetries and first integrals to that already available in the literature. We also present and discuss a new form of Lagragian (may be taken as standard Lagragian) for DHO system which has not been reported so far. Though the new form of Lagrangian is completely different yet it yields same results as derived from the Bateman Lagrangian. So this new form of optimal Lagragian function may conceal something interesting for the physicists.\\
A brief outline of the present work is as follows: In Section II, we outline our treatment on Noether symmetries and Noether's theorem. In Section III, we give a classification of DHO with respect to Noether symmetries and construct correspondent conserved quantities for three cases. We provide Lie algebras of Noether operators and also discuss new form of Lagrangian for DHO. Finally, we present a brief summary of our work in Section IV.
\section{Preliminaries on Noether Symmetries and First integrals}
In order to present results of the current study in a straightforward manner, we start with reviewing expressions on finding Noether symmetry generators and associated conservation laws for DHO. We shall use the salient features of this section to derive important results in the next section. \\
Consider the point type vector field
\begin{eqnarray}
Z=\xi(t,u)\partial_{t}+\eta(t,u)\partial_{u},
\end{eqnarray}
then its first prolongation can be expressed as
\begin{equation}
Z^{[1]}=Z+(\eta_{t}+\eta_{u}u^{\prime}-\xi_{t}u^{\prime}-\xi_{u}u^{\prime
	2})\partial_{u^{\prime}}.
\end{equation}
Further assume that we have a first-order Lagrangian function $L(t,u,u')$ whose insertion in the Euler-Lagrange equation
\begin{eqnarray}
\frac{d}{dt}(\frac{\partial L}{\partial u^{\prime}})-\frac{\partial
	L}{\partial u}=0,
\end{eqnarray}
yields the following equation of motion
\begin{eqnarray}
u^{\prime\prime}=w(t,u,u^{\prime}).
\end{eqnarray}
We recall that the vector field $Z$ is called a Noether point symmetry generator corresponding to the Lagrangian $L(t,u,u')$ if a gauge function $B(t,u)$ exists such that
\begin{eqnarray}
Z^{[1]}(L)+D(\xi)L=D(B),
\end{eqnarray}
holds, where
\begin{eqnarray}
D=\partial_{t}+u^{\prime}\partial_{u}+u^{\prime\prime}\partial_{u^{\prime}}+...  .
\end{eqnarray}
The significance of a Noether symmetry is that it provides the first integral as well as double reduction of order of underlying equation of motion.\\
\textbf{Noether's Theorem:} Suppose that a Lie point symmetry generator $Z$ satisfies the the condition given in Eq. (5), then the following expression leads to the conserved quantity or constant of motion of (4) related to $Z$.
\begin{eqnarray}
I=\xi L+(\eta-u^{\prime}\xi)\frac{\partial L}{\partial
	u^{\prime}}-B.
\end{eqnarray}
\section{Noether symmetries and first integrals of damped harmonic oscillator}
In this section we wish to find Noether symmetries and then utilizing them in the Noether's theorem to determine first integrals of a famous dynamical equation depicting motion of DHO. We shall classify the problem with respect to Noether's symmetries and construct the associated conserved quantities for over-damped, under-damped and critically damped cases and also determine their Lie algebras. So for this we proceed as:\\
Consider the physical system of DHO whose dynamics is governed by the following second-order ordinary differential equation
\begin{eqnarray}
u^{\prime\prime}+\frac{c}{m}u^{\prime}+\frac{k}{m}u=0.\label{eq1}
\end{eqnarray}
The well known time dependent Lagrangian of above Eq. (8) is given by the so called Bateman Lagrangian \cite{Bat,Yan,Edw,Ray}
\begin{eqnarray}
L=\frac{1}{2}\exp\left({\frac{ct}{m}}\right)(mu^{\prime 2}-ku^{2}),\label{lag1}.
\end{eqnarray}
From Eq. (9), it is evident that the time dependency of Lagrangian does not reveal the energy conservation of the system as one observes in the usual harmonic oscillator case. After inserting the above Lagrangian in eq. (5) and solving the resulting set of partial differential equations, we obtain the following Noether symmetries
\begin{equation}\label{p13}
\begin{aligned}
\textbf{X}_{1}&=\frac{\partial}{\partial t}-\frac{cu}{2m}\frac{\partial}{\partial u},\\
\textbf{X}_{2}&=\sin\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)\frac{\partial}{\partial t}
+\Big(\frac{u\sqrt{4km-c^2}}{2m}\cos\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)-\frac{cu}{2m}\sin\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)\Big)\frac{\partial}{\partial u},\\
\textbf{X}_{3}&=\cos\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)\frac{\partial}{\partial t}-\Big(\frac{u\sqrt{4km-c^2}}{2m}\sin\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)+\frac{cu}{2m}\cos\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)\Big)\frac{\partial}{\partial u},\\
\textbf{X}_{4}&=\exp\left(\frac{(-c-\sqrt{c^2-4km})t}{2m}\right)\frac{\partial}{\partial u},\\
\textbf{X}_{5}&=\exp\Big(\frac{(-c+\sqrt{c^2-4km})t}{2m}\Big)\frac{\partial}{\partial u}.
\end{aligned}
\end{equation}
Now by employing these operators in the Noether's theorem gives rise respectively the following 5-first integrals
\begin{equation}
\begin{aligned}
I_{1}&=-\exp\Big({\frac{ct}{m}}\Big)(mu^{\prime 2}+cuu^{\prime}+ku^{2}),\\
I_{2}&=-\frac{\exp\Big({\frac{ct}{m}}\Big)}{4m}(2m^2u^{\prime 2}+2mu(cu^{\prime}-ku)+c^2u^2)\sin\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)\\ &+\frac{\exp\Big({\frac{ct}{m}}\Big)}{4m}u\sqrt{4km-c^2}(cu+2mu^{\prime})\cos\Big(\frac{\sqrt{4km-c^2}t}{m}\Big),\\
I_{3}&=-\frac{\exp\Big({\frac{ct}{m}}\Big)}{4m}(2m^2u^{\prime 2}+2mu(cu^{\prime}-ku)+c^2u^2)\cos\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)\\
&-\frac{\exp\Big({\frac{ct}{m}}\Big)}{4m}u\sqrt{4km-c^2}(cu+2mu^{\prime})\sin\Big(\frac{\sqrt{4km-c^2}t}{m}\Big),\\
I_{4}&=\exp\Big({\frac{ct}{m}}\Big)\exp\Big({\frac{-(c+\sqrt{c^2-4km})t}{2m}}\Big)mu^{\prime}+\frac{1}{2}(c+\sqrt{c^2-4km})u\exp\Big({\frac{(c-\sqrt{c^2-4km})t}{2m}}\Big),\\
I_{5}&=\exp\Big({\frac{ct}{m}}\Big)\exp\Big({\frac{(-c+\sqrt{c^2-4km})t}{2m}}\Big)mu^{\prime}+\frac{1}{2}(c-\sqrt{c^2-4km})u\exp\Big({\frac{(c+\sqrt{c^2-4km})t}{2m}}\Big).
\end{aligned}
\end{equation}
Henceforth, we shall discuss 3-forms of DHO which include over-damped, under-damped and critical damped. We shall use conditions for each form and retrieve some known as well as some new results.\\
\textbf{Case 1. \textbf{Over-damped oscillator:}}\\ If the damping force is much stronger than the restoring force of the system, then it is called the over-damped harmonic oscillator. In this case $4km<c^2$, and the set of generators given in Eq. (10) takes the form
\begin{equation}\label{genset2}
\begin{aligned}
\textbf{X}_{1}&=\frac{\partial}{\partial t}-\frac{cu}{2m}\frac{\partial}{\partial u},\\
\textbf{X}_{2}&=\iota\sinh\Big(\frac{\sqrt{c^2-4km}t}{m}\Big)\frac{\partial}{\partial t}+\iota\Big(\frac{u\sqrt{c^2-4km}}{2m}\cosh\Big(\frac{\sqrt{c^2-4km}t}{m}\Big)-\frac{cu}{2m}\sinh\Big(\frac{\sqrt{c^2-4km}t}{m}\Big)\Big)\frac{\partial}{\partial u},\\
\textbf{X}_{3}&=\cosh\Big(\frac{\sqrt{c^2-4km}t}{m}\Big)\frac{\partial}{\partial t}-\Big(\frac{u\sqrt{c^2-4km}}{2m}\sinh\Big(\frac{\sqrt{c^2-4km}t}{m}\Big)+\frac{cu}{2m}\cosh\Big(\frac{\sqrt{c^2-4km}t}{m}\Big)\Big)\frac{\partial}{\partial u},\\
\textbf{X}_{4}&=\exp\Big(\frac{(-c-\sqrt{c^2-4km})t}{2m}\Big)\frac{\partial}{\partial u},\\
\textbf{X}_{5}&=\exp\Big(\frac{(-c+\sqrt{c^2-4km})t}{2m}\Big)\frac{\partial}{\partial u}.
\end{aligned}
\end{equation}
On closely viewing these operators, we note that the Noether operators  $\textbf{X}_{2}$ and $\textbf{X}_{3}$ attain different forms to those given in Eq. (10) while the other three generators $\textbf{X}_{1}$, $\textbf{X}_{4}$, $\textbf{X}_{5}$ and the associated conserved quantities ${I}_{1}$, ${I}_{4}$, ${I}_{5}$ remained unchanged. Hence, from the constants of the motion ${I}_{4}$ and ${I}_{5}$ we obtain the following solution for the over-damped oscillator
\begin{equation}
u(x)=\frac{2m}{c^2-4km}\Big[I_4e^\frac{-c+\sqrt{c^2-4km}}{2m}-I_5e^\frac{-c-\sqrt{c^2-4km}}{2m}\Big].
\end{equation}
Moreover, by utilizing the Noether symmetries $\textbf{X}_{2}$ and $\textbf{X}_{3}$ along with the Lagrangian (9) in Noether's theorem, we  determine respectively the following two new correspondent conserved quantities
\begin{equation}
\begin{aligned}
I_{2}&=\frac{1}{16km^2}\Big(-4km\exp\Big(\frac{ct}{m}\Big)(2m^2u^{\prime 2}+2m(cuu^{\prime}-ku^{2})+c^2u^2)\sinh\Big(\frac{\sqrt{c^2-4km}t}{m}\Big)+\\&\exp\Big(\frac{ct}{m}\Big)(cu^2+2muu^{\prime})(c^2\sqrt{c^2-4km}-(c^2-4km)^{3/2})\cosh\Big(\frac{\sqrt{c^2-4km}t}{m}\Big)\Big),\\
I_{3}&=\frac{1}{16km^2}\Big(-4km\exp\Big(\frac{ct}{m}\Big)(2m^2u^{\prime 2}+2m(cuu^{\prime}-ku^{2})+c^2u^2)\cosh\Big(\frac{\sqrt{c^2-4km}t}{m}\Big)+\\&\exp\Big(\frac{ct}{m}\Big)(cu^2+2muu^{\prime})(c^2\sqrt{c^2-4km}-(c^2-4km)^{3/2})\sinh\Big(\frac{\sqrt{c^2-4km}t}{m}\Big)\Big).
\end{aligned}
\end{equation}
Furthermore, the Lie algebra of these operators is shown in the following table
\begin{table}[h!]
	\centering
	\begin{tabular}{c| c c c c c}
		&$\textbf{X}_{1}$&$\textbf{X}_{2}$ & $\textbf{X}_{3}$ & $\textbf{X}_{4}$ & $\textbf{X}_{5}$ \\
		\hline
		\hline
		$\textbf{X}_{1}$ & $0$ & $\frac{\sqrt{c^2-4km}}{m}\textbf{X}_{3}$ & $\frac{\sqrt{c^2-4km}}{m}\textbf{X}_{2}$ & $-\frac{\sqrt{c^2-4km}}{2m}\textbf{X}_{4}$ &$\frac{\sqrt{c^2-4km}}{2m}\textbf{X}_{5}$ \\

		$\textbf{X}_{2}$ &$-\frac{\sqrt{c^2-4km}}{m}\textbf{X}_{3}$&$0$& $-\frac{\sqrt{c^2-4km}}{m}\textbf{X}_{1}$ & $-\frac{\sqrt{c^2-4km}}{2m}\textbf{X}_{5}$ &$-\frac{\sqrt{c^2-4km}}{2m}\textbf{X}_{4}$ \\

		$\textbf{X}_{3}$ &$-\frac{\sqrt{c^2-4km}}{m}\textbf{X}_{2}$ & $\frac{\sqrt{c^2-4km}}{m}\textbf{X}_{1}$ & $0$ &$-\frac{\sqrt{c^2-4km}}{2m}\textbf{X}_{5}$ & $\frac{\sqrt{c^2-4km}}{2m}\textbf{X}_{4}$\\

		$\textbf{X}_{4}$ & $\frac{\sqrt{c^2-4km}}{2m}\textbf{X}_{4}$ & $\frac{\sqrt{c^2-4km}}{2m}\textbf{X}_{5}$ & $\frac{\sqrt{c^2-4km}}{2m}\textbf{X}_{5}$ &$0$ &$0$ \\

		$\textbf{X}_{5}$ & $-\frac{\sqrt{c^2-4km}}{2m}\textbf{X}_{5}$ & $\frac{\sqrt{c^2-4km}}{2m}\textbf{X}_{4}$ & $-\frac{\sqrt{c^2-4km}}{2m}\textbf{X}_{4}$ &$0$ & $0$\\ [1ex]
	\end{tabular}
	\caption{Commutators Table}
	\label{table:1}
\end{table}\\
The list of generators provided in Eq. (12) except $X_1$ is entirely different from that presented in \cite{Cer}. Consequently, the associated first integrals are also non-identical to those mentioned in \cite{Cer}. Moreover, from above we see that if we put $c=0$, the operator $X_1$ becomes time translational symmetry which yields integral $I_1$ as energy conservation of the usual harmonic oscillator. Furthermore, under the assumption that the system is free of damping term, all the above symmetry generators and corresponding first integrals become identical to what determine by Lutzky \cite{Lut}.\\
\textbf{Case 2. Under-damped Oscillator:}\\ If the restoring force is large compared to the damping force term, then the system is termed as the under-damped. This situation occurs when $4km>c^2$, so the resulting Noether symmetry generators are expressed as
\begin{equation}\label{genset3}
\begin{aligned}
\textbf{X}_{1}&=\frac{\partial}{\partial t}-\frac{cu}{2m}\frac{\partial}{\partial u},\\
\textbf{X}_{2}&=\sin\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)\frac{\partial}{\partial t}+\Big(\frac{u\sqrt{4km-c^2}}{2m}\cos\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)-\frac{cu}{2m}\sin\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)\Big)\frac{\partial}{\partial u},\\
\textbf{X}_{3}&=\cos\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)\frac{\partial}{\partial t}-\Big(\frac{u\sqrt{4km-c^2}}{2m}\sin\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)+\frac{cu}{2m}\cos\Big(\frac{\sqrt{4km-c^2}t}{m}\Big)\Big)\frac{\partial}{\partial u},\\
\textbf{X}_{4}&=\exp\Big(\frac{-ct}{2m}\Big)\Big(\cos\Big(\frac{\sqrt{4km-c^2}t}{2m}\Big)-\iota \sin\Big(\frac{\sqrt{4km-c^2}t}{2m}\Big)\Big)\frac{\partial}{\partial u},\\
\textbf{X}_{5}&=\exp\Big(\frac{-ct}{2m}\Big)\Big(\cos\Big(\frac{\sqrt{4km-c^2}t}{2m}\Big)+\iota \sin\Big(\frac{\sqrt{4km-c^2}t}{2m}\Big)\Big)\frac{\partial}{\partial u}.
\end{aligned}
\end{equation}
In this case the operators $\textbf{X}_{4}$ and $\textbf{X}_{5}$, are different from that shown in (10). Interestingly, the splitting one of these two symmetries $\textbf{X}_{4}$ and $\textbf{X}_{5}$ into their real and imaginary parts results in two new Noether symmetries as shown
\begin{equation}
\begin{aligned}
\textbf{G}_{4}&=\exp\Big(\frac{-ct}{2m}\Big)\cos\Big(\frac{\sqrt{4km-c^2}t}{2m}\Big)\frac{\partial}{\partial u},\\
\textbf{G}_{5}&=\exp\Big(\frac{-ct}{2m}\Big)\sin\Big(\frac{\sqrt{4km-c^2}t}{2m}\Big)\frac{\partial}{\partial u}.
\end{aligned}
\end{equation}
Now insertion of the Lagrangian (9) along with Noether symmetries mentioned in Eq. (16), in Noether's theorem (7) leads towards following pair of conserved quantities
\begin{equation}
\begin{aligned}
I_{4}&=\exp\Big(\frac{ct}{2m}\Big)\Big[(mu^{\prime}+\frac{cu}{2})\cos\Big(\frac{\sqrt{4km-c^2}t}{2m}\Big)+\frac{1}{2}u\sqrt{4km-c^2}\sin\Big(\frac{\sqrt{4km-c^2}t}{2m}\Big)\Big],\\
I_{5}&=\exp\Big(\frac{ct}{2m}\Big)\Big[(mu^{\prime}+\frac{cu}{2})\sin\Big(\frac{\sqrt{4km-c^2}t}{2m}\Big)-\frac{1}{2}u\sqrt{4km-c^2}\cos\Big(\frac{\sqrt{4km-c^2}t}{2m}\Big)\Big].
\end{aligned}
\end{equation}
Interestingly, though forms of both pairs of generators $X_4, X_5$ given in Eq. (15) and $G_4, G_5$ in (16) are different, even then they provide the same first integrals. Moreover, the Lie algebra of the operators $G_4$ and $G_5$ is Abelian, i.e., $[G_4, G_5]=0$ and the elimination of $u'$ from the Eq. (17) yields an analytical solution for the equation of under-damped oscillator as
\begin{equation}
u(t)=\frac{1}{\sqrt{4km-c^2}}\exp\Big(-\frac{ct}{2m}\Big)\Big[C_1\sin\Big(\frac{\sqrt{4km-c^2}t}{2m}\Big)-C_2\cos\Big(\frac{\sqrt{4km-c^2}t}{2m}\Big)\Big].
\end{equation}
\\ The commutation relations in view of operators (14) are displayed in TABLE II
\\
\\
\begin{table}[h!]
	\centering
	\begin{tabular}{c| c c c c c}
		&$\textbf{X}_{1}$&$\textbf{X}_{2}$ & $\textbf{X}_{3}$ & $\textbf{X}_{4}$ & $\textbf{X}_{5}$ \\
		\hline
		\hline
		$\textbf{X}_{1}$ & $0$ & $\frac{\sqrt{4km-c^2}}{m}\textbf{X}_{3}$ & $-\frac{\sqrt{4km-c^2}}{m}\textbf{X}_{2}$ & $-\frac{\sqrt{4km-c^2}}{2m}\textbf{X}_{5}$ &$\frac{\sqrt{4km-c^2}}{2m}\textbf{X}_{4}$ \\

		$\textbf{X}_{2}$ &$-\frac{\sqrt{4km-c^2}}{m}\textbf{X}_{3}$&$0$& $-\frac{\sqrt{4km-c^2}}{m}\textbf{X}_{1}$ & $-\frac{\sqrt{4km-c^2}}{2m}\textbf{X}_{4}$ &$\frac{\sqrt{4km-c^2}}{2m}\textbf{X}_{5}$ \\

		$\textbf{X}_{3}$ &$\frac{\sqrt{4km-c^2}}{m}\textbf{X}_{2}$ & $\frac{\sqrt{4km-c^2}}{m}\textbf{X}_{1}$ & $0$ &$\frac{\sqrt{4km-c^2}}{2m}\textbf{X}_{5}$ & $\frac{\sqrt{4km-c^2}}{2m}\textbf{X}_{4}$\\

		$\textbf{X}_{4}$ & $\frac{\sqrt{4km-c^2}}{2m}\textbf{X}_{5}$ & $\frac{\sqrt{4km-c^2}}{2m}\textbf{X}_{4}$ & $-\frac{\sqrt{4km-c^2}}{2m}\textbf{X}_{5}$ &$0$ &$0$ \\

		$\textbf{X}_{5}$ & $-\frac{\sqrt{4km-c^2}}{2m}\textbf{X}_{4}$ & $-\frac{\sqrt{4km-c^2}}{2m}\textbf{X}_{5}$ & $-\frac{\sqrt{4km-c^2}}{2m}\textbf{X}_{4}$ &$0$ & $0$\\ [1ex]
	\end{tabular}
	\caption{Commutators Table}
	\label{table:2}
\end{table}\\
\textbf{Case 3. Critical-damped oscillator:}\\ If the restoring force and damping force are comparable in effect then the dynamical system behaves as critical damping oscillator. In this situation, we have $4km-c^2=0$ and the Eq. (5) provides the following 5-Noether symmetries
\begin{equation}
\begin{aligned}
\textbf{X}_{1}&=\exp\left(-\sqrt{\frac{k}{m}}t\right)\frac{\partial}{\partial u},\\
\textbf{X}_{2}&=t\exp\left(-\sqrt{\frac{k}{m}}t\right)\frac{\partial}{\partial u},\\
\textbf{X}_{3}&=\frac{\partial}{\partial t}-\sqrt{\frac{k}{m}}u\frac{\partial}{\partial u},\\
\textbf{X}_{4}&=t\frac{\partial}{\partial t}-\frac{1}{2}\Big(2t\sqrt{\frac{k}{m}}-1\Big)u\frac{\partial}{\partial u},\\
\textbf{X}_{5}&=\frac{1}{2}t^2\frac{\partial}{\partial t}-\frac{1}{2}\left(t^2\sqrt{\frac{k}{m}}-t\right)u\frac{\partial}{\partial u}.
\end{aligned}
\end{equation}
We mention here that the forms of Noether symmetry generators in the critical damped harmonic case are entirely different from what we have seen for over-damped and under-damped cases.\\ The Lie algebra of these operators is shown below
\begin{table}[h!]
	\centering
	\begin{tabular}{c| c c c c c}
		&$\textbf{X}_{1}$&$\textbf{X}_{2}$ & $\textbf{X}_{3}$ & $\textbf{X}_{4}$ & $\textbf{X}_{5}$ \\
		\hline
		\hline
		$\textbf{X}_{1}$ & $0$ & $0$ & $0$ & $\frac{1}{2}\textbf{X}_{1}$ &$\frac{1}{2}\textbf{X}_{2}$ \\
		$\textbf{X}_{2}$ &$0$&$0$& $-\textbf{X}_{1}$ & $-\frac{1}{2}\textbf{X}_{2}$ &$0$ \\
		$\textbf{X}_{3}$ &$0$ & $\textbf{X}_{1}$ & $0$ &$\textbf{X}_{3}$ & $\textbf{X}_{4}$\\
		$\textbf{X}_{4}$ & $-\frac{1}{2}\textbf{X}_{1}$ & $\frac{1}{2}\textbf{X}_{2}$ & $-\textbf{X}_{3}$ &$0$ &$\textbf{X}_{5}$ \\
		$\textbf{X}_{5}$ & $-\frac{1}{2}\textbf{X}_{2}$ & $0$ & $-\textbf{X}_{4}$ &$-\textbf{X}_{5}$ & $0$\\ [1ex]
	\end{tabular}
	\caption{Commutators Table}
	\label{table:3}
\end{table}\\
Now considering these Noether symmetries and employing the classical Noether's theorem provide respectively, the five first integrals
\begin{equation}
\begin{aligned}
I_{1}&=\exp\Big(\sqrt{\frac{k}{m}}t\Big)\Big(\sqrt{\frac{k}{m}}u+u^{\prime}\Big),\\
I_{2}&=\exp\Big(\sqrt{\frac{k}{m}}t\Big)\Big(\sqrt{\frac{k}{m}}tu+tu^{\prime}-u\Big),\\
I_{3}&=\Big(\exp\Big(\sqrt{\frac{k}{m}}t\Big)\Big)^{2}\Big(\frac{ku^{2}}{2m}+\sqrt{\frac{k}{m}} uu^{\prime}+\frac{1}{2}u^{\prime 2}\Big),\\
I_{4}&=\Big(\exp\Big(\sqrt{\frac{k}{m}}t\Big)\Big)^{2}\Big(\frac{t}{2}\Big(\sqrt{\frac{k}{m}}u+u^{\prime}\Big)^{2}-\Big(\sqrt{\frac{k}{m}}tu-\frac{u}{2}+tu^{\prime}\Big)\Big(\sqrt{\frac{k}{m}}u+u^{\prime}\Big)\Big),\\
I_{5}&=\frac{u^{2}}{4}\exp\Big(2t\sqrt{\frac{k}{m}}\Big)+\frac{t}{4m}(2mu(tu^{\prime}-u)\sqrt{\frac{k}{m}}+tmu^{\prime 2}+tku^{2}-2muu^{\prime})\Big(\exp\Big(\sqrt{\frac{k}{m}}t\Big)\Big)^{2}.
\end{aligned}
\end{equation}
We stress here that the forms of Noether symmetry generators $X_1, X_2, X_3, X_4$ given in Eq. (20) are similar to that found in \cite{Cer} and the first integrals as well. However, the form of operator $X_5$ is different and the resulting first integral is also new and fulfill the condition to becomes first integral of the critical damped harmonic oscillator. Furthermore, from the invariant quantities $I_1, I_2$ given in Eq. (20), we obtain following solution for the critical damped harmonic oscillator
\begin{equation}
u(x)=e^{-\sqrt{k/m}t}\Big[I_1 t-I_2\Big].
\end{equation}
Interestingly, we point out here that the insertion of the following Lagragian in Euler-Lagrange equation (3) also yields the same damped harmonic oscillator equation (8)
\begin{equation}
L=\frac{e^{{\frac{c}{m}}t}}{4m^2}\Big[2m^2u'^2+2cmuu'+(c^2-2km)u^2\Big].
\end{equation}
Remarkably, if we employ this new form of Lagragian for the DHO in Noether symmetry condition and classical Noether's theorem, we obtain same set of Noether symmetries and constants of motion as we have determined for over-, under- and critical-damped oscillators. So this Lgargian can also be considered as the standard Lagragian. Though the new form of Lagrangian is completely different yet it yields same results as we derived using the Bateman Lagrangian. So this new form of optimal Lagragian function given in (22) may conceal something interesting for the physicists. We conclude here that this new form of Lagrangian may lead to entirely different quantum mechanical properties related to damped harmonic oscillator. Moreover, some interesting results related to DHO employing Eq. (22) will be discussed in future work.
\section{Conclusion}
In this paper we have classified the harmonic oscillator problem with respect to Noether symmetries by constructing related conserved quantities. We have discussed three cases of damped harmonic oscillator namely, over-damped, under-damped and critical-damped. We observe that when we impose over-damped, under-damped and critical-damped conditions, the forms of few generators given in Eq. (10) get changed and provide new first integrals compared to that mentioned in Eq (11). Being an old and classical problem, several authors have attempted to study the group theoretic properties of variational DHO equation. In the current study we have also investigated the DHO equation of motion from group theoretic view point. In this paper we have determined some new results along with old one which have already been reported in the literature. For instance, in the over-damped harmonic oscillator case, the forms of all generators except $X_1$ are entirely different to what mentioned in \cite{Cer} and consequently the corresponding conserved quantities are also different. Similarly, for the under-damped harmonic oscillator case the real and imaginary parts of operator $X_4$ provide two Noether symmetry generators. These resulting operators form an Abelian algebra and yield two functionally independent first integrals. Furthermore, these two integrals suffice to provide the exact solution of under-damped oscillator. In the critical case the four Noether symmetry generators are similar to that mentioned in \cite{Cer}, while the operator $X_5$ is new and distinct from the list provided in \cite{Cer} and the resulting first integral also fulfills the condition to become the first integral of critical damped harmonic oscillator. In all three cases of damped harmonic oscillator the algebra of resulting Noether symmetry generators is closed. \\
We have also presented a new form of a Lagrangian function for the DHO system. We observe that for DHO on using new form of Lagrangian provides same set of Noether symmetries and associated conservation Laws as previously discussed for over-, under- and critical-damped oscillators. So this Lagrangian can also be considered as the standard Lagrangian. Though the new form of Lagrangian is completely different to that Bateman Lagrangian but it yields same results as we derive employing Bateman Lagrangian. So this new form of optimal Lagrangian function given in (22) may be usefully applied to explore some mechanical properties related to DHO.


\section{Introduction}

Entanglement is the clearest nonclassical signature of quantum physics. A composite system is considered to be entangled when the quantum state that describes it is inseparable, \textit{i.e.},  it is impossible to write it as a product of the states of each subsystem \cite{horodecki2009}. In the last decades, entangled states have been the subject of great interest, presenting themselves as a resource for several quantum schemes and applications such as quantum communication \cite{ekert1991,bennett1993,bennett2000}, quantum computation \cite{gottesman1999}, metrology \cite{giovannetti2006} and quantum sensing \cite{degen2017}.

The success of the aforementioned applications and tests often depends on the ability to generate long-lived entangled states. However, in a realistic situation, the system will interact with the environment. This will inevitably lead to the deterioration of the entangled state, which is sensitive to decoherence \cite{horodecki2009}. For this reason, entanglement preservation schemes have gained great prominence. Among the proposed methods to minimize unwanted decoherence, we should mention the use of decoherence-free subspaces \cite{lidar1998,viola1998}, quantum error correction codes \cite{steane1996,knill1997}, weak measurements \cite{sun2010,kim2012} and the quantum Zeno effect \cite{facchi2004,maniscalco2008}.

Instead of aiming to prevent decoherence, a different strategy involves engineering the system-environment interaction to \textit{generate} entangled states, these are called dissipation-assisted entanglement generation methods \cite{diehl2008,kraus2008,verstraete2009}. Since this idea was presented, numerous implementations have been proposed and experimentally realized using several physical platforms, such as cavity QED \cite{plenio1999,nicolosi2004,wang2010,kastoryano2011,del2011}, superconducting qubits \cite{shankar2013,kimchi2016,liu2016,doucet2020}, macroscopic atomic ensembles \cite{krauter2011,santos2022}, Rydberg atoms \cite{rao2013,carr2013,shao2017} and trapped ions \cite{barreiro2011,lin2013,bentley2014,horn2018,cole2021,cole2022}. Another widely-studied approach is the creation of long-lived entangled states via stimulated Raman adiabatic passage (STIRAP) \cite{vitanov2017,bergmann2019}, since it offers robustness against decoherence by not populating lossy states.

In this paper we propose two new schemes for producing highly-entangled states in a system of two non-degenerate qubits. It is known that, in some systems, it's possible to manipulate the degree of entanglement between two qubits via Stark shifts of their electronic levels \cite{hettich2002,wang2010,trebbia2021}. Here, however, we show that the symmetry in the energy shifts between the emitters with respect to the bosonic mode actually allows us to achieve a stronger entanglement. By considering their effective interaction through a bosonic mode, such as an optical cavity or a motional mode, we are able to achieve a stronger coupling between the qubits without needing to place the qubits particularly close together. This interaction with the quantized mode provides a coupling regime strong enough so that the time scales of the effective interactions are much faster than the qubit relaxation, leading to higher degrees of entanglement. This shows that the distinguishability between the quantum emitters can be an advantage in the quest for the producing highly-entangled states. We are able to achieve a maximally-entangled steady state, which is maintained by injecting power via a pump field.

We show that the dynamics leading to the highly-entangled two-qubit state can be understood by comparison with electromagnetically induced transparency (EIT) regime \cite{fleischhauer2005}. In that same analogy, but restricting ourselves to the subspace of just a single excitation, we are able to drive the system into an entangled state via a STIRAP-like process.

\section{Model}

Let us consider a system of two qubits, with different resonance frequencies $\omega_e^{(1)}$ and $\omega_e^{(2)}$, which are coupled to the same bosonic mode with frequency $\omega_m$, as illustrated in Fig.~\ref{fig:1}(a). Here, the bosonic mode is symmetrically detuned from each of the  qubits, so that $\omega_e^{(1)}=\omega_m-\Delta$ and $\omega_e^{(2)}=\omega_m+\Delta$. In the Schr\"odinger picture, the system Hamiltonian reads ($\hbar=1$):
\begin{eqnarray}
\hat{H} = &\omega_{e}^{(1)}\hat{\sigma}_{ee}^{(1)}+\omega_{e}^{(2)}\hat{\sigma}_{ee}^{(2)}+\omega_{m}\hat{a}^{\dagger}\hat{a}\nonumber\\
&+g\left[\hat{a}\left(\hat{\sigma}_{+}^{(1)}+\hat{\sigma}_{+}^{(2)}\right)+h.c.\right],\label{eq:1}
\end{eqnarray}
where $\hat{\sigma}_{+}^{(k)}=\vert e \rangle\langle g\vert$, $\hat{\sigma}_{-}^{(k)}=\vert g\rangle\langle e\vert$ and $\hat{\sigma}_{ee}^{(k)}=\vert e\rangle\langle e\vert$ are the raising, lowering and excited-state population operators, respectively, acting on the $k$-th qubit (with $k\in\{1,2\}$). Without loss of generality, we define the ground state energy to be zero. $\hat{a}$ ($\hat{a}^\dagger$) is the annihilation (creation) operator of the bosonic mode, $g$ is the coupling strength between the bosonic mode and each of the qubits, and $h.c.$ stands for the Hermitian conjugate.

\begin{figure}[ht]
\begin{centering}
\includegraphics[width=0.9\columnwidth]{FIG1_final3.png}
\par\end{centering}
\caption{a)~Two non-degenerate qubits, with both ground states coupled to the bosonic mode of frequency $\omega_{m}$, detuned by $\pm \Delta$. b) Level scheme of the same system, but now in the basis up to one excitation:  $\vert G\rangle \otimes \vert 0\rangle$, $\vert E\rangle=\vert G\rangle \otimes \vert 1\rangle$, and $\vert\Psi_{\pm}\rangle=\vert\Phi_{\pm}\rangle\otimes\vert 0 \rangle$, where $\vert G\rangle=\vert g,g \rangle$ and $\vert \Phi_{\pm}\rangle = \left(\vert e,g\rangle\pm\vert g,e\rangle\right)/\sqrt{2}$. Here, $g$ promotes transitions from $\vert E\rangle$ to $\vert\Psi_{+}\rangle$, while $\Delta$ from $\vert\Psi_{+}\rangle$ to $\vert\Psi_{-}\rangle$. We consider a pump field, of strength $\varepsilon$, continuously injecting energy into the mode to combat decay from the mode with rate $\kappa$. The inset in b) shows the entangled steady state partial density matrix, where the bosonic mode has been traced out.}
\label{fig:1}
\end{figure}

For convenience, we move to the interaction picture, make the rotating wave approximation and move to a rotating referential of relative coordinates in which both qubits are stationary, thus eliminating the Hamiltonian time dependence. Then, Eq. \eqref{eq:1} becomes
\begin{equation}
\mathcal{\hat{H}}=\Delta\left(\hat{\sigma}_{ee}^{(1)}-\hat{\sigma}_{ee}^{(2)}\right)+g\left[\hat{a}\left(\hat{\sigma}_{+}^{(1)}+\sigma_{+}^{(2)}\right)+h.c.\right].\label{eq:3}
\end{equation}

To account for decoherence, we consider our system to be a weak system-environment coupling regime, which allows us to use the Lindblad master equation \cite{petruccione2002} at temperature $T=0\,$K. The assumption of zero temperature is reasonable since we work within the optical regime, where the number of thermal photons remain negligible even for room temperatures. Thus, we obtain the dynamical equations for the density matrix $\rho$
\begin{equation}
\dot{\hat{\rho}}=-i\left[\hat{\mathcal{H}},\hat{\rho}\right]+\hat{\mathscr{L}}_{\mathrm{q}}^{(1)}+\hat{\mathscr{L}}_{\mathrm{q}}^{(2)}+\hat{\mathscr{L}}_{m},\label{eq:5}
\end{equation}
where
\begin{equation}
\hat{\mathscr{L}}_{\mathrm{q}}^{(k)}=\Gamma\left(2\hat{\sigma}_{-}^{(k)}\hat{\rho}\hat{\sigma}_{+}^{(k)}-\hat{\sigma}_{ee}^{(k)}\hat{\rho}-\hat{\rho}\hat{\sigma}_{ee}^{(k)}\right)
\end{equation}
is the Lindblad term that accounts for the spontaneous decay from the excited state of the $k$-th qubit, with $k\in\{1,2\}$ and $\Gamma$ is the qubit decay rate, here assumed the same for both qubits. The Lindblad term
\begin{equation}
\hat{\mathscr{L}}_{m}=\kappa\left(2\hat{a}\hat{\rho} \hat{a}^{\dagger}-\hat{a}^{\dagger}\hat{a}\hat{\rho}-\hat{\rho} \hat{a}^{\dagger}\hat{a}\right)
\end{equation}
accounts for decay of the bosonic mode, where $\kappa$ is the decay rate.

\section{Steady state entanglement production}\label{sec2}
To describe the main mechanism responsible for the generation of entanglement, we restrict ourselves, for the moment, to the single-excitation subspace, which is composed of the following three states: $\vert E\rangle=\vert G\rangle \otimes \vert 1\rangle$ and $\vert\Psi_{\pm}\rangle=\vert\Phi_{\pm}\rangle\otimes\vert 0 \rangle$, where $\vert G\rangle=\vert g,g \rangle$ and $\vert \Phi_{\pm}\rangle = \left(\vert e,g\rangle\pm\vert g,e\rangle\right)/\sqrt{2}$. $\vert\Phi_{\pm}\rangle$ are maximally-entangled two-qubit states. In this subspace, the reduced Hamiltonian is:
\begin{equation}
    \hat{\bar{\mathcal{H}}}=\Delta\vert\Psi_{-}\rangle\langle\Psi_{+}\vert +\sqrt{2}g\vert\Psi_{+}\rangle\langle E\vert +h.c.
    \label{eq:4}
\end{equation}
In analogy with the typical three-level $\Lambda$ systems, we consider that $\vert E\rangle$ and $\vert \Psi_-\rangle$ play the roles of the two ground states, while $\vert \Psi_+\rangle$ is the excited state, as depicted in Fig.~\ref{fig:1}(b). According to Eq.~\eqref{eq:4}, the transitions $\vert E\rangle \leftrightarrow \vert \Psi_+\rangle$ and $\vert \Psi_-\rangle \leftrightarrow \vert \Psi_+\rangle$ have effective coupling strengths $\sqrt{2}g$ and $\Delta$, respectively. The Hamiltonian in Eq.~\eqref{eq:4} has a \textit{dark eigenstate} (an eigenstate without a $\vert \Psi_{+}\rangle$ component) given by
\begin{equation}
    \vert D\rangle = -\frac{\Delta}{\sqrt{\Delta^2+2g^2}}\vert E\rangle+\frac{\sqrt{2}g}{\sqrt{\Delta^2+2g^2}}\vert\Psi_{-}\rangle. \label{eq:darkstate}
\end{equation}
Just as in other three-level systems \cite{fleischhauer2005}, when the condition $g\gg\Delta$ is fulfilled, the dark state $\vert D\rangle$ transforms to $\vert\Psi_{-}\rangle$, which is our maximally-entangled target state.

A significant difference between our system and a typical three-level $\Lambda$-system is the fact that neither $\vert E \rangle$ or $\vert\Psi_-\rangle$ are actually ground states. Because of this, the system keeps spontaneously decaying to the true zero energy ground state $\vert G\rangle \otimes \vert 0\rangle$. For this reason, the analogy becomes more accurate as the mode dissipation rate $\kappa$ and the spontaneous decay rate $\Gamma$ become negligible ($\kappa,\Gamma \ll g,\Delta$). Then the system effectively remains in the one-excitation subspace. One way to circumvent the decay and maintain the system in the one-excitation subspace is to keep injecting energy into the mode. We consider this energy injection as an additional term
\begin{equation}
    \hat{H}_{\mathrm{pump}}=\varepsilon\left(\hat{a}^{\dagger}+\hat{a}\right)
\end{equation}
of the system's Hamiltonian, where $\varepsilon$ is the the pump strength. In our simulations, $\varepsilon$ was kept constant at $\varepsilon=\kappa$.

To characterize the steady-state entanglement between the qubits, we choose the monotone quantifier concurrence \cite{Wootters1997,Wootters1998}. We numerically simulate the full system dynamics, taking into account the decoherence and also higher excited states. The concurrence is derived from the steady-state density matrix, which we obtain using the Quantum Toolbox in Python (QuTiP) \cite{qutip2013}, after the mode is traced out.
\begin{figure}[ht]
\begin{centering}
\includegraphics[width=1\columnwidth]{Fig2_matplotlib2.png}
\par\end{centering}
\caption{Entanglement generation using the steady-state method. (a) Colormap showing the concurrence of the steady state as a function of $\Delta/\kappa$ and $g/\kappa$. Points A to D are referred to in panel (c). (b) Concurrence as a function of the spontaneous decay rate $\Gamma$ of each of the two qubits (top x-axis) and as a function of the detuning (bottom x-axis). (c) Concurrence as a function of time for different parameter sets, given in panel (a). The concurrence generally grows and then stabilises. The insets in panel (c) show the steady state partial density matrix of each curve, where the mode has been traced out. In all panels, the pump strength was set to $\varepsilon=\kappa$, and for (a) and (c) the spontaneous decay rate of each qubit was $\Gamma=10^{-5}\kappa$.}
\label{fig:2}
\end{figure}

As shown in Fig.~\ref{fig:2}(a), we obtain a strong entanglement within a large region of parameters $\Delta$ and $g$. Moreover, we see that satisfying the EIT condition $g\gg \Delta$ is a necessary but not sufficient condition to reach our maximally entangled target state. We have chosen four sets of parameters to illustrate the system dynamics: $A$ shows the situation where $g \ll \Delta$ and the population is lead to the respective dark state, which has a $\vert G\rangle$ character [see Eq.~\eqref{eq:darkstate}], and no entanglement is observed; in $B$, we have $g = \Delta$ which leads to weakly entangled dark state given by a mixture of $\vert g,g\rangle$ and $\vert \Phi_- \rangle$; in $C$ we have an interesting situation where the EIT condition $g \gg \Delta$ is fulfilled, the dark state a $\vert\Phi_- \rangle$ character, but the detuning is so small that it takes too long to populate $\vert \Phi_- \rangle$, and since the system keeps decaying to the true ground state $\vert G \rangle \otimes \vert 0\rangle$, the entanglement is affected; $D$ shows a near-optimal situation, with $g \gg \Delta$ and $\Delta \gg \Gamma$, driving approximately all the population to the dark state $\vert \Phi_- \rangle$.

In Fig.~\ref{fig:2}(b) we show the behavior of the concurrence as a function of the detuning for a fixed value of $g=\kappa$. We observe that the concurrence decreases significantly when we move away from the optimal point $\Delta\approx0.1g$. This finding is in accordance with what was previously discussed for the parameters sets $A$ and $C$. This optimal ratio between $g$ and $\Delta$ is influenced by the system's decay rate $\Gamma$. A reduced value of $\Gamma$ allows us to achieve entanglement for smaller detunings. For $\Gamma=0$, even an infinitesimal detuning would eventually take the system to $\vert \Phi_- \rangle$. Fig.~\ref{fig:2}(b) shows that the spontaneous decay reduce the entanglement even for the optimal case $D$; the concurrence decreases exponentially for $\Gamma>5\times10^{-4}\kappa$.

Focusing on the entanglement dynamics, we show in Fig.~\ref{fig:2}(c) the concurrence over time for the cases $A$, $B$, $C$ and $D$, where we observe different time scales to achieve the maximum entanglement for each curve, respectively. We also include a visual representation of the time evolved partial density matrix, where the bosonic mode has been traced out, showing that we achieve the maximally entangled target state $\vert\Phi_-\rangle$ for the near-optimal set of parameters $D$.

\section{Adiabatic process}\label{sec:stirap}

Inspired by the schemes to counteract decoherence developed for multilevel atoms, we propose a STIRAP-like process to efficiently populate the entangled state $\vert\Psi_{-}\rangle$. To this end, we consider a trapped ion implementation, where two ions are confined in a harmonic potential and coupled to the same motional mode, as depicted in Fig.~\ref{fig:3alt}(a). The ions are subjected to a magnetic field gradient, thus will experience different energy shifts of their excited states, resembling the system illustrated in Fig.~\ref{fig:1}(a). To perform the adiabatic population transfer, the system must be initially prepared in the state $\vert E\rangle$, which consists of the two ions in the ground state and the motional state with one excitation; the latter can be prepared by exciting one of the ions and then by letting it exchange energy with the vibration mode via a red sideband interaction. The initial state $\vert E\rangle$ corresponds to the dark state when $\Delta\gg g$. By reversing this condition to $g\gg\Delta$ the dark state adiabatically transforms to $\vert\Psi_-\rangle$. By manipulating the coupling strength $g$ and the detuning $\Delta$, which can be accomplished by changing the power of a laser resonant to the red sideband transition and by magnetic field gradient varying in time, respectively, we can control the adiabaticity throughout the process so that all the population is transferred coherently to the final maximally entangled state $\vert\Psi_{-}\rangle$, with a negligible population in the state $\vert\Psi_{+}\rangle$ (see Appendix A).

\begin{figure}[htb]
\begin{centering}
\includegraphics[width=1\columnwidth]{STIRAP_g_const.png}
\par\end{centering}
\caption{Entanglement generation - adiabatic method. (a) Two two-level ions are trapped in a harmonic potential and coupled to the same phonon mode. The ions are subjected to a magnetic field gradient, which promotes different energy shifts to their excited states, recovering the system illustrated in Fig.~\ref{fig:1}(a). (b) Pictorial representation of the energy shifts in the ions excited states due to a magnetic field gradient, as a function of their position in the trap. (c) Time evolution of $g$ and $\Delta$, with $\Delta$ varying accordingly to Eq.~\eqref{eq:D} and a constant $g$. (d) Population changes over time. We considered $\Gamma=10^{-3}\kappa$, $\Delta_{\mathrm{max}}=2\times10^{5}\kappa$, $g=2.5\times10^{2}\kappa$, $t_0=4\times10^{-3}/\kappa$, $\lambda=5\times10^2\kappa$.}
\label{fig:3alt}
\end{figure}

This scheme can also be applied to other experimental platforms. However, these might not allow a perfect control over both parameters at the same time, e.g., a fixed $g$ and a tunable $\Delta$, as often observed in quantum circuits. To show that even in this situation we can achieve a high degree of entanglement, we assume $g$ to be constant and we vary $\Delta$ in time, according to
\begin{equation}
\Delta\left(t\right) = \frac{\Delta_{\mathrm{max}}}{2}\left[1-\tanh{\lambda(t-t_0)}\right], \label{eq:D}
\end{equation}
as shown in Fig.~\ref{fig:3alt}(c). Here, $\Delta_{\mathrm{max}}$ corresponds to the maximum value of the detuning and $t_0$ is the time at which this function reaches its respective half maximum value $\Delta_{\mathrm{max}}/2$. The adiabaticity of the process is controlled by $\lambda$, which determines the time scale of the parameter swap and, consequently, how fast the parameters are changed. The specific choice of function is not crucial to the method, as long as it guarantees an near-adiabatic inversion of parameters. Setting the value of $g$ that ensures the initial and final conditions $\Delta \gg g$ and $g \gg \Delta$, respectively, we make the parameter swapping as smooth as possible. On the other hand, the lack of control over $g$ restricts the time scale of the population transfer to maintain the adiabaticity during the parameter swap. Nevertheless, we show that one can still perform a STIRAP-like process, obtaining a highly entangled final state, as shown in \ref{fig:3alt}(d), with a negligible population of the state $\vert \Psi_+\rangle$. In the case where both $g$ and $\Delta$ can be simultaneously controlled, a near-perfect STIRAP process can be achieved (see Appendix A).

\section{Conclusion}

In conclusion, here we present two novel strategies for producing maximally entangled states in a system of two non-degenerate qubits coupled to a single bosonic mode. In both cases, we use a direct analogy with ordinary three-level atomic systems in the $\Lambda$-configuration and the EIT phenomenon, which allows us to draw a parallel with the processes of optical pumping and adiabatic population transfer.

In the first proposed scheme, we show, that it is possible to generate steady state entanglement with concurrence $C>0.99$. Moreover, we show that the symmetry between the qubits with respect to the bosonic mode is beneficial for the generation of entanglement. As for the second scheme, we generate a highly entangled state by means of an adiabatic process, and we achieve a population over $99\%$ in the state $\vert \Psi_-\rangle$. We emphasize that this was done even considering non-ideal situations, where there is no complete control over $g$ and $\Delta$. In ideal cases, where both parameters can be controlled simultaneously, the adiabatic process can be controlled more efficiently, which leads to a perfectly coherent transfer of populations, as shown in detail in Appendix A.

The results presented in this work, besides indicating new ways to generate highly entangled states in a simple way, have potential application in several experimental platforms, such as trapped ions and quantum dots molecules coupled to a cavity mode.

\acknowledgments
C.J.V.-B., R.B. and M.H.O. thank the support from the National Council for Scientific and Technological Development (CNPq) grants 311612/2021-0, 141247/2018-5, 409946/2018-4 and 313886/2020-2, and from the S\~ao Paulo Research Foundation (FAPESP) through Grants No.~2020/00725-9, 2019/13143-0, 2019/11999-5 and 2018/15554-5. The authors also thank the Joint Brazilian-Swedish Research Collaboration (CAPES-STINT), grant 88887.304806/2018-00 and BR2018-8054. M.H., G.H., C.Z. thank the support by the Knut \& Alice Wallenberg Foundation (through the Wallenberg Centre for Quantum Technology [WACQT]), by Olle Enkvists Stiftelse (Project No 204-0283), by Carl Trygger Stiftelse Grant No. CTS 19-139, and by the Swedish Research Council (Trapped Rydberg Ion Quantum Simulator).

\bibliographystyle{unsrt}