# How to program mathematica (simple example) [closed]

How do I program mathematica to perform basic second quantization operation of creation and annihilation operators such as c† and commutation rules for c†c. I know that out there there are libraries (.m files) but I am looking for simple programming examples (.m files ?) so that I can do my own.

I am trying to create operators such as creation operator and annihilation operator and Symbol dagger. How to create new function and rules in mathematica. Also Green function advanced or (G^a) etc.

Just few and simple example

Thanks

## closed as too broad by MarcoB, m_goldberg, corey979, Feyre, ubpdqnFeb 24 '17 at 7:18

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

This is something I made mostly to learn the language rather than have something that I would use heavily. I don't claim it to be an efficient or even correct implementation of some simple quantum operators, but nevertheless I hope it can help.

I want to define a new head Operator, and a multiplication operation such that all the terms which do not contain the head Operator explicitly are factored out, leaving only the Operator objects inside. This could be done as follows.

Remove[ourTimes, Operator];

Format[Operator[a_]] := OverHat[a];
Format[ourTimes[a__]] := AngleBracket[Row@List[a]];

ourTimes[a___, b__, c___] /; FreeQ[{b}, Operator] :=
Times[b, ourTimes[a, c]];

ourTimes[a___, b_ + c__, d___] :=
ourTimes[a, b, d] +
ourTimes[a, Plus[c], d];
ourTimes[a___, b_.*(x_Operator + c__), d___] :=
ourTimes[a, b*x, d] + ourTimes[a, b*Plus[c], d];
ourTimes[a___, b_.*x_Operator + c__, d___] :=
ourTimes[a, b*x, d] + ourTimes[a, Plus[c], d];

ourTimes[a___, b__*x__Operator, d___] /; FreeQ[{b}, Operator] :=
Times[b, ourTimes[a, x, d]];

ourTimes[a___, ourTimes[b_, c__], d___] :=
ourTimes[a, b, c, d];

ourTimes[a___] /; FreeQ[{a}, Operator] := Times[a];

Attributes[ourTimes] = {Flat};


We can try it out by defining the following operators (results omitted, you can try these yourself).

q = (Operator[a] + Operator[a^\[Dagger]])/Sqrt[2 \[Omega]];
p = -I (Operator[a] - Operator[a^\[Dagger]])*Sqrt[\[Omega]/2];

ourTimes[q, p]
ourTimes[q, q] - Expand[ourTimes[q]*ourTimes[q]]


To get a numeric value as a final result, you could interpret expressions where ourTimes contains only operator arguments as expectation values, then define a list of rules that transform the expectation values. You could even define different rule sets for different states to get different results depending on which states you assume to have.

I made the example above when reading Leonid's book, section 4.3.6., which should also offer some insight into how the code above works. I recommend the book to you, it's a good read.