# I-CTRL: Imitation to Control Humanoid Robots Through Constrained
Reinforcement Learning

Yashuai Yan1, Esteve Valls Mascaro1, Tobias Egle2, Dongheui Lee1,3

###### Abstract

This paper addresses the critical need for refining robot motions that,
despite achieving a high visual similarity through human-to-humanoid
retargeting methods, fall short of practical execution in the physical realm.
Existing techniques in the graphics community often prioritize visual fidelity
over physics-based feasibility, posing a significant challenge for deploying
bipedal systems in practical applications. Our research introduces a
constrained reinforcement learning algorithm to produce physics-based high-
quality motion imitation onto legged humanoid robots that enhance motion
resemblance while successfully following the reference human trajectory. We
name our framework: I-CTRL. By reformulating the motion imitation problem as a
constrained refinement over non-physics-based retargeted motions, our
framework excels in motion imitation with simple and unique rewards that
generalize across four robots. Moreover, our framework can follow large-scale
motion datasets with a unique RL agent. The proposed approach signifies a
crucial step forward in advancing the control of bipedal robots, emphasizing
the importance of aligning visual and physical realism for successful motion
imitation.

Figure 1: Overview of our task of mimicking dynamic human motions into four
bipedal humanoid robots with constrained goal-conditioned reinforcement
learning. Each block shows diverse behaviors for different robots: Bruce,
Atlas, H1, and JVRC-1.

††footnotetext: 1Yashuai Yan, Esteve Valls Mascaro and Dongheui Lee are with
Autonomous Systems Lab, Technische Universität Wien (TU Wien), Vienna, Austria
(e-mail: {yashuai.yan, esteve.valls.mascaro,
dongheui.lee}@tuwien.ac.at).††footnotetext: 2Tobias Egle is with Robotics
Systems Lab, Technische Universität Wien (TU Wien), Vienna, Austria (e-mail:
{tobias.egle}@tuwien.ac.at).††footnotetext: 3Dongheui Lee is also with the
Institute of Robotics and Mechatronics (DLR), German Aerospace Center,
Wessling, Germany.

##  I Introduction

The task of endowing humanoid robots with human-like motion capabilities has
been a longstanding challenge in robotics, driven by the vision of integrating
machines into populated environments. A critical bottleneck in achieving this
vision lies in the ability to impart dynamic movements to robots that resemble
the robustness and sophistication inherent in human motor skills. This paper
addresses this challenge through constrained deep reinforcement learning to
control bipedal robots to accurately follow a reference human motion.

Recently, the animation and robotics communities have increasingly focused on
translating human motions into various agents, each with its distinct focus.
Motion imitation in animation has mainly addressed retargeting to preserve the
style and nuances of dynamics of human motions, while roboticists have focused
on the physically plausible executions in real robots. Given the inherent
difference in the tasks, two divergent approaches were taken by each
community. On one hand, in the animation community, deep-learning-based models
have been widely adopted lately to accurately translate human motions into
simulated characters [22, 36, 41, 1, 39]. However, most works in animation
often do not satisfy the constraints existing in our physical world: foot
penetrations, joint-limits violations, or instability, among others. On the
other hand, control-based optimization algorithms [26, 6, 34, 27, 42] have
been used by the robots community to ensure the proper robot control, but
either fail to generalize to different kinematics and scale to motion
diversities or do not correctly preserve the style with high quality.

Thanks to the recent advancement in DRL, some works [28, 29, 24] are now
making strides to ensure physically-based motions to animated characters.
Early works [8, 13] showed the potential of DRL to produce complex behaviors,
but exhibit numerous artifacts. Later, Peng et al. [28, 29] aimed to improve
the quality of the produced motions by reformulating the problem through goal-
conditioned RL [30], where the task is to follow a goal, considered in
imitation as a reference motion, to guide the RL agent. Finally, Luo et al.
[24] proposed a unique agent for physically-based avatar imitation learned
over thousands of human motions, thus opening the possibility of
generalization in motion imitation using RL. Despite the great evolution in
recent years, these works [28, 29, 24] rely on human data that identically
matches the kinematic of their simulated avatar, while not constrained to the
different mass distributions, joint limits, and kinematic structures of
robots. In the field of robotics, DRL has been applied mainly to quadrupedal
robots, which led to significant improvements in the agility and robustness of
the produced motions [2, 10], but only a few works apply DRL for bipedal
humanoid robots [19, 21]. Still, Li et al. [19, 21] trained a DRL agent to
accomplish specific tasks (walk, jump, run), but without generalization to
several different motions. We believe that quadrupedal robots are less
susceptible to losing balance due to different mass distributions than bipedal
robots. This complicates the physics-based motion imitation problem in legged
humanoid robots.

In this work, we demonstrate how the less stable nature inherent in bipedal
robots can be tackled efficiently by Imitating humans through ConsTrained deep
Reinforcement Learning, named I-CTRL. First, we reformulate the problem of
motion imitation as a goal-conditioned RL task that minimally refines a
reference motion produced by ImitationNet [39]. ImitationNet serves as the
baseline of our work as it learns to retarget human poses into robots using
unsupervised deep learning. However, we observed that the motions produced by
Yan et al. [39] for bipedal robots are often infeasible and the robots fall in
the real world. Therefore, I-CTRL refines the reference motion from Yan et al.
[39] to ensure feasible robot control and scalability to a large diversity of
motions. Additionally, we propose a general training policy that allows to
control of four bipedal humanoid robots without tuning the reward signals:
Bruce, Unitree H1, Atlas, or JVRC-1. The versatility of I-CTRL is crucial to
improving the rapid adaptation to new robots. Finally, to overcome the high-
constrained and complex nature of the different kinematics while ensuring
efficient training, we adapt contained RL to reduce the exploration space and
effectively preserve the motion style. Overall, our physics-based constrained
RL algorithm, named I-CTRL, produces highly accurate human motion imitation
onto various bipedal robots that can follow around 10.000 different dynamic
motions with a single learned policy. Our efforts resulted in the following
contributions.

  1. 1.

Learning physics-based human-like motions for various bipedal humanoid robots.

  2. 2.

Designing a novel sample-efficient constrained reinforcement learning
algorithm that better preserves the style of the targeted motion and
generalizes to around 10.000 motions using a single policy and shared rewards.

  3. 3.

Evaluating I-CTRL quantitatively and qualitatively on four bipedal robots and
a variety of motions to showcase the versatility of our model without reward
tuning.

##  II Related Work

Existing literature on motion imitation and robot control is reviewed next,
highlighting the importance of achieving physics-based motion to control
bipedal robots while preserving the style of reference human motion.

###  II-A Motion Imitation for Animated Characters

Producing highly realistic character animations that resemble human motions is
a persistent challenge in computer graphics. The goal is to translate human
motions onto numerous characters in a versatile, accurate, and generalizable
manner, that visually seems natural while being cost-efficient.

Historically, the first approaches handcrafted kinematic constraints to
properly map poses between the human and character domains. The emergence of
large-scale motion-capture (MoCap) datasets led to data-driven approaches that
could reduce this engineering problem. Initial works [7, 15] still required
pairing data between the source and target characters to build a shared
representation space that overcame the gap between kinematics encountered in
the retargeting process. To reduce the implementation cost and enhance
generalization to different agents, Villegas et al. [35] proposed an
unsupervised method with adversarial training and cycle consistency, which was
recently extended in [39] with a contrastive approach. Yan et al. [39] adopted
a synergy of different losses to ensure a high-quality unsupervised motion
retargeting system that is generalizable to various animated robots while
preserving the nuances of the source motion. However, the poses produced by
Yan et al. [39] do not ensure its plausibility in the physical world, and
therefore a bipedal robot could only be controlled for animation. Our goal is
to extend Yan et al. [39] to ensure that the highly accurate retarget motions
satisfy the physics laws of the real world.

The task of ensuring physically-based motions in the animated character has
recently gained followers as an additional condition to ensure realism.
Previous works [36, 41] adopted additional constraints for the deep learning
optimization to reduce interpenetrations in the character meshes. However, to
satisfy the plausibility of the motion, it is required to use a physics engine
that emulates the real world and allows the algorithms to extend animation to
reality. Inspired by the success of Deep Reinforcement Learning (DRL) in
numerous tasks [25, 32, 18], the animation field has adopted DRL showing
promising results in controlling physics-based animated characters [8, 23, 28,
29, 24]. Initial works [8, 13] tackled the physics-based motion synthesis from
an uncontrolled scenario, producing complex motion behaviors, but resulting in
unrealistic poses and gaits. To generate more human-like dynamics, Peng et al.
[28, 29] reformulated the task as a goal-conditioned DRL problem [30], where
human motion served as a reference to guide the DRL in learning. However, Peng
et al. [28] required manual tuning of the rewards for different motions, and
Peng et al. [29] trained a discriminant network over a small-scale
distribution of similar motions to act as a reward function, which highly
increased the complexity of the training. Recently, [23, 24] have proposed the
task of imitating large-scale motion datasets using a single policy. To
achieve a high success rate, Luo et al. [23] artificially added residual force
control [40] to prevent the animated character from falling and improve
balance. On the contrary, Luo et al. [24] adopted progressive RL [4] to permit
imitation of large-scale motions through a mixture of experts (MoE). Opposite
to Luo et al. [24], our work adopts a single network to follow a large
diversity of motions, which also results in simpler and more efficient
training. Additionally, the aforementioned works [28, 29, 24, 23] rely on
human data that identically matches the kinematic of their simulated avatar,
while not being challenged by the different mass distributions, joint limits,
and kinematic structures of robots. On the contrary, we propose a constrained
deep RL policy with simple rewards that generalize to various full-size
humanoid robots and large-scale motions from a generated reference motion from
Yan et al. [39]. Thanks to the reformulation of motion imitation as a
refinement task, I-CTRL can significantly improve the learning efficiency of
RL while preserving a high resemblance to the reference human motion.

###  II-B Model- and Learning-based Bipedal Locomotion

Whole-body control of bipedal robots is an ongoing challenge in the robotics
community. Historically, recurring motions such as bipedal walking and running
have been addressed using biologically inspired template models, with
particular emphasis on the linear inverted pendulum model (LIPM) [16] for
walking and the spring-loaded inverted pendulum [37] for running. Later, model
predictive control (MPC) [38, 14] showed improvements in generating stable
center-of-mass (COM) trajectories based on future footsteps. The current
state-of-the-art in robot control adopts divergent component of motion (DCM)
trajectory generation [33] with MPC-based step adaptation [17, 9] for the
locomotion of full-scale humanoid robots. However, those frameworks require
high configuration tuning and only focus on simple walking behaviors, without
allowing for versatile human motion imitation.

Lately, MPC-based approaches have been combined with deep learning to produce
robot motions that resemble humans or animals from casual videos [42].
However, Zhang et al. [42] required eight hours and eight GPUs to imitate a
single one-minute video, thus making it infeasible for real-time imitation and
generalization to large-scale motions. On the contrary, deep RL has already
shown promising potential in producing agile and robust locomotion for
quadrupedal robots [2, 10] which can work in real-time and deployed to outdoor
unconstrained environments. In regards to bipedal locomotion, Li et al. [19,
21] have recently adopted Deep RL to produce specific motions: walk, run, or
jump. To this end, Li et al. [21] archived real-world robustness through
multi-goal training, but only considered a few motions as goals to follow. On
the contrary, we propose imitating large-scale diverse human motions for
various bipedal robots, showing a promising direction to scalability and
generalization of the DeepRL for robot control.

###  II-C Constrained Reinforcement Learning

Robot control is typically formulated as an optimization task. On one hand,
model-based control roboticists typically use knowledge of the system dynamics
to create a mathematical model, set objective functions, and then use
optimization techniques to find control inputs that achieve those objectives
while adhering to system constraints. These constraints ensure feasibility and
safety but also shape the behavior of the system by restricting the solution
search space to the region that aligns with the engineer’s intent. DRL
reformulates the problem so that the engineer designs the objectives to
achieve a task, named rewards, and optimizes a neural network through trial
and error to maximize these rewards. Inspired by Constrained Markov Decision
Processes (CMDPs) [3], constrained DRL (CDRL) aims to find a policy that
maximizes the rewards while satisfying the constraints. CDRL has been used to
enforce velocities further from limits or avoid self-collisions [11], to
promote safety when the agent experiments in the real world, or to limit
feasible regions of exploration [11]. In this paper, we use goal-conditioned
DRL with constraints to impose that the generated robot motion is not far from
the reference motion, thus encouraging better alignment.

##  III Methodology

Figure 2: Overview of our I-CTRL framework. First, we train ImitationNet [39]
to retarget human motions into different robots in an unsupervised manner.
[39] provides an accurate robot pose $\mathbf{q}^{ref}_{t}$ that resembles the
original human but is not feasible for executing in the physical world. Then,
we train our constrained reinforcement learning agent to refine the next
reference pose $\mathbf{q}^{ref}_{t+1}$ based on our current robot joints
$\mathbf{q}_{t}$. IsaacGym is used as the physical simulation.

In this section, we present our whole-body imitation algorithm for humanoid
robots. An overview of our framework can be observed in Fig. 2

###  III-A Problem Formulation

Let $\mathbf{h}_{t}\in\mathbb{R}^{J\times n}$ be a human pose at time $t$
composed by $J$ joints in local quaternion representation ($n=4$), and
$\mathbf{H}=[\mathbf{h}_{1},\cdots,\mathbf{h}_{T}]\in\mathbb{R}^{T\times
J\times n}$ a human motion. Similarly, $\mathbf{q}_{t}\in\mathbb{R}^{J\times
1}$ represents a robot pose described in joint angles. As a notation
convention, we use $\mathbf{q}^{ref}_{t}$ and $\mathbf{q}^{ref}_{1:T}$ to
denote the robot reference pose and trajectory, respectively. $\mathbf{q}_{t}$
denotes the robot pose observed from the physics simulation. In addition, we
define the robot’s state at time $t$ as $\mathbf{s}_{t}$, consisting of its
joint angles $\mathbf{q}_{t}$, joint velocities $\dot{\mathbf{q}}_{t}$, root
position $\mathbf{p}^{root}_{t}\in\mathbb{R}^{3}$ and orientation
$\boldsymbol{\theta}^{root}_{t}\in\mathbb{R}^{4}$ in quaternion and its linear
and angular velocity in
$\mathbf{v}_{t},\boldsymbol{\omega}_{t}\in\mathbb{R}^{3}$. In our work,
robots’ root frame is defined as their base link frame. The reference state
$\mathbf{s}_{t}^{ref}$ is identical to $\mathbf{s}_{t}$, with the only
difference being the absence of joint velocities in $\mathbf{s}_{t}^{ref}$.
The velocities of reference motions are computed through finite differences.
Overall, our motion imitation task consists of first translating human motions
$\mathbf{H}$ into robot reference motions $\mathbf{q}^{ref}_{1:T}$ and then
refining these reference motions to robot motions $\mathbf{q}_{1:T}$ that
satisfy the physical laws of the real world.

###  III-B Human-to-robot pose retargeting

We define our physics-based motion imitation as two steps. Our task initially
requires a retargeting process from human $\mathbf{h}_{t}$ to a robot pose
$\mathbf{q}^{ref}_{t}$ that we use as reference. The goal is to achieve
$\mathbf{q}^{ref}_{t}$ that highly resembles a human pose in appearance, thus
the reference captures the style and nuances of the motion, despite not being
physically possible. For that, we make use of Yan et al. [39], which is a
simple but effective deep-learning model that uses unsupervised learning to
translate a human pose to a robot pose with high accuracy. We adopt Yan et al.
[39] as it does not require robot data for the retargeting, which allows us to
later evaluate the generalization of our RL algorithm to various robots. We
train a model $f$ based on Yan et al. [39] for each robot, and the model is
applied to map a human pose $\mathbf{h}_{t}$ to a robot pose
$\mathbf{q}^{ref}_{t}$ ($f:\mathbf{h}_{t}\longmapsto\mathbf{q}^{ref}_{t}$).

###  III-C Constrained Reinforcement Learning

In standard RL, the agents learn to maximize the expected return through trial
and error, and the objective function is only shaped by reward functions. Most
RL approaches are usually not sample-efficient and very sensible to the
designed reward functions. To cope with those limitations, we propose to add
constraints to the DRL for the human-to-robot imitation task, inspired by the
Constrained Markov Decision Process (CMDP) [3]. A CMDP is defined as
$\left\langle\mathcal{S,A},T,R,C,\rho,\gamma\right\rangle$, where
$\mathcal{S,A}$ denote the state and action space, $T:\mathcal{S\times
A}\rightarrow\mathcal{S}$ is the transition function and $R,C$ are reward and
cost functions. Lastly, $\gamma\in(0,1)$ and $\rho$ are the discount rate and
distribution of the initial state. As an extension to the Markov Decision
Process (MDP), the CMDP maximizes the expected discounted return while
satisfying certain constraints. Let

$J(\pi)=\mathop{\mathbb{E}}_{\rho,\pi,T}\left[\sum_{t=0}^{\infty}\gamma^{t}R(s_{t},a_{t},s_{t+1})\right]$
(1)

denote the expected return. To this end, we formulate the imitation task as an
optimization task,

$\begin{gathered}\underset{x}{\text{maximize}}\;\;\;J(\pi)\\\ \text{s.
t.}\;\;\;\;D(T(s_{t},\pi(s_{t})),f(h_{t+1}))<\delta\end{gathered}$ (2)

where $D(\cdot,\cdot)$ indicates the distance between two states, $\pi$ is the
policy to optimize. and $f$ is any function that maps human motions into state
space $\mathcal{S}$. As discussed in Sec. III-B, ImitationNet [39] is employed
as $f$ in this work due to the high resemblance of retargeted motions.
$\delta$ constrains the exploration region for searching optimal solutions and
serves as an upper limit on the allowable deviation between robot and human
motion. The robot states and the transition function $T$ are determined by the
real-world environment or a physics simulator while our policy $\pi$ computes
per-step action $\mathbf{a}_{t}\in\mathcal{A}$ conditioned on the state
$\mathbf{s}_{t}$ and the reference $\mathbf{s}_{t}^{ref}$. The action
$\mathbf{a}_{t}$ specifies target joint positions for PD controllers at each
joint. The reward function $R$ computing per-step reward
$r_{t}=R(\mathbf{s}_{t},\mathbf{a}_{t},\mathbf{s}_{t+1})$ is defined in
Sec.III-D. Similar to prior works [28, 29, 20, 24], the proximal policy
optimization [31] is used to learn our imitation policy $\pi$.

###  III-D Rewards

The reward functions define the optimization objective in reinforcement
learning. They directly influence the behavior of the agent and have a
significant impact on the agent’s performance. In our scenario, we aim to
design simple but general task rewards that can be applied to any robot
without finetuning. First, we adopt the task reward $r_{t}^{task}$ shown in
Eq. 3 to encourage the robot’s root to follow the trajectory of the reference
motion. Here, $\mathbf{p}_{t}^{root},\boldsymbol{\theta}_{t}^{root}$ and
$\mathbf{p}_{t}^{ref},\boldsymbol{\theta}_{t}^{ref}$ denote the position and
orientation of the robot’s root at time $t$ for robot and reference motion,
respectively.
$\boldsymbol{\theta}_{t}^{root}\ominus\boldsymbol{\theta}_{t}^{ref}$ denotes
distance between two quaternions.

$\displaystyle r_{t}^{task}=$ $\displaystyle
0.5\cdot\exp(-10\cdot||\mathbf{p}_{t}^{root}-\mathbf{p}_{t}^{ref}||)+$
$\displaystyle
0.5\cdot\exp(-10\cdot||\boldsymbol{\theta}_{t}^{root}\ominus\boldsymbol{\theta}_{t}^{ref}||)$
(3)

Secondly, we use the style reward $r_{t}^{style}$ defined in Eq. 4 to push the
imitated motion to the reference one by penalizing the difference of joint
angles between imitated and reference motions.

$r_{t}^{style}=\exp(-||\mathbf{q}_{t}-\mathbf{q}_{t}^{ref}||)$ (4)

Thirdly, we apply a smooth reward $r_{t}^{sm}$ as shown in Eq. 5 to penalize
high joint accelerations. This reward helps to avoid jitter motions during
imitation.

$r_{t}^{sm}=\exp(-0.1\cdot||\ddot{\mathbf{q}}_{t}||)$ (5)

Finally, we adopt a reward function $r_{t}^{c}$ to encourage robots’ feet to
follow the status of contacting the floor as in human motions.
$\mathbf{c}_{t},\mathbf{c}_{t}^{ref}\in\\{0,1\\}^{2}$, describing the contact
status, are binary values of two feet: $0$ for no contact, $1$ for contact.
For instance, $\mathbf{c}_{t}^{ref}=\left[1,0\right]^{T}$ indicates that the
left foot is touching the floor while the right foot is not. $r_{t}^{c}$ is
shown in Eq. 6. While $r_{t}^{c}$ helps enforce synchronization in the
stepping phase between human reference and robot, we only use $\mathbf{c}_{t}$
as an auxiliary reward but not in the robot state $\mathbf{s}_{t}$ defined in
Sec. III-A, to avoid depending on values which are difficult to estimate in
real scenarios.

$r_{t}^{c}=-||\mathbf{c}_{t}-\mathbf{c}_{t}^{ref}||$ (6)

Our total reward is a weighted sum of these four rewards, where
$w^{task}=1.0,w^{style}=1.0,w^{sm}=0.2,w^{c}=0.1$ are the weighting factors
presented in Eq. 7. All weighting values presented in each of the exponential
functions are chosen to normalize the average value of each reward, to later
on facilitate distributing the importance through the weighting factors. We
consider $w^{sm}$ and $w^{c}$ as soft auxiliary rewards to improve the
imitation performance, but with lower importance than $w^{task}$ and
$w^{style}$.

$r_{t}^{total}=w^{task}\cdot r_{t}^{task}+w^{style}\cdot
r_{t}^{style}+w^{sm}\cdot r_{t}^{sm}+w^{c}\cdot r_{t}^{c}$ (7)

Figure 3: Performance of our RL agent in translating a dynamic human dance
motion into four different humanoid robots. The depicted robots are, from top
to bottom, Bruce, H1, JVRC-1, and Atlas.

###  III-E Initialization and Early Termination

During the training, we use the first pose of the reference motion,
$\mathbf{q}_{1}^{ref}$, to initialize our robot joint angles in the physics
simulator. Similar to prior works [28, 24, 29], we also adopt the early
termination technique to boost our training process. We terminate the episode
when the robot’s root position $\mathbf{p}_{t}^{root}$ is more than 0.5 meters
away from the reference position $\mathbf{p}_{t}^{ref}$ or the root rotation
distance is more than 45 degrees.

##  IV Experiments

In this section, we evaluate the ability to control different bipedal humanoid
robots to imitate human motions.

###  IV-A Human MoCap Dataset

Following Yan et al. [39], we make use of the HumanML3D dataset [12] to train
our human-to-robot retargeting algorithm. HumanML3D consists of 14616 high-
quality human motions captured through motion capture (MoCap) systems with a
duration of 28.59 hours. As training dataset, we used 10,000 motions,
including _walking_ , _cooking_ , _punching_ , _dancing_ , _jumping_ ,
_running in a circle_ , _strong turns while dancing_ , etc. We excluded
motions that are unfeasible in physics-based simulators without the necessary
environmental data, such as _swimming_ without water or _walking up stairs_
without staircases. Also, we excluded highly dynamic motions, which are hard
to perform even for ordinary human, such as _parkour_ , _break-dance_ , or
_handstand_. Each motion in the training dataset is also annotated with
textual inputs, which allows for the connection of our work with text-to-
motion generators. Contrary to Yan et al. [39] that imitates humans’ upper
body only, we retarget the whole-body poses from human to humanoid robots.

###  IV-B Humanoid Robots

To demonstrate the generalization capabilities of our algorithm, we assess our
motion imitation system across four humanoid robots of varying sizes and
weights, and distinct configurations, eliminating the necessity for fine-
tuning rewards or specified values. Fig. 3 offers a qualitative view of our RL
approach, where each of the humanoid robots imitates the same human dancing
motion in the physics simulator.

  * •

Bruce is a small-sized humanoid robot with 16 actuated joints: 5 DoFs for each
leg and 3 DoFs for each arm.

  * •

H1 is a full-sized but lightweight (about 47 kg) humanoid robot with 19
actuated joints in total. H1 robots have 5 DoFs per leg and 4 DoFs per arm.
They also have 1 DoF for the trunk that allows it to rotate in the z-axis.

  * •

Atlas is a full-sized and heavy robot with 23 actuated joints in total. Atlas
has 6 DoFs for each leg and 4 DoFs for each arm, and it has 3 DoFs for the
body trunk.

  * •

JVRC-1 is a virtual humanoid robot created for the Japan Virtual Robotics
Challenge. Although there is no corresponding physical robot for it, JVRC-1 is
close to the real humanoid robots in terms of kinematic structure, link
masses, inertias, armature values, joint ranges, etc. The original JVRC-1
robot has 34 actuated joints. For our imitation task, we remove unused joints
(i.e. joints for the grippers and head rotation) to keep similar settings with
other robots. We actuate 23 joints for the legs, arms, and trunk, which is the
same as the Atlas robot.

As the sizes of different robots can vary a lot from the human model, we
normalize the human root trajectories by multiplying a scale factor on the
root’s position and linear velocity. The scaling factor is computed as the
ratio between the robot’s height and the human’s height.

###  IV-C Experimental setting

For reinforcement learning, both actor and critic are implemented as
feedforward neural networks with 8 fully connected linear layers. Each linear
layer has 256 neurons and is followed by an activation function Exponential
Linear Unit (ELU) [5]. We utilize the physics simulator NVIDIA’s Isaac Gym to
parallelize our training on an NVIDIA A4000 GPU. During the training, the
control policy is running at 60 Hz, while the simulation runs at 120 Hz. The
experimental settings are shared among different robots.

###  IV-D Evaluation

####  IV-D1 Metrics

We evaluate our model for human-like daily motions such as walking, running,
jumping, and arm movements. Next, we define four metrics used for our
evaluation:

  * •

Living Rate indicates whether a robot has fallen. Robots are considered to be
falling when their root height is less than a certain margin (e.g., 0.3 m).

$\displaystyle R_{t}^{living}=\begin{cases}1&\text{if }||z||>0.3\\\
0&\text{otherwise}\end{cases}$

  * •

Root Tracking Success Rate evaluates how well our agents are following the
reference root trajectory. The root tracking is considered to fail when
robot’s root deviates more than $D$ meters from the reference root position.
This metric is also used in [23, 24] for human character’s imitation with
$D=0.5$m. In this work we adapt $D$ depending on the different sizes of the
robots. For instance, given that $height_{Bruce}=0.4*height_{human}$, then
$D_{Bruce}=D*0.4=0.2m$.

$\displaystyle R_{t}^{tracking}=\begin{cases}1&\text{if
}||\mathbf{p}_{t}^{root}-\mathbf{p}_{t}^{ref}||<D\\\
0&\text{otherwise}\end{cases}$

  * •

Style Distance, similarly to our style reward, assesses the visual similarity
of imitated motions to reference ones. Style distance is defined as the Mean
Squared Error (MSE) of the $xyz$-euclidean positions of joints in meter units.
It seems a fair metric for motion resemblance, considering that ours and Peng
et al. [28] use a different representation of the style rewards (joint angles
versus body-link quaternions).

  * •

Root Position Distance, similarly to our task reward, showcases the ability to
follow the root trajectories in a global coordinate system. The Root Position
Distance is also considered as the MSE of the root position between the
reference and the robot motion.

  * •

Root Rotation Distance evaluates the ability of the robot to follow the
orientation of the human motion in a global coordinate system. The rotation
distance is computed as the quaternion distance between the orientation of the
robot’s root and the reference motion.

####  IV-D2 Generalization to various robots

We evaluate I-CTRL trained for a single-motion policy for Bruce, JVRC, H1, and
Atlas robots in the straight-line walking motion illustrated in Fig. 5 and
more dynamic dancing motions as shown in Fig. 3. We have reported the results
in Table I, which showcase the perfect success rate in all motions, with low
divergence from the reference motions both in terms of task and style.

|  |  Success Rate (%) $\uparrow$ |  Style Dist. (m) $\downarrow$ |  Root Pos. Dist. (m) $\downarrow$ |  Root Rot. Dist. (°) $\downarrow$
---|---|---|---|---|---
Robots | Policy | Walk | Dance | Walk | Dance | Walk | Dance | Walk | Dance
Bruce | DeepMimic, [28] | 72.07 | 97.94 | 0.031 | 0.050 | 0.136 | 0.112 | 0.557 | 2.51
AMP, [29] | 59.62 | 98.92 | 0.094 | 0.089 | 0.169 | 0.098 | 9.710 | 19.54
Ours | 100.00 | 100.00 | 0.010 | 0.049 | 0.036 | 0.039 | 4.263 | 1.036
JVRC | Ours | 100.00 | 100.00 | 0.047 | 0.073 | 0.034 | 0.037 | 0.855 | 0.952
H1 | Ours | 100.00 | 100.00 | 0.408 | 0.075 | 0.058 | 0.044 | 0.405 | 1.296
Atlas | Ours | 100.00 | 100.00 | 0.111 | 0.137 | 0.132 | 0.043 | 8.70 | 13.39

TABLE I: Quantitative evaluation of our different robots in the single policy
training. Here, bold denotes the best result per metric and motion.

####  IV-D3 Generalization of our policy

Following recent tendencies in the computer graphics community [24, 23], we
train a universal policy that can generalize to a large variety of different
human motions. To achieve this generalization, Luo et al. [23] adds residual
forces to avoid falling, which are inexistent in real robots, and Luo et al.
[24] considers a complex mixture of neural networks trained with high
parameter tuning to achieve generalization in humanoid characters. On the
contrary, we focus on a single network and policy with simple rewards and no
parameter tuning that can generalize to multiple complex robotic agents. We
evaluate the generalization of our policy in the Bruce, JVRC and H1 robots by
training our agent over 10.000 motions. The results presented in Table II
showcase that our RL agent can fulfill 89.64% of the 10.000 motions for Bruce,
88.70% for JVRC, and 76.69% for H1. Additionally, we compared our universal
models with the overfitted agents in Table I to illustrate the robustness and
high quality of the agents to dynamic motions such as walking and dancing.

####  IV-D4 Comparison with baselines

Due to the non-existence of a proper benchmark for human imitation into
robots, we train DeepMimic [28] and AMP [29] with the reference motions from
Yan et al. [39] for the Bruce robot. Given the one-policy-per-motion nature of
Peng et al. [28], we choose a walking motion and a dancing motion for the
evaluation. The reported results are indicated in Table I. We observe, as also
pointed out by the authors of [28], that DeepMimic is very sensible to the
reward weights compared to our model, which generalizes across robots and
motion variance without tuning. We conducted substantial hyperparameter tuning
for Peng et al. [28] to achieve high-quality results for the walking motion
and trained an agent with the same configuration for the dancing motion.
Still, the obtained performance of DeepMimic underperforms compared to our
behaviors. We also compared I-CTRL with [29], who designed a discriminator
network to predict a higher reward when it can not differentiate between the
reference and the produced motions, thus pushing the network to achieve
faithful motions. However, AMP underperforms for both walking and dancing
compared to I-CTRL as shown in Table I. We observe that AMP, as well as
DeepMimic, were trained originally with highly accurate and realistic motions.
In this case, the use of ImitationNet’s output causes a discrepancy of the
reference motion with the physics-based motion, which the discriminant of AMP
can easily identify. Some indicators that might ease the discriminant task are
foot penetrations or both feet not in contact with the ground.

| Bruce | JVRC | H1
---|---|---|---
Success Rate (%) $\uparrow$ | 89.64 | 88.70 | 76.69
Style Distance (m) $\downarrow$ | 0.029 | 0.056 | 0.080
Root Pos. Dist. (m) $\downarrow$ | 0.082 | 0.149 | 0.314
Root Rot. Dist. (°) $\downarrow$ | 3.427 | 2.502 | 18.122

TABLE II: Quantitative evaluation of our universal agent for a large variety
of human motions and three different bipedal robots. The living rate has not
been reported as it reaches 100% on both robots for the 10K motions.

$0$$0.5$$1$Living Rate (%) $\uparrow$ImitationNet [39]Ours (w/
constraints)Ours (w/o constraints)DeepMimic [28]AMP [29]$0$$0.5$$1$Root
Tracking Success Rate (%) $\uparrow$$0$$200$$400$$0$$5\cdot
10^{-2}$$0.1$$0.15$Training StepsStyle Distance (m)
$\downarrow$$0$$200$$400$$1\cdot 10^{-1}$$0.2$$0.3$Training StepsRoot Distance
(m) $\downarrow$

Figure 4: Quantitative evaluation of DeepMimic [28] and AMP [29] versus our
proposed approach, with and without constraints. Figure 5: Comparison of the
performance of DeepMimic [28] and AMP [29] versus our approach. Here, the
ground truth refers to the animated robot movement retargeted from a real
human walking using Yan et al. [39]. The ground-truth robot has the root
position fixed for visualization, and therefore no physics laws are applied.
Figure 6: Control of different bipedal humanoid robots for various everyday
motions.

####  IV-D5 Ablation Study

We evaluate the impact of constrained RL to enforce a high-quality imitation
in terms of style and efficiency in training. These results are shown in Fig.
4, where our constrained system converges faster in all metrics while
achieving lower style error.

####  IV-D6 Robustness to noisy reference

We address the question of the effectiveness of I-CTRL under low-quality
reference motions.

As shown in Sec. IV-D4, Peng et al. [28, 29] require a high-precision MoCAP
reference motion to achieve high-quality imitation, and underperform when the
reference motion is obtained from ImitationNet [39]. In order to evaluate the
robustness of I-CTRL under poor reference motions, we artificially add
Gaussian noise at different scales to the reference motion from ImitationNet.
Precisely, $\mathbf{\tilde{q}}_{1:T}^{ref}=\mathbf{q}_{1:T}^{ref}+G$ where
$G\sim\mathcal{N}(0,\,\sigma^{2})$. We report our quantitative results in
Table III for different standard deviations ($\sigma$) . As we expected, as
the noise level increases, it also augments the uncertainty of I-CTRL in the
task performance. However, Bruce robot is still able to successfully follow
the walking trajectory in 73.1% of the cases when the body has been affected
with an additional gaussian noise of $\sigma=5\degree$ on top of the errors of
the imitation algorithm.

Noise - $\sigma$ (°) | 0 | 5 | 10
---|---|---|---
Living Rate (%) $\uparrow$ | 100 | 98.0 | 86.5
Success Rate (%) $\uparrow$ | 91.1 | 73.1 | 60.6
Root Pos. Dist. (m) $\downarrow$ | $0.10\pm 0.06$ | $0.13\pm 0.10$ | $0.15\pm 0.11$
Root Rot. Dist. (°) $\downarrow$ | $0.90\pm 0.42$ | $1.50\pm 2.22$ | $3.01\pm 4.80$

TABLE III:  Quantitative evaluation of the robustness of our model in
following a walking forward motion given different noise rates in the
reference style motion for the Bruce robot.

###  IV-E Future work

Despite the high performance of our robotic agents in following the desired
motions, we observed that controlling heavyweight robots, such as Atlas, is
very challenging, as it requires manually tuning the stiffness and damping per
joint for all movements. Investigating an agent that produces torque control
directly based on the motion to follow will be considered in the future to
solve this manual tuning.

##  V Conclusion

In this paper, we introduced a novel constrained reinforcement learning
algorithm to imitate human motion on various legged humanoid robots, named
I-CTRL. We consider the imitation task as a two-stage problem: first, we
retarget human motions to robots and then refine those non-feasible motions
using I-CTRL. Therefore, we constrain the exploration space for the RL agent
to a smaller $\delta$-ball subspace, which improves motion resemblance and
training efficiency. Additionally, we focus on a versatile imitation system
that generalizes to different bipedal robots without reward function fine-
tuning, from small-sized Bruce to full-sized H1, and from lightweight JVRC-1
to heavyweight Atlas. Under consistent settings, our method enables these
robots to imitate complex dynamic motions such as walking, running, jumping,
or miscellaneous dances, among others. Furthermore, we expand the capabilities
of a single agent on Bruce, JVRC, and H1 robots to approximately 10,000
diverse dynamic motions and achieve an average of 84.9% success rate. The
average style and root distance to the reference motion are less than 3 cm and
8 cm, respectively, for the Bruce robot, showcasing that the produced physics-
based motions preserve visual similarity and replicate root trajectories with
high quality. These behaviors are shown through various quantitative and
qualitative experiments for different motions and robots.

## References

  * Aberman et al. [2020] Kfir Aberman, Peizhuo Li, Dani Lischinski, Olga Sorkine-Hornung, Daniel Cohen-Or, and Baoquan Chen.  Skeleton-aware networks for deep motion retargeting.  _ACM Transactions on Graphics (TOG)_ , 39(4):62–1, 2020.
  * Agarwal et al. [2023] Ananye Agarwal, Ashish Kumar, Jitendra Malik, and Deepak Pathak.  Legged locomotion in challenging terrains using egocentric vision.  In _Conference on Robot Learning_ , pages 403–415. PMLR, 2023.
  * Altman [1999] E. Altman.  _Constrained Markov Decision Processes_.  Chapman and Hall, 1999.
  * Bergamin et al. [2019] Kevin Bergamin, Simon Clavet, Daniel Holden, and James Richard Forbes.  Drecon: data-driven responsive control of physics-based characters.  _ACM Trans. Graph._ , 38(6), nov 2019.  ISSN 0730-0301.  doi: 10.1145/3355089.3356536.  URL https://doi.org/10.1145/3355089.3356536.
  * Clevert et al. [2016] Djork-Arné Clevert, Thomas Unterthiner, and Sepp Hochreiter.  Fast and accurate deep network learning by exponential linear units (elus), 2016.
  * Darvish et al. [2019] Kourosh Darvish, Yeshasvi Tirupachuri, Giulio Romualdi, Lorenzo Rapetti, Diego Ferigo, Francisco Javier Andrade Chavez, and Daniele Pucci.  Whole-body geometric retargeting for humanoid robots.  In _2019 IEEE-RAS 19th International Conference on Humanoid Robots (Humanoids)_ , pages 679–686. IEEE, 2019.
  * Delhaisse et al. [2017] Brian Delhaisse, Domingo Esteban, Leonel Rozo, and Darwin Caldwell.  Transfer learning of shared latent spaces between robots with similar kinematic structure.  02 2017.  doi: 10.1109/IJCNN.2017.7966379.
  * Duan et al. [2016] Yan Duan, Xi Chen, Rein Houthooft, John Schulman, and Pieter Abbeel.  Benchmarking deep reinforcement learning for continuous control.  In Maria Florina Balcan and Kilian Q. Weinberger, editors, _Proceedings of The 33rd International Conference on Machine Learning_ , volume 48 of _Proceedings of Machine Learning Research_ , pages 1329–1338, New York, New York, USA, 20–22 Jun 2016. PMLR.  URL https://proceedings.mlr.press/v48/duan16.html.
  * Egle et al. [2023] Tobias Egle, Johannes Englsberger, and Christian Ott.  Step and timing adaptation during online dcm trajectory generation for robust humanoid walking with double support phases.  In _2023 IEEE-RAS 22nd International Conference on Humanoid Robots (Humanoids)_ , pages 1–8, 2023.  doi: 10.1109/Humanoids57100.2023.10375152.
  * Fu et al. [2023] Zipeng Fu, Xuxin Cheng, and Deepak Pathak.  Deep whole-body control: Learning a unified policy for manipulation and locomotion.  In Karen Liu, Dana Kulic, and Jeff Ichnowski, editors, _Proceedings of The 6th Conference on Robot Learning_ , volume 205 of _Proceedings of Machine Learning Research_ , pages 138–149. PMLR, 14–18 Dec 2023.  URL https://proceedings.mlr.press/v205/fu23a.html.
  * Gangapurwala et al. [2020] Siddhant Gangapurwala, Alexander Mitchell, and Ioannis Havoutis.  Guided constrained policy optimization for dynamic quadrupedal robot locomotion.  _IEEE Robotics and Automation Letters_ , 5(2):3642–3649, 2020.  doi: 10.1109/LRA.2020.2979656.
  * Guo et al. [2022] Chuan Guo, Shihao Zou, Xinxin Zuo, Sen Wang, Wei Ji, Xingyu Li, and Li Cheng.  Generating diverse and natural 3d human motions from text.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)_ , pages 5152–5161, June 2022.
  * Heess et al. [2017] Nicolas Heess, Dhruva Tb, Srinivasan Sriram, Jay Lemmon, Josh Merel, Greg Wayne, Yuval Tassa, Tom Erez, Ziyu Wang, SM Eslami, et al.  Emergence of locomotion behaviours in rich environments.  _arXiv preprint arXiv:1707.02286_ , 2017.
  * Herdt et al. [2010] Andrei Herdt, Nicolas Perrin, and Pierre-Brice Wieber.  Walking without thinking about it.  In _2010 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , pages 190–195, 2010.  doi: 10.1109/IROS.2010.5654429.
  * Jang et al. [2018] Hanyoung Jang, Byungjun Kwon, Moonwon Yu, Seong Kim, and Jongmin Kim.  A variational u-net for motion retargeting.  pages 1–2, 12 2018.  doi: 10.1145/3283289.3283316.
  * Kajita et al. [2001] S. Kajita, F. Kanehiro, K. Kaneko, K. Yokoi, and H. Hirukawa.  The 3d linear inverted pendulum mode: a simple modeling for a biped walking pattern generation.  In _Proceedings 2001 IEEE/RSJ International Conference on Intelligent Robots and Systems. Expanding the Societal Role of Robotics in the the Next Millennium (Cat. No.01CH37180)_ , volume 1, pages 239–246 vol.1, 2001.  doi: 10.1109/IROS.2001.973365.
  * Khadiv et al. [2020] Majid Khadiv, Alexander Herzog, S. Ali. A. Moosavian, and Ludovic Righetti.  Walking control based on step timing adaptation.  _IEEE Transactions on Robotics_ , 36(3):629–643, 2020.  doi: 10.1109/TRO.2020.2982584.
  * Levine et al. [2018] Sergey Levine, Peter Pastor, Alex Krizhevsky, Julian Ibarz, and Deirdre Quillen.  Learning hand-eye coordination for robotic grasping with deep learning and large-scale data collection.  _The International journal of robotics research_ , 37(4-5):421–436, 2018.
  * Li et al. [2021] Zhongyu Li, Xuxin Cheng, Xue Bin Peng, Pieter Abbeel, Sergey Levine, Glen Berseth, and Koushil Sreenath.  Reinforcement learning for robust parameterized locomotion control of bipedal robots.  In _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , pages 2811–2817. IEEE, 2021.
  * Li et al. [2023a] Zhongyu Li, Xue Bin Peng, Pieter Abbeel, Sergey Levine, Glen Berseth, and Koushil Sreenath.  Robust and versatile bipedal jumping control through reinforcement learning, 2023a.
  * Li et al. [2023b] Zhongyu Li, Xue Bin Peng, Pieter Abbeel, Sergey Levine, Glen Berseth, and Koushil Sreenath.  Robust and versatile bipedal jumping control through multi-task reinforcement learning.  _arXiv preprint arXiv:2302.09450_ , 2023b.
  * Lim et al. [2019] Jongin Lim, Hyung Jin Chang, and Jin Young Choi.  Pmnet: learning of disentangled pose and movement for unsupervised motion retargeting.  In _Proceedings of the 30th British Machine Vision Conference (BMVC 2019)_. British Machine Vision Association, BMVA, September 2019.  30th British Machine Vision Conference (BMVC 2019) ; Conference date: 09-09-2019 Through 12-09-2019.
  * Luo et al. [2021] Zhengyi Luo, Ryo Hachiuma, Ye Yuan, and Kris Kitani.  Dynamics-regulated kinematic policy for egocentric pose estimation.  _Advances in Neural Information Processing Systems_ , 34:25019–25032, 2021.
  * Luo et al. [2023] Zhengyi Luo, Jinkun Cao, Kris Kitani, Weipeng Xu, et al.  Perpetual humanoid control for real-time simulated avatars.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 10895–10904, 2023.
  * Mnih et al. [2013] Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra, and Martin Riedmiller.  Playing atari with deep reinforcement learning.  _arXiv preprint arXiv:1312.5602_ , 2013.
  * Ott et al. [2008] Christian Ott, Dongheui Lee, and Yoshihiko Nakamura.  Motion capture based human motion recognition and imitation by direct marker control.  In _Humanoids 2008 - 8th IEEE-RAS International Conference on Humanoid Robots_ , pages 399–405, 2008.  doi: 10.1109/ICHR.2008.4755984.
  * Penco et al. [2018] Luigi Penco, B. Clement, V. Moduano, Enrico Mingo Hoffman, Gabriele Nava, Daniele Pucci, Nikos Tsagarakis, J. Mourert, and Serena Ivaldi.  Robust real-time whole-body motion retargeting from human to humanoid.  pages 425–432, 11 2018.  doi: 10.1109/HUMANOIDS.2018.8624943.
  * Peng et al. [2018] Xue Bin Peng, Pieter Abbeel, Sergey Levine, and Michiel van de Panne.  Deepmimic: Example-guided deep reinforcement learning of physics-based character skills.  _ACM Trans. Graph._ , 37(4):143:1–143:14, July 2018.  ISSN 0730-0301.  doi: 10.1145/3197517.3201311.  URL http://doi.acm.org/10.1145/3197517.3201311.
  * Peng et al. [2021] Xue Bin Peng, Ze Ma, Pieter Abbeel, Sergey Levine, and Angjoo Kanazawa.  Amp: Adversarial motion priors for stylized physics-based character control.  _ACM Trans. Graph._ , 40(4), July 2021.  doi: 10.1145/3450626.3459670.  URL http://doi.acm.org/10.1145/3450626.3459670.
  * Plappert et al. [2018] Matthias Plappert, Marcin Andrychowicz, Alex Ray, Bob McGrew, Bowen Baker, Glenn Powell, Jonas Schneider, Josh Tobin, Maciek Chociej, Peter Welinder, et al.  Multi-goal reinforcement learning: Challenging robotics environments and request for research.  _arXiv preprint arXiv:1802.09464_ , 2018.
  * Schulman et al. [2017] John Schulman, Filip Wolski, Prafulla Dhariwal, Alec Radford, and Oleg Klimov.  Proximal policy optimization algorithms, 2017.
  * Silver et al. [2016] David Silver, Aja Huang, Chris J Maddison, Arthur Guez, Laurent Sifre, George Van Den Driessche, Julian Schrittwieser, Ioannis Antonoglou, Veda Panneershelvam, Marc Lanctot, et al.  Mastering the game of go with deep neural networks and tree search.  _nature_ , 529(7587):484–489, 2016.
  * Takenaka et al. [2009] Toru Takenaka, Takashi Matsumoto, and Takahide Yoshiike.  Real time motion generation and control for biped robot -1st report: Walking gait pattern generation-.  In _2009 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , pages 1084–1091, 2009.  doi: 10.1109/IROS.2009.5354662.
  * Viceconte et al. [2022] Paolo Viceconte, Raffaello Camoriano, Giulio Romualdi, Diego Ferigo, Stefano Dafarra, Silvio Traversaro, Giuseppe Oriolo, Lorenzo Rosasco, and Daniele Pucci.  Adherent: Learning human-like trajectory generators for whole-body control of humanoid robots.  _IEEE Robotics and Automation Letters_ , PP:1–1, 01 2022.  doi: 10.1109/LRA.2022.3141658.
  * Villegas et al. [2018] Ruben Villegas, Jimei Yang, Duygu Ceylan, and Honglak Lee.  Neural kinematic networks for unsupervised motion retargetting.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 8639–8648, 2018.
  * Villegas et al. [2021] Ruben Villegas, Duygu Ceylan, Aaron Hertzmann, Jimei Yang, and Jun Saito.  Contact-aware retargeting of skinned motion.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 9720–9729, 2021.
  * Wensing and Orin [2013] Patrick M. Wensing and David E. Orin.  High-speed humanoid running through control with a 3d-slip model.  In _2013 IEEE/RSJ International Conference on Intelligent Robots and Systems_ , pages 5134–5140, 2013.  doi: 10.1109/IROS.2013.6697099.
  * Wieber [2006] Pierre-Brice Wieber.  Trajectory free linear model predictive control for stable walking in the presence of strong perturbations.  _2006 6th IEEE-RAS International Conference on Humanoid Robots_ , pages 137–142, 2006.  URL https://api.semanticscholar.org/CorpusID:2154754.
  * Yan et al. [2023] Yashuai Yan, Esteve Valls Mascaro, and Dongheui Lee.  Imitationnet: Unsupervised human-to-robot motion retargeting via shared latent space.  In _2023 IEEE-RAS 22nd International Conference on Humanoid Robots (Humanoids)_ , pages 1–8, 2023.  doi: 10.1109/Humanoids57100.2023.10375150.
  * Yuan and Kitani [2019] Ye Yuan and Kris Kitani.  Ego-pose estimation and forecasting as real-time pd control.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 10082–10092, 2019.
  * Zhang et al. [2023a] Jiaxu Zhang, Junwu Weng, Di Kang, Fang Zhao, Shaoli Huang, Xuefei Zhe, Linchao Bao, Ying Shan, Jue Wang, and Zhigang Tu.  Skinned motion retargeting with residual perception of motion semantics & geometry.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 13864–13872, 2023a.
  * Zhang et al. [2023b] John Z Zhang, Shuo Yang, Gengshan Yang, Arun L Bishop, Swaminathan Gurumurthy, Deva Ramanan, and Zachary Manchester.  Slomo: A general system for legged robot motion imitation from casual videos.  _IEEE Robotics and Automation Letters_ , 2023b.


# Quantum Discord for Multi-qubit Systems

Bo Li<EMAIL_ADDRESS>School of Mathematics and Computer science,
Shangrao Normal University, Shangrao 334001, China  Chen-Lu Zhu
<EMAIL_ADDRESS>Department of Mathematics, East China University of
Technology, Nanchang 330013, China  Xiao-Bin Liang<EMAIL_ADDRESS>School of Mathematics and Computer science, Shangrao Normal University,
Shangrao 334001, China  Biao-Liang Ye<EMAIL_ADDRESS>Quantum
Information Research Center, Shangrao Normal University, Shangrao 334001,
China  Shao-Ming Fei<EMAIL_ADDRESS>School of Mathematical Sciences,
Capital Normal University, Beijing 100048, China  Max-Planck-Institute for
Mathematics in the Sciences, 04103, Leipzig, Germany

###### Abstract

We evaluate analytically the quantum discord for a large family of multi-qubit
states. It is interesting to note that the quantum discord of three-qubits and
five-qubits is the same, as is the quantum discord of two-qubits and six-
qubits. We discover that the quantum discord of this family states can be
concluded into three categories. The level surfaces of the quantum discord in
the three categories is shown through images. Furthermore, we investigated the
dynamic behavior of quantum discord under decoherence. For the odd partite
systems, we prove the frozen phenomenon of quantum discord doesn’t exist under
the phase flip channel, while it can be found in the even partite systems.

###### pacs:

03.67.Mn, 03.65.Ud, 03.65.Yz

##  I Introduction

Quantum correlations are essential features of quantum mechanics which
distinguish the quantum from the classical world and play very important roles
in quantum information processing. The quantum correlated states are shown to
be more useful than the classically correlated ones in performing
communication and computation tasks. Understanding and quantifying various
quantum correlations are the primary goals in quantum information theory. The
quantum entanglement and nonlocal correlations can be considered the most
fundamental resources in quantum information processing Adesso ; Ma ;
Baumgratz ; Streltsov ; Streltsov1 ; Styliaris ; Xi ; Xu ; Xiong ; Huang ;
Oszmaniec , which are tightly related to quantum coherence Theurer ; Yao ; Yu
; Jin ; zhuhj ; xiya .

The quantum discord is of the most famous quantum correlations proposed by
Ollivier and Zurek Ollivier  and Henderson and Vedral Henderson , which
quantifies the quantum correlations in bipartite systems without quantum
entanglement. It is defined as the minimum difference between the quantum
versions of two classically equivalent expressions of mutual information under
projective measurements Rulli ; Luo . Due to the complexity of the
minimization process, the computation of quantum discord is a hard task and
seldom analytic results are know only for some restricted families of states
Ali ; Vinjanampathy ; Xie ; luo1 ; lang ; ali . For a bipartite state $\rho$
in systems $A$ and $B$Ollivier ; Henderson , the quantum discord
$D_{A;B}(\rho)$ is defined by
$D_{A;B}(\rho)=\min\limits_{\Pi^{A}}\Big{[}S_{B|\Pi^{A}}(\rho)-S_{B|A}(\rho)\Big{]},$
where the conditional entropy $S_{B|A}(\rho)=S(\rho)-S(\rho_{A})$ with
$S(X)=-TrX\log_{2}X$ the von Neumann entropy of a state $X$, $\rho_{A}$ is the
reduced state associated to the system $A$.
$S_{B|\Pi^{A}}(\rho)=\sum\limits_{j}p_{j}^{A}S(\Pi^{A}_{j}\rho\Pi^{A}_{j}/p_{j}^{A}),$
where $\Pi^{A}_{j}$ is the von Neumann projection operator on the subsystem A
and $p_{j}^{A}=Tr(\Pi^{A}_{j}\rho\Pi^{A}_{j})$ is the probability with respect
to the measurement outcome $j$.

Very recently, Radhakrishnan _et.al_ Radhakrishnan  introduced a
generalization of discord for tripartite and multipartite states. One of the
main features of this approach is the use of conditional measurements, where
each successive measurement is conditionally related to the previous
measurements. The $(N-1)$-partite measurement is written as

$\displaystyle\Pi_{j_{1}\cdots j_{N-1}}^{A_{1}\cdots
A_{N-1}}=\Pi_{j_{1}}^{A_{1}}\otimes\Pi_{j_{2}|j_{1}}^{A_{2}}\cdots\otimes\Pi_{j_{N-1}|j_{1}\cdots
j_{N-2}}^{A_{N-1}},$

where $\Pi_{j_{1}|j_{2}}^{A_{2}}$ is a projector on subsystem $A_{2}$
conditioned on the measurement outcome of $A_{1}$. Here the measurements take
place in the order $A_{1}\rightarrow A_{2}\rightarrow\cdots A_{N-1}.$ Such
conditioned measurements are essential to take into account all the classical
correlations that may exist among the subsystems. Viewing the measurements as
operations to break the quantum correlations, the optimization over all such
measurements allows one to recover the pure quantum contributions. Moreover,
there is an obvious asymmetry due the fixed ordering of the measurements. This
asymmetry has similarities with the quantum steering where one also considers
measurements on part of a system, while the aims are somewhat different in
that for discord, one minimizes the disturbance due to measurements rather
than compares it to a local hidden state theory. In deed, in some quantum
information processing such as one-way quantum computing, there is a definite
ordering of measurements, in compatible with the multipartite discord.

The quantum discord of an N-partite state $\rho$ is defined by

$\displaystyle D_{A_{1};A_{2};\dots;A_{N}}(\rho)=$
$\displaystyle\min_{\Pi^{A_{1}\dots A_{N-1}}}\Big{[}-S_{A_{2}\dots
A_{N}|A_{1}}(\rho)$
$\displaystyle+S_{A_{2}|\Pi^{A_{1}}}(\rho)\dots+S_{A_{N}|\Pi^{A_{1}\dots
A_{N-1}}}(\rho)\Big{]}$ (1)

for the measurement ordering $A_{1}\rightarrow A_{2}\rightarrow\dots A_{N-1}$.
Here we have defined $S_{A_{k}|\Pi^{A_{1}\dots
A_{k-1}}}(\rho)=\sum_{j_{1}\dots j_{k-1}}p_{j}^{(k-1)}S_{A_{1}\dots
A_{k}}(\Pi^{(k-1)}_{j}\rho\Pi^{(k-1)}_{j}/p_{j}^{(k-1)})$ with
$\Pi^{(k)}_{{j}}\equiv\Pi^{A_{1}\dots A_{k}}_{j_{1}\dots j_{k}}$,
$p_{j}^{(k)}=\text{Tr}(\Pi^{(k)}_{j}\rho\Pi^{(k)}_{j})$.

In general, it is difficult to evaluate the quantum discord (1) due to the
complexity of the optimization. We analyze and evaluate this quantum discord
for a family of multi-qubit states, and graphically show the level surfaces of
the quantum discord of this family. Moreover, due to the interaction with the
environments, bipartite quantum discord may decrease asymptomatically with
time Ban , and may be also frozen Mazzola ; Maziero ; Haikka ; Bromley  and
decoherence free for certain time. We also study the dynamic behavior of
quantum discord for a family of three-qubit and four-qubit states under
decoherence. We discover that the multi-qubit quantum discord of some states
cannot be destroyed by decoherence in finite time.

The rest of this article is organized as follows. In Sec. II, we calculate
analytically the multi-qubit discord for a family of quantum states. We shown
that the quantum discord can be classified into three categories. In Sec. III,
we investigated the dynamical behavior of the discord for a family of three-
qubit and four-qubit states. We discuss and summarize the results in Sec. IV.

##  II Quantum discord for multi-qubit systems

Consider the following family of $N$-qubit states,

$\rho=\frac{1}{2^{N}}(I+\sum\limits_{j=1}^{3}c_{j}\sigma_{j}\otimes\cdots\otimes\sigma_{j}),$
(2)

where $\sigma_{j}$, $j=1,2,3$, are the standard Pauli matrices, $I$ stands for
the corresponding identity operator. The motivation to consider the states (2)
is that, for $N=2$ (2) reduces to the well-known Bell diagonal states whose
famous analytical formulae of quantum discord have been provided by Luo luo1 ,
which attracted much attention and resulted in further vital results. For
general $N$, these states are highly symmetric and include some generalized
GHZ or W states as special ones.

We consider the family of three-qubit state, associated with systems $A$, $B$
and $C$,

$\rho=\frac{1}{8}(I+\sum\limits_{j=1}^{3}c_{j}\sigma_{j}\otimes\sigma_{j}\otimes\sigma_{j}).$
(3)

From (1) the quantum discord is given by

$\displaystyle D_{A;B;C}(\rho)=$
$\displaystyle\min_{\Pi^{AB}}\Big{[}-S_{BC|A}(\rho)+S_{B|\Pi^{A}}(\rho)+S_{C|\Pi^{AB}}(\rho)\Big{]}.$
(4)

Since $\rho_{A}=Tr_{BC}(\rho)=\frac{I}{2}$, we have the entropy
$S(\rho_{A})=1$. Set $\xi=\sqrt{c_{1}^{2}+c_{2}^{2}+c_{3}^{2}}$. One can
verify that

$S(\rho)=-4\times\frac{1+\xi}{8}\log_{2}\frac{1+\xi}{8}-4\times\frac{1-\xi}{8}\log_{2}\frac{1-\xi}{8}.$

Hence,

$\displaystyle-S_{BC|A}(\rho)$ $\displaystyle=-(S(\rho)-S(\rho_{A}))$
$\displaystyle=\frac{1+\xi}{2}\log_{2}(1+\xi)+\frac{1-\xi}{2}\log_{2}(1-\xi)-2.$
(5)

Denote $\\{\Pi_{k}=|k\rangle\langle k|:k=0,1\\}$. The von Neumann measurement
on subsystem $A$ is given by $\\{A_{k}=V_{A}\Pi_{k}V_{A}^{\dagger}:k=0,1\\}$,
where $V_{A}=t_{A}I+i\overrightarrow{y_{A}}\cdot\overrightarrow{\sigma}$ is
the unitary operator with $t_{A}\in\mathbb{R}$,
$\overrightarrow{y_{A}}=(y_{A1},y_{A2},y_{A3})\in\mathbb{R}^{3}$ and
$t_{A}^{2}+y_{A1}^{2}+y_{A2}^{2}+y_{A3}^{2}=1$. After the measurement
${A_{k}}$, the state $\rho$ is going to become the ensemble
$\\{\rho_{k},p_{k}\\}$ with $\rho_{k}:=\frac{1}{p_{k}}(A_{k}\otimes
I)\rho(A_{k}\otimes I)$ and $p_{k}=Tr(A_{k}\otimes I)\rho(A_{k}\otimes I).$
Then we obtain $p_{0}=p_{1}=\frac{1}{2}$,

$\displaystyle\rho_{0}=$
$\displaystyle\frac{1}{4}V_{A}\Pi_{0}V_{A}^{\dagger}\otimes(I\otimes
I+c_{1}z_{1}\sigma_{1}\otimes\sigma_{1}+c_{2}z_{2}\sigma_{2}\otimes\sigma_{2}$
$\displaystyle+c_{3}z_{3}\sigma_{3}\otimes\sigma_{3})$ (6)

and

$\displaystyle\rho_{1}=$
$\displaystyle\frac{1}{4}V_{A}\Pi_{1}V_{A}^{\dagger}\otimes(I\otimes
I-c_{1}z_{1}\sigma_{1}\otimes\sigma_{1}-c_{2}z_{2}\sigma_{2}\otimes\sigma_{2}$
$\displaystyle-c_{3}z_{3}\sigma_{3}\otimes\sigma_{3}),$ (7)

with

$\displaystyle z_{1}=2(-t_{A}y_{A2}+y_{A1}y_{A3}),$ $\displaystyle
z_{2}=2(t_{A}y_{A1}+y_{A2}y_{A3}),$ $\displaystyle
z_{3}=t_{A}^{2}-y_{A1}^{2}-y_{A2}^{2}+y_{A3}^{2}.$

Thus, we have $Tr_{C}(\rho_{0})=\frac{1}{2}V_{A}\Pi_{0}V_{A}^{\dagger}\otimes
I$ and $Tr_{C}(\rho_{1})=\frac{1}{2}V_{A}\Pi_{1}V_{A}^{\dagger}\otimes I$. The
average entropy of the subsystem $B$ after measure $\Pi^{A}$ is given by

$S_{B|\Pi^{A}}(\rho)=\frac{1}{2}\times 1+\frac{1}{2}\times 1=1.$ (8)

To evaluate $S_{C|\Pi^{AB}}(\rho)$, one needs to measure the subsystem $B$
under the conditions of the outcomes on measuring $A$. Let

$\\{B_{k}^{j}=V_{B^{j}}\Pi_{k}V_{B^{j}}^{\dagger}:k=0,1\\},~{}j=0,1,$

be the von Neumann measurement on the subsystem $B$ when the outcome of the
measurement on $A$ is $j$ ($j=0,1$), where
$V_{B^{j}}=t_{B^{j}}I+i\overrightarrow{y_{B^{j}}}\cdot\overrightarrow{\sigma}$
is the unitary operator with $t_{B^{j}}\in\mathbb{R}$,
$\overrightarrow{y_{B^{j}}}=(y_{B^{j}1},y_{B^{j}2},y_{B^{j}3})\in\mathbb{R}^{3}$
and $t_{B^{j}}^{2}+y_{B^{j}1}^{2}+y_{B^{j}2}^{2}+y_{B^{j}3}^{2}=1$.

If the measurement outcome on system $A$ is $0$,the state after the
measurement will be reduced to $\rho_{0}$ given in (6). Notice that the
subsystems $B$ and $C$ in (6) is still in a Bell-diagonal state. After
performing the measurement $\\{B_{k}^{0}:k=0,1\\}$, the state reduces to

$\displaystyle\rho_{00}$ $\displaystyle=$
$\displaystyle\frac{1}{2}V_{A}\Pi_{0}V_{A}^{\dagger}\otimes
V_{B^{0}}\Pi_{0}V_{B^{0}}^{\dagger}\otimes(I+c_{1}z_{1}l_{1}\sigma_{1}$
$\displaystyle+c_{2}z_{2}l_{2}\sigma_{2}+c_{3}z_{3}l_{3}\sigma_{3}),$
$\displaystyle\rho_{01}$ $\displaystyle=$
$\displaystyle\frac{1}{2}V_{A}\Pi_{0}V_{A}^{\dagger}\otimes
V_{B^{0}}\Pi_{1}V_{B^{0}}^{\dagger}\otimes(I-c_{1}z_{1}l_{1}\sigma_{1}$
$\displaystyle-c_{2}z_{2}l_{2}\sigma_{2}-c_{3}z_{3}l_{3}\sigma_{3}),$

with the probability $p_{00}=p_{01}=\frac{1}{4}$,where

$\displaystyle l_{1}=2(-t_{B^{0}}y_{B^{0}2}+y_{B^{0}1}y_{B^{0}3}),$
$\displaystyle l_{2}=2(t_{B^{0}}y_{B^{0}1}+y_{B^{0}2}y_{B^{0}3}),$
$\displaystyle
l_{3}=t_{B^{0}}^{2}-y_{B^{0}1}^{2}-y_{B^{0}2}^{2}+y_{B^{0}3}^{2}.$

If the measurement outcome on system $A$ is $1$, performing the measurement
$\\{B_{k}^{1}:k=0,1\\}$ on the subsystem $B$ of the state $\rho_{1}$, we
obtain

$\displaystyle\rho_{10}$ $\displaystyle=$
$\displaystyle\frac{1}{2}V_{A}\Pi_{1}V_{A}^{\dagger}\otimes
V_{B^{1}}\Pi_{0}V_{B^{1}}^{\dagger}\otimes(I-c_{1}z_{1}m_{1}\sigma_{1}$
$\displaystyle-c_{2}z_{2}m_{2}\sigma_{2}-c_{3}z_{3}m_{3}\sigma_{3}),$
$\displaystyle\rho_{11}$ $\displaystyle=$
$\displaystyle\frac{1}{2}V_{A}\Pi_{1}V_{A}^{\dagger}\otimes
V_{B^{1}}\Pi_{1}V_{B^{1}}^{\dagger}\otimes(I+c_{1}z_{1}m_{1}\sigma_{1}$
$\displaystyle+c_{2}z_{2}m_{2}\sigma_{2}+c_{3}z_{3}m_{3}\sigma_{3}),$

with the probability $p_{10}=p_{11}=\frac{1}{4}$, where

$\displaystyle m_{1}=2(-t_{B^{1}}y_{B^{1}2}+y_{B^{1}1}y_{B^{1}3}),$
$\displaystyle m_{2}=2(t_{B^{1}}y_{B^{1}1}+y_{B^{1}2}y_{B^{1}3}),$
$\displaystyle
m_{3}=t_{B^{1}}^{2}-y_{B^{1}1}^{2}-y_{B^{1}2}^{2}+y_{B^{1}3}^{2}.$

The state $\rho_{\Pi^{AB}}$ is given by
$\rho_{\Pi^{AB}}=p_{00}\rho_{00}+p_{01}\rho_{01}+p_{10}\rho_{10}+p_{11}\rho_{11}$.
Set
$\alpha=\sqrt{c_{1}^{2}z_{1}^{2}l_{1}^{2}+c_{2}^{2}z_{2}^{2}l_{2}^{2}+c_{3}^{2}z_{3}^{2}l_{3}^{2}}$
and
$\beta=\sqrt{c_{1}^{2}z_{1}^{2}m_{1}^{2}+c_{2}^{2}z_{2}^{2}m_{2}^{2}+c_{3}^{2}z_{3}^{2}m_{3}^{2}}$.
Then

$\displaystyle S_{C|\Pi^{AB}}(\rho)=$
$\displaystyle-\frac{1+\alpha}{4}\log_{2}(1+\alpha)-\frac{1-\alpha}{4}\log_{2}(1-\alpha)$
$\displaystyle-\frac{1+\beta}{4}\log_{2}(1+\beta)-\frac{1-\beta}{4}\log_{2}(1-\beta)$
$\displaystyle+1.$

It can be directly verified that $z_{1}^{2}+z_{2}^{2}+z_{3}^{2}=1$,
$l_{1}^{2}+l_{2}^{2}+l_{3}^{2}=1$, $m_{1}^{2}+m_{2}^{2}+m_{3}^{2}=1$. Denote

$c:=max\\{|c_{1}|,|c_{2}|,|c_{3}|\\}.$ (9)

Then

$\alpha\leq\sqrt{|c^{2}|(|z_{1}|^{2}|l_{1}|^{2}+|z_{2}|^{2}|l_{2}|^{2}+|z_{3}|^{2}|l_{3}|^{2})}=c,$
(10)

and

$\beta\leq\sqrt{|c^{2}|(|z_{1}|^{2}|m_{1}|^{2}+|z_{2}|^{2}|m_{2}|^{2}+|z_{3}|^{2}|m_{3}|^{2})}=c.$
(11)

The equality holds in (10) for the following cases: (1) if $c=|c_{1}|$, then
$|z_{1}|=|l_{1}|=1,z_{2}=z_{3}=l_{2}=l_{3}=0$. For instance,
$|t_{A}|=|y_{A2}|=|t_{B^{0}}|=|y_{B^{0}2}|=\frac{1}{\surd 2}$ and
$y_{A1}=y_{A3}=y_{B^{0}1}=y_{B^{0}3}=0;$ (2) if $c=|c_{2}|$, then
$|z_{2}|=|l_{2}|=1,z_{1}=z_{3}=l_{1}=l_{3}=0$. For example,
$|t_{A}|=|y_{A1}|=|t_{B^{0}}|=|y_{B^{0}1}|=\frac{1}{\surd 2}$ and
$y_{A2}=y_{A3}=y_{B^{0}2}=y_{B^{0}3}=0;$ (3) if $c=|c_{3}|,$ then
$|z_{3}|=|l_{3}|=1,z_{1}=z_{2}=l_{1}=l_{2}=0,$ e.g.,
$y_{A1}=y_{A2}=y_{B^{0}1}=y_{B^{0}2}=0.$ Similarly, one can prove that the
equality holds in (11) too for the above cases. Therefore, we obtain

$\displaystyle min(S_{C|\Pi^{AB}}(\rho))$ $\displaystyle=$
$\displaystyle-\frac{1+c}{2}\log_{2}(1+c)-\frac{1-c}{2}\log_{2}(1-c)+1.$ (12)

From (II), (8) and (II), we get the quantum discord

$\displaystyle D_{A;B;C}(\rho)$
$\displaystyle=\min_{\Pi^{AB}}\Big{[}-S_{BC|A}(\rho)+S_{B|\Pi^{A}}(\rho)+S_{C|\Pi^{AB}}(\rho)\Big{]}$
$\displaystyle=\frac{1+\xi}{2}\log_{2}(1+\xi)+\frac{1-\xi}{2}\log_{2}(1-\xi)$
$\displaystyle\quad-\frac{1+c}{2}\log_{2}(1+c)-\frac{1-c}{2}\log_{2}(1-c).$
(13)

We now consider the family of four-qubit case,

$\rho=\frac{1}{16}(I+\sum\limits_{j=1}^{3}c_{j}\sigma_{j}\otimes\sigma_{j}\otimes\sigma_{j}\otimes\sigma_{j})$
(14)

in systems $A_{1}$, $A_{2}$, $A_{3}$ and $A_{4}$. The four-qubit quantum
discord is given by

$\displaystyle D_{A_{1};A_{2};A_{3};A_{4}}(\rho)=$
$\displaystyle\min_{\Pi^{A_{1}A_{2}A_{3}}}\Big{[}-S_{A_{2}A_{3}A_{4}|A_{1}}(\rho)$
$\displaystyle+S_{A_{2}|\Pi^{A_{1}}}(\rho)+S_{A_{3}|\Pi^{A_{1}A_{2}}}(\rho)$
$\displaystyle+S_{A_{4}|\Pi^{A_{1}A_{2}A_{3}}}(\rho)\Big{]}.$ (15)

For (14) we have $\rho_{A_{1}}=Tr_{A_{2}A_{3}A_{4}}(\rho)=\frac{I}{2}$ and the
entropy of the subsystem $A_{1}$ is $S(\rho_{A_{1}})=1.$ It can be directly
verified that

$\displaystyle S(\rho)=$
$\displaystyle-\frac{1}{4}[(1+c_{1}-c_{2}-c_{3})\log_{2}(1+c_{1}-c_{2}-c_{3})$
$\displaystyle+(1-c_{1}+c_{2}-c_{3})\log_{2}(1-c_{1}+c_{2}-c_{3})$
$\displaystyle+(1-c_{1}-c_{2}+c_{3})\log_{2}(1-c_{1}-c_{2}+c_{3})$
$\displaystyle+(1+c_{1}+c_{2}+c_{3})\log_{2}(1+c_{1}+c_{2}+c_{3})]+4.$

Therefore,

$\displaystyle-S_{A_{2}A_{3}A_{4}|A_{1}}(\rho)$ $\displaystyle=$
$\displaystyle\frac{1}{4}[(1+c_{1}-c_{2}-c_{3})\log_{2}(1+c_{1}-c_{2}-c_{3})$
$\displaystyle+(1-c_{1}+c_{2}-c_{3})\log_{2}(1-c_{1}+c_{2}-c_{3})$
$\displaystyle+(1-c_{1}-c_{2}+c_{3})\log_{2}(1-c_{1}-c_{2}+c_{3})$
$\displaystyle+(1+c_{1}+c_{2}+c_{3})\log_{2}(1+c_{1}+c_{2}+c_{3})]-3.$ (16)

The von Neumann measurement on the subsystem $A_{1}$ is given by
$\\{A_{1k}=V_{A_{1}}\Pi_{k}V_{A_{1}}^{\dagger}:~{}k=0,1\\}$, where
$V_{A_{1}}=t_{A_{1}}I+i\overrightarrow{y_{A_{1}}}\cdot\overrightarrow{\sigma}$,
with $t_{A_{1}}\in\mathbb{R}$,
$\overrightarrow{y_{A_{1}}}=(y_{{A_{1}}1},y_{{A_{1}}2},y_{{A_{1}}3})\in\mathbb{R}^{3}$,
and $t_{A_{1}}^{2}+y_{{A_{1}}1}^{2}+y_{{A_{1}}2}^{2}+y_{{A_{1}}3}^{2}=1.$

The state $\rho_{\Pi^{A_{1}}}$ is given by
$\rho_{\Pi^{A_{1}}}=p_{0}\rho_{0}+p_{1}\rho_{1}$, where
$p_{0}=p_{1}=\frac{1}{2}$, and

$\displaystyle\rho_{0}=$
$\displaystyle\frac{1}{8}V_{A_{1}}\Pi_{0}V_{A_{1}}^{\dagger}\otimes(I\otimes
I\otimes I+c_{1}d_{1}\sigma_{1}\otimes\sigma_{1}\otimes\sigma_{1}$
$\displaystyle+c_{2}d_{2}\sigma_{2}\otimes\sigma_{2}\otimes\sigma_{2}+c_{3}d_{3}\sigma_{3}\otimes\sigma_{3}\otimes\sigma_{3}),$
$\displaystyle\rho_{1}=$
$\displaystyle\frac{1}{8}V_{A_{1}}\Pi_{1}V_{A_{1}}^{\dagger}\otimes(I\otimes
I\otimes I-c_{1}d_{1}\sigma_{1}\otimes\sigma_{1}\otimes\sigma_{1}$
$\displaystyle-
c_{2}d_{2}\sigma_{2}\otimes\sigma_{2}\otimes\sigma_{2}-c_{3}d_{3}\sigma_{3}\otimes\sigma_{3}\otimes\sigma_{3}),$

where

$\displaystyle d_{1}=2(-t_{A_{1}}y_{{A_{1}}2}+y_{{A_{1}}1}y_{{A_{1}}3}),$
$\displaystyle d_{2}=2(t_{A_{1}}y_{{A_{1}}1}+y_{{A_{1}}2}y_{{A_{1}}3}),$
$\displaystyle
d_{3}=t_{A_{1}}^{2}-y_{{A_{1}}1}^{2}-y_{{A_{1}}2}^{2}+y_{{A_{1}}3}^{2}.$

Thus, we have
$Tr_{A_{3}A_{4}}(\rho_{0})=\frac{1}{2}V_{A_{1}}\Pi_{0}V_{A_{1}}^{\dagger}\otimes
I$ and
$Tr_{A_{3}A_{4}}(\rho_{1})=\frac{1}{2}V_{A_{1}}\Pi_{1}V_{A_{1}}^{\dagger}\otimes
I$. The average entropy of the subsystem $A_{2}$ after the measurement
$\Pi^{A_{1}}$ is given by $S_{{A_{2}}|\Pi^{A_{1}}}(\rho)=1$.

To evaluate $S_{A_{3}|\Pi^{A_{1}A_{2}}}(\rho)$ and
$S_{A_{4}|\Pi^{A_{1}A_{2}A_{3}}}(\rho)$, we need to measure the subsystem
$A_{2}$ based on the measurement outcomes on $A_{1}$. We obtain

$\displaystyle\rho_{00}=$
$\displaystyle\frac{1}{4}V_{A_{1}}\Pi_{0}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{0}}\Pi_{0}V_{A_{2}^{0}}^{\dagger}\otimes(I\otimes
I+c_{1}d_{1}e_{1}\sigma_{1}\otimes\sigma_{1}$
$\displaystyle+c_{2}d_{2}e_{2}\sigma_{2}\otimes\sigma_{2}+c_{3}d_{3}e_{3}\sigma_{3}\otimes\sigma_{3}),$
$\displaystyle\rho_{01}=$
$\displaystyle\frac{1}{4}V_{A_{1}}\Pi_{0}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{0}}\Pi_{1}V_{A_{2}^{0}}^{\dagger}\otimes(I\otimes
I-c_{1}d_{1}e_{1}\sigma_{1}\otimes\sigma_{1}$ $\displaystyle-
c_{2}d_{2}e_{2}\sigma_{2}\otimes\sigma_{2}-c_{3}d_{3}e_{3}\sigma_{3}\otimes\sigma_{3}),$
$\displaystyle\rho_{10}=$
$\displaystyle\frac{1}{4}V_{A_{1}}\Pi_{1}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{1}}\Pi_{0}V_{A_{2}^{1}}^{\dagger}\otimes(I\otimes
I-c_{1}d_{1}f_{1}\sigma_{1}\otimes\sigma_{1}$ $\displaystyle-
c_{2}d_{2}f_{2}\sigma_{2}\otimes\sigma_{2}-c_{3}d_{3}f_{3}\sigma_{3}\otimes\sigma_{3}),$
$\displaystyle\rho_{11}=$
$\displaystyle\frac{1}{4}V_{A_{1}}\Pi_{1}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{1}}\Pi_{1}V_{A_{2}^{1}}^{\dagger}\otimes(I\otimes
I+c_{1}d_{1}f_{1}\sigma_{1}\otimes\sigma_{1}$
$\displaystyle+c_{2}d_{2}f_{2}\sigma_{2}\otimes\sigma_{2}+c_{3}d_{3}f_{3}\sigma_{3}\otimes\sigma_{3}),$

where the $k$ in the unitary $\\{V_{A_{2}^{k}}:k=0,1\\}$ is the outcome of the
measurement of $A_{1}$ , and $V_{A_{2}^{k}}$ can be written as
$V_{A_{2}^{k}}=t_{A_{2}^{k}}I+i\overrightarrow{y_{A_{2}^{k}}}\cdot\overrightarrow{\sigma}$,
with $t_{A_{2}^{k}}\in\mathbb{R}$,
$\overrightarrow{y_{A_{2}^{k}}}=(y_{{A_{2}^{k}}1},y_{{A_{2}^{k}}2},y_{{A_{2}^{k}}3})\in\mathbb{R}^{3}$,
and
$t_{A_{2}^{k}}^{2}+y_{{A_{2}^{k}}1}^{2}+y_{{A_{2}^{k}}2}^{2}+y_{{A_{2}^{k}}3}^{2}=1,$

$\displaystyle
e_{1}=2(-t_{A_{2}^{0}}y_{A_{2}^{0}2}+y_{A_{2}^{0}1}y_{A_{2}^{0}3}),$
$\displaystyle
e_{2}=2(t_{A_{2}^{0}}y_{A_{2}^{0}1}+y_{A_{2}^{0}2}y_{A_{2}^{0}3}),$
$\displaystyle
e_{3}=t_{A_{2}^{0}}^{2}-y_{A_{2}^{0}1}^{2}-y_{A_{2}^{0}2}^{2}+y_{A_{2}^{0}3}^{2},$
$\displaystyle
f_{1}:=2(-t_{A_{2}^{1}}y_{A_{2}^{1}2}+y_{A_{2}^{1}1}y_{A_{2}^{1}3}),$
$\displaystyle
f_{2}:=2(t_{A_{2}^{1}}y_{A_{2}^{1}1}+y_{A_{2}^{1}2}y_{A_{2}^{1}3}),$
$\displaystyle
f_{3}:=t_{A_{2}^{1}}^{2}-y_{A_{2}^{1}1}^{2}-y_{A_{2}^{1}2}^{2}+y_{A_{2}^{1}3}^{2}.$

The state $\rho_{\Pi^{A_{1}A_{2}}}$ is given by
$\rho_{\Pi^{A_{1}A_{2}}}=p_{00}\rho_{00}+p_{01}\rho_{01}+p_{10}\rho_{10}+p_{11}\rho_{11}$.
Thus, we have
$Tr_{A_{4}}(\rho_{00})=\frac{1}{2}V_{A_{1}}\Pi_{0}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{0}}\Pi_{0}V_{A_{2}^{0}}^{\dagger}\otimes I,$
$Tr_{A_{4}}(\rho_{01})=\frac{1}{2}V_{A_{1}}\Pi_{0}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{0}}\Pi_{1}V_{A_{2}^{0}}^{\dagger}\otimes I,$
$Tr_{A_{4}}(\rho_{10})=\frac{1}{2}V_{A_{1}}\Pi_{1}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{1}}\Pi_{0}V_{A_{2}^{1}}^{\dagger}\otimes I,$
$Tr_{A_{4}}(\rho_{11})=\frac{1}{2}V_{A_{1}}\Pi_{1}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{1}}\Pi_{1}V_{A_{2}^{1}}^{\dagger}\otimes I$. The average entropy of
the sub system $A_{3}$ after the measurement $\Pi^{A_{1}A_{2}}$ is given by
$S_{{A_{3}}|\Pi^{A_{1}A_{2}}}(\rho)=4\times(\frac{1}{4}\times 1)=1$.

To evaluate $S_{A_{4}|\Pi^{A_{1}A_{2}A_{3}}}(\rho)$, one needs to continue to
measure the subsystem $A_{3}$ based on the measurement outcomes on $A_{1}$ and
$A_{2}$. We obtain

$\displaystyle\rho_{000}=$
$\displaystyle\frac{1}{2}V_{A_{1}}\Pi_{0}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{0}}\Pi_{0}V_{A_{2}^{0}}^{\dagger}\otimes
V_{A_{3}^{00}}\Pi_{0}V_{A_{3}^{00}}^{\dagger}\otimes(I$
$\displaystyle+c_{1}d_{1}e_{1}g_{1}\sigma_{1}+c_{2}d_{2}e_{2}g_{2}\sigma_{2}+c_{3}d_{3}e_{3}g_{3}\sigma_{3}),$
$\displaystyle\rho_{001}=$
$\displaystyle\frac{1}{2}V_{A_{1}}\Pi_{0}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{0}}\Pi_{0}V_{A_{2}^{0}}^{\dagger}\otimes
V_{A_{3}^{00}}\Pi_{1}V_{A_{3}^{00}}^{\dagger}\otimes(I$ $\displaystyle-
c_{1}d_{1}e_{1}g_{1}\sigma_{1}-c_{2}d_{2}e_{2}g_{2}\sigma_{2}-c_{3}d_{3}e_{3}g_{3}\sigma_{3}),$
$\displaystyle\rho_{010}=$
$\displaystyle\frac{1}{2}V_{A_{1}}\Pi_{0}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{0}}\Pi_{1}V_{A_{2}^{0}}^{\dagger}\otimes
V_{A_{3}^{01}}\Pi_{0}V_{A_{3}^{01}}^{\dagger}\otimes(I$ $\displaystyle-
c_{1}d_{1}e_{1}h_{1}\sigma_{1}-c_{2}d_{2}e_{2}h_{2}\sigma_{2}-c_{3}d_{3}e_{3}h_{3}\sigma_{3}),$
$\displaystyle\rho_{011}=$
$\displaystyle\frac{1}{2}V_{A_{1}}\Pi_{0}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{0}}\Pi_{1}V_{A_{2}^{0}}^{\dagger}\otimes
V_{A_{3}^{01}}\Pi_{1}V_{A_{3}^{01}}^{\dagger}\otimes(I$
$\displaystyle+c_{1}d_{1}e_{1}h_{1}\sigma_{1}+c_{2}d_{2}e_{2}h_{2}\sigma_{2}+c_{3}d_{3}e_{3}h_{3}\sigma_{3}),$
$\displaystyle\rho_{100}=$
$\displaystyle\frac{1}{2}V_{A_{1}}\Pi_{1}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{1}}\Pi_{0}V_{A_{2}^{1}}^{\dagger}\otimes
V_{A_{3}^{10}}\Pi_{0}V_{A_{3}^{10}}^{\dagger}\otimes(I$ $\displaystyle-
c_{1}d_{1}f_{1}n_{1}\sigma_{1}-c_{2}d_{2}f_{2}n_{2}\sigma_{2}-c_{3}d_{3}f_{3}n_{3}\sigma_{3}),$
$\displaystyle\rho_{101}=$
$\displaystyle\frac{1}{2}V_{A_{1}}\Pi_{1}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{1}}\Pi_{0}V_{A_{2}^{1}}^{\dagger}\otimes
V_{A_{3}^{10}}\Pi_{1}V_{A_{3}^{10}}^{\dagger}\otimes(I$
$\displaystyle+c_{1}d_{1}f_{1}n_{1}\sigma_{1}+c_{2}d_{2}f_{2}n_{2}\sigma_{2}+c_{3}d_{3}f_{3}n_{3}\sigma_{3}),$
$\displaystyle\rho_{110}=$
$\displaystyle\frac{1}{2}V_{A_{1}}\Pi_{1}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{1}}\Pi_{1}V_{A_{2}^{1}}^{\dagger}\otimes
V_{A_{3}^{11}}\Pi_{0}V_{A_{3}^{11}}^{\dagger}\otimes(I$
$\displaystyle+c_{1}d_{1}f_{1}r_{1}\sigma_{1}+c_{2}d_{2}f_{2}r_{2}\sigma_{2}+c_{3}d_{3}f_{3}r_{3}\sigma_{3}),$
$\displaystyle\rho_{111}=$
$\displaystyle\frac{1}{2}V_{A_{1}}\Pi_{1}V_{A_{1}}^{\dagger}\otimes
V_{A_{2}^{1}}\Pi_{1}V_{A_{2}^{1}}^{\dagger}\otimes
V_{A_{3}^{11}}\Pi_{1}V_{A_{3}^{11}}^{\dagger}\otimes(I$ $\displaystyle-
c_{1}d_{1}f_{1}r_{1}\sigma_{1}-c_{2}d_{2}f_{2}r_{2}\sigma_{2}-c_{3}d_{3}f_{3}r_{3}\sigma_{3}),$

where the $k$ in the unitary $\\{V_{A_{3}^{ku}}:k=0,1;u=0,1\\}$ is the outcome
of the measurement of $A_{1}$ , and $u$ is the outcome of the measurement of
$A_{2}$. Denote

$\displaystyle\mu_{1}=\sqrt{c_{1}^{2}d_{1}^{2}e_{1}^{2}g_{1}^{2}+c_{2}^{2}d_{2}^{2}e_{2}^{2}g_{2}^{2}+c_{3}^{2}d_{3}^{2}e_{3}^{2}g_{3}^{2}},$
$\displaystyle\mu_{2}=\sqrt{c_{1}^{2}d_{1}^{2}e_{1}^{2}h_{1}^{2}+c_{2}^{2}d_{2}^{2}e_{2}^{2}h_{2}^{2}+c_{3}^{2}d_{3}^{2}e_{3}^{2}h_{3}^{2}},$
$\displaystyle\mu_{3}=\sqrt{c_{1}^{2}d_{1}^{2}f_{1}^{2}n_{1}^{2}+c_{2}^{2}d_{2}^{2}f_{2}^{2}n_{2}^{2}+c_{3}^{2}d_{3}^{2}f_{3}^{2}n_{3}^{2}},$
$\displaystyle\mu_{4}=\sqrt{c_{1}^{2}d_{1}^{2}f_{1}^{2}r_{1}^{2}+c_{2}^{2}d_{2}^{2}f_{2}^{2}r_{2}^{2}+c_{3}^{2}d_{3}^{2}f_{3}^{2}r_{3}^{2}}.$

We have

$\displaystyle S_{A_{4}|\Pi^{A_{1}A_{2}A_{3}}}(\rho)$ $\displaystyle=$
$\displaystyle-\frac{1+\mu_{1}}{8}\log_{2}(1+\mu_{1})-\frac{1-\mu_{1}}{8}\log_{2}(1-\mu_{1})$
$\displaystyle-\frac{1+\mu_{2}}{8}\log_{2}(1+\mu_{2})-\frac{1-\mu_{2}}{8}\log_{2}(1-\mu_{2})$
$\displaystyle-\frac{1+\mu_{3}}{8}\log_{2}(1+\mu_{3})-\frac{1-\mu_{3}}{8}\log_{2}(1-\mu_{3})$
$\displaystyle-\frac{1+\mu_{4}}{8}\log_{2}(1+\mu_{4})-\frac{1-\mu_{4}}{8}\log_{2}(1-\mu_{4})+1.$

It can be directly verified that $d_{1}^{2}+d_{2}^{2}+d_{3}^{2}=1,$
$e_{1}^{2}+e_{2}^{2}+e_{3}^{2}=1,$ $f_{1}^{2}+f_{2}^{2}+f_{3}^{2}=1,$
$g_{1}^{2}+g_{2}^{2}+g_{3}^{2}=1,$ $h_{1}^{2}+h_{2}^{2}+h_{3}^{2}=1,$
$n_{1}^{2}+n_{2}^{2}+n_{3}^{2}=1$ and $r_{1}^{2}+r_{2}^{2}+r_{3}^{2}=1.$ Since
$\mu_{1}\leq c$, $\mu_{2}\leq c$, $\mu_{3}\leq c$ and $\mu_{4}\leq c$, we
obtain

$\displaystyle min(S_{A_{4}|\Pi^{A_{1}A_{2}A_{3}}}(\rho))$ $\displaystyle=$
$\displaystyle-\frac{1+c}{2}\log_{2}(1+c)-\frac{1-c}{2}\log_{2}(1-c)+1.$ (17)

By the definition of the four-qubit quantum discord (15), we get

$\displaystyle D_{A_{1};A_{2};A_{3};A_{4}}(\rho)$ $\displaystyle=$
$\displaystyle\frac{1}{4}[(1+c_{1}-c_{2}-c_{3})\log_{2}(1+c_{1}-c_{2}-c_{3})$
$\displaystyle+(1-c_{1}+c_{2}-c_{3})\log_{2}(1-c_{1}+c_{2}-c_{3})$
$\displaystyle+(1-c_{1}-c_{2}+c_{3})\log_{2}(1-c_{1}-c_{2}+c_{3})$
$\displaystyle+(1+c_{1}+c_{2}+c_{3})\log_{2}(1+c_{1}+c_{2}+c_{3})]$
$\displaystyle-\frac{1+c}{2}\log_{2}(1+c)-\frac{1-c}{2}\log_{2}(1-c).$ (18)

From the results of three-qubit and four-quibt states, we can prove the
following conclusion for general $N$-qubit case.

###### Theorem 1

For the family of $N$-qubit states (2), we have the quantum discord:

(1) if $N=2v+1$, $v\in\mathbf{N^{+}}$,

$\displaystyle D_{A_{1};A_{2};\dots;A_{2v+1}}(\rho)$ $\displaystyle=$
$\displaystyle\frac{1+\xi}{2}\log_{2}(1+\xi)+\frac{1-\xi}{2}\log_{2}(1-\xi)$
$\displaystyle-\frac{1+c}{2}\log_{2}(1+c)-\frac{1-c}{2}\log_{2}(1-c),$ (19)

where $\xi=\sqrt{c_{1}^{2}+c_{2}^{2}+c_{3}^{2}},$
$c=max\\{|c_{1}|,|c_{2}|,|c_{3}|\\};$

(2) if $N=4v-2$, $v\in\mathbf{N^{+}}$,

$\displaystyle D_{A_{1};A_{2};\dots;A_{4v-2}}(\rho)$ $\displaystyle=$
$\displaystyle\frac{1}{4}[(1-c_{1}-c_{2}-c_{3})\log_{2}(1-c_{1}-c_{2}-c_{3})$
$\displaystyle+(1-c_{1}+c_{2}+c_{3})\log_{2}(1-c_{1}+c_{2}+c_{3})$
$\displaystyle+(1+c_{1}-c_{2}+c_{3})\log_{2}(1+c_{1}-c_{2}+c_{3})$
$\displaystyle+(1+c_{1}+c_{2}-c_{3})\log_{2}(1+c_{1}+c_{2}-c_{3})]$
$\displaystyle-\frac{1+c}{2}\log_{2}(1+c)-\frac{1-c}{2}\log_{2}(1-c);$ (20)

(3) if $N=4v$, $v\in\mathbf{N^{+}}$,

$\displaystyle D_{A_{1};A_{2};\dots;A_{4v}}(\rho)$ $\displaystyle=$
$\displaystyle\frac{1}{4}[(1+c_{1}-c_{2}-c_{3})\log_{2}(1+c_{1}-c_{2}-c_{3})$
$\displaystyle+(1-c_{1}+c_{2}-c_{3})\log_{2}(1-c_{1}+c_{2}-c_{3})$
$\displaystyle+(1-c_{1}-c_{2}+c_{3})\log_{2}(1-c_{1}-c_{2}+c_{3})$
$\displaystyle+(1+c_{1}+c_{2}+c_{3})\log_{2}(1+c_{1}+c_{2}+c_{3})]$
$\displaystyle-\frac{1+c}{2}\log_{2}(1+c)-\frac{1-c}{2}\log_{2}(1-c).$ (21)

###### Proof 1

If $N=2v+1$, we have $Tr_{A_{2}A_{3}\dots A_{2v+1}}(\rho)=\frac{I}{2}$, and
$S_{A_{1}}(\rho)=1.$ Let $\lambda$ be the eigenvalues of $\rho$. From the
characteristic equation $det|\rho_{2v+1}-\lambda I|=0$, we get

$\displaystyle[\frac{(1-2^{2v+1}\lambda)^{2}-c_{1}^{2}-c_{2}^{2}-c_{3}^{2}}{2^{(4v+2)}}]^{2^{2v}}=0.$

The $2^{2v}$ eigenvalues are given by
$\frac{1}{2^{2v+1}}(1-\sqrt{c_{1}^{2}+c_{2}^{2}+c_{3}^{2}})$ and
$\frac{1}{2^{2v+1}}(1+\sqrt{c_{1}^{2}+c_{2}^{2}+c_{3}^{2}})$, respectively.
One can verify that

$\displaystyle-S_{A_{2};A_{3};\dots;A_{2v+1}|A_{1}}(\rho)$ $\displaystyle=$
$\displaystyle\frac{1+\xi}{2}\log_{2}(1+\xi)+\frac{1-\xi}{2}\log_{2}(1-\xi)-2v.$
(22)

We obtain $S_{A_{k}|\Pi^{A_{1}A_{2}\dots A_{k-1}}}(\rho)=1$, where
$k=2,\cdots,2v$, namely,

$\displaystyle
S_{A_{2}|\Pi^{A_{1}}}(\rho)=S_{A_{3}|\Pi^{A_{1}A_{2}}}(\rho)=\dots$
$\displaystyle=$ $\displaystyle S_{A_{2v}|\Pi^{A_{1}A_{2}\dots
A_{2v-1}}}(\rho)=1.$ (23)

Hence,

$\displaystyle min(S_{A_{2v+1}|\Pi^{A_{1}A_{2}\dots A_{2v}}}(\rho))$
$\displaystyle=$
$\displaystyle-\frac{1+c}{2}\log_{2}(1+c)-\frac{1-c}{2}\log_{2}(1-c)+1.$ (24)

By the definition (1), we obtain (1).

If $N=4v-2$, we have $Tr_{A_{2}A_{3}\dots A_{4v-2}}(\rho)=\frac{I}{2}$ and
$S_{A_{1}}(\rho)=1.$ $\rho$ has $2^{4v-4}$ eigenvalues given by
$\frac{1}{2^{4v-2}}(1-c_{1}-c_{2}-c_{3})$,
$\frac{1}{2^{4v-2}}(1-c_{1}+c_{2}+c_{3})$,
$\frac{1}{2^{4v-2}}(1+c_{1}-c_{2}+c_{3})$ and
$\frac{1}{2^{4v-2}}(1+c_{1}+c_{2}-c_{3})$, respectively. Then we obtain
$-S_{A_{2};A_{3};\dots;A_{4}v-2}(\rho)=\frac{1}{4}[(1-c_{1}-c_{2}-c_{3})\log_{2}(1-c_{1}-c_{2}-c_{3})+(1-c_{1}+c_{2}+c_{3})\log_{2}(1-c_{1}+c_{2}+c_{3})+(1+c_{1}-c_{2}+c_{3})\log_{2}(1+c_{1}-c_{2}+c_{3})+(1+c_{1}+c_{2}-c_{3})\log_{2}(1+c_{1}+c_{2}-c_{3})]-4v+3.$
The entropy after the measurement is the same as (23) and (24). Therefore, we
obtain (1). (1) is similarly proved.

From the Theorem one has that for five-qubit states,

$\displaystyle D_{A_{1};A_{2};A_{3};A_{4};A_{5}}(\rho)$ $\displaystyle=$
$\displaystyle\frac{1+\xi}{2}\log_{2}(1+\xi)+\frac{1-\xi}{2}\log_{2}(1-\xi)$
$\displaystyle-\frac{1+c}{2}\log_{2}(1+c)-\frac{1-c}{2}\log_{2}(1-c),$ (25)

and for six-qubit states

$\displaystyle D_{A_{1};A_{2};A_{3};A_{4};A_{5};A_{6}}(\rho)$ $\displaystyle=$
$\displaystyle\frac{1}{4}[(1-c_{1}-c_{2}-c_{3})\log_{2}(1-c_{1}-c_{2}-c_{3})$
$\displaystyle+(1-c_{1}+c_{2}+c_{3})\log_{2}(1-c_{1}+c_{2}+c_{3})$
$\displaystyle+(1+c_{1}-c_{2}+c_{3})\log_{2}(1+c_{1}-c_{2}+c_{3})$
$\displaystyle+(1+c_{1}+c_{2}-c_{3})\log_{2}(1+c_{1}+c_{2}-c_{3})]$
$\displaystyle-\frac{1+c}{2}\log_{2}(1+c)-\frac{1-c}{2}\log_{2}(1-c).$ (26)

Interestingly, the result (25) is equivalent to the three-qubit quantum
discord (13), while (26) is equivalent to the two-qubit quantum discord given
by Luo luo1 .

Fig. 1 shows the level surfaces of discord for $D(\rho)=0.03,0.15$ and $0.55$.
The three figures $(F_{21})$, $(F_{22})$ and $(F_{23})$ in the second row of
Fig. 1 are for $4v-2$-qubit states, which are in consist with the ones given
in lang  for two-qubit states. For small discord, $D(\rho)=0.03$ and $0.15$,
the level surfaces are centrally symmetric, consisting of three intersecting
”tubes” along the three coordinate axes. For larger discord value $0.55$,
these intersecting tubes expand until only a few vertices remained, where
$(F_{13})$ has level surfaces in eight corners, while $(F_{23})$ and
$(F_{33})$ have only four corners left.

Figure 1: Level surfaces of constant discord. $N=2v+1$ for figures $(F_{11})$,
$(F_{12})$ and $(F_{13})$ with $D(\rho)=0.03,~{}0.15$ and $0.55$,
respectively. $N=4v-2$ for figures $(F_{21})$, $(F_{22})$ and $(F_{23})$ with
$D(\rho)=0.03,~{}0.15$ and $0.55$, respectively. $N=4v$ for figures
$(F_{31})$, $(F_{32})$ and $(F_{33})$ with $D(\rho)=0.03,~{}0.15$ and $0.55$,
respectively.

##  III Dynamics of quantum discord under local nondissipative channels

It has been discovered that for two-qubit states, the quantum discord is
invariant under some decoherence channels in a finite time intervalMazzola ;
Li . To verify if such phenomena still exist in multi-qubit systems, we
consider that the states $\rho$ (3) and (14) under the phase flip channel,
with the Kraus operators $\Gamma_{0}^{(A_{1})}=$
diag$(\sqrt{1-p/2},\sqrt{1-p/2})\otimes I\otimes\dots\otimes I$,
$\Gamma_{1}^{(A_{1})}=$ diag$(\sqrt{p/2},-\sqrt{p/2})\otimes
I\otimes\dots\otimes I$, $\dots$, $\Gamma_{0}^{(A_{N})}=I\otimes\dots\otimes
I\otimes$ diag$(\sqrt{1-p/2},\sqrt{1-p/2})$,
$\Gamma_{1}^{(A_{N})}=I\otimes\dots\otimes I\otimes$
diag$(\sqrt{p/2},-\sqrt{p/2})$, where $N=3,4$, $p=1-\exp(-\gamma t)$, $\gamma$
is the phase damping rate.

Let $\varepsilon(\cdot)$ represent the operator of decoherence. For the three-
qubit state (3) under the phase flip channel, we have

$\displaystyle\varepsilon(\rho)=$ $\displaystyle\frac{1}{8}(I\otimes I\otimes
I+(1-p)^{3}c_{1}\sigma_{1}\otimes\sigma_{1}\otimes\sigma_{1}$
$\displaystyle+(1-p)^{3}c_{2}\sigma_{2}\otimes\sigma_{2}\otimes\sigma_{2}+c_{3}\sigma_{3}\otimes\sigma_{3}\otimes\sigma_{3}).$
(27)

From (13), we obtain

$\displaystyle D_{A_{1};A_{2};A_{3}}(\varepsilon(\rho))$ $\displaystyle=$
$\displaystyle\frac{1+\delta}{2}\log_{2}(1+\delta)+\frac{1-\delta}{2}\log_{2}(1-\delta)$
$\displaystyle-\frac{1+\theta}{2}\log_{2}(1+\theta)-\frac{1-\theta}{2}\log_{2}(1-\theta),$
(28)

where $\delta=\sqrt{(1-p)^{6}c_{1}^{2}+(1-p)^{6}c_{2}^{2}+c_{3}^{2}},$
$\theta=max\\{|(1-p)^{3}c_{1}|,|(1-p)^{3}c_{2}|,|c_{3}|\\}.$ Notice that the
derivative of $D_{A_{1};A_{2};A_{3}}(\varepsilon(\rho))$ with respect to $p$
is always less than $0$.

[Proof] From (28), in general, the derivative of
$D_{A_{1};A_{2};A_{3}}(\varepsilon(\rho))$ can be cast as

$\displaystyle
D^{\prime}_{A_{1};A_{2};A_{3}}(\varepsilon(\rho))=\frac{1}{2}[\theta^{\prime}\log_{2}(\frac{1-\theta}{1+\theta})-\delta^{\prime}\log_{2}(\frac{1-\delta}{1+\delta})].$

In particular for

$\displaystyle\delta=\sqrt{(1-p)^{6}c_{1}^{2}+(1-p)^{6}c_{2}^{2}+c_{3}^{2}}>0,$
$\displaystyle\theta=max\\{|(1-p)^{3}c_{1}|,|(1-p)^{3}c_{2}|,|c_{3}|\\}>0.$

This implies that

$\displaystyle\log_{2}(\frac{1-\theta}{1+\theta})<0\quad
and\quad\log_{2}(\frac{1-\delta}{1+\delta})<0.$

we have that

$\displaystyle\delta^{\prime}=-\frac{3(c_{1}^{2}+c_{2}^{2})}{\delta}(1-p)^{5},$

given that $0<(1-p)^{5}<1$ then $\delta^{\prime}<0.$

If $\theta=|c_{3}|$ then

$\displaystyle
D^{\prime}_{A_{1};A_{2};A_{3}}(\varepsilon(\rho))=\frac{-\delta^{\prime}}{2}\log_{2}(\frac{1-\delta}{1+\delta})<0.$

And for $\theta=|(1-p)^{3}c_{1}|$, in this case $\theta=(1-p)^{3}|c_{1}|$
because $0<1-p<1$ then

$\displaystyle\theta^{\prime}=-3(1-p)^{2}|c_{1}|<0.$

Then

$\displaystyle D^{\prime}_{A_{1};A_{2};A_{3}}(\varepsilon(\rho))=$
$\displaystyle\frac{3}{2}[\frac{(c_{1}^{2}+c_{2}^{2})}{\delta}(1-p)^{5}\log_{2}(\frac{1-\delta}{1+\delta})$
$\displaystyle-(1-p)^{2}|c_{1}|\log_{2}(\frac{1-\theta}{1+\theta})].$

Here, we assume that $D_{A_{1};A_{2};A_{3}}(\varepsilon(\rho))$ is
monotonically decreasing, i.e.,
$D^{\prime}_{A_{1};A_{2};A_{3}}(\varepsilon(\rho))<0.$ In order to show these,
it must satisfy that

$\displaystyle\frac{(c_{1}^{2}+c_{2}^{2})}{\delta}(1-p)^{5}\log_{2}(\frac{1-\delta}{1+\delta})<(1-p)^{2}|c_{1}|\log_{2}(\frac{1-\theta}{1+\theta}),$

by multiplying by the positive numbers $(1-p)$ and $\delta$

$\displaystyle(c_{1}^{2}+c_{2}^{2})(1-p)^{6}\log_{2}(\frac{1-\delta}{1+\delta})<\delta\theta\log_{2}(\frac{1-\theta}{1+\theta}).$

Therefore, $D^{\prime}_{A_{1};A_{2};A_{3}}(\varepsilon(\rho))<0,$ which means
that $D_{A_{1};A_{2};A_{3}}(\varepsilon(\rho))$ is monotonically decreasing.
$\Box$

Hence, the frozen phenomenon of quantum discord dos not exist for three-qubit
states under the phase flip channel. Since the quantum discord of three-qubit
and $(2v+1)$-qubit are the same, the odd-qubit systems do not exhibit frozen
phenomenon of quantum discord under the phase flip channel. For instance, take
$c_{1}=\frac{4}{5}$, $c_{2}=\frac{c_{1}}{2}$ and $c_{3}=\frac{1}{2}$ in the
initial state, the dashed line in Fig. 2 shows the dynamic behavior of the
quantum discord under the phase flip channel.

Figure 2: Quantum discord of the three-qubit state (dashed line) and quantum
discord of the four-qubit state (solid line) under phase flip channel for
$c_{1}=\frac{4}{5}$, $c_{2}=\frac{c_{1}}{2}$, $c_{3}=\frac{1}{2}$.

For $N=4$, the state $\rho$ under the phase flip channel is give by

$\displaystyle\varepsilon(\rho)=$ $\displaystyle\frac{1}{16}(I\otimes I\otimes
I\otimes
I+(1-p)^{4}c_{1}\sigma_{1}\otimes\sigma_{1}\otimes\sigma_{1}\otimes\sigma_{1}$
$\displaystyle+(1-p)^{4}c_{2}\sigma_{2}\otimes\sigma_{2}\otimes\sigma_{2}\otimes\sigma_{2}$
$\displaystyle+c_{3}\sigma_{3}\otimes\sigma_{3}\otimes\sigma_{3}\otimes\sigma_{3}).$
(29)

Noting that $c_{3}$ is independent on time, we consider the case that
$c_{2}=c_{1}c_{3},~{}-1\leq c_{3}\leq 1$. Then we have $|c_{2}|\leq|c_{1}|$
for any $p$. From (II) we obtain the quantum discord

$\displaystyle D_{A_{1};A_{2};A_{3};A_{4}}(\varepsilon(\rho))$
$\displaystyle=$
$\displaystyle\frac{1+c_{3}}{2}\log_{2}(1+c_{3})+\frac{1-c_{3}}{2}\log_{2}(1-c_{3})$
$\displaystyle+\frac{1+(1-p)^{4}c_{1}}{2}\log_{2}(1+(1-p)^{4}c_{1})$
$\displaystyle+\frac{1-(1-p)^{4}c_{1}}{2}\log_{2}(1-(1-p)^{4}c_{1})$
$\displaystyle-\frac{1+\sigma}{2}\log_{2}(1+\sigma)-\frac{1-\sigma}{2}\log_{2}(1-\sigma),$

where $\sigma=max\\{|(1-p)^{4}c_{1}|,|c_{3}|\\}.$

When $max\\{|(1-p)^{4}c_{1}|,|c_{3}|\\}=|(1-p)^{4}c_{1}|$, we have

$\displaystyle D_{A_{1};A_{2};A_{3};A_{4}}(\varepsilon(\rho))$
$\displaystyle=$
$\displaystyle\frac{1+c_{3}}{2}\log_{2}(1+c_{3})+\frac{1-c_{3}}{2}\log_{2}(1-c_{3}),$

$D_{A_{1};A_{2};A_{3};A_{4}}(\varepsilon(\rho))$ is constant under the
decoherence channel during the time interval. Otherwise,

$\displaystyle D_{A_{1};A_{2};A_{3};A_{4}}(\varepsilon(\rho))$
$\displaystyle=$
$\displaystyle\frac{1+(1-p)^{4}c_{1}}{2}\log_{2}(1+(1-p)^{4}c_{1})$
$\displaystyle+\frac{1-(1-p)^{4}c_{1}}{2}\log_{2}(1-(1-p)^{4}c_{1}),$

which monotonically decreases to zero.

Therefore, to calculate the quantum discord, we need to determine the
magnitude of $|(1-p)^{4}c_{1}|$ and $|c_{3}|$. If for $|c_{1}|>|c_{3}|$ there
exist $0\leq p_{0}\leq 1$ such that
$max\\{|(1-p)^{4}c_{1}|,|c_{3}|\\}=|(1-p)^{4}c_{1}|$ for $0\leq p\leq p_{0},$
and $max\\{|(1-p)^{4}c_{1}|,|c_{3}|\\}=|c_{3}|$ for $p_{0}\leq p\leq 1$, then
$D_{A_{1};A_{2};A_{3};A_{4}}(\varepsilon(\rho))$ remains unchanged first, and
then monotonicity goes down to zero. As an example, set $c_{1}=\frac{4}{5}$,
$c_{2}=\frac{c_{1}}{2}$ and $c_{3}=\frac{1}{2}$. The solid line in Fig. 2
shows the dynamic behavior of quantum discord under the phase flip channel. A
sudden transition of quantum discord happens at $p=0.11086.$ The frozen
phenomenon of quantum discord exists for the four-qubit states under the phase
flip channel, while for the case of three-qubit states, such phenomenon does
not exist.

In Li , it has been shown that the frozen phenomenon of quantum discord also
exists when the phase noise acts on two-qubit states. Since the quantum
discord of two-qubit and $(4v-2)$-qubit are the same, and the four-qubit
quantum discord is equal to that of $(4v)$-qubit states, the even-qubit
systems exhibit frozen phenomenon of quantum discord under the phase flip
channel, while the odd-qubit systems not.

##  IV summary

We have studied the quantum discord for a family to multi-qubit states.
Analytical formulae have been derived in detail for $(2v+1)$, $(4v-2)$ and
$(4v)$-qubit states. The level surfaces of quantum discord have been depicted.
It has been shown that under the phase flip channel the quantum discord could
still keep constant in a certain time interval for the even-qubit systems, but
not for odd-qubit systems. Our results may highlight further investigations on
multipartite quantum discord and their applications in quantum information
processing.

Acknowledgments This work is supported by NSFC under numbers 11765016,
12075159, 12065021, 11905131, Key Project of Beijing Municipal Commission of
Education (KZ201810028042), Beijing Natural Science Foundation (Z190005),
Academy for Multidisciplinary Studies, Capital Normal University, and Shenzhen
Institute for Quantum Science and Engineering, Southern University of Science
and Technology, Shenzhen 518055, China (No. SIQSE202001). X.B.L. was supported
by the Natural Science Foundation of Jiangxi Province (Grant No.
20202BAB201010).

B.L and C.L.Z contributed equally to this work.

## References

  * (1) G. Adesso, T.R. Bromley, and M. Cianciaruso, J. Phys. A 49, 473001 (2016).
  * (2) Z. Ma, J. Cui, Z. Cao, S. Fei, V. Vedral, T. Byrnes, and C. Radhakrishnan, Europhys. Lett 125, 50005 (2019).
  * (3) H. Ollivier and W. H. Zurek, Phys. Rev. Lett. 88, 017901 (2001).
  * (4) L. Henderson and V. Vedral, J. Phys. A 34, 6899 (2001).
  * (5) C.C. Rulli and M.S. Sarandy, Phys. Rev. A 84, 042109 (2011).
  * (6) S. Luo and S.Fu, Phys. Rev. A 82, 034302 (2010).
  * (7) M. Ali, J. Phys. A 43, 495303 (2010).
  * (8) S. Vinjanampathy and A. R. P. Rau, J. Phys. A 45, 095303 (2012).
  * (9) C. Xie, Z. Zhang, J. Chen, and X. Yin, Entropy 22, 147 (2020).
  * (10) S. Luo, Phys. Rev. A 77, 042303 (2008).
  * (11) M.D. Lang, and C.M. Caves, Phys. Rev. Lett. 105, 150501 (2010).
  * (12) M. Ali, A.R.P. Rau, and G. Alber, Phys. Rev. A 81, 042105 (2010); M. Ali, A.R.P. Rau, and G. Alber, ibid 82, 069902 (2010).
  * (13) C. Radhakrishnan, M. Laurière, and T. Byrnes, Phys. Rev. Lett. 124, 110401 (2020).
  * (14) T. Baumgratz, M. Cramer, and M.B. Plenio, Phys. Rev. Lett. 113, 140401 (2014).
  * (15) A. Streltsov, U. Singh, H.S. Dhar, M.N. Bera, and G. Adesso, Phys. Rev. Lett. 115, 020403 (2015).
  * (16) A.Streltsov, G.Adesso, and M.B. Plenio, Rev. Mod. Phys. 89, 041003 (2017).
  * (17) G. Styliaris, N. Anand, L.C. Venuti, and P. Zanardi, Phys. Rev. B 100, 224204 (2019).
  * (18) Z. Xi, and S. Yuwen, Phys. Rev. A 99, 022340 (2019).
  * (19) J. Xu, L. Shao, and S. Fei, Phys. Rev. A 102, 012411 (2020).
  * (20) C. Xiong, A, Kumar, M. Huang, S. Das, U. Sen, and J.Wu, Phys. Rev. A 99, 032305 (2019).
  * (21) Y. Huang, New J. Phys. 16 033027( 2014).
  * (22) M. Oszmaniec, P. Suwara and A. Sawicki, J. Math. Phys. 55, 062204 (2014).
  * (23) T. Theurer, S. Satyajit, and M.B. Plenio, Phys. Rev. Lett. 125, 130401 (2020).
  * (24) Y. Yao, D. Li, and C.P. Sun, Phys. Rev. A 100, 032324 (2019).
  * (25) D. Yu, L. Zhang, and C. Yu, Phys. Rev. A 101, 062114 (2020).
  * (26) Z. Jin and S. Fei, Phys. Rev. A 97, 062342 (2018).
  * (27) H. Zhu, Z. Ma, Z. Cao, S.M. Fei, V. Vedral, Phys. Rev. A 96, 032316 (2017).
  * (28) Y. Xi, T. Zhang, Z. Zheng, X. Li-Jost, S.M. Fei, Phys. Rev. A 100, 022310 (2019).
  * (29) M. Ban, Phys. Rev. A 99, 012116 (2019).
  * (30) L. Mazzola, J. Piilo, and S. Maniscalco, Phys. Rev. Lett. 104, 200401 (2010).
  * (31) J. Maziero, L. C. Céleri, R. M. Serra, and V. Vedral, Phys. Rev. A. 80, 044102 (2009).
  * (32) P. Haikka, T. H. Johnson, and S. Maniscalco, Phys. Rev. A 87, 010103 (2013).
  * (33) T.R. Bromley, M. Cianciaruso, and G. Adesso, Phys. Rev. Lett. 114, 210401 (2015).
  * (34) B. Li, Z. Wang, S. Fei, Phys. Rev. A 83, 022321 (2011)