Message Boards Message Boards

0
|
5617 Views
|
12 Replies
|
1 Total Likes
View groups...
Share
Share this post:

Solve equation of motion with Dirac-fermions?

Posted 6 years ago

Dear Wolfram team:

I am a beginner of Mathematica.

My Problem is that I want solve a System with n equation of Motion in first order. In this equation of motion are creation and annihilation operators of Dirac-fermions. I don't know and don't find a it, how I can describes the creation and annihilation operators of Dirac-fermions in Mathematica. The equation of motion have the form:

$$\dot{c}_i^\dagger [t]=f*c_i^\dagger[t]+g[t]*c_{i+1}[t]-g[t]*c_{i+1}^\dagger[t]+h[t]*c_{i-1}[t]-h[t]*c_{i+1}^\dagger[t]\\ \dot{c}_i [t]=f*c_i[t]+g[t]*c_{i+1}^\dagger[t]-g[t]*c_{i+1}[t]-h[t]*c_{i-1}^\dagger[t]+h[t]*c_{i+1}[t],$$

where $c_i^\dagger,c_i $ are creation and annihilation operators and f,g,h are functions.

Then I want use DSolve or NDSolve to solve the equation of motion.

Thanks, for your help.

POSTED BY: C. H.
12 Replies
Posted 6 years ago

Thanks for your answer.

I found a lot of nice commands but Mathematica don't work. My code is for the equation are

eoMkitaevchaindiracnormal = {
   -I c1'[t] == -\[Omega] (c2[t] + ConjugateTranspose[c2[t]]) + 
     Exp[-((x1 - v*t)^2/(2 \[Sigma]^2))]/
      2 c1[t], -I ConjugateTranspose[c1]'[
      t] == \[Omega] (c2[t] + ConjugateTranspose[c2[t]]) + 
     Exp[-((x1 - v*t)^2/(2 \[Sigma]^2))]/
      2 ConjugateTranspose[c1[t]], -I c2'[
      t] == -\[Omega] (c1[t] - ConjugateTranspose[c1[t]]) + 
     Exp[-((x2 - v*t)^2/(2 \[Sigma]^2))]/
      2 c2[t], -I ConjugateTranspose[c2]'[
      t] == \[Omega] (-c1[t] + ConjugateTranspose[c1[t]]) + 
     Exp[-((x2 - v*t)^2/(2 \[Sigma]^2))]/2 ConjugateTranspose[c2[t]], ,
   c1[-100] == 1, c2[-100] == 1, ConjugateTranspose[c1][-100] == 1, 
   ConjugateTranspose[c2][-100] == 1};

Now, I want use NDSolve to solve the equation and use this code

solutionkitaevchaindiracnormal = 
         NDSolve[eoMkitaevchaindiracnormal /. {v0 -> 1, 
            v -> 1, \[Sigma] -> 0.75, x1 -> 1, x2 -> 2, \[Omega] -> 1}, {c1[
            t], c2[t], ConjugateTranspose[c1][t], 
           ConjugateTranspose[c2][t]}, {t, -100, 100}];

After the second step Mathematics give me this Output:

NDSolve::deqn: Equation or list of equations expected instead of Null in the first argument {-I c1^\[Prime](t)==1/2 E^(-0.888889 Plus[<<2>>]^2) c1(t)-c2(t)-c2(t)^\[ConjugateTranspose],-I ((c1^\[ConjugateTranspose])^\[Prime])\[InvisibleApplication](t)==c2(t)+1/2 E^(-0.888889 Plus[<<2>>]^2) c1(t)^\[ConjugateTranspose]+c2(t)^\[ConjugateTranspose],-I c2^\[Prime](t)==-c1(t)+1/2 E^(-0.888889 Plus[<<2>>]^2) c2(t)+c1(t)^\[ConjugateTranspose],-I ((c2^\[ConjugateTranspose])^\[Prime])\[InvisibleApplication](t)==-c1(t)+c1(t)^\[ConjugateTranspose]+1/2 E^(-0.888889 Plus[<<2>>]^2) c2(t)^\[ConjugateTranspose],Null,c1(-100)==1,c2(-100)==1,(c1^\[ConjugateTranspose])\[InvisibleApplication](-100)==1,(c2^\[ConjugateTranspose])\[InvisibleApplication](-100)==1}. >>

What is my mistake? Can I use ConjugateTranspose too describe a creation operator or I have use an other command ?

POSTED BY: C. H.

NDSolve is for functions, not operators.

POSTED BY: Frank Kampas
Posted 6 years ago

Ok, which command can I use for operators?

POSTED BY: C. H.

I don't think there is one in Mathematica. That's why I suggested FeynCalc.

POSTED BY: Frank Kampas
Posted 6 years ago

Ok, I have the package FeynCalc installed but I don't find in the list of FeynCalc an other command which solve my system of equations.

POSTED BY: C. H.

List subscribers can write emails by sending an email to feyncalcATfeyncalc.org (replacing "AT" with "@"). Others can use the link "mail a new topic".

from

https://www.feyncalc.org/forum/subject.html

POSTED BY: Frank Kampas

Umm, did you ever actually look at eoMkitaevchaindiracnormal? Because there is a Null in it, and the error message is telling you that, and you should find and remove it before appealing to an entire forum.

Next there is the business of double naming. Some places use ConjugateTranspose[c1[t]], others use ConjugateTranspose[c1][t]. I will propose getting rid of both and just using ctc1[t].NDSolve will handle the variant below.

eoMkitaevchaindiracnormal = {-I c1'[
     t] == -\[Omega] (c2[t] + ctc2[t]) + 
    Exp[-((x1 - v*t)^2/(2 \[Sigma]^2))]/2 c1[t], -I ctc1'[
     t] == \[Omega] (c2[t] + ctc2[t]) + 
    Exp[-((x1 - v*t)^2/(2 \[Sigma]^2))]/2 ctc1[t], -I c2'[
     t] == -\[Omega] (c1[t] - ctc1[t]) + 
    Exp[-((x2 - v*t)^2/(2 \[Sigma]^2))]/2 c2[t], -I ctc2'[
     t] == \[Omega] (-c1[t] + ctc1[t]) + 
    Exp[-((x2 - v*t)^2/(2 \[Sigma]^2))]/2 ctc2[t], c1[-100] == 1, 
  c2[-100] == 1, ctc1[-100] == 1, ctc2[-100] == 1}

Bottom line: this could and should have been debugged off-line. If the result from NDSolve is not satisfactory, that's the sort of issue to bring to the forum.

POSTED BY: Daniel Lichtblau
Posted 6 years ago

Thank,s for the hint with double naming. That was the problem

It is irrelevant to use ctc1[t] or ConjugateTranspose[c1][t].

When I solve and Plot the System with this code :

eoMkitaevchaindiracnormal = {
   -I c1'[t] == -\[Omega] (c2[t] + ConjugateTranspose[c2][t]) + 
     Exp[-((x1 - v*t)^2/(2 \[Sigma]^2))]/
      2 c1[t], -I ConjugateTranspose[c1]'[
      t] == \[Omega] (c2[t] + ConjugateTranspose[c2][t]) + 
     Exp[-((x1 - v*t)^2/(2 \[Sigma]^2))]/
      2 ConjugateTranspose[c1][t], -I c2'[
      t] == -\[Omega] (c1[t] - ConjugateTranspose[c1][t]) + 
     Exp[-((x2 - v*t)^2/(2 \[Sigma]^2))]/
      2 c2[t], -I ConjugateTranspose[c2]'[
      t] == \[Omega] (-c1[t] + ConjugateTranspose[c1][t]) + 
     Exp[-((x2 - v*t)^2/(2 \[Sigma]^2))]/2 ConjugateTranspose[c2][t],
   c1[-100] == 1, c2[-100] == 1, ConjugateTranspose[c1][-100] == 1, 
   ConjugateTranspose[c2][-100] == 1};
solutionkitaevchaindiracnormal = 
         NDSolve[eoMkitaevchaindiracnormal /. {v0 -> 1, 
            v -> 1, \[Sigma] -> 0.75, x1 -> 1, x2 -> 2, \[Omega] -> 1}, {c1[
            t], c2[t], ConjugateTranspose[c1][t], 
           ConjugateTranspose[c2][t]}, {t, -100, 100}]
ParametricPlot[{Re[#], Im[#]} &@
  Evaluate[{ConjugateTranspose[c1][t]*c1[t]} /. 
    solutionkitaevchaindiracnormal], {t, -10, 10}, PlotRange -> All, 
 PlotStyle -> {Blue,Thick}]

I become this solution for $c_1^\dagger c_1$

enter image description here

But $c_1^\dagger c_1$ is a density Operator.

I think that the imaginary part of $c_1^\dagger c_1$ should be 0 and

ConjugateTranspose[c1] and ctc1

is not $c_1^\dagger$. But what is the correct description of $c_1^\dagger$?

POSTED BY: C. H.
Posted 6 years ago

Thanks you for the hint with the double Name. Now Mathematica solve the System. I use ctc1[t} und ConjugateTranspose[c1[t]] and both have the same solution. But I think that Mathematica don't know that ConjugateTranspose[c1[t]] should be $c_1^\dagger$ and >Mathematica interpret c1 and ConjugateTranspose[c1[t]] as operator because if I plot $c_1^\dagger c_1$ then I get

enter image description here

where horizontal axis is the real part and vertical axis is the imaginary part. I think that ConjugateTranspose[c1[t]] is not $c_1^\dagger$. I think that the imaginary part are Zero. I use this Code:

eoMkitaevchaindiracnormal = {
   -I c1'[t] == -\[Omega] (c2[t] + ctc2[t]) + 
     v0 Exp[-((x1 - v*t)^2/(2 \[Sigma]^2))]/2 c1[t],
   -I ctc1'[t] == \[Omega] (c2[t] + ctc2[t]) + 
     v0 Exp[-((x1 - v*t)^2/(2 \[Sigma]^2))]/2 ctc1[t],
   -I c2'[t] == -\[Omega] (c1[t] - ctc1[t]) + 
     v0 Exp[-((x2 - v*t)^2/(2 \[Sigma]^2))]/2 c2[t],
   -I ctc2'[t] == \[Omega] (-c1[t] + ctc1[t]) + 
     v0 Exp[-((x2 - v*t)^2/(2 \[Sigma]^2))]/2 ctc2[t],
   c1[-0] == 1, c2[-0] == 1, ctc1[-0] == 1, ctc2[-0] == 1};
solutionkitaevchaindiracnormal = 
  NDSolve[eoMkitaevchaindiracnormal /. {v0 -> 0.1, 
     v -> 0, \[Sigma] -> 0.75, x1 -> 1, x2 -> 2, \[Omega] -> 1}, {c1[
     t], c2[t], ctc1[t], ctc2[t]}, {t, -0, 10}];
ParametricPlot[{Re[#], Im[#]} &@
           Evaluate[{ctc1[t]*c1[t]} /. solutionkitaevchaindiracnormal], {t, -0,
            10}, PlotRange -> All, 
          PlotStyle -> {{Blue, Thick}, {Black, Dashed, Thick}}]

When I include

ctc1 = ComplexConjugate[c1];
ctc2 = ComplexConjugate[c2];

then dont solve Mathematica the System.

What is my mistake and who can I describe the operators in Mathematica korrect?

POSTED BY: C. H.
Posted 6 years ago

The answer from https://www.feyncalc.org/forum/subject.html:

Hi, solving equations of motion is not what FeynCalc is used for. I also have no experience in that. Perhaps someone of mathematica.stackexchange.com can advice you a suitable package. Cheers, Vladyslav

Am 13.08.2018 um 18:14 schrieb Constantin:

Hey guys,

My Problem is that I want solve a System with n equation of Motion in first order. In this equation of motion are creation and annihilation operators of Dirac-fermions. I don't know and don't find a it, how I can describes the creation and annihilation operators of Dirac-fermions in Mathematica and then I want to solve the equation of motion.

I have asked the Question also by community.wolfram.com.

Thanks, for your help.

Other ideas to solve a system with operators of Dirac fermions?

POSTED BY: C. H.
Posted 6 years ago

I think I have it. My solution:

eoMkitaevchaindiracnormal = {
   -I c1'[t] == -\[Omega] (c2[t] + ctc2[t]) - 
     v0 Exp[-((x1 - v*t)^2/(2 \[Sigma]^2))]/2 c1[t],
   -I ctc1'[t] == \[Omega] (c2[t] + ctc2[t]) + 
     v0 Exp[-((x1 - v*t)^2/(2 \[Sigma]^2))]/2 ctc1[t],
   -I c2'[t] == -\[Omega] (c1[t] - ctc1[t]) - 
     v0 Exp[-((x2 - v*t)^2/(2 \[Sigma]^2))]/2 c2[t],
   -I ctc2'[t] == \[Omega] (-c1[t] + ctc1[t]) + 
     v0 Exp[-((x2 - v*t)^2/(2 \[Sigma]^2))]/2 ctc2[t],
   c1[-10] == 1, c2[-10] == 1, ctc1[-10] == 1, ctc2[-10] == 1};
solutionkitaevchaindiracnormal = 
  NDSolve[eoMkitaevchaindiracnormal /. {v0 -> 1, 
     v -> 1, \[Sigma] -> 0.75, x1 -> 1, x2 -> 2, \[Omega] -> 1}, {c1[
     t], c2[t], ctc1[t], ctc2[t]}, {t, -100, 100}];
Plot[Evaluate[{ctc1[t]*c1[t], ctc2[t]*c2[t], 
    ctc1[t]*c1[t] + ctc2[t]*c2[t]} /. 
   solutionkitaevchaindiracnormal], {t, -10, 10}, PlotRange -> All, 
 PlotStyle -> {{Blue, Thick}, {Orange, Thick}, {Black, Dashed, 
    Thick}}]

enter image description here

Thanks for your help.

Now, I have another problem. I want that the absolute value of $c_1, c_1^\dagger,c_2, c_2^\dagger$ is between 0 and 1. Assumptions and Assuming don't work with NDSolve. I know why I can't use Assumptions and Assuming. But I don't know what is the command for NDSolve?

POSTED BY: C. H.
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract