\section{Introduction}

\begin{figure}[t!]
    \centering
    \includegraphics[width=1.0\columnwidth]{{6371.problem_of_previous_few_shot}.pdf}
    \caption{
    Examples of identity preservation failures and improved results generated by the proposed method. Each row shows (a) driver shape interference, (b) losing details of target identity, and (c) failure of warping at large poses.
    }
    \label{fig:problem_of_previous_few_shot}
 \end{figure}

\begin{figure*}[t]
    \includegraphics[width=0.9\linewidth]{{6371.architecture}.pdf}
    \centering
    \caption{
    The overall architecture of MarioNETte.
    }
    \label{fig:marionette_architecture}
 \end{figure*}

Given a \textit{target} face and a \textit{driver} face, face reenactment aims to synthesize a \textit{reenacted} face which is animated by the movement of a driver while preserving the identity of the target.

Many approaches make use of generative adversarial networks (GAN) which have demonstrated a great success in image generation tasks.
\citet{xu-arxiv-2017-facetransfer,wu-eccv-2018-reenactgan} achieved high-fidelity face reenactment results by exploiting CycleGAN~\cite{zhu-cvpr-2017-cyclegan}.
However, the CycleGAN-based approaches require at least a few minutes of training data for each target and can only reenact predefined identities, which is less attractive in-the-wild where a reenactment of unseen targets cannot be avoided.

The few-shot face reenactment approaches, therefore, try to reenact any unseen targets by utilizing operations such as adaptive instance normalization (AdaIN)~\cite{zakharov-arxiv-2019-samsung} or warping module~\cite{wiles-eccv-2018-x2face,siarohin-cvpr-2019-monkeynet}.
However, current state-of-the-art methods suffer from the problem we call \textit{identity preservation problem}: the inability to preserve the identity of the target leading to defective reenactments.
As the identity of the driver diverges from that of the target, the problem is exacerbated even further.

Examples of flawed and successful face reenactments, generated by previous approaches and the proposed model, respectively, are illustrated in Figure~\ref{fig:problem_of_previous_few_shot}.
The failures of previous approaches, for the most part, can be broken down into three different modes~\footnote{Additional example images and videos can be found at the following URL: \url{http://hyperconnect.github.io/MarioNETte}}:

\begin{enumerate}
    \item Neglecting the identity mismatch may lead to a identity of the driver interfere with the face synthesis such that the generated face resembles the driver (Figure~\ref{fig:problem_of_previous_few_shot}a).
    \item Insufficient capacity of the compressed vector representation (e.g., AdaIN layer) to preserve the information of the target identity may lead the produced face to lose the detailed characteristics (Figure~\ref{fig:problem_of_previous_few_shot}b).
    \item Warping operation incurs a defect when dealing with large poses (Figure~\ref{fig:problem_of_previous_few_shot}c).
\end{enumerate}

We propose a framework called \textit{MarioNETte}, which aims to reenact the face of unseen targets in a few-shot manner while preserving the identity without any fine-tuning.
We adopt \textit{image attention block} and \textit{target feature alignment}, which allow MarioNETte to directly inject features from the target when generating image.
In addition, we propose a novel \textit{landmark transformer} which further mitigates the identity preservation problem by adjusting for the identity mismatch in an unsupervised fashion.
Our contributions are as follows:

\begin{itemize}
    \item We propose a few-shot face reenactment framework called MarioNETte, which preserves the target identity even in situations where the facial characteristics of the driver differs widely from those of the target. Utilizing image attention block, which allows the model to attend to relevant positions of the target feature map, together with target feature alignment, which includes multiple feature-level warping operations, proposed method improves the quality of the face reenactment under different identities.
    \item We introduce a novel method of landmark transformation which copes with varying facial characteristics of different people. The proposed method adapts the landmark of a driver to that of the target in an unsupervised manner, thereby mitigating the identity preservation problem without any additional labeled data.
    \item We compare the state-of-the-art methods when the target and the driver identities coincide and differ using VoxCeleb1~\cite{nagrani-interspeech-2017-voxceleb} and CelebV~\cite{wu-eccv-2018-reenactgan} dataset, respectively. Our experiments including user studies show that the proposed method outperforms the state-of-the-art methods.

\end{itemize}

\section{MarioNETte Architecture}

Figure~\ref{fig:marionette_architecture} illustrates the overall architecture of the proposed model.
A conditional \textit{generator} $G$ generates the reenacted face given the driver $\mathbf{x}$ and the target images $\{\mathbf{y}^i\}_{i=1 \ldots K}$, and the \textit{discriminator} $D$ predicts whether the image is real or not.
The generator consists of following components:

\begin{itemize}
    \item The \textbf{preprocessor} $P$ utilizes a 3D landmark detector~\cite{bulat-iccv-2017-facealignment} to extract facial keypoints and renders them to landmark image, yielding $\mathbf{r}_x = P(\mathbf{x})$ and $\mathbf{r}_{y}^{i} = P(\mathbf{y}^{i})$, corresponding to the driver and the target input respectively. Note that proposed landmark transformer is included in the preprocessor. Since we normalize the scale, translation and rotation of landmarks before using them in a landmark transformer, we utilize 3D landmarks instead of 2D ones.

    \item The \textbf{driver encoder} $E_x(\mathbf{r}_x)$ extracts pose and expression information from the driver input and produces driver feature map $\mathbf{z}_x$.

    \item The \textbf{target encoder} $E_y(\mathbf{y}, \mathbf{r}_y)$ adopts a U-Net architecture to extract style information from the target input and generates target feature map $\mathbf{z}_y$ along with the warped target feature maps $\hat{\mathbf{S}}$.

    \item The \textbf{blender} $B(\mathbf{z}_x, \{\mathbf{z}_{y}^{i}\}_{i=1 \ldots K})$ receives driver feature map $\mathbf{z}_x$ and target feature maps $\mathbf{Z}_y = [\mathbf{z}_y^{1}, \ldots ,\mathbf{z}_y^{K}]$ to produce mixed feature map $\mathbf{z}_{xy}$. Proposed image attention block is basic building block of the blender.

    \item The \textbf{decoder} $Q(\mathbf{z}_{xy}, \{\hat{\mathbf{S}}^{i}\}_{i=1 \ldots K})$ utilizes warped target feature maps $\hat{\mathbf{S}}$ and mixed feature map $\mathbf{z}_{xy}$ to synthesize reenacted image. The decoder improves quality of reenacted image exploiting proposed target feature alignment.
\end{itemize}

\noindent For further details, refer to Supplementary Material A1.

\subsection{Image attention block}

\begin{figure}[t]
    \includegraphics[width=1.0\columnwidth]{{6371.image_attention_block}.pdf}
    \caption{
    Architecture of the image attention block.
    Red boxes conceptually visualize how each position of $\mathbf{z}_x$ and $\mathbf{Z}_y$ are associated. Our attention can attend different position of each target feature maps with different importance.
    }
    \label{fig:image_attention_block}
 \end{figure}

To transfer style information of targets to the driver, previous studies encoded target information as a vector and mixed it with driver feature by concatenation or AdaIN layers~\cite{liu-arxiv-2019-funit,zakharov-arxiv-2019-samsung}.
However, encoding targets as a spatial-agnostic vector leads to losing spatial information of targets.
In addition, these methods are absent of innate design for multiple target images, and thus, summary statistics (e.g. mean or max) are used to deal with multiple targets which might cause losing details of the target.

We suggest image attention block (Figure~\ref{fig:image_attention_block}) to alleviate aforementioned problem.
The proposed attention block is inspired by the encoder-decoder attention of transformer~\cite{vaswani-nips-2017-transformer}, where the driver feature map acts as an attention query and the target feature maps act as attention memory.
The proposed attention block attends to proper positions of each feature (red boxes in Figure~~\ref{fig:image_attention_block}) while handling multiple target feature maps (i.e., $\mathbf{Z}_y$).

Given driver feature map $\mathbf{z}_x \in \mathbb{R}^{h_x \times w_x \times c_x}$ and target feature maps $\mathbf{Z}_y = [\mathbf{z}_y^{1},\ldots,\mathbf{z}_y^{K}] \in \mathbb{R}^{K \times h_y \times w_y \times c_y}$, the attention is calculated as follows:
\begin{equation}
\begin{aligned}
 \mathbf{Q} &= \mathbf{z}_{x} \mathbf{W}_{q} + \mathbf{P}_{x} \mathbf{W}_{qp} & \in &\; \mathbb{R}^{h_x \times w_x \times c_a} \\
 \mathbf{K} &= \mathbf{Z}_{y} \mathbf{W}_{k} + \mathbf{P}_{y} \mathbf{W}_{kp} & \in &\; \mathbb{R}^{K \times h_y \times w_y \times c_a} \\
 \mathbf{V} &= \mathbf{Z}_{y} \mathbf{W}_{v} & \in &\; \mathbb{R}^{K \times h_y \times w_y \times c_x} \\
\end{aligned}
\end{equation}

\begin{equation}
     A(\mathbf{Q}, \mathbf{K}, \mathbf{V}) = \text{softmax}\left(\frac{f(\mathbf{Q}) f(\mathbf{K})^T}{\sqrt{c_a}}\right)f(\mathbf{V}),
\end{equation}

\noindent where $f: \mathbb{R}^{d_1 \times \ldots \times d_k \times c} \xrightarrow{} \mathbb{R}^{(d_1 \times \ldots \times d_k) \times c}$ is a flattening function, all $\mathbf{W}$ are linear projection matrices that map to proper number of channels at the last dimension, and $\mathbf{P}_x$ and $\mathbf{P}_y$ are sinusoidal positional encodings which encode the coordinate of feature maps (further details of sinusoidal positional encodings we used are described in Supplementary Material A2).
Finally, the output $A(\mathbf{Q}, \mathbf{K}, \mathbf{V}) \in \mathbb{R}^{(h_x \times w_x) \times c_x}$ is reshaped to $\mathbb{R}^{h_x \times w_x \times c_x}$.

Instance normalization, residual connection, and convolution layer follow the attention layer to generate output feature map $\mathbf{z}_{xy}$.
The image attention block offers a direct mechanism of transferring information from multiple target images to the pose of driver.

\subsection{Target feature alignment}

\begin{figure}[t]
    \includegraphics[width=1.0\columnwidth]{{6371.style_feature_flow}.pdf}
    \caption{Architecture of target feature alignment.}
    \label{fig:style_feature_flow}
\end{figure}

The fine-grained details of the target identity can be preserved through the warping of low-level features~\cite{siarohin-cvpr-2019-monkeynet}.
Unlike previous approaches that estimate a warping flow map or an affine transform matrix by computing the difference between keypoints of the target and the driver~\cite{balakrishnan-cvpr-2018-synthesizing,siarohin-cvpr-2018-deformablegan,siarohin-cvpr-2019-monkeynet}, we propose a target feature alignment (Figure~\ref{fig:style_feature_flow}) which warps the target feature maps in two stages: (1) \textit{target pose normalization} generates pose normalized target feature maps and (2) \textit{driver pose adaptation} aligns normalized target feature maps to the pose of the driver.
The two-stage process allows the model to better handle the structural disparities of different identities.
The details are as follows:

\begin{enumerate}
    \item \textbf{Target pose normalization.} In the target encoder $E_y$, encoded feature maps $\{\mathbf{S}_{j}\}_{j=1 \ldots n_{y}}$ are processed into $\hat{\mathbf{S}} = \{\mathcal{T}(\mathbf{S}_{1}; \mathbf{f}_{y}), \ldots ,\mathcal{T}(\mathbf{S}_{n_y}; \mathbf{f}_{y})\}$ by estimated normalization flow map $\mathbf{f}_{y}$ of target and warping function $\mathcal{T}$ (\textcircled{\raisebox{-0.9pt}{1}} in Figure~\ref{fig:style_feature_flow}).
    The following \textit{warp-alignment block} at decoder treats $\hat{\mathbf{S}}$ in a target pose-agnostic manner.

    \item \textbf{Driver pose adaptation.} The warp-alignment block in the decoder receives $\{\hat{\mathbf{S}}^i\}_{i=1 \ldots K}$ and the output $\mathbf{u}$ of the previous block of the decoder.
    In a few-shot setting, we average resolution-compatible feature maps from different target images (i.e., $\hat{\mathbf{S}}_j = \sum_i\hat{\mathbf{S}}^i_j/K$).
    To adapt pose-normalized feature maps to the pose of the driver, we generate an estimated flow map of the driver $\mathbf{f}_u$ using $1 \times 1$ convolution that takes $\mathbf{u}$ as the input.
    Alignment by $\mathcal{T}(\hat{\mathbf{S}}_j;\mathbf{f}_{u})$ follows (\textcircled{\raisebox{-0.9pt}{2}} in Figure~\ref{fig:style_feature_flow}). Then, the result is concatenated to $\mathbf{u}$ and fed into the following residual upsampling block.
\end{enumerate}

\section{Landmark Transformer}

Large structural differences between two facial landmarks may lead to severe degradation of the quality of the reenactment.
The usual approach to such a problem has been to learn a transformation for every identity~\cite{wu-eccv-2018-reenactgan} or by preparing a paired landmark data with the same expressions~\cite{zhang-arxiv-2019-faceswapnet}.
However, these methods are unnatural in a few-shot setting where we handle unseen identities, and moreover, the labeled data is hard to be acquired.
To overcome this difficulty, we propose a novel \textit{landmark transformer} which transfers the facial expression of the driver to an arbitrary target identity.
The landmark transformer utilizes multiple videos of unlabeled human faces and is trained in an unsupervised manner.

\subsection{Landmark decomposition}

Given video footages of different identities, we denote $\mathbf{x}(c,t)$ as the $\mathit{t}$-th frame of the $\mathit{c}$-th video, and $\mathbf{l}(c,t)$ as a 3D facial landmark.
We first transform every landmark into a normalized landmark $\mathbf{\bar{l}}(c,t)$ by normalizing the scale, translation, and rotation.
Inspired by 3D morphable models of face~\cite{blanz-siggraph-1999-3dmm}, we assume that normalized landmarks can be decomposed as follows:

\begin{equation}
\begin{aligned}
    \mathbf{\bar{l}}(c,t) &= \mathbf{\bar{l}}_m + \mathbf{\bar{l}}_{id}(c) + \mathbf{\bar{l}}_{exp}(c, t),
\end{aligned}
\label{eq:landmark_decompose}
\end{equation}
where $\mathbf{\bar{l}}_m$ is the average facial landmark geometry computed by taking the mean over all landmarks, $\mathbf{\bar{l}}_{id}(c)$ denotes the landmark geometry of identity $c$, computed by $\mathbf{\bar{l}}_{id}(c) = \sum_{t} \mathbf{\bar{l}}(c,t)/{T_c} - \mathbf{\bar{l}}_m$ where $T_c$ is the number of frames of $c$-th video, and $\mathbf{\bar{l}}_{exp}(c, t)$ corresponds to the expression geometry of $\mathit{t}$-th frame.
The decomposition leads to $\mathbf{\bar{l}}_{exp}(c,t) = \mathbf{\bar{l}}(c,t) - \mathbf{\bar{l}}_m - \mathbf{\bar{l}}_{id}(c)$.

Given a target landmark $\mathbf{\bar{l}}(c_y,t_y)$ and a driver landmark $\mathbf{\bar{l}}(c_x,t_x)$ we wish to generate the following landmark:
\begin{equation}
    \mathbf{\bar{l}}(c_x \xrightarrow{} c_y,t_x) = \mathbf{\bar{l}}_m + \mathbf{\bar{l}}_{id}(c_y) + \mathbf{\bar{l}}_{exp}(c_x, t_x),
\label{eq:landmark_transformation}
\end{equation}
i.e., a landmark with the identity of the target and the expression of the driver.
Computing $\mathbf{\bar{l}}_{id}(c_y)$ and $\mathbf{\bar{l}}_{exp}$ is possible if enough images of $c_y$ are given, but in a few-shot setting, it is difficult to disentangle landmark of unseen identity into two terms.

\subsection{Landmark disentanglement}

To decouple the identity and the expression geometry in a few-shot setting, we introduce a neural network to regress the coefficients for linear bases.
Previously, such an approach has been widely used in modeling complex face geometries~\cite{blanz-siggraph-1999-3dmm}.
We separate expression landmarks into semantic groups of the face (e.g., mouth, nose and eyes) and perform PCA on each group to extract the expression bases from the training data:
\begin{equation}
    \mathbf{\bar{l}}_{exp}(c,t) = \sum_{k=1}^{n_{exp}} \alpha_k(c,t) \mathbf{b}_{exp, k} = \mathbf{b}_{exp}^{T} \bm{\alpha}(c,t),
    \label{eq:expression_regress}
\end{equation}
where $\mathbf{b}_{exp, k}$ and $\alpha_k$ represent the basis and the corresponding coefficient, respectively.

The proposed neural network, a \textit{landmark disentangler} $M$, estimates $\bm{\alpha}(c,t)$ given an image $\mathbf{x}(c,t)$ and a landmark $\mathbf{\bar{l}}(c,t)$.
Figure~\ref{fig:landmark_transform_architecture} illustrates the architecture of the landmark disentangler.
Once the model is trained, the identity and the expression geometry can be computed as follows:
\begin{equation}
\begin{aligned}
    \bm{\hat{\alpha}}(c,t) &= M\left(\mathbf{x}(c,t), \;  \mathbf{\bar{l}}(c,t)\right) \\
    \mathbf{\hat{l}}_{exp}(c,t) &= \lambda_{exp} \mathbf{b}_{exp}^T \bm{\hat{\alpha}}(c,t) \\
    \mathbf{\hat{l}}_{id}(c) &= \mathbf{\bar{l}}(c,t) - \mathbf{\bar{l}}_m - \mathbf{\hat{l}}_{exp}(c,t),
\end{aligned}
\label{eq:estimate_id_and_exp}
\end{equation}
where $\lambda_{exp}$ is a hyperparameter that controls the intensity of the predicted expressions from the network.
Image feature extracted by a ResNet-50 and the landmark, $\mathbf{\bar{l}}(c,t) - \mathbf{\bar{l}}_m$, are fed into a 2-layer MLP to predict $\bm{\hat{\alpha}}(c,t)$.

During the inference, the target and the driver landmarks are processed according to Equation~\ref{eq:estimate_id_and_exp}.
When multiple target images are given, we take the mean value over all $\mathbf{\hat{l}}_{id}(c_y)$.
Finally, landmark transformer converts landmark as:
\begin{equation}
    \mathbf{\hat{l}}(c_x \xrightarrow{} c_y,t_x) = \mathbf{\bar{l}}_m + \mathbf{\hat{l}}_{id}(c_y) + \mathbf{\hat{l}}_{exp}(c_x, t_x).
\end{equation}
Denormalization to recover the original scale, translation, and rotation is followed by the rasterization that generates a landmark adequate for the generator to consume.
Further details of landmark transformer are described in Supplementary Material B.

 \begin{figure}[t]
    \includegraphics[width=1.0\columnwidth]{{6371.marionette_fig_landdis}.pdf}
    \caption{
    Architecture of landmark disentangler.
    Note that $\mathbf{\bar{l}}(c, t)$ is a set of landmark points but visualized as an image in the figure.
    }
    \label{fig:landmark_transform_architecture}
 \end{figure}

\section{Experimental Setup}

\subsubsection{Datasets}

We trained our model and the baselines using VoxCeleb1~\cite{nagrani-interspeech-2017-voxceleb}, which contains $256\times256$ size videos of 1,251 different identities.
We utilized the test split of VoxCeleb1 and CelebV~\cite{wu-eccv-2018-reenactgan} for evaluating self-reenactment and reenactment under a different identity, respectively.
We created the test set by sampling 2,083 image sets from randomly selected 100 videos of VoxCeleb1 test split, and uniformly sampled 2,000 image sets from every identity from CelebV.
The CelebV data includes the videos of five different celebrities of widely varying characteristics, which we utilize to evaluate the performance of the models reenacting unseen targets, similar to in-the-wild scenario.
Further details of the loss function and the training method can be found at Supplementary Material A3 and A4.

\begin{figure*}[t]
    \includegraphics[width=1.0\linewidth]{{6371.main_figure_final}.pdf}
    \centering
    \caption{Images generated by the proposed method and baselines, reenacting different identity on CelebV in one-shot setting.}
    \label{fig:main_result}
 \end{figure*}

\subsubsection{Baselines}

MarioNETte variants, with and without the landmark transformer (\textit{MarioNETte+LT} and \textit{MarioNETte}, respectively), are compared with state-of-the-art models for few-shot face reenactment.
Details of each baseline are as follows:

\begin{itemize}
    \item \textbf{X2Face}~\cite{wiles-eccv-2018-x2face}. X2face utilizes direct image warping. We used the pre-trained model provided by the authors, trained on VoxCeleb1.

    \item \textbf{Monkey-Net}~\cite{siarohin-cvpr-2019-monkeynet}. Monkey-Net adopts feature-level warping. We used the implementation provided by the authors. Due to the structure of the method, Monkey-Net can only receive a single target image.

    \item \textbf{NeuralHead}~\cite{zakharov-arxiv-2019-samsung}.
    NeuralHead exploits AdaIN layers. Since a reference implementation is absent, we made an honest attempt to reproduce the results.
    Our implementation is a feed-forward version of their model (\textit{NeuralHead-FF}) where we omit the meta-learning as well as fine-tuning phase, because we are interested in using a single model to deal with multiple identities.
\end{itemize}

\subsubsection{Metrics}

We compare the models based on the following metrics to evaluate the quality of the generated images.
Structured similarity (\textbf{SSIM})~\cite{wang-tip-2004-ssim} and peak signal-to-noise ratio (\textbf{PSNR}) evaluate the low-level similarity between the generated image and the ground-truth image.
We also report the masked-SSIM (\textbf{M-SSIM}) and masked-PSNR (\textbf{M-PSNR}) where the measurements are restricted to the facial region.

In the absence of the ground truth image where different identity drives the target face, the following metrics are more relevant.
Cosine similarity (\textbf{CSIM}) of embedding vectors generated by pre-trained face recognition model~\cite{deng-cvpr-2019-csim} is used to evaluate the quality of identity preservation.
To inspect the capability of the model to properly reenact the pose and the expression of the driver, we compute \textbf{PRMSE}, the root mean square error of the head pose angles, and \textbf{AUCON}, the ratio of identical facial action unit values, between the generated images and the driving images.
OpenFace~\cite{baltrusaitis-ieee-2018-openface} is utilized to compute pose angles and action unit values.

\section{Experimental Results}

Models were compared under self-reenactment and reenactment of different identities, including a user study.
Ablation tests were conducted as well.
All experiments were conducted under two different settings: one-shot and few-shot, where one or eight target images were used respectively.

\subsection{Self-reenactment}

\begin{table*}[t]
    \centering
    \begin{tabular}{ccccccccc}
        \toprule
Model (\# target) & CSIM$\uparrow$ & SSIM$\uparrow$ & M-SSIM$\uparrow$ & PSNR$\uparrow$ & M-PSNR$\uparrow$ & PRMSE$\downarrow$ & AUCON$\uparrow$ \\
\midrule
X2face (1) & 0.689 & 0.719 & 0.941 & 22.537 & 31.529 & 3.26 & 0.813 \\
Monkey-Net (1) & 0.697 & 0.734 & 0.934 & \bf{23.472} & 30.580 & 3.46 & 0.770 \\
NeuralHead-FF (1) & 0.229 & 0.635 & 0.923 & 20.818 & 29.599 & 3.76 & 0.791 \\
MarioNETte (1) & \bf{0.755} & \bf{0.744} & \bf{0.948} & 23.244 & \bf{32.380} & \bf{3.13} & \bf{0.825} \\
\midrule
X2face (8) & 0.762 & 0.776 & 0.956 & 24.326 & 33.328 & 3.21 & 0.826 \\
NeuralHead-FF (8) & 0.239 & 0.645 & 0.925 & 21.362 & 29.952 & 3.69 & 0.795 \\
MarioNETte (8) & \bf{0.828} & \bf{0.786} & \bf{0.958} & \bf{24.905} & \bf{33.645} & \bf{2.57} & \bf{0.850} \\

        \bottomrule
    \end{tabular}
    \caption{Evaluation result of self-reenactment setting on VoxCeleb1. Upward/downward pointing arrows correspond to metrics that are better when the values are higher/lower.}
    \label{tab:self_reenact_result}
\end{table*}

Table~\ref{tab:self_reenact_result} illustrates the evaluation results of the models under self-reenactment settings on VoxCeleb1.
\textit{MarioNETte} surpasses other models in every metric under few-shot setting and outperforms other models in every metric except for PSNR under the one-shot setting.
However, \textit{MarioNETte} shows the best performance in M-PSNR which implies that it performs better on facial region compared to baselines.
The low CSIM yielded from \textit{NeuralHead-FF} is an indirect evidence of the lack of capacity in AdaIN-based methods.

\subsection{Reenacting Different Identity}

\begin{table}[t]
    \centering
    \begin{tabular}{cccc}
        \toprule
Model (\# target) & CSIM$\uparrow$ & PRMSE$\downarrow$ & AUCON$\uparrow$ \\
\midrule
X2face (1) & 0.450 & 3.62 & 0.679 \\
Monkey-Net (1) & 0.451 & 4.81 & 0.584 \\
NeuralHead-FF (1) & 0.108 & \bf{3.30} & \bf{0.722} \\
MarioNETte (1) & \underline{0.520} & \underline{3.41} & \underline{0.710} \\
MarioNETte+LT (1) & \bf{0.568} & 3.70 & 0.684 \\
\midrule
X2face (8) & 0.484 & \bf{3.15} & 0.709 \\
NeuralHead-FF (8) & 0.120 & \underline{3.26} & \bf{0.723} \\
MarioNETte (8) & \underline{0.608} & \underline{3.26} & \underline{0.717} \\
MarioNETte+LT (8) & \bf{0.661} & 3.57 & 0.691 \\
        \bottomrule
    \end{tabular}
    \caption{Evaluation result of reenacting a different identity on CelebV. \textbf{Bold} and \underline{underlined} values correspond to the best and the second-best value of each metric, respectively.}
    \label{tab:diff_reenact_result}
\end{table}

\begin{table}[ht]
    \centering
    \begin{tabular}{cccc}
        \toprule

        Model (\# target) & \begin{tabular}[c]{@{}c@{}}vs. \\Ours\end{tabular} & \begin{tabular}[c]{@{}c@{}}vs. \\Ours+LT\end{tabular} & Realism $\uparrow$\\

        \midrule
        X2Face (1) & 0.07 & 0.09 & 0.093 \\
        Monkey-Net (1) & 0.05 & 0.09 & 0.100 \\
        NeuralHead-FF (1)& 0.17 & 0.17 & 0.087 \\
        MarioNETte (1) & - & 0.51 & 0.140 \\
        MarioNETte+LT (1) & - & - & \bf{0.187} \\
        \midrule
        \midrule
        X2Face (8) & 0.09 & 0.07 & 0.047 \\
        NeuralHead-FF (8) & 0.15 & 0.16 & 0.080 \\
        MarioNETte (8)& - & 0.52 & 0.147 \\
        MarioNETte+LT (8) & - & - & \bf{0.280} \\
        \bottomrule
    \end{tabular}
    \caption{User study results of reenacting different identity on CelebV. \textit{Ours} stands for our proposed model, \textit{MarioNETte}, and \textit{Ours+LT} stands for \textit{MarioNETte+LT}.}

    \label{tab:user_study}

\end{table}

\begin{table}[ht]
    \centering
    \begin{tabular}{cccc}
        \toprule
Model (\# target) & CSIM$\uparrow$ & PRMSE$\downarrow$ & AUCON$\uparrow$ \\
\midrule
AdaIN (1) & 0.063 & 3.47 & \underline{0.724} \\
+Attention (1) & 0.333 & \bf{3.17} & \bf{0.729} \\
+Alignment (1) & \bf{0.530} & 3.44 & 0.700 \\
\midrule
MarioNETte (1) & \underline{0.520} & \underline{3.41} & 0.710 \\
\midrule
\midrule
AdaIN (8) & 0.069 & 3.40 & \underline{0.723} \\
+Attention (8) & 0.472 & \bf{3.22} & \bf{0.727} \\
+Alignment (8) & \underline{0.605} & 3.27 & 0.709 \\
\midrule
MarioNETte (8) & \bf{0.608} &\underline{3.26} & 0.717 \\
\bottomrule
    \end{tabular}
    \caption{Comparison of ablation models for reenacting different identity on CelebV.}
    \label{tab:ablation_result}
\end{table}

Table~\ref{tab:diff_reenact_result} displays the evaluation result of reenacting a different identity on CelebV, and Figure~\ref{fig:main_result} shows generated images from proposed method and baselines.
\textit{MarioNETte} and \textit{MarioNETte+LT} preserve target identity adequately, thereby outperforming other models in CSIM.
The proposed method alleviates the identity preservation problem regardless of the driver being of the same identity or not.
While \textit{NeuralHead-FF} exhibits slightly better performance in terms of PRMSE and AUCON compared to \textit{MarioNETte}, the low CSIM of \textit{NeuralHead-FF} portrays the failure to preserve the target identity.
The landmark transformer significantly boosts identity preservation at the cost of a slight decrease in PRMSE and AUCON. The decrease may be due to the PCA bases for the expression disentanglement not being diverse enough to span the whole space of expressions. Moreover, the disentanglement of identity and expression itself is a non-trivial problem, especially in a one-shot setting.

\subsection{User Study}

Two types of user studies are conducted to assess the performance of the proposed model:

\begin{itemize}
\item \textbf{Comparative analysis.} Given three example images of the target and a driver image, we displayed two images generated by different models and asked human evaluators to select an image with higher quality.
The users were asked to assess the quality of an image in terms of (1) identity preservation, (2) reenactment of driver's pose and expression, and (3) photo-realism.
We report the winning ratio of baseline models compared to our proposed models.
We believe that user reported score better reflects the quality of different models than other indirect metrics.

\item \textbf{Realism analysis.} Similar to the user study protocol of \citet{zakharov-arxiv-2019-samsung}, three images of the same person, where two of the photos were taken from a video and the remaining generated by the model, were presented to human evaluators.
Users were instructed to choose an image that differs from the other two in terms of the identity under a three-second time limit.
We report the ratio of deception, which demonstrates the identity preservation and the photo-realism of each model.
\end{itemize}

\noindent For both studies, 150 examples were sampled from CelebV, which were evenly distributed to 100 different human evaluators.

Table~\ref{tab:user_study} illustrates that our models are preferred over existing methods achieving realism scores with a large margin.
The result demonstrates the capability of \textit{MarioNETte} in creating photo-realistic reenactments while preserving the target identity in terms of human perception.
We see a slight preference of \textit{MarioNETte} over \textit{MarioNETte+LT}, which agrees with the Table~\ref{tab:diff_reenact_result}, as \textit{MarioNETte+LT} has better identity preservation capability at the expense of slight degradation in expression transfer.
Since the identity preservation capability of \textit{MarioNETte+LT} surpasses all other models in realism score, almost twice the score of even \textit{MarioNETte} on few-shot settings, we consider the minor decline in expression transfer a good compromise.

\subsection{Ablation Test}

\begin{figure}[t]
    \centering
    \includegraphics[width=1.0\columnwidth]{{6371.ablation_qualitive}.pdf}
    \caption{(a) Driver and target images overlapped with attention map. Brightness signifies the intensity of the attention.
    (b) Failure case of \textit{+Alignment} and improved result generated by \textit{MarioNETte}.}
    \label{fig:ablation_qualitive}
 \end{figure}

We performed ablation test to investigate the effectiveness of the proposed components.
While keeping all other things the same, we compare the following configurations reenacting different identities: (1) \textbf{MarioNETte} is the proposed method where both image attention block and target feature alignment are applied. (2) \textbf{AdaIN} corresponds to the same model as MarioNETte, where the image attention block is replaced with AdaIN residual block while the target feature alignment is omitted. (3) \textbf{+Attention} is a MarioNETte where only the image attention block is applied. (4) \textbf{+Alignment} only employs the target feature alignment.

Table~\ref{tab:ablation_result} shows result of ablation test.
For identity preservation (i.e., CSIM), \textit{AdaIN} has a hard time combining style features depending solely on AdaIN residual blocks.
\textit{+Attention} alleviates the problem immensely in both one-shot and few-shot settings by attending to proper coordinates.
While \textit{+Alignment} exhibits a higher CSIM compared to \textit{+Attention}, it struggles in generating plausible images for unseen poses and expressions leading to worse PRMSE and AUCON.
Taking advantage of both attention and target feature alignment, \textit{MarioNETte} outperforms \textit{+Alignment} in every metric under consideration.

Entirely relying on target feature alignment for reenactment, \textit{+Alignment} is vulnerable to failures due to large differences in pose between target and driver that \textit{MarioNETte} can overcome.
Given a single driver image along with three target images (Figure~\ref{fig:ablation_qualitive}a), \textit{+Alignment} has defects on the forehead (denoted by arrows in Figure~\ref{fig:ablation_qualitive}b).
This is due to (1) warping low-level features from a large-pose input and (2) aggregating features from multiple targets with diverse poses.
\textit{MarioNETte}, on the other hand, gracefully handles the situation by attending to proper image among several target images as well as adequate spatial coordinates in the target image.
The attention map, highlighting the area where the image attention block is focusing on, is illustrated with white in Figure~\ref{fig:ablation_qualitive}a.
Note that \textit{MarioNETte} attends to the forehead and  adequate target images (Target 2 and 3 in Figure~\ref{fig:ablation_qualitive}a) which has similar pose with driver.

\section{Related Works}

The classical approach to face reenactment commonly involves the use of explicit 3D modeling of human faces~\cite{blanz-siggraph-1999-3dmm} where the 3DMM parameters of the driver and the target are computed from a single image, and blended eventually~\cite{thies-siggraphasia-2015-realtimeexpression,thies-cvpr-2016-face2face}.
Image warping is another popular approach where the target image is modified using the estimated flow obtained form 3D models~\cite{cao-ieee-2013-facewarehouse} or sparse landmarks~\cite{averbuch-2017-siggraph-portraitstolife}.
Face reenactment studies have embraced the recent success of neural networks exploring different image-to-image translation architectures~\cite{isola-cvpr-2017-pix2pix} such as the works of \citet{xu-arxiv-2017-facetransfer} and that of \citet{wu-eccv-2018-reenactgan}, which combined the cycle consistency loss~\cite{zhu-cvpr-2017-cyclegan}.
A hybrid of two approaches has been studied as well.
\citet{kim-siggraph-2018-deepvideoportraits} trained an image translation network which maps reenacted render of a 3D face model into a photo-realistic output.

Architectures, capable of blending the style information of the target with the spatial information of the driver, have been proposed recently.
AdaIN~\cite{huang-cvpr-2017-adain,huang-eccv-2018-munit,liu-arxiv-2019-funit} layer, attention mechanism~\cite{zhu-2019-cvpr-progressiveattention,lathuiliere-arxiv-2019-attentionfusion,park-arxiv-2019-styleattention}, deformation operation~\cite{siarohin-cvpr-2018-deformablegan,dong-2018-nips-softgatedgan}, and GAN-based method~\cite{bao2018towards} have all seen a wide adoption.
Similar idea has been applied to few-shot face reenactment settings such as the use of image-level~\cite{wiles-eccv-2018-x2face} and feature-level~\cite{siarohin-cvpr-2019-monkeynet} warping, and AdaIN layer in conjuction with a meta-learning~\cite{zakharov-arxiv-2019-samsung}.
The identity mismatch problem has been studied through methods such as CycleGAN-based landmark transformers~\cite{wu-eccv-2018-reenactgan} and landmark swappers~\cite{zhang-arxiv-2019-faceswapnet}.
While effective, these methods either require an independent model per person or a dataset with image pairs that may be hard to acquire.

\section{Conclusions}

In this paper, we have proposed a framework for few-shot face reenactment.
Our proposed image attention block and target feature alignment, together with the landmark transformer, allow us to handle the identity mismatch caused by using the landmarks of a different person.
Proposed method do not need additional fine-tuning phase for identity adaptation, which significantly increases the usefulness of the model when deployed in-the-wild.
Our experiments including human evaluation suggest the excellence of the proposed method.

One exciting avenue for future work is to improve the landmark transformer to better handle the landmark disentanglement to make the reenactment even more convincing.

\bibliographystyle{aaai}
\fontsize{9.0pt}{10.0pt} \selectfont


\section{Final instructions}

\begin{acknowledgments}
The authors thank Prof. Xavier Serra for insightful discussions on the problem.
\end{acknowledgments}


\section{Introduction}

A system in contact with a thermal reservoir relaxes when initially prepared in a state out of thermal equilibrium. This fundamental process is ubiquitous in all fields of science and is particularly interesting if, after a transient period of time, the dynamics is completely governed by only a single characteristic time scale, the inverse of the relaxation rate. The archetypical situation is that of an ensemble of particles initially confined in a potential well which is separated from an adjacent well or a continuum by a sufficiently high energy barrier. Then, the local dynamics within the well can be assumed to be much faster than barrier escape processes due to either thermal activation or quantum tunneling \cite{RevModPhys.62.251,Weiss2012}.
 The situation becomes even richer when a system in contact with a thermal bath is externally driven periodically to approach a stationary state: First, these states can be of very diverse nature, for example, following the periodicity of the external source or not (broken time-translational symmetry) \cite{PhysRevLett.117.090402,PhysRevLett.109.160401,PhysRevLett.114.251603,PhysRevB.93.174305,PhysRevLett.116.250401,PhysRevB.93.245146,PhysRevLett.118.030401,Sacha_2017};

 second, the relaxation dynamics towards these states occurs not by traversing energy barriers but rather dynamical barriers in phase space. Accurate control of the parameters of the external drive (amplitude, frequency) allows one not only to precisely access different stationary states but also to explore the relaxation dynamics towards these states in different domains of phase space \cite{Kitagawa2012,PhysRevLett.106.220402,10.1038/nature21413,nature21426}.

A particularly fascinating class of systems are periodically driven nonlinear oscillators that have received substantial attention in the last decade. Despite their putative simplicity, they reveal a wealth of dynamical features due to the subtle interplay of driving, nonlinearity, and dissipation and allow for a wide range of experimental implementations from superconducting circuits to nanomechanical systems, and  cold atomic gases \cite{Dykmanbook2012,PootvanderZant2012}.
Mostly the conventional cases of linear \cite{DykmanKrivoglaz1980,PeanoThorwart20062,SerbanWilhelm2007,GuoZhengLi2010, GuoZhengLiZhengLiYan2011,andreguopeanoschoen2012} and parametric driving \cite{DykmanMaloneySilverstein1998,PhysRevB.83.224506,DykmanMarthalerPeano2011,PhysRevB.87.184501,doi10106315116533,PhysRevApplied.8.024018} of a weakly anharmonic oscillator have been studied. These systems are conveniently described in a rotating frame by quasi-energy Hamiltonians , where quasienergy levels are associated with Floquet states. Typically, the symmetry of  the Hamiltonian in the laboratory frame with respect to time translations is reflected  in phase space by rotational symmetries. Stationary states may then appear in the rotating frame as either localized or delocalized accessible by tuning drive parameters. The transition between them is generally associated with the occurrence of slow modes, bifurcations, and the existence of multiple orbits.

At low temperatures, quantum fluctuations are the dominant source to induce switching between classical fixed points.
Accordingly, when the quantum oscillator is initially prepared close to one of the classical fixed points, relaxation to the stationary state occurs either via quantum tunneling \cite{PhysRevA.76.010102, PhysRevLett.109.090401}, in absence or for very weak dissipation, or by quantum activation, for stronger dissipation \cite{MarthalerDykmanswitchin2006, DykmanMarthalerPeano2011,PhysRevLett.110.047001, Peano_2014}. The latter phenomenon is a manifestation of the fact that in a rotating frame excitation and relaxation processes between local quasi-energy levels behave very different from the situation for Fock states in the laboratory frame. Likewise, a stationary state is in general not determined by a detailed balance condition \cite{MarthalerDykmanswitchin2006, Peano_2014} in contrast to a thermal equilibrium. Additionally, time translation symmetry breaking has been found for period-two vibrations where the state of full symmetry (the unbroken state) can merge with the broken states \cite{PhysRevE.92.022105}.  In any case, driven nonlinear quantum oscillators may serve as testbeds to explore features of phase-transition-like phenomena far from equilibrium.

In contrast to the cases of linear and parametric driving, much less attention has been paid to the study of period tripling \cite{nature21426,ZhangGosnerDykman2017,zhangDykman2019,1742-6596-681-1-012018,loerchDykman2019}, where the oscillator is periodically driven with three times its fundamental frequency. In fact, this is not just a generalization of the conventional situation but rather displays unique dynamical features, absent for linear and parametric drive, but generic for all higher order resonances \cite{PhysRevLett.111.205303,1367-2630-18-2-023006,doi10106315116533}. Classically, stationary orbits oscillating with an integer multiple of the drive period do {\em not} emerge continuously out of period-1 orbits (oscillating with the frequency of the drive). Quantum mechanically, tunneling between broken-symmetry states depends on a rotational geometric phase in phase space \cite{ZhangGosnerDykman2017}, very different from the case of the parametric oscillator \cite{PhysRevA.76.010102, PhysRevLett.109.090401}. Even weak dissipation destroys coherent tunneling and quantum activation is expected to set in.

While period tripling in the pure quantum case has been studied in \cite{ZhangGosnerDykman2017} and the quantum activation in \cite{zhangDykman2019}, here, we complement these findings by exploring in detail the relaxation dynamics towards stationary states and the dissipation-induced transition between these states being either localized or delocalized in phase space. Indeed, experimentally higher order photon resonances have already been observed in a set-up including superconducting resonators \cite{SvenssonShumeikoDelsing2017,180209259,SandboWilson2019}. They may be of relevance for the creation of higher order cat-states for quantum simulations, new types of quantum limited detectors, or to explore the fundamental physics of phase transitions far from equilibrium. Here, we lay the basis to explore the deep quantum regime and provide detailed predictions for future experiments, e.g. with superconducting circuits including Josephson junctions.

To set the stage, in Sec.~\ref{sec:clorbits} we introduce the generic model and briefly discuss the classical behavior in both the laboratory frame and the rotating frame. For the generalization to the quantum regime, we study the relaxation dynamics starting from a localized state in phase space numerically in Sec.~\ref{sec:simulation} based on a Lindblad-Master equation and in Sec.~\ref{escaperatelambda} in Liouville space for the dissipative generator of the dynamics. Section \ref{sc:semiclassical} collects some results for the switching in the semiclassical regime which allows in Sec.~\ref{sec:comparison} for a comparison of the various extracted relaxation rates and in Sec.~\ref{sec:tunneling} a comparison between ground state  quantum tunneling and quantum activation. Finally, Sec.~\ref{sec:transition} addresses the phase transition between localized and delocalized phase. Main findings are summarized in the Conclusions.

\section{Preliminaries: Classical nonlinear oscillator}
\label{sec:clorbits}
We start by briefly addressing the classical regime which besides introducing the setting and the basic notation provides a physical picture to better understand the quantum problem.

\subsection{Steady state orbits}

The system consists of a mechanical model, where a weakly anharmonic oscillator of the generic form
\begin{equation}
\ddot{q} + 2 \gamma \dot{q} + \omega_0^2 q + \alpha q^3 = F_0 q^{2} \cos(\omega_F t)
\label{eqofmot}
\end{equation}
is subject to weak damping with rate $\gamma$ and to an external driving with amplitude $F_0$ and frequency $\omega_F$. The nonlinearity is parametrized by $\alpha$ which is supposed to be weak, i.e., $(\alpha/\omega_0^2) A^2\ll 1$, where $A$ is a typical amplitude of steady state orbits (see Fig.~\ref{bifdiag}). For the driving frequency we will specifically focus on  $\omega_F \approx 3 \omega_0$ with a quadratic coupling to the oscillator degree of freedom, in contrast to the conventional situations of linear [$\omega_F \approx \omega_0$ with $F_0\cos(\omega_F t)$] and parametric [$\omega_F \approx 2 \omega_0$ with $F_0 q\cos(\omega_F t)$] driving, respectively. As we have demonstrated already in \cite{ZhangGosnerDykman2017}, this case is of particular interest as it shows features that are absent in these conventional situations but generic for all higher than second order resonances.

By assuming periodic steady state solutions of the form $q(t)=A \cos(\Omega t+\varphi)$ possible orbits can easily be obtained from (\ref{eqofmot}), for details see App.~\ref{appendix}. Two types of orbits are found, namely,  those oscillating with the  frequency $\Omega=\omega_F$ of the external drive, and those oscillating with $\Omega=\omega_F/3 \approx \omega_0$, the fundamental frequency of the bare oscillator; the former are termed period-1 orbits, the latter period-3 orbits. A stability analysis reveals that there is a stable and an unstable branch of period-3 orbits, see Fig.\ref{bifdiag}(a). Period-3 orbits only exist beyond a threshold for the driving strength $F_0$ and do not grow continuously out of period-1 solutions, in contrast to the situation for parametric amplification, where period-2 solutions smoothly emerge out of period-1 orbits.
While the amplitudes for stable period-3 orbits grow with increasing driving, those for unstable period-3 orbits shrink and asymptotically approach the stable branch of the period-1 orbits. Amplitudes for the latter always remain small due to the far-off-resonance driving, see Fig.\ref{bifdiag}(b). %
According to the time translational symmetry of (\ref{eqofmot}), period-3 solutions appear (up to an overall off-set) with three different phases $\varphi=0, 2\pi/3, 4\pi/3$, see Fig.\ref{bifdiag}(c). Note that if initially one prepares the oscillator in a period-1 state for weak driving and adiabatically increases the driving, in absence of external perturbations, the oscillator never approaches period-3 orbits.
	 \begin{figure}[htp]
\centering
\includegraphics[width=1\linewidth]{Fig1.pdf}
 \caption{(a): Classical bifurcation diagram in rotating frame depicting steady state amplitudes [in units of $C=\sqrt{8 \omega_0 \delta \omega/3 \alpha}$] versus driving strength [scaled as in (\ref{scaleddriving})] for $\gamma/\omega_0=0.01$, $\omega_F/\omega_0=3.02$; the threshold beyond which multiple orbits coexist is $f_{c, \rm cl}\approx 0.49$. (b),(c): Steady state orbits in the laboratory frame solving \eqref{eqofmot} for $f=0.76$. The period-1 orbit (b)  is oscillating with the driving frequency $\omega_F$ (note the scale of the amplitude). Stable period-3 orbits (c) oscillating with $\omega_F/3 \approx \omega_0$  occur with relative phase-shifts $2\pi/3$, $4\pi/3$;
 other parameters are $\alpha L^2/\omega_0^2=0.003$, $F_0 L/\omega_0^2=0.05$ with arbitrary length scale $L$.}
 \label{bifdiag}
\end{figure}

\subsection{Quasi-energy in the rotating frame}
The non-dissipative part of (\ref{eqofmot}) results from the time-dependent Hamiltonian
\begin{equation}
H=\frac{p^2}{2}+\frac{\omega_0^2 q^2}{2}+ \frac{\alpha q^4}{4} - \frac{F_0 q^3}{3} \cos(\omega_F t)\, .
\end{equation}
If the driving is not too strong, so that the driving term, the anharmonic term, and $\omega_0 \delta\omega q^2$ with  de-tuning $|\delta\omega|=|\omega_F/3-\omega_0|\ll \omega_0$ are small compared to the bare harmonic part, it is convenient to map this Hamiltonian in the laboratory frame via a canonical transformation of the form
\begin{equation}
q  = C Q \cos\left(\frac{\omega_F t}{3}\right) + C P \sin\left(\frac{\omega_F t}{3}\right)
\label{unitarytrafo}
\end{equation}
to a rotating frame Hamiltonian $H_{RWA}$. Here, we introduced dimensionless quadratures $Q$, $P$ and  a scaling factor $C=\sqrt{\frac{8 \omega_0 \delta \omega}{3 \alpha}}$ (we assume $\delta\omega>0$ in the sequel). This way, one finds the dimensionless quasi-energy
\begin{equation}
\label{quasienergy}
g(Q,P)\equiv \frac{8 H_{RWA}}{3 \alpha C^4}=\frac{1}{4} (Q^2+P^2-1)^2-  f (Q^3-3 P^2 Q),
\end{equation}
with a dimensionless driving strength
\begin{equation}
\label{scaleddriving}
f=\frac{F_0}{3\sqrt{24 \omega_0 \alpha \delta \omega}}\, .
\end{equation}
	 \begin{figure}[htp]
\centering
\includegraphics[width=1\linewidth]{Fig2.pdf}
\caption{Classical quasi-energy $g(Q,P)$ according to (\ref{quasienergy}) for different driving strengths $f=0.1$ (a) and $f=0.5$ (b). The quasi-energy has a $C_3$ symmetry in $QP-$ space, where the local maximum at $Q_0=0, P_0=0$ corresponds to the period-1 orbit in Fig.~\ref{bifdiag}(b) and the three wells located at $(Q_m,P_m)$, $m=1,2,3$ are related to the period-3 orbits in Fig.~\ref{bifdiag}(c). }
\label{classpotg}
\end{figure}
Note that here the definition of $f$ differs by a numerical factor from \cite{ZhangGosnerDykman2017, zhangDykman2019}.
This quasi-energy has a three-fold symmetry in $QP-$ space ($C_3$ symmetry), see Fig.~\ref{classpotg}, that reflects the time-translational symmetry of the three steady state period-3 orbits discussed in the previous section. Accordingly, a local maximum at $Q_0=0, P_0=0$ associated with the period-1 orbit with low amplitude is complemented by three wells located at $(Q_m,P_m)$, $m=1,2,3$ corresponding to the period-3 orbits with
\begin{equation}
\begin{split}
&P_1=0 \text{,  } Q_1=Q_{\rm min}\\
&P_{2,3}= \pm \sin(2\pi/3)\,  Q_{\rm min} \text{,  } Q_{2,3}=\cos(2\pi/3)\,  Q_{\rm min}\, ,
\end{split}
\end{equation}
where $Q_{\rm min}=\frac{3}{2}f + \sqrt{\frac{9}{4}f^2+1}$  (see \cite{ZhangGosnerDykman2017}). Rotations in the $QP$-plane  by $2\pi/3$ correspond  to time translations in the laboratory frame $t\to t+t_F, t_F=2\pi/\omega_F$, so that rotations by $2\pi$ map onto $t\to t+3t_F$, see \cite{ZhangGosnerDykman2017} for further details. The three wells are separated by saddle points with quasi-energy $g_{\rm saddle}<0$. The equations of motion in the rotating frame in presence of dissipation easily follow from Hamilton's equations augmented by friction, i.e.,
\begin{equation}
\begin{split}
&\dot{Q}= \partial_P g - \kappa\,  Q\\
&\dot{P} = - \partial_Q g - \kappa\,  P
\label{eq:classical}
\end{split}
\end{equation}
with the dimensionless friction constant $\kappa=\frac{\gamma}{\delta \omega}$. In App.~\ref{appendix} we discuss a protocol how to realize experimentally the situation in (\ref{eqofmot}) via a conventional Duffing oscillator with linear drive. We note in passing that time translational symmetry breaking in the classical realm has  recently been studied experimentally in a two-oscillator set-up in \cite{heugel2019}.

\section{Dynamics towards steady state in the quantum regime}
\label{sec:simulation}

We will now turn to the main subject of this paper, namely, the quantum induced switching out of domains in $QP$-space around classical fixed points. The quantum equivalent to the canonical transformation \eqref{unitarytrafo} is provided by the unitary operator ${U}=\exp({-i \frac{\omega_F t}{3} {a}^{\dagger} {a}})$ with the standard ladder operators of the bare harmonic system obeying $[a, a^\dagger]=1$.   The quantized rotating frame Hamiltonian is then obtained as
\begin{equation}
\begin{split}
\hat{g}=& - \Lambda n + \frac{1}{4}\left(1-\Lambda\right)^2 + \Lambda^2 (n+n^2) - 4 f \left(\frac{\Lambda}{2}\right)^{\frac{3}{2}} (a^3+ a^{\dagger 3} )\,
\end{split}
\label{quasi-operator}
\end{equation}
with the effective Planck constant
\begin{equation}
\label{lambda}
\Lambda=\frac{ \hbar}{M \omega_0 C^2}\equiv \frac{3\hbar\alpha}{M 8\omega_0^2\delta\omega} \,,
\end{equation}
including the mass $M$ of the oscillator. The effective Planck constant serves as an externally tunable parameter, for example via the de-tuning $\delta\omega$.
The operator $\hat{g}$ has an equivalent representation in terms of phase space operators $\hat{Q}, \hat{P}$ with canonical commutation relations  $[\hat{Q},\hat{P}]=i \Lambda$ according to ${a} = \frac{1}{\sqrt{2 \Lambda}} ({\hat{Q}} + {i \hat{P}})$ and ${a}^{\dagger} = \frac{1}{\sqrt{2 \Lambda}} (\hat{Q} - {i \hat{P}})$. This way, $\hat{g}(\hat{Q}, \hat{P})$ results from the classical quasi-energy $g(Q,P)$ in (\ref{quasienergy}) by replacing classical variables by operators and $P^2 Q \to \hat{Q}\hat{P}\hat{Q}$.

For the dissipative dynamics of the quantum oscillator we invoke a weak coupling approximation which leads to a standard
Lindblad master equation
\begin{equation}
\begin{split}
\dot{\rho}(t) =& \frac{1}{i \Lambda} [g,\rho] +  \kappa (\bar{n}+1) \left(2a \rho a^\dagger -  a^\dagger a \rho -  \rho  a^\dagger a\right)\\
&+ \kappa \bar{n} \left(-2a^\dagger \rho a + a a^\dagger \rho +  \rho a a^\dagger\right)\, ,
\end{split}
\label{Lindblad}
\end{equation}
where time is scaled with $\delta \omega$ and $\bar{n}=1/[\exp(\hbar \omega_0\beta)-1]$ at inverse temperature $\beta=1/k_{\rm B}T$ denotes the Bose distribution.

With the above ingredients and starting from a specific initial state, one can now explore the quantum dynamics analogous to the discussion of the classical steady state orbits above. Here, we are particularly interested in the quantum fluctuations of one of the wells around $(Q_m, P_m), m>0$ in the limit of vanishing temperature $\bar{n}=0$. For that purpose, we consider the dynamics of a locally relaxed state namely the stationary solution of the local Fokker-Planck equation in the harmonic well around $Q_1$, $P_1$ \cite{Dykmanbook2012,PhysRevE.75.011101}. This way we obtain in the Wigner representation
\begin{equation}
\rho_W(Q, P,t=0) \propto \exp\left[- \frac{2}{\Lambda} \left(\frac{g_{QQ} \delta Q^2 + g_{PP} \delta P^2} {g_{PP}+g_{QQ}} \right) \right]\, ,
\end{equation}
 where $\delta Q=Q-Q_1$, $\delta P=P-P_1$, and $g_{XY}$ is the second derivative with respect to $X, Y$ taken at $Q_1, P_1$.

Monitoring the dynamics of this initial state according to (\ref{Lindblad}) we show snapshots of corresponding Wigner distributions at different time steps in Fig.~\ref{Wignerlindz}. Apparently, in the long time limit the distribution approaches a steady state that is delocalized predominantly among the three wells with minimal contributions around the origin. The set of parameters is chosen such that classically the system is beyond the threshold (see Fig.~\ref{bifdiag}), where period-1 and period-3 orbits coexist. This in turn implies that quantum fluctuations induce the decay of classically stable fixed points towards delocalized distributions with $C_3$-symmetry. Note the slight distortions of the arms of the distribution due to friction.
 %
	 \begin{figure}[htp]
\centering
\includegraphics[width=1\linewidth]{Fig3.pdf}
 \caption{Wigner distribution $\rho_W(Q,P,t)$ starting initially from a locally relaxed state in the well $Q_1, P_1$ of the quasi-energy surface at different time steps and for $\kappa=1.5$, $\Lambda=0.1688$, $f=0.9$: $t=33$ (a), $t=1636$ (b), $t=4975$ (c), and asymptotic state (d). The population dynamics of a closed area indicated by a white box in (a) is depicted in Fig.~\ref{decayswitching}, such that the population $p_1(0)\approx 1$ at time $t=0$.}
 \label{Wignerlindz}
\end{figure}

After a transient period of time, this relaxation is very accurately characterized by a single time scale $\Gamma_L$ when following the population dynamics
\begin{equation}
p_1(t) = \int_{B_1} dQ dP \, \rho_W(Q,P,t)\,
\end{equation}
according to $\dot{p}_1(t)=- \Gamma_L [p_1(t) - p_1(\infty)]$, see Fig.~\ref{decayswitching}. Here, $B_1$ indicates a closed area in $QP$-space around $Q_1, P_1$ (white box in Fig.~\ref{Wignerlindz}(a)) such that initially $p_1(0)\approx 1$. Asymptotically, one has $p_1(\infty)<1/3$ due to the portion of the distribution that is still located around $Q=P=0$.
\begin{figure}[htp]
\centering
\includegraphics[width=1\linewidth]{Fig4.pdf}
 \caption{Decay of the population in the well around $Q_1, P_1$ (in the white box in Fig.~\ref{Wignerlindz}(a)) versus dimensionless times for $\kappa=1.5$, $\Lambda=0.1688$, $f=0.9$. On long time scales it is governed by a single time constant $1/\Gamma_L$. The black line is a fit function, where $p_1(\infty)=0.32$. The inset shows the transient behavior on very short time scales (blue); also shown is the decay of a state initially located around the origin at $Q=P=0$ (green). The initial distribution at $Q=P=0$ decays much faster than the initial population around the well $Q_1, P_1$.}
 \label{decayswitching}
\end{figure}
Apparently, the observed dynamics is purely relaxational in contrast to coherent oscillations between the wells that are expected to occur in complete absence of dissipation and on different time scales (cf.\ also below in Sec.~\ref{sec:tunneling}). As we discussed in \cite{ZhangGosnerDykman2017}, due to the energy level structure of (\ref{quasi-operator}) sufficiently beyond the threshold, only exponentially weak friction is needed to induce decoherence (associated with transitions between tunnel-splitted energy levels) and mixing of corresponding eigenstates. For friction that  substantially exceeds the tunnel splitting (as considered here) the coherent dynamics turns into uni-directional decay. Thus, on the time scale associated with $\Gamma_L$ there occur transitions between the localized wells  (interwell transitions) that, as we will discuss in detail in Secs.~\ref{sc:semiclassical} and \ref{sec:tunneling}, are attributed to the phenomenon of quantum activation. A condition for the long time dynamics to be governed by this single time scale is a time scale separation between interwell and intrawell processes,  an issue that we will explore more carefully in the next section.

Furthermore, we show in Fig.~\ref{decayswitching} that starting initially from a ground state distribution localized around $Q=P=0$ (classical period-1 orbit) provides a much faster decay towards the delocalized asymptotic steady state.
The time scale associated with this decay turns out to play an important role near those domains in parameter space, where the nature of the steady state changes from localized to delocalized, i.e. the quantum analogue of the classical bifurcation as will be discussed in Sec.~\ref{sec:transition}.

\section{Time scale separation between intrawell and interwell processes}
\label{escaperatelambda}

In the previous section, the time scale which governs the long time dynamics was obtained from the time evolution of the Lindblad master equation $\dot{\rho}(t)=\mathcal{L} \rho(t)$ with the super-operator $\mathcal{L}$ as given in (\ref{Lindblad}). Here, we extend this analysis to access also shorter time scales which capture not only inter-well but also intra-well dynamics. The most convenient way to do so, is to switch to a Liouville representation, where the density is considered as a vector in Liouville space on which  $\mathcal{L}$ operates.  The goal is then to find eigenfunctions and eigenvalues of  $\mathcal{L}$ according to
\begin{equation}
\mathcal{L} \ket{\ket{v_k}} = \lambda_k \ket{\ket{v_k}} ,
\end{equation}
where the density matrix at time $t=0$ is represented as
\begin{equation}
\ket{\ket{\rho(0)}}=\sum_{k\geq 0} a_k \ket{\ket{v_k}}\, .
\end{equation}
Of course, due to the existence of a steady state $\mathcal{L} \rho_{ss}=0 $ there is always a vanishing eigenvalue $\lambda_0=0$; all other eigenvalues have negative real parts ${\rm Re}\{\lambda_k\}<0$ and are either real or appear in complex conjugate pairs. The Lindblad time evolution in Liouville space
\begin{equation}
\ket{\ket{\rho(t)}}= e^{\mathcal{L} t} \ket{\ket{\rho(0)}}=a_0 \ket{\ket{v_0}}+\sum_{k\geq 1} a_k e^{\lambda_k t} \ket{\ket{v_k}}  ,
\end{equation}
is completely determined by the initial state projections $a_k$ and the eigenvalues $\lambda_k$,  and asymptotically approaches $\rho_{ss}= \lim\limits_{t \rightarrow \infty}{\rho(t)}\sim \ket{\ket{v_0}}$ .

In the sequel, we focus at zero temperature $\bar{n}=0$ on the lowest eigenvalues. The Liouville operator $\mathcal{L}$ is converted accordingly to a $N^2\times N^2$ matrix which is diagonalized using the Lanczos/Arnoldi algorithm.
This provides eigenvectors  $\ket{\ket{\nu_l}}, l\geq 0$ to each eigenvalue $\lambda_l$ which are linear independent, but not orthogonal.

Upon representing eigenvectors $\ket{\ket{\nu_l}}$ again in $N\times N$ matrices, i.e.\ $\ket{\ket{\nu_l}}\to \rho^{(l)}$, one can classify the eigen-densities into three classes according to the following structure in the Fock state basis
\begin{eqnarray}
\ket{\ket{\nu_{3k}}}\to \rho^{(3k)} &=& \sum_{n \ge 0}^N p_{nn}^{(k)} \ket{n} \bra{n} \nonumber\\
&&+ \sum_{n,m\ge 0}^N \left(\alpha_{nm}^{(k)} \ket{n} \bra{3m+3+n}+h.c.\right) \nonumber\\
\ket{\ket{\nu_{3k+1}}}\to \rho^{(3k+1)} &=&  \sum_{n,m\ge 0}^N \Big(\gamma_{nm}^{(k)} \ket{n} \bra{3m+1+n}\nonumber\\&& + \delta_{nm}^{(k)}  \ket{3m+2+n} \bra{n}\Big)\nonumber\\
\ket{\ket{\nu_{3k+2}}}\to  \rho^{(3k+2)} &=&{\rho^{(3k+1)}}^\dagger \ \ , \ \ k\geq 0
\end{eqnarray}
with real-valued coefficients $p_n^{(k)}$ and complex-valued $\alpha_{nm}^{(k)}, \gamma_{nm}^{(k)}, \delta_{nm}^{(k)}$. Eigen-density matrices have matrix elements $\rho_{nm}$ where $n-m$ is identical $(\text{mod } 3)$ within the same class. This structure is explained by the fact that both the Hamiltonian \eqref{quasi-operator} part as well as the dissipative part of the Lindblad master equation \eqref{Lindblad} only couple such matrix elements. Diagonal elements of $\rho$ (populations) only appear in the class $\rho^{(3k)}$ associated with real eigenvalues $\lambda_{3k}$, while $\lambda_{3k+1}=\lambda_{3k+2}^*$. As we will see in Sec.~\ref{sec:transition}, for sufficiently weak driving the dominating contribution to the steady state density $\rho^{(0)}$ is the ground state density of the oscillator around $Q_0, P_0$, while for stronger driving beyond a threshold strongly delocalized  components prevail.
	 \begin{figure}[htp]
\centering
\includegraphics[width=0.5\textwidth]{Fig5.pdf}
 \caption{Real (a) and imaginary (b) parts of the lowest lying eigenvalues of the Lindbladian in (\ref{Lindblad}) versus the driving for $\Lambda=0.1688$ and $\kappa=1.5$. Shown
 are eigenvalues corresponding to the class $\rho^{(3k)}$ (blue), namely,
  $\lambda_0$ (\textcolor{hblue}{$\diamond$}) and $\lambda_3$ (\textcolor{hblue}{$\bullet$}), and eigenvalues corresponding to the classes $\rho^{(3k+1)}, \rho^{(3k+2)}$ (black), namely, ${\rm Re}\{\lambda_1\}={\rm Re}\{\lambda_2\}$ ($\nabla$). For very weak driving
  ${\rm Re}\{\lambda_{1, 2}\}$ describe relaxation of a harmonic oscillator near $Q=P=0$, while beyond a threshold $f_c$ (here at the maximum of $\lambda_3$) the eigenvalues become exponentially small (see inset of (a) with a log-scale) and capture interwell transition processes. The imaginary part (b) of $\lambda_{1,2}$ also decreases towards $\lambda_0$.}
\label{realimageigliou}
\end{figure}

Here, we proceed with the analysis of the subset of eigenvalues with least negative real parts.
Figure \ref{realimageigliou}(a) depicts together with $\lambda_0=0$ the three eigenvalues ${\rm Re}\{\lambda_1\}={\rm Re}\{\lambda_2\}, \lambda_3$ versus the driving strength. For very weak driving these eigenvalues reflect  the expected results for a harmonic oscillator located at $Q_0=P_0=0$, i.e.\ $\lambda_{1,2}=- \kappa$, $\lambda_3 = - 2 \kappa$, according to the low energy properties of $\hat{g}$ in (\ref{quasi-operator}) for $f=0$.
 With increasing driving, however, near a critical driving $f_c$, $ {\rm Re}\{\lambda_{1, 2}\}$ turn into an exponential decrease towards $\lambda_0=0$ for $f>f_c$ [see Fig.~\ref{realimageigliou}(a), inset]; a similar behavior is found for their imaginary parts [Fig.~\ref{realimageigliou}(b)]. The appearance of a new (exponentially small) time scale is associated with the appearance of a new dynamical process, namely, the switching between wells separated by quasi-energy barriers (dynamical barriers in phase space). For the present case, it is also accompanied by the tendency of closing of the gap between $\lambda_0=0$ and the next lowest lying eigenvalue $\lambda_3$, a signature typically associated with the emergence of slow modes well-known for phase transition phenomena \cite{RevModPhys.69.315}.

 It becomes an exact closing only for $\Lambda\to 0$ which may thus be interpreted as the thermodynamic limit for this system, see also Sec.~\ref{sec:transition} and Fig.~\ref{changeinlambda}. For any finite $\Lambda$ the transition is smeared out in a range around $f_c$.

Thus, we identify a critical driving $f_c\approx 0.63$  as that driving strength, where $|\lambda_3|$ approaches a minimum; as we will see in Sec.~\ref{sec:transition}, around this point the steady state changes from being localized around the origin of the $QP$-plane to being predominantly localized around $Q_m, P_m, m=1, 2, 3$. Of course, there is no sharp transition for finite $\Lambda$.  All other eigenvalues $\lambda_{k>3}$ remain sufficiently separated from $\lambda_3$ near $f_c$. Beyond this threshold, $|\lambda_3|$ increases again to become well-separated from ${\rm Re}\{\lambda_{1, 2}\}$ and dives down into the pattern of the other eigenvalues $\lambda_{k>3}$. For increasing $f>f_c$ this pattern approximately captures that of a harmonic oscillator localized in one of the wells $Q_m P_m, m=1, 2, 3$ (not shown). We thus conclude that for $f>f_c$, one has a clear time scale separation between intrawell and interwell processes so that the long time relaxation behavior is governed by a single time scale, i.e.\ $|{\rm Re}\{\lambda_{1, 2}\}|=\Gamma_L$. Near the critical driving $f_c$, another dynamical role is assigned to the eigenvalues $\lambda_{1,2}$: While for very weak driving they are related to the relaxation dynamics of the local harmonic oscillator at the origin of the $QP$-plane, for $f$ sufficiently beyond $f_c$ their real parts correspond to interwell transition processes. How this switching between classical period-3 states happens is in more detail elucidated in Secs.~\ref{sc:semiclassical} and \ref{sec:tunneling}. A general theory connecting metastability and separations in the eigenvalue spectrum of Markovian open quantum systems has been put forward in \cite{PhysRevE.94.052132,PhysRevLett.116.240404}.

We found a slightly modified pattern of the lowest lying eigenvalues for parameters $\kappa<\Lambda$, which, however, does not change the conclusions about the transition. Specifically for those parameters the eigenvalues exhibit a more complex dependence on the driving below $f_c$ and the critical behavior is observed for driving values slightly below the minimum of $|\lambda_3|$. For driving strengths near and above $f_c$ one finds a similar behavior as for $\kappa>\Lambda$ though.

\section{Semiclassical switching rate}
\label{sc:semiclassical}
In order to obtain a deeper insight into the interwell relaxation process, we employ a semiclassical treatment applicable when formally Planck's constant $\Lambda \ll 1$. Practically, this implies a large number of states in the potential wells of $g(Q, P)$ such that  $|g(Q_m, P_m)|/\Lambda \sim f^4/\Lambda\gg 1$ and also that $\kappa$ by far exceeds the tunnel splitting and switching rate to ensure a time scale separation (see below).
Technically, we follow the methodology explained in \cite{MarthalerDykmanswitchin2006}: One represents the master equation  \eqref{Lindblad} around one of the well basins in corresponding intrawell eigenstates (Wannier basis) $|\mu\rangle$ of $\hat{g}$ that are linear combinations of the global eigenstates, see \cite{ZhangGosnerDykman2017}. Since we are interested in the population dynamics, only the diagonal part $\rho_\mu\equiv \langle \mu|\rho|\mu\rangle$ is considered. We assume a quasi-continuum of energy levels in the wells so that quantum mechanical matrix elements can be obtained from classical orbits. Within this semiclassical approximation, the steady state distribution is obtained based on an exponential ansatz. For a system in thermal equilibrium this would lead to a Gibbs distribution, while here we find an exponent that depends nonlinearily on quasi-energies. Following these lines, we start from
\begin{equation}
\label{eq:lindblad2}
\dot{\rho}_{\mu}  = -2\kappa \sum_\nu \left(W_{\mu \nu}\rho_\mu - W_{\nu \mu}\rho_\nu\right)\,
\end{equation}
with transition rates
\begin{equation}
W_{\mu+\nu, \mu}=(1+\bar{n}) |\langle \mu|a|\mu+\nu\rangle|^2 +\bar{n} |\langle\mu+\nu| a|\mu\rangle|^2\, .
\end{equation}
For vanishing temperature $\bar{n}=0$ only transitions from $\mu+\nu \rightarrow \mu$ take place. Semiclassically, the relevant matrix elements follow from
\begin{equation}
\begin{split}
a_\nu(g_\mu) & \equiv \sqrt{2 \Lambda} \bra{\mu} a \ket{\mu+\nu}\\
&=\frac{1}{T(g)} \int_{P.O.} dt\,   {\rm e}^{- i \omega_\mu  \nu t}  [Q(g_\mu,t)+i P(g_\mu,t)]
\end{split}
\end{equation}
 where we set $(g_{\mu+\nu}-g_\mu)/\Lambda\approx \nu \omega(g_\mu)$ with $g_\mu=\langle \mu|\hat{g}|\mu\rangle$. The integration is taken along periodic orbits (P.O.) in the well area with position $Q(g,t)$ and momentum  $P(g, t)$ at fixed energy $g$ with period $T(g)=2\pi/\omega(g)$;  these are period-3 orbits in the language of Sec.~\ref{sec:clorbits}. Note that since the right hand side of (\ref{eq:lindblad2}) is of order $\kappa$, the orbits are obtained from (\ref{eq:classical}) in absence of friction.

The above matrix elements carry interesting information about intrawell processes such as relaxation and quasi-energy diffusion in the rotating frame as we will discuss now.
Figure \ref{asymmetricmatrixel}(a) depicts $a_\nu(g)$ versus $\nu$  at different values of $g$, where in accordance with the semiclassical limit we assume a quasi-continuum of energies $\{g_\mu\}\to g_{\rm min}\leq g\leq g_{\rm saddle}$. The general tendency is that the matrix elements decrease exponentially with growing $\nu$ and increase as the energy grows towards the saddle point.  Further, $|a_\nu|<|a_{-\nu}| (\nu>0)$ with a slight asymmetry in the decay as $\nu$ increases (for details and analytical results cf.~\cite{zhangDykman2019}). This unusual behavior at zero temperature is a manifestation of the fact that relaxation processes in the laboratory frame appear in the rotating frame as both excitation and relaxation processes which gives rise to the phenomenon of quantum activation \cite{MarthalerDykmanswitchin2006}.
Further, the asymmetry in the decay with growing $\nu$  has direct consequences on the detailed balance condition
	 \begin{figure}[htp]
\centering
\includegraphics[width=1\linewidth]{Fig6.pdf}
 \caption{(a): Transition matrix elements $|a_\nu(g)|$ at two different quasi-energies $g=-1.2$ (diamonds), $g=-0.2$ (dots) for $f=0.66$ [$g_\text{min}=-3.55$, $g_\text{saddle}=0.22$]. Shown are $|a_\nu|$ (blue) and $|a_{-\nu}|$ (red), where $|a_\nu|<|a_{-\nu}| (\nu>0)$. The matrix elements decrease exponentially with growing $\nu$ and increase as energy grows towards $g_\text{saddle}$. (b): Detailed balance ratio $r_{\mu\nu}$ according to (\ref{detailedratio}) for $r_{42}$ (blue) and $r_{31}$ (red) at $f=0.5$. The ratio $r_{\mu\nu}$ differs from $1$ in both cases so that the detailed balance condition is not fulfilled.}
\label{asymmetricmatrixel}
\end{figure}
\begin{equation}
\frac{W_{\mu\nu}}{W_{\nu\mu}}=\frac{W_{\mu\mu'} W_{\mu'\nu}}{W_{\nu\mu'} W_{\mu'\mu}}\, .
\end{equation}
As seen in Fig. \ref{asymmetricmatrixel}(b), it is not fulfilled , even for $\bar{n}=0$, since the ratio
\begin{equation}
\label{detailedratio}
r_{\mu\nu}=\frac{|a_\nu|^2 |a_{-\nu+\mu} a_{-\mu}|^2}{|a_{-\nu}|^2 | a_{\nu-\mu} a_{\mu}|^2}
\end{equation}
clearly differs from 1. This in turn is a signature of the fact that in the rotating frame dissipation induced relaxation and excitation processes behave quite differently to the situation in the laboratory frame. Notably, while for the Duffing oscillator with parametric drive, detailed balance is at least guaranteed at zero temperature, we here find detailed balance to be violated for {\em all} temperatures.

Now, in steady state according to an eikonal ansatz the density is of the form $\rho_\mu\sim \exp(-R_\mu/\Lambda)$ with $R_\mu=R(g_\mu)$ so that $\rho_{\mu+\nu}\approx \rho_\mu \exp[-\nu \omega_\mu R'(g_\mu)]$,  $R'=dR/dg$. The balance equation (\ref{eq:lindblad2}) can then be cast in the form
\begin{equation}
\label{stat-master}
\sum_{\nu=-\mu}^{\mu_\text{saddle}}= W_{\mu+\nu,\mu} \left(1-\xi^\nu_{(\mu)} \right)=0,
\end{equation}
where the parameter $\xi_{(\mu)} = e^{-R'_\mu\omega_\mu}$ is independent of $\Lambda$ in leading order.
The above equation can easily be solved numerically for the parameter $\xi$ from which $R'(g)$ is obtained. This way,
  by assuming a simple relation between steady state density and interwell relaxation rate (similar to the relation between the Boltzmann distribution and the thermal escape rate for barrier escape problems in absence of driving), one finds
\begin{equation}
\Gamma_{\rm scl} = D_0 \,  {\rm e}^{-\frac{R(g_{\rm saddle})}{\Lambda}}
\label{escsem}
\end{equation}
with an unknown prefactor $D_0$ and the action $R(g)=\int_{g_{\rm min}}^{g} dx R'(x) $, known for quantum activated processes \cite{MarthalerDykmanswitchin2006}. The calculation of the prefactor is much more challenging. According the known procedure for static barriers, one would seek in the low viscosity limit for a flux solution which describes deviations from the stationary distribution for energies close to the barrier top. The matching between this solution and the stationary one then provides the prefactor. It is not clear yet how such a procedure can be extended to phase space barriers and the deep quantum regime. The classical Kramers result \cite{RevModPhys.62.251} suggests $D_0\approx \kappa$ for sufficiently weak friction, while for somewhat larger friction
$D_0\approx 1/T(g_{\rm min})$ is the expected transition state theory (TST) result. This issue will be discussed in the next section.

One important remark is in order here: The basic assumption underlying the above treatment, particularly (\ref{stat-master}), is that state populations remain stationary due to transitions between neighboring states. As a recent study revealed \cite{zhangDykman2019}, however, this is in a strict sense no longer true for the present case, but only applies approximately for values of $f$ sufficiently away from $f_c$. In fact, as shown by Zhang and Dykman, stationarity is maintained by non-local transitions so that, at least in principle, one must start from the full balance equation (\ref{eq:lindblad2}) to determine the stationary distribution. Nevertheless, in case of a sufficient time scale separation between intrawell processes (fast) and interwell decay (much slower), the dominating contribution to the interwell rate (i.e.\ on a logarithmic scale) is still expected to be fairly accurately captured by (\ref{escsem}).

Figure \ref{RandRprime}(a) illustrates the dependence of  the action $R(g)$ on $g$ for various temperatures $\bar{n}$.
	 \begin{figure}[htp]
\centering
\includegraphics[width=1\linewidth]{Fig7.pdf}
\caption{(a): Action $R(g)$ at $f=0.66$ [$g_{\rm min}=-3.55, g_{\rm saddle}=0.22$] and for different temperatures: $\bar{n}=0$ (black), $\bar{n}=0.1$ (blue), $\bar{n}=0.2$ (red). The action increases linearly with $g$. (b): Saddle point action $R(g_{\rm saddle})$ determining the rate for quantum activation versus the driving strength at $\bar{n}=0$ (black), $\bar{n}=0.1$ (blue), $\bar{n}=0.2$ (red). $R(g_{\rm saddle})$ increases almost linearly with increasing driving beyond $f_c$. }
\label{RandRprime}
\end{figure}
The action grows basically linearly with $g$ with slight deviations for elevated temperatures. The action for interwell transitions $R(g_{\rm saddle})$ is shown in Fig. \ref{RandRprime}(b) versus the driving strength and reveals a monotonously growing, nonlinear behavior which turns into an almost linear one beyond $f_c$.

\section{Comparison of switching rates}
\label{sec:comparison}

	 \begin{figure}[htp]
\centering
\includegraphics[width=1\linewidth]{Fig8.pdf}
\caption{Semiclassical rate for interwell transitions via quantum activation $\Gamma_{\rm scl}$ according to (\ref{escsem}) with $D_0= 1/T(g_{\rm min})$ (purple and blue) and $D_0= \kappa$ (orange and turquoise) for $\Lambda=0.1688$ and $\Lambda=0.1125$ together with the numerically obtained switching rates $\Gamma_L=|{\rm Re}\{\lambda_{1, 2}\}|$  at $\Lambda=0.1688$ (red crosses) and  at $\Lambda=0.1125$ (green crosses) versus the friction constant $\kappa$  at $f=0.8$.}
\label{kappadepf}
\end{figure}

We now collect the results of the previous sections with respect to the relaxation time scale in the regime $f>f_c$. First, as to be expected, one finds that numerically to very high accuracy
\begin{equation}
\Gamma_L\approx | {\rm Re}\{\lambda_{1, 2}\}|,
\end{equation}
for driving strengths sufficiently beyond the threshold, cf.~Fig.~\ref{realimageigliou}(a). The rate exponent depends roughly linearly on $f$ for $f>f_c$ [see inset of Fig.~\ref{realimageigliou}(a)], a behavior that is approximately also found for the exponent in $\Gamma_{\rm scl}$, see  Fig.~\ref{RandRprime}(b). However, while $\Gamma_L$ exhibits a dependence on the friction strength $\kappa$, this is by construction not the case for $\Gamma_{\rm scl}$ (see Fig.~\ref{kappadepf}) if $D_0 \propto 1/T(g_{\rm min})$ is chosen. Instead for $D_0 \propto \kappa$ we find a decreasing rate with decreasing $\kappa$ with substantial deviations, however, for growing $\kappa > 1$ (see Fig.~\ref{kappadepf}). Note that for the parameters chosen in Fig.~\ref{kappadepf} one has $V_b/\Lambda \omega_1 \approx 12$ ($\Lambda=0.1688$) or $\approx 18$ ($\Lambda=0.1125$) levels in the well ($V_b$ is the well depth and $\omega_1$ the frequency near the well bottom). The diagonalization of the full Lindbladian according to Sec.~\ref{escaperatelambda} is limited to values of the friction with $\kappa>0.9$. Nonetheless, the dominating exponential in $\Gamma_L$ is captured by the simple semiclassical expression $\Gamma_{\rm scl}$, while a more refined description is required for the prefactor. This must also include the impact of the period-1 oscillator located around $Q=P=0$ and the dynamics close to barrier energies.

\section{Quantum activation versus quantum tunneling}
\label{sec:tunneling}
The question whether switching between the wells happens via quantum activation only or to which extent direct low-energy quantum tunneling between the wells (ground state tunneling between period-3 orbits) may play a role, has recently been addressed in \cite{zhangDykman2019}. Here, we recall the main results to keep the line of reasoning self-contained.
Accordingly, we consider also the dominant exponential factor for coherent tunneling from one well to the others (corresponding to the maximal  level splitting between adjacent energy levels of $\hat{g}$). In a semiclassical limit, the latter is given by $\propto e^{-2S_{tun}/\Lambda}$, where up to corrections of order $\Lambda$, the action at a given quasi-energy $g<g_{\rm saddle}$ reads
\begin{equation}
S_{\rm tun}(g)=\int_{Q_1}^{-Q_1/2} {\rm Im}\{P(Q)\} dQ\,
\end{equation}
Note that since the integration contour includes not only areas, where $g<g(Q,P)$, the momentum carries apart from imaginary also real parts. Accordingly, the full action has also a real part that induces a phase dependence of the energy level splittings, for details see \cite{ZhangGosnerDykman2017}.

Now, the probability density at a given $g$ for times short compared to interwell processes but long compared to local relaxation is $\sim {\rm e}^{-R(g)/\Lambda}$, where $R(g)$ monotonously increases with $g$ from $R(g_{\rm min})=0$.
Tunneling can now happen from the bottom of the well or at a higher quasi-energy $g$. If the product $e^{-2S_{tun}(g)/\Lambda} \, {\rm e}^{-R(g)/\Lambda}$ monotonously grows with $g$, this implies that the system has the tendency to tunnel towards increasing $g \to g_{\rm saddle}$, and eventually will go over the barrier.
This phenomenon has been termed quantum activation. The condition for quantum activation is thus
\begin{equation}
2 \partial_g S_{\rm tun}(g)+\partial_g R(g)<0 \ \ \mbox{for}\ \ g_{\rm min}\leq g\leq g_{\rm saddle}\, .
\end{equation}
Since the derivative of the tunnel action provides the (dimensionless) tunneling time
\begin{equation}
\tau(g)\equiv \partial_g S_{\rm tun}(g)= {\rm Im}\int \frac{dQ}{\partial_P g(Q,P)}\, ,
\end{equation}
the above condition can also be written as $\partial_g R(g)< 2|\tau(g)|$ with $\tau(g)<0$ in the respective energy range.  As Fig.~\ref{tunnell} reveals, this condition is always fulfilled so that we can conclude that the relaxation dynamics  towards a steady state  with rate $\Gamma_L$, as for example  seen in Fig.~\ref{Wignerlindz}, occurs dominantly via quantum activation.

In essence, the picture that emerges from these and the previous sections is the following: Classically, the switching between period-3 states occurs via thermally activated processes, where a trajectory diffuses up in energy until it approaches the saddle points, from where it switches to another period-3 state via the period-1 oscillator as an intermediate state (for a detailed discussion see also \cite{zhangDykman2019}). Quantum mechanically, for sufficiently large friction, the system diffusively climbs up the quasi-energy ladder within a period-3 well until it approaches the energy range around the saddle point energy $g_{\rm saddle}$. In this energy range, direct tunneling through the quasi-energy barrier towards one of the adjacent period-3 states occurs.  The period-1 oscillator looses its role as an intermediate state since tunneling between period-3 and period-1 states is less likely for energetic reasons. However, it is still an open question how in detail the tunneling  between the period-3 states near the saddle points happens, a problem that we intend to consider elsewhere.

	 \begin{figure}[htp]
\centering
\includegraphics[width=1\linewidth]{Fig9.pdf}
 \caption{Condition for quantum activation: The tunneling time $|\tau(g)|$ exceeds the slope of the action $R(g)$ for all energies $g$ at $f=0.66$ [$g_{\rm min}=-3.55, g_{\rm saddle}=0.22$] so that the condition is always fulfilled. Consequently, the relaxation dynamics towards a steady state is dominated by quantum activation.}
 \label{tunnell}
\end{figure}

\section{Dissipative induced phase transition}
\label{sec:transition}

	 \begin{figure*}[htb]
\includegraphics[width=1\linewidth]{Fig10.pdf}
 \caption{(a)-(d): Steady state Wigner density $\rho_W^{(0)}(Q,P)$ of $\ket{\ket{\nu_0}}$; (e)-(h): Wigner density $\rho_W^{(3,c)}(Q,P)$ of $\ket{\ket{\nu_3}}_{c}$ for various driving strengths $f=0.27$, $f=0.47$, $f=0.6$, $f=0.67$ (from (a) to (d) and (e) to (h)) at $\bar{n}=0$, $\kappa=1.5$ and $\Lambda=0.1688$, where $f_c\approx 0.63$ [cf.~Fig.~\ref{realimageigliou}(a)]. $\rho_W^{(0)}(Q,P)$ makes a transition from a localized state (ground state) around $Q=P=0$  for small driving to a delocalized distribution predominantly localized in the wells $Q_m, P_m, m=1, 2, 3$ corresponding to classical period-3 fixed points. $\rho_W^{(3,c)}(Q,P)$ exhibits the complementary behavior. }
 \label{wignerofphysexc}
\end{figure*}

 \begin{figure*}[htb]
\includegraphics[width=1\textwidth]{Fig11.pdf}
 \caption{Fock state representation of the steady state density $\langle n|\rho^{(3, c)}|m\rangle$ corresponding to the
 eigenvector $\ket{\ket{\nu_3}}_{c}$ orthogonal to the steady state for $f=0.27$, $f=0.47$, $f=0.6$, $f=0.67$ (from (a) to (d)) and $\kappa=1.5$, $\Lambda=0.1688$. The density has a dominant contribution in the first excited Fock state for small driving which turns into a dominant contribution in the ground state for strong driving.}
 \label{wignerofphysexch}
\end{figure*}

So far we have discussed the time scales for relaxation towards the steady state and its dependence on the driving strength. Now we turn to the steady state itself to explore to what extent the emergence of multiple steady state orbits at a critical driving in the classical regime yields signatures in the quantum realm. As we discussed in Sec.~\ref{escaperatelambda}, near the critical driving $f_c$, a new time scale appears that is associated with interwell processes. This coincides  with a fundamental change in the nature of the steady state from a localized to a delocalized one. We recall that for finite $\Lambda$ the transition is smeared out due to fluctuations.

To analyse this changeover in more detail, we come back to the eigenvectors $\ket{\ket{\nu_0}}$ and $\ket{\ket{\nu_3}}$ corresponding to the steady state density $\rho_{ss}\equiv \rho^{(0)}$ and to the eigen-density $\rho^{(3)}$, respectively. While the former is asymptotically approached for long times, the latter includes those populations in Fock space contained in the full time-dependent density that survive the longest.
It turns out that the above transition leaves direct signatures in the Wigner and the Fock state representations of these two densities. As already mentioned above, the associated eigenvectors are not orthogonal. Consequently, in order to distill that part of $\ket{\ket{\nu_3}}$ that has no overlap with $\ket{\ket{\nu_0}}$, we consider
\begin{equation}
\ket{\ket{\nu_3}}_c= N_3\left(\ket{\ket{\nu_3}}-\bra{\bra{\nu_0}}\ket{\ket{\nu_3}}\, \ket{\ket{\nu_0}}\right)\,
\end{equation}
with normalization $N_3$. Of course, $\ket{\ket{\nu_3}}_c$ is also an eigenvector of $\mathcal{L}$ with eigenvalue $\lambda_3$.

Then, the corresponding density $\rho^{(3)}_c$ provides information about
those populations in Fock space, complementary to those in $\rho^{(0)}$,  that disappear with the least negative eigenvalue $\lambda_3$ during the relaxation process of the full density matrix $\rho(t)$ (the densities associated with $\lambda_{1, 2}$ carry off-diagonal elements only). One has to keep in mind though that it is not directly related to a physical density, but rather serves here as a diagnostic tool to analyse the (approximate) phase transition.

Figure \ref{wignerofphysexc} illustrates how the corresponding Wigner densities evolve with increasing driving. While the steady state displays the expected transition from localized to delocalized (from strong overlap with the oscillator ground state at $Q=P=0$ to strong overlap with those around $Q_m, P_m; m>0$), the density $\rho^{(3)}_c$ shows the opposite behavior, from delocalized to localized. Accordingly, the period-1 oscillator at the origin turns from a stationary into a transient state, while period-3 oscillators are associated with steady states beyond the critical driving $f_c\approx 0.63$.

The same behavior can be more precisely detected in the Fock state representation, see Fig. \ref{wignerofphysexch}. Below $f_c$ the distribution to $\rho^{(3)}_c$ occupies predominantly the {\em first excited} state of the oscillator at $Q=P=0$, i.e.\ $\rho^{(3)}_c\sim |1\rangle \langle1|$,  while beyond $f_c$ the {\em ground} state dominates, i.e.\ $\rho^{(3)}_c\sim |0\rangle \langle 0|$, and the excited state is absent. At that value of the driving, where this transition happens, i.e.\ $f_c$, we find empirically that  $|\lambda_3|$ takes its minimal value [see Fig.~\ref{realimageigliou}(a)]. This resembles behavior such as the closing of energy gaps and the emergence of slow modes that can be observed in phase transitions. For the ranges of parameters explored here, particularly small $\Lambda$, and $\kappa>1$, this coincides with a minimum in $|\lambda_3|$.
This still holds for $\Lambda\sim O(1)$ and stronger friction.
 Outside of these domains, one observes a shift of  the critical driving towards smaller values (not shown).

	 \begin{figure}[htp]
\centering
\includegraphics[width=1\linewidth]{Fig12.pdf}
 \caption{
The overlap $|\bra{\bra{\nu_0}}\ket{\ket{\nu_3}}|$ for $\kappa=1.5$, $\Lambda=0.1688$ (blue line) exhibits a maximum followed by a sharp drop to almost zero when $\ket{\ket{\nu_3}}$ becomes least unstable ($|\lambda_3|$ has a minimum); it sharply increases again with growing driving. Most of this behavior results from the components of the eigenvectors $\ket{\ket{\nu_0}}$ and $\ket{\ket{\nu_3}}$ in the sub-space spanned by the two lowest lying  Fock states $\{|0\rangle, |1\rangle\}$ ($|\bra{\bra{\nu_{0,r}}}\ket{\ket{\nu_{3,r}}}|$, black line).}
\label{realimageiglio}
\end{figure}

In terms of the original (non-orthogonal) eigenvectors, somewhat below the critical driving strength the overlap $|\bra{\bra{\nu_0}}\ket{\ket{\nu_3}}|$ (see Fig.~\ref{realimageiglio} blue line) exhibits a maximum followed by a sharp drop to almost zero when $\ket{\ket{\nu_3}}$ becomes least unstable ($|\lambda_3|$ has a minimum); it sharply increases again with growing driving. Most of this behavior results indeed from the components of the eigenvectors $\ket{\ket{\nu_0}}$ and $\ket{\ket{\nu_3}}$ in the sub-space spanned by the two lowest lying  Fock states $\{|0\rangle, |1\rangle\}$, i.e.\ $\ket{\ket{\nu_{0, r}}}$, $\ket{\ket{\nu_{3, r}}}$, see Fig.~\ref{realimageiglio} black line.

How the location of the minimum of $|\lambda_{3}|$ changes with the driving $f$ and the quantum scale $\Lambda$ is depicted in Fig.~\ref{changeinlambda}. For $\Lambda\to 0$, one  approaches the classical threshold $f_{c,cl}\approx 0.49$ (cf.~Fig.~\ref{bifdiag}), where the classical limit cannot be fully resolved (cf.~Fig.~\ref{changeinlambda}). With growing impact of quantum fluctuations $f_c(\Lambda)$ increases. Fig.~\ref{changeinlambda} also shows how $\text{min}(|\rm Re(\lambda_3)|)$ diminishes as the quantum fluctuations $\Lambda$ go to zero and the gap closes.

\begin{figure}
\centering
\includegraphics[width=1\linewidth]{Fig13.pdf}
 \caption{Critical driving $f_c$ (blue line) located at the minimum of the eigenvalue $|\lambda_3|$ with changing $\Lambda$ ($\kappa=1.5$ and $\bar{n}=0$). For decreasing $\Lambda$ the systems turns into the semiclassical regime and approaches the classical threshold $f_{c,cl}\approx 0.49$ for $\Lambda\to 0$. $\text{min}(|\rm Re(\lambda_3)|)$ (green line) decreases as $\Lambda$ goes to zero.}
 \label{changeinlambda}
\end{figure}

 \begin{figure}
\centering
\includegraphics[width=1\linewidth]{Fig14.pdf}
 \caption{Phase diagram (de-tuning versus driving amplitude) of a dissipative quantum oscillator driven by three times its fundamental frequency. The black line separates the classical regions of single (left) and multiple (right) solutions. Quantum mechanically, the green ($\bar{n}=0$)/ blue ($\bar{n}=0.1$) line separates the domain, where the steady state is localized around the origin of the $QP$-plane (left) from the domain, where it is localized with dominant contributions around classical period-3 fixed points (right). In comparison to the classical prediction for the bifurcation line in absence of noise which separates the regime of period-1 orbits from the one of period-3 orbits, the quantum system exhibits ranges with a re-entrant behavior when increasing the de-tuning for fixed driving. Parameters are $\gamma/\omega_0=0.01$, $\alpha q_0^2/\omega_0^2=0.0015$ with the new quantum parameter $q_0^2=\hbar/(2M\omega_0)=1/2$. The red point corresponds to the transition point in Fig.~\ref{bifdiag} ($f_{c, \rm cl}= 0.49$).  }
 \label{quphaset}
\end{figure}
One unique feature of periodically driven oscillators lies in the fact that they allow exploring dissipative phase transitions in steady state within set-ups which are easily tunable by the parameters of the external drive. In our situation, the corresponding parameter space is the two-dimensional space $(\delta\omega, F_0)$. In order to characterize in this space the domains of localized and delocalized steady states, respectively, and thus to provide predictions for experimental realizations, one has to work with dimensionless quantities where the scaling does {\em not} include $\delta\omega$, as we have done so far. Accordingly, we here switch to another scaling and consider physical parameters.
Figure \ref{quphaset} displays the corresponding phase diagram of the quantum system in the parameter space $([\omega_F-3\omega_0]/\omega_0, F_0/\omega_0^2)$.
The green/ blue solid line separates for vanishing/ finite temperature the domain, where the steady state is localized (left area), from the domain, where it is delocalized with most of the weight sitting in the three well regions of the quasi-energy surface (right area). Apparently, tuning $[\omega_F-3\omega_0]/\omega_0$ for fixed driving $F_0/\omega_0^2$, one finds ranges for the driving, where a re-entrant behavior can be seen: one starts at a given $F_0/\omega_0^2$ for low de-tuning from a localized state that turns into  a delocalized one for growing $[\omega_F-3\omega_0]/\omega_0$, but ends again in a localized state for large de-tunings. In comparison we show the classical situation at $\text{T}=0$, where the re-entrant property is absent (black line).

\section{Conclusion}

In this work we studied quantum-noise induced switching between stable states of an oscillator that displays stationary period-three vibrations. The oscillator's quasi-energy consists of three localized wells in phase space when a sufficiently strong driving is applied. Starting initially from a locally relaxed state in one of the wells in phase space, transitions to the other wells occur via quantum activation rather than through quantum tunneling. The relaxation dynamics is characterized by a single time scale that separates interwell from intrawell processes. A standard semiclassical treatment does not fully capture the complexity of the dynamics and, thus, can at best yield the order of magnitude for this time scale. It also reveals a breaking of detailed balance for all temperatures down to $\text{T}=0$.
The appearance of a time scale separation is associated with an approximate closing of the gap between the lowest lying eigenvalues for the relaxation of populations and in turn a dissipation-induced transition in the nature of the stationary state from being localized to being delocalized. In contrast to the classical situation, this phase transition exhibits a re-entrant behavior in the parameter space of the external drive (amplitude, frequency). These findings may further stimulate on-going experimental investigations based on superconducting circuits.

\section*{acknowledgments}
We are grateful for financial support by the German Science Foundation through DFG AN336/11-1 and the Center for Integrated Quantum Science and Technology (IQST). We acknowledge fruitful discussions with M. I. Dykman and Y. Zhang.

\clearpage