††thanks: These two authors contributed equally††thanks: These two authors
contributed equally

# Longitudinal elastic wave control by pre-deforming semi-linear materials

Dengke Guo  School of Aerospace Engineering, Beijing Institute of Technology,
100081, Beijing, China  Yi Chen  School of Aerospace Engineering, Beijing
Institute of Technology, 100081, Beijing, China  Zheng Chang  College of
Science, China Agricultural University, Beijing 100083, China  Gengkai Hu
<EMAIL_ADDRESS>School of Aerospace Engineering, Beijing Institute of
Technology, 100081, Beijing, China

###### Abstract

An incremental wave superimposed on a pre-deformed hyper-elastic material
perceives an elastic media with the instantaneous modulus of the current
material. This offers a new route with a broadband feature to control elastic
waves by purposely creating finite deformation field. This study proves that
the governing equation of a semi-linear material under a symmetric pre-
deformation condition maintains the form invariance for longitudinal wave, so
the longitudinal wave control can be made by transformation method without the
constraint condition on principle stretches, but this is not the case for
shear waves. Therefore pre-deforming a semi-linear material provides a
potential method for treating longitudinal and shear waves differently.
Examples with elastic wave control and band structure shift through pre-
deforming a semi-linear material are provided to illustrate this finding.
Finally, a one-dimensional spring lattice is proposed to mimic a semi-linear
material, and the dispersion relation for longitudinal waves in a sandwich
structure with such spring lattice is shown to be invariant during elongation,
confirming the result found based on a homogeneous semi-linear material. These
results may stimulate researches on designing new hyper-elastic
microstructures as well as designing new devices based on pre-deformed hyper-
elastic materials.

##  I INTRODUCTION

Wave steering by carefully distributing material in space is an interesting
topic for both scientific and engineering communities. A recent significant
progress along this line is the proposition of transformation methodPendry
_et al._ (2006); Milton  _et al._ (2006); Cummer and Schurig (2007); Norris
(2008), which in essence maps physical fields from one region to another based
on the form-invariance of the governing equation. Within the framework of
conventional elasticity, the elastic wave equation cannot retain its form
under a general mappingMilton  _et al._ (2006); Norris and Shuvalov (2011). An
approximate methodChang  _et al._ (2010) with conventional elastic materials
or an exact control strategyBrun  _et al._ (2009) with hypothesized elastic
materials without the elastic minor symmetry are proposed to control elastic
waves. Another interesting method is to look at the degenerated elastic
materials, for example, pentamode materials are shown to satisfy the form-
invarianceNorris (2008), but they are used to control specifically the pseudo
pressure wave instead of a fully coupled elastic wave. The transformation
method based on pentamode materials generally requires a symmetric or quasi-
symmetric mappingChen  _et al._ (2016). Realization of the materials to meet
the design requirement is of a great challenge and needs meticulous
microstructure designBückmann  _et al._ (2014); Chen  _et al._ (2017). In
addition, they may suffer from limited frequency band due to inevitable
dispersion at high frequency.

To overcome these shortcomings, an alternative promising route is proposed to
steer elastic waves by pre-deforming a hyper-elastic material. Wave on the
pre-deformed material perceives a new media with the modified local density
and instantaneous modulus of the current materialOgden (2007). This effect
purely originates from the finite deformation and hyper-elastic constitutive
relation of the hyper-elastic material. The non-dispersive local density and
instantaneous modulus may provide a new way to control elastic waves with a
broadband efficiency. By noticing the similarity between the material
parameters required by the asymmetric transformationBrun  _et al._ (2009) and
the local density and instantaneous modulusOgden (2007) of a pre-deformed
hyper-elastic material, the hyper-elastic transformation theoryParnell (2012);
Norris and Parnell (2012); Chang  _et al._ (2015) is proposed to control the
superimposed incremental waves by pre-deforming a hyper-elastic material. Neo-
Hookean1Parnell (2012); Parnell  _et al._ (2012) hyper-elastic materials have
been investigated to design a cylindrical cloak for anti-plane shear waves,
the cloaking effect is demonstrated from scattering analysis. This hyper-
elastic material is also used to design a phononic crystal with invariant band
gapsBarnwell  _et al._ (2016) for shear waves. Chang et al.Chang  _et al._
(2015) further studied the in-plane longitudinal and shear waves in a
compressible neo-Hookean material under a pre-shear deformation, and they
found that the in-plane shear waves will follow the design by the hyper-
elastic transformation while the longitudinal waves will not be affected.
However, it is only rigorously proved that incompressible neo-Hookean
materials maintain the form invariance for anti-plane shear wavesParnell
(2012); Parnell  _et al._ (2012); Barnwell  _et al._ (2016). Norris and
Parnell proposed to use a semi-linear hyper-elastic material with a symmetric
pre-deformation to control fully coupled elastic wavesNorris and Parnell
(2012), and a cylindrical cloak is designed for both in-plane longitudinal and
shear waves. Based on the semi-linear energy function in terms of three
principal stretches, they find that the out of plane extensions should be
adapted accordingly with the in-plane stretches in order to control the in-
plane longitudinal and shear wavesNorris and Parnell (2012), simultaneously.
In this paper, we observe that the second order derivatives of the semi-linear
energy function are composed of the deformation gradients and a fourth order
tensor with the minor anti-symmetry. This forth order tensor with the minor
anti-symmetry is cancelled out with a divergence field, implying that a semi-
linear material, under a general symmetric deformation without the out of
plane extensions constraint, still maintains the form invariance for a
divergence field. Therefore the transformation method can be used to control
longitudinal waves. This finding provides a useful tool to design devices
interacting differently with the longitudinal and shear components of a fully
coupled elastic wave.

This research is arranged as follows. First, we will prove that the divergence
component of elastic waves on a symmetrically pre-deformed semi-linear
material can be strictly mapped from the initial configuration. These results
indicate that, the governing equation for longitudinal waves in a semi-linear
material under a symmetric pre-deformation maintain the form invariance, so
they can be controlled within the strategy of transformation method by pre-
deformation. Second, numerical simulations, demonstrating wavelength
magnifying and dispersion relation invariance for longitudinal waves on a pre-
deformed semi-linear material, will be conducted to validate the
transformation method. Finally, a spring lattice is proposed as a one-
dimensional (1D) prototype for a homogeneous semi-linear material, the
invariant dispersion relation for longitudinal waves under a finite elongation
will be demonstrated through a 1D sandwich structure to confirm the
theoretical finding.

##  II METHODS

###  II.1 Wave propagation on a deformed hyper-elastic material

A stress free hyper-elastic material occupying an initial configuration
$\Omega_{0}$ is deformed to the current configuration $\Omega$ through the
boundary and body loads. Mathematically, this process is represented by a
mapping between the two configurations
$\mathbf{x}:\Omega_{0}\rightarrow\Omega$, which maps any point $\mathbf{X}$ in
$\Omega_{0}$ to another point $\mathbf{x}$ in $\Omega$,
$\mathbf{x}=\mathbf{x}(\mathbf{X})$. This mapping can also be denoted by the
deformation gradient $\mathbf{F}=\partial\mathbf{x}/\partial\mathbf{X}$ in
nonlinear elasticityOgden (2007), or in component form $F_{iJ}=\partial
x_{i}/\partial X_{J}$. Here, lowercase and uppercase subscripts refer to the
initial configuration and current configuration, respectively. The deformation
tensor $\mathbf{F}$ in general can take any form compatible with physical
deformation, e.g., identity, diagonal or symmetric form. A small displacement
perturbation $\mathbf{u}(\mathbf{x},t)$ superimposed on the deformed body
$\Omega$ is governed by the dynamic equation in the current configurationOgden
(2007),

$\rho u_{j,tt}=\nabla_{i}(C_{ijkl}\nabla_{k}u_{l})$ (1)

Here, $i$ with a lowercase subscript representing the partial derivative with
respect to the current coordinate $\mathbf{x}_{i}$, the displacements ul with
a lowercase subscript mean the fields physically occurred in the current
configuration, and subscript $t$ denotes time derivative. The used subscripts
range from 1 to 3 and duplicated indexes should be understood as summation.
The local density $\rho$ and instantaneous elastic tensor $\mathbf{C}$ in the
current configuration are,

$\rho=J^{-1}\rho_{0}$ (2a)
$C_{ijkl}=J^{-1}F_{iM}F_{kN}A_{MjNl}=J^{-1}F_{iM}F_{kN}\frac{\partial^{2}W}{\partial
F_{jM}\partial F_{lN}}$ (2b)

In which, $J$ is determinant of the deformation tensor $\mathbf{F}$, $\rho$ is
derived from the mass conservation law and $\rho_{0}$ represents the density
in the initial configuration, $W$ is the strain energy function of the hyper-
elastic material. The instantaneous elastic tensor $\mathbf{C}$ is a forth
order tensor without the minor symmetry, and the incremental wave is likely to
be governed by an elastic material with asymmetric incremental stress
$C_{ijkl}\nabla_{k}u_{l}$. The true Cauchy stress in the hyper-elastic
material is indeed still symmetric, only the incremental stress is asymmetric.
This observation stimulates researches on the transformation theory for
elastic waves through a finite deformation in a hyper-elastic materialParnell
(2012); Norris and Parnell (2012); Chang  _et al._ (2015). It should be noted,
although the actual perturbation occurs on the current configuration $\Omega$,
its governing equation can be transformed to the initial configuration
$\Omega_{0}$. Following the finite deformation elasticity theoryOgden (2007),
Eq. (1) can be easily pulled back to the initial configuration,

$\rho_{0}u_{j,tt}=J\nabla_{i}(J^{-1}F_{iM}F_{kN}A_{MjNl}\nabla_{k}u_{l})=\nabla_{M}(A_{MjNl}\nabla_{N}u_{l})$
(3)

In above equation, two mathematic identities, $\nabla_{i}(J^{-1}F_{iM})=0$,
$F_{iM}\nabla_{i}=\nabla_{M}$ are used. $\nabla_{M}$ with an uppercase
subscript represents the partial derivative with respect to the initial
coordinate $\mathbf{X}_{M}$. It is seen that the modulus,
$A_{MjNl}=\partial^{2}W/\partial F_{jM}/\partial F_{lN}$, functions as a
nominal elastic modulus in the initial configuration, and therefore is called
the pull-backed modulus in the following.

For the same hyper-elastic material without the pre-deformation, wave
propagation is governed by the equation in the initial configuration
$\Omega_{0}$,

$\begin{split}\rho_{0}v_{J,tt}=&\nabla_{M}(C_{0mjnl}\nabla_{N}v_{l})\\\
=&(\lambda_{0}+\mu_{0})\nabla_{J}(\nabla_{L}v_{L})+\mu_{0}\nabla_{M}\nabla_{M}v_{J}\end{split}$
(4)

Here, $\mathbf{v}(\mathbf{X},t)$ is used for the displacement fields on the
hyper-elastic material without the pre-deformation to distinguish from the
displacement field $\mathbf{u}(\mathbf{x},t)$ with the pre-deformation. The
displacements $v_{J}$ with an uppercase subscript mean the fields physically
occurred in the initial configuration. Apparently, the above equation is a
special case of Eq. (3) with the deformation tensor being an identity matrix
$\mathbf{F}=\mathbf{I}$, and the elastic tensor for a conventional isotropic
elastic material is recovered, i.e.,
$C_{0mjnl}=\lambda_{0}\delta_{mj}\delta_{nl}+\mu_{0}(\delta_{mn}\delta_{jl}+\delta_{ml}\delta_{jn})$
with $\lambda_{0}$ and $\mu_{0}$ being the lamé constants.

###  II.2 Transformation method for longitudinal waves with semi-linear
materials

In the following, we will prove the governing equation for longitudinal waves
is form-invariant for semi-linear materials under a symmetric deformation. The
strain energy function of a semi-linear material writesNorris and Parnell
(2012),

$W=\frac{\lambda_{0}}{2}(U_{KK}-3)^{2}+\mu_{0}(U_{KL}-\delta_{KL})(U_{KL}-\delta_{KL})$
(5)

In which, $\mathbf{U}=\mathbf{R}^{-1}\mathbf{F}$ is the right stretch tensor
for the finite deformation gradient $\mathbf{F}$ and $\mathbf{R}$ is an
orthogonal matrix $\mathbf{R}^{\mathbf{T}}\mathbf{R}=\mathbf{I}$. The strain
energy function $W$ is only related to the stretch tensor $\mathbf{U}$ as
required by objectivity. The pull-backed modulus $A_{MjNl}$ for a general
deformation gradient $\mathbf{F}$ is given in AppendixA, and can be further
simplified when $\mathbf{F}$ is symmetric with $\mathbf{R}=\mathbf{I}$,

$\begin{split}A_{MjNl}|_{\mathbf{F}=\mathbf{F}^{\mathbf{T}}}=&\lambda_{0}\delta_{jM}\delta_{lN}+2\mu_{0}\delta_{ij}\delta_{MN}\\\
+&(\lambda_{0}(U_{KK}-\delta_{KK})-2\mu_{0})\frac{\partial R_{jM}}{\partial
F_{lN}}\end{split}$ (6)

In addition to the first two constant terms, the pull-backed modulus has a
third term depending on the deformation gradient $\mathbf{F}$. From the
geometric point of view, the symmetry of the deformation tensor $\mathbf{F}$
means the deformation field is irrotational. If $\mathbf{F}$ further becomes
an identity matrix $\mathbf{I}$, the modulus $A_{MjNl}$ will degenerate to the
conventional isotropic elasticity tensor (see more details in AppendixA).
Substituting the above pull-backed modulus Eq. (6) into Eq. (3) and taking the
divergence operation on both sides, one will obtain the following equation,

$\begin{split}\rho\nabla_{J}u_{j,tt}&=\nabla_{J}\nabla_{M}(A_{MjNl}\nabla_{N}u_{l})\\\
&=(\lambda_{0}+2\mu_{0})\nabla_{M}\nabla_{M}(\nabla_{L}u_{l})\end{split}$ (7)

In deriving the above equation, the skew-symmetry of $\partial R_{jM}/\partial
F_{lN}$ with respect to index $j$ and $M$ is used. This equation is exactly
the same as the governing equation for the divergence field on the initially
un-deformed configuration, i.e.,
$\rho_{0}\nabla_{J}v_{J,tt}=(\lambda_{0}+2\mu_{0})\nabla_{M}\nabla_{M}(\nabla_{J}v_{J})$.
In this case, once the boundary and loading conditions for the two
configurations are set accordinglyChen  _et al._ (2016), the divergence field
on the pre-deformed body can be mapped from that on the initial configuration
$\nabla_{L}u_{l}=F_{kJ}(\nabla_{k}u_{j})=\nabla_{J}v_{J}$. This observation
implies that, if one region of a semi-linear material is pre-deformed without
altering its boundary, then an impinged longitudinal wave inside the pre-
deformed region will follow the mapping, while outside the pre-deformed region
it will not be influenced at all and just perceives a homogeneous semi-linear
material. It can be proved similarly that the curl field on the pre-deformed
semi-linear material doesn’t follow the same equation as on the initial
configuration, and therefore cannot be mapped. This property can be used to
design wave controlling devices interacting differently with longitudinal and
shear waves.

A more explicit expression for $A_{MjNl}$ is not available for a general
symmetric gradient, however, a very simple formula for the in-plane components
is possible for the special symmetric gradients with $F_{13}=F_{23}=0$,
$F_{33}\not=0$,

$\begin{split}A_{MjNl}|_{\mathbf{F}=\mathbf{F}^{\mathbf{T}}}=&\lambda_{0}\delta_{jM}\delta_{lN}+2\mu_{0}\delta_{ij}\delta_{MN}\\\
+&(\lambda_{0}(U_{kk}+U_{33}-3)-2\mu_{0})\frac{1}{F_{kk}}\epsilon_{jM}\epsilon_{lN}\end{split}$
(8)

Here, all the index $M$, $j$, $N$, $l$ and $k$ range from 1 to 2,
$\epsilon_{jM}$ is the 2D permutation tensor. Further if the out of plane
stretch is constraint accordingly with the in-plane stretches,

$U_{33}=1-\frac{\lambda_{0}+\mu_{0}}{\lambda_{0}}(U_{kk}-2)$ (9)

Equation (9) is the same as (4.14) by Norris and ParnellNorris and Parnell
(2012). Substituting Eq. (9) into Eq. (8), and noticing $U_{ηη}=F_{ηη}$, one
will obtain an isotropic in-plane elastic tensor,

$A_{MjNl}|_{\mathbf{F}=\mathbf{F}^{\mathbf{T}}}=\lambda_{0}\delta_{jM}\delta_{lN}+\mu_{0}(\delta_{ij}\delta_{MN}+\delta_{jN}\delta_{Ml})$
(10)

So with the constraint condition of Eq. (9), both longitudinal and in-plane
shear waves can be controlled simultaneously with transformation method, as
observed by Norris and ParnellNorris and Parnell (2012). However we
demonstrate here, for longitudinal waves the constraint condition can be
removed.

##  III NUMERICAL EXAMPLES

###  III.1 Wave propagation on a symmetrically pre-deformed semi-linear
material

Numerical simulation based on the hyper-elastic transformation can be
performed in three steps: First, a stress free body in the initial
configuration $\Omega_{0}$ is deformed to the current configuration $\Omega$
by loading or an imaginary mapping, and then the deformation tensor
$\mathbf{F}$ is obtained. Second, the density $\rho$ and instantaneous elastic
tensor $\mathbf{C}$ in the current configuration are calculated using Eq.
(2b). Third, the wave equation will be solved on the current configuration
with the derived local material parameters. The numerical simulations are all
conducted using PDE (Partial Differential Equation) and Nonlinear Mechanics
module in the commercial software COMSOL Multiphysics.

Figure 1: Wave propagation on the stress-free and pre-deformed semi-linear
materials. (a), (b) Divergence fields for the stress-free and pre-deformed
cases with longitudinal wave excitation. (c),(d) Curl fields in the stress-
free and pre-deformed cases with shear wave excitation. Figure 2: Band
structure of the stress-free and pre-deformed semi-linear materials. (a)
Stress free semi-linear material case and the first four modes correspondent
to a wave vector $\mathbf{k}=\boldsymbol{\Gamma}\mathbf{X}/6$. (b) Semi-linear
material with pre-deformed circular region and the first four mode
correspondent to a wave vector $\mathbf{k}=\boldsymbol{\Gamma}\mathbf{X}/6$.

We simulate the wave propagation in a pre-deformed semi-linear material that
can magnify the wavelength. As shown in FIG. 1(a), the initial configuration
consists of two regions, an inner circular region ($R<b=11$ m), where $b$ is
the radius of the outer boundary (black dashed line) for the pre-deformed
region, and an outer annular region ($b<R<24$ m). Perfectly Matched Layers
adjacent to the outer circular boundary is employed to mimic an infinite large
space without reflection. With the mapping $r=R^{4}/b^{3}$ for the inner
circular region ($R<b$), this region is transformed to another circular region
with the same radius b as shown in FIG. 1(b). The radial mapping ensures the
symmetric deformation gradient $\mathbf{F}$ and therefore the condition of
transformation method for longitudinal wave holds. A tiny circle placed at
$r=20$ m can excite longitudinal or shear waves ($f=6.5$ Hz) by expansion or
rotation motions, respectively. Material parameters of the semi-linear
material are $\lambda_{0}=770$ Pa, $\mu_{0}=260$ Pa, $\rho_{0}=1$ kg/m3. In
the simulation, the longitudinal and shear wave excitations are modeled
through expanding and rotating boundary of the small circle in the domain.
Figures 1(a) and 1(b) show divergence fields for the stress free and pre-
deformed cases with the longitudinal wave excitation, respectively. Exactly
the same fields are observed outside the pre-deformed region as predicted
theoretically. Further, the divergence fields in the pre-deformed region
($r<b$) can be perfectly mapped from the fields in the initial configuration
with $\nabla_{i}u_{i}=F^{-1}_{kL}\nabla_{K}v_{L}$. While, the curl fields for
the stress free and pre-deformed cases (FIG. 1(c),(d)) with the shear wave
excitation are completely different, and significant scattering caused by the
pre-deformed region is observed. This interesting phenomenon indicates that
the pre-deformation region is invisible for longitudinal waves, but is a
strong scatter for shear waves.

Transformation method for longitudinal waves with a pre-deformed semi-linear
material is further demonstrated with calculation of dispersion relation. Two
square unit cells of the same sizes with a lattice constant $a=1$ m are
considered, and the periodic medium are constructed from tiling unit cells
along two orthogonal directions. One is composed of a homogeneous stress free
semi-linear material, while the other has its inner region ($r<0.35$ m) pre-
deformed from the same circular region ($R<0.35$ m) with a mapping function
$r=R(R/0.35)^{1/2}$. The Bloch wave frequency can be obtained by solving the
eigen-frequency problem on one unit cell with the Bloch wave boundary
condition. The dispersion relation is obtained by sweeping the Bloch wave
vector along the edges of the first irreducible Brillouin Zone
($\mathbf{X}-\boldsymbol{\Gamma}-\mathbf{M}-\mathbf{X}$). Figure 2(a) is the
band structure of the stress-free homogenous unit cell, and FIG. 2(b) presents
the result of the partly pre-deformed semi-linear material unit cell. Wave
polarizations for different branches can be identified from the modes shown
below the band structures. The second branch is longitudinal and the mode
shows an overall displacement parallel to the wave vector, while all other
three branches represent shear waves. From both band structures, it is shown
that the purely longitudinal branches (marked by red dots) are exactly the
same for the stress-free and pre-deformed unit cells, while the shear related
branches become quite different for both cases. This is expected from the
longitudinal wave invariance property of the symmetrically pre-deformed semi-
linear material.

Figure 3: Spring lattice as a 1D semi-linear material.(a) Unit cell of the
U-shaped spring and sandwich spring composite. (b) The deformation field
($\nabla\cdot\mathbf{u}$) of the composite with a finite strain
($\lambda=1.5$). (c) Mechanical property of 10 connected unit cells in (a).
(d) Enlarged view of the indicated cyan region in (c) with a strain range
($\lambda=2.3\sim 3.0$).

###  III.2 Prototype of a 1D semi-linear material with spring lattice

Inspired by the Hooke’s law of a perfect spring, a 1D spring lattice is
proposed to mimic the semi-linear material. As shown in FIG. 3(a), a unit cell
of the proposed spring lattice is composed of connected beams with the
following geometry parameters, $t=0.012$ m, $h=0.2$ m and $w=0.04$ m, and all
the beams are made of the same material. A sandwiched spring composite (SSC)
can be constructed by filling the springs into two separated homogenous
materials. A unit cell of the SSC consists of multiple spring unit cells
($5\times 10$, $l=10w$) and two layers of the homogeneous materials
($l_{1}=l$, $l_{2}=5h$).

The hyper-elastic strain energy function $W(F_{11})$ of the proposed 1D semi-
linear material is obtained numerically by solving the stored strain energy in
10 connected unit cells with the applied displacement boundary condition,
$u=(F_{11}-1)10w$, and shows quadratic relation (marked by squares in FIG.
3(c)) with respect to the deformation gradient $F_{11}=\lambda$ for a finite
deformation range $F_{11}=1.0\sim 3.0$ as compared to the ideal semi-linear
case (red line). Two important points should be noted in the numerical
simulation. First, since finite rotations occur in this problem during
elongation, the geometric nonlinear effect must be turned on during the
numerical simulation. Second, although the overall elongation ratio of the
unit cell is no longer small, the local strains in the springs and the
homogeneous layers are still quite small (the Green strain is of the order of
magnitude of $10^{-3}$ for a elongation ratio $\lambda=3$), so a linear
elastic constitutive relation is set for the spring material. Here, the
constitutive relation is taken to be the default setting in COMSOL with
geometry nonlinearity, i.e., St-Venant-Kirchohoff constitutive relation. In
short, the hyper-elastic behavior of the proposed 1D spring lattice comes from
the finite rotation but small strain of the integrated microstructure. The
numerically computed pull-backed modulus $A_{1111}=\partial^{2}W/\partial
F_{11}/\partial F_{11}$ (marked by triangles) is nearly constant for the
investigated finite strain range. The highlighted region in FIG. 3(c) is shown
more clearly in FIG. 3(d), where the modulus changes less than 4 percent from
a constant value. The numerical result shows the proposed spring lattice has a
nearly constant pull-backed modulus under the elongation $F_{11}$, and
therefore can be used with the hyper-elastic transformation method to control
longitudinal wave along the elongation direction. As for its impact on the
longitudinal wave in the other directions, the hyper-elastic strain energy
function should be studied for a general deformation gradient and this is
beyond the main scope of this paper. Finally, since the pull-backed modulus is
obtained from the static homogenization or under long wave approximation, so
in the following analysis, we consider only the low frequency case where the
static homogenization condition holds.

Figure 4: Band structure of the proposed 1D semi-linear material. (a) Band
structure of the SSC unit cell without (line) and with (discrete square) pre-
deformation. (b) Modes of the first four branches correspondent to a
normalized wave number $k_{x}L_{0}=\pi/2$, left panel is for the stress-free
case and right panel is for the pre-deformed case.

Figure 3(b) shows the deformed state of the SSC with a finite strain
$\lambda=1.5$, where the deformation mainly concentrates in the turning points
of the springs and the strain in the homogenous material is extremely small,
therefore the springs can be regarded as pre-deformed while the layered
materials remain in its initial configuration. The dispersion relations for
the un-deformed and pre-deformed SSC unit cells (FIG. 3(b)) are shown in FIG.
4(a) together with their corresponding deformation modes in FIG. 4(b). Since
the unit cells of the un-deformed and pre-deformed cases have different
lengths, the horizontal axis in FIG. 4(a) represents a non-dimensional wave
number $k_{x}L_{0}$ with $L_{0}$ being the length of the corresponding unit
cell. The true deformation is rather small and has been exaggerated in the
figures for clarity. The modes indicate that the second and forth branches are
mainly related to the longitudinal waves and therefore these two branches of
the pre-deformed case are almost the same as the stress-free case. While the
first and third branches are dominated by shear waves as shown by the modes,
and are changed during stretching, especially for the first branch. The above
results show that the SSC can be used as a robust semi-linear material for
controlling longitudinal waves with the hyper-elastic transformation method.

##  IV CONCLUSIONS

The hyper-elastic transformation method with semi-linear materials is
revisited. It is proved that semi-linear materials under a general symmetric
pre-deformation maintain the form invariance for longitudinal elastic waves,
while not for shear waves. This implies that the longitudinal wave control can
be made with the hyper-elastic transformation method without the additional
constraint on pre-stretches. Wave simulations on a pre-deformed semi-linear
material show that, the longitudinal waves follow exactly the designed path by
the hyper-elastic transformation method without any scattering, while the
shear waves are strongly scattered. The dispersion relation of a semi-linear
material, with a symmetric pre-deformation, is shown invariant for the
longitudinal wave branches, but not for the shear related branches. These
numerical results confirm the theoretical finding. Finally, a 1D semi-linear
material realized by a spring lattice is proposed to mimic a semi-linear
material. Numerical simulations reveal that the longitudinal wave branches
remain intact during stretching. The demonstrated property of the semi-linear
material may find applications where longitudinal and shear waves are expected
to be controlled differently with external deformation.

###### Acknowledgements.

This work was supported by National Natural Science Foundation of China (Grant
Nos. 11472044, 11521062, 11632003).

*

##  Appendix A DERIVATIVE OF STRAIN ENERGY FUNCTION

Following the polar decomposition of a deformation tensor
$\mathbf{F}=\mathbf{R}\mathbf{U}$, we have
$\mathbf{F}^{\mathbf{T}}\mathbf{F}=\mathbf{U}^{2}$ or in component form
$U_{IK}U_{KJ}=F_{kI}F_{kJ}$. Taking the derivative of both sides with respect
to $F_{lN}$ leads to,

$\frac{\partial U_{IK}}{\partial F_{lN}}U_{KJ}+\frac{\partial U_{KJ}}{\partial
F_{lN}}U_{IK}=F_{lI}\delta_{Nj}+F_{lJ}\delta_{Nl}$ (11)

From det($\mathbf{U}$)=det($\mathbf{F}$)>0, one can deduce that $\mathbf{U}$
is reversible. Multiplying both sides with $\delta_{IJ}$ and $U^{-1}_{IJ}$
gives, respectively,

$\frac{\partial U_{IK}}{\partial F_{lN}}U_{KI}=F_{lN},\quad\frac{\partial
U_{KK}}{\partial F_{lN}}=F_{lI}U^{-1}_{IN}=R_{lN}$ (12)

Further differentiating Eq. (11) with respect to $F_{rS}$ and multiplying both
sides with $U^{-1}_{IJ}$ gives,

$\frac{\partial R_{lN}}{\partial F_{rS}}+\frac{\partial U_{IK}}{\partial
F_{lN}}\frac{\partial U_{KJ}}{\partial
F_{rS}}U^{-1}_{IJ}=\delta_{lr}U^{-1}_{SN}$ (13)

If $\mathbf{F}$ becomes an identity matrix, Eq. (13) can be simplified as,

$\frac{\partial R_{lN}}{\partial
F_{rS}}|_{\mathbf{F}=\mathbf{I}}=\frac{1}{2}(\delta_{lr}\delta_{SN}-\delta_{lS}\delta_{rN})$
(14)

With above formulas in Eq. (12), the second derivative of the semi-linear
strain energy function can be derived,

$\begin{split}A_{MjNl}&=\frac{\partial^{2}W}{\partial F_{jM}\partial
F_{lN}}=\frac{\partial}{\partial F_{jM}}\frac{\partial}{\partial
F_{lN}}(\frac{\lambda_{0}}{2}(U_{KK}-\delta_{KK})^{2}+\mu_{0}(U_{KL}U_{KL}-2U_{KK}))\\\
&=\frac{\partial}{\partial
F_{jM}}(\lambda_{0}(U_{KK}-\delta_{KK})\frac{\partial U_{KK}}{\partial
F_{lN}}+2\mu_{0}U_{KL}\frac{\partial U_{KL}}{\partial
F_{lN}}-2\mu_{0}\frac{\partial U_{KK}}{\partial F_{lN}})\\\
&=\frac{\partial}{\partial
F_{jM}}((\lambda_{0}U_{KK}-\lambda_{0}\delta_{KK}-2\mu_{0})R_{lN}+2\mu_{0}F_{lN})\\\
&=\lambda_{0}R_{jM}R_{lN}+2\mu_{0}\delta_{lj}\delta_{MN}+(\lambda_{0}U_{KK}-\lambda_{0}\delta_{KK}-2\mu_{0})\frac{\partial
R_{lN}}{\partial F_{jM}}\end{split}$ (15)

If the deformation tensor becomes symmetric
$\mathbf{F}=\mathbf{F}^{\mathbf{T}}$, we can simplify the pull-backed modulus
Eq. (15) as,

$\begin{split}A_{MjNl}|_{\mathbf{F}=\mathbf{F}^{\mathbf{T}}}=&\lambda_{0}\delta_{jM}\delta_{lN}+2\mu_{0}\delta_{ij}\delta_{MN}\\\
+&(\lambda_{0}(U_{KK}-\delta_{KK})-2\mu_{0})\frac{\partial R_{lN}}{\partial
F_{jM}}\end{split}$ (16)

Notice that, $\partial R_{lN}/\partial F_{jM}=\partial R_{jM}/\partial F_{lN}$
is skew-symmetric with respect to index $j$ and $M$ when $\mathbf{F}$ is
symmetric. It can be easily proved from this mathematical identity
$\mathbf{R}\mathbf{R}^{\mathbf{T}}=\mathbf{I}$ by differentiation,

$\frac{\partial R_{iK}}{\partial F_{lN}}R_{jK}+R_{iK}\frac{\partial
R_{jK}}{\partial F_{lN}}=0$ (17a) $\frac{\partial R_{iJ}}{\partial
F_{lN}}|_{\mathbf{F}=\mathbf{F}^{\mathbf{T}}}+\frac{\partial R_{jI}}{\partial
F_{lN}}|_{\mathbf{F}=\mathbf{F}^{\mathbf{T}}}=0$ (17b)

Further if $\mathbf{F}$ becomes an identity matrix, the pull-backed modulus
becomes the conventional isotropic one,

$\begin{split}A_{MjNl}|_{\mathbf{F}=\mathbf{F}^{\mathbf{T}}}=&\lambda_{0}\delta_{jM}\delta_{lN}+2\mu_{0}\delta_{lj}\delta_{MN}-2\mu_{0}\frac{\partial
R_{lN}}{\partial F_{jM}}|_{\mathbf{F}=\mathbf{I}}\\\
=&\lambda_{0}\delta_{jM}\delta_{lN}+\mu_{0}\delta_{lj}\delta_{MN}+\mu_{0}\delta_{lM}\delta_{jN}\end{split}$
(18)

## References

  * Pendry  _et al._ (2006) J. B. Pendry, D. Schurig, and D. R. Smith, science 312, 1780 (2006).
  * Milton  _et al._ (2006) G. W. Milton, M. Briane, and J. R. Willis, New Journal of Physics 8, 248 (2006).
  * Cummer and Schurig (2007) S. A. Cummer and D. Schurig, New Journal of Physics 9, 45 (2007).
  * Norris (2008) A. N. Norris, Proc. R. Soc. A 464, 2411 (2008).
  * Norris and Shuvalov (2011) A. Norris and A. Shuvalov, Wave Motion 48, 525 (2011).
  * Chang  _et al._ (2010) Z. Chang, X. Zhou, J. Hu, and G. Hu, Optics express 18, 6089 (2010).
  * Brun  _et al._ (2009) M. Brun, S. Guenneau, and A. B. Movchan, Applied Physics Letters 94, 061903 (2009).
  * Chen  _et al._ (2016) Y. Chen, X. Liu, and G. Hu, The Journal of the Acoustical Society of America 140, EL405 (2016).
  * Bückmann  _et al._ (2014) T. Bückmann, M. Thiel, M. Kadic, R. Schittny, and M. Wegener, Nature communications 5, 4130 (2014).
  * Chen  _et al._ (2017) Y. Chen, M. Zheng, X. Liu, Y. Bi, Z. Sun, P. Xiang, J. Yang, and G. Hu, Physical Review B 95, 180104 (2017).
  * Ogden (2007) R. W. Ogden, COURSES AND LECTURES 495, 1 (2007).
  * Parnell (2012) W. J. Parnell, Proc. R. Soc. A 468, 563 (2012).
  * Norris and Parnell (2012) A. N. Norris and W. J. Parnell, Proc. R. Soc. A 468, 2881 (2012).
  * Chang  _et al._ (2015) Z. Chang, H.-Y. Guo, B. Li, and X.-Q. Feng, Applied Physics Letters 106, 161903 (2015).
  * Parnell  _et al._ (2012) W. J. Parnell, A. N. Norris, and T. Shearer, Applied Physics Letters 100, 171907 (2012).
  * Barnwell  _et al._ (2016) E. G. Barnwell, W. J. Parnell, and I. D. Abrahams, Wave Motion 63, 98 (2016).


# High-precision force sensing using a single trapped ion

Peter A. Ivanov  Department of Physics, St. Kliment Ohridski University of
Sofia, James Bourchier 5 blvd, 1164 Sofia, Bulgaria  Nikolay V. Vitanov
Department of Physics, St. Kliment Ohridski University of Sofia, James
Bourchier 5 blvd, 1164 Sofia, Bulgaria  Kilian Singer  Experimentalphysik I,
Universität Kassel, Heinrich-Plett-Str. 40, D-34132 Kassel, Germany

(August 27, 2024)

###### Abstract

We introduce quantum sensing schemes for measuring very weak forces with a
single trapped ion. They use the spin-motional coupling induced by the laser-
ion interaction to transfer the relevant force information to the spin-degree
of freedom. Therefore, the force estimation is carried out simply by observing
the Ramsey-type oscillations of the ion spin states. Three quantum probes are
considered, which are represented by systems obeying the Jaynes-Cummings,
quantum Rabi (in 1D) and Jahn-Teller (in 2D) models. By using dynamical
decoupling schemes in the Jaynes-Cummings and Jahn-Teller models, our force
sensing protocols can be made robust to the spin dephasing caused by the
thermal and magnetic field fluctuations. In the quantum-Rabi probe, the
residual spin-phonon coupling vanishes, which makes this sensing protocol
naturally robust to thermally-induced spin dephasing. We show that the
proposed techniques can be used to sense the axial and transverse components
of the force with a sensitivity beyond the yN $/\sqrt{\text{Hz}}$ range, i.e.
in the xN$/\sqrt{\text{Hz}}$ (xennonewton, $10^{-27}$). The Jahn-Teller
protocol, in particular, can be used to implement a two-channel vector
spectrum analyzer for measuring ultra-low voltages.

###### pacs:

03.67.Ac, 03.67.Bg, 03.67.Lx, 42.50.Dv

##  I Introduction

Over the last few years, research of mechanical systems coupled to quantum
two-level systems has attracted great deal of experimental and theoretical
interest Treutlein2015 ; Kurizki2015 . Micro- and nano-mechanical oscillators
can respond to very weak electric, magnetic and optical forces, which allows
one to use them as highly sensitive force detectors Stowe1997 . For example,
the cantilever with attonewton ($10^{-18}$ N) force sensitivity can be used to
test the violation of Newtonian gravity at sub-millimeter length scale
Geraci2008 . With current quantum technologies coupling between a
nanomechanical oscillator and a single spin can be achieved experimentally by
using strong magnetic-field gradient. Such a coupling paves the way for
sensing the magnetic force associated with the single electron spin Rugar2004
. To this end, a recent experiment demonstrated that the coherent evolution of
the electronic spin of an individual nitrogen vacancy center can be used to
detect the vibration of a magnetized mechanical resonator Kolkowitz2012 .

Another promising quantum platform with application in high-precision sensing
is the system of laser-cooled trapped ions, which allows excellent control
over the internal and motional degrees of freedom Blatt2008 . Force
sensitivity of order of 170 yN $\rm{Hz}^{-1/2}$ ($10^{-24}$ N) was reported
recently with an ensemble of ions in a Penning trap Biercuk2010 . Force
measurement down to $5$ yN has been demonstrated experimentally using the
injection-locking technique with a single trapped ion Knunz2010 . Moreover,
force detection with sensitivity in the range of 1 yN $\rm{Hz}^{-1/2}$ is
possible for single-ion experiments based on the measurement of the ion’s
displacement amplitude Maiwald2009 .

In this work, we propose ion-based sensing schemes for measuring very rapidly
varying forces, which follow an earlier proposal Ivanov2015  wherein the
relevant force information is mapped into the spin degrees of freedom of the
single trapped ion. In contrast to Ivanov2015 , the techniques proposed here
do not require specific adiabatic evolution of the control parameters but
rather they rely on using Ramsey-type oscillations of the ion’s spin states,
which are detected via state-dependent fluorescence measurements. Moreover, we
show that by using dynamical decoupling schemes, the sensing protocols become
robust against dephasing of the spin states caused by thermal and magnetic-
field fluctuations.

We consider a quantum system described by the Jaynes-Cummings (JC) model which
can be used as a highly sensitive quantum probe for sensing of the axial force
component. By applying an additional strong driving field Zheng ; Solano2003
the dephasing of the spin states induced by the residual spin-phonon
interaction can be suppressed such that the sensing protocol does not require
initial ground-state cooling of the ion’s vibrational state. We show that the
axial force sensing can be implemented also by using a probe represented by
the quantum Rabi (QR) model. Because of the absence of residual spin-motional
coupling in this case, the force estimation is robust to spin dephasing
induced by the thermal motion fluctuations.

Furthermore, we introduce a sensing scheme capable to extract the two-
dimensional map of the applied force. Here the quantum probe is represented by
the Jahn-Teller (JT) model, in which the spin states are coupled with phonons
in two spatial directions. We show that the two transverse components of the
force can be measured by observing simply the coherent evolution of the spin
states. In order to protect the spin coherence during the force estimation we
propose a dynamical decoupling sequence composed of phonon phase-shift
operators, which average to zero the residual spin-phonon interaction.

We estimate the optimal force sensitivity in the presence of motional heating
and find that with current ion trap technologies force sensitivity better than
1 yN $\rm{Hz}^{-1/2}$ can be achieved. Thus, a single trapped ion may serve as
a high-precision sensor of very weak electric fields generated by small needle
electrodes with sensitivity as low as 1 $\mu{\rm V}/{\rm m}$ $\rm{Hz}^{-1/2}$.

This paper is arranged as follows. In Sec. II we describe the sensing protocol
for detection of the axial component of very weak forces using a quantum probe
represented by the Jaynes-Cummings and quantum Rabi models. It is shown that
by using dynamical decoupling technique the sensing protocol using Jaynes-
Cummings system is immune to thermal spin dephasing. In Sec. III we introduce
a sensing scheme, which is able to detect the two components of the external
force. Finally, in Sec. IV we summarize our findings.

##  II 1D Force Sensing

###  II.1 Jaynes-Cummings quantum probe

In our model we consider a single two-state ion with a transition frequency
$\omega_{0}$, in a linear Paul trap with an axial trap frequency $\omega_{z}$.
The small axial oscillation of the ion is described by the vibrational
Hamiltonian $\hat{H}_{\rm ax}=\hbar\omega_{z}\hat{a}^{{\dagger}}\hat{a}$,
where $\hat{a}^{{\dagger}}$ ($\hat{a}$) creates (annihilates) a phonon
excitation. We assume that the ion interacts with a laser field with a
frequency $\omega_{\rm L}=\omega_{0}-\omega_{z}+\delta$, tuned near the red-
sideband resonance with a detuning $\delta$. The interaction Hamiltonian in
the Lamb-Dicke limit and the rotating-wave approximation reads Wineland1998 ;
Haffner2008 ; Schneider2012

$\hat{H}_{\rm
JC}=\hbar\omega\hat{a}^{{\dagger}}\hat{a}+\hbar\Delta\sigma_{z}+\hbar
g(\sigma^{-}\hat{a}^{{\dagger}}+\sigma^{+}\hat{a}),$ (1)

with $\delta=\Delta-\omega$, where $\Delta$ is the effective spin frequency
and $\omega$ is the effective phonon frequency. Here, $\sigma_{x,y,z}$ are the
Pauli matrices, $\sigma^{\pm}$ are the respective raising and lowering
operators for the effective spin system, and $g$ determines the strength of
the spin-phonon coupling.

The external time-varying force with a frequency
$\omega_{d}=\omega_{z}-\omega$, e.g., $F(t)=F\cos(\omega_{d}t)$, displaces the
motional amplitude of the ion oscillator along the axial direction, as
described by the term

$\hat{H}_{F}=\frac{z_{\rm ax}F}{2}(\hat{a}^{{\dagger}}+\hat{a}).$ (2)

Here $z_{\rm ax}=\sqrt{\hbar/2m\omega_{z}}$ is the spread of the zero-point
wavefunction along the axial direction and $F$ is the parameter we wish to
estimate. The origin of the oscillating force can be a very weak electric
field, an optical dipole force, spin-dependent forces created in a magnetic-
field gradient or a Stark-shift gradient, etc. With the term (2) the total
Hamiltonian becomes

$\hat{H}_{\rm T}=\hat{H}_{\rm JC}+\hat{H}_{F}.$ (3)

In the following, we consider the weak-coupling regime $g\ll\omega$, in which
the phonon degree of freedom can be eliminated from the dynamics. This can be
carried out by applying the canonical transformation $\hat{U}=e^{\hat{S}}$ to
$\hat{H}_{\rm T}$ (3) such that $\hat{H}_{\rm eff}^{\rm
JC}=e^{-\hat{S}}\hat{H}_{\rm T}e^{\hat{S}}$ with
$\hat{S}=(g/\omega)(\sigma^{+}\hat{a}-\sigma^{-}\hat{a}^{{\dagger}})+(\Omega_{F}/g)(\hat{a}-\hat{a}^{{\dagger}})$.
Keeping only the terms of order of $g/\omega$ we arrive at the following
effective Hamiltonian (see the Appendix),

$\displaystyle\hat{H}_{\rm eff}^{\rm{JC}}$
$\displaystyle=\hbar\widetilde{\Delta}\sigma_{z}-\hbar\Omega_{F}\sigma_{x}-\hat{H}_{\rm
JC}^{\prime},$ (4a) $\displaystyle\hat{H}^{\prime}_{\rm JC}$
$\displaystyle=\frac{\hbar
g^{2}}{\omega}\sigma_{z}\hat{a}^{{\dagger}}\hat{a}.$ (4b)

This result indicates that the spin-motional interaction in Eq. (3) shifts the
effective spin frequency by the amount
$\widetilde{\Delta}=\Delta-g^{2}/2\omega$, while the effect of the force term
is to induce transitions between the spin states. The strength of the
transition is quantified by the Rabi frequency $\Omega_{F}=gz_{\rm
ax}F/2\hbar\omega$, which is proportional to the applied force $F$. Hence the
force estimation can be carried out by observing the coherent evolution of the
spin population that can be read out via state-dependent fluorescence.

Figure 1: (Color online) a) Time-evolution of the probability to find the
system in spin state $\left|\uparrow\right\rangle$ for the JC system. We
compare the probabilities derived from the original Hamiltonian (3) (dots) and
the effective Hamiltonian (4a) (solid lines). We assume an initial thermal
distribution with a mean phonon number $\bar{n}=1.2$. The parameters are set
to $g=4$ kHz, $\omega=170$ kHz, $\Delta=g^{2}/2\omega$, $z_{\rm ax}=14.5$ nm,
$F=20$ yN and $\Omega=10$ kHz. For the same initial state but in the absence
of driving field ($\Omega=0$), the signal loses contrast (blue dashed line).
b) Contrast of the Rabi oscillations defined as
$S=P_{\uparrow}(t_{2})-P_{\uparrow}(t_{1})$ with $t_{1}=\pi/2\Omega_{F}$ and
$t_{2}=\pi/\Omega_{F}$ with $\Omega_{F}=60$ kHz as a function of the mean
phonon number $\bar{n}$.

The last term $\hat{H}_{\rm JC}^{\prime}$ in Eq. (4a) is the residual spin-
motional coupling. This term affects the force estimation because it can be a
source of pure spin dephasing Viola1999 . Indeed, the $\sigma_{z}$ factor in
$\hat{H}_{\rm JC}^{\prime}$ induces transitions between the eigenstates
$\left|\pm\right\rangle$ of the operator $\sigma_{x}$ depending on the
vibrational state of the oscillator. As long as the oscillator is prepared
initially in an incoherent vibrational state at a finite temperature this
would lead to a random component in the spin energy. As we will see below, by
using dynamical decoupling the effect of the pure spin dephasing can be
reduced.

The sensing protocol starts by preparing the system in state
$\hat{\rho}(0)=\left|\uparrow\right\rangle\left\langle\uparrow\right|\otimes\hat{\rho}_{\rm
osc}$, where $\hat{\rho}_{\rm osc}$ stands for the initial density operator of
the oscillator. According to Eq. (4a), the evolution of the system is driven
by the unitary propagator $\hat{U}_{\rm JC}(t,0)=e^{-{\rm i}\hat{H}_{\rm
eff}^{\rm JC}t/\hbar}$. Assuming for the moment that $\hat{\rho}_{\rm
osc}=|0\rangle\langle 0|$ where $|n\rangle$ is the harmonic oscillator Fock
state with $n$ phonon excitations, the probability to find the system in state
$\left|\uparrow\right\rangle$ is $P_{\uparrow}(t)=\cos^{2}(\Omega_{F}t)$,
where for simplicity we set $\Delta=g^{2}/2\omega$, hence
$\widetilde{\Delta}=0$. In this case, the effect of $\hat{H}_{\rm
JC}^{\prime}$ automatically vanishes such that the signal exhibits a cosine
behavior according to the effective Hamiltonian (4a). An initial thermal
phonon distribution, however, would introduce dephasing on the spin
oscillations caused by thermal fluctuations. The spin coherence can be
protected, for example, by applying a sequence of fast pulses, which flip the
spin states and average the residual spin-motional interaction to zero during
the force estimation Yang2010 . On the other hand, because the relevant force
information is encoded in the $\sigma_{x}$ term in Eq. (4a), continuously
applying an additional strong driving field $\hat{H}_{\rm
d}=\hbar\Omega\sigma_{x}$ in the same basis Zheng ; Solano2003 , such that
$\hat{H}_{\rm T}\rightarrow\hat{H}_{\rm T}+\hat{H}_{\rm d}$, would not affects
the force estimation but rather will suppress the effect of the residual spin-
motional coupling. Indeed, going in the interaction frame with respect to
$\hat{H}_{\rm d}$, the residual spin-motional coupling becomes

$\hat{H}_{\rm JC}^{\prime}(t)=\frac{\hbar g^{2}}{\omega}(e^{2{\rm i}\Omega
t}|+\rangle\langle-|+e^{-2{\rm i}\Omega
t}|-\rangle\langle+|)\hat{a}^{{\dagger}}\hat{a}.$ (5)

The latter result indicates that the off-resonance transitions between states
$\left|\pm\right\rangle$ induced by $\hat{H}_{\rm JC}^{\prime}$ are suppressed
if $g^{2}/2\omega\ll\Omega$. By separating the pulse sequences from $t=0$ to
$t/2$ with a Hamiltonian $\hat{H}_{\rm T}+\hat{H}_{\rm d}$, and then from
$t/2$ to $t$ with a Hamiltonian $\hat{H}_{\rm T}-\hat{H}_{\rm d}$, the spin
states are protected from the thermal dephasing and the signal depends only on
the Rabi frequency $\Omega_{F}$ at the final time $t$. Note that the effect of
the magnetic field fluctuations of the spin states is described by an
additional $\sigma_{z}$ term in Eq. (4a), therefore the strong driving field
used here suppresses the spin dephasing caused by the magnetic-field
fluctuations, as was experimentally demonstrated Timoney2011 ; Webster2013 .

In Fig. 1(a) we show the time evolution of the probability $P_{\uparrow}(t)$
for an initial thermal vibrational state. Applying the driving field during
the force estimation leads to reduction of the spin dephasing and hence
protecting the contrast of the Rabi oscillations, see Fig. 1(b). We note that
a similar technique using a strong driving carrier field for dynamical
decoupling was proposed for the implementation of a high-fidelity phase gate
with two trapped ions Tan2013 ; Bermudez2012 .

Figure 2: (Color online) The sensitivity of the force measurement versus time
$t$ for various values of $\omega$. We assume an initial thermal vibrational
state with a mean phonon number $\bar{n}=1$. The solid lines represent the
analytical result given by Eq. (7) while the dots are the exact numerical
solution with the Hamiltonian (3) including the strong driving term. The other
parameters are set to $g=4$ kHz and $\Omega=7$ kHz.

The shot-noise-limited sensitivity for measuring $\Omega_{F}$ is

$\delta\Omega_{F}=\frac{\Delta P_{\uparrow}(t)}{\frac{\partial
P_{\uparrow}(t)}{\partial\Omega_{F}}\sqrt{\nu}},$ (6)

where $\Delta P_{\uparrow}(t)$ stands for the variance of the signal and
$\nu=T/\tau$ is the repetition number. Here $T$ is the total experimental
time, and the time $\tau$ includes the evolution time as well as the
preparation and measurement times. Because our technique relies on state-
projective detection, such that the preparation and measurement times are much
smaller than the other time scale, we assume $\tau\approx t$. From Eq. (6) we
find that the sensitivity, which characterizes the minimal force difference
that can be discriminated within a total experimental time of 1 s, is

$F_{\rm min}\sqrt{T}=\frac{\hbar\omega}{gz_{\rm ax}\sqrt{t}}.$ (7)

In Fig. 2 we show the sensitivity of the force estimation versus time $t$ for
different frequencies $\omega$ assuming an initial thermal vibrational state.
For an evolution time of $20$ ms, force sensitivity of $2$ yN $\rm{Hz}^{-1/2}$
can be achieved.

Let us now estimate the effect of the motional heating which limits the force
estimation. Indeed, the heating of the ion motion causes damping of the
signal, which leads to Wineland1998

$P_{\uparrow}(t)=\frac{1}{2}[1+e^{-\gamma t}\cos(2\Omega_{F}t)],$ (8)

where $\gamma$ is the decoherence rate. We assume that
$\gamma\sim\langle\dot{n}_{\rm ax}\rangle$ where $\langle\dot{n}_{\rm
ax}\rangle$ stands for the axial ion’s heating rate. Thus, the optimal force
sensitivity is Huelga1997

$F_{\rm min}\sqrt{T}=\frac{\hbar\omega}{gz_{\rm ax}}\sqrt{2\langle\dot{n}_{\rm
ax}\rangle e}.$ (9)

Using the parameters in Fig. 2 with $\omega=180$ kHz and assuming
$\langle\dot{n}_{\rm ax}\rangle=0.01$ $\rm{ms}^{-1}$ we estimate force
sensitivity of 2.4 yN $\rm{Hz}^{-1/2}$. For a cryogenic ion trap with heating
rate in the range of $\langle\dot{n}_{\rm ax}\rangle=1$ ${\rm s}^{-1}$ and
evolution time of $t=500$ ms, the force sensitivity would be 0.8 yN
$\rm{Hz}^{-1/2}$.

Figure 3: (Color online) Time-evolution of the probability to find the system
in spin state $\left|\uparrow\right\rangle$ for the QR system. We assume an
initial thermal vibrational state with a mean phonon number $\bar{n}=1.2$. Due
to the absence of residual spin-motion coupling the Rabi oscillations are
robust with respect to the spin dephasing caused by the thermal fluctuations.
We compare the probability derived from the Hamiltonian $\hat{H}_{\rm
T}=\hat{H}_{\rm QR}+\hat{H}_{F}$ with the analytical solution
$P_{\uparrow}=\cos^{2}(2\Omega_{F})$. The parameters are set to $g=4$ kHz,
$\omega=170$ kHz, $z_{\rm ax}=14.5$ nm, $F=20$ yN.

###  II.2 Quantum Rabi model

An alternative approach to sense the axial component of the force is to use a
probe described by the quantum Rabi model,

$\hat{H}_{\rm QR}=\hbar\omega\hat{a}^{{\dagger}}\hat{a}+\hbar
g\sigma_{x}(\hat{a}^{{\dagger}}+\hat{a}),$ (10)

which includes it the counter-rotating wave terms. This Hamiltonian can be
implemented by using a bichromatic laser field along the axial direction
Pedernales2015 . In the weak-coupling regime, $g\ll\omega$, we find by using
the unitary transformation $\hat{U}=e^{\hat{S}}$ with
$\hat{S}=-(g/\omega)\sigma_{x}(\hat{a}^{{\dagger}}-\hat{a})-(2\Omega_{F}/g)(\hat{a}^{{\dagger}}-\hat{a})$
that (see the Appendix)

$\hat{H}_{\rm eff}^{\rm QR}=-2\hbar\Omega_{F}\sigma_{x}.$ (11)

In contrast to Eq. (4a), now the effective Hamiltonian (11) does not contain
an additional residual spin-motional coupling, which implies that the spins
are immune to dephasing caused by the thermal motion fluctuations, see Fig. 3.
Thereby the force estimation can be carried out without additional strong
driving field. We find that the optimal force sensitivity is similar to Eq.
(9) but with extra factor of 2 in the denominator,

$F_{\rm min}\sqrt{T}=\frac{\hbar\omega}{2gz_{\rm
ax}}\sqrt{2\langle\dot{n}_{\rm ax}\rangle e}.$ (12)

Up to now we have considered probes that are responsive only to the axial
component of the force. In the following we propose a sensing technique that
can be used to detect the two transverse components of the time-varying
external force.

##  III Jahn-Teller quantum probe

In conventional ion trap sensing methods, the information on the force
direction can be extracted by using the three spatial vibrational modes of the
ion Maiwald2009 ; Munro2003 . Such an experiment requires an independent
measurement of the displacement amplitudes in each vibrational mode, which,
however, increases the complexity of the measurement procedure and can lead to
longer total experimental times. Here we show that by utilizing the laser-
induced coupling between the spin states and the transverse ion oscillation we
are able to detect the transverse components of the force by observing simply
the coherent evolution of the spin states.

Indeed, let us consider the case in which the small transverse oscillations of
the ion with a frequency $\omega_{\rm t}$ described by the Hamiltonian
$\hat{H}_{\rm t}=\hbar\omega_{\rm
t}(\hat{a}_{x}^{{\dagger}}\hat{a}_{x}+\hat{a}_{y}^{{\dagger}}\hat{a}_{y})$ are
coupled with the spin states via Jahn-Teller interaction. Such a coupling can
be achieved by using bihromatic laser fields with frequencies
$\omega_{b,r}=\omega_{0}\pm(\omega_{\rm t}-\omega)$ tuned respectively near
the blue- and red-sideband resonances, with a detuning $\omega$, which excite
the transverse $x$ and $y$ vibrational modes of the trapped ion. The
interaction Hamiltonian of the system is given by Porras2012 ; Ivanov2013

$\hat{H}_{\rm
JT}=\hbar\omega(\hat{a}^{{\dagger}}_{x}\hat{a}_{x}+\hat{a}^{{\dagger}}_{y}\hat{a}_{y})+\hbar
g\sigma_{x}(\hat{a}^{{\dagger}}_{x}+\hat{a}_{x})+\hbar
g\sigma_{y}(\hat{a}^{{\dagger}}_{y}+\hat{a}_{y}).$ (13)

Here $\hat{a}^{{\dagger}}_{\beta}$ and $\hat{a}_{\beta}$ are the creation and
annihilation operators of phonon excitations along the transverse direction
($\beta=x,y$) with an effective frequency $\omega$. The last two terms in Eq.
(13) describe the Jahn-Teller $E\otimes e$ spin-phonon interaction with a
coupling strengths $g$. In the following, we assume that a classical
oscillating force with a frequency $\omega_{d}=\omega_{\rm t}-\omega$
displaces the vibrational amplitudes along the transverse $x$ and $y$
directions of the quantum oscillator described by

$\hat{H}_{\vec{F}}=\frac{z_{\rm
t}F_{x}}{2}(\hat{a}^{{\dagger}}_{x}+\hat{a}_{x})+\frac{z_{\rm
t}F_{y}}{2}(\hat{a}^{{\dagger}}_{y}+\hat{a}_{y}),$ (14)

where $z_{\rm t}=\sqrt{\hbar/2m\omega_{\rm t}}$ is the size of the transverse
ion’s harmonic oscillator ground-state wavefunction. $F_{x}$ and $F_{y}$ are
the two transverse components of the force we wish to estimate. With the
perturbation term (14) the total Hamiltonian becomes

$\hat{H}_{\rm T}=\hat{H}_{\rm JT}+\hat{H}_{\vec{F}}.$ (15)

Assuming the weak-coupling regime, $g\ll\omega$, the two phonon modes are only
virtually excited. After performing the canonical transformation
$\hat{U}=e^{\hat{S}}$ of $\hat{H}_{\rm T}$ (15), where

$\hat{S}=(\hat{a}_{x}-\hat{a}_{x}^{{\dagger}})\left(\frac{g}{\omega}\sigma_{x}+\frac{\Omega_{x}}{g}\right)+(\hat{a}_{y}-\hat{a}_{y}^{{\dagger}})\left(\frac{g}{\omega}\sigma_{y}+\frac{\Omega_{y}}{g}\right),$
(16)

we obtain the following effective Hamiltonian (see the Appendix)

$\hat{H}_{\rm
eff}^{\rm{JT}}=-\hbar\Omega_{x}\sigma_{x}-\hbar\Omega_{y}\sigma_{y}+\hat{H}^{\prime}_{\rm
JT}.$ (17)

Here $\Omega_{x,y}=gz_{\rm t}F_{x,y}/\hbar\omega$ are the respective driving
Rabi frequencies of the transition between spin states
$\left|\uparrow\right\rangle$ and $\left|\downarrow\right\rangle$. The last
term in Eq. (17) is the residual spin-phonon interaction described by

$\hat{H}^{\prime}_{\rm JT}=2{\rm i}\frac{\hbar
g^{2}}{\omega}\sigma_{z}(\hat{a}_{x}^{{\dagger}}\hat{a}_{y}-\hat{a}_{x}\hat{a}_{y}^{{\dagger}}),$
(18)

which can be a source of thermal spin dephasing as long as the two phonon
modes are prepared in initial thermal vibrational states.

Figure 4: (Color online) a) Time-evolution of the probability to find the
system in spin state $\left|\uparrow\right\rangle$ for the JT system. We
compare the probability calculated from the Hamiltonian (15) assuming the
initial states
$\left|\psi(0)\right\rangle=\left|\uparrow\right\rangle\left|0_{x},0_{y}\right\rangle$
(red dots) and
$\left|\psi(0)\right\rangle=2^{-1/2}(\left|\uparrow\right\rangle+\left|\downarrow\right\rangle)|0_{x},0_{y}\rangle$
(blue triangles) with those given by the effective Hamiltonian (17) (solid
lines). The parameters are set to $g=4$ kHz, $\omega=170$ kHz, $z_{\rm t}=12$
nm, $F_{x}=20$ yN and $F_{y}=15$ yN. b) Oscillations of the signal for fixed
$t$ as a function of the phase $\phi$ for an initial superposition spin state.

The two-dimensional force sensing protocol starts by preparing the system in
state
$|\psi(0)\rangle=(c_{\uparrow}(0)\left|\uparrow\right\rangle+c_{\downarrow}(0)\left|\downarrow\right\rangle)\otimes\left|0_{x},0_{y}\right\rangle$,
where $c_{\uparrow,\downarrow}(0)$ are the respective initial spin probability
amplitudes and $\left|n_{x},n_{y}\right\rangle$ stands for the Fock state with
$n_{\beta}$ excitations in each phonon mode. According to the effective
Hamiltonian (17) the evolution of the system is driven by the free propagator
$\hat{U}_{\rm JT}=e^{-{\rm i}\hat{H}_{\rm eff}^{\rm JT}t/\hbar}$. Neglecting
the residual spin-motional coupling (18) the propagator reads

$\hat{U}_{\rm JT}^{0}(t,0)=\left[\begin{array}[]{cc}a&b\\\
-b^{*}&a^{*}\end{array}\right].$ (19)

Here $a=\cos(\widetilde{\Omega}t)$ and $b={\rm
i}e^{-\rm{i}\xi}\sin(\widetilde{\Omega}t)$ are the Cayley-Klein parameters,
which depend on the rms Rabi frequency $\widetilde{\Omega}=\frac{gz_{\rm
t}}{\hbar\omega}|\vec{F}_{\perp}|$, which is proportional to the magnitude of
the force $|\vec{F}_{\perp}|=\sqrt{F_{x}^{2}+F_{y}^{2}}$. In addition to
$|\vec{F}_{\perp}|$, we introduce the relative amplitude parameter
$\xi=\tan^{-1}\left(\frac{F_{y}}{F_{x}}\right)$. Assuming an initial state
with $c_{\uparrow}(0)=1$, $c_{\downarrow}(0)=0$, the respective probability to
find the system in state $\left|\uparrow\right\rangle$ is
$P_{\uparrow}(t)=\cos^{2}(\widetilde{\Omega}t)$, which implies that the Rabi
oscillations depends only on the magnitude of the force, see Fig. 4(a). Using
Eq. (6) we find that the shot-noise-limited sensitivity for measuring the
magnitude of the force is given by

$|\vec{F}_{\perp}|_{\rm min}\sqrt{T}=\frac{\hbar\omega}{2gz_{\rm t}\sqrt{t}}.$
(20)

In the presence of motional heating of both vibrational modes, the signal is
damped with decoherence rate
$\gamma\sim\langle\dot{n}_{x}\rangle+\langle\dot{n}_{y}\rangle$, where
$\langle\dot{n}_{\beta}\rangle$ is the heating rate along the $\beta$ spatial
direction. Therefore we find that the optimal force sensitivity is

$|\vec{F}_{\rm\perp}|_{\rm min}\sqrt{T}=\frac{\hbar\omega}{2gz_{\rm
t}}\sqrt{2(\langle\dot{n}_{x}\rangle+\langle\dot{n}_{y}\rangle)e}.$ (21)

It is important that due to the strong transverse confinement the sensing
scheme for measuring $|\vec{F}_{\perp}|$ is less sensitive to the ion’s
heating Turchette2000 ; Zhu2006 . Using the parameters in Fig. 4 and assuming
$\langle\dot{n}_{x}\rangle=\langle\dot{n}_{y}\rangle=1$ $\rm{s}^{-1}$ we
estimate force sensitivity of 0.6 yN $\rm{Hz}^{-1/2}$.

In order to detect the parameter $\xi$ we prepare the spin state in an initial
superposition state with $c_{\uparrow}(0)=1/\sqrt{2}$ and
$c_{\downarrow}(0)=e^{{\rm i}\phi}/\sqrt{2}$. Then the probability oscillates
with time as

$\displaystyle
P_{\uparrow}(t)=\frac{1}{2}\left[1+\sin(\xi-\phi)\sin(2\widetilde{\Omega}t)\right].$
(22)

Hence, for fixed evolution time $t$, the Ramsey oscillations versus the phase
$\phi$ provide a measure of the relative phase $\xi$, see Fig. 4(b).

In fact, Eq. (22) allows one to determine both the magnitude of the force
$|\widetilde{\Omega}|$ and the mixing parameter $\xi$ from the same signal
when plotted vs the evolution time $t$: $|\widetilde{\Omega}|$ is related to
the oscillation frequency and $\xi$ to the oscillation amplitude. The
parameter $\xi$ can be determined also by varying the externally controlled
superposition phase $\phi$, until the oscillation amplitude vanishes at some
value $\phi_{0}$; this signals the value $\xi=\phi_{0}$ (modulo $\pi$).

Finally, we discuss the dynamical decoupling schemes, which can be used to
suppress the effects of the term $\hat{H}_{\rm JT}^{\prime}$ (18) during the
force estimation. In that case, applying continuous driving field, e.g., along
the $\sigma_{x}$ direction, would reduce the thermal fluctuation induced by
$\hat{H}_{\rm JT}^{\prime}$, but additionally, the relevant force information,
which is encoded in the $\sigma_{y}$ term in (17), will be spoiled. Here we
propose an alternative dynamical decoupling scheme, which follows the Carr-
Purcell-Meiboom-Gill (CPMG) pulse sequence Carr1954 ; Meiboom1958 , in which,
however, the single instantaneous $\pi$ pulse is replaced by the phonon phase-
flip operator $\hat{R}_{\pi}=e^{{\rm
i}\pi\hat{a}_{x}^{{\dagger}}\hat{a}_{x}}$. Such a phonon phase shift
$\Delta\omega_{x}\tau=\pi$ can be achieved by switching the RF potential of
the trap by the fixed amount $\Delta\omega_{x}$ for a time $\tau$ Singer2010 .
The effect of $\hat{R}_{\pi}$ is to change the sign of the $\hat{H}_{\rm
JT}^{\prime}$ such that $\hat{R}^{{\dagger}}_{\pi}\hat{H}_{\rm
JT}^{\prime}\hat{R}_{\pi}=-\hat{H}_{\rm JT}^{\prime}$ but it leaves the other
part of the Hamiltonian (17) unaffected. Using that the pulse sequence
$\hat{U}_{1}=\hat{R}_{\pi}\hat{U}_{\rm JT}\hat{R}_{\pi}\hat{U}_{\rm JT}$
eliminates the residual spin-phonon coupling in the first order of the
interaction time $t$, a high-order reduction can be achieved by the recursion
$\hat{U}_{n}=\hat{R}_{\pi}\hat{U}_{n-1}\hat{R}_{\pi}\hat{U}_{n-1}$, which
eliminates the spin-phonon coupling up to $n$th order in $t$.

##  IV Summary and Outlook

We have proposed quantum sensing protocols, which rely on mapping the relevant
force information onto the spin degrees of freedom of the single trapped ion.
The force sensing is carried out by observing the Ramsey-type oscillations of
the spin states, which can be detected via state-dependent fluorescence. We
have considered quantum probes represented by the JC and QR systems, which can
be used to sense the axial component of the force. We have shown that when
using a JC system as a quantum probe, one can apply dynamical decoupling
schemes to suppress the effect of the spin dephasing during the force
estimation. When using a QR system as a probe, the absence of a residual spin-
phonon coupling makes the sensing protocol robust to thermally-induced spin
dephasing. Furthermore, we have shown that the transverse-force direction can
be measured by using a system described by the JT model, in which the spin
states are coupled with the two spatial phonon modes. Here the information of
the magnitude of the force and the relative ratio can be extracted by
observing the time evolution of the respective ion’s spin states, which
simplify significantly the experimental procedure.

Tuning the trap frequencies over the broad range, the force sensing methods
proposed here can be employed to implement a spectrum analyzer for ultra-low
voltages. Moreover, because in the force-field direction sensing the mutual
ratio can be additionally estimated our method can be used to implement a two-
channel vector spectrum analyzer. Finally, the realization of the proposed
force sensing protocols are not restricted only to trapped ions but could be
implemented with other quantum optical setups such as cavity-QED Dimer2007  or
circuit-QED systems Ballester2015 .

##  Appendix A Elimination of the vibrational degree of freedom

Let us make the canonical transformation of Hamiltonian
$\hat{H}=\hat{H}_{0}+\hat{H}_{\rm int}$,

$\displaystyle\hat{H}_{\rm eff}$ $\displaystyle=$ $\displaystyle
e^{-\hat{S}}\hat{H}e^{\hat{S}}=\hat{H}_{0}+\hat{H}_{\rm
int}+[\hat{H}_{0},\hat{S}]+[\hat{H}_{\rm int},\hat{S}]$ (23)
$\displaystyle+\tfrac{1}{2}[[\hat{H}_{0},\hat{S}],\hat{S}]+\tfrac{1}{2}[[\hat{H}_{\rm
int},\hat{S}],\hat{S}]+\ldots.$

Our goal is to choose $\hat{S}$ in a such a way that all terms of order $g$ in
$\hat{H}_{\rm eff}$ are canceled and the first term describing the spin-boson
interaction is of order $g^{2}/\omega$. If we determine $\hat{S}$ by the
condition

$\hat{H}_{\rm int}+[\hat{H}_{0},\hat{S}]=0,$ (24)

then the effective Hamiltonian becomes

$\hat{H}_{\rm eff}\approx\hat{H}_{0}+\tfrac{1}{2}[\hat{H}_{\rm int},\hat{S}].$
(25)

Let us consider the time-dependent operator $\hat{S}(t)=e^{{\rm
i}\hat{H}_{0}t/\hbar}\hat{S}e^{-{\rm i}\hat{H}_{0}t/\hbar}$, which obeys the
Heisenberg equation ${\rm i}\hbar\dot{\hat{S}}(t)=[\hat{S}(t),\hat{H}_{0}]$.
Using Eq. (24) we arrive at the equation

${\rm i}\hbar\dot{\hat{S}}(t)=\hat{H}_{\rm int}(t),$ (26)

where $\hat{H}_{\rm int}(t)=e^{{\rm i}\hat{H}_{0}t/\hbar}\hat{H}_{\rm
int}e^{-{\rm i}\hat{H}_{0}t/\hbar}$. Solving Eq. (26) we determine the desired
operator $\hat{S}$.

###  A.1 Jaynes-Cummings model

We identify $\hat{H}_{0}=\hbar\omega\hat{a}^{{\dagger}}\hat{a}$ and
$\hat{H}_{\rm int}=\hbar
g(\sigma^{-}\hat{a}^{{\dagger}}+\sigma^{+}\hat{a})+\frac{z_{\rm
ax}F}{2}(\hat{a}^{{\dagger}}+\hat{a})$. Using Eq. (26) we obtain

$\hat{S}=\frac{g}{\omega}(\sigma^{+}\hat{a}-\sigma^{-}\hat{a}^{{\dagger}})+\frac{z_{\rm
ax}F}{2\hbar\omega}(\hat{a}-\hat{a}^{{\dagger}}),$ (27)

which fulfills the condition (24). For the effective Hamiltonian we derive

$\displaystyle\hat{H}_{\rm eff}$ $\displaystyle=$
$\displaystyle\hbar\omega\hat{a}^{{\dagger}}\hat{a}+\hbar\left(\Delta-\frac{g^{2}}{2\omega}\right)\sigma_{z}-\hbar\Omega_{F}\sigma_{x}$
(28) $\displaystyle-\frac{\hbar
g^{2}}{\omega}\sigma_{z}\hat{a}^{{\dagger}}\hat{a}-\frac{\hbar
g^{2}}{2\omega}-\frac{z_{\rm ax}^{2}F^{2}}{4\hbar\omega}+\hat{H}^{\prime},$

where $\Omega_{F}=gz_{\rm ax}F/2\hbar\omega$ is the Rabi frequency and
$\hat{H}^{\prime}=\frac{1}{3}[[\hat{H}_{\rm int},\hat{S}],\hat{S}]+\ldots$
contains the higher-order terms in (23). We find

$\displaystyle\frac{1}{3}[[\hat{H}_{\rm int},\hat{S}],\hat{S}]$
$\displaystyle=$ $\displaystyle\frac{2g^{2}z_{\rm
ax}F}{3\omega^{2}}\sigma_{z}(\hat{a}^{{\dagger}}+\hat{a})-\frac{4\hbar
g^{3}}{3\omega^{2}}(\sigma^{-}\hat{a}^{{\dagger}}+\sigma^{+}\hat{a})$ (29)
$\displaystyle-\frac{4\hbar
g^{3}}{3\omega^{2}}(\sigma^{-}\hat{a}^{{\dagger}}\hat{a}^{{\dagger}}\hat{a}+\sigma^{+}\hat{a}^{{\dagger}}\hat{a}\hat{a}).$

As long as $g/\omega\ll 1$ the higher-order terms can be neglected and thus
the lowest-order effective Hamiltonian is given by Eq. (4a).

###  A.2 Quantum Rabi Model

Here the interaction Hamiltonian is $\hat{H}_{\rm int}=\hbar
g\sigma_{x}(\hat{a}^{{\dagger}}+\hat{a})+\frac{z_{\rm
ax}F}{2}(\hat{a}^{{\dagger}}+\hat{a})$ and the canonical transformation is
given by the operator

$\hat{S}=\frac{g}{\omega}\sigma_{x}(\hat{a}-\hat{a}^{{\dagger}})+\frac{z_{\rm
ax}F}{2\hbar\omega}(\hat{a}-\hat{a}^{{\dagger}}).$ (30)

The effective Hamiltonian is

$\hat{H}_{\rm
eff}=\hbar\omega\hat{a}^{{\dagger}}\hat{a}-2\hbar\Omega_{F}\sigma_{x}-\frac{\hbar
g^{2}}{\omega}-\frac{(z_{\rm ax}F)^{2}}{4\hbar\omega}.$ (31)

Remarkably, due to the equality $[[\hat{H}_{\rm int},\hat{S}],\hat{S}]=0$ all
higher-order terms in Eq. (23) vanish.

###  A.3 Jahn-Teller Model

Following the same procedure we have

$\displaystyle\hat{H}_{0}=\hbar\omega(\hat{a}_{x}^{{\dagger}}\hat{a}_{x}+\hat{a}_{y}^{{\dagger}}\hat{a}_{y}),$
$\displaystyle\hat{H}_{\rm int}=\hbar
g\sigma_{x}(\hat{a}_{x}^{{\dagger}}+\hat{a}_{x})+\hbar
g\sigma_{y}(\hat{a}_{y}+\hat{a}_{y})+\frac{z_{\rm
t}F_{x}}{2}(\hat{a}_{x}^{{\dagger}}+\hat{a}_{x})$ $\displaystyle+\frac{z_{\rm
t}F_{y}}{2}(\hat{a}^{{\dagger}}_{y}+\hat{a}_{y}).$ (32)

In this case the canonical transformation is represented by the operator

$\displaystyle\hat{S}$ $\displaystyle=$
$\displaystyle\frac{g}{\omega}\sigma_{x}(\hat{a}_{x}-\hat{a}_{x}^{{\dagger}})+\frac{g}{\omega}\sigma_{y}(\hat{a}_{y}-\hat{a}_{y}^{{\dagger}})+\frac{z_{\rm
t}F_{x}}{2\hbar\omega}(\hat{a}_{x}-\hat{a}_{x}^{{\dagger}})$ (33)
$\displaystyle+\frac{z_{\rm
t}F_{y}}{2\hbar\omega}(\hat{a}_{y}-\hat{a}_{y}^{{\dagger}}).$

Using Eq. (33) we obtain the following effective Hamiltonian

$\displaystyle\hat{H}_{\rm eff}$ $\displaystyle=$
$\displaystyle\hbar\omega(\hat{a}_{x}^{{\dagger}}\hat{a}_{x}+\hat{a}_{y}^{{\dagger}}\hat{a}_{y})-\hbar\Omega_{x}\sigma_{x}-\hbar\Omega_{y}\sigma_{y}-2{\rm
i}\frac{\hbar g^{2}}{\omega}$ (34)
$\displaystyle\times\sigma_{z}(\hat{a}_{x}\hat{a}_{y}^{{\dagger}}-\hat{a}_{x}^{{\dagger}}\hat{a}_{y})-\frac{2\hbar
g^{2}}{\omega}-\frac{z_{\rm
t}^{2}|\vec{F}_{\perp}|^{2}}{4\hbar\omega}+\hat{H}^{\prime},$

where $\Omega_{x,y}=gz_{\rm t}F_{x,y}/\hbar\omega$ are the respective Rabi
driving frequencies. The next higher-order terms in $\hat{H}^{\prime}$ (34)
are given by

$\displaystyle\frac{1}{3}[[\hat{H}_{\rm int},\hat{S}],\hat{S}]$
$\displaystyle=$ $\displaystyle 2{\rm i}\frac{g^{2}z_{\rm
t}F_{x}}{\omega^{2}}\sigma_{z}(\hat{a}_{y}^{{\dagger}}-\hat{a}_{y})-2{\rm
i}\frac{g^{2}z_{\rm
t}F_{y}}{\omega^{2}}\sigma_{z}(\hat{a}_{x}^{{\dagger}}-\hat{a}_{x})$ (35)
$\displaystyle-\frac{4\hbar
g^{3}}{\omega^{2}}\sigma_{y}\\{(\hat{a}_{y}^{{\dagger}}+\hat{a}_{y})(1+2\hat{n}_{x})-2\hat{a}_{x}^{{\dagger}2}\hat{a}_{y}$
$\displaystyle-2\hat{a}_{x}^{2}\hat{a}_{y}^{{\dagger}}\\}-\frac{4\hbar
g^{3}}{\omega^{2}}\sigma_{x}\\{(\hat{a}_{x}^{{\dagger}}+\hat{a}_{x})(1+2\hat{n}_{y})$
$\displaystyle-2\hat{a}_{y}^{{\dagger}2}\hat{a}_{x}-2\hat{a}_{y}^{2}\hat{a}_{x}^{{\dagger}}\\}.$

## References

  * (1) P. Treutlein, C. Genes, K. Hammerer, M. Poggio and P. Rabl, _Cavity Optomechanics_ ed. M. Aspelmeyer, T. Kippenberg and F. Marquardt (Berlin: Springer).
  * (2) G. Kurizki, P. Bertet, Y. Kubo, K. Molmer, D. Petrosyan, P. Rabl, J. Schmiedmayer, Proc. Natl. Acad. Sci. USA 112, 3866 (2015)
  * (3) T. D. Stowe, K. Yasumura, T. W. Kenny, D. Botkin, K. Wago and D. Rugar, Appl. Phys. Lett. 71, 288 (1997).
  * (4) A. A. Geraci, S. J. Smullin, D. M. Weld, J. Chiaverini, and A. Kapitulnik, Phys. Rev. D 78, 022002 (2008).
  * (5) D. Rugar, R. Budakian, H. Mamin, B. Chui, Nature 430, 329 (2004).
  * (6) S. Kolkowitz, A. C. B. Jayich, Q. P. Unterreithmeier, S. D. Bennett, P. Rabl, J. G. E. Harris, M. D. Lukin, Science 335, 1603 (2012).
  * (7) R. Blatt and D. Wineland, Nature 453, 1008 (2008).
  * (8) M. J. Biercuk, H. Uys, J. W. Britton, A. P. VanDevender and J. J. Bollinger, Nat. Nanotechnol. 5, 646 (2010).
  * (9) S. Knünz, M. Herrmann, V. Batteiger, G. Saathoff, T. W. Hänsch, K. Vahala, and Th. Udem, Phys. Rev. Lett. 105, 013004 (2010).
  * (10) R. Maiwald, D. Leibfried, J. Britton, J. C. Bergquist, G. Leuchs and D. J. Wineland, Nat. Phys. 5, 551 (2009).
  * (11) P. A. Ivanov, K. Singer, N. V. Vitanov and D. Porras, Phys. Rev. Appl. 4, 054007 (2015).
  * (12) Shi-Biao Zheng, Phys. Rev. A 66, 060303(R) (2002).
  * (13) E. Solano, G. S. Agarwal, and H. Walther, Phys. Rev. Lett. 90, 027903 (2003).
  * (14) D. J. Wineland, C. Monroe, W. M. Itano, D. Leibfried, B. E. King, and D. M. Meekhof, J. Res. Natl. Inst. Stand. Technol. 103, 259 (1998).
  * (15) H. Häffner, C. F. Roos, and R. Blatt, Phys. Rep. 469, 155 (2008).
  * (16) C. Schneider, D. Porras, and T. Schaetz, Rep. Prog. 75, 024401 (2012)
  * (17) L. Viola, E. Knill, and S. Lloyd, Phys. Rev. Lett. 82, 2417 (1999).
  * (18) W. Yang, Z.-Y. Wang, and R.-B. Liu, Front. Phys. 6, 1 (2010).
  * (19) N. Timoney, I. Baumgart, M. Johanning, A. F. Varon. M. B. Plenio, A. Retzker, and Ch. Wunderlich, Nature 476, 185 (2011).
  * (20) S. C. Webster, S. Weidt, K. Lake, J. J. McLoughlin, and W. K. Hensinger, Phys. Rev. Lett. 111, 140501 (2013).
  * (21) T. R. Tan, J. P. Gaebler, R. Bowler, Y. Lin, J. D. Jost, D. Leibfried, and D. J. Wineland, Phys. Rev. Lett. 110, 263002 (2013).
  * (22) A. Bermudez, P. O. Schmidt, M. B. Plenio, and A. Retzker, Phys. Rev. A 85, 040302(R) (2012).
  * (23) S. F. Huelga, C. Macchiavello, T. Pellizzari, A. K. Ekert, M. B. Plenio, J. I. Cirac, Phys. Rev. Lett. 79, 3865 (1997).
  * (24) J. S. Pedernales, I. Lizuain, S. Felicetti, G. Romero, L. Lamata and E. Solano, Sci. Rep. 5, 15472 (2015).
  * (25) W. J. Munro, K. Nemoto, G. J. Milburn, and S. L. Braunstein, Phys. Rev. A 66, 023819 (2002).
  * (26) D. Porras, P. A. Ivanov, and F. Schmidt-Kaler, Phys. Rev. Lett. 108, 235701 (2012).
  * (27) P. A. Ivanov, D. Porras, S. S. Ivanov and F. Schmidt-Kaler, J. Phys. B: At. Mol. Opt. Phys. 46, 104003 (2013).
  * (28) Q. A. Turchette, D. Kielpinski, B. E. King, D. Leibfried, D. M. Meekhof, C. J. Myatt, M. A. Rowe, C. A. Sackett, C. S. Wood,cW. M. Itano, C. Monroe, and D. J. Wineland, Phys. Rev. A 61, 063418 (2000).
  * (29) Shi-Liang Zhu, C. Monroe, and K.-M. Duan, Phys. Rev. Lett. 97, 050505 (2006).
  * (30) H. Carr and E. M. Purcell, Phys. Rev. 94, 630 (1954).
  * (31) S. Meiboom and D. Gill, Rev. Sci. Instrum. 29, 688 (1958).
  * (32) K. Singer, U. Poschinger, M. Murphy, P. Ivanov, F. Ziesel, T. Calarco, F. Schmidt-Kaler, Rev. Mod. Phys. 82, 2609 (2010).
  * (33) F. Dimer, B. Estienne, A. S. Parkins, H. J. Carmichael, Phys. Rev. A 75, 013804 (2007).
  * (34) D. Ballester, G. Romero, J. J. Garcia-Ripoll, F. Deppe, E. Solano, Phys. Rev. X 2, 021007 (2012).