structure: this namespace comprises classes wherefrom all elements must be derived to present the necessary interfaces to trajectory drivers (cf. the quantumtrajectory namespace)
quantumoperator: here, such operator classes are found that represent special operator structures ubiquitous in atomic, molecular, and optical physics.

This set of tools constitutes the Level-2 interface of C++QED, below Level-1, which is the interface for writing scripts from already available elements (cf. userguide).

Note: The files referred to in this Section are found in directory examples in the distribution. The code examples are expected to compile, cf. examples/Jamfile.

Implementation of a class representing a harmonic-oscillator mode

We demonstrate how to implement an element representing a pumped lossy mode in a truncated Fock space. In the frame rotating with the pump frequency, it is described by the Hamiltonian:

$H=-\delta a^\dagger a+\lp\eta a^\dagger+\hermConj\rp,$

where $\delta$ is the detuning between the oscillator and the pump frequencies. The Liouvillean:

$\Liou\rho=\kappa\lp(n_\text{Th}+1)\lp2a\rho a^\dagger-\comm{a^\dagger a}{\rho}_+\rp+n_\text{Th}\lp2a^\dagger\rho a-\comm{a\,a^\dagger}{\rho}_+\rp\rp$

The frequency-like parameters are $\delta$ , $\kappa$ , and $\eta$ , representing the mode frequency, loss rate, and pumping Rabi frequency, respectively. A dimensionless parameter is $n$ (the average number of thermal photons in the heat bath) and the cutoff of the Fock space.

Using the notation of Sec. Description of the MCWF method

$J_0=\sqrt{2\kappa(n_\text{Th}+1)}a,\quad J_1=\sqrt{2\kappa n_\text{Th}}a^\dagger.$

The non-Hermitian Hamiltonian for the Monte Carlo wave-function method reads:

$\HnH=\lp-\delta-i\kappa(2n_\text{Th}+1)\rp a^\dagger a+\lp\eta a^\dagger+\hermConj\rp\equiv-iz\,a^\dagger a+\lp\eta a^\dagger+\hermConj\rp,$

where we have introduced the complex frequency

$z\equiv\kappa(2n_\text{Th}+1)-i\delta.$

The element has to be represented by a class which inherits publicly from the necessary classes in the structure namespace. In this simple case, it is basically two helper functions returning quantumoperator::Tridiagonal instances, a constructor, and two virtual functions inherited from structure::ElementAveraged that have to be written.

See also: The description of quantumoperator::Tridiagonal

Consider the file ExampleMode.h:

#include "Free.h"
#include "ElementLiouvillean.h"
#include "ElementAveraged.h"
#include "TridiagonalHamiltonian.h"
using namespace structure;
using quantumoperator::TridiagonalHamiltonian;
using namespace freesystem;
const Tridiagonal aop(size_t dim); // ladder operator of the given dimension
const Tridiagonal nop(size_t dim); // number operator "
namespace basic {
class PumpedLossyMode 
  : public Free, public TridiagonalHamiltonian<1,false>, public ElementLiouvilleanStrategies<1,2>, public ElementAveraged<1>
{
public:
  typedef ElementAveraged<1>::LazyDensityOperator LazyDensityOperator;
  PumpedLossyMode(double delta, double kappa, dcomp eta, double nTh, size_t cutoff);
private:
  const Averages average_v(NoTime, const LazyDensityOperator&) const;
};
} // basic
inline const Tridiagonal aop(const basic::PumpedLossyMode& mode) {return aop(mode.getDimension());} // just for convenience

Though the structure::NoTime tagging class as a function argument creates some redundancy, it is necessary because of the design of the structure-bundle. Note that structure::ElementLiouvilleanStrategies (like structure::ElementLiouvillean) assumes that the number of Lindblads is known @ compile time, as is the case here. If this is not the case, structure::Liouvillean has to be used instead.

This will suffice here. Let us look at the implementations in ExampleMode.cc:

 // Copyright András Vukics 2006–2014. Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE.txt)
 #include "ExampleMode.h"
 
 #include "ElementLiouvillean.tcc"
 #include "Tridiagonal.tcc"
 #include "TridiagonalHamiltonian.tcc"
 
 #include <boost/bind.hpp>
 
 using boost::bind;
 
 
 void aJump   (StateVectorLow&, double kappa_nPlus1);
 void aDagJump(StateVectorLow&, double kappa_n     );
 
 double aJumpRate   (const LazyDensityOperator&, double kappa_nPlus1);
 double aDagJumpRate(const LazyDensityOperator&, double kappa_n     );
 
 
 basic::PumpedLossyMode::PumpedLossyMode(double delta, double kappa, dcomp eta, double nTh, size_t cutoff)
   : Free(cutoff,
          {
            CF{"(kappa*(2*nTh+1),delta)",dcomp(kappa*(2*nTh+1),delta),     cutoff },
            CF{"eta"                    ,eta                         ,sqrt(cutoff)}
          }),

We construct the structure::Free base with the dimension of the system and the name-value-multiplier tuples for the frequency-like parameters of the system, which in this case are all complex (cf. the explanation @ structure::DynamicsBase). We use C++11 initializer lists and the structure::DynamicsBase::CF typedef.

     TridiagonalHamiltonian<1,false>(dcomp(-kappa*(2*nTh+1),delta)*nop(cutoff)
                                     +
                                     tridiagPlusHC_overI(conj(eta)*aop(cutoff))),

We construct the time-independent quantumoperator::TridiagonalHamiltonian base. This is greatly facilitated by the algebra and helpers of the quantumoperator::Tridiagonal class.

Warning: When implementing the Hamiltonian, not itself but $\frac Hi$ has to supplied!

     ElementLiouvilleanStrategies<1,2>(JumpStrategies(bind(aJump   ,_1,kappa*(nTh+1)),
                                                      bind(aDagJump,_1,kappa* nTh   )),
                                       JumpRateStrategies(bind(aJumpRate   ,_1,kappa*(nTh+1)),
                                                          bind(aDagJumpRate,_1,kappa* nTh   )),
                                       "Mode",{"photon loss","photon absorption"}),

We construct the structure::ElementLiouvilleanStrategies base, whose second template argument denotes the number of different quantum jumps, which is 2 in this case. The constructor takes the strategies for calculating the impact of a jump on a structure::freesystem::StateVectorLow, and for calculating the rate from a structure::freesystem::LazyDensityOperator. These strategy functions are produced from the free-standing helpers in Lines 10-14 above through argument binding. The strategies are followed by a description of the lossy element and the decay channels. The number of descriptive labels must agree with the number of strategies.

1 ElementAveraged<1>("PumpedLossyMode",{"<number operator>","real(<ladder operator>)","imag(\")"})

We construct the structure::ElementAveraged base, with parameters necessary to produce a simple key of the quantum averages that are communicated towards the user. Here we calculate only three such averages, the expectation value of the number operator, and the real and imaginary parts of that of the ladder operator.

 {
   getParsStream()<<"# Pumped lossy mode";
 }

With the structure::DynamicsBase::getParsStream function we obtain a stream whereon we can write more information about the object that gets communicated towards the user in that part of the output which summarizes the parameters of the actual run.

Next, the inherited function structure::Averaged::average is implemented (according to the non-virtual interface idiom, its virtual counterpart is called average_v):

 auto basic::PumpedLossyMode::average_v(NoTime, const LazyDensityOperator& matrix) const -> const Averages
 {
   auto averages(initializedAverages());
 
   averages(0)=aJumpRate(matrix,1);

the expectation value of the photon number is calculated (where we can reuse our function aJumpRate, with unit loss rate).

 
   for (int n=1; n<int(matrix.getDimension()); n++) {
     double sqrtn=sqrt(n);
     dcomp offdiag(matrix(n)(n-1));
     averages(1)+=sqrtn*real(offdiag);
     averages(2)+=sqrtn*imag(offdiag);

the expectation value of the ladder operator is calculated (real & imaginary parts)

   }
 
   return averages;
 
 }

The implementation of the helpers is also quite straightforward. It may come to a separate file ExampleModeImpl.cc:

 // Copyright András Vukics 2006–2014. Distributed under the Boost Software License, Version 1.0. (See accompanying file LICENSE.txt)
 #include "ExampleMode.h"
 #include "Tridiagonal.tcc"
 
 void aJump   (StateVectorLow& psi, double kappa_nPlus1)
 {
   double fact=sqrt(2.*kappa_nPlus1);
   int ubound=psi.ubound(0);
   for (int n=0; n<ubound; ++n)
     psi(n)=fact*sqrt(n+1)*psi(n+1);
   psi(ubound)=0;
 }
 
 
 void aDagJump(StateVectorLow& psi, double kappa_n     )
 {
   double fact=sqrt(2.*kappa_n);
   for (int n=psi.ubound(0); n>0; --n)
     psi(n)=fact*sqrt(n)*psi(n-1);
   psi(0)=0;
 }
 
 
 namespace {
 
 double photonNumber(const LazyDensityOperator& matrix)
 {
   double res=0;
   for (size_t n=1; n<matrix.getDimension(); ++n)
     res+=n*matrix(n);
   return res;
 }
 
 }
 
 
 double aJumpRate   (const LazyDensityOperator& matrix, double kappa_nPlus1)
 {
   return 2.*kappa_nPlus1*photonNumber(matrix);
 }
 
 
 double aDagJumpRate(const LazyDensityOperator& matrix, double kappa_n     )
 {
   return 2.*kappa_n*(photonNumber(matrix)+matrix.trace());
 }
 
 
 const Tridiagonal aop(size_t dim)
 {
   typedef Tridiagonal::Diagonal Diagonal;
   Diagonal diagonal(dim-1);
   return Tridiagonal(Diagonal(),1,Diagonal(),diagonal=sqrt(blitz::tensor::i+1.));

Warning: In situations like this, it is extremely important to write blitz::tensor::i+1., otherwise the blitz::sqrt function will operate within the integers, which is perhaps a shortcoming of the Blitz++ library in this respect.

 }
 
 
 const Tridiagonal::Diagonal mainDiagonal(const dcomp& z, size_t dim)
 {
   Tridiagonal::Diagonal res(dim);
   res=blitz::tensor::i;
   return res*=z;
 }
 
 
 const Tridiagonal nop(size_t dim)
 {
   return Tridiagonal(mainDiagonal(1.,dim));
 }

Exploiting interaction picture

In many situations, it pays to transfer to interaction picture defined by the first term of the Hamiltonian $\HnH$ above. The ladder operator in interaction picture is

$a\Int(t)=a e^{-zt},$

so that the Hamiltonian reads

$H\Int(t)=\lp\eta a^\dagger e^{zt}+\eta^*ae^{-zt}\rp.$

See also: These notes on how to treat interaction pictures defined by non-unitary transition operators in a consistent way.

In this case, the class representing the element has to be derived from structure::Exact as well, which provides an interface allowing for the transformation between the two pictures. In addition, instead of quantumoperator::TridiagonalHamiltonian <1,false>, we need to derive from quantumoperator::TridiagonalHamiltonian <1,true> because the Hamiltonian is now time-dependent.

Note: In general usage, the jump and the averages are calculated in the normal picture also in this case (cf. explanation of classes structure::Hamiltonian, structure::Exact, structure::Liouvillean, and structure::Averaged, and furthermore quantumtrajectory::MCWF_Trajectory). This allows for reusing the same code in both pictures. (Incidentally, here the Liouvillean remains unchanged anyway.)

 #include "FreeExact.h" // Further inculdes and using directives the same as above
 
 namespace basic {
 
 class PumpedLossyModeIP
   : public Free, public FreeExact<false>, public TridiagonalHamiltonian<1,true>, public ElementLiouvilleanStrategies<1,2>, public ElementAveraged<1>
 {
 public:
   typedef ElementAveraged<1>::LazyDensityOperator LazyDensityOperator;
 
   PumpedLossyModeIP(double delta, double kappa, dcomp eta, double nTh, size_t cutoff);
 
   const dcomp get_z() const {return z_;}
 
 private:
   void updateU(OneTime) const;
 
   bool applicableInMaster_v() const {return true;}
 
   const Averages average_v(NoTime, const LazyDensityOperator&) const;
 
   const dcomp z_; // Needed for updateU
 
 };
 
 } // basic

The structure::OneTime tagging class at the same time carries the information about the time instant to which the (diagonal) transformation operator has to be updated. structure::OneTime can be implicitly converted into a double.

In the implementation, the only difference from the previous case will be the constructor, because the Hamiltonian now also requires furnishing with frequencies (cf. quantumoperator::furnishWithFreqs), and the implementation of the virtual function structure::FreeExact::updateU.

Note: In connection with the structure::Exact::applicableInMaster virtual function implemented here, cf. also this note.

When “furnished with frequencies”, a quantumoperator::Tridiagonal object will internally take care about the time-dependent phases appearing in $a\Int(t)$ and $H\Int(t)$ :

 const Tridiagonal::Diagonal mainDiagonal(const dcomp& z, size_t cutoff);
 
 
 basic::PumpedLossyModeIP::PumpedLossyModeIP(double delta, double kappa, dcomp eta, double nTh, size_t cutoff)
   : Free(cutoff,
          {
            RF{"kappa*(2*nTh+1)",kappa*(2*nTh+1),cutoff},
            RF{"delta",delta,1}
          },
          CF{"eta",eta,sqrt(cutoff)}),
     FreeExact<false>(cutoff),
     TridiagonalHamiltonian<1,true>(furnishWithFreqs(tridiagPlusHC_overI(conj(eta)*aop(cutoff)),
                                                     mainDiagonal(dcomp(kappa*(2*nTh+1),-delta),cutoff))),
     ElementLiouvilleanStrategies<1,2>(JumpStrategies(bind(aJump   ,_1,kappa*(nTh+1)),
                                                      bind(aDagJump,_1,kappa* nTh   )),
                                       JumpRateStrategies(bind(aJumpRate   ,_1,kappa*(nTh+1)),
                                                          bind(aDagJumpRate,_1,kappa* nTh   )),
                                       "Mode",{"photon loss","photon absorption"}),
     ElementAveraged<1>("PumpedLossyMode",{"<number operator>","real(<ladder operator>)","imag(\")"}),
     z_(kappa*(2*nTh+1),-delta)
 {}
 
 
 void basic::PumpedLossyModeIP::updateU(OneTime dtDid) const
 {
   getDiagonal()=exp(-z_*(dtDid*blitz::tensor::i));
 }
 
 
 // PumpedLossyModeIP::average_v exactly the same as PumpedLossyMode::average_v above

structure::FreeExact assumes that the operator transforming between the two pictures is diagonal, and the factors to update are simply its diagonal elements. If this is not the case, Exact has to be used instead. Here, since there are also real frequency-like parameters, we have to use structure::DynamicsBase::RF as well.

Note: Since a lot of the code from the previous case can be reused here, one will usually adopt an inheritence- or class-composition-based solution to implement classes like PumpedLossyMode and PumpedLossyModeIP (for an inheritance-based solution, cf. below; for one based on class-composition, cf. the actual implementation of a harmonic-oscillator mode in the framework in elements/frees/Mode_.h).

Implementing an X-X interaction

Let us consider the interaction described by the Hamiltonian

$H_\text{X-X}=g(a+a^\dagger)(b+b^\dagger).$

The class implementing this interaction has to be derived from structure::Interaction <2> because it is a binary interaction, and quantumoperator::TridiagonalHamiltonian <2,...> (note that quantumoperator::Tridiagonal is capable to represent direct products of tridiagonal matrices in the case of RANK>1).

The only thing requiring some care is that once we transform some elements into interaction picture, the whole Hamiltonian is transformed, that is, $a$ or $b$ or both may be in interaction picture. Here, for the sake of simplicity, we assume that both constituents are of the type PumpedLossyMode. Hence, the Hamiltonian is in fact

$H_{\text{X-X;I}}(t)=g(ae^{-z_at}+a^\dagger e^{z_at})(be^{-z_bt}+b^\dagger e^{z_bt}).$

Consider ExampleInteraction.h:

 #include "ExampleMode.h"
 
 #include "Interaction.h"
 #include "TridiagonalHamiltonian.tcc"
 #include "StateVector.tcc"
 
 namespace basic {
 
 class InteractionX_X
   : public Interaction<2>, public TridiagonalHamiltonian<2,false>
 {
 public:
   InteractionX_X(const PumpedLossyMode&, const PumpedLossyMode&, double g);
 
 };
 
 } // basic

ExampleInteraction.cc then reads

 #include "ExampleInteraction.h"
 
 basic::InteractionX_X::InteractionX_X(const PumpedLossyMode& m0, const PumpedLossyMode& m1, double g)
   : Interaction<2>(m0,m1,RF{"g",g,sqrt(m0.getDimension()*m1.getDimension())}),
     TridiagonalHamiltonian<2,false>(g*
                                     (aop(m0)+aop(m0).dagger())*
                                     (aop(m1)+aop(m1).dagger())
                                     )
 {}

As we see, the Hamiltonian can be written in a rather straightforward way, and it internally takes care about the time-dependent phases appearing in $H_{\text{X-X;I}}(t)$ , which result from the use of interaction picture.

The free elements are stored as shared pointers in the interaction element, and the constant references supplied to the constructor are turned into (non-owning) shared pointers (cf. cpputils::sharedPointerize). Of course, the free elements have to live in a larger scope than the interaction, otherwise we may run into trouble with dangling pointers.

Using class inheritance

For an inheritance-based solution, it pays to define a base class collecting all the common services. Consider the following snippet from ExampleMode.h:

 namespace hierarchical {
 
 // All inculdes and using directives the same as above
 
 class ModeBase : public Free, public ElementLiouvillean<1,2>, public ElementAveraged<1>
 {
 public:
   typedef ElementLiouvillean<1,2>::StateVectorLow StateVectorLow     ;
   typedef ElementAveraged<1>::LazyDensityOperator LazyDensityOperator;
 
 protected:
   ModeBase(double kappa, double nTh, size_t cutoff);
   
 private:
   void doActWithJ(NoTime, StateVectorLow&, LindbladNo<0>) const;
   void doActWithJ(NoTime, StateVectorLow&, LindbladNo<1>) const;
   
   double rate(NoTime, const LazyDensityOperator&, LindbladNo<0>) const;
   double rate(NoTime, const LazyDensityOperator&, LindbladNo<1>) const;
   
   const Averages average_v(NoTime, const LazyDensityOperator&) const;
 
   const double kappa_, nTh_; // needed for calculating jumps & rates 
   
 };
 
 } // hierarchical

Here, instead of structure::ElementLiouvilleanStrategies, we can rather use structure::ElementLiouvillean, which has as many virtual functions doActWithJ and rate as there are jumps (indicated by the second template argument), distinguished by the tagging classes structure::lindblad::Base::LindbladNo. It results in a compile-time error to instantiate such a class with an argument not smaller than the number of Lindblads (since the numbering of jumps begins with 0). Via this solution we can get around the awkwardness of specifying the jump and rate strategies for structure::ElementLiouvilleanStrategies, while retaining a way of controlling the number of Lindblads @ compile time.

Deriving from ModeBase, the definition of PumpedLossyMode is trivial, while for PumpedLossyModeIP, we have to define the virtual functions inherited from structure::FreeExact:

 namespace hierarchical {
 
 
 class PumpedLossyMode
   : public ModeBase, public TridiagonalHamiltonian<1,false>
 {
 public:
   PumpedLossyMode(double delta, double kappa, dcomp eta, double nTh, size_t cutoff);
 
 };
 
 
 class PumpedLossyModeIP
   : public ModeBase, public FreeExact<false>, public TridiagonalHamiltonian<1,true >
 {
 public:
   PumpedLossyModeIP(double delta, double kappa, dcomp eta, double nTh, size_t cutoff);
 
   const dcomp get_z() const {return z_;}
 
 private:
   void updateU(OneTime) const;
 
   bool applicableInMaster_v() const {return true;}
 
   const dcomp z_; // Needed for updateU
 
 };
 
 
 } // hierarchical

Note: The correct treatment of frequency-like parameters would require more care in this case, since ModeBase does not know about delta & eta

The implementations come in ExampleMode.cc:

 hierarchical::ModeBase::ModeBase(double kappa, double nTh, size_t cutoff)
   : Free(cutoff,RF{"kappa*(2*nTh+1)",kappa*(2*nTh+1),cutoff}),
     ElementLiouvillean<1,2>("Mode",{"photon loss","photon absorption"}),
     ElementAveraged<1>("PumpedLossyMode",{"<number operator>","real(<ladder operator>)","imag(\")"}),
     kappa_(kappa), nTh_(nTh)
 {}
 
 
 void hierarchical::ModeBase::doActWithJ(NoTime, StateVectorLow& psi, LindbladNo<0>) const
 {
   aJump(psi,kappa_*(nTh_+1));
 }
 
 void hierarchical::ModeBase::doActWithJ(NoTime, StateVectorLow& psi, LindbladNo<1>) const
 {
   aDagJump(psi,kappa_*nTh_);
 }
 
 double hierarchical::ModeBase::rate(NoTime, const LazyDensityOperator& matrix, LindbladNo<0>) const
 {
   return aJumpRate(matrix,kappa_*(nTh_+1));
 }
 
 double hierarchical::ModeBase::rate(NoTime, const LazyDensityOperator& matrix, LindbladNo<1>) const
 {
   return aDagJumpRate(matrix,kappa_*nTh_);
 }
 
 // ModeBase::average_v exactly the same as PumpedLossyMode::average_v above

and the derived classes:

 hierarchical::PumpedLossyMode::PumpedLossyMode(double delta, double kappa, dcomp eta, double nTh, size_t cutoff)
   : ModeBase(kappa,nTh,cutoff),
     TridiagonalHamiltonian<1,false>(dcomp(-kappa*(2*nTh+1),delta)*nop(cutoff)
                                     +
                                     tridiagPlusHC_overI(conj(eta)*aop(cutoff)))
 {}
 
 
 hierarchical::PumpedLossyModeIP::PumpedLossyModeIP(double delta, double kappa, dcomp eta, double nTh, size_t cutoff)
   : ModeBase(kappa,nTh,cutoff),
     FreeExact<false>(cutoff),
     TridiagonalHamiltonian<1,true>(furnishWithFreqs(tridiagPlusHC_overI(conj(eta)*aop(cutoff)),
                                                     mainDiagonal(dcomp(kappa*(2*nTh+1),-delta),cutoff))),
     z_(kappa*(2*nTh+1),-delta)
 {}
 
 
 // PumpedLossyModeIP::updateU exactly the same as above

We may define a function using runtime-dispatching for a ladder operator either furnished or not furnished with frequencies, depending on the actual type. It should be implemented via dynamic casting:

 const Tridiagonal aop(const hierarchical::ModeBase& mode)
 {
   using namespace hierarchical;
   size_t dim=mode.getDimension();
   if (const auto modeIP=dynamic_cast<const PumpedLossyModeIP*>(&mode))
     return furnishWithFreqs(aop(dim),mainDiagonal(modeIP->get_z(),dim));
   else
     return aop(dim);
 }

Interaction element in the inheritance-based design

Based on the above design of the mode-element and the dispatching ladder-operator function, we may define an interaction element that works correctly with either free element, if the constructor expects arguments of type ModeBase:

 namespace hierarchical {
 
 class InteractionX_X
   : public Interaction<2>, public TridiagonalHamiltonian<2,true>
 {
 public:
   InteractionX_X(const ModeBase&, const ModeBase&, double g);
 
 };
 
 } // hierarchical

The implementation is formally equivalent to the previous:

 hierarchical::InteractionX_X::InteractionX_X(const ModeBase& m0, const ModeBase& m1, double g)
   : Interaction<2>(m0,m1,RF{"g",g,sqrt(m0.getDimension()*m1.getDimension())}),
     TridiagonalHamiltonian<2,true>(g*
                                    (aop(m0)+aop(m0).dagger())*
                                    (aop(m1)+aop(m1).dagger())
                                    )
 {}

Warning: Here, it would cause a hard-to-detect physical error to use quantumoperator::TridiagonalHamiltonian <2,false> instead of quantumoperator::TridiagonalHamiltonian <2,true>, because in the former case, the time update of the binary tridiagonal would not occur even with PumpedLossyModeIP.

Other uses of interaction elements

In the language of the framework, every element that operates on more than one quantum numbers is an interaction.

Hence, if we need for instance an element calculating correlations between two free subsystems (or, two quantum numbers in general), it has to be derived from Interaction because only an interaction element has a chance to access more than quantum numbers.

Assume we need an element calculating the dynamics

Since now we need to operate

at this point, the averages of an X-X interaction between two modes as above, but it also lets the user monitor the correlations $\avr{XQ},\;\avr{XP}\;\avr{YQ},$ and $\avr{YP}$ between the modes, where $X,\;Y$ and $Q,\;P$ are the quadratures of the two modes, respectively. The element can be derived from the former interaction element in the inheritance-based design:

1 namespace hierarchical {

2

3 class InteractionX_X_Correlations : public InteractionX_X, public ElementAveraged<2>

4 {

5 public:

6 InteractionX_X_Correlations(const ModeBase& m0, const ModeBase& m1, double g)

7 : InteractionX_X(m0,m1,g),

8 ElementAveraged<2>("ModeCorrelations",{"<XQ>","<XP>","<YQ>","<YP>"})

9 {}

10

11 private:

12 const Averages average_v(NoTime, const LazyDensityOperator&) const;

13

14 };

15

16 } // hierarchical

on two quantum numbers to calculate the quantum averages, we derived from structure::ElementAveraged <2>, which operates on a binary quantumdata::LazyDensityOperator. The implementation of the averaging function may read

1 auto hierarchical::InteractionX_X_Correlations::average_v(NoTime, const LazyDensityOperator& matrix) const -> const Averages

2 {

3 auto averages(initializedAverages());

4

5 for (int n=0; n<int(matrix.getDimension(0)); n++) for (int m=1; m<int(matrix.getDimension(1)); m++) {

6 if(n<int(matrix.getDimension(0))-1) {

7 dcomp temp=sqrt(m*(n+1))*matrix(n,m)(n+1,m-1);

8 averages(0)+=real(temp);

9 averages(1)+=imag(temp);

10 }

11 if(n>0) {

12 dcomp temp=sqrt(m*n)*matrix(n,m)(n-1,m-1);

13 averages(2)+=real(temp);

14 averages(3)+=imag(temp);

1 }

array contains the real and imaginary parts of $\avr{a^\dagger b}$ and $\avr{a b},$ respectively.

Now the desired set of quantum averages can be obtained via linear operations:

 
   double
     xq= averages(0)+averages(2),
     xp= averages(1)+averages(3),
     yq=-averages(1)+averages(3),
     yp= averages(0)-averages(2);
 
   averages=xq,xp,yq,yp;
 
   return averages;
 
 }

Table of Contents