Initial conditions for calculation

[aleonau]

Hi. I would like to consider the evolution of the Tavis-Cummings (or Jaynes-Cummings for N = 1) model, which starts from the ground state of the two-level atoms and some COHERENT state of the driving field.

How one should introduce this initial COHERENT state in the code, which is not the vacuum state? For instance, I would like to have n_average = 25 photons.

Thank you

[modmido]

Hello, there is no ready made density matrix setup function for this case, but it should not be too difficult to implement:

You could e.g. start with the code of the PetscErrorCode System::DMWriteDiagThermal(Vec dm, PetscReal beta) function in src/sutil.cpp make a copy of that and then use index methods like inline PetscInt Index::IsMLSGroundState() and inline PetscInt Index::IsModeDensity() to check whether the index corresponds to the two-level system ground state and a diagonal element of the bosonic mode.

Then you just need to write the code for the coherent state distribution which is a poissonian if i remember correctly.

[modmido]

while ( index->ContinueLocal()) {
      if( index->IsModeDensity() && index->IsMLSGroundState() ) {
    value   = 1;
    i = num_mlsdims; // start with this index, this corresponds to the first bosonic mode
    while( i < num_dims ) {
      ierr = Energies(i,&energy); CHKERRQ(ierr);
      value *= here enter Poissonian as a function of the number state index of the bosonic mode: index->ModeQN(i)
      i += 2;
    }
    ierr    = VecSetValues(dm,1,&locindex,&value,INSERT_VALUES); CHKERRQ(ierr);
      }
      locindex = index->Increment();
    }

[modmido]

so if you replace the while loop in DMWriteDiagThermal with this one, and provide your coherent state population function then this should do the trick ^^

[modmido]

Hi there, did it work? If you want to, you could also open a pull request to include your initial state setup function to the repo. Anyways i would like to close the issue

[aleonau]

Hello. I introduced the changes and tried it but the results of calculation were not affected at all. I got the same curves as running the code without it. Probably this happened because I was running the particular version of the code from "example/ex1a" folder, which addresses some other functions. I decided to postpone it for a while to test some other modes. I ran the code for the vacuum state, which is both Fock and coherent at the same time with initially excited atomic states. In this regime the results of my calculation coincide with the ones obtained via "psiquasp". So I am happy to test my code and have at least this comparison :)

Anyway, thank you very much for the comments above! Probably, I will do one more attempt later. You can close the issue.

P.S. I also wanted to ask you about the "biggest" system you managed to consider. I tried running N = 100 (initially excited) with the field in the vacuum initial state and the following parameters:

ntls = 100, m0=100, dm0=10; qnumbers [5] = {100,0,0,0,0}; decay mode: only spontaneous emission with some given gamma (other modes disabled)

I got "out of memory" error , although I was running the code on the cluster. The requested memory size was insane...

[modmido]

Hi there, if the setup function is correct, then it should affect the results, at least for the initial state or time equal to zero.

Concerning memory: yes, also a polynomial scaling can explode quite quickly. However the part

//add the degrees of freedom
TLSAdd(ntls,ntls,ntls,tlsenergy);   
ModeAdd(m0+1,dm0,modeenergy);   

Allows for a lot of finetuning: dm0 is the number of offdiagonal terms in the bosonic mode, the second and third argument in TLSAdd are the offdiagonal elements in the TLS part of the density matrix, so setting this to e.g. 10 like

TLSAdd(ntls,10,10,tlsenergy);

reduces the storage and computation time from ~100^3=10^6 to~100*10*10 =10^4 e.g. by a factor of 100, which is a lot.

Generally it is always advisable to use offdiagonal truncation since the dephasing/damping of these elements scales with the distance to the diagonal in the exponent! Thus having a quantum coherence in a state like |100><0| is very, very unlikely to occur, that would be an extreme case of a very strongly coupled system and you should use different basis states altogether.

The number of offidagonal elements must be chosen depending on the system parameters. 10 is actually quite something, sometimes you need less.

[modmido]

Also I found that 100 TLS is actually already the regime where it is best to use the coherent state representation of the quantum master equation and the resulting Fokker-Planck Equation, probably best introduced in the Book of Carmichael:

H. J. Carmichael. Statistical Methods in Quantum Optics I: Master Equations and Fokker- Planck Equations. Springer, 2002.

The permutation symmetric method is best applied in the regime between 1 and 100. For one the system is very simple for 100 and above the system essentially behaves classically and should thus be treated with a classical equation, as expected from the quantum to classical transition.

[aleonau]

Allows for a lot of finetuning: dm0 is the number of offdiagonal terms in the bosonic mode, the second and third argument in TLSAdd are the offdiagonal elements in the TLS part of the density matrix, so setting this to e.g. 10 like

TLSAdd(ntls,10,10,tlsenergy);

reduces the storage and computation time from ~100^3=10^6 to~100*10*10 =10^4 e.g. by a factor of 100, which is a lot.

Oh, that's really a good advice to reduce the memory consumption. I will definitely try tuning it.

[aleonau]

Also I found that 100 TLS is actually already the regime where it is best to use the coherent state representation of the quantum master equation and the resulting Fokker-Planck Equation, probably best introduced in the Book of Carmichael:

H. J. Carmichael. Statistical Methods in Quantum Optics I: Master Equations and Fokker- Planck Equations. Springer, 2002.

The permutation symmetric method is best applied in the regime between 1 and 100. For one the system is very simple for 100 and above the system essentially behaves classically and should thus be treated with a classical equation, as expected from the quantum to classical transition.

Actually, this is almost the thing, which I am currently simulating :) That's why I need some solid results to compare with :) I am working in the positive P-representation, leading to the stochastic differential equations in the Ito form.

[modmido]

Happy to help. If you want to dive deeper into the usage of the library or generally the method, numerics on multi level system quantum master equations I recommend you to read the respective chapters in my dissertation:

https://api-depositonce.tu-berlin.de/server/api/core/bitstreams/859bbb2b-aa00-41fe-a7c8-ae1c8d241402/content

I tried my best to write it all down, since I left academia after that.

[aleonau]

Thank you!