\section{introduction}

	The advent of intelligent reflecting surface (IRS) reshapes the wireless propagation environment between transceivers by tuning the reflecting elements, which is able to improve not only the efficiency of wireless energy transmission (WET) but also that of wireless information transmission (WIT) for future Internet-of-Things (IoT) applications \cite{clerckx2021wireless}. This thus has motivated recent research on IRS-aided wireless powered communication networks (WPCNs) and other related applications \cite{lyubin2021IRS,zheng2020joint,wu2021irs,bai2021resource,pan2020intelligent}.
	In the meanwhile, to further improve the spectral efficiency (SE) as well as user (UE) fairness, non-orthogonal multiple access (NOMA) has been recently employed for uplink (UL) WIT in WPCNs, which allows multiple users to utilize the same spectrum simultaneously by adopting successive interference cancellation (SIC) at the receiver \cite{diamantoulakis2016wireless}. A recent related work proved that dynamic IRS beamforming (adjusting the IRS phase-shift vector over downlink (DL) and UL) was not needed for the  IRS-aided WPCN with a hybrid access point (HAP), which helps reduce the number of IRS phase shifts to be optimized \cite{wu2021irs}.

	\begin{figure}
		\setlength{\abovecaptionskip} {0.cm}
		\setlength{\belowcaptionskip} {-1cm}
		\centering
		\includegraphics[width=0.45\textwidth]{figure1.eps}
		\caption{An IRS-aided WPCN.}
		\label{fig:Channel}
	\end{figure}

		Nevertheless, most prior works on WPCNs focused on maximizing the sum throughput \cite{lyubin2021IRS,zheng2020joint,wu2021irs}, which neglect the fairness and quality-of-service (QoS) requirements among users. Moreover, they assumed an oversimplified linear energy harvesting (EH) model, which may cause severe resource allocation mismatches and lead to significant performance degradation for practical systems with highly non-linear EH circuits \cite{valenta2014harvesting}. To our best knowledge, the QoS guaranteed resource allocation and beamforming design in the IRS-aided wireless powered NOMA with non-linear EH model have not been investigated yet, which may yield new insights and conclusions different from those in existing works. In particular, it remains unclear whether the dynamic IRS beamforming can improve the system performance or not for power-efficient WPCNs.

	Motivated by the above considerations, we consider an IRS-aided WPCN where an IRS is deployed to assist the DL WET and UL WIT.
	By employing the separated power station (PS) and data sink (DS) as shown in Fig. \ref{fig:Channel}, the “doubly near-far phenomenon” in WPCNs with a HAP can be avoided \cite{ju2013throughput}. Specifically, our objective is to minimize the system transmit energy consumption at the PS via joint optimization of resource allocation and IRS phase shifts. In particular, minimum throughput requirements are enforced for all the devices and the practical non-linear EH model is incorporated. By decomposing the coupled joint optimization into two subproblems, we solve the energy minimization problem by applying the alternating optimization (AO) technique and propose efficient algorithms for each subproblem respectively. Specifically, we first transform the subproblem of resource allocation into the difference of convex functions (DC) programming and solve it using the successive convex approximation (SCA) technique. On the other hand, for the IRS beamforming optimization, a penalty based algorithm is invoked.
	Numerical results demonstrate the superiority of the proposed algorithms compared to benchmark schemes and reveal the significant benefits of integrating IRS into WPCNs. Moreover, it is found that: 1) More transmit energy consumption can be saved by adopting dynamic IRS reconfigurations over DL and UL compared to the static IRS beamforming, which is different from the result for the throughput maximization problem in \cite{wu2021irs}; 2) Increasing the number of low-cost IRS reflecting elements instead of raising the emitted power at the PS is a more energy-efficient approach for achieving sustainable and green WPCNs.

	\emph{Notations:} Throughout the paper, superscripts $(\cdot)^{T}$, $(\cdot)^{H}$ and \text{diag}$(\cdot)$ represent the transpose,  Hermitian transpose and diagonalization operator, respectively. $\mathbb{C}^{a \times b}$ denotes the space of $a \times b$ complex matrices. $|\cdot|$ denotes the absolute value of a complex scalar.  $\|\mathbf{X}\|$, $\text{Tr}(\mathbf{X})$ and \text{Rank}$(\mathbf{X})$ denote the Euclidean norm, trace and rank of the matrix $\mathbf{X}$, respectively. ${\bf 1}_{N}$ denotes an all-one $N \times 1$ vector.

	\vspace{-0.25cm}
	\section{System Model and Problem Formulation}

	\vspace{-0.08cm}
	\subsection{Channel Model}
	\vspace{-0.08cm}
	As shown in Fig. \ref{fig:Channel}, we consider an IRS-aided WPCN, which consists of one PS, one IRS, one DS, and $K$ IoT devices. We assume that the PS, DS, and all the devices are equipped with a single antenna, while the IRS has $N$ reflecting elements. The total available transmission time is denoted by $T_{\max}$.
	To characterize the maximum achievable performance, it is assumed that the channel state information (CSI) of all channels is perfectly known at the PS and the DS (see \cite{wu2019towards} and its follow-up works for practical channel acquisition methods).

	Denote $\bs{g}_\text{E} \in \mathbb{C}^{N \times 1}$, $\bs{h}_{r,k} \in \mathbb{C}^{N \times 1}$ and $h_{d,k} \in \mathbb{C}$, $\forall k \in \mathcal{K}\triangleq\{1,\cdots,K\}$, as the equivalent baseband channels from the PS to the IRS, from the IRS to device $k$, and from the PS to device $k$, respectively. Likewise, the equivalent baseband channels from the IRS to the DS, from device $k$ to the IRS, and from device $k$ to the DS are denoted by $\bs{g}_\text{I} \in \mathbb{C}^{N \times 1}$, $\bs{q}_{r,k} \in \mathbb{C}^{N \times 1}$ and $q_{d,k} \in \mathbb{C}$, $\forall k \in \mathcal{K} $, respectively.
	The reflection phase-shift matrix for DL WET and UL WIT are denoted by $\bm{\Theta}_0 = \text{diag}(e^{j\phi_1},\cdots,e^{j\phi_N})$ and $\bm{\Theta}_1 = \text{diag}(e^{j\varphi_1},\cdots,e^{j\varphi_N})$, respectively, where $\phi_n$ and $\varphi_n  \in [0,2\pi), \forall n \in \mathcal{N}\triangleq\{1,\cdots,N\}$.
	%
	We assume that the WPCN adopts the typical ``harvest and then transmit'' protocol, where each device first collects energy from the radio frequency (RF) signal  broadcasted by the PS with constant transmit power $P_\text{E}$ for a time duration of $\tau_0$ in DL WET. Then, all the devices transmit their information signals to the DS simultaneously for a duration of $\tau_1$ with transmit power $p_k$, $\forall k \in \mathcal{K}$, in UL WIT. The RF power received at device $k$ is then given by
	\begin{sequation}
		\begin{aligned}
			P_{{\text{EH}}_k} &= P_\text{E} | h_{d,k} + \bs{h}_{r,k}^{H} \bm{\Theta}_0 \bs{g}_\text{E} |^2 = P_\text{E} |h_{d,k} + \bs{h}_k^H \bs{v}_0|^2.
		\end{aligned}
		\label{eq:P_EHk}
	\end{sequation}\noindent
	where $\bs{h}_k^H \triangleq \bs{h}_{r,k}^{H} \text{diag} (\bs{g}_\text{E})$, and $\bs{v}_0 \triangleq [e^{j\phi_1},\cdots,e^{j\phi_N}]^{T}$.

	\subsubsection{EH Model}
	To avoid the resource allocation mismatches associated with the traditional linear EH model, we adopt a more general and practical non-linear EH model, which is given by \cite{boshkovska2015practical}
	\begin{sequation}
		\begin{aligned}
			{\Phi}^{\text{non-linear}}_{{\text{EH}}_k} &= \frac{\Psi^{\text{non-linear}}_{{\text{EH}}_k} - M_k \Omega_k}{1-\Omega_k},  \\
			{\Psi}^{\text{non-linear}}_{{\text{EH}}_k} &= \frac{M_k}{1+e^{-a_k(P_{\text{EH}_k}-b_k)}}, \hspace{0.2cm}
			\Omega_k = \frac{1}{1+e^{a_k b_k}},
		\end{aligned}
		\label{eq:nonlinear_P}
	\end{sequation}\noindent
	where ${\Phi}^{\text{non-linear}}_{{\text{EH}}_k}$ denotes the amount of harvested energy at the $k^{th}$ device, and $\Psi^{\text{non-linear}}_{{\text{EH}}_k}$ is the conventional logistic function with respect to the received RF power $P_{{\text{EH}}_k}$. By adjusting the parameters $M_k$, $a_k$ and $b_k$, such a non-linear EH model is able to characterize the joint effects of various non-linear phenomena caused by hardware limitations \cite{boshkovska2015practical}.

	\subsubsection{NOMA Strategy}
	To improve SE and user fairness, NOMA technology is adopted for UL WIT
	with the help of SIC at the DS.
	For exposition purpose, we assume that the DS decodes devices' signals based on the indexes from 1 to $K$ in this paper, where the order of indexes are sorted according to the channel gain of direct link, i.e., $q_{d,k}$,  while the proposed algorithm and results are applicable to any given SIC order in practice.
	Thus, the achievable throughput of the $k$-th device in bits/Hz is given by \cite{wu2021irs}
	\begin{sequation}
		\begin{aligned}
			r_k = \tau_1 \log_2 \bigg( 1+ \frac{|q_{d,k} + \bs{q}_k^H \bs{v}_1|^2 p_k}{\sum\limits^K_{i=k+1} |q_{d,i} + \bs{q}_i^H \bs{v}_1|^2 p_i + \sigma^2} \bigg) ,
		\end{aligned}
		\label{eq:rk}
	\end{sequation}where $\bs{q}_k^H \triangleq \bs{g}_\text{I}^{H} \text{diag} (\bs{q}_{r,k})$,  $\bs{v}_1 \triangleq [e^{j\varphi_1},\cdots,e^{j\varphi_N}]^{T}$, and $\sigma^2$ is the additive white Gaussian noise power at the DS.

	\vspace{-0.4cm}
	\subsection{Problem Formulation}
	Different from existing works aiming at maximizing the sum throughput, we pursue a greener transmission policy, which leads to the objective of minimizing the total emitted energy at the PS, i.e., $\tau_0 P_\text{E}$, subject to the total available transmission time and user QoS constraints.
	Note that, since $P_\text{E}$ is usually a given constant due to practical limitation, the objective is equivalent to minimize the time duration $\tau_0$ for the DL WET. Accordingly, the time allocation, the transmit power at the IoT devices and the IRS phase shifts need to be jointly optimized, which leads to the following problem

		\vspace{-0.4cm}
		\begin{small}
		\begin{align}
			& { \min_{\scriptstyle \tau_0,\tau_1,\{ p_k \}, \atop \scriptstyle \bs{v}_0,\bs{v}_1} }  {\tau_0} \label{pro:A}\\
			& \hspace{0.5cm} {\text{s.t.}} \hspace{0.5cm}  { p_k \tau_1 \leq \tau_0 \Phi^{\text{non-linear}}_{{\text{EH}}_k}, \forall k, } \tag{\ref{pro:A}{a}} \label{pro:Aa}\\
			& \hspace{1.35cm} \tau_{0}+\tau_{1} \leq T_{\max } , \tag{\ref{pro:A}{b}} \label{pro:Ab} \\
			& \hspace{1.35cm}  r_k  \geq \Gamma_k, \forall k. \tag{\ref{pro:A}{c}} \label{pro:Af} \\
			& \hspace{1.35cm} \tau_{0} \geq 0, \tau_{1} \geq 0, p_k \geq 0, \forall k , \tag{\ref{pro:A}{d}} \label{pro:Ac} \\
			& \hspace{1.35cm} \left|\left[\bs{v}_0\right]_{n}\right|=1, \forall n \in \mathcal{N}, \tag{\ref{pro:A}{e}} \label{pro:Ad} \\
			& \hspace{1.35cm} \left|\left[\bs{v}_1\right]_{n}\right|=1, \forall n \in \mathcal{N}. \tag{\ref{pro:A}{f}} \label{pro:Ae}
		\end{align}
		\end{small}In problem (\ref{pro:A}), (\ref{pro:Aa}) and (\ref{pro:Ab}) represent the energy causality and total time constraints, respectively, (\ref{pro:Af}) are the QoS constraints of IoT devices, (\ref{pro:Ac}) are non-negativity constraints, (\ref{pro:Ad}) and (\ref{pro:Ae}) are unit-modulus constraints of the IRS phase shifts employed for DL WET and UL WIT, respectively. It can be shown that, problem (\ref{pro:A}) is a challenging non-convex optimization problem and difficult to solve optimally in general due to the highly coupled optimization variables in (\ref{pro:Aa}), as well as the non-convex unit-modulus and QoS constraints. In particular, when $\bs{v}_0$ is forced to be the same as $\bs{v}_1$ in problem (\ref{pro:A}), i.e., $\bs{v}_0=\bs{v}_1$, it means that static IRS beamforming is adopted, which will also be investigated for comparison.

	\vspace{-0.3cm}
	\section{Proposed Solution}

	In this section, we first decompose the original problem into two subproblems, namely the resource allocation problem and IRS passive beamforming problem. Then, we propose efficient solutions to the two subproblems respectively.

	\vspace{-0.34cm}
	\subsection{Resource Allocation Optimization}
	\vspace{-0.1cm}
	First, we design the resource allocation strategy for fixed reflection phase-shift vectors $\bs{v}_0$ and $\bs{v}_1$.
	Note that, different from the conclusion in \cite{wu2021irs} where the energy constraint (\ref{pro:Aa}) is met with {\it equality}, letting each device exhaust all its harvested energy may no longer be the optimal solution since each device may not need to use up all the harvested energy to meet with the minimum throughput constraint and a higher transmit power of one device may also cause stronger interference to other devices due to the use of NOMA for UL WIT. Therefore, $\{ p_k \}$ cannot be eliminated in the optimization as in \cite{wu2021irs}, which are coupled with the time variable $\tau_1$, thus making the problem more difficult to solve. Therefore, we introduce a new set of variables $\{ f_k \}$, where $f_k = p_k \tau_1, \forall k$. For notation simplicity, define $\widetilde{\bs{q}}_i^H \triangleq [\bs{q}_i^H \hspace{0.1cm} q_{d,i}]$, $ \widetilde{\bs{v}}_1^H \triangleq [\bs{v}_1^H \hspace{0.1cm} 1]$, $\bs{Q}_i \triangleq \widetilde{\bs{q}}_i \widetilde{\bs{q}}_i^H$, $\bs{V}_1 \triangleq \widetilde{\bs{v}}_1 \widetilde{\bs{v}}_1^H$, and $\boldsymbol{f} = [f_1, \cdots, f_K]^T$. The subproblem is reformulated as

	\vspace{-0.4cm}
	\begin{small}
		\begin{align}
			& \min_{\scriptstyle \tau_0,\tau_1,\{ f_k \}} \tau_0 \label{pro:1.1}\\
			&\hspace{0.45cm} {\text{s.t.}} \hspace{0.2cm}  { f_k  \leq \tau_0 \Phi^{\text{non-linear}}_{{\text{EH}}_k}, \forall k,} \tag{\ref{pro:1.1}{a}} \label{pro:1.1a}\\
			& \hspace{1.cm} (\ref{pro:Ab}), \hspace{0.2cm} \tau_{0} \geq 0, \tau_{1} \geq 0, f_k \geq 0, \forall k , \tag{\ref{pro:1.1}{b}} \label{pro:1.1b} \\
			& \hspace{1.cm}  \tau_1 \log_2 \bigg( 1+ \frac{\text{Tr} \big(\bs{Q}_k \bs{V}_1 \big) f_k}{\sum\limits^K_{i=k+1} \text{Tr} \big(\bs{Q}_i \bs{V}_1 \big) f_i + \sigma^2 \tau_1} \bigg) \notag  \geq \Gamma_k, \forall k. \tag{\ref{pro:1.1}{c}} \label{pro:1.1d}
		\end{align}
	\end{small}\noindent
	Problem (\ref{pro:1.1}) is still intractable due to the non-convexity of constraint (\ref{pro:1.1d}). To proceed, first we rewrite constraint (\ref{pro:1.1d}) as:
	\begin{sequation}
		\begin{aligned}
			h_k \big(\boldsymbol{f}, \tau_1 \big) - g_k \big(\boldsymbol{f}, \tau_1 \big) \geq {\Gamma_k}/{\tau_1}, \forall k ,
		\end{aligned}
		\label{eq:c1}
	\end{sequation}\noindent
	where
	\vspace{-0.4cm}
	\begin{small}
		\begin{align}
			&h_k \big(\boldsymbol{f}, \tau_1 \big) \triangleq \log_2 \bigg( \sum^K_{i=k} \frac{\text{Tr} \big(\bs{Q}_i \bs{V}_1 \big)}{\sigma^2} f_i + \tau_1 \bigg), \label{eq:c1.1.1}\\
			&g_k \big(\boldsymbol{f}, \tau_1 \big) \triangleq \log_2 \bigg( \sum^K_{i=k+1} \frac{\text{Tr} \big(\bs{Q}_i \bs{V}_1 \big)}{\sigma^2} f_i + \tau_1 \bigg) , \label{eq:c1.2.2}
		\end{align}
	\end{small}\noindent
	which yield to the classic DC form \cite{zargari2021max}. Note that $g_k \big(\boldsymbol{f}, \tau_1 \big)$ is a differentiable concave function with respect to both $\boldsymbol{f}$ and $\tau_1$, and is globally upper-bounded by its first-order Taylor expansion at any point. Therefore, we can obtain the upper bound of $g_k \big(\boldsymbol{f}, \tau_1 \big)$ at fixed local points ${\bs{f}}^r$ and ${\tau}^r_1$ as
	\begin{sequation}
		\begin{aligned}
			g_k \big(\boldsymbol{f}, \tau_1 \big) \leq \hspace{0.2cm} & g_k \big({\bs{f}}^r, {\tau}^r_1 \big) + \sum^K_{i=k+1}  \frac{\partial g_k \big({\bs{f}}^r, {\tau}^r_1 \big)}{\partial f_i}  (f_i - {f}^r_i) \\
			& + \frac{\partial g_k \big({\bs{f}}^r, {\tau}^r_1 \big)}{\partial \tau_1}  (\tau_1 - {\tau}^r_1) \triangleq \widehat{g}_k \big(\boldsymbol{f}, \tau_1 \big).
		\end{aligned}
	\end{sequation}\noindent
	By applying the SCA method, problem (\ref{pro:1.1}) can be reformulated in the following convex form

	\vspace{-0.5cm}
	\begin{small}
		\begin{align}
			& \min_{\scriptstyle \tau_0,\tau_1,\{ f_k \}} \tau_0  \label{pro:1.2}\\
			& \hspace{0.5cm} \text{s.t.} \hspace{0.5cm}   \text{(\ref{pro:1.1a}), (\ref{pro:1.1b})}, \hspace{0.2cm}   h_k \big(\boldsymbol{f}, \tau_1 \big) - \widehat{g}_k \big(\boldsymbol{f}, \tau_1 \big) \geq {\Gamma_k}/{\tau_1}, \forall k , \tag{\ref{pro:1.2}{b}} \label{pro:1.2b}
		\end{align}
	\end{small}\noindent
	which can be solved by convex optimization solvers such as CVX \cite{GRANT2014CVX}. After obtaining the optimal $\{ f_k \}$, $\tau_0$ and $\tau_1$, we can readily derive $\{ p_k \}$ as $p_k = f_k / \tau_1, \forall k$.

	\vspace{-0.3cm}
	\subsection{IRS Passive Beamforming Optimization}
	\vspace{-0.1cm}
	Next, we optimize the IRS phase-shift vectors $\bs{v}_0$ and $\bs{v}_1$ for the given $\tau_0$, $\tau_1$, and $\{ p_k \}$. Note that, since the optimization variables $\bs{v}_0$ and $\bs{v}_1$ are not coupled in the objective and constraints, we can separate this subproblem into two parallel feasibility check problems.

	\subsubsection{DL WET}
	In this stage, our target is to optimize $\bs{v}_0$, subject to (\ref{pro:Aa}) and (\ref{pro:Ad}). Define $\bs{H}_k \triangleq \widetilde{\bs{h}}_k \widetilde{\bs{h}}_k^H$, $\widetilde{\bs{h}}_k^H \triangleq [\bs{h}_k^H \hspace{0.2cm} h_{d,k}], \forall k$, $\bs{V}_0 \triangleq \widetilde{\bs{v}}_0 \widetilde{\bs{v}}_0^H$ and $\widetilde{\bs{v}}_0^H \triangleq [\bs{v}_0^H \hspace{0.2cm} 1]$. We recast (\ref{pro:Aa}) in a more explicit form which is given by
	\begin{sequation}
		\begin{aligned}
			\text{Tr} ( \bs{H}_k \bs{V}_0 ) \geq & \frac{b_k}{P_\text{E}} \! -  \! \frac{1}{a_k P_\text{E}} \ln \! \bigg( \! \frac{M_k}{\frac{1-\Omega_k}{\tau_0}p_k \tau_1+M_k \Omega_k} \! - \! 1 \! \bigg) , \! \forall k.
		\end{aligned}
		\label{eq:c1.4_v0}
	\end{sequation}\noindent
	To guarantee the rank and unit modulus constraints, the following condition should be satisfied as well.
	\begin{sequation}
		\begin{aligned}
			\text{Rank} ( \bs{V}_0 ) \leq 1, \hspace{0.2cm}
			\bs{V}_0 \succeq 0, \hspace{0.2cm} \text{diag} ( \bs{V}_0 ) = {\bf 1}_{N+1}.
		\end{aligned}
		\label{eq:v0_semi}
	\end{sequation}\noindent
	%
	Then, the optimal $\bs{V}_0$ can be obtained by solving
	\begin{small}
		\begin{align}
			&  \text{Find}  \hspace{0.8cm} \bs{V}_0 \label{pro:2.1.3}\\
			& \hspace{0.2cm} {\text{s.t.}} \hspace{0.7cm}  \text{ (\ref{eq:c1.4_v0}), (\ref{eq:v0_semi})} . \tag{\ref{pro:2.1.3}{a}} \label{pro:2.1.3a}
		\end{align}
	\end{small}\noindent
	Note that most prior works handle such kind of problems by discarding the rank-one constraint and then constructing the approximate solution by applying Gaussian randomization on the high rank solution \cite{wu2019intelligent, wu2019towards, wu2021irs}. However, such a scheme performs far from satisfactory within the AO process. Therefore, we exploit the penalty-based method as in \cite{zargari2021max} instead. Specifically, first we reformulate problem (\ref{pro:2.1.3}) as

	\vspace{-0.5cm}
	\begin{small}
		\begin{align}
			& \ { \min_{\scriptstyle \bs{V}_0} } \hspace{0.8cm} \mu\bigg( \text{Tr} ( \bs{V}_0 ) - \| \bs{V}_0 \|_2 \bigg) \label{pro:2.1.4}\\
			& \hspace{0.3cm} {\text{s.t.}} \hspace{0.7cm}   (\ref{eq:c1.4_v0}), \hspace{0.2cm} \bs{V}_0 \succeq 0, \hspace{0.2cm} \text{diag} ( \bs{V}_0 ) = {\bf 1}_{N+1}, \tag{\ref{pro:2.1.4}{a}} \label{pro:2.1.4a}
		\end{align}
	\end{small}\noindent
	where $\mu$ is the penalty factor. It is worth noting that although the rank-one constraint is relaxed in (\ref{pro:2.1.4}), the solution yields a rank-one solution when $\mu$ is sufficiently large. Nevertheless, (\ref{pro:2.1.4}) is not convex yet due to the DC form of the objective function. To tackle such non-convexity, we utilize the first-order Taylor expansion to obtain the lower bound of $\| \bs{V}_0 \|_2$ and rewrite the objective of (\ref{pro:2.1.4}) as

	\vspace{-0.5cm}
	\begin{small}
		\begin{align}
			 z(\bs{V}_0) \triangleq \mu\bigg( \text{Tr} ( \bs{V}_0 ) - \| {\bV}^r_0 \|_2 - \text{Tr} \big[ {\bv}^r_{\max}{\bv}^r_{\max}^H( \bs{V}_0 - {\bV}^r_0) \big] \bigg), \label{pro:2.1.5}
		\end{align}
	\end{small}\noindent
	where ${\bv}^r_{\max}$ is the eigenvector corresponding to the maximum eigenvalue of matrix ${\bV}^r_0$. Afterwards, we can solve this problem with CVX tool \cite{GRANT2014CVX}.

	 To achieve better converged solution, we introduce a ``residual'' variable $\delta$.
	 The rationale is that optimizing the phase shift to enforce the user's equivalent cascaded channel gain to be larger than the original channel gain target in (\ref{pro:2.1.DC.delta.b}) improves the channels' condition, which
	 leads to the reduction of the DL WET time duration $\tau_0$ in problem (\ref{pro:A}). For this consideration, the slack variable $\delta$ can be interpreted as the “channel gain residual” of the users in phase shift optimization, which is beneficial in terms of the convergence manner \cite{wu2019intelligent}. Thus, problem (\ref{pro:2.1.4}) can be further transformed as

	\vspace{-0.4cm}
	\begin{small}
		\begin{align}
			&  \min_{\scriptstyle \bs{V}_0, \delta}  \hspace{0.3cm} z(\bs{V}_0) - \delta  \label{pro:2.1.DC.delta}\\
			& \hspace{0.3cm} {\text{s.t.}} \hspace{0.3cm}  \bs{V}_0 \succeq 0, \hspace{0.2cm} \text{diag} ( \bs{V}_0 ) = {\bf 1}_{N+1},  \tag{\ref{pro:2.1.DC.delta}{a}} \label{pro:2.1.DC.delta.a} \\
			& \hspace{1cm} \text{Tr} ( \bs{H}_k \bs{V}_0 ) \geq \delta + \frac{b_k}{P_\text{E}} - \notag \\
			& \hspace{1.6cm} \frac{1}{a_k P_\text{E}} \ln \bigg( \frac{M_k}{\frac{1-\Omega_k}{\tau_0}p_k \tau_1+M_k \Omega_k} -1 \bigg), \forall k. \tag{\ref{pro:2.1.DC.delta}{b}} \label{pro:2.1.DC.delta.b}
		\end{align}
	\end{small}\noindent
	Finally, the optimal $\bs{V}_0$ can be obtained by solving problem (\ref{pro:2.1.DC.delta}) with penalty base method iteratively until the algorithm converges \cite{zargari2021max}.

	\subsubsection{UL WIT}
	To optimize $\bs{V}_1$, first we rewrite the constraint (\ref{pro:Af}) as
	\vspace{-0.2cm}
	\begin{sequation}
		\begin{aligned}
			& p_k \text{Tr} ( \bs{Q}_k \bs{V}_1 )  \geq  \left( \! 2^{\frac{\Gamma_k}{\tau_1}} \! - \! 1 \! \right) \! \bigg(\! \sum\limits^K_{i=k+1} p_i  \text{Tr} ( \bs{Q}_i \bs{V}_1 )  + \sigma^2 \! \bigg) , \forall k.
		\end{aligned}
		\label{eq:c6.1_v1}
	\end{sequation}\noindent
	Then, the optimal $\bs{V}_1$ can be obtained similarly as that for problem (\ref{pro:2.1.4}), except replacing the constraint (\ref{eq:c1.4_v0}) by (\ref{eq:c6.1_v1}). Thus, the detail is omitted here for simplicity.

	\vspace{-0.3cm}
	\subsection{Static IRS Beamforming}
	To reduce the algorithmic computations and lower the feedback signalling overhead from the PS to the IRS for sending the optimized IRS phase shifts, we consider a special case of employing static IRS beamforming where the IRS phase shifts remain unchanged in DL WET and UL WIT, i.e., $\bs{V}_0=\bs{V}_1\triangleq \bs{V}$. The optimal IRS passive beamforming can be obtained by solving

	\vspace{-0.4cm}
	\begin{small}
		\begin{align}
			  \text{Find}\hspace{0.1cm}  \hspace{0.2cm} \bs{V} \hspace{6.5cm} & \label{pro:2.2}\\
			  {\text{s.t.}}\hspace{0.2cm}   \text{Tr} ( \bs{H}_k \bs{V} ) \geq  \frac{b_k}{P_\text{E}} -  \hspace{4.5cm} & \notag \\
			  \hspace{1.8cm}  \frac{1}{a_k P_\text{E}} \ln \! \bigg( \! \frac{M_k}{\frac{1-\Omega_k}{\tau_0}p_k \tau_1+M_k \Omega_k} \! - \! 1 \! \bigg) , \! \forall k, & \tag{\ref{pro:2.2}{a}} \label{pro:2.2a} \\
			 p_k  \text{Tr} ( \bs{Q}_k \bs{V} )  \geq  \left( \! 2^{\frac{\Gamma_k}{\tau_1}} \! - \! 1 \! \right) \! \bigg(\! \sum\limits^K_{i=k+1} p_i  \text{Tr} ( \bs{Q}_i \bs{V} )  + \sigma^2 \! \bigg) , \forall k, & \tag{\ref{pro:2.2}{b}} \label{pro:2.2b} \\
				\text{Rank} ( \bs{V} ) \leq 1, \hspace{0.2cm}
			\bs{V} \succeq 0, \hspace{0.2cm} \text{diag} ( \bs{V} ) = {\bf 1}_{N+1}, \hspace{1.7cm} & \tag{\ref{pro:2.2}{c}}
		\end{align}
	\end{small}\noindent
	which can be tackled by the aforementioned penalty-based method algorithm as well. Note that the number of optimization variables is reduced greatly in this case, especially when the number of IRS elements is practically large, which might be a compromise choice when the computation complexity is a more critical issue in some certain scenarios, striking a balance between the system performance and signaling overhead as well as implementation complexity.

	To sum up, the overall joint design problem is solved by alternately optimizing the resource allocation and IRS beamforming until convergence.
	More rigorously, the convergence of such AO-based algorithm is ensured since the objective value of problem (\ref{pro:A}) is non-increasing over the AO iterations \cite{wu2019intelligent}.
	In addition, the complexity of the AO-based algorithm is analyzed as follows. For the resource allocation subproblem, i.e., problem (\ref{pro:1.1}), the approximate computational complexity is given by $\mathcal{O}(K^{3.5})$. While for the subproblem of IRS passive beamforming optimization, the computational complexity is given by $\mathcal{O}(N^{6.5})$. Thus, the computational complexity of the overall AO algorithm is $\mathcal{O}(I_{AO} (K^{3.5}+N^{6.5}))$, where $I_{AO}$ denotes the number of iterations required for convergence.

	\vspace{-0.2cm}
	\section{Numerical Results}

	In this section, we consider an IRS-aided WPCN where the PS, IRS, and DS are located at (-5,0,0) meter (m), (0,0,1) m, and (45,0,0) m, respectively. The IoT devices are randomly and uniformly distributed within a radius of 1 m centered at (0,0,0) m.
	For large scale fading, the distance-dependent path loss model with respect to the link distance $d$ is given by $PL(d) = c_0(d/d_0)^{-\alpha}$, where $c_0=-30$ dB is the signal attenuation at a reference distance of $d_0 = 1$ m, and $\alpha$ is the pathloss exponent, which is set to be 2.2 for the PS-IRS, IRS-DS and IRS-device channels, while 2.8 for the PS-device and device-DS channels.
	For small scale fading, Rayleigh fading is adopted for PS-device and device-DS direct link, $h_{d,k}$ and $q_{d,k}$, $\forall k$. While for $\bs{g}_\text{E}$, $\bs{g}_\text{I}$, $\bs{h}_{r,k}$ and $\bs{q}_{r,k}$, $\forall k$, Rician fading is adopted with a Rician factor of $10$. The parameters for non-linear EH model are set as $a_k = 150, b_k = 0.014$, and $M_k = 0.024$ \cite{zargari2021max}. The other system parameters are set as follows: $P_{\text{E}}=40$ dBm, $\sigma^2=-
	90$ dBm, $T_{\max}=1$ s, $\Gamma =0.1$ bits/Hz, $K=3$ and $N=100$ \cite{wu2021irs}, if not specified otherwise.

	For comparison, we consider the following schemes.  1) \textbf{Proposed AO with dynamic IRS:} employ different IRS phase shifts in DL and UL with the proposed solutions;
	2) \textbf{GR with dynamic IRS:} where Gaussian randomization is applied to recover the rank-one $\bs{V}_0$ and $\bs{V}_1$ based on the SDR solution \cite{wu2019intelligent}; 3) \textbf{ETA with dynamic IRS:} i.e., $\tau_0=\tau_1$ based on the solution of the scheme in 1);
	4) \textbf{Proposed AO with static IRS:} use the same IRS phase shifts in DL and UL, i.e., $\bs{V}_0 = \bs{V}_1$ with the proposed solutions; 5) \textbf{Fixed IRS phase shift:} Only optimize resources allocation, i.e., $\tau_0,\tau_1,\{ p_k \}$, with randomly initialized IRS phase shifts; 6) \textbf{Without IRS.}

	\begin{figure*}[!t]
		\centering
		\subfigure[Optimized $\tau_0$ versus $N$.]{
			\label{fig:N_t0}
			\includegraphics[width=0.4885\columnwidth,height=0.487\columnwidth]{figure2_a.eps}}
		\subfigure[Optimized $\tau_0$ versus $P_\text{E}$.]{
			\label{fig:PA_t0}
			\includegraphics[width=0.4885\columnwidth,height=0.4885\columnwidth]{figure2_b.eps}}
		\subfigure[Transmit energy  $\tau_0 P_\text{E}$ versus $P_\text{E}$.]{
			\label{fig:PA_E}
			\includegraphics[width=0.4885\columnwidth,height=0.4885\columnwidth]{figure2_c.eps}}
		\subfigure[Effect of imperfect CSI.]{
			\label{fig:imCSI}
			\includegraphics[width=0.4885\columnwidth,height=0.5\columnwidth]{figure2_d.eps}}
			\caption{Simulation results.}
	\end{figure*}

	First, we investigate the impact of the number of IRS elements. As shown in Fig. \ref{fig:N_t0}, the optimized $\tau_0$ keeps decreasing as the number of IRS elements increases for all the strategies, except the benchmark of without IRS. This is expected since for the IRS-aided WPCN, the increasing number of IRS elements brings more cascade channels, which enables the higher channel gain. Thus, the DL WET duration $\tau_0$ can be lowered, which shows the benefits of IRS for WPCNs since we can reduce the PS's consumed transmit power for DL WET by deploying IRS with a large number of low-cost reflection elements.
	Besides, it can be seen that when employing dynamic IRS, the optimized $\tau_0$ obtained by the proposed algorithm is much lower compared with that of the system with "GR", "ETA" and "Fixed IRS phase shifts", which reveals the necessity of properly designed joint resource allocation and phase shift optimization.
	Furthermore, it is observed that the "Proposed AO with dynamic IRS" outperforms "Proposed AO with static IRS", indicating the benefits of taking different IRS configurations over the DL WET and UL WIT, which is different from the conclusion in \cite{wu2021irs}, where using same IRS phase shifts suffices to achieve the optimal performance with a HAP.
	This results from the considered system model and problem formulation.
	In detail, employing separated PS and DS architecture leads to different UL/DL channels. Besides, using the non-linear energy harvesting model in DL WPT further makes DL channels asymmetric as the UL channels.
	This result highlights the advantage of dynamic IRS beamforming as well as the superiority of the proposed design among the benchmarks.

	Fig. \ref{fig:PA_t0} and Fig. \ref{fig:PA_E} present respectively the optimized DL WET duration and the corresponding transmit energy consumption as the transmit power $P_\text{E}$ at the PS varies. Specifically, it is observed from Fig. \ref{fig:PA_t0} that the optimized $\tau_0$ for DL WET declines as the transmit power increases, because increasing the transmit power is able to shorten the energy charging time to meet with the QoS constraint in UL WIT. Nevertheless, we find from Fig. \ref{fig:PA_E} that increasing the transmit power leads to an increase in total transmit energy consumption $\tau_0 P_\text{E}$, indicating that increasing the transmit power is not a favorable choice to reduce the overall transmit energy consumption in practical WPCNs.
	This is due to the fact that we adopt non-linear EH model in the paper, the harvest energy is not linearly proportional to the increase of the transmit power, which leads to that the reduction of the optimized DL WET time duration $\tau_0$ is not linearly proportional to the increase of the transmit power. From Fig. \ref{fig:PA_E}, it is easy to infer that the increasing of the transmit power is more prominent than the reduction of time $\tau_0$. Therefore, the total transmit energy is increasing with large transmit power.
	By comparing Fig. \ref{fig:N_t0} and Fig. \ref{fig:PA_E}, it is not difficult to conclude that we can save the transmit energy consumption by deploying more low-cost reflecting elements instead of raising the emitted power at the PS, which further demonstrates the benefits of integrating properly designed IRS into WPCNs.

	In addition, we also examine the effect of imperfect CSI on the system performance as shown in Fig.\ref{fig:imCSI}, where $\epsilon$ denotes the value of channel estimation error.
	It can be seen that, when the available CSIs are imperfect, the performance of the proposed algorithm decreases as expected. And this degradation gets worsen as the channel estimation errors $\epsilon$ are increasing. However, the observations regarding to the comparison between dynamic IRS beamforming and static IRS beamforming in this letter still hold even with estimation errors of CSIs. More importantly, we also observe that as $N$ increases, the performance gap between the perfect CSI and imperfect CSI cases becomes smaller, which implies that the performance loss caused by imperfect CSI can be effectively alleviated by using a large number of reflecting elements and thus validates the practically of the proposed method.

	\vspace{-0.2cm}
	\section{Conclusions}

	In this paper, we investigated the joint passive beamforming and resource allocation optimization for an IRS-aided WPCN employing NOMA for UL WIT. In particular, the total emitted energy at the PS was minimized subject to the total available transmission time and devices' QoS constraints under a non-linear EH model. The joint optimization problem was decomposed into two subproblems and solved respectively in an alternative manner by using the SCA technique and penalty based algorithm. Numerical results not only validated the effectiveness of the proposed AO algorithm but also unveiled the following useful insights. 1) Using dynamic IRS beamforming in DL and UL outperforms employing static IRS beamforming for separated PS and DS deployment. 2) The transmit energy consumption can be effectively reduced by deploying more low-cost IRS reflecting elements instead of raising the emitted power at the PS.

	\vspace{0.05cm}

	\bibliographystyle{IEEEtran}


\section{Introduction}
\label{introduction}
Text summarization is studied in two paradigms: extractive and abstractive. Different from extractive models which directly select text spans from source content, abstractive approaches can generate summaries more concisely with flexible information aggregation and paraphrasing, but raise the higher requirement for context understanding and language generation.
In the past few years, various neural sequence-to-sequence models are proposed and applied to abstractive summarization, from the vanilla attention-based summarizer \cite{rush-etal-2015-neuralSumm}, pointer-generator networks \cite{see-2017-pgnet}, to the recent large-scale pre-trained language models such as BERT and BART \cite{devlin-2019-BERT,lewis-etal-2020-bart}. These models have brought significant improvement on the generation quality especially in fluency and readability \cite{liu-lapata-2019-PreSumm}.

While most prior work focuses on documents such as news articles \cite{hermann-2015-cnnDaily}, summarizing conversations starts to raise more research interest \cite{goo2018-gatedDiaSumm,zhu2020-SummHNet}. One typical multi-party conversation scenario is meeting \cite{mccowan2005ami}. Since meetings often aim to cover several sub-topic discussions, the average number of their dialogue turns is much larger than that of short conversations such as social chat \cite{gliwa-etal-2019-samsum} and inquiring-answering \cite{liu-2019-topic}.
As the source content becomes longer, the summarizer may fall distracted on the massive contextual information, and the quality of neural generation tends to degrade when producing long text \cite{fan2018StoryGeneration}. Thus it is challenging for neural models to achieve the overall high quality for the long transcript summarization. Moreover, if adopting Transformer-based pre-trained language models like BART \cite{lewis-etal-2020-bart}, the transcript length often exceeds their maximum positional embedding limitations, and conducting text truncation will cause loss of some context.

Meetings usually consist of a set of agenda items \cite{mccowan2005ami}, and each item focuses on discussing a specific topic. Thus the whole meeting conversation is inherently organized at the topic level. This is also reflected in the summaries written by humans, which condense the key information on each discussion part.
Therefore, we postulate the divide-and-conquer strategy can be useful to tackle the aforementioned challenges of meeting summarization. More specifically, our target is to split the long transcript into multiple segments, and obtain the final output by aggregating the summary snippets of each segment. To construct the ground-truth data for training, we first investigate the linguistic characteristics of meeting transcripts and reference summaries on the representative corpus AMI \cite{mccowan2005ami}.
Then, we propose to enhance the sequence-to-sequence summarizer with a dynamic sliding window strategy, to tackle the length limitation of Transformer-based language models. Unlike the conventional fixed window sliding \cite{devlin-2019-BERT}, the dynamic method makes a model learn to decide the start position of the next segment during the generation process.
Experimental results show that the summary generation can benefit from adopting the dynamic window sliding, and achieves state-of-the-art performance. Further example analysis demonstrates that the enhanced framework performs better than the base model considering factual consistency.

In the contemporary studies, efficient attention schemes of the Transformer \cite{beltagy2020longFormer} and multi-block context aggregation \cite{grail2021globalizing} are proposed to support long text processing, and \citet{schuller2020windowing} investigate the fixed and automatic sliding window for processing long documents, and they focused on summarizing Wikipedia and news articles. Our proposed method is inspired by the organizational features of meeting transcripts, and is applicable for the language backbones with input length limitations.

\begin{figure}
    \begin{center}
    \includegraphics[width=0.49\textwidth]{images/Meeting_Example.png}
    \end{center}
    \caption{One example of the meeting transcript and reference summary in the AMI corpus. Here we only show several extractive utterances in three sub-topic chunks and their corresponding abstractive sentences.}
    \label{fig:meet-example}
\vspace{-0.3cm}
\end{figure}

\begin{table}
\linespread{1.1}
    \centering
    \small
    \begin{tabular}{lcccc}
    \hline
         & CNN & DailyMail & NYT & AMI \\
    \hline
         \multicolumn{4}{l}{Average Source Content Length:} \\
         Sentence Level & 33.98 & 29.33 & 35.55 & 288.7 \\
         Word Level & 760.5 & 653.3 & 800.1 & 4757 \\
    \hline
    \multicolumn{4}{l}{Average Reference Summary Length:} \\
        Sentence Level & 3.59 & 3.86 & 2.44 & 17.55 \\
         Word Level & 45.8 & 54.65 & 45.54 & 323.3 \\
    \hline
    \end{tabular}
    \caption{Data statistics of the news summarization benchmarks and the AMI meeting corpus.}
    \label{tab:dataset}
\vspace{-0.2cm}
\end{table}

\section{Meeting Transcript Analysis}
\label{sec:meeting_analysis}
In this section, we conduct data analysis on the AMI meeting corpus \cite{mccowan2005ami}, in which the participants work in a team and conduct meetings to discuss product design, development, and planning. There are four main speaker roles: a project manager (PM), a marketing expert (ME), an industrial designer (ID), and a user interface designer (UI). Following previous work \cite{shang-etal-2018-unsupervised}, we split the whole dataset into train (97 transcripts), validation (20 transcripts), and test (20 transcripts) sets. The data statistics are shown in Table \ref{tab:dataset}, where we count one utterance in the conversation as one sentence, and conduct word-level tokenization. Compared with the news summarization benchmarks, the average length of meeting transcripts as well as the reference summaries are much larger (In our settings, we use the long version abstracts in the AMI corpus, as previous work \cite{feng2020-DDAGCN}). Moreover, human-written summaries of news articles often concentrate on the first few parts of source content \cite{liu2020conditional}, thus the truncation with a fixed length does not affect the final performance significantly \cite{jung2019earlier}. However, summarizing meetings requires grasping useful contextual information across the whole conversation, in this case, simple text truncation will lead to certain information loss.

\begin{figure}
    \begin{center}
    \includegraphics[width=0.43\textwidth]{images/Meeting_Coref.png}
    \end{center}
    \caption{One example of meeting summary in the AMI corpus. The colored spans are named entities or pronouns in a coreference chain.}
    \label{fig:coref-example}
\vspace{-0.3cm}
\end{figure}

\begin{figure*}[ht!]
    \begin{center}
    \includegraphics[width=0.82\textwidth]{images/Meeting_framework.png}
    \end{center}
    \caption{Overview of the meeting summarization framework with the dynamic sliding window strategy.}
    \label{fig:framework}
\vspace{-0.2cm}
\end{figure*}

\subsection{Organization Analysis}
\label{ssec:organization_analysis}
Since meetings are usually taken with a set of sub-topics, we analyze the organizational co-relation of these sub-topics and sentences in meeting summaries.
In the AMI corpus, for each meeting transcript, aside from a human-written abstractive summary, there is an additional extractive annotation. As shown in Figure \ref{fig:meet-example}, human annotators were asked to choose sentences from the source conversation, as the supporting segments for each abstractive summary sentence. Based on this extractive annotation, we obtained the start/end supporting utterance indices in the conversation of each summary sentence. Then we sorted them by their occurrence order in the source content, and observed that 80\% sentences in reference summaries match the same agenda order as in the source content. This indicates that when humans write a summary, they often read the source content and record the key points sequentially.

\subsection{Segmentation Analysis}
\label{ssec:segment_analysis}
Summarization is a process to write a shorter version of a text span, and context integrity is essential. When adopting a divide-and-conquer strategy, only a part of the conversation will be extracted. To assess information integrity, we calculate the word level coverage with the extractive annotations (Figure \ref{fig:meet-example}) in the AMI. With the start/end supporting utterance indices described in Section \ref{ssec:organization_analysis}, we split the original transcripts into a number of segments. Then the word-level recall is calculated between each summary sentence and its corresponding segment. Summaries are often written via paraphrasing and introducing novel words, and the overall word-level coverage in our observation is 74\%.

\subsection{Summary Conversion}
\label{ssec:summary_conversion}
As the final output under a divide-and-conquer strategy is produced by aggregating summary snippets, to build the training ground-truth, we split each reference summary into multiple parts based on the analysis in Section \ref{ssec:organization_analysis} and Section \ref{ssec:segment_analysis}. For each summary sentence, we constructed its context segment by using the corresponding extractive annotation. Then we ordered the summary sentences by their occurrence indices in the source conversation, and merged adjacent summary sentences to summary snippets if their context segments had a certain overlap. For summary sentences without extractive annotation, we reorganized them via calculating word-level overlap with existing content segments.
Furthermore, we observed that some summary snippets are started with a pronoun that refers to a precedent personal named entity, as shown in Figure \ref{fig:coref-example}. To maximize the semantic integrity when producing one snippet, we used coreference resolution on the original summary, and if one summary snippet is started with a personal pronoun, it will be replaced by its referring.

\section{Neural Meeting Summarizer}
In this section, we elaborate our framework for meeting summarization with the dynamic sliding window strategy.

\subsection{Dynamic Sliding Window}
The sliding window is to split a long input into a number of shorter spans, process the spans in order (e.g., from left to right), and adopt aggregation for final output. It is a straightforward but useful method that is commonly applied to long sequence encoding \cite{devlin-2019-BERT}, and it is generally controlled by two parameters: \textit{window size} denotes the context span size, and \textit{stride size} is the amount of movement at each sliding step.
In previous work, these two parameters are often fixed, here we propose a dynamic sliding window strategy, where the stride size at each sliding step is predicted by the model. More specifically, given the window size $k$, at each sliding step, one index $j$ is selected (equals stride size is $j$) in the range of $[0, k-1]$ as the start position of next window.

\subsection{Base Neural Architecture}
For the neural summarizer, we use a Transformer-based sequence-to-sequence architecture \cite{vaswani-2017-Transformer}, and select the large-scale pre-trained language backbone BART \citep{lewis-etal-2020-bart} for initialization. The encoder consists of 6 stacked Transformer layers, and each layer has two sub-components: a multi-head self-attention layer and a position-wise feed-forward layer. Between the two sub-components, residual connection and layer normalization are used. The $u$-th encoding layer is formulated as:
\begin{equation}
\label{eq-attn-head}
    \widetilde{h}^u=\mathrm{LN}(h^{u-1}+\mathrm{MHAttention}(h^{u-1}))
\end{equation}
\vspace{-0.5cm}
\begin{equation}
\label{eq-attn-residual}
    h^u=\mathrm{LN}(\widetilde{h}^u+ \mathrm{FFN}(\widetilde{h}^u))
\end{equation}
where $h^0$ is the first layer input. $\mathrm{MHAttention}$, $\mathrm{FNN}$, $\mathrm{LN}$ are multi-head attention, feed-forward, and layer normalization, respectively.

The decoder consists of 6 stacked Transformer layers as well. In addition to the two encoding sub-components, the decoder performs another multi-head attention over the previous decoding hidden states and all the encoded representations. Then, the decoder generates tokens in an auto-regressive manner from left to right.

\subsection{Adopting Dynamic Sliding Window for Meeting Summarization}
\label{ssec:adopting_sliding}
The proposed dynamic sliding window is a general design that can be applied to various neural architectures.
In our summarization setting, the source content is a sequence with $n$ tokens $C=\{w_1, w_2,..., w_n\}$, and the summary is composed of $m$ snippets $S = \{s_1, s_2,..., s_m\}$. In a traditional sequence-to-sequence process, the whole source content $C$ is fed to a model as input, and the output summary is generated at one decoding stage. With the dynamic sliding window strategy, the encoding-decoding process will be conducted in a number of steps, as shown in Figure \ref{fig:framework}. Given the window size $k$, at $i$-th sliding step, the start token index is $i_{left}$, and the end token index $i_{right}$ is $i_{left} + k$. Then the input sequence of encoder at $i$-th step is $c_i=\{w_{i_{left}},..., w_{i_{right}}\}$.
After encoding, the contextualized hidden representation $h_i=\{h_{i_{left}},..., h_{i_{right}}\}$ is fed to the decoder to generate a summary snippet, and all snippets will be merged to form the final output.

To obtain the dynamic sliding prediction, one way is to directly predict the value of stride size as in \cite{gong-etal-2020-recurrentChunk}, where reinforcement learning is used to decide the stride value, and \citet{schuller2020windowing} proposed to generate a special indicator for the moving operation. Here, based on the analysis in Section \ref{sec:meeting_analysis}, we introduce a learning-efficient method named \textit{Retrospective}. More specifically, as shown in Figure \ref{fig:framework}, at $i$-th step, the decoder will not only generate the summary snippet $s_i$, but also the last supporting utterance of its context segment that is described in Section \ref{ssec:organization_analysis} (concatenated with a ``[SEP]'' token). Since the predicted supporting utterance provides the context boundary of the current generation step on the right, it is used to determine the amount of sliding movement. After all sliding steps, the generated snippets are concatenated as the output summary.

\section{Experiment Results and Analysis}

\subsection{Configuration}
The proposed framework was implemented using PyTorch \citep{paszke2019pytorch} and Hugging Face \citep{wolf-etal-2020-transformers}. The learning rate was set at 2e-5, and AdamW \cite{loshchilov2017-adamW} optimizer was applied. We trained each model for 20 epochs, and selected the best checkpoints on the validation set based on ROUGE-2 score \cite{lin-2004-rouge}. All input sequences were processed with the sub-word tokenization scheme as in \cite{lewis-etal-2020-bart}, and repeated output sentences were removed.
Based on the summary conversion described in Section \ref{ssec:summary_conversion}, we obtained 598/131/143 snippet-level samples as the train, validation, and test set, respectively. In the training stage, to simulate the input noise during inference time, we randomly added $k$ adjacent utterances ($5<k<15$) in each context chunk.
In the testing stage, for each example, the generation process of our \textbf{BART-SW-Dynamic} model was started with the first chunk of a default window size, which is initialized as 1024.\footnote{The token-level maximum input length of BART \cite{lewis-etal-2020-bart} is 1024 by default.} Then we used the sliding prediction described in Section \ref{ssec:adopting_sliding} to form the next chunk, until the whole transcript was processed. We also assessed the \textbf{BART-SW-Fixed} model by fixing the window size at 1024, the \textbf{BART-Truncate} model by truncating the transcript of 1024 tokens, and the \textbf{BART-SW-GoldSeg} model by using the gold context segmentation.

\begin{table}[t!]
\linespread{1.0}
\centering
\begin{tabular}{p{3.2cm}ccc}
\hline
                & \textbf{R-1} & \textbf{R-2} & \textbf{R-L}  \\
\hline
\multicolumn{2}{l}{Extractive baselines*:} \\
TextRank   &  35.19   &  6.13 &  15.70  \\
SummaRunner   &  30.98  &  5.54 &  13.91  \\
\hline
\hline
\multicolumn{2}{l}{Abstractive baselines*:} \\
Seq2Seq+Attention  &  36.42   &  11.97 &  21.78 \\
Pointer-Generator   &  42.60   &  14.01 &  22.62  \\
Sentence-Gated   &  49.29   & 19.31  &  24.82  \\
TopicSeg   &  51.53   & 12.23  &  25.47  \\
HMNet   & 53.02 & 18.57 & 24.85 \\
DDA-GCN  & 53.15  & 22.32 & 25.67 \\
\hline
\hline
\multicolumn{3}{l}{Our models:} \\
BART-SW-GoldSeg & 53.42  & 22.57 & 28.52 \\
BART-Truncate & 48.31  & 17.52 & 20.13 \\
BART-SW-Fixed & 50.02  & 19.86 & 23.98 \\
BART-SW-Dynamic & 52.83  & 21.77 & 26.01 \\
\hline
\end{tabular}
\caption{\label{output-ROUGE-table} ROUGE F1 scores on the AMI test set from baseline models and our framework. * Results of baselines are reported as in \cite{feng2020-DDAGCN}.}
\vspace{-0.2cm}
\end{table}

\subsection{Results on the AMI Corpus}
Following previous work \cite{chen-yang-2020-MultiView,feng2020-DDAGCN}, we used the ROUGE score \citep{lin-2004-rouge} for generation assessment, and reported ROUGE-1 (R-1), ROUGE-2 (R-2) and ROUGE-L (R-L) scores. We selected a set of strong baseline models for extensive comparison including Pointer-Generator \cite{see-2017-pgnet}, Sentence-Gated \cite{goo2018-gatedDiaSumm}, TopicSeg \cite{li-etal-2019-keepMeeting}, HMNet \cite{chen-yang-2020-MultiView}, and DDA-GCN \cite{feng2020-DDAGCN}.
As shown in Table \ref{output-ROUGE-table}, for meeting summarization, abstractive approaches generally perform much better than extractive ones. In our settings, the model with the dynamic sliding window strategy (BART-SW-Dynamic) outperforms the model with the fixed window (BART-SW-Fixed), and the model with text truncation (BART-Truncate). This shows the effectiveness of the sliding window approach. Moreover, our proposed framework achieves the comparable performance of the contemporary state-of-the-art models.

\subsection{Sample Analysis}
We first conducted a text quality analysis on the generated summaries across models. As shown in Figure \ref{fig:compare-example}, based on the strong generation capability of the language backbone, all BART-based models can produce fluent and grammatically correct summaries. While the model with text truncation achieves relatively acceptable ROUGE scores (as reported in Table \ref{output-ROUGE-table}), it produces sentences that are factually inconsistent with the source content, as shown in Figure \ref{fig:compare-example}. We speculate that this is caused by over-fitting the training samples. In contrast, since our proposed framework can produce the final summary based on relevant context segments without truncation, it performs better than the base model considering the factual consistency.

We then conducted a stride prediction assessment. For the \textit{Retrospective} method described in Section \ref{ssec:adopting_sliding}, the predicted context boundaries are expected to be located closely to the ground-truth. With the best checkpoint, we observed that the average utterance-level distance between gold boundary span and model prediction is 2.7 (48 characters), and this shows that the model is able to predict the correct start position at each sliding step.

\begin{figure}[t!]
    \begin{center}
    \includegraphics[width=0.48\textwidth]{images/Meeting_Compare.png}
    \end{center}
    \caption{Generated summaries. Spans in blue and pink are consistent with the gold summary and source content. Spans in yellow are factually inconsistent.}
    \label{fig:compare-example}
\vspace{-0.3cm}
\end{figure}

\section{Conclusion}
In this paper, to tackle the challenges from lengthy meeting transcript inputs, we proposed a dynamic sliding window strategy for abstractive summarization. Experimental results demonstrate that the neural sequence-to-sequence models can benefit from the proposed method, and suggest that the long transcript summarizing can be conducted in a divide-and-conquer manner. One of the future works can be extending the proposed method to other corpora with larger sample sizes.

\bibliographystyle{acl_natbib}


\section{Introduction}
Nanomechanical vibrational systems provide a natural platform for studying a broad range of classical and quantum phenomena in a well characterized setting, cf. Refs. \cite{Cattiaux2021,Siskins2021,Tepsic2021,Keskekler2021,Ari2020,MacCabe2020} for recent examples. An important advantageous feature of such systems is a small decay rate, with the ratio $Q$ of the vibration eigenfrequency $\omega_0$ to the decay rate $2\Gamma$ reaching  $5\times 10^{10}$ \cite{MacCabe2020} for localized acoustic modes and $8\times 10^8$ for flexural modes \cite{Ghadimi2018}. This makes nanomechanical modes highly sensitive to a resonant force, which underlies many of their applications. A consequence of this sensitivity is that a comparatively weak resonant force can drive the vibrations into a regime where their nonlinearity comes into play. This enables using resonantly driven nanomechanical vibrations for studying various nonlinear phenomena far from thermal equilibrium, cf. Refs.~ \cite{Aldridge2005,Chan2008a,Segev2008,Defoort2015,Dolleman2019} and papers cited therein.

In many cases, the nonlinear response of nanomechanical vibrations to a comparatively weak resonant drive is well described by the Duffing model \cite{Kozinsky2007}. In this model the nonlinearity comes from the term in the
potential energy, which is quartic in the mode coordinate $q$. This is the lowest-order anharmonic term  for a mode with the potential that has inversion symmetry, and in this sense, the Duffing model is minimalistic. The major effect of this term for comparatively small vibration amplitudes comes from making the vibration frequency amplitude-dependent \cite{Landau2004a}. For a small decay rate, the frequency change due to this dependence can exceed the frequency uncertainty due to the decay, making the nonlinearity significant. At the same time, the vibrations remain close to sinusoidal as long as the drive is comparatively weak. The response, in this case,  is often analyzed using the Bogoliubov-Krylov averaging method \cite{Bogoliubov1961}, which in this context is equivalent to the rotating wave approximation (RWA) of quantum optics \cite{Mandel1995}.

Flexural modes, which are most frequently studied in nanomechanics, do not necessarily have inversion symmetry. Such symmetry implies that the nanoresonator lies in a plane and the vibrations occur transverse to this plane. Typically, nanoresonators like nanobeams, nanomembranes, or carbon nanotubes, are bent because of an applied gate voltage \cite{Eichler2013} or the asymmetry of the clamping \cite{Schmid2016} or, as in the system  studied here, because the asymmetry imposed by the dielectric transduction electrodes~\cite{Unterreithmeier2009}, see Fig.~\ref{fig:f1}.  For asymmetric modes, along with the term $\propto q^4$ in the mode potential energy, it is necessary to take into account the lower-order term $\propto q^3$. However, in the standard analysis based on the RWA, the effect of this  term on the vibrations at the drive frequency comes to renormalizing the Duffing parameter \cite{Landau2004a}, such that the standard Duffing model remains applicable with an effective Duffing parameter.

Here we demonstrate that, for a  nanoresonator with a broken symmetry, the resonant response can significantly deviate from the standard Duffing response already for a moderately strong driving. In the studied system this happens for the vibration amplitudes where the mode frequency differs from its zero-amplitude value $\omega_0$ by $\lesssim 10^{-4}\omega_0$.

The physics of the effect can be understood from the following argument. In our system the Duffing nonlinearity is hardening,  the frequency increases with the vibration amplitude. On the other hand, the RWA-change of the effective Duffing parameter due to the broken symmetry is negative \cite{Landau2004a}. This means that the symmetry breaking  results in a decrease of the Duffing parameter compared to its value in the symmetric system. For a weak driving it is the decreased effective value that determines the resonant response. However, it is clear that, for sufficiently large amplitudes, the quartic term in the  potential becomes ``stronger'' than the cubic term. Therefore for such amplitudes, the frequency dependence on the amplitude should be different from the small-amplitude range. Remarkably, this happens where the amplitude is still small compared to the scale where the change of the vibration frequency becomes comparable to $\omega_0$.

Another effect of a broken symmetry is the onset of vibrations at the  even multiples of the frequency of the resonant drive (whereas odd multiples arise from the regular Duffing model), and in particular at twice the drive frequency \footnote{The \addMD{onset of vibrations} \sout{peak} at twice the drive frequency has a counterpart at zero frequency, which was studied for nanoresonators in Ref.~\cite{Eichler2013}.}.  The occurrence of vibrations at twice the drive frequency is often referred to as the second harmonic generation. Such vibrations, which are also referred to as temporal harmonics or overtones, were seen earlier in microscale vibrational systems, cf. Ref.~\cite{Asadi2021} and references therein. Here we measure the amplitude and phase of these vibrations directly, and by comparing them to the amplitude and phase of the main tone establish that they are indeed due to the nonlinearity of the mode potential. This allows quantifying this potential experimentally.

To describe the observations, the theoretical analysis should go beyond the standard RWA. A simplifying factor is the small decay rate of the nanoresonator studied in the experiment. This suggests extending the methods of the Hamiltonian  nonlinear dynamics to the problem at hand \cite{Arnold1989}. Such an extension should allow for both the nonlinearity and the weak damping. A theory should also address the observation that, even though the signal at twice the drive frequency has an appreciable amplitude in the studied range, higher-order overtones remain very small.

Below in Sec.~II we discuss the setup of the experiment and in Sec.~III summarize the experimental observations. In Sec.~IV we outline the theory. Section~V provides a discussion of the results and a comparison between the theory and the experiment. Section~VI contains concluding remarks. The Appendices describe auxiliary experimental observations, including the observed dependence of the nonlinearity parameters on the applied DC control voltage, a discussion of potential other mechanisms leading to a vibration  at twice the drive frequency, and further details of the theory.

\section{Setup and Characterization}
We investigate a nanomechanical doubly clamped string resonator fabricated from pre-stressed silicon nitride on a fused silica substrate, similar to the one depicted in Fig.~\ref{fig:f1}.

\begin{figure}[t!]
	\centering
	\includegraphics[width=0.7\linewidth]{F1.pdf}
	\caption{ Scanning electron micrograph of the doubly clamped silicon nitride sting resonator (green) and two adjacent gold electrodes (yellow) for dielectric drive and detection. Depicted is one of the clamping pads (right), the onset of the string and the two control electrodes.}
	\label{fig:f1}
\end{figure}
The string is 270\,nm wide, 100\, nm thick and 55\,$\mu$m long.  It is flanked by two adjacent gold electrodes, enabling the dielectric transduction combined with a microwave cavity-enhanced heterodyne detection scheme, discussed in more detail in \cite{Unterreithmeier_2009,Faust_2012,Rieger_2012}.
The two gold electrodes are placed asymmetrically with respect to the string, leading to an inhomogeneous electrical field when a DC voltage is applied between the electrodes. As a consequence, the dielectric string resonator gets polarized and experiences a gradient force which displaces it from its original equilibrium position as it is getting pulled towards the electrodes where the field is strongest. The corresponding change in the field gradient alters the eigenfrequency, enabling frequency tuning with the applied DC voltage \cite{Rieger_2012}. Concurrently, it also leads to the breaking of the symmetry of the restoring potential (see discussion of Sec.~\ref{sec:theory}).

For the measurements shown in the following, the DC voltage is fixed to $5$\,V, such that the fundamental flexural out-of-plane mode, which will be considered in the following, is well-separated from the corresponding in-plane mode.
The experiment is performed at room temperature of $293$\,K and under vacuum at a pressure below $10^{-4}$\,mbar.

The response of the fundamental out-of-plane mode (referred to as mode M in the following) to a comparatively weak resonant field is described by the Duffing model. In this model the equation of motion for the mode coordinate $q(t)$ reads
\begin{align}
	\label{eq:Duffing}
	\ddot{q} +2 \Gamma \dot{q} + \omega_0^2 q + \gamma_\mathrm{eff} q^3 = F_{d} \cos(\omega_{d} t)  \, .
\end{align}
Here, $\omega_0=2\pi f_{0}$ is the angular mode eigenfrequency, $\Gamma$ is the damping rate, and $\gamma_\mathrm{eff}$ is the effective Duffing nonlinearity parameter. The parameter $\gamma_\mathrm{eff}$ describes the resonant response in the range of the driving amplitude $F_d$ where the RWA applies, i.e., this is the Duffing parameter renormalized by the cubic nonlinearity of the potential. The driving amplitude $F_d$ in Eq.~(\ref{eq:Duffing})  is scaled by the mass. The driving angular frequency $\omega_d=2\pi f_{d}$ is assumed to be close to $\omega_0$, with $|\omega_d - \omega_0|\ll \omega_0$.

In the experiment, both the drive tone and the measured signal are voltage signals. Therefore we calibrate the system in units of volts, as discussed in detail in the Supplemental Material of Ref.~\onlinecite{Huber2020}, and apply the driving amplitude as an RF voltage $V_d$.
For $V_d < 1$\,mV the mode dynamics is linear. The spectrum of the linear response in this range,  including the Lorentzian fit, is shown in Fig.~\ref{fig:af1} in the Appendix. We find $f_0=6.528$\,MHz and $2\Gamma/(2\pi)$=20\,Hz, which gives the quality factor $Q\approx 325 000$.

\section{Experimental Observations}

When the drive voltage is increased but remains close to the linear regime, the measured response is well described by the solution of Eq.~(\ref{eq:Duffing}). The corresponding bidirectional scan at a drive voltage of $9$\,mV is shown in Fig.\,\ref{fig:f2}  as gray dots. A fit  (red line) yields the effective nonlinear Duffing parameter. Since the signals are measured in volts, in what follows we use the superscript $(\mathrm{V})$ to indicate that the appropriate nonlinearity parameters obtained by fitting the signals are also in volts. Using the measured value of the signal $q(t)$ in volts, we obtain $\gamma_\mathrm{eff}^{(\rm V)}/(2\pi)^2 = 2.48 \cdot 10^{15}$\,V$^{-2}$ s$^{-2}$.

The measured response for a stronger drive voltage, $V_d=40$\,mV, is shown in Fig.~\ref{fig:f2} by the black dots. The red curve in this figure, on the other hand, shows the response for this voltage calculated using the above value of $\gamma_\mathrm{eff}^{(\rm V)}$. Equation~(\ref{eq:Duffing}) does not fit the data, indicating that the Duffing model no longer applies for $V_d=40$~mV.

\begin{figure}[t!]
	\centering
	\includegraphics[width=0.8\linewidth]{F2.pdf}
	\caption{Duffing response and non-Duffing response. Measured nonlinear response curves at drive voltages of $V_d= 9$\,mV (gray) and $40$\,mV (black). Both traces display unprocessed raw data; the noise level remains below the size of the dots. A fit to the Duffing model at $V_d= 9$\,mV as well as a theory curve using the obtained $\gamma_\mathrm{eff}^{(\rm V)}/(2\pi)^2 = 2.48 \cdot 10^{15}$ V$^{-2}$ s$^{-2}$ for $V_d= 40$\,mV are included (red lines). The deviation of the data from the Duffing model is clearly visible for the $V_d=40$\,mV data. A better agreement is achieved by taking the influence of the cubic nonlinearity of the potential $\beta$ into account (green line). The curve is calculated for $\beta^{(\rm V)}/(2\pi)^2=2.08\cdot 10^{15}$ V$^{-1}$ s$^{-2}$
	and is truncated at a detuning of $2$\,kHz.}
	\label{fig:f2}
\end{figure}

Along with the deviation from the Duffing response curve we have also observed the onset of a signal at twice the drive frequency. This is depicted for a drive voltage $V_d=100$\,mV applied at the eigenfrequency of the mode, i.e., $f_d = f_0$, in Fig.~\ref{fig:f3}. The spectrum clearly shows the forced vibrations at the eigenfrequency of the driven mode (labelled M), along with a pronounced peak at the  overtone frequency (O) $2f_d =2f_0= 13.056$\,MHz. We note that the second spatial eigenmode of out-of-plane vibrations of the nanostring has the frequency $13.2$\,MHz, almost $200$\,kHz above $2f_d$. A frequency response measurement, indicating the frequency separation of the two features, is shown in Fig.~\ref{fig:af2} in the Appendix. With a frequency separation much larger than the damping rate, both features can clearly be distinguished. This in itself allows us to unambiguously associate the signal overtone at $2f_d$ with the (non-sinusoidal) oscillation of the resonantly driven fundamental mode.

The dependence on the drive amplitude $V_d$ of the amplitudes of the fundamental mode as well as the signals at $2f_d$ and $3f_d$ are shown in Fig.\,\ref{fig:f4}\,(a). The data refers to the drive frequency $f_d=f_0$. The drive voltage $V_d$ is swept between $0$ and $400$\,mV.
The three vibration amplitudes are measured simultaneously with a high-frequency lock-in amplifier by using multiple demodulators.
As already seen in the spectral measurements in Fig.~\ref{fig:f3}, the signal at  $2f_d$ is significantly stronger than the signals at $3f_d$ and $4f_d$ (the signal at $4f_d$ is not shown).

\begin{figure}[b!]
	\centering
	\includegraphics[width=0.7\linewidth]{F3}
	\caption{Spectrum of the resonantly driven fundamental mode (M). For a resonant drive at $f_d =f_0=6.528$\,MHz with $V_d=100$\,mV, the  overtone (O) at $13.056$\,MHz is clearly visible. Higher-order overtones are barely discerned in this representation. A constant noise background has been subtracted from the data.}
	\label{fig:f3}
\end{figure}

\begin{figure*}[t!]
	\centering
	\includegraphics[width=0.8\linewidth]{F4}
	\caption{(a) Drive voltage dependence of the resonant amplitude of the fundamental mode (M, black) as well as the amplitudes of the vibrations at twice (O, blue) and three times the drive frequency (dark blue); the amplitudes of the vibrations at the overtones have been rescaled to account for the different frequency-dependent displacement-to-voltage conversion factor.
	The drive frequency $f_d$ is fixed at the resonance frequency of the fundamental mode $f_0$.
		(b) The amplitude of the vibrations at $2f_d$ [marked by O in panel (a)] depends quadratically on the amplitude of the fundamental mode. Data shown in blue. The red line is a quadratic fit. }
	\label{fig:f4}
\end{figure*}

\section{Theory}
\label{sec:theory}

To account for the broken inversion symmetry, we have to include the cubic nonlinearity in the mode potential $U(q)$ along with the quartic nonlinearity. The equation of motion then reads
\begin{align}
	\label{eq:eom_coordinate}
	\ddot{q} &= - 2 \Gamma \dot{q} -\partial_q U(q)+ F_{d} \cos(\omega_{d} t) \,
\end{align}
where
\begin{align}
	\label{eq:potential}
	U(q) = \frac{1}{2}\omega_0^2q^2  + \frac{1}{3}\beta q^3 + \frac{1}{4}\gamma q^4.
\end{align}
Here $\beta$ and $\gamma$ are the parameters of the cubic and quartic nonlinearity, respectively.
In what follows, we consider a comparatively weak nonlinearity, so that the nonlinear part of the vibration energy remains smaller than the harmonic part $\sim \omega_0^2\langle q^2\rangle$. For concreteness we set $\beta>0$; the sign of $\beta$ can be changed by changing the sign of the coordinate $q$ and incrementing the phase of the drive by $\pi$. We note that in the engineering literature Eq.~(\ref{eq:eom_coordinate}) in the absence of the driving is sometimes called the Helmholtz-Duffing equation. In the context of elastic cables a numerical analysis of Eq.~(\ref{eq:eom_coordinate}) was done  in Ref.~\onlinecite{Benedettini1987}.

For a resonantly driven nonlinear oscillator, a major consequence of the broken symmetry is the occurrence of vibrations at even multiples of the drive frequency beyond the odd multiples expected for a Duffing oscillator. For a weak nonlinearity, forced vibrations are almost sinusoidal,  $q(t) \simeq  q^{(1)}(t) $ with $q^{(1)}(t)  = A\cos(\omega_d t +\varphi)$. To the leading order, the cubic nonlinearity $\beta$  leads to the onset of vibrations at $2{}f_d$. They are described by the expression
\begin{align}
	\label{eq:second_harmonic}
	q^{(2)}(t) =& \frac{1}{6\omega_0^2}\beta A^2 \cos(2{}\omega_d t + 2 \varphi).
\end{align}
We use that  $|\omega_d-\omega_0|\ll \omega_0$ and $\Gamma \ll \omega_0$. An observation of such vibrations is an unambiguous signature of the broken inversion symmetry of the vibrational mode. There are, however, other contributions to the signal at  $2{}f_d$, which are also related to the broken symmetry, as discussed in Appendix~\ref{sec:second_harmonic}.

Besides leading to the onset of vibrations at the even multiples of the drive frequency, the cubic nonlinearity modifies the dependence of the amplitude of forced vibrations on the drive amplitude compared to the Duffing response. The Duffing model (\ref{eq:Duffing}) has been very successful in describing many observations in  nanomechanical systems, and as mentioned in the Introduction, in the majority of cases the analysis was based on the RWA. In the RWA, one changes from the fast oscillating coordinate $q(t)$ and momentum $p(t)$ to slowly varying in-phase and quadrature components, $q(t)  - i\omega_d^{-1}p(t) = [Q(t) - iP(t)] \, \exp(i\omega_dt)$. In the equations for $Q,P$ one then disregards the terms  that oscillate at the frequency $\omega_d$ and its overtones. Then the major effect of the Duffing nonlinearity is the dependence of the mode frequency on the vibration amplitude $A$ \cite{Landau2004a}
\begin{align}
	\label{eq:frequency_renormalized}
	\omega\to\omega_\mathrm{eff} \approx \omega_0 + \frac{3\gamma}{8\omega_0} \, A^2 \, ,
\end{align}
with $A$ given by the value of $(Q^2 + P^2)^{1/2}$ in the stable vibrational state.

The RWA is often applied also to a vibrational system with additional cubic nonlinearity $\beta$. In this approximation the response to the resonant field is mapped onto that of the Duffing model (\ref{eq:Duffing}) with the renormalized nonlinearity parameter $\gamma_\mathrm{eff}$ replacing the bare Duffing parameter $\gamma$ \cite{Landau2004a}
\begin{align}
	\label{eq:renormalized_gamma}
	\gamma\to \gamma_\mathrm{eff} =\gamma-\frac{10\beta^2}{9\omega_0^2}\,.
\end{align}

It is seen from Eq.~(\ref{eq:renormalized_gamma}) that the cubic nonlinearity can strongly affect the amplitude dependence of the mode frequency (\ref{eq:frequency_renormalized}). Indeed, if $\gamma>0$, but $\gamma_\mathrm{eff}<0$, even the sign of $d\omega_\mathrm{eff}/dA^2$ changes. This leads to the so-called zero-dispersion behavior \cite{Dykman1990d} (see Ref.~\cite{Soskin2003} for a comprehensive review). In what follows we consider the case $\gamma$ and $\gamma_\mathrm{eff}>0$, which is relevant for the experiment described in this paper.

The calculation outlined below shows, however, that the strong change of $d\omega_\mathrm{eff}/dA^2$ occurs only in the region of comparatively small amplitudes $A$. Indeed, for large amplitudes the term $\propto q^4$ in $U(q)$ becomes more important than the term $\propto q^3$.
Simple dimensional arguments show that for amplitudes $ A^2 \gtrsim \omega_0^2\gamma_\mathrm{eff}/\gamma^2$, the RWA approximation (\ref{eq:renormalized_gamma}) becomes inapplicable (see Appendix~\ref{app:theory}). For small $\gamma_\mathrm{eff}/\gamma$ this happens where the nonlinear part of the energy $\sim \gamma A^4$ is still small compared to the harmonic part $\sim \omega_0^2A^2$.

The analysis of the resonant  response beyond the renormalization~(\ref{eq:renormalized_gamma}) is significantly simplified in the case of weak damping, where the decay rate $\Gamma\ll \omega_0$. Here it is convenient to change from the coordinate and momentum of the mode to its action-angle variables. This is a canonical transformation. For an isolated mode with the Hamiltonian
\begin{align}
	\label{eq:Hamiltonian}
	H=\frac{1}{2}p^2 + U(q)
\end{align}
the action  $I$ and the angle (phase) $\phi$ are defined as $I=(2\pi)^{-1}\oint p\,dq $ and $\phi = \partial_I  \int p\,dq $ \cite{Landau2004a}. The natural vibration frequency of the mode is
\[\omega(I)=(\partial I/\partial E)^{-1},\]
where $E$ is the mode energy. The coordinate and momentum are functions of $I$ and $\phi$ and are periodic in $\phi$,
\[q(I, \phi+2\pi) = q(I,\phi), \quad p(I, \phi+2\pi) = p(I,\phi).\]

In the presence of friction and dissipation, the second-order equation of motion (\ref{eq:eom_coordinate}) becomes a set of two first-order equations for $I$ and $\phi$,
\begin{align}
	\label{eq:eom_I_phi_full}
	&\dot I = R \, \partial_\phi q , \quad   \dot\phi = \omega(I) - R\, \partial_Iq \nonumber\\
	&R=-2\Gamma p + F_d\cos\omega_dt \, .
\end{align}
In the stationary regime forced vibrations occur at the drive frequency $\omega_d$. This means that $\dot\phi \approx \omega_d $, if we neglect  terms oscillating at $\omega_d$. The action $I$ in this regime has a time-independent component and components oscillating at $\omega_d$. However, as seen from Eq.~(\ref{eq:eom_I_phi_full}), keeping oscillating terms in the right-hand sides of the equations for $\dot I$ and $\dot \phi$  lead to small corrections to $I$ and $\phi$ for a comparatively weak drive; in particular, the corrections to $I$ are $\sim \Gamma I/\omega_d$ and $~|F_d \partial_\phi q  \, /\omega_d|$. The smallness of these corrections is the condition of the applicability of the analysis.

If we disregard the fast-oscillating corrections, the right-hand sides of the equations for $\dot I$ and $\dot\phi$ can be averaged over the period $2\pi/\omega_d$.  In this approximation the equations for the stationary states  read
\begin{align}
	\label{eq:working}
	d\overline{I}/dt = \overline{R\, \partial_\phi q}=0, \quad d\overline{\phi}/dt = \omega(\overline{I}) -\overline{R\, \partial_Iq}=\omega_d
\end{align}
where the overbar implies period averaging. It is clear from Eq.~(\ref{eq:working}) in particular that only the component of $q(I,\phi)$ that oscillates as $\cos\phi$, i.e., the main tone, contributes to the terms that multiply $F_d\cos\omega_dt$.

Equation~(\ref{eq:working}) gives two parameters of the stationary vibrational state, the action $\overline{I}=I_\mathrm{st}$ and the time-independent part $\varphi\equiv \varphi_\mathrm{st}$ of the phase, $\overline{\phi(t)} = \omega_d t + \varphi_\mathrm{st}$. The solution of Eq.~(\ref{eq:working}) is provided in Appendix \ref{app:theory}. The calculation is significantly simplified by the fact that the functions $q(I,\phi), p(I,\phi)$ can be expressed in terms of the Jacobi elliptic functions. This property allows one to find the frequency $\omega(I)$ as well as the amplitudes of  vibrations at $\omega(I)$ and its overtones as functions of $I$ in terms of the elliptic integrals. Inversely, for not too strong nonlinearity, it allows one to express the action $I$ in terms of the amplitude $A$ of the main tone, i.e., of the vibrations at frequency $\omega(I)$.

In Fig.~\ref{fig:f5}\,(a) we plot the frequency $\omega(I)$ vs the square of the amplitude $A$ scaled by  the typical displacement $(2\omega_0^2/\gamma)^{1/2}$ at which the Duffing nonlinearity becomes pronounced. The plots refer to different values of the scaled cubic nonlinearity  $\beta/\omega_0\sqrt{\gamma}$. For $\beta=0$ (brown line) the frequency is linear in $A^2$ for small amplitudes.
In contrast, for the critical value $\beta_\mathrm{cr}/\omega_0  \sqrt{\gamma} = \sqrt{0.9}$ (light green line), where the effective Duffing parameter $\gamma_\mathrm{eff}=0$, the frequency is parabolic in $A^2$ for small $A^2$. For intermediate $0 <\beta < \beta_\mathrm{cr}$ (petrol line) the frequency displays a significant curvature as function of $A^2$ in the small-$A^2$ range.
The curve with the experimental value of $\beta$ obtained in Section~\ref{results} (dark green line) is very close to the line for the critical $\beta$, as the effective Duffing parameter is much smaller than the bare Duffing parameter.
A zoom into the small amplitude regime is shown in Fig.~\ref{fig:f5}\,(b) where
the experimental range of the solution for the experimental value of $\beta$ (dark green line) is
compared to the solution of the renormalized Duffing model, Eq.~(\ref{eq:renormalized_gamma}). The two curves clearly disagree.

We  note that a simple way to think of the cubic term in the potential of a nanoresonator mode is to relate it to a linear bias. Such bias can come, for example, from a gate voltage that ``pulls'' the nanoresonator. The  potential of the Duffing oscillator with linear bias is
\begin{align}
	\label{eq:bias}
	U_B(q) = -\xi_\mathrm{bias} q + \frac{1}{2}\omega_0^2 q^2  + \frac{1}{4}\gamma q^4,
\end{align}
where $\xi_\mathrm{bias}$ is the bias strength.

The potential $U_B$ has the same form as  the potential~(\ref{eq:potential}) for $|\beta|<(4\gamma\omega_0^2)^{1/2}$. This is seen if one shifts the equilibrium position $q\to q+\delta q$ to compensate the linear term; in the limit of small nonlinearity $\delta q \approx \xi_\mathrm{bias}/\omega_0^2$. After the shift $U_B(q)$ becomes of the same form as $U(q)$ with
\[\beta = 3\gamma \delta q, \quad \omega_0^2 \to \omega_0^2 + 3\gamma (\delta q)^2.\]
The condition  $\beta^2 < 4 \gamma \omega_0^2$ corresponds to the experimental situation where the potential $U(q)$ has a single minimum.

	\begin{figure*}[t!]
		\centering
		\includegraphics[width=0.8\linewidth]{F5}

				\caption  {Relative variation of the frequency $\omega(I)$ vs the square of the amplitude $A$ scaled by $q_0 = \sqrt{2 \omega_0^2 /\gamma}$.
				(a) The brown, petrol, and light green curves refer, respectively, to $\beta /\omega_0\sqrt{\gamma}= 0, 0.6$, and $\sqrt{0.9}$ (this is the critical value for which $\gamma_\mathrm{eff}=0$). The dark green curve is computed using the value of  $\beta /\omega_0\sqrt{\gamma}$ used in the comparison with the experiment in Fig.~\ref{fig:f2}.
				(b) Closeup showing  the relative variation of $\omega(I)$ using the same $\beta$ as in Fig.~\ref{fig:f2} (dark green) compared with the  Duffing model for the experimental value of $\gamma_\mathrm{eff}$ (gray) for the upper vibrational branch.
						}
		\label{fig:f5}
	\end{figure*}

\section{Results}\label{results}

The response curve shown in Fig.~\ref{fig:f2} is anomalous in the sense that it significantly differs from the conventional Duffing curve. The theoretical model of Sec.~\ref{sec:theory}, which takes into account the cubic nonlinearity of the potential, allows us to describe this curve quantitatively and to find the nonlinearity parameter $\beta^{(\rm V)}$. The analysis is based on Eqs.~(\ref{eq:working}) and uses the experimentally determined eigenfrequency $\omega_0$ and linewidth $2\Gamma$ (see  Appendix \ref{app:theory}).

The effective Duffing parameter $\gamma_\mathrm{eff}^{(\rm V)}$ (in volts) is extracted from the measurements for a moderately weak driving, where the effective Duffing model applies, cf. the data for $V_d=9$~mV in Fig.~\ref{fig:f2}. The cubic nonlinearity parameter $\beta^{(\rm V)}$ is then chosen to match the theoretical result to the experimental curve for a stronger drive. The underlying theory goes beyond the conventional RWA approximation.  The green line in Fig.~\ref{fig:f2} displays the theoretical curve for $\beta^{(\rm V)}/(2\pi)^2=2.08 \cdot 10^{15}$\,V$^{-1}$s$^{-2}$. With this parameter value, we find good agreement between the theory and the experiment.

However, it is seen from Fig.~\ref{fig:f2} that the upper branch of the experimental response curve for $V_d=40$~mV ends at a detuning of approximately $1.5$\,kHz, whereas the theoretical branch extends much further and is truncated at a detuning of $2$\,kHz. We attribute this discrepancy to a comparatively short lifetime of the large-amplitude state. Our data acquisition system does not allow us to observe state with a short lifetime, and therefore we do not observe the stable large-amplitude state beyond a certain frequency.

We now comment on the lifetime of the large-amplitude state.  As seen from Fig.~\ref{fig:f2}, the theoretical values of the amplitudes of this state and the unstable state  are very close. In the considered very weakly damped system these amplitudes are determined by the quasienergies (Floquet eigenvalues) of the periodically driven mode in the corresponding states, which indicates that the corresponding quasienergies are also close. Thermal noise, which is invariably present in the system, leads to escape from a dynamically stable vibrational state \cite{Dykman1979b}. In the system investigated here such escape was seen earlier \cite{Huber2020}.

In the weak-damping regime, periodically driven systems escape via diffusion over quasienergy. Generically, the escape rate increases exponentially with the decreasing distance between the quasienergies of the stable and unstable (saddle-type) states \cite{Dykman2005d}. Therefore we expect it to be comparatively large where the amplitude of the stable state is close to that of the unstable state. A full calculation of the escape rate is beyond the scope of the present paper.

The values of $\gamma_\mathrm{eff}^{\rm (V)}$ and $\beta^{\rm (V)}$ allow us to find the ``bare'' Duffing parameter $\gamma^{(\rm V)}$ using Eq.~(\ref{eq:renormalized_gamma}). The value of this parameter  $\gamma^{(\rm V)}/(2\pi)^2 \approx 1.16 \cdot 10^{17} \, \rm V^{-2} \rm s^{-2} $  is two orders of magnitude larger than the effective Duffing coefficient $\gamma_\mathrm{eff}^{(\rm V)}/(2\pi)^2 = 2.48 \cdot 10^{15}$ V$^{-2}$s$^{-2}$ measured for the moderately weak driving regime.

The signal observed in the experiment at twice the drive frequency arises from the broken inversion symmetry of the nanostring under investigation. The symmetry breaking is caused primarily by the asymmetric arrangement of the dielectric transduction electrodes (Fig.~\ref{fig:f1}). However, besides the mechanism described in the previous section and associated with the cubic nonlinearity of the potential $\propto \beta$, other mechanisms originating from the broken inversion symmetry could also contribute. They include the nonlinear excitation of the second spatial harmonic eigenmode of the nanostring and the nonlinear coupling to the driving force resulting in a direct as well as a parametric drive at $2f_d$. These mechanisms are discussed in more detail in Appendix~\ref{sec:second_harmonic}. Notably, some of them lead to a different dependence of the signal at $2f_d$ on the amplitude and phase of the fundamental mode M. This can be used to identify the origin of the signal. In the following, we compare our experimental observations with the theoretical predictions for all these mechanisms and indeed find that the nonlinearity of the potential energy of the nanoresonator characterized by the parameter $\beta$ is the dominant source of the signal at $2f_d$.

Figure\,\ref{fig:f4}\,(b) displays the amplitude of the  overtone signal (in volts) $V_{\rm out,O}$ at $2{}f_d$ as a function of the amplitude $V_{\rm out, M}$ of the resonant response of the mode M. A fit with $V_{\rm out,O}=cV^2_{\rm out, M}$ clearly shows that the overtone amplitude scales quadratically with that of the main tone. In addition, we have determined the phase of the signal measured at $2{}f_d$ with respect to the phase $\varphi$ of the  mode M with the lock-in amplifier. The measurement indeed reveals that the phase of the overtone coincides with $2\varphi$ for a comparatively weak driving  used to derive Eq.~(\ref{eq:second_harmonic}). Both observations are in agreement with Eq.~(\ref{eq:second_harmonic}).
This suggests that the nonlinear driving terms (\ref{eq:2nd_overtone_direct}) and (\ref{eq:nonlinear_suscept}) discussed in Appendix \ref{sec:second_harmonic} make a small contribution to the  signal  at $2{}f_d$ at most, since they exhibit a different dependence on the amplitude and phase of the main tone.

However, this is not sufficient to establish the nonlinearity of the mode potential $U(q)$  as the only source of the signal at $2{}f_d$, as the experimental observations are also compatible with the nonlinear coupling to the second spatial harmonic  of the nanostring, see  Eq.~(\ref{eq:mode_2_vibrations}). In order to estimate the role of this mechanism as the remaining alternative source of the overtone signal, we repeat the experiment for a drive frequency at half the eigenfrequency of the fundamental mode $\omega_0/2$. In this regime, according to Eq.~(\ref{eq:eom_coordinate}),
the force $F_d\cos(\omega_dt)$ still excites forced  vibrations with a displacement $q^{(0)}(t)\approx [F_d/(\omega_0^2-\omega_d{}^2)]\cos\omega_d t$, even when nonresonant as it is the case for $\omega_d \approx \omega_0/2$. Because of the nonlinearity of the mode, which is characterized by the parameter $\beta$ in Eq.~(\ref{eq:potential}), these vibrations resonantly excite vibrations at $2{}f_d$, which is the effect of resonant second harmonic generation,
	%
	\begin{align}
		\label{eq:half_eigenfrequency}

		q^{(2)}(t) \approx \frac{\beta F_d{}^2}{2(\omega_0^2-\omega_d{}^2)^2}

		\times\mathrm{Re}\frac{\exp(2i\omega_d t)}{\omega_0^2 - 4\omega_d{}^2 +4i\Gamma\omega_d}, \nonumber \\ \quad \omega_d\approx \omega_0/2.
	\end{align}
	%

As seen from Fig.~\ref{fig:f6}, we clearly observe the corresponding overtone at $2f_d\approx f_0$, whereas the response at the drive frequency $f_d$ could not be resolved for bandwidth limitations of the experimental setup.
Figure~\ref{fig:f6}\,(a) displays the overtone amplitude for a sweep of the drive frequency around $f_0/2$ along with a Lorentzian fit. The fit is in full agreement with Eq.~(\ref{eq:half_eigenfrequency}). Figure~\ref{fig:f6}\,(b) shows the scaling of the overtone amplitude for $f_d=f_0/2$ with the drive voltage fitted with a quadratic function.

The response of the second spatial harmonic of the nanostring is nonresonant for $f_d\approx f_0/2$ and may not lead to an appreciable signal. This demonstrates that the cubic nonlinearity of the potential $U(q)$  of the mode M is the major contributor to the response at twice the drive frequency and suggests that it is also a major contributor in the case of the driving at $f_d\approx f_0$. In Appendix \ref{sec:second_harmonic} we provide an argument why the nonlinear coupling to the second spatial harmonic of the nanostring should be weak in addition to being just nonlinear.

\begin{figure*}[t!]
	\centering
	\includegraphics[width=0.8\linewidth]{F6}
	\caption  {Resonant excitation of the overtone  for a drive at half the eigenfrequency. (a) Amplitude of the resonantly excited vibrations at $2f_d$ as a function of the drive frequency $f_d$ swept around $f_0/2$ with $V_d=100$\,mV. The red line displays a Lorentzian fit with halfwidth $\Gamma$. (b) Peak amplitude of the resonantly excited vibrations at $f_0$ for $f_d=f_0/2$ as a function of the drive voltage. The red line correponds to a quadratic fit.
	}
	\label{fig:f6}
\end{figure*}

The finding that the observed signal at $2f_d \approx 2f_0$ is predominantly caused by the nonlinearity $\beta$ suggests that the ratio $c=V_{\rm out,O}/V^2_{\rm out, M}$, which can be extracted from the quadratic fit to the data in Fig.~\ref{fig:f4}\,(b), could be employed to quantitatively determine the cubic nonlinearity parameter in units of volts, $\beta^{(\rm V)}$, from Eq. (\ref{eq:second_harmonic}). However, the displacement-to-voltage conversion factor of our read-out apparatus at $2f_0$ does not coincide with the calibrated one at $f_0$ \cite{Gajo_2020}, such that $\beta^{(\rm V)}$ cannot be directly extracted from the data.

\section{Conclusions}
Conventionally, the nonlinear response of vibrational modes to moderately strong resonant driving is described by the Duffing model. In contrast, in this work we have observed and explained a non-Duffing resonant nonlinear response of  the fundamental mode of an underdamped  nanomechanical resonator. We found that, even though the dependence of the amplitude of forced vibrations on the frequency of the drive is of the familiar Duffing form for moderately weak driving, its shape changes significantly as the driving becomes stronger. This happens in the range where the driving still remains not too strong, so that the forced vibrations are still close to sinusoidal. Our explanation is based on taking into account the broken inversion symmetry of the resonator. The symmetry breaking leads to the onset of a term in the potential energy of the mode, which  is cubic in the mode coordinate.

We show that, while the cubic term does not change the form of the response to a moderately weak drive, it significantly changes the response for a stronger drive. The analysis required us to go beyond the standard rotating-wave approximation. The obtained results allow describing the spectrum of the nonlinear response in a broad range of amplitudes of the resonant drive, where the shape of the response curve changes significantly. The comparison between the experimental data and the theoretical model allows us to determine the parameters of the quartic and cubic nonlinearity of the potential of the mode.

Along with the change of the response form, the broken inversion symmetry leads to the onset of response at even multiples of the drive frequency. For the drive at frequency close to the mode eigenfrequency, we have observed vibrations at twice the drive frequency, i.e., second harmonic generation, in optics terms. We have discussed several microscopic mechanisms that lead to the onset of such vibrations. The dependence of the vibration amplitude and phase on the amplitude and phase of the vibrations at the drive frequency suggests that the major contribution to the frequency doubling comes from the cubic nonlinearity of the mode potential. Characteristically, higher overtones have a very small amplitude, consistent with the model.

We have also observed resonant second harmonic generation when the driving frequency was close to half the mode eigenfrequency. The spectrum of the vibrations at twice the drive frequency has a characteristic Lorentzian shape, while the vibration amplitude is quadratic in the driving amplitude.

Our results demonstrate that weakly damped vibrations of nanomechanical systems display very rich nonlinear dynamics.  It comes from the interplay and competition of different nonlinearity parameters even in simple cases where the vibrations remain close to sinusoidal.  The significant compensation of the nonlinearity that we have found in a certain range of the vibration amplitudes, which is controlled by how strongly the inversion symmetry is broken, can be used in applications as it extends the practically important regime of almost linear behavior of the mode.  At the same time, the deviation from the traditionally assumed Duffing behavior should be generic for weakly damped nanomechanical modes, as in many cases such modes lack inversion symmetry.

\section{Acknowledgements}
We are grateful to H. Yamaguchi for the discussion of the mechanisms of overtone generation. J.\,S.\,O. and E.\,M.\,W. gratefully acknowledge financial support from the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) through Project-ID 425217212 - SFB 1432, the European Union’s Horizon 2020 Research and Innovation Programme under Grant
Agreement No 732894 (FET Proactive HOT), and the German Federal Ministry of Education and Research
(contract no. 13N14777) within the European QuantERA cofund project QuaSeRT.
M.\,I.\,D. acknowledges support from the National Science Foundation, Grants No. DMR-1806473 and CMMI 1661618. M.\,I.\,D. is a senior fellow of the Zukunftskolleg of the University of Konstanz; he is grateful for the warm hospitality at the University of Konstanz where this work was started.