Optimal estimation of the optomechanical coupling strength
József Zsolt Bernád Claudio Sanavio André Xuereb
quant-ph





zsolt.bernad@um.edu.mt
Department of Physics, University of Malta, Msida MSD 2080, Malta

Department of Physics, University of Malta, Msida MSD 2080, Malta

Department of Physics, University of Malta, Msida MSD 2080, Malta




We apply the formalism of quantum estimation theory to obtain information about the value of the nonlinear optomechanical coupling strength. 
In particular, we discuss the minimum mean-square error estimator and a quantum Cramér–Rao-type inequality for the estimation of the coupling strength. 
Our estimation strategy reveals some cases where quantum statistical inference is inconclusive and merely result in the reinforcement of prior expectations. 
We show that these situations also involve the highest expected information losses. We demonstrate that interaction times in the order of one time period of mechanical oscillations 
are the most suitable for our estimation scenario, and compare situations involving different photon and phonon excitations.


Optimal estimation of the optomechanical coupling strength
    André Xuereb
    December 30, 2023
==========================================================





§ INTRODUCTION

Quantum estimation theory attempts to find the best strategy for learning the value of one or more parameters of the density matrix of a quantum mechanical system <cit.>, 
expressed as a positive-operator valued measure (POVM). An estimation protocol consists in identifying this POVM, applying its elements in repeated measurements of the system, and finally 
estimating the unknown parameters from the data set. The optimum strategy consists of those POVMs which minimize an average cost functional, typically considered for the maximum likelihood or 
mean-square error estimators. The mathematical framework for studying the conditions under which solutions of the optimization problem exists was established by Holevo <cit.>. 
Subsequently, considerable theoretical and experimental developments in quantum statistical inference have led to various applications in quantum tomography and metrology <cit.>.

We explore the implementation of this methodology to the specific case of an optomechanical system. It has been known for a long time that electromagnetic radiation exerts “radiation pressure” on 
any surface exposed to it <cit.>. The resulting momentum transfer is particularly notable on mirrors forming a cavity for the electromagnetic radiation. This subject has been brought 
to the forefront in recent years because of its impact on the design and operation on laser-based interferometric gravitational wave observatories <cit.>, where it imposes limits on the 
continuous detection of the mirror positions <cit.>. Under the guise of optomechanics, it has also drawn much attention both theoretically and experimentally <cit.>, 
motivated by its applications in sensitive optical detection of weak forces, mechanical motion in the quantum regime, and coherent light–matter interfaces that could form the backbone of future 
quantum information devices.

The simplest optomechanical interaction is deceptively simple to describe. Starting from two uncoupled harmonic oscillators, one representing a single mode of the electromagnetic field and one representing 
the mechanical 
oscillator, the only other free parameter in this model is the coupling strength that quantifies the extent by which the mirror rest position moves when a single photon is added to the electromagnetic field. 
The aim of this paper is to obtain the optimal strategy for measuring the value of this coupling strength. Our approach is based on quantum inference techniques, which have been successfully applied to phase 
estimations of quantum states <cit.>. In our case the parameter to be estimated is not a simple phase parameter, but rather a parameter which appears both in the 
spectrum <cit.> and the eigenvalues of the quantum state. If the motion of the mechanical oscillator is adiabatically slow compared to the frequency separation of the optical 
modes <cit.>, an analytical solution can be found for the evolution of a system described by this model <cit.>. The resulting density matrix describes the joint state of the 
field and the mechanical oscillator. Since measurements are, in practice, usually performed on the state of the electromagnetic field emerging from the system, we trace out the mechanical degrees of freedom. 
We therefore concentrate on the resulting optical state, which will be subject of our quantum estimation procedure.

In this paper we shall focus on a mean-square error estimator and assume that the prior probability density function of the coupling constant is a Gaussian distribution, which will allow us to keep our 
calculations analytic as far as possible. We set the mean and the standard deviation of this distribution function to values obtained by a canonical quantization procedure with a high frequency cut-off of 
the radiation field and adiabatically slow motion of the mechanical oscillator. In order to illustrate the basic features of our proposal, we consider the mechanical oscillator to be initially in 
(i) a coherent state, (ii) a thermal state, and (iii) a squeezed state. For the sake of simplicity and to keep our numerical calculations tractable, we will assume that the initial optical field state has 
only a few excitations. We determine the mean-square error estimator which minimizes the cost functional and study a quantum Cramér–Rao-type inequality for the mean-squared error of a biased estimator. 
For our model, a right logarithmic or a symmetrized logarithmic derivative operator of the electromagnetic field state with respect to the coupling constant is very challenging to construct since the Hilbert space 
of the harmonic oscillator is infinite-dimensional. Therefore, we explore the possibilities of deriving an analytically accessible lower bound formula of the mean-squared error of the estimator. Obviously, 
we will employ the standard techniques for the derivation of the unbiased quantum Cramér–Rao inequality <cit.>.

This paper is organized as follows. In model we discuss the model and its solutions of the single mode radiation coupled via radiation pressure to a vibrational mode of a mechanical oscillator. 
In mst we introduce quantum estimation theory for minimizing mean-square error estimators, and study its properties when applied to the optomechanical model. We then address the mean-squared error of the 
biased estimator, in qCR, deriving a lower bound and employing this result to the optomechanical model. A discussion about our analytical and numerical findings is the subject of a summary in 
conclusion.



§ MODEL


We consider a system composed of two harmonic oscillators, a single mode of the radiation field and a vibrational mode of a mechanical oscillator. Provided that the electromagnetic field 
describes a high-finesse cavity field mode, and that one of the mirrors is movable, one is able to derive a radiation-pressure interaction Hamiltonian <cit.> by using time-varying 
boundary conditions in the quantization procedure, resulting in (ħ=1)

    Ĥ=ω_câ^†â + ω_mb̂^†b̂+ g â^†â (b̂^†+ b̂),

where â (â^†) is the annihilation (creation) operator of the single mode electromagnetic field, with frequency ω_c, and b̂ (b̂^†) is the 
annihilation (creation) operator of the mirror motion, with frequency ω_m. The strength of the optomechanical interaction, denoted g, depends on the specifics of the realization in question.

Starting from a joint field–mechanics state |Ψ(0)⟩, the time evolution of the system is given by the Schrödinger equation and can be rephrased as

    |Ψ(t)⟩=e^-i Ĥ t|Ψ(0)⟩.

We are interested in the case where there are no initial correlations between the field and the mechanical oscillator. Therefore, we choose an initial state of the form

    |Ψ(0)⟩=∑^∞_n=0 a_n |n⟩_c|ψ⟩_m,

where the exact form of |ψ⟩_m depends on the desired initial conditions. In the following subsections we consider initial coherent, thermal, and squeezed states.



 §.§ Coherent state

We set
We start off by setting the initial mechanical oscillator to a coherent state <cit.>

    |ψ⟩_m=|α⟩_m=∑_n=0^∞
     e^-|α|^2/2α^n/√(n!)|n⟩_m, with α=|α| e^iϕ,

which we write in terms of the field number states |n⟩_m (n∈ℕ_0). Here, |α| is the amplitude of the coherent state and ϕ its phase. 
We allow the coefficients a_n of the photon-number states to be general and only impose the normalization condition ∑_n|a_n|^2=1. The choice of the initial state in cohinitial 
is our basic approach for determining the time evolution of the system, and will eventually be extended to cover initial thermal and squeezed states of the mechanical oscillator.

The interaction Hamiltonian g â^†â (b̂^†+ b̂) commutes with the free Hamiltonian of the radiation field, ω_câ^†â, which yields

    _c⟨n|Ĥ|m⟩_c=[n ω_cÎ + ω_mb̂^†b̂+ n g (b̂^†+ b̂)] δ_n,m

with δ_n,m being the Kronecker delta and Î the identity operator on the Hilbert space of the mechanical oscillator. Thus, the Hamiltonian (<ref>) is block-diagonal with 
respect to photon-number states |n⟩_c.

In order to evaluate the expression exp{-i ω_mb̂^†b̂t-i n g (b̂^†+ b̂)t} we employ the Baker–Campbell–Hausdorff formula and obtain 
(see, for example, Ref. <cit.>)

    e^-i ω_mb̂^†b̂ t - i n g (b̂^†+ b̂)t e^i ω_mb̂^†b̂ t=
     e^iΦ_n(t) e^α_n(t)b̂^†-α^∗_n(t)b̂,

where we have introduced the parameters

    α_n(t)   =n g/ω_m(e^-iω_m t-1) and
    Φ_n(t)   =n^2 g^2/ω^2_m[ω_m t-sin(ω_m t)].

This implies that the full time evolution can be viewed as photon-number dependent displacements of the mechanical oscillator; with the help of calc1 and (<ref>), we find


    |Ψ(t)⟩   =∑^∞_n=0 a_n e^i φ_n(t)|n⟩_c|β_n(t)⟩_m,
    φ_n(t)   =-n ω_c t + n^2 g^2/ω^2_m[ω_m t-sin(ω_m t)]
           +ng/ω_mα^∗(1-e^iω_m t)-α(1-e^-iω_m t) /2i, and
    β_n(t)   =n g/ω_m(e^-iω_m t-1)+ α e^-i ω_m t,



where we have also used a corollary of the Baker–Campbell–Hausdorff formula, which states that the product of two displacement operators is also a displacement operator with an overall phase factor.

The quantum state of mainr yields a complete description of the interaction between the single mode of the radiation field and the single vibration mode of the mechanical oscillator, i.e., 
neglecting all losses and sources of decoherence. In the subsequent sections we will be interested in possible measurement scenarios of the field, which are capable to estimate the couple constant g. 
Therefore, the field to be measured for estimation reads

    ρ̂_F=_mech{|Ψ(t)⟩⟨Ψ(t)|}=∑^∞_n,m=0 A_n,m|n⟩_c⟨m|

with

    A_n,m   =a_n a^∗_m e^i φ_n(t)-iφ_m(t)-(|β_n(t)|^2+|β_m(t)|^2)/2+β_n(t)β^∗_m(t)
       =a_n a^∗_m e^-g^2 f^(2)_n,m(t)+g f^(1)_n,m(t)-f^(0)_n,m(t),

where


    f^(0)_n,m(t)   =i ω_c t (n-m),
    
    f^(1)_n,m(t)   =α^∗(1-e^iω_m t)-α(1-e^-iω_m t) /ω_m(n-m), and
    
    f^(2)_n,m(t)   =1-cos(ω_m t)/ω^2_m (n-m)^2
           -i ω_m t-sin(ω_m t)/ω^2_m (n^2-m^2).



We note in passing that these results imply that the coefficient of the linear term in g contributes to A_n,m only when the initial state of the mechanical oscillator is not in the ground state, 
i.e., α≠0. Before moving on to developing the quantum error theory we require, let us briefly extend our considerations to thermal and squeezed states.



 §.§ Thermal state


If the initial state of the mechanical oscillator is thermal, as a first step we must switch from discussing state vectors to density matrices. In this case, the uncorrelated initial state of the 
optomechanical system has the form

    ρ̂(t=0)=∑^∞_n,m=0 a_n a^∗_m∫e^-|γ|^2/n_th/π n_th|n⟩_c|γ⟩_m⟨m|_c⟨γ|_m ^2 γ,

where we have used the Glauber–Sudarshan representation <cit.> of the mechanical oscillator thermal state with the average phonon number

    n_th = [exp( ħω_m/k_BT)-1 ]^-1,

where k_B is the Boltzmann constant and T the thermodynamic temperature of the initial state of the mechanical system. The time evolution of this system is given by

    ρ̂(t)=e^-i Ĥ tρ̂(0) e^i Ĥ t,

which yields

    ρ̂(t)=∑^∞_n,m=0 a_n a^∗_m∫e^-|γ|^2/n_th/π n_th e^i φ_n(t)-i φ_m(t)|n⟩_c⟨m|_c
    ×|β_n(t)⟩_m⟨β_m(t)|_m ^2 γ,

with

    φ_n(t)   =-n ω_c t + n^2 g^2/ω^2_m[ω_m t-sin(ω_m t)]
           +ng/ω_mγ^∗(1-e^iω_m t)-γ(1-e^-iω_m t) /2i, and
    β_n(t)   =n g/ω_m(e^-iω_m t-1)+ γ e^-i ω_m t.

In the next step we trace out the mechanical degrees of freedom, as before, obtaining

    ρ̂_F=_mech{ρ̂(t)}=∑^∞_n,m=0 a_n a^∗_m|n⟩_c⟨m|_c
    ×∫e^-|γ|^2/n_th/π n_th B_n,m(γ,γ^∗)  ^2γ,

with

    B_n,m(γ,γ^∗)= e^-h_0+h_1,

where

    h_0   =i ω_c t (n-m)+g^2/ω^2_m{[1-cos(ω_m t)](n-m)^2
           -i [ ω_m t-sin(ω_m t) ] (n^2-m^2) }, and
    
    h_1   =g/ω_m[γ^∗(1-e^iω_m t)-γ(1-e^-iω_m t) ] (n-m).

Now, we perform the Gaussian integral by using ^2γ={γ} {γ} and obtain a density matrix in the form of densitym. Employing the notation of Aform, we find

    f^(0)_n,m(t)   =i ω_c t (n-m),
    
    f^(1)_n,m(t)   =0, and
    
    f^(2)_n,m(t)   =(2 n_th+1)1-cos(ω_m t)/ω^2_m (n-m)^2
           -i ω_m t-sin(ω_m t)/ω^2_m (n^2-m^2).




 §.§ Squeezed state


Let us now consider the case where the initial state of the mechanical system is a displaced squeezed state; we write, therefore,

    |Ψ(0)⟩=∑^∞_n=0 a_n |n⟩_c|α,ζ⟩_m

with the mechanical oscillator state being defined as <cit.>

    |α,ζ⟩_m=D̂(α)Ŝ(ζ)|0⟩_m

where D̂(α)=exp(αb̂^†-α^∗b̂), with α=|α|e^iϕ, is the displacement operator, and Ŝ(ζ)=exp[12(ζ^∗b^2-ζb^† 2)], with ζ=|ζ|e^iθ, is the squeezing operator.

We employ the squeezed state of squeezed; in passing, however, we note that it is possible to invert the order of the displacement and squeezing operator. This
results in a generalized squeezed state, which differs from the original state by the displacement parameter:

    Ŝ(ζ)D̂(α)=D̂[αcosh(|ζ|)-α^∗ e^iθsinh(|ζ|)]Ŝ(ζ).

Exploiting the block-diagonal structure of the Hamiltonian with respect to the photon-number states |n⟩_c, and calc2, we find

    _c⟨n| e^-iĤt|Ψ(0)⟩=a_n e^i φ_n(t)|β_n(t),ζe^-2iω_m t⟩_m

where φ_n(t) and β_n(t) are defined in mainr. Next, tracing out the mechanical degrees of freedom yields the state of the field in the form of densitym, with

    A_n,m=a_n a^∗_m e^i φ_n(t)-iφ_m(t)
    ×{|β_n(t),ζe^-2iω_m t⟩⟨β_m(t),ζe^-2iω_m t|}.


The trace in this equation can be evaluated with the help of the Glauber–Sudarshan representation, which allows us to write

    {|β_n(t),ζe^-2iω_m t⟩⟨β_m(t),ζe^-2iω_m t|}=∫^2γ/π⟨γ|β_n(t),ζe^-2iω_m t⟩⟨β_m(t),ζe^-2iω_m t|γ⟩.

First, we note that

    ⟨γ|β_n(t),ζe^-2iω_m t⟩_m=e^-(γβ_n^∗(t)-γ^∗β_n(t))/2⟨0|D̂^†(γ-β_n(t))Ŝ(ζ
    e^-2iω_m t)|0⟩,

where we have used the relation D̂^†(-γ)=D̂(γ). The overlap integral between the coherent state |γ⟩ and the squeezed state |0,ζ⟩ is

    ⟨γ|0,ζ⟩=√(e^-|γ|^2/cosh(|ζ|))exp[-12γ^∗2e^iθtanh(|ζ|)].

For the purposes of todosq we thus obtain

    ⟨γ|β_n(t),ζe^-2iω_m t⟩_m=√(e^-|γ-β_n(t)|^2-γβ_n^∗(t)+γ^∗β_n(t)/cosh(|ζ|))exp{-12[γ^∗-β^∗_n(t)]^2e^i(θ-2 ω_m t)tanh(|ζ|)}.

Substituting this result into tracetodosq and performing the integral by using ^2γ={γ} {γ}
we obtain the coefficients in Aform:


    f^(0)_n,m   =i ω_c t(n-m)+|α|^2[ 1+tanh(|ζ|)cos(θ-2ϕ)]-I^(0)+ln[cosh(|ζ|)√(1-tanh^2(|ζ|))],
    
    f^(1)_n,m(t)   =α^∗(1-e^iω_m t)/ω_m(n-m)+I^(1)_n,m/ω_m+tanh(|ζ|)α^∗(1-e^-iω_m t)
    e^iθn+α(1-e^iω_m t)e^-iθm/2ω_m
    
    f^(2)_n,m(t)   =-i ω_m t-sin(ω_m t)/ω^2_m (n^2-m^2)+tanh(|ζ|)(e^-iω_m t-1)^2e^iθn^2+(e^iω_m t-1)^2
    e^-iθm^2/2ω^2_m
           +1-cosω_m t/ω^2_m (n^2+m^2)+I^(2)_n,m/ω^2_m,



where, for simplicity of presentation, we relegate the explicit form of the coefficients I^(0), I^(1)_n,m, and I^(2)_n,m to Calc.




§ QUANTUM MINIMUM MEAN-SQUARE ERROR ESTIMATION


Quantum estimation theory attempts to find the best strategy for estimating one or more parameters of the density matrix <cit.>. In our case, the density matrix of the field in densitym 
depends on the parameter g to be estimated. Any outcome of a measurement on the field is a variable with probability depending on the estimanda g, the parameter to be estimated. 
As our knowledge of g is limited, we assume that the 
estimanda is a random variable with prior probability density function

    p(g)=1/√(2πσ^2) e^-(g-g_0)^2/2σ^2

with mean g_0 and variance σ^2. We shall return to these parameters and their physical meanings later on.

Our estimation problem is to thus find the best measurements on ρ̂_F(g) to estimate g. In practice, we are looking for a POVM the elements Π̂(Δ) of which are defined on a compact 
interval Δ⊂ℝ of the real line, i.e., the set of all possible values for g, and satisfy


    0 ⩽Π̂(Δ) ⩽Î,

where Î is the identity operator. We also suppose that the infinitesimal operators Π̂(g) can be formed, thus yielding

    Π̂(Δ)=∫_ΔΠ̂(g),

and

    Î=∫^∞_-∞Π̂(g).



In order to solve the estimation problem we have to also provide a cost function, a measure of the cost incurred upon making errors in the estimate of g. Here, we wish to minimize the average squared cost of 
error, which is encoded in the cost function

    C(g̃, g)=(g̃ -g)^2,

where g̃ is the estimate of g, and therefore a function of the measurement data.

Now, we are able to formulate the quantum estimation problem. We are looking for Π̂(g̃), which minimizes the average cost of this estimation strategy

    C̅={∫^∞_-∞∫^∞_-∞ p(g) C(g̃, g) ρ̂_F(g)Π̂(g̃) g },

under the constraints embodied in cons. This is a variational problem for the functional C̅. To proceed, we consider each possible estimate g̃ to be an eigenvalue of the Hermitian operator

    M̂=∫^∞_-∞g̃ Π̂(g̃) = ∫^∞_-∞g̃|g̃⟩⟨g̃| g̃

with eigenstates |g̃⟩. Thus, the average cost functional in averagecost, calculated using the cost function in costf, can be written as

    C̅ [M̂]={∫^∞_-∞ p(g) (M̂-g Î)^2 ρ̂_F(g) g }.

For convenience, we now define the following operators (k=0,1,2):

    Γ̂_k= ∫^∞_-∞ g^k p(g) ρ̂_F(g) g.

Now, let ϵ be a real number and N̂ any Hermitian operator. Let M̂_min be the Hermitian operator which minimizes C̅ [M̂]. Then, we have

    C̅ [M̂_min] ⩽C̅ [M̂_min+ϵN̂],

because the sum of Hermitian operators is itself a Hermitian operator. Evaluating the right-hand side of the inequality and using the operators defined in Gamma, we obtain

    C̅ [M̂_min+ϵN̂]=C̅ [M̂_min]
    
    +ϵ {N̂(Γ̂_0 M̂_min +
    M̂_minΓ̂_0 -2Γ̂_1 )}
    
    +ϵ^2 {Γ̂_0 N̂^2 }.

By differentiating this relation with respect to ϵ and equating the result to zero, one is able to show that the unique Hermitian operator M̂_min minimizing C̅ must 
satisfy <cit.>

    Γ̂_0 M̂_min +M̂_minΓ̂_0=2Γ̂_1.

The average minimum cost of error for this measurement is

    C̅_min   ={Γ̂_0 M̂^2_min-2 Γ̂_1 M̂_min + Γ̂_2}
       ={Γ̂_2-M̂_minΓ̂_0 M̂_min},

where we have used the relation in tosolve to simplify the result. In order to determine M̂_min we thus have to solve the operator equation tosolve. It has been shown in
Ref. <cit.> that the unique solution of this equation can be written as

    M̂_min=2 ∫^∞_0 exp(-Γ̂_0 x) Γ̂_1 exp(-Γ̂_0 x)  x.

A comment about this solution is in order. The operator that we have found does not necessarily represent the best estimator of g, but rather the measurement operator which protects best 
against information loss, no matter what the true value of g is. Further discussion about this subtlety and its relation to biased estimators can be found in the illuminating monograph by Jaynes <cit.>.

We can evaluate all the Γ̂_k by using the form of ρ̂_F(g) in densitym, thus obtaining (k=0,1,2)

    Γ̂_k=∑^∞_n,m=0 a_n a^∗_m A^(k)_n,mexp(-γ_n,m) |n⟩_c⟨m|,

with


    A^(0)_n,m   = 1/σ',
    
    A^(1)_n,m   = g_0+f^(1)_n,m(t) σ^2/σ'^3, and
    
    A^(2)_n,m   = (g_0+f^(1)_n,m(t) σ^2)^2+σ^2 σ'^2/σ'^5,


where we have also introduced

    γ_n,m   ={2 g^2_0 f^(2)_n,m(t)- 2g_0 f^(1)_n,m(t)+2f^(0)_n,m(t)σ'^2
           - [f^(1)_n,m(t)]^2σ^2}/(2σ'^2), and
    σ'^2   =2f^(2)_n,m(t) σ^2 +1.


Written in this form, our results are very general. In the following section we will investigate some simple cases, related to the optomechanical model introduced previously.



§ A CASE STUDY: OPTOMECHANICS


We shall now put the results obtained in the two preceding sections together, thus allowing us to study the process of quantum mean-square error estimation as it applies to an optomechanical system.

The simplest non-trivial case results when a_n=0 for n>1 in initial. In order to maximize the absolute values of the off-diagonal elements of the density matrix we choose a_0=a_1=1/√(2). 
This specific choice is due to the fact that the unknown parameter g is only present in the off-diagonal elements, as can be seen from example1. Here, the estimate g̃ is simply one of the 
two eigenvalues of M̂_min, which turn up as a result of applying the two projective measurements defined by their accompanied eigenvectors.

Furthermore, we ought to define g_0 and σ in apriori, the prior probability density function of the estimanda g. We set


    g_0   =ω_c/L√(⟨x̂^2 ⟩_0)=ω_c/L√(ħ/2 m ω_m), and
    σ^2   =(ω_c/L)^2 √(⟨x̂^4 ⟩_0-⟨x̂^2 ⟩^2_0)=(ω_c/L)^2
     ħ/√(2) m ω_m,



where L is the length of the cavity, m is the mass of the mechanical oscillator, and ⟨Â⟩_0 is the expectation value of operator Â, acting only on the Hilbert space of 
the mechanical oscillator, in the ground state <cit.>. For the sake of simplicity we perform our calculations in the rotating frame of the single-mode field, i.e., ρ̂_F→Ûρ̂_FÛ^† with Û=exp{-iω_c t  â^†â}.



 §.§ Coherent state


We determine M̂_min from the Γ̂_k in gammak by using intx. One can obtain analytical results; however, due to their complex structure we omit their explicit 
presentation here. Instead, we focus on numerical solutions. First, we investigate the average minimum cost of error C̅_min; costfig shows C̅_min as a function of time, 
which decreases until it reaches its minimum and then returns asymptotically to its initial value, which is equal to σ^2. At t=0, where no interaction occurred, the eigenvalues of M̂_min 
are g_0 and zero. The probability of
measuring the eigenvalue zero is zero and therefore the estimate is g_0. It is immediate from the form of the prior probability distribution p(g) in apriori that the average minimum cost of error 
is σ^2, or simply the variance of p(g), at t=0.



In the opposite limit, t →∞, the average minimum cost of error C̅_min is also σ^2; however, the estimates or the eigenvalues of M̂_min are g_0, 
as can be seen in eigenvM. This means that for long interaction times the inference of the parameter g from the measurement data only yields the mean g_0 of the probability distribution p(g). 
Since the average minimum cost of error attains its maximum at both t=0 and t→∞, we are going to neglect these situations and focus on intermediate times, when C̅_min decreases. 
The fact that average minimum cost of error reaches a minimum at a finite time implies the existence of a particular duration for the interaction that yields the greatest amount of information on g. 
For each set of parameters, we can determine the time t^∗ as the time when C̅_min reaches its minimum value. We can work backwards to obtain the specific 
M̂^∗_min=M̂_min(t=t^∗) to be measured, which is the measurement that best protects against information loss.



The value of α, the amplitude of the initial mechanical oscillator coherent state, has a strong influence on the value of t^∗. We show in costfig that the limit α→∞, 
with α∈ℝ, results in t^∗=0 and the lowest observed value for C̅_min≈0.636σ^2. However, the eigenvalues of M̂^∗_min are still zero and g_0 
at t=0, from which it follows that highly excited initial states of the mechanical oscillator result in an estimation scenario where measuring M^∗_min merely reinforces prior knowledge and 
yields only the mean g_0 of the prior probability distribution p(g). In the next step, we investigate the position of the minimum for α∈ℂ, to deduce its dependence on the phase of α. 
costimalpha shows a shift of t^∗ towards higher values and an increase of the minimum value of C̅_min as the imaginary part of α gets larger. We see that the case with very 
large |α| may lead to inconclusive measurement scenarios because C̅_min is only significantly smaller than its maximum for a short time period. This observation is of significance 
in the discussion of initial thermal states, since it implies that higher initial temperatures will degrade the quality of the estimation procedure.



Let us turn our attention to M̂^∗_min, which has already been defined as the optimal measurement, made at the time t^∗ that minimizes the average minimum cost of error. Every outcome of 
the measurement of M̂^∗_min is an estimate of g. The most important quantity for a possible experimental implementation is the average estimate at t=t^∗

    h(g)={M̂^∗_minρ̂_F(g)}.

Thus, measurement data determine the value of h(g). From this, one may deduce the value of g. In avestimator, we show the curves of h(g) for different values of the real parameter 
α. t^∗ is independently calculated for each specific initial state. When α=0, the average estimator is an even function of g. This is a direct consequence of our particular 
choice of the cost function (<ref>), which is also an even function.

Before turning our attention to initial thermal states, let us conclude this section by summarizing the measurement procedure. Given a specific initial state, the system is allowed to evolve for a time t^∗. 
At this point in time, one would conduct a measurement of M̂^∗_min. This process is repeated, obtaining an average measurement, h(g). Using calculations of the kind shown in avestimator 
allows one to work backward and obtain g.





 §.§ Thermal state


Similar to the case for an initial coherent state for the mechanical oscillator, an initial thermal state exhibits an average minimum cost of error C̅_min that is equal to σ^2 for t=0 
and t→∞. In these limits, the eigenvalues of M̂_min have the same values as for the coherent state, so our earlier observations hold for the present case as well. In thermal, we 
show the time dependence of C̅_min and the average estimator h(g) for different average phonon numbers n_th obtained from the density matrix (<ref>) with the help of the 
expressions in example2. One can observe that an increase in the value of n_th increases C̅_min for most times, while it does not induce any significant change in the average 
estimator
h(g). Furthermore, the oscillations seen in costfig for longer interaction times are damped by the increase of n_th. Thus, in the context of this optimal estimation scenario the lower the 
temperature T of the mechanical oscillator, the less sensitive is the average minimum cost of error. This provides additional impetus to one of the central pillars of optomechanical experiments, which is to 
cool down the mechanical oscillator to temperatures as low as possible <cit.>.





 §.§ Squeezed state


We make use of example3 to construct the density matrix in Aform. The properties of C̅_min and M̂_min for t=0 and t→∞ are essentially the same as in the 
two cases discussed above. Let us recall that, in our discussion above, we showed that for an initial coherent state |α⟩ of the mechanical oscillator large |α|^2 reduces the
average minimum cost of error C̅_min, but at the expense of pushing the minimum towards very short interaction times, which leads to inconclusive measurement scenarios. In our discussion above we 
also identified a preferable scenario where |α|^2 is large but with approximately equal real and imaginary parts. In the case of initial squeezed states, squeezed shows an interesting effect, 
namely the squeezing parameter ζ may also reduce the average minimum cost of error. The average estimator h(g) is again an even function, this time because we have chosen two squeezed states without 
displacement. The cost of error C̅_min attains a minimum when the squeezing angle lies between the position and momentum quadratures of the oscillator. This can be understood as an effective 
continuous sampling of the noise ellipse during the evolution of the system for the first fraction of a mechanical time-period. Squeezing along either position or momentum quadrature will result in a greater 
uncertainty, whereas squeezing at an angle half-way between these two quadratures allows the measurement to take place with the least possible uncertainty.





 §.§ Different initial photonic states


So far we have discussed in detail the estimation problem of the optomechanical coupling g for the simplest initial state of the single-mode field. In this section, we consider the situation where the optical 
field may have more than one photon, and where the mechanical oscillator is initially in the ground state. Since ρ̂_F(g) depends on g only in its off-diagonal elements, we therefore set the 
amplitude of all participating photon number states to be equal. This ensures the maximum allowed absolute value for the off-diagonal elements in the density matrix. Due to the added complexity of dealing with 
intx we restrict our comparison to the following family of initial states of the optical field, indexed by the parameter N=2,3,4:

    |ψ_N⟩_c=∑_n=0^N-1a_n |n⟩_c= 1/√(N)∑_n=0^N-1|n⟩_c.

Our earlier investigations consider exclusively the case N=2.



costdim shows that the average minimum cost of error is reduced as the number of the photons in the initial state increases. This can be understood by examining carefully the Hamiltonian in optoH,
which reveals that the interaction between the single mode field and the mechanical oscillator gets stronger with increased number of participating photons. Thus, we have a better chance to estimate the 
optomechanical coupling g. The time t^∗ when C̅_min attains its minimum does not change markedly with N. We have also calculated the average estimator h(g) for t^∗; 
estimatordim shows the three curves obtained. Since the mechanical oscillator is initially in the ground state in every case, all curves are even.





§ QUANTUM CRAMÉR–RAO-TYPE INEQUALITY


In the preceding sections we discussed the properties of the optimum Hermitian operator M̂_min which minimizes the average cost in averagecost, and the eigenvalues of which are the estimates 
of the unknown optomechanical coupling strength g. An important task is to find out the accuracy with which g can be estimated. We would like to employ here the quantum Cramér–Rao inequality, which is 
widely used in the case of unbiased estimators <cit.>. In the present case, however, we have a biased estimator

    {ρ̂_F(g) (M̂_min-g Î)}=f(g),

where f(g) is the bias of the estimation and is not necessarily equal to zero. To properly account for this situation, we have to review the derivation of the Cramér–Rao inequality.

Let us first, however, deal with an extra issue regarding the derivative of the density matrix ρ̂_F(g) with respect to the parameter g. For concreteness, let us recall the density matrix 
ρ̂_F(g) from densitym, together with example1, and observe that

    ρ̂_F(g)=∑^∞_n,m=0 a_n a^∗_m e^-a_1 (n-m)^2+a_2 (n^2-m^2)-a_3 (n-m)|n⟩⟨m|,

where

    a_1=   g^2/ω^2_m[1-cos(ω_m t) ],
    
    a_2=   i g^2/ω^2_m[ω_m t-sin(ω_m t) ], and
    
    a_3=   i ω_c t-g/ω_m[α^∗(1-e^iω_m t)-α(1-e^-iω_m t) ].

Therefore,

    ∂ρ̂_F(g)/∂ g   =-∂ a_1/∂ g[â^†â,[â^†â,ρ̂_F(g)]]
           +∂ a_2/∂ g[(â^†â)^2,ρ̂_F(g)]
           -∂ a_3/∂ g[â^†â,ρ̂_F(g)]=ℒ[ρ̂_F(g)],

which demonstrates that ℒ[ρ̂_F(g)] does not have the form of either a right logarithmic or a symmetrized logarithmic derivative of the density matrix 
ρ̂_F(g) by default (see the definitions in RSLD). Therefore, we need the spectral decomposition of ρ̂_F(g) to construct at least the symmetrized logarithmic 
derivative operator, which is very challenging due to the fact that we have to deal with states defined on an infinite dimensional Hilbert space. Although this problem can be easily circumvented in 
numerical simulations, here we are motivated to derive an analytically expressible lower bound. This situation will result in a departure from the standard analysis <cit.>. In the standard proof, 
a Cauchy–Schwartz–Bunyakovsky inequality is employed, which suggests that in our new situation we would have to introduce the operator ρ̂^-1/2_F(g). This operator does not exist when the 
spectrum of ρ̂_F(g) contains zero (e.g., a pure state). We avoid this situation by following a different path.

In order to derive a lower bound for the mean-squared error,

    MSE(M̂_min)={ρ̂_F(g) (M̂_min-g Î)^2},

we define

    x_1(g)={ρ̂^2_F(g) M̂_min},

and then we differentiate both sides with respect to the parameter g,

    {(∂ρ̂_F(g)/∂ gρ̂_F(g)+ρ̂_F(g)
     ∂ρ̂_F(g)/∂ g)
     M̂_min}=x'_1(g).

We also define x_2(g)={ρ̂^2_F(g)}, and find

    {(∂ρ̂_F(g)/∂ gρ̂_F(g)+ρ̂_F(g)
     ∂ρ̂_F(g)/∂ g) g Î}=g x'_2(g).

Subtracting three from two, we obtain

    {(∂ρ̂_F(g)/∂ gρ̂_F(g)+ρ̂_F(g)
     ∂ρ̂_F(g)/∂ g) (M̂_min-g Î)}
    
    =x'_1(g)-g x'_2(g)=x(g),

where the last equality defines the function x(g). We make use of Liou and write four as

    {(ℒ[ρ̂_F]ρ̂_F+ρ̂_Fℒ[ρ̂_F]) (M̂_min-g Î)}=x(g),

where for the sake of notational simplicity we have omitted the argument g of ρ̂_F(g).



Before continuing, we discuss an issue connected with the boundedness of ρ̂_F. The Banach space of the Hilbert–Schmidt operators is defined as

    ℬ_2(ℋ):={X̂∈ℬ(ℋ): {X̂^†X̂}<∞},

where ℬ(ℋ) is Banach space of all bounded operators defined on the Hilbert space ℋ. The space ℬ_2(ℋ) with the inner product

    ⟨ A, B ⟩={A^† B},

where A,B ∈ℬ_2(ℋ), is a Hilbert space <cit.>. The Cauchy–Schwartz–Bunyakovsky inequality reads

    |{A^† B}|^2⩽{A^† A}{B^† B}.

In our case the Hilbert space is the symmetric Fock space, i.e., ℋ=Γ_s(ℂ), and ℒ contains powers of â^†â, which is an unbounded operator. 
This clearly shows that our proof is limited to density matrices which fulfill the conditions ρ̂^1/2_F (â^†â)^2 ρ̂_F, 
ρ̂^1/2_Fâ^†âρ̂_Fâ^†â∈ℬ_2(Γ_s(ℂ)). These conditions, together with the cyclic property 
of the trace, imply that ρ̂^1/2_Fℒ(ρ̂_F) is a Hilbert–Schmidt operator. Similarly, the condition ρ̂^1/2_FM̂_min∈ℬ_2(Γ_s(ℂ)) may restrict further the set of the density matrices. In other words, there are restrictions on the choice of the a_n in the 
initial state initial. In the case of finite dimensional examples, i.e, if there exists an N>0 such that a_n=0 for n ⩾ N, these complications do not arise, because all matrices are 
Hilbert–Schmidt operators. This is the typical case encountered in numerical simulations.







Now, provided that ρ̂^1/2_Fℒ[ρ̂_F] and ρ̂^1/2_FM̂_min are Hilbert–Schmidt operators, fh implies

    |x(g)|=|{ℒ(ρ̂_F)ρ̂^1/2_Fρ̂^1/2_F(M̂_min-g Î)}
    
     +{ρ̂^1/2_Fℒ(ρ̂_F) (M̂_min-g Î)ρ̂^1/2_F}|.

Applying first the subadditivity of the absolute value and then the Cauchy–Schwartz–Bunyakovsky inequality (<ref>) twice, we find

    |x(g)|^2 ⩽ 2{ρ̂_F(ℒ[ρ̂_F])^2 }MSE(M̂_min),

where we have used the fact that (ℒ[ρ̂_F])^†=ℒ[ρ̂_F], as can be deduced from Liou.

Finally, we obtain a lower bound for the mean-squared error

    MSE(M̂_min) ⩾|x(g)|^2/2{ρ̂_F(ℒ[ρ̂_F])^2 }.

The quantity on the right of this inequality is very similar to the standard quantum Cramér–Rao bound. In this expression, the function of x(g) in the numerator represents the fact that the estimator is 
biased, and includes information about the purity of the density matrix ρ̂_F(g). The denominator has a similar but slightly more complicated structure than the quantum Fischer 
information <cit.>, due to our approach to finding the derivative of ρ̂_F(g) with respect to the optomechanical coupling g.

We shall now apply the technique we just described to study the same cases as we did above, allowing us to study how the the mean-squared error behaves in each case.



 §.§ Coherent state

We investigate numerically the lower bound of the mean-squared error. We consider g_0 and σ to be the same as in mst. In cr1, we recall the results of costfig
and <ref>, and show the behavior of the lower bound of the mean-squared error as a function of g/ω_m. The most interesting feature occurs when |α| grows, where the lower 
bound is the smallest. This may seem to suggest that measurement strategies perform better under these conditions. However, this in apparent contrast with our findings in OM. What we can 
deduce is that measurements made with large |α| may simply return g_0, i.e., our prior expectation, for the value of the coupling strength. In such circumstances, we gain no information 
about the system; these scenarios are therefore to be avoided.



 §.§ Thermal state

Let us consider again the parameters of thermal, where we have seen that the average estimator is insensitive to the change of the average phonon number n_th, i.e., the change in the 
temperature of the mechanical oscillator. Here, we observe something different for the lower bound of the mean-squared error, (see thq). These findings indicate that the accuracy of the 
measurements is slightly worsened with the increase of n_th. However, we have also demonstrated an increase in the average minimum cost of error as n_th is increasing. Therefore, 
in accordance with intuition, high temperatures once again lead to inconclusive estimation results.



 §.§ Squeezed state

As we have seen in squeezed for the average minimum cost of error, squeezing is beneficial in the sense of reducing the mean-squared error. In sqq we see the lower bound of the mean-squared error 
may also be reduced by squeezing, in a manner that depends highly on the squeezing angle as well as the magnitude of the squeezing. This is, once again, in accordance with our earlier arguments and 
with intuition.



 §.§ Different initial photon states

Finally, we compare the lower bound of the mean-squared error for the three different initial single-mode field states given in initialstate, with the mechanical oscillator again assumed to be in 
its ground state. The time t^∗, when the average minimum cost of error C̅_min attains its minimum is taken to be the same as in estimatordim. The lower bound of the mean-squared error, 
as shown in cr2, generally decreases with the photon number states in the initial state. This is in agreement with our findings in costdim, namely that the average minimum cost of 
error is reduced by the increase of the photon number states. This also suggests that the initial preparation of the optical field is crucial to the outcome of the estimation procedure, with an equally 
weighted superposition of many photon number states being preferable. For ranges of values of g, however, either this improvement with increasing N is not seen, or in some cases, the situation worsens 
as N increases.

We conclude from this qualitative assessment that preparation of initial states of both the optical field and mechanical oscillator is crucial to obtaining more precise measurement outcomes, and 
consequently better estimations of the optomechanical coupling strength.



§ CONCLUDING REMARKS


In this paper we have investigated the simplest optomechanical model, consisting of a single mode of an optical field interacting with a single vibrational mode of a mechanical oscillator, from the 
perspective of quantum estimation theory. The object of our analysis was to determine the optomechanical coupling strength optimally, based on measurements made on the optical field. We have discussed 
this problem by introducing a quantum estimation scenario, in which one seeks for the best estimator, which minimizes the mean-square error cost functional. This Bayesian-inference approach requires a 
prior probability density function of the coupling strength, which represents the limited prior information held about the system. In particular, we have considered a normal distribution, where the 
mean and the standard deviation have been set to values emerging from the derivation of the radiation pressure Hamiltonian <cit.>. This derivation motivates our analysis, which develops an 
estimation procedure that results in a updated posterior probability density function for the coupling strength.

We have concentrated on the average mean-square error estimator, where the measurements occur at those interaction times where the average minimum cost of error reaches a minimum. The estimates are the 
eigenvalues of this estimator, with the eigenvectors determining a projective POVM that implements the measurement strategy. Our analysis has shown that highly excited initial coherent states of the 
mechanical oscillator limit the efficiency of this estimation procedure, unless the imaginary and the real parts of the displacement amplitude are approximately equal. We have demonstrated that the most 
promising estimates involve measurements being made during the first time period of the mechanical oscillation. We have, moreover, explored the effect of increasing the photon number states involved in the 
state of the optical field, sticking to the case of an equally weighted superposition of photon number states, initialstate; we find that increasing photon numbers reduces the average information loss. 
Furthermore, we have investigated scenarios where the mechanical oscillator is initially in a thermal state or a squeezed state. In general, thermal states lead to inconclusive measurement outcomes, 
where the updated posterior probability density function is the same as the prior one. The situation with an initial squeezed state is different, because we find that for certain choices of squeezing angle, 
squeezing reduces the average minimum cost of error.

Third, we have investigated the accuracy of the mean-square error estimator by means of a lower bound for the mean-squared error. The quantum Cramér–Rao inequality, defining this lower 
bound, is derived with the help of a symmetrized logarithmic derivative operator. In our situation, this operator is demanding to construct due to the infinite dimensionality of the Hilbert space on which 
the states to be measured are defined. Therefore, we have derived a new lower bound, lbound, for the mean-squared error of our biased estimator. In fact, we have reproduced the derivation of the 
quantum Cramér–Rao inequality by applying its standard methods to our case. Our numerical investigations here largely corroborate our previous conclusions. However, the lowest bounds for the estimation 
accuracy have been found for those limiting cases when the eigenvalues of the estimator are either zero or the mean of the prior normal distribution, where measurement yields no further information about the 
system. In particular, we have found that the initial state of the mechanical oscillator has to be carefully prepared, otherwise the outcome of the measurement process will be to simply reinforce prior 
expectations of the optomechanical coupling strength.

Finally let us make some comments on our approach. The analysis clearly indicates a characteristic set of parameters when the estimation of the optomechanical coupling can be done with minimal loss of 
information. Despite the fact that our results pinpoint some important results for a scenario of much experimental relevance, the question of how to implement the optimal detection strategy or to 
compare with less optimal but implementable measurement setups (see Ref. <cit.>) has not been answered, and is the subject of ongoing investigations. Another critical point is the preparation 
of the initial state of the optical field; in this paper we have considered this state to be an equally-weighted superposition of photon number states, but we have not tackled the question of whether this 
family of states is optimal. These questions define the direction of our future investigations. As a final word we think that the present paper may offer interesting perspective and viewpoint, which 
provides a different way of thinking about optomechanical systems.



§ ACKNOWLEDGEMENT

J.Z.B. is grateful to Matteo G. A. Paris for stimulating discussions. This paper is supported by the European Union's Horizon 2020 research and innovation program under Grant Agreement 
No. 732894 (FET Proactive HOT––Hybrid Optomechanical Technologies).




§ MECHANICAL OSCILLATOR IN AN INITIAL SQUEEZED STATE


In this appendix, we present the full expressions of I^(0), I^(1)_n,m, and I^(2)_n,m, which appear in example3. First, we introduce the following notation:

    χ^(1)_n,m   =n(e^-iω_m t-1)-m(e^iω_m t-1),
    χ^(2)_n,m   =n(e^-iω_m t-1)+m(e^iω_m t-1),
    χ^(3)_n,m   =n(1-e^-iω_m t)e^i(θ-ω_m t)-m(1-e^iω_m t)e^-i(θ-ω_m t), and
    χ^(4)_n,m   =n(1-e^-iω_m t)e^i(θ-ω_m t)+m(1-e^iω_m t)e^-i(θ-ω_m t),

as well as

    ξ^(0)   =1/4[1-tanh(|ζ|)cos(z)],
    ξ^(1)   =1-tanh(|ζ|)cos(z)/4[1-tanh^2(|ζ|)], and
    ξ^(2)   =itanh(|ζ|)sin(z)/1-tanh(|ζ|)cos(z),

with z=θ-2ω_m t. Finally, we can write

    I^(2)_n,m   =ξ^(0)[χ^(1)_n,m+tanh^2(|ζ|)χ^(3)_n,m]^2-ξ^(1){χ^(2)_n,m+tanh^2(|ζ|)χ^(4)_n,m+ξ^(2)[χ^(1)_n,m+
    tanh^2(|ζ|)χ^(3)_n,m]}^2,
    
    I^(1)_n,m   =4|α|(ξ^(0)[χ^(1)_n,m+tanh^2(|ζ|)χ^(3)_n,m]i[sin(z_1)+tanh(|ζ|)sin(z_2)]+ξ^(1){χ^(2)_n,m+
    tanh^2(|ζ|)χ^(4)_n,m+ξ^(2)[χ^(1)_n,m
           +tanh^2(|ζ|)χ^(3)_n,m]}{cos(z_1)+tanh(|ζ|)cos(z_2)-iξ^(2)[sin(z_1)+tanh(|ζ|)sin(z_2)]}), and
    
    I^(0)   =4|α|^2(ξ^(0)[sin(z_1)+tanh(|ζ|)sin(z_2)]^2+ξ^(1){[cos(z_1)+tanh(|ζ|)cos(z_2)]
           -iξ^(2)[sin(z_1)+tanh(|ζ|)sin(z_2)]}^2),

where z_1=ω_m t -ϕ and z_2=ω_m t -ϕ-θ.




§ THE SYMMETRIZED AND THE RIGHT LOGARITHMIC DERIVATIVE OPERATORS


In this appendix, we present some well-known material in order to support the arguments of this paper. In a single parameter estimation scenario, the symmetrized logarithmic derivative L̂ of the
density matrix ρ̂(x) is defined by

    ∂ρ̂/∂ x=12(L̂ρ̂+ρ̂L̂).

The operator L̂ is Hermitian <cit.>. If we consider the spectral decomposition 

    ρ̂=∑_i p_i |i⟩⟨i|,

then

    L̂=∑_i,j2⟨i|∂ρ̂/∂ x |j⟩/p_i+p_j|i⟩⟨j|

satisfies the above definition. However, in order to construct L̂ one must know the exact eigenvalues and the eigenvectors of ρ̂.

Another way of defining the derivative of ρ̂(x) with respect to x involves the non-Hermitian operator 𝕃 that is the solution of the equation

    ∂ρ̂/∂ x= ρ̂𝕃=𝕃^†ρ̂,

this being called the right logarithmic derivative operator. This operator may not exist for many density matrices, and in particular for those representing pure states <cit.>.

99

Helstrom1969 C. W. Helstrom, J. Stat. Phys. 1, 231 (1969).

Holevo1973 A. S. Holevo, J. Multivar. Anal. 3, 337 (1973). 

Holevo2011 A. S. Holevo, Probabilistic and Statistical  Aspects of Quantum  Theory (Edizioni della Normale, Pisa, 2011).

Paris2004 Quantum State Estimation, edited by M. G. A. Paris and J. Řeháček (Springer-Verlag, Berlin, 2004).

Aspelmeyer2014 M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Rev. Mod. Phys. 86, 1391 (2014).

Abbott2016 B. P. Abbott et al., Phys. Rev. Lett. 116, 061102 (2016).

Braginski1975 V. B. Braginsky and Yu. I. Vorontsov, Sov. Phys. Usp. 17, 644 (1975).

Caves1981 C. M. Caves, Phys. Rev. D 23, 1693 (1981).

Giovannetti2011 V. Giovannetti, S. Lloyd, and L. Maccone, Nature Photon. 5, 222 (2011).

Berni2015 A. A. Berni, T. Gehring, B. M. Nielsen, V. Händchen, M. G. A. Paris, and U. L. Andersen, Nat. Photon. 9, 577 (2015).

Hayashi2006  M. Hayashi, Quantum Information (Springer-Verlag, Berlin, 2006).

Seveso2017 L. Seveso, M. A. C. Rossi, and M. G. A. Paris, Phys. Rev. A 95, 012111 (2017).

Law1995 C. K. Law, Phys. Rev. A 51, 2537 (1995).

Bernad2013 J. Z. Bernád, H. Frydrych, and G. Alber, J. Phys. B 46, 235501 (2013); Appendix A.

Bose1997 S. Bose, K. Jacobs, and P. L. Knight, Phys. Rev. A 56, 4175 (1997).

Helstrom1974 C. W. Helstrom, Int. J. Theor. Phys. 8, 361 (1974).

Law1994 C. K. Law, Phys. Rev. A 49, 433 (1994).

Perelomov1986 A. Perelomov, Generalized Coherent States and Their Applications (Springer-Verlag, Berlin, 1986).

Glauber1963 R. J. Glauber, Phys. Rev. Lett. 10, 84 (1963).

Sudarshan1963 E. C. G. Sudarshan, Phys. Rev. Lett. 10, 277 (1963).

Schleich2001 W. P. Schleich, Quantum Optics in Phase Space (Wiley-VCH, Weinheim, 2001).

Helstrom1976 C. W. Helstrom, Quantum Detection and Estimation
Theory (Academic Press, New York, 1976).

Personick1971 S. D. Personick, IEEE Trans. Inf. Theory 17, 240 (1971).

Jaynes2003 E. T. Jaynes, Probability Theory: The Logic of Science (Cambridge University Press, Cambridge, 2003).

Helstrom1968 C. W. Helstrom, IEEE Trans. Inf. Theory 14, 234 (1968).

Reed1972 M. Reed and B. Simon, Methods of Modern Mathematical Physics I: Functional Analysis (Academic Press, San Diego, 1980).

Braunstein1994 S. L. Braunstein and C. M. Caves, Phys. Rev. Lett. 72, 3439 (1994).



