# AirCrab: A Hybrid Aerial-Ground Manipulator with An Active Wheel

Muqing Cao*, Jiayan Zhao*, Xinhang Xu, and Lihua Xie1 *Equal ContributionAll
authors are with the School of Electrical and Electronic Engineering, Nanyang
Technological University, 50 Nanyang Avenue, Singapore. 1 Corresponding author

###### Abstract

Inspired by the behavior of birds, we present AirCrab, a hybrid aerial ground
manipulator (HAGM) with a single active wheel and a 3-degree of freedom
(3-DoF) manipulator. AirCrab leverages a single point of contact with the
ground to reduce position drift and improve manipulation accuracy. The single
active wheel enables locomotion on narrow surfaces without adding significant
weight to the robot. To realize accurate attitude maintenance using propellers
on the ground, we design a control allocation method for AirCrab that
prioritizes attitude control and dynamically adjusts the thrust input to
reduce energy consumption. Experiments verify the effectiveness of the
proposed control method and the gain in manipulation accuracy with ground
contact. A series of operations to complete the letters ‘NTU’ demonstrates the
capability of the robot to perform challenging hybrid aerial-ground
manipulation missions.

## Supplementary Material

A video summarizing the approach and experiments is available at
https://youtu.be/Q1n-IiIt400.

##  I INTRODUCTION

In recent years, new aerial robots have emerged that perform tasks requiring
physical interaction with the surroundings. These include aerial robots
equipped with (1) robot arms to pick and place objects [1], (2) contact-based
sensing devices for data collection [2], and (3) ground-driving wheels for
hybrid air-ground locomotion [3]. These emerging aerial robots open up many
application scenarios, such as item manipulation and transportation [4],
contact-based façade and bridge inspection [2], and long-duration exploration
of unknown regions. However, the need for high robustness and accuracy in the
interaction poses a significant challenge for aerial robots, whose control and
localization accuracy are affected by complex aerodynamic effects, wind
disturbances, and the limitation of the onboard sensors and computational
resources [5]. The highly coupled position and attitude of the traditional
multi-rotor robots also make it challenging to perform accurate contact.

Various solutions have been proposed in the context of contact-based
inspection. In [6, 7], fully actuated UAVs are designed to achieve decoupled
position and attitude control for precise sensing using fixed-mounted contact
probes. In a different research direction, multi-rotor robots are equipped
with several top-mounted wheels for inspection of ceiling or bridge floor [8,
9]. During the inspection, thrust from rotors pushes the robot against the
surface while the wheels in contact with the top surface realize the
locomotion. Such designs reduce the degree of freedom of the robots and
effectively alleviate the position drift and attitude oscillation during the
inspection. However, adding numerous active and passive wheels increases the
weight and the power consumption.

In the context of aerial manipulation (AM), a new class of robots called
hybrid aerial-ground manipulators (HAGM) has been proposed [1]. Incorporating
design principles from hybrid aerial-ground vehicles (HAGVs) [3, 10, 11],
HAGMs are equipped with ground-driving mechanisms to reduce energy consumption
and improve manipulation accuracy by exploiting surface contact. In [12], a
hexarotor UAV is equipped with a crawler and roller for locomotion on pipes
while inspecting the pipes by a 3-DoF arm, but the heavy driving mechanisms
compromise the aerial efficiency. In [13, 14], a HAGM design is proposed with
two active wheels and two passive wheels but no physical robot is built. Chat-
PM[15] combines a dual-passive-wheel HAGV [16] with a lightweight 4-DoF arm to
achieve the state-of-the-art operation accuracy. However, the use of passive
wheels means that attitude and position are still highly coupled in ground
motion, reducing its flexibility in manipulation tasks. In general, the
existing HAGMs employ a multi-wheel design requiring multiple contact points,
hindering applications on narrow surfaces.

Figure 1: The AirCrab mimics the nesting behavior of egrets, standing on a
single active wheel, maintaining attitude with propellers, and performing
precise operation.

In this work, we propose AirCrab, a design of HAGM with only one active wheel.
Such a design is inspired by the behavior of egrets who can stand on a single
leg to provide an anchor for precise manipulation and flap wings to maintain
balance, as shown in Fig. 1. The lightweight active wheel provides a single
contact point on narrow terrain and does not add much redundant weight in the
aerial operation. The single contact point facilitates rotation around the
robot’s vertical axis, adding a degree of freedom for manipulation.
Furthermore, the active wheel enables decoupled position and attitude control,
where the propellers are used to maintain attitude while the translational
motion is realized by wheel rotation. To reduce energy consumption and improve
control performance on the ground, we design a control allocation approach for
the propellers that prioritizes attitude maintenance and dynamically selects
the minimum necessary propeller thrust. The main contributions of this paper
are summarized as:

  * •

A design of a HAGM with an active wheel and a 3-DoF manipulator is proposed
and built;

  * •

Detailed dynamics modeling and controller design are carried out, and a
control allocation strategy is proposed to prioritize attitude control and
reduce propeller thrust;

  * •

Experimental results verify the high accuracy and energy efficiency of the
proposed control allocation strategy and the trajectory tracking accuracy of
the end effector in different modes;

  * •

Challenging missions showcase the capability of AirCrab in grasping an object
on the ground, taking off with payload, landing on a table, and placing the
object to complete the characters ‘NTU’ (Fig. 7).

The rest of the paper is organized as follows: Section II describes the
mechanical design of AirCrab. The detailed dynamics model is shown in the
section III. Section IV explains the attitude controller design and presents
the proposed control allocation strategy. Section V details the experimental
verification in terms of attitude control, power consumption, manipulator
tracking performance, and the execution of complex hybrid aerial-ground
missions. Finally, section VI concludes the paper.

##  II Mechanical Design

We construct a prototype based on a general quadrotor design with a diagonal
length of $450$ mm, equipped with four brushless HS$2216$ $920$KV motors and
$10$-inch propellers. The onboard Cube Orange flight controller runs a
customized Ardupilot firmware for attitude and position control. The total
weight of the robot including the wheel, the manipulator, the quadcopter and a
4S 5200 mAh battery is 2.655 kg.

Figure 2: Coordinate frames for AirCrab. The inertial frame is labeled with a
subscript $I$ and the body frame is labeled with a subscript $b$.

A servo-driven rubber wheel is connected to the bottom plate of the quadrotor
chassis through a carbon fiber tube, as shown in Fig 1. A lightweight
servomotor STS3020 with 24.3 kg·cm maximum torque capability is chosen. To
ensure the wheel lies directly below the geometric center of the quadrotor, we
use a belt and two pulleys to drive the wheel instead of connecting the servo
directly to the wheel. The total weight of the wheel, the servomotor, and the
mechanical connections is $305$ g.

The $3$-DoF manipulator arm utilizes three rotational joints to offer 2 DoF in
position and 1 DoF in rotation for the end effector. The rotation axes for the
joints are shown in Figure 2. A gripper is mounted as the end effector for
object pick and place. Each joint and gripper is driven by a servomotor and
connected to an ESP32-based driver board for control. Two extra carbon fiber
plates are mounted on the end effector to provide two widely spaced support
points, forming a triangular support for AirCrab when powered off.

##  III Dynamics Modelling

The translational and rotational dynamics of AirCrab can be expressed as:

$\displaystyle
m\ddot{\bm{p}}=\bm{f}_{g}+\bm{R}(\bm{f}+\bm{f}_{w}+\bm{f}_{d}),$ (1)
$\displaystyle\bm{J}\dot{\bm{\omega}}+[\bm{\omega}]_{\times}\bm{J}\bm{\omega}=\bm{\tau}-\bm{\tau}_{w}+\bm{\tau}_{d}+[\bm{d}_{w}]_{\times}\bm{f}_{w}+[\bm{d}_{g}]_{\times}\bm{f}_{g},$
(2)

where $\bm{p}=[p_{x},p_{y},p_{z}]^{\top}$ denotes the position of the robot in
the inertial frame; $\bm{\omega}\in\mathbb{R}^{3}$ denotes the angular
velocity in the body frame; $\bm{R}$ denotes the rotation matrix from the
robot body frame to the inertial frame; $m\in\mathbb{R}$ and
$\bm{J}\in\mathbb{R}^{3\times 3}$ denote the mass and the moment of inertial
about the body frame of the robot; $\bm{f}=[0,0,T]^{\top}$ and
$\bm{\tau}\in\mathbb{R}^{3}$ denote the collective thrust and torque generated
by the rotors in the body frame; $\bm{f}_{w}$ denotes the force acting on the
wheel due to contact with the ground and $\bm{\tau}_{w}=[0,\tau_{w},0]^{\top}$
represents the torque exerted by the servomotor on the wheel; $\bm{f}_{d}$ and
$\bm{\tau}_{d}$ represent the force and torque arising from disturbances due
to the movements of the manipulator and aerodynamic effects; $\bm{d}_{w}$ and
$\bm{d}_{g}$ denote the vector from the robot’s center of gravity (CG) and
ground contact point to the origin of its body frame, expressed in body and
inertial frame, respectively. $[\cdot]_{\times}$ denotes the screw-symmetric
cross product of a vector.

When the robot is in the aerial mode, the wheel is not in contact with the
ground, $\bm{f}_{w}$ and $\bm{\tau}_{w}$ are equal to zero. Hence the dynamics
model (1) and (2) becomes identical to those of a traditional multi-rotor UAV
[17].

When the robot is in contact with the ground, $\bm{f}_{w}$ and $\bm{\tau}_{w}$
are nonnegligible. We specifically consider the case where the attitude of the
robot is maintained near the level condition, i.e., roll and pitch angles
$\phi,\theta\approx 0$. In practice, this can be achieved by controlling the
rotor speed, as will be shown in Section V. In this case, $\bm{R}$ can be
approximated as a yaw-only rotation around the inertial $z$-axis,
$\bm{R}_{z}(\gamma)$. As a result, $\bm{R}\bm{f}\approx[0,0,T]^{\top}$, the
effect of rotor thrust on the x and y movement of the robot is negligible. The
2D translation of the robot on the ground can be approximated as

$\displaystyle m\begin{bmatrix}\ddot{p}_{x}\\\
\ddot{p}_{y}\end{bmatrix}\approx\begin{bmatrix}\cos{\gamma}&-\sin\gamma\\\
\sin\gamma&\cos\gamma\end{bmatrix}\left(\begin{bmatrix}f_{wx}\\\
f_{wy}\end{bmatrix}+\begin{bmatrix}f_{dx}\\\ f_{dy}\end{bmatrix}\right),$ (3)

where $f_{wx},f_{wy},f_{dx},f_{dy}$ are the $x$ and $y$ components of
$\bm{f}_{w}$ and $\bm{f}_{d}$, respectively. $f_{wx}$ and $f_{wy}$ indicate
the frictional force at the ground contact point. When the wheel is
stationary, small disturbances acting on the robot are compensated by the
static frictional forces, reducing the translational drifts due to
disturbances.

When the robot maintains a level attitude, the ground contact point lies close
to the $z$-axis of the body frame, i.e., $\bm{d}_{w}\approx[0,0,-l_{z}]$ for a
$l_{z}>0$. Referring to Equation (2), the effect of $\bm{f}_{w}$ and
$\bm{f}_{g}$ on the yaw rotation can be neglected. Furthermore, given that the
torque input $\bm{\tau}_{w}$ does not apply on the body $z$-axis (Figure 2),
the yaw rotation in ground mode is similar to that in aerial mode, dominated
by the torques from rotors and disturbances. Clearly, the influence of ground
contact on the roll and pitch dynamics is more significant than on the yaw
dynamics. Finally, the wheel dynamics and kinematics can be expressed as:

$\displaystyle J_{w}\dot{\omega}_{w}=\tau_{w}-rf_{wx},$ (4)
$\displaystyle\omega_{w}r=\dot{p}_{wx},$ (5)

where $\omega_{w}$ and $r$ are the rotation rate and the radius of the wheel,
$p_{wx}$ is the position of the wheel in the body $x$-axis, $J_{w}$ indicates
the moment of inertia of the wheel about its rotation axis.

##  IV Controller Design

###  IV-A Attitude Control

The control inputs for AirCrab are the rotor thrust $T$, the rotor torques
$\bm{\tau}$, and the servomotor torque $\tau_{w}$. In the aerial mode, desired
thrust $T$ and the desired attitude $\bm{R}_{d}$ are computed from the high-
level position controller similar to [17]. In the ground mode, the desired
attitude is set to a level condition $\bm{R}_{d}=\bm{R}_{z}(\gamma_{d})$,
where $\gamma_{d}$ is the desired yaw angle obtained from pilot input or high-
level planner. The thrust in ground mode is dynamically adjusted in the
control allocation as will be explained in Section IV-B. Given the desired
attitude, the input torque is computed using a cascade PID controller:

$\displaystyle\bm{\omega}_{d}=\bm{K}_{p}\frac{1}{2}(\bm{R}^{\top}_{d}\bm{R}-\bm{R}^{\top}\bm{R}_{d})^{\vee},$
(6) $\displaystyle\bm{\omega}_{e}=\bm{\omega}_{d}-\bm{\omega},$ (7)
$\displaystyle\bm{\tau}=\bm{K}^{\omega}_{p}\bm{\omega}_{e}+\bm{K}^{\omega}_{i}\int\bm{\omega}_{e}+\bm{K}^{\omega}_{d}\dot{\bm{\omega}}_{e}+\bm{\tau}_{c},$
(8)

where $\bm{\omega}_{d}$ and $\bm{\omega}_{e}$ are the desired angular velocity
and the angular velocity error,
$\bm{K}_{p},\bm{K}^{\omega}_{p},\bm{K}^{\omega}_{i},\bm{K}^{\omega}_{d}$ are
the diagonal matrices consisting of the PID gains, ∨ is the mapping from
$so(3)$ to $\mathbb{R}^{3}$. $\bm{\tau}_{c}\in\mathbb{R}^{3}$ is added to
compensate for the off-center CG. For simplicity, we set $\bm{\tau}_{c}$ to be
constant in this work, however, techniques from existing literature [18] can
be employed to compensate for changing disturbances. The torque to drive the
wheel is computed from the desired rotation rate $\omega_{wd}$:

$\displaystyle\tau_{w}=K\omega_{wd},$ (9)

where $K$ is a proportion gain. $\omega_{wd}$ is obtained from the high-level
planner or the pilot input.

###  IV-B Control Allocation in Ground Mode

Among the computed control inputs, $\tau_{w}$ can be converted directly to a
signal and sent to the servomotor. On the other hand, $T$ and $\bm{\tau}$ have
to be collectively realized by four propellers due to the coupling between
row, pitch, yaw, and thrust. For a frame configuration as AirCrab (commonly
called an X frame), there exists a linear mapping from quadratic motor speeds
of the four motors $\bm{v}=[v_{1}^{2},v_{2}^{2},v_{3}^{2},v^{2}_{4}]^{\top}$
to control inputs $\bm{u}=[T,\tau^{\top}]^{\top}\in\mathbb{R}^{4}$, more
specifically, $\bm{u}=\bm{M}\bm{v}$ can be written as

$\displaystyle\begin{bmatrix}T\\\ \tau_{x}\\\ \tau_{y}\\\
\tau_{z}\end{bmatrix}=\begin{bmatrix}C_{t}&C_{t}&C_{t}&C_{t}\\\
-\frac{\sqrt{2}}{2}C_{t}l&-\frac{\sqrt{2}}{2}C_{t}l&\frac{\sqrt{2}}{2}C_{t}l&\frac{\sqrt{2}}{2}C_{t}l\\\
-\frac{\sqrt{2}}{2}C_{t}l&\frac{\sqrt{2}}{2}C_{t}l&-\frac{\sqrt{2}}{2}C_{t}l&\frac{\sqrt{2}}{2}C_{t}l\\\
-C_{d}&C_{d}&C_{d}&-C_{d}\end{bmatrix}\begin{bmatrix}v_{1}^{2}\\\ v_{2}^{2}\\\
v_{3}^{2}\\\ v_{4}^{2}\end{bmatrix}$ (10)

where $C_{t}>0$ and $C_{d}>0$ are the thrust coefficient and the drag
coefficient, respectively, and $l$ is the distance between propellers and
geometry center of frame.

It may seem straightforward to apply the inverse mapping
$\bm{v}=\bm{M}^{-1}\bm{u}$ to attain the motor speeds; however, the results
may exceed the mechanical limits of the motors, as each motor is subject to
velocity constraints, $0\leqslant v_{min}\leqslant v_{i}\leqslant v_{max}$.
Therefore, a control allocation problem is solved to find a modified control
input $\hat{\bm{u}}=\bm{g}(\bm{u})$ such that
$\bm{v}=\bm{M}^{-1}\hat{\bm{u}}\in\mathbb{V}$, where
$\mathbb{V}=\left\\{\bm{v}\in\mathbb{R}^{4}|v_{min}^{2}\leqslant
v^{2}_{i}\leqslant v_{max}^{2}\right\\}$. Common control allocation approaches
are weighted least square (WLS) [19], direct control allocation (DCA)[20] and
partial control allocation (PCA) [21]. Our approach is based on the PCA which
groups the inputs into different priorities and preserves the inputs with the
higher priority as much as possible. We choose to classify $\bm{u}$ into the
tilt input ($\tau_{x},\tau_{y}$), yaw input ($\tau_{z}$), and the thrust input
($T$), in the order of decreasing priority. The motivation is that more
control effort should be allocated to the roll and pitch, which faces larger
disturbances from the ground contact as analyzed in Section III. Furthermore,
with ground support, the thrust from the motors is not necessary to maintain
the position and should be minimized to reduce energy consumption.

The control allocation strategy is presented in Algorithm 1. By mapping the
feasible set of motor velocity $\mathbb{V}$ using Equation (10), we obtain the
feasible set of control inputs $\mathbb{U}$. The algorithm first checks if the
tilt control inputs $\tau_{x},\tau_{y}$ stay within the feasible subset
$\mathbb{B}$, which is obtained by projecting the feasible set $\mathbb{U}$ on
$T=0,\tau_{y}=0$ using Fourier-Motzkin elimination[22]. Specifically,
$\tau_{x},\tau_{y}$ should satisfy

$\displaystyle\sqrt{2}C_{t}l(v_{min}^{2}-v_{max}^{2})\leqslant\tau_{x}\pm\tau_{y}\leqslant\sqrt{2}C_{t}l(v_{max}^{2}-v_{min}^{2}).$
(11)

If (11) is not satisfied, a scale factor $\alpha\in(0,1)$ is computed such
that $\alpha[\tau_{x},\tau_{y}]\in\partial\mathbb{B}$, i.e.,
$\alpha[\tau_{x},\tau_{y}]$ lies on the boundary of $\mathbb{B}$ to preserve
the magnitude of the original input as much as possible. Then, we check if the
scaled tilt input and the original yaw input,
$[\alpha\tau_{x},\alpha\tau_{y},\tau_{z}]$, lie in the feasible set
$\mathbb{A}$ obtained by projecting the feasible set $\mathbb{U}$ on $T=0$:

$\displaystyle v_{min}^{2}-v_{max}^{2}\leqslant\left\\{\begin{aligned}
&\frac{\sqrt{2}\alpha}{2C_{t}l}\tau_{x}\pm\frac{1}{2C_{d}}\tau_{z}\\\
&\frac{\sqrt{2}\alpha}{2C_{t}l}\tau_{y}\pm\frac{1}{2C_{d}}\tau_{z}\end{aligned}\right\\}\leqslant
v_{max}^{2}-v_{min}^{2}.$ (12)

Similarly, if (12) is not satisfied, we scale $\tau_{z}$ with a factor $\beta$
such that
$[\alpha\tau_{x},\alpha\tau_{y},\beta{\tau}_{z}]\in\partial\mathbb{A}$. Next,
we compute the feasible range for $T$ given the torque inputs
$\alpha\tau_{x},\alpha\tau_{y},\beta{\tau}_{z}$:

$\displaystyle\begin{aligned} &T\geq
4C_{t}v_{min}^{2}-\min\left\\{-\frac{\sqrt{2}}{l}\alpha|\tau_{x}|\pm\frac{\sqrt{2}}{l}\alpha\tau_{y}\pm\frac{C_{t}}{C_{d}}\beta\tau_{z}\right\\}\\\
&T\leq
4C_{t}v_{max}^{2}-\min\left\\{-\frac{\sqrt{2}}{l}\alpha|\tau_{x}|\pm\frac{\sqrt{2}}{l}\alpha\tau_{y}\pm\frac{C_{t}}{C_{d}}\beta\tau_{z}\right\\}.\end{aligned}$
(13)

In theory, we can simply set $T$ to the lower bound of Equation (13) to
minimize thrust input. In practice, we add a small bias $T_{\text{ground}}$ to
$T$ because the thrust coefficient $C_{t}$, obtained by fitting the full
thrust curve, does not approximate the thrust well when the motors operate at
a very low velocity, causing larger control error as will be shown in the
experiments (Section V-A). Finally, the motor velocities are obtained from the
inverse mapping of the modified control inputs, which are guaranteed to
satisfy the velocity constraints.

1if  _$[\tau_{x},\tau_{y}]^{T}\in\mathbb{B}$_ then

2 $\alpha=1$

3else

4 find $\alpha\in(0,1)$ s.t.
$\alpha[\tau_{x},\tau_{y}]^{\top}\in\partial\mathbb{B}$

5if  _$[\alpha\tau_{x},\alpha\tau_{y},\tau_{z}]^{T}\in\mathbb{A}$_ then

6 $\beta=1$

7else

8 find $\beta$ s.t.
$[\alpha\tau_{x},\alpha\tau_{y},\beta{\tau}_{z}]^{\top}\in\partial\mathbb{A}$

9set $T$ to the lower bound of Eqn. (13)

$T=T+T_{\text{ground}}$ // add small bias

$\bm{v}=\bm{M}^{-1}[T,\alpha\tau_{x},\alpha\tau_{y},\beta{\tau}_{z}]^{\top}$
// inverse map

10

Algorithm 1 Control allocation algorithm

##  V Experiments

###  V-A Attitude Control and Power Consumption

Figure 3: The response of proposed method with
${T_{\text{ground}}}/{T_{\text{max}}}=7.5\%$, the settling time
$t_{s}=0.9195s$. In this figure, $0^{\circ}$ is the desired state and the
region of $\pm 3^{\circ}$ error is colored in green. Figure 4: Pitch angle
error and power consumption under various throttle levels.

We first conduct experiments to validate the attitude control performance of
AirCrab in ground mode. During the experiment, the robot starts at a rest
position with a large pitch angle and tries to reach the level condition
($\phi,\theta=0$). The attitude and power consumption of the robot are
recorded with varying values of $T_{\text{ground}}$. For comparison, we also
conduct the same experiment using the default attitude control algorithm of
the widely used open-source Ardupilot firmware
111https://ardupilot.org/copter/ and vary the hover throttle level (by
adjusting the parameter MOT_THST_HOVER). Figure 3 shows the attitude of the
robot for ${T_{\text{ground}}}/{T_{\text{max}}}=7.5\%$. The system takes
$0.9195$ s to reach and stay within $\pm 3^{\circ}$, and the pitch error
continues to reduce and maintain within $\pm 1^{\circ}$ after $5.26$ s. The
root mean square error (RMSE) of the pitch angle for $10$ s is computed for
each experiment. Figure 4 plots the trend of the power consumption and the
RMSE for the proposed method and the Ardupilot method. As we increase
$T_{\text{ground}}$, the pitch RMSE decreases while the power consumption
increases for the proposed method. As explained in Section IV-B, the
improvement in accuracy is due to a better approximation of the thrust using
the experimentally determined thrust coefficient at higher motor velocity. In
comparison, using the Ardupilot controller, the pitch angle could not
stabilize within $\pm 3^{\circ}$ for a hovel throttle smaller than $17.5\%$.
This is because the control allocation prioritizes maintaining the altitude
over attitude, which is reasonable for aerial operation. However, when the
thrust input is small, the control allocation scales down the rotation torque
to maintain a small total thrust, causing insufficient torque to bring the
robot to a level condition. When the hover throttle is above $17.5\%$, the
pitch RMSE reduces but is still higher than the proposed method, and the
energy consumption is higher. Compared with the proposed method at
${T_{\text{ground}}}/{T_{\text{max}}}=7.5\%$, the Ardupilot method produces
$180\%,96.2\%$, and $12.1\%$ more error and consumes $11.7\%,23.5\%$, and
$47.1\%$ more energy at hover throttle of $17.5\%,20\%$, and $22.5\%$,
respectively. Hence, the proposed attitude controller and the control
allocation method can achieve good control performance and high energy
efficiency. For subsequent experiments, we choose
${T_{\text{ground}}}/{T_{\text{max}}}=7.5\%$ to exploit the good energy
efficiency without much compromise on the control accuracy.

Figure 5: The time-lapse photos of the manipulation experiment. The yellow
points and red arrows illustrate the movement of the end effector. Figure 6:
The attitude of the AirCrab and the displacement of the object during the
mission. Different phases of the mission are colored and numbered as Fig.5.

We then conduct an experiment in which the robot is manually controlled to
move forward using the wheel, pick up an object of $90$ gram, move back and
rotate $90^{\circ}$, and then place the object. Textured mats are used as
uneven ground surfaces with bumps and ridges. The attitude of the robot and
the displacement of the object are recorded using the Vicon motion capture
system and plotted in Figure 6. The time-lapse images of the experiment are
shown in Figure 5. It is observed that the roll and pitch angles maintained
within $\pm 3^{\circ}$ for most of the operation, and the yaw angle stayed
within $3^{\circ}$ error to the desired yaw. A deviation of $3.9^{\circ}$ is
observed when the robot picks up the object, but the error reduces to less
than $1^{\circ}$ in $0.5$ s. Oscillations of pitch angles are observed when
the robot moves (phases $4$ and $5$) because the disturbance caused by the
wheel rotation and uneven ground is not compensated in the controller. The
effect of the disturbance on roll and yaw is less significant as the wheel
rotation overcomes friction forces mainly in the body $X$-axis.

###  V-B Manipulator Accuracy

Figure 7: Overview of the hybrid aerial-ground manipulation experiments. In
each experiment, one wooden is placed to complete one of the letters in ‘NTU’.
Figure 8: Comparison of tracking performance of the end effector in different
modes.

We measure the tracking accuracy of the end effector in three cases: (1)
completely static, (2) ground mode, holding level attitude, and (3) aerial
mode, holding position. In aerial mode, the position estimate is obtained from
an extended Kalman filter that fuses the onboard IMU and optical flow sensor.
We use onboard localization to simulate a realistic manipulation task where
accurate external localization systems are unavailable. In each case, the
robot arm is commanded to execute a trajectory such that the end effector
reaches $2$ setpoints sequentially and holds for $7$ seconds at each setpoint.
The trajectory setpoints are given in the robot arm’s base frame (labeled with
the subscript $a0$ in Fig. 2), and inverse kinematics is solved to obtain the
target angles for each joint at $5$ Hz. Figures 8 show the tracking
performance for all modes. As expected, the tracking accuracy is best in
static mode, with $0.90$ cm RMSE and $2.11$ cm max error. The tracking error
in static mode arises due to the limitations of the robot arm, specifically,
the low angular resolution of the low-cost servomotor joint and imperfect
parameters used to compute kinematics. The RMSE and maximum error in ground
mode are $1.02$ cm and $2.14$ cm, respectively, slightly increased from those
in static mode. On the other hand, the RMSE and maximum error in aerial mode
are significantly higher at $2.00$ cm and $5.53$ cm, respectively. From Figure
8, it is observed that the tracking performance in ground mode is similar to
that in the static mode. Much larger position oscillation is observed in
aerial mode, mainly due to drift in the position estimate and air disturbance.
Therefore, ground contact is indeed beneficial in improving the manipulator’s
accuracy. The high accuracy in ground mode enables challenging manipulation
tasks as described in the next section.

###  V-C Hybrid Aerial-Ground Operation

We conduct a series of hybrid aerial-ground manipulation tasks to demonstrate
the capabilities of AirCrab. In each experiment, the robot is manually
controlled to pick up a wooden block, take off and fly over a table, land, and
place the block to complete one of the letters of the word NTU. A pilot
remotely controls the ground and aerial locomotion, while another operator
controls the joint angles of the arm. As a safety precaution, the robot is
connected to a long rod through a safety tether, which does not exert force on
the robot except during emergencies. Despite the manual operation, high
accuracy and stability of the end effector are required to place the block in
the designated location. Otherwise, the block may fall from the other blocks
or fail to form the desired pattern. Figure 7 shows an overview of the tasks;
a video of the experiments can be viewed in the supplementary document.

##  VI CONCLUSIONS

In this work, we present the dynamics modeling and control design of AirCrab,
a HAGM with a single active wheel. Experiments have demonstrated the attitude
control performance and the improvement in manipulation accuracy when using
the proposed control strategy on the ground. This work serves as a starting
point for our subsequent research. First, the control performance of the robot
on various types of terrain could be enhanced by estimating the terrain slope
and roughness. Second, active disturbance rejection techniques could improve
control accuracy under moving wheels and changing payload. Eventually, we
would like to realize a fully autonomous hybrid aerial-ground manipulation
mission with onboard sensors and computational resources.

## References

  * [1] A. Ollero, M. Tognon, A. Suarez, D. Lee, and A. Franchi, “Past, present, and future of aerial robotic manipulators,” IEEE Transactions on Robotics, vol. 38, no. 1, pp. 626–645, 2021.
  * [2] A. Jimenez-Cano, J. Braga, G. Heredia, and A. Ollero, “Aerial manipulator for structure inspection by contact from the underside,” in 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 1879–1884, 2015.
  * [3] M. Cao, X. Xu, S. Yuan, K. Cao, K. Liu, and L. Xie, “Doublebee: A hybrid aerial-ground robot with two active wheels,” in 2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 6962–6969, IEEE, 2023.
  * [4] A. X. Appius, E. Bauer, M. Blöchlinger, A. Kalra, R. Oberson, A. Raayatsanati, P. Strauch, S. Suresh, M. von Salis, and R. K. Katzschmann, “Raptor: Rapid aerial pickup and transport of objects by robots,” in 2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 349–355, IEEE, 2022.
  * [5] A. Suarez, V. M. Vega, M. Fernandez, G. Heredia, and A. Ollero, “Benchmarks for aerial manipulation,” IEEE Robotics and Automation Letters, vol. 5, no. 2, pp. 2650–2657, 2020.
  * [6] R. Rashad, D. Bicego, R. Jiao, S. Sanchez-Escalonilla, and S. Stramigioli, “Towards vision-based impedance control for the contact inspection of unknown generically-shaped surfaces with a fully-actuated uav,” in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 1605–1612, IEEE, 2020.
  * [7] K. Bodie, M. Brunner, M. Pantic, S. Walser, P. Pfändler, U. Angst, R. Siegwart, and J. Nieto, “Active interaction force control for contact-based inspection with a fully actuated aerial vehicle,” IEEE Transactions on Robotics, vol. 37, no. 3, pp. 709–722, 2021.
  * [8] R. Ladig and K. Shimonomura, “High precision marker based localization and movement on the ceiling employing an aerial robot with top mounted omni wheel drive system,” in 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 3081–3086, IEEE, 2016.
  * [9] H.-F. Wang, L. Zhai, H. Huang, L.-M. Guan, K.-N. Mu, and G.-p. Wang, “Measurement for cracks at the bottom of bridges based on tethered creeping unmanned aerial vehicle,” Automation in Construction, vol. 119, p. 103330, 2020.
  * [10] N. Pan, J. Jiang, R. Zhang, C. Xu, and F. Gao, “Skywalker: A compact and agile air-ground omnidirectional vehicle,” IEEE Robotics and Automation Letters, vol. 8, no. 5, pp. 2534–2541, 2023.
  * [11] Z. Zheng, J. Wang, Y. Wu, Q. Cai, H. Yu, R. Zhang, J. Tu, J. Meng, G. Lu, and F. Gao, “Roller-quadrotor: A novel hybrid terrestrial/aerial quadrotor with unicycle-driven and rotor-assisted turning,” in 2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 6927–6934, IEEE, 2023.
  * [12] A. Lopez-Lora, P. Sanchez-Cuevas, A. Suarez, A. Garofano-Soldado, A. Ollero, and G. Heredia, “Mhyro: Modular hybrid robot for contact inspection and maintenance in oil & gas plants,” in 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 1268–1275, 2020.
  * [13] M. Saeed, M. Fanni, and A. M. Mohamed, “Modeling and control of a new aerial manipulation system of a hybrid ground aerial robot,” in 2018 IEEE International Conference on Industrial Technology (ICIT), pp. 234–238, IEEE, 2018.
  * [14] M. Salman, A. Sameh, M. Fanni, S. Sugano, and A. M. Mohamed, “Design, control, and dynamic simulation of securing and transformation mechanisms for a hybrid ground aerial robot,” International Journal of Mechanical and Mechatronics Engineering, vol. 20, no. 2, pp. 100–107, 2020.
  * [15] Y. Ding, X. Ji, L. Zhang, Y. Zhu, Y. Dong, T. Wu, C. Han, and J. Yang, “Chat-pm: A class of composite hybrid aerial/terrestrial precise manipulator,” IEEE Robotics and Automation Letters, 2024.
  * [16] R. Zhang, Y. Wu, L. Zhang, C. Xu, and F. Gao, “Autonomous and adaptive navigation for terrestrial-aerial bimodal vehicles,” IEEE Robotics and Automation Letters, vol. 7, no. 2, pp. 3008–3015, 2022.
  * [17] D. Mellinger and V. Kumar, “Minimum snap trajectory generation and control for quadrotors,” in 2011 IEEE International Conference on Robotics and Automation, pp. 2520–2525, 2011.
  * [18] H. Cao, Y. Li, C. Liu, and S. Zhao, “Eso-based robust and high-precision tracking control for aerial manipulation,” IEEE Transactions on Automation Science and Engineering, pp. 1–17, 2023.
  * [19] Q.-L. Zhou, Y. Zhang, C.-A. Rabbath, and D. Theilliol, “Reconfigurable control allocation technology using weighted least squares for nonlinear system in unmanned aerial vehicle,” in 48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, p. 1138, 2010.
  * [20] W. C. Durham, “Constrained control allocation,” Journal of Guidance, control, and Dynamics, vol. 16, no. 4, pp. 717–725, 1993.
  * [21] J. C. Monteiro, F. Lizarralde, and L. Hsu, “Optimal control allocation of quadrotor uavs subject to actuator constraints,” in 2016 American Control Conference (ACC), pp. 500–505, IEEE, 2016.
  * [22] G. M. Ziegler, Lectures on polytopes, vol. 152.  Springer Science & Business Media, 2012.


††thanks: ,† First three authors contributed equally.††thanks: ,† First three
authors contributed equally.

# Demonstration of multi-time quantum statistics without measurement back-
action

Pengfei Wang  State Key Laboratory of Low Dimensional Quantum Physics,
Department of Physics, Tsinghua University, Beijing 100084, China  Beijing
Academy of Quantum Information Sciences, Beijing 100193, China  Hyukjoon Kwon
<EMAIL_ADDRESS>Korea Institute for Advanced Study, Seoul 02455, Korea
Chun-Yang Luan  State Key Laboratory of Low Dimensional Quantum Physics,
Department of Physics, Tsinghua University, Beijing 100084, China  Wentao Chen
State Key Laboratory of Low Dimensional Quantum Physics, Department of
Physics, Tsinghua University, Beijing 100084, China  Mu Qiao  State Key
Laboratory of Low Dimensional Quantum Physics, Department of Physics, Tsinghua
University, Beijing 100084, China  Zinan Zhou  State Key Laboratory of Low
Dimensional Quantum Physics, Department of Physics, Tsinghua University,
Beijing 100084, China  Kaizhao Wang  State Key Laboratory of Low Dimensional
Quantum Physics, Department of Physics, Tsinghua University, Beijing 100084,
China  M. S. Kim<EMAIL_ADDRESS>Korea Institute for Advanced Study,
Seoul 02455, Korea  QOLS, Blackett Laboratory, Imperial College London, London
SW7 2AZ, United Kingdom  Kihwan Kim<EMAIL_ADDRESS>State Key
Laboratory of Low Dimensional Quantum Physics, Department of Physics, Tsinghua
University, Beijing 100084, China  Beijing Academy of Quantum Information
Sciences, Beijing 100193, China  Frontier Science Center for Quantum
Information, Beijing 100084, China

###### Abstract

It is challenging to obtain quantum statistics of multiple time points due to
the principle of quantum mechanics that a measurement disturbs the quantum
state. We propose an ancilla-assisted measurement scheme that does not suffer
from the measurement-induced back-action and experimentally demonstrate it
using dual-species trapped ions. By ensemble averaging the ancilla-measurement
outcomes with properly chosen weights, quantum statistics, such as quantum
correlation functions and quasi-probability distributions can be
reconstructed. We employ ${}^{171}\rm{Yb}^{+}$-${}^{138}\rm{Ba}^{+}$ ions as
the system and the ancilla to perform multi-time measurements that consist of
repeated initialization and detection of the ancilla state without effecting
the system state. The two- and three-time quantum correlation functions and
quasi-probability distributions are clearly revealed from experimental data.
We successfully verify that the marginal distribution is unaffected by the
measurement at each time and identify the nonclassicality of the reconstructed
distribution. Our scheme can be applied for any $N$-time measurements of a
general quantum process, which will be an essential tool for exploring
properties of various quantum systems.

A striking difference between quantum and classical statistics arises as there
exist observables in a quantum mechanics that cannot be precisely determined
simultaneously. Correlations between these incompatible observables do not
follow the description of classical joint probability distributions but can be
explained by introducing quasi-probability distributions (QPDs), which allows
negative Wigner (1932); Margenau and Hill (1961) or even non-real values
Kirkwood (1933); Dirac (1945). Such nonclassical features of QPDs have been
widely studied in the fields of quantum foundations BELL (1966); Leggett and
Garg (1985); Feynman (1987) and thermodynamics Allahverdyan (2014); Perarnau-
Llobet  _et al._ (2017); Yunger Halpern (2017); Lostaglio (2018); Levy and
Lostaglio (2020); Kwon and Kim (2019), as well as being considered as a useful
resource for quantum computing Veitch  _et al._ (2012); Mari and Eisert
(2012); Howard  _et al._ (2014); Pashayan  _et al._ (2015); Rahimi-Keshari
_et al._ (2016) and metrology Kwon  _et al._ (2019); Arvidsson-Shukur  _et
al._ (2020); Lostaglio (2020). At the same time, quantum correlation functions
and QPDs serve as essential tools for exploring both static and dynamic
properties of various quantum systems Kubo (1957); Mandel and Wolf (1995);
Gardiner  _et al._ (2004); Clerk  _et al._ (2010); Krumm  _et al._ (2016);
Yunger Halpern  _et al._ (2018); Dowling  _et al._ (2021).

Meanwhile, accessing the quantum statistics in experiments is way more
challenging than the classical system. Quantum theory prohibits one from
directly obtaining the statistics of incompatible observables as a fundamental
trade-off relation holds between information gain and measurement-induced
disturbance Groenewold (1971); Fuchs and Peres (1996); Ozawa (2003); Buscemi
_et al._ (2014); Fan  _et al._ (2015); Lee  _et al._ (2021); Hong  _et al._
(2022). In particular, projection measurements are not sufficient to fully
capture the quantum nature as all the off-diagonal elements are washed out
after the measurement and can not contribute to the subsequent measurement
statistics. To detour this problem, various indirect methods to obtain quantum
correlation functions and QPDs have been theoretically proposed Buhrman  _et
al._ (2001); Somma  _et al._ (2002); Johansen (2007); Buscemi  _et al._
(2013); Pedernales  _et al._ (2014); Perarnau-Llobet  _et al._ (2017) and
experimentally demonstrated Souza  _et al._ (2011); Piacentini  _et al._
(2016); Xin  _et al._ (2017); Ringbauer  _et al._ (2018); Wu  _et al._ (2019);
Del Re  _et al._ (2022). While most of these approaches focus on quantum
statics of two-time points, only a few experimental realizations of quantum
correlation functions beyond two-time points are reported Xin  _et al._
(2017); Ringbauer  _et al._ (2018).

In this Letter, we propose and experimentally verify that quantum statistics
of multiple time points can be reconstructed from ancilla-assisted
measurements. Remarkably, the reconstructed quantum statistics do not suffer
back-action from the measurement-induced disturbance. To this end, we interact
the system with an ancilla state at each time followed by the ancilla
measurement with respect to a certain basis set. The quantum statistics is
then reconstructed by taking an ensemble average of the sequential measurement
outcomes with properly weights. By increasing the number of the ancilla
measurement bases, richer quantum statistics, such as quantum correlation
function and QPDs, can be obtained.

Our approach has several advantages: (i) the measurement scheme is constructed
independent of the system’s dynamics, (ii) the protocol does not require
preparing multiple copies of quantum states at the same time, and (iii)
coherence between the system and ancilla states is required to be maintained
only during a short interaction time. These advantages make the proposed
scheme easier to be applied to any multi-time measurements with a longer time
scale. On the other hand, the ancilla state should be measured and used
repeatedly without affecting the system state. Such in-circuit detection (ICD)
and in-circuit initialization (ICI) of the ancilla are demanded for
experimental realization of the scheme.

In trapped-ion systems, ICD and ICI of the ancilla can be achieved by adopting
ion shuttling Kielpinski  _et al._ (2002); Wan  _et al._ (2019); Kaushal  _et
al._ (2020); Pino  _et al._ (2021) or multiple types of qubits including
multi-species systemsHome (2013); Tan  _et al._ (2015); Ballance  _et al._
(2015); Inlek  _et al._ (2017); Negnevitsky  _et al._ (2018); Bruzewicz  _et
al._ (2019); Wang  _et al._ (2022). We adopt a dual-species
${}^{171}\rm{Yb}^{+}$-${}^{138}\rm{Ba}^{+}$ trapped-ion system to reconstruct
two- and three-time quantum correlation functions and QPDs. The reconstructed
statistics follow quantum mechanical prediction with full contributions of
coherence. We also verify that the experimentally obtained QPDs have negative
or non-real values and their marginal distributions are preserved without
being disturbed by measurements. This is, to the best of our knowledge, the
first direct experimental realization of multi-time QPDs without process
tomography.

Quantum correlation function and QPDs.— A general form of the $N$-time quantum
correlation function Gardiner  _et al._ (2004) between observables
$\bm{A}=(A^{(1)},A^{(2)},\cdots,A^{(N)})$ at multiple times
$\bm{t}=(t_{1},\cdots,t_{N})$ can be expressed as

${\cal C}_{\bm{A}}(\bm{t}):=\mbox{Tr}[\rho A^{(1)}(t_{1})A^{(2)}(t_{2})\cdots
A^{(N)}(t_{N})].$ (1)

Here, we denote an observable $A$ after time $t$ in the Heisenberg picture,
$A(t)=U_{0\rightarrow t}^{\dagger}AU_{0\rightarrow t}$, where $U_{t\rightarrow
t^{\prime}}$ is a unitary operator describing the system’s time-evolution from
time $t$ to $t^{\prime}$. As $A^{(1)}(t_{1})A^{(2)}(t_{2})\cdots
A^{(N)}(t_{N})$ is not a Hermitian operator, the quantum correlation becomes,
in general, complex-valued, which can not be directly obtained from a single
observable in experiments.

By decomposing the observables into
$A^{(n)}=\sum_{i_{n}}a_{i_{n}}^{(n)}\Pi_{i_{n}}^{(n)}$ using the projection
operators $\Pi_{i_{n}}^{(n)}$, we also define a $N$-time QPDs as

${\cal
Q}_{\bm{i}}(\bm{t}):=\mbox{Tr}\left[\rho\Pi_{i_{1}}^{(1)}(t_{1})\Pi_{i_{2}}^{(2)}(t_{2})\cdots\Pi_{i_{N}}^{(N)}(t_{N})\right],$
(2)

where $\bm{i}=(i_{1},i_{2},\cdots,i_{N})$. For two-time points, such a
distribution is known as the Kirkwood-Dirac distribution Kirkwood (1933);
Dirac (1945), which has been widely adopted in the field of quantum
thermodynamics Allahverdyan (2014); Yunger Halpern (2017); Lostaglio (2018);
Levy and Lostaglio (2020) and quantum metrology Arvidsson-Shukur  _et al._
(2020); Lostaglio (2020) (See Ref. Lostaglio  _et al._ (2022) for a recent
overview). While being a complex-valued similarly to the quantum correlation
function, the QPDs preserves the marginal distribution of the statistics at
each time. We note that the QPDs contain more information than the correlation
function, as the latter can always be calculated from the former:

${\cal C}_{\bm{A}}(\bm{t})=\sum_{\bm{i}}a_{\bm{i}}{\cal Q}_{\bm{i}}(\bm{t}),$
(3)

where $a_{\bm{i}}=\prod_{n=1}^{N}a^{(n)}_{i_{n}}$.

As quantum coherence contributes to both ${\cal C}_{\bm{A}}(\bm{t})$ and
${\cal Q}_{\bm{i}}(\bm{t})$, these quantum statistics cannot be directly
derived from the classical distributions of projective measurement outcomes.
This is due to the fact that the projective measurement destroys all the off-
diagonal elements so that the state is disturbed after the measurement, i.e.,
$\sum_{i}\Pi_{i}\rho\Pi_{i}\neq\rho$ unless $[\rho,\Pi_{i}]=0$, which is often
referred to as measurement back-action.

Extracting quantum statics from ancilla-assisted measurements.—

Figure 1: (a) Diagram for sequential ancilla-assisted measurements. (b) By
adding more measurement bases on top of $\\{\left|0\right>,\left|1\right>\\}$
measurements, the quantum correlation function and the QPDs can be
reconstructed from the ancilla measurement outcomes.

Our main goal is to extract quantum correlation function and QPDs without
inducing measurement back-action. The main idea is finding a set of positive
operator-valued measurements (POVMs) $\\{M_{m}\\}$ satisfying

$\rho A=\sum_{m}\gamma_{m}(A)M_{m}\rho M_{m}^{\dagger},$ (4)

with complex coefficients $\gamma_{m}(A)$. To meet the completeness relation,
we also impose the condition $\sum_{m}M^{\dagger}_{m}M_{m}=\mathbb{1}$. In
principle, any POVM element $M_{m}$ can be realized by introducing an ancilla
state interacting with the system followed by the measurement of the ancilla
state, as described in Fig. 1(a). Our main observation is that the $N$-time
correlation function can be obtained as

${\cal
C}_{\bm{A}}(\bm{t})=\mathbb{E}\left[\gamma_{\bm{m}}(\bm{A})\right]=\sum_{\bm{m}}{\cal
P}_{\bm{m}}(\bm{t})\left[\prod_{n=1}^{N}\gamma_{m_{n}}(A^{(n)})\right],$ (5)

as an ensemble average over sequential outcome trajectories
$\bm{m}:=(m_{1},m_{2},\cdots,m_{N})$. The probability distribution of the
outcome $\bm{m}$ can explicitly be expressed as

$\displaystyle{\cal P}_{\bm{m}}(\bm{t})$ (6) $\displaystyle\quad:={\rm
Tr}\left[{\cal M}_{m_{N}}\circ\cdots\circ{\cal U}_{t_{2}\rightarrow
t_{3}}\circ{\cal M}_{m_{2}}\circ{\cal U}_{t_{1}\rightarrow t_{2}}\circ{\cal
M}_{m_{1}}(\rho)\right],$

where ${\cal U}_{t_{i}\rightarrow t_{j}}(\rho):=U_{t_{i}\rightarrow t_{j}}\rho
U_{t_{i}\rightarrow t_{j}}^{\dagger}$ and ${\cal M}_{m_{n}}:=M_{m_{n}}\rho
M_{m_{n}}^{\dagger}$. We note that ${\cal P}_{\bm{m}}(\bm{t})$ is an
experimentally accessible non-negative probability distribution, while the
quantum correlation function can be reconstructed by taking complex-valued
weights to it. As the measurement $\\{M_{m}\\}$ and the weights
$\gamma_{m}(A)$ at each time can be constructed independent of the system’s
dynamics, one can obtain joint distribution of various quantum dynamics
without changing the measurement setting. Following the same argument, the
QPDs can also be obtained from the sequential POVM outcomes as

${\cal
Q}_{\bm{i}}(\bm{t})=\mathbb{E}\left[\gamma_{\bm{m}}(\bm{\Pi}_{\bm{i}})\right],$
(7)

where
$\gamma_{\bm{m}}(\bm{\Pi}_{\bm{i}}):=\prod_{n=1}^{N}\gamma_{m_{n}}\left(\Pi^{(n)}_{i_{n}}\right)$.

Figure 2: Ancilla-assisted sequential measurement realized in a
${}^{171}\rm{Yb}^{+}$-${}^{138}\rm{Ba}^{+}$ trapped-ion system. The ancilla-
assisted measurement is performed as follows: (a) The system qubit involves
the evolution ${\cal U}_{t_{n-1}\rightarrow t_{n}}$, then (b) initialize the
${}^{138}\rm{Ba}^{+}$ qubit to $\left|0\right>$ by optical pumping (OPT). (c)
Perform a CNOT gate between two qubits through entangling operation. (d)
Measure the ancilla qubit with fluorescence detection (DET). The whole round
of measurement (b)–(d) is then completed.

While the proposed protocol extracts quantum statistics without measurement
back-action, we highlight that it does not violate the information
gain–disturbance trade-off relation Groenewold (1971); Fuchs and Peres (1996);
Ozawa (2003); Buscemi  _et al._ (2014); Fan  _et al._ (2015); Lee  _et al._
(2021); Hong  _et al._ (2022). This can be observed from the fact that
obtaining quantum statistics in Eqs. (5) and (7) requires more resources than
classical statistics, as a larger number of trajectories should be collected
to compensate the information disturbance caused by ancilla-assisted
measurements. More precisely, the expectation value
$\mathbb{E}\left[\gamma_{\bm{m}}(\bm{A})\right]$ fluctuates with
$|\gamma|_{\rm
max}:=\prod_{n=1}^{N}\left[\max_{m_{n}}\\{|\gamma_{m_{n}}(A^{(n)})|\\}\right]$
so that the number of trajectory for a desired precision scales with
$|\gamma|^{2}_{\rm max}$ from Hoeffding’s inequality Hoeffding (1963).

We also note that our approach is free from the restrictions of the
conventional approaches, which require either preparation of multiple copies
of quantum states Perarnau-Llobet  _et al._ (2017); Wu  _et al._ (2019) or
maintaining long time entanglement between the system and the ancilla Somma
_et al._ (2002); Souza  _et al._ (2011); Pedernales  _et al._ (2014); Xin  _et
al._ (2017).

Explicit protocol for a two-level system.— For an given observable $A$,
various choices of $\\{M_{m},\gamma_{m}\\}$ can be taken to satisfy Eq. (4).
In order to focus on the statistics of a two-level system, we consider a POVM
of the following form:

$M_{m}=\frac{\langle\phi_{m}|0\rangle}{\sqrt{\alpha}}\left|0\right>\langle
0|+\frac{\langle\phi_{m}|1\rangle}{\sqrt{\alpha}}\left|1\right>\langle 1|,$
(8)

where $\sum_{m}\left|\phi_{m}\right>\langle\phi_{m}|=\alpha\mathbb{1}$ and
$\left|\phi_{m}\right>$ are not necessarily orthogonal to each other. Such a
POVM can be implemented using ancilla-assisted measurement as follows: (i)
prepare an ancilla state at $\left|0\right>$, (ii) interact the system and the
ancilla via the CNOT-gate, (iii) measure the ancilla onto states
$\left|\phi_{m}\right>$ (see Fig. 1). After the measurement, the ancilla qubit
is re-initialized to $\left|0\right>$ so that it can be used for the
measurement at the next time. We note that this scheme readily incorporates
the classical projection measurements by performing the $z$-basis measurement
of the ancilla state, i.e.,
$\\{\left|\phi_{m}\right>\\}=\\{\left|0\right>,\left|1\right>\\}$.

Quantum statistics can be obtained by taking more measurements on top of the
$z$-basis measurement. For example, the Pauli-$Z$ operator,
$Z=\left|0\right>\langle 0|-\left|1\right>\langle 1|$, acting on the right
side of the quantum state, i.e., $\rho Z$ can be realized by a set of
measurements
$\\{\left|\phi_{m}\right>\\}=\\{\left|0\right>,\left|1\right>,\left|+_{y}\right>,\left|-_{y}\right>\\}$,
where $\left|\pm_{y}\right>=(\left|0\right>\pm i\left|1\right>)/\sqrt{2}$ (see
Fig. 1(b)). In this case, the weight can be calculated as
$\\{\gamma_{m}(Z)\\}=\\{2,-2,2i,-2i\\}$. On the other hand, such a measurement
set is not enough to obtain QPDs described in Eq. (7). This can be intuitively
understood by the fact that QPDs contain more information than correlation
function as seen in Eq. (3). We find that by adding $x$-basis measurements
onto $\left|\pm_{x}\right>=(\left|0\right>\pm\left|1\right>)/\sqrt{2}$, which
compose a set of measurements
$\\{\left|\phi_{m}\right>\\}=\\{\left|0\right>,\left|1\right>,\left|+_{y}\right>,\left|-_{y}\right>,\left|+_{x}\right>,\left|-_{x}\right>\\}$
(see Fig. 1(c)), is sufficient to obtain QPDs with various choices of
$\gamma_{m}(\Pi_{0})$ and $\gamma_{m}(\Pi_{1})$ (see Supplemental Material Sup
for details and further generalizations).

Experimental realization.— We demonstrate the proposed protocol with two
different species of ${}^{171}\rm{Yb}^{+}$ and ${}^{138}\rm{Ba}^{+}$ ions
trapped in a single trap Wang  _et al._ (2017, 2021, 2022). A crucial part of
the protocol is repeated times of ancilla measurement and initialization
without disturbing the system. The dual-species trapped-ion system provides a
promising direction for this as each ion is controlled by lasers with
different wavelengths. With minimal influence on each other, ICD and ICI on
two different ions can independently performed Wang  _et al._ (2022). In the
experiment, two hyperfine levels in the $S_{1/2}$ of the ${}^{171}\rm{Yb}^{+}$
ion, $\left|F=0,m_{F}=0\right>\equiv\left|0\right>_{\mathrm{Yb}}$ and
$\left|F=1,m_{F}=0\right>\equiv\left|1\right>_{\mathrm{Yb}}$ with the energy
splitting of 12.6428 GHz, serves as the system qubit. The qubit is insensitive
to environmental noise and demonstrated a long coherence time. For the
${}^{138}\rm{Ba}^{+}$ ion, two Zeeman levels in the $S_{1/2}$ manifold are
encode as the ancilla qubit,
$\left|s_{j}=1/2\right>\equiv\left|0\right>_{\mathrm{Ba}}$ and
$\left|s_{j}=-1/2\right>\equiv\left|1\right>_{\mathrm{Ba}}$ with the energy
splitting of 16.2 MHz. Raman transitions are used to manipulate the
${}^{171}\rm{Yb}^{+}$ and ${}^{138}\rm{Ba}^{+}$ qubits with 355 nm and 532 nm
lasers, respectively, as shown in Fig. 2. For the entangling operations for
${}^{171}\rm{Yb}^{+}$ and ${}^{138}\rm{Ba}^{+}$ qubits, we simultaneously
apply the 355 nm and 532 nm laser beams with properly chosen frequencies.

(a)

(c)

(e)

(b)

(d)

(f)

Figure 3: (a) Two- and (b) three-time correlation functions reconstructed from
the ancilla-assisted measurements of the
${}^{171}\rm{Yb}^{+}$-${}^{138}\rm{Ba}^{+}$ system. The real (blue) and
imaginary (red) parts of the correlation functions obtained in the experiments
are represented in dots, while solid and dashed curves refer to imaginary and
real parts of the theoretical values, respectively. (c) Two- and (d) three-
time QPDs reconstructed from the outcome trajectories for $\theta=0.3\pi$ and
$\theta=0.74\pi$, respectively. The horizontal index
$\bm{i}=(i_{1},i_{2},\cdots)$ indicates the projection operator $\Pi_{i_{n}}$
at time $t_{n}$ described in Eq. (2). For the three-time QPDs, negativity is
verified for (0,1,0) and (1,1,0), pointed in red-dashed-line boxes. (e) The
marginal distributions of the $z$-basis measurement at time $t_{2}$. Dots are
reconstructed from the experimental data where three-time ancilla-assisted
measurements are performed. Solid and dashed lines refer to the exact
$z-$population when no measurement is performed before $t_{2}$. dash-dotted
line refers to the $z-$population when projection measurements are performed
at time $t_{1}$ or $t_{2}$, which is balanced for every $\theta$. (f) Leggett-
Garg inequality test with three-time QPDs (blue) and projection measurements
(gray). Dots are experimental data. Solid curves are theoretical prediction.
At the point $\theta=0.74\pi$ from three-time QPDs, we obtain $K\approx 1.17$,
which violates the classical bound $1$. For all the data, the dynamics are
given as $U_{t_{1}\rightarrow t_{2}}=R_{x}(\theta)$ and $U_{t_{2}\rightarrow
t_{3}}=R_{y}(\theta^{2})$, and error bars represent the standard error of the
mean (SEM). For each measurement configuration, $100$ outcome trajectories are
collected.

We first prepare the system qubit to
$\left|+\right>_{\mathrm{Yb}}=\frac{1}{\sqrt{2}}(\left|0\right>_{\mathrm{Yb}}+\left|1\right>_{\mathrm{Yb}})$
and the ancilla qubit to $\left|0\right>_{\mathrm{Ba}}$ (at $t_{1}$ in Fig.
1(a)). The initialization to $\left|0\right>_{\mathrm{Yb(Ba)}}$ states is
performed by the standard optical pumping methods for both
${}^{171}\rm{Yb}^{+}$ and ${}^{138}\rm{Ba}^{+}$ ions Wang  _et al._ (2022). We
then use 355 nm laser beams to prepare $\left|+\right>_{\mathrm{Yb}}$ state.
The system state then evolves under a single-qubit $x$-axis rotation
$U_{t_{1}\rightarrow t_{2}}=R_{x}(\theta)=R(\theta,0)$ (from $t_{1}$ to
$t_{2}$) followed by $y$-axis rotation $U_{t_{2}\rightarrow
t_{3}}=R_{y}(\theta^{2})=R(\theta^{2},-\pi/2)$ with $\theta\in[0,\pi]$ (from
$t_{2}$ to $t_{3}$), both of which are performed by applying 355 nm laser
beams. Here a general single-qubit rotation is defined as
$R(\theta,\phi)=\left(\begin{array}[]{cc}\cos(\frac{\theta}{2})&-ie^{-i\phi}\sin(\frac{\theta}{2})\\\
-ie^{i\phi}\sin(\frac{\theta}{2})&\cos(\frac{\theta}{2})\end{array}\right)$.
At each time, the ancilla-assisted measurement is realized by the CNOT gate
between the system (${}^{171}\rm{Yb}^{+}$) and ancilla (${}^{138}\rm{Ba}^{+}$)
states followed by the projection measurement of the ancilla. The Mølmer-
Sørensen (M-S) gate Sørensen and Mølmer (1999) is used to construct the CNOT
gate. The average fidelity of the M-S gate is $94\pm 2\%$ in the experiment.
Details of the M-S gate and its imperfection are described in Supplemental
Material Sup .

The measurement of the ${}^{138}\rm{Ba}^{+}$ ion is realized by applying
detection laser after shelving $\left|0\right>_{\mathrm{Ba}}$ state to
$D_{5/2}$ manifold with the help of 1762 nm narrow linewidth laser, where
detection fidelity is $98.9\%$. In order to simplify the experimental steps,
the final measurement is replaced by the projection measurement on the
${}^{171}\rm{Yb}^{+}$ ion as we do not need to worry about the system
disturbance. The detection fidelity of ${}^{171}\rm{Yb}^{+}$ ion is $98.4\%$.
The detection infidelities are improved by using the detection-error correct
method introduced in Ref. Shen and Duan (2012) for both ${}^{171}\rm{Yb}^{+}$
and ${}^{138}\rm{Ba}^{+}$ ions.

We reconstruct correlation functions from the experimental measurement
results. Fig. 3(a) and (b) show the two- and three-time correlation functions
${\cal C}_{ZZ}(t_{1},t_{2})$ and ${\cal C}_{ZZZ}(t_{1},t_{2},t_{3})$
reconstructed from the experimental measurements. The two-time correlation
function is in a good agreement with the theoretical predictions, while the
three-time correlation function has relatively large deviations. This is
mainly because of that adding more measurement steps induces a larger
fluctuation of $\gamma_{\bm{m}}=\prod_{n=1}^{N}\gamma_{m_{n}}$ for each
realization. Another reason is that durations for the ancilla measurements are
relatively long to collect enough fluorescence, which introduces large phase
fluctuations and degrades the performance of the CNOT gates.

Fig. 3(c) and (d) show the two- and three-time QPDs ${\cal
Q}_{\bm{i}}(t_{1},t_{2})$ and ${\cal Q}_{\bm{i}}(t_{1},t_{2},t_{3})$ obtained
by collecting additional trajectories including the $x$-basis measurement at
$t_{1}$ and $t_{2}$. To our knowledge, it is the first experimental
reconstruction of the QPDs beyond two-time points. The negativity in real
parts as well as non-vanishing imaginary parts are clearly observed in the
experimentally reconstructed three-time QPDs (see Fig. 3(d)). The contribution
of coherence can be clearly seen by comparing the marginal distribution at
time $t_{3}$ of the reconstructed QPDs with the projection measurement case.
When the projective measurement is performed to the initial state
$\rho=\left|+\right>\langle+|$, the state after the measurement becomes the
maximally mixed state without any coherence,
$\rho^{\prime}=\sum_{i}\Pi_{i}\rho\Pi_{i}=\mathbb{1}/2$, so that the
population is always balanced at $t_{3}$ regardless of the unitary dynamics
$U_{t_{1}\rightarrow t_{2}}$ and $U_{t_{2}\rightarrow t_{3}}$. Fig. 3(e) shows
that the marginal distribution at $t_{3}$ is not always balanced but varies by
taking different unitary dynamics, following the quantum mechanical prediction
without measurement back-action.

While the negative and non-real values can already be interpreted as a
nonclassical feature of QPDs, a stronger notion of nonclassicality can be
witnessed by Leggett-Garg inequality violation Leggett and Garg (1985). For
any classical joint probability distribution $P(i_{1},i_{2},i_{3})$ and a
bounded function $E_{t_{n}}(i_{n})\in\\{1,-1\\}$ at each time $t_{n}$, the
following inequality should hold:

$K=\langle E_{t_{1}}E_{t_{2}}\rangle+\langle E_{t_{2}}E_{t_{3}}\rangle-\langle
E_{t_{1}}E_{t_{3}}\rangle\leq 1,$ (9)

which was introduced by Leggett and Garg Leggett and Garg (1985) to test
microscopic realism. By taking $E_{t_{1}}(i_{1})=(-1)^{i_{1}}$,
$E_{t_{2}}(i_{2})=1$, and $E_{t_{3}}(i_{3})=(-1)^{i_{3}}$, we observe the
violation of Eq. (9) with $K\approx 1.17$ from the real part of the
experimentally reconstructed three-time QPDs (see Fig. 3(f)). This strongly
supports that the reconstructed QPDs can be sharply distinguished from
classical distributions.

Remarks.— We have proposed and experimentally demonstrated a measurement
scheme to overcome the measurement back-action in obtaining multi-time quantum
statistics. Using the ${}^{171}\rm{Yb}^{+}$-${}^{138}\rm{Ba}^{+}$ trapped-ion
system, we have experimentally demonstrated that two- and three-time quantum
correlation functions and QPDs reconstructed from the ancilla-assisted
measurements follow the quantum mechanical prediction with full contributions
of coherence.

We highlight that our method is not limited to the two-level system and can be
generalized to any $d$-dimensional quantum system by replacing the CNOT gate
with the CSUM gate, $U_{\rm CSUM}=\sum_{i,j=0}^{d-1}\left|i,i\oplus
j\right>\langle i,j|$ and taking informationally complete basis measurements
on a $d$-dimensional ancilla state Sup . Furthermore, it is straightforward to
extend the protocol to be applied to non-unitary dynamics, by simply replacing
${\cal U}_{t_{i}\rightarrow t_{j}}$ in Eq. (6) to any completely positive
trace-preserving quantum channel. As the ancilla-assisted measurement protocol
is independent of the system’s dynamics and only requires a single copy of the
system state with a short time system-ancilla interaction, our approach can be
a useful experimental tool for exploring quantum statistics of both open and
closed quantum systems.

Note added.– While preparing the manuscript, we recognized a recent work Del
Re  _et al._ (2022) which deals with the multi-time correlation function.

###### Acknowledgements.

This work was supported by the National Key Research and Development Program
of China under Grants No.2016YFA0301900 and No.2016YFA0301901, the National
Natural Science Foundation of China Grants No.92065205, and No.11974200. H.K.
is supported by the KIAS Individual Grant No.CG085301 at the Korea Institute
for Advanced Study. MK acknowledges the KIST Open Research program.

## References

  * Wigner (1932) E. Wigner, Phys. Rev. 40, 749 (1932).
  * Margenau and Hill (1961) H. Margenau and R. N. Hill, Prog. Theor. Exp. Phys. 26, 722 (1961).
  * Kirkwood (1933) J. G. Kirkwood, Phys. Rev. 44, 31 (1933).
  * Dirac (1945) P. A. M. Dirac, Rev. Mod. Phys. 17, 195 (1945).
  * BELL (1966) J. S. BELL, Rev. Mod. Phys. 38, 447 (1966).
  * Leggett and Garg (1985) A. J. Leggett and A. Garg, Phys. Rev. Lett. 54, 857 (1985).
  * Feynman (1987) R. P. Feynman, Quantum implications: essays in honour of David Bohm , 235 (1987).
  * Allahverdyan (2014) A. E. Allahverdyan, Phys. Rev. E 90, 032137 (2014).
  * Perarnau-Llobet  _et al._ (2017) M. Perarnau-Llobet, E. Bäumer, K. V. Hovhannisyan, M. Huber, and A. Acin, Phys. Rev. Lett. 118, 070601 (2017).
  * Yunger Halpern (2017) N. Yunger Halpern, Phys. Rev. A 95, 012120 (2017).
  * Lostaglio (2018) M. Lostaglio, Phys. Rev. Lett. 120, 040602 (2018).
  * Levy and Lostaglio (2020) A. Levy and M. Lostaglio, PRX Quantum 1, 010309 (2020).
  * Kwon and Kim (2019) H. Kwon and M. S. Kim, Phys. Rev. X 9, 031029 (2019).
  * Veitch  _et al._ (2012) V. Veitch, C. Ferrie, D. Gross, and J. Emerson, New J. Phys 14, 113011 (2012).
  * Mari and Eisert (2012) A. Mari and J. Eisert, Phys. Rev. Lett. 109, 230503 (2012).
  * Howard  _et al._ (2014) M. Howard, J. Wallman, V. Veitch, and J. Emerson, Nature (London) 510, 351 (2014).
  * Pashayan  _et al._ (2015) H. Pashayan, J. J. Wallman, and S. D. Bartlett, Phys. Rev. Lett. 115, 070501 (2015).
  * Rahimi-Keshari  _et al._ (2016) S. Rahimi-Keshari, T. C. Ralph, and C. M. Caves, Phys. Rev. X 6, 021039 (2016).
  * Kwon  _et al._ (2019) H. Kwon, K. C. Tan, T. Volkoff, and H. Jeong, Phys. Rev. Lett. 122, 040503 (2019).
  * Arvidsson-Shukur  _et al._ (2020) D. R. Arvidsson-Shukur, N. Yunger Halpern, H. V. Lepage, A. A. Lasek, C. H. Barnes, and S. Lloyd, Nat. Commun. 11, 1 (2020).
  * Lostaglio (2020) M. Lostaglio, Phys. Rev. Lett. 125, 230603 (2020).
  * Kubo (1957) R. Kubo, J. Phys. Soc. Japan 12, 570 (1957).
  * Mandel and Wolf (1995) L. Mandel and E. Wolf,  _Optical coherence and quantum optics_ (Cambridge University Press, Cambridge, England, 1995).
  * Gardiner  _et al._ (2004) C. Gardiner, P. Zoller, and P. Zoller,  _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_ (Springer, Berlin, 2004).
  * Clerk  _et al._ (2010) A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Rev. Mod. Phys. 82, 1155 (2010).
  * Krumm  _et al._ (2016) F. Krumm, J. Sperling, and W. Vogel, Phys. Rev. A 93, 063843 (2016).
  * Yunger Halpern  _et al._ (2018) N. Yunger Halpern, B. Swingle, and J. Dressel, Phys. Rev. A 97, 042105 (2018).
  * Dowling  _et al._ (2021) N. Dowling, P. Figueroa-Romero, F. A. Pollock, P. Strasberg, and K. Modi, arXiv:2108.07420 (2021).
  * Groenewold (1971) H. J. Groenewold, Int. J. Theor. Phys. 4, 327 (1971).
  * Fuchs and Peres (1996) C. A. Fuchs and A. Peres, Phys. Rev. A 53, 2038 (1996).
  * Ozawa (2003) M. Ozawa, Phys. Rev. A 67, 042105 (2003).
  * Buscemi  _et al._ (2014) F. Buscemi, M. J. W. Hall, M. Ozawa, and M. M. Wilde, Phys. Rev. Lett. 112, 050401 (2014).
  * Fan  _et al._ (2015) L. Fan, W. Ge, H. Nha, and M. S. Zubairy, Phys. Rev. A 92, 022114 (2015).
  * Lee  _et al._ (2021) S.-W. Lee, J. Kim, and H. Nha, Quantum 5, 414 (2021).
  * Hong  _et al._ (2022) S. Hong, Y.-S. Kim, Y.-W. Cho, J. Kim, S.-W. Lee, and H.-T. Lim, Phys. Rev. Lett. 128, 050401 (2022).
  * Buhrman  _et al._ (2001) H. Buhrman, R. Cleve, J. Watrous, and R. de Wolf, Phys. Rev. Lett. 87, 167902 (2001).
  * Somma  _et al._ (2002) R. Somma, G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, Phys. Rev. A 65, 042323 (2002).
  * Johansen (2007) L. M. Johansen, Phys. Rev. A 76, 012119 (2007).
  * Buscemi  _et al._ (2013) F. Buscemi, M. Dall’Arno, M. Ozawa, and V. Vedral, arXiv:1312.4240 (2013).
  * Pedernales  _et al._ (2014) J. S. Pedernales, R. Di Candia, I. L. Egusquiza, J. Casanova, and E. Solano, Phys. Rev. Lett. 113, 020505 (2014).
  * Souza  _et al._ (2011) A. Souza, I. Oliveira, and R. Sarthour, New J. Phys. 13, 053023 (2011).
  * Piacentini  _et al._ (2016) F. Piacentini, A. Avella, M. P. Levi, M. Gramegna, G. Brida, I. P. Degiovanni, E. Cohen, R. Lussana, F. Villa, A. Tosi, F. Zappa, and M. Genovese, Phys. Rev. Lett. 117, 170402 (2016).
  * Xin  _et al._ (2017) T. Xin, J. S. Pedernales, L. Lamata, E. Solano, and G.-L. Long, Sci. Rep. 7, 1 (2017).
  * Ringbauer  _et al._ (2018) M. Ringbauer, F. Costa, M. E. Goggin, A. G. White, and A. Fedrizzi, npj Quantum Inf. 4, 1 (2018).
  * Wu  _et al._ (2019) K.-D. Wu, Y. Yuan, G.-Y. Xiang, C.-F. Li, G.-C. Guo, and M. Perarnau-Llobet, Sci. Adv. 5, eaav4944 (2019).
  * Del Re  _et al._ (2022) L. Del Re, B. Rost, M. Foss-Feig, A. Kemper, and J. Freericks, arXiv:2204.12400 (2022).
  * Kielpinski  _et al._ (2002) D. Kielpinski, C. Monroe, and D. J. Wineland, Nature (London) 417, 709 (2002).
  * Wan  _et al._ (2019) Y. Wan, D. Kienzler, S. D. Erickson, K. H. Mayer, T. R. Tan, J. J. Wu, H. M. Vasconcelos, S. Glancy, E. Knill, D. J. Wineland, _et al._ , Science 364, 875 (2019).
  * Kaushal  _et al._ (2020) V. Kaushal, B. Lekitsch, A. Stahl, J. Hilder, D. Pijn, C. Schmiegelow, A. Bermudez, M. Müller, F. Schmidt-Kaler, and U. Poschinger, AVS Quantum Sci. 2, 014101 (2020).
  * Pino  _et al._ (2021) J. M. Pino, J. M. Dreiling, C. Figgatt, J. P. Gaebler, S. A. Moses, M. Allman, C. Baldwin, M. Foss-Feig, D. Hayes, K. Mayer, _et al._ , Nature (London) 592, 209 (2021).
  * Home (2013) J. P. Home, Adv. At. Mol. Opt. Phys. 62, 231 (2013).
  * Tan  _et al._ (2015) T. R. Tan, J. P. Gaebler, Y. Lin, Y. Wan, R. Bowler, D. Leibfried, and D. J. Wineland, Nature (London) 528, 380 (2015).
  * Ballance  _et al._ (2015) C. Ballance, V. Schäfer, J. P. Home, D. Szwer, S. C. Webster, D. Allcock, N. M. Linke, T. Harty, D. Aude Craik, D. N. Stacey, _et al._ , Nature (London) 528, 384 (2015).
  * Inlek  _et al._ (2017) I. V. Inlek, C. Crocker, M. Lichtman, K. Sosnova, and C. Monroe, Phys. Rev. Lett. 118, 250502 (2017).
  * Negnevitsky  _et al._ (2018) V. Negnevitsky, M. Marinelli, K. K. Mehta, H.-Y. Lo, C. Flühmann, and J. P. Home, Nature (London) 563, 527 (2018).
  * Bruzewicz  _et al._ (2019) C. Bruzewicz, R. McConnell, J. Stuart, J. Sage, and J. Chiaverini, npj Quantum Inf. 5, 1 (2019).
  * Wang  _et al._ (2022) P. Wang, J. Zhang, C.-Y. Luan, M. Um, Y. Wang, M. Qiao, T. Xie, J.-N. Zhang, A. Cabello, and K. Kim, Sci. Adv. 8, eabk1660 (2022).
  * Lostaglio  _et al._ (2022) M. Lostaglio, A. Belenchia, A. Levy, S. Hernández-Gómez, N. Fabbri, and S. Gherardini, arXiv:2206.11783 (2022).
  * Hoeffding (1963) W. Hoeffding, J. Am. Stat. Assoc. 58, 13 (1963).
  * (60) See Supplemental Material at [url] for more details on the theory of demonstration of multi-time quantum statistics without measurement back-action, which includes Ref. [57] .
  * Wang  _et al._ (2017) Y. Wang, M. Um, J. Zhang, S. An, M. Lyu, J.-N. Zhang, L.-M. Duan, D. Yum, and K. Kim, Nat. Photonics 11, 646 (2017).
  * Wang  _et al._ (2021) P. Wang, C.-Y. Luan, M. Qiao, M. Um, J. Zhang, Y. Wang, X. Yuan, M. Gu, J. Zhang, and K. Kim, Nat. Commun. 12, 1 (2021).
  * Sørensen and Mølmer (1999) A. Sørensen and K. Mølmer, Phys. Rev. Lett. 82, 1971 (1999).
  * Shen and Duan (2012) C. Shen and L. Duan, New J. Phys 14, 053053 (2012).

Supplemental Material for “Demonstration of multi-time quantum statistics
without measurement back-action”

##  I Experimental Methods

###  I.1 Experimental circuit

Multi-time measurement circuit is shown in Fig. 4(a). The system qubit
information is measured with the assistance of ancilla qubit. The ancilla
qubit is connected to the system qubit via a CNOT gate consisting of a Mølmer-
Sørensen (M-S) gate and four single qubit rotations, as shown in Fig. 4(b).
Before applying the M-S gate, the average phonon number in the axial out-of-
phase (in-phase) mode is cooled to 0.04 (0.11) by three steps: Doppler
cooling, electromagnetically-induced-transparency (EIT) cooling, and Raman
sideband cooling. State fidelity of the Bell state created by the M-S gates is
$0.94\pm 0.02$, which is limited by parameter shift in long-sequence
operations and the axial in-phase mode cooling imperfection.

Figure 4: Experimental circuit for multi-time measurement. (a) The system and
ancilla qubits are first prepared to $\rho$ and $\left|0\right>$,
respectively. The first-time measurement starts with a CNOT gate between the
system and the ancilla, followed by a measurement of the ancilla qubit. After
a unitary evolution ${\cal U}_{t_{1}\rightarrow t_{2}}$ on the system qubit, a
second-time measurement is performed with a similarly process to the first
one. This scheme can be directly extended to multi-time measurement. Red and
green boxes indicate the two-time and three-time measurement, respectively.
(b) The CNOT gate is composed of a M-S gate and four-single qubit rotations.
The M-S gate can be described as $\exp(i\frac{\pi}{4}X\otimes X)$, where $X$
is the Pauli operator. The single-qubit rotation is defined as
$R(\theta,\phi)=\left(\begin{array}[]{cc}\cos(\frac{\theta}{2})&-ie^{-i\phi}\sin(\frac{\theta}{2})\\\
-ie^{i\phi}\sin(\frac{\theta}{2})&\cos(\frac{\theta}{2})\end{array}\right)$.
(c) The last-time measurement can be performed by direct measurement in basis
$m_{2}\in\\{\left|0\right>_{z},\left|1\right>_{z}\\}$ on the system qubit
since coherence is no more needed for further measurement.

###  I.2 Experimental setup for ancilla-assisted measurement

As shown in Fig. 5(a), a dual-species trapped-ion system, which traps one
${}^{171}\rm{Yb}^{+}$ ion and one ${}^{138}\rm{Ba}^{+}$ ion in a four-rod
trap, is used to realize the ancilla-assisted measurements.
${}^{171}\rm{Yb}^{+}$ and ${}^{138}\rm{Ba}^{+}$ ions serve as the system and
the ancilla respectively. The two ions have different energy structures, and
require different initialization and detection lasers. Therefore, the
operation on one ion does not affect the other ion. It is thus possible to
perform in-circuit detection (ICD) and initialization (ICI), where one qubit
is detected or initialized without affecting the other qubit.Wang  _et al._
(2022). This special property is essential for the ancilla-assisted
measurement. Experimentally, we also need single-qubit rotations and two-qubit
M-S gate. As shown in Fig. 5, we use 355 nm and 532 nm Raman transitions to
realize gate operations of ${}^{171}\rm{Yb}^{+}$ and ${}^{138}\rm{Ba}^{+}$
ions, respectively.

Figure 5: Experimental setup. (a) ${}^{171}\rm{Yb}^{+}$ and
${}^{138}\rm{Ba}^{+}$ ions trapped in a four-rod trap. The magnetic field of
5.8 Gauss is applied along the $x$-axis. In order to realize M-S gate between
two ions, the ${}^{171}\rm{Yb}^{+}$ ion is controlled by illuminating a pair
of 355 nm lasers (purple) with a frequency of $f_{355}$ and
$f_{355}+f_{\mathrm{Yb}}\pm(\delta+f_{z})$, where $f_{\mathrm{Yb}}$ is the
qubit splitting of ${}^{171}\rm{Yb}^{+}$ qubit. The ${}^{138}\rm{Ba}^{+}$ ion
is controlled by illuminating a pair of 532 nm lasers (green) with a frequency
of $f_{532}$ and $f_{532}+f_{\mathrm{Ba}}\pm(\delta+f_{z})$, where
$f_{\mathrm{Ba}}$ is the qubit splitting of ${}^{138}\rm{Ba}^{+}$ qubit,
$f_{z}$ is the axial OOP (out-of-phase) mode frequency, and $\delta$ is the
laser detuning from the sideband of OOP mode. Raman laser directions are
represented by thick arrows, and polarization is represented by thin arrows
and dots. (b) Energy level and related lasers for ${}^{171}\rm{Yb}^{+}$ and
${}^{138}\rm{Ba}^{+}$ ions. The ${}^{171}\rm{Yb}^{+}$ ion’s hyperfine qubit
and ${}^{138}\rm{Ba}^{+}$ ion’s Zeeman qubit serve as the system qubit and the
ancilla qubit, respectively. Single-qubit rotations are realized by resonant
Raman transitions. Two-qubit entangling operation is realized by applying
bichromatic Raman laser for both ions. The 1762 nm laser is the shelving laser
used for ${}^{138}\rm{Ba}^{+}$ qubit detection.

###  I.3 Post-selection detection

Fig. 6 shows the detailed process of the ${}^{138}\rm{Ba}^{+}$ ion
fluorescence detection. In this detection process, we first use 1762 nm laser
to shelve $\left|0\right>_{\mathrm{Ba}}$ population to $\mathrm{D}_{5/2}$, and
then apply 493 nm laser to drive the transition between $\mathrm{S}_{1/2}$ and
$\mathrm{P}_{1/2}$. If the qubit state is $\left|1\right>_{\mathrm{Ba}}$, the
shelving operation has no effect and the later 493 nm laser will produce a
large number of photons. But if the qubit state is
$\left|0\right>_{\mathrm{Ba}}$, the population is shelved to
$\mathrm{D}_{5/2}$ and the 493 nm laser does not produce any photons. Then
qubit states can be distinguished by the number of photons scattered.

In our multi-time measurements circuit, the ancilla (${}^{138}\rm{Ba}^{+}$ )
needs to be repeatedly detected and used. However, as shown in Fig. 6(a), the
detecting of bright state produces a large number of photons, thereby heats
the ion chain and further degrades the performance of subsequent CNOT
operations. To solve this problem, we adopt a post-selection approach for all
ICD, which use only dark state data. For bright state data, we transfer the
bright state to dark state by $\pi$ flip pulse and then measure.

Figure 6: Fluorescence detection of ${}^{138}\rm{Ba}^{+}$ ion. The detection
process involves two steps: first, use a 1762 nm laser to shelve the
$\left|0\right>_{\mathrm{Ba}}$ population to $\mathrm{D}_{5/2}$, and then use
a 493 nm laser to drive the transition between $\mathrm{S}_{1/2}$ and
$\mathrm{P}_{1/2}$. (a) Bright state detection. when the ancilla is at state
$\left|1\right>_{\mathrm{Ba}}$, the shelving operation has no effect on the
qubit. The later 493 nm laser produces a large number of photons. (b) Dark
state detection. When the ancilla is at state $\left|0\right>_{\mathrm{Ba}}$,
its population will be shelved to $\mathrm{D}_{5/2}$. Therefore, the 493 nm
laser does not produce any photons.

##  II Obtaining quantum statistics from POVM measurements

###  II.1 General conditions

For a $d$-dimensional quantum system. Let us consider a slightly more general
scenario such that operators are acting both on left and right-side of the
density matrix:

$B\rho A=\sum_{m}\gamma_{m}(B,A)M_{m}\rho M_{m}^{\dagger},$ (10)

where the case discussed in the main text is a special case by taking $B=I$.
We also generalize the POVM measurements by introducing an ancilla state
having the same dimension as the system, initially prepared at
$\left|0\right>_{R}$. After applying the CSUM gate, a generalized CNOT gate,
$U_{\rm CSUM}=\sum_{i,j=0}^{d-1}\left|i,i\oplus j\right>\langle i,j|$ followed
by the ancilla measurement onto the states $\\{\left|\phi_{m}\right>_{R}\\}$
satisfying $\sum_{m}\left|\phi_{m}\right>\langle\phi_{m}|=\alpha\mathbb{1}$,
we can construct each POVM element as,

$M_{m}=\frac{\langle\phi_{m}|U_{\rm
CSUM}\left|0\right>_{R}}{\sqrt{\alpha}}=\sum_{i=0}^{d-1}\left(\frac{\langle\phi_{m}|i\rangle_{R}}{\sqrt{\alpha}}\right)\left|i\right>\langle
i|.$ (11)

We then find a sufficient condition for a set of measurements
$\\{\left|\phi_{m}\right>\\}$ to reconstruct the operation in Eq. (10) for
diagonal operators $A$ and $B$ as follows:

######  Proposition 1 (Sufficient condition of a measurement set).

For operators $A=\sum_{i}a_{i}\left|i\right>\langle i|$ and
$B=\sum_{i}b_{i}\left|i\right>\langle i|$ that are diagonal in the basis
$\left|i\right>$, there always exists $\gamma_{m}(B,A)$ satisfying Eq. (10)
for a measurement set $\\{\left|\phi_{m}\right>\\}$ which is informationally
complete.

###### Proof.

By definition of a set of informationally complete states, one can express any
operator $O$ (no matter Hermitian or not) as
$O=\sum_{m}c_{m}(O)\left|\phi_{m}\right>\langle\phi_{m}|$ with some complex
coefficients $c_{m}(O)$. Now let’s take an operator $O(B,A)$ as follows:

$O(B,A)=\sum_{i,j=0}^{d-1}a_{i}b_{j}\left|i\right>\langle
j|=\sum_{m}c_{m}(O(B,A))\left|\phi_{m}\right>\langle\phi_{m}|.$ (12)

We then observe that taking $\gamma_{m}(B,A)=\alpha c_{m}(O(B,A))$ leads to

$\displaystyle\sum_{m}\gamma_{m}(B,A)M_{m}\rho M_{m}^{\dagger}$ (13)
$\displaystyle=\sum_{m}\sum_{i,j=0}^{d-1}c_{m}(O(B,A))\langle
i|\phi_{m}\rangle\langle\phi_{m}|j\rangle\left|j\right>\langle
j|\rho\left|i\right>\langle i|$
$\displaystyle=\sum_{i,j=0}^{d-1}a_{i}b_{j}\left|j\right>\langle
j|\rho\left|i\right>\langle i|$ $\displaystyle=B\rho A,$

which completes the proof. ∎

As an informationally complete basis has at least $d^{2}$ elements for a
$d$-dimensional system, we note that obtaining the quantum statistics in the
form of Eq. (10) always requires more measurement setting than the projection
measurements $\\{\Pi_{i}\\}_{i=0}^{d-1}$. However, there could be an
informationally incomplete measurement set even with $d^{2}$ elements. For
example, for the measurement sets discussed in the main text,
$\\{\left|\phi_{m}\right>\\}=\\{\left|0\right>,\left|1\right>,\left|+_{y}\right>,\left|-_{y}\right>\\}$
is not informationally complete, while
$\\{\left|\phi_{m}\right>\\}=\\{\left|0\right>,\left|1\right>,\left|+_{y}\right>,\left|-_{y}\right>,\left|+_{x}\right>,\left|-_{x}\right>\\}$
is informationally complete. Nevertheless, such the measurement set based on
the Pauli basis is easier to be realized in experiments.

###  II.2 Optimizing the weights

When a set of POVM $\\{M_{m}\\}$ to obtain the statistics of hermitian
operators $A$ and $B$ is overcomplete, there could be a various choice of
weight vector $\gamma_{m}(B,A)$ satisfying Eq. (10). As the expectation value
$\mathbb{E}\left[\gamma_{m}(B,A)\right]$ fluctuates with $|\gamma(B,A)|_{\rm
max}:=\max_{m}\\{|\gamma_{m}(B,A)|\\}$, the optimal choice would be to
minimize $|\gamma|_{\rm max}$. Such a optimization problem can be expressed
as:

$\displaystyle{\rm for~{}given:}$ $\displaystyle\qquad A,~{}B,~{}\\{M_{m}\\}$
(14) $\displaystyle{\rm minimize:}$ $\displaystyle\qquad|\gamma|_{\rm
max}=\max_{m}\\{|\gamma_{m}|\\}$ (15) $\displaystyle{\rm subject~{}to:}$
$\displaystyle\qquad B\otimes A^{T}=\sum_{m}\gamma_{m}M_{m}\otimes M_{m}^{*}.$
(16)

For $A$ and $B$ that are diagonal in the same basis $\\{\left|i\right>\\}$ and
the POVM measurement described in Eq. (11), the condition in Eq. (16) reduces
to

$\bm{T}\bm{\gamma}=\alpha\bm{y},$ (17)

where $[\bm{T}]_{i+jd,m}=\langle\phi_{m}|j\rangle\langle i|\phi_{m}\rangle$,
$[\bm{\gamma}]_{m}=\gamma_{m}$, and $[\bm{y}]_{i+jd}=a_{i}b_{j}$. For the
correlation function discussed in the main text, the weight vector
$\\{\gamma_{m}(I,Z)\\}=\\{2,-2,2i,-2i\\}$ with the measeurement set
$\\{\left|\phi_{m}\right>\\}=\\{\left|0\right>,\left|1\right>,\left|+_{y}\right>,\left|-_{y}\right>,\left|+_{x}\right>,\left|-_{x}\right>\\}$
yields the minimum value of $|\gamma(I,Z)|_{\rm max}=2$. For the QPD
reconstruction, the numerical optimization for the measurement set
$\\{\left|\phi_{m}\right>\\}=\\{\left|0\right>,\left|1\right>,\left|+_{y}\right>,\left|-_{y}\right>,\left|+_{x}\right>,\left|-_{x}\right>\\}$
leads to $|\gamma(I,\Pi_{i})|_{\rm max}\approx 1.775$.