Author: Xiang Gao
For more information about QuTiP see http://qutip.org
If you have questions about this part, feel free to email qasdfgtyuiop@gmail.com.
This notebook correspond to the module "qutip.perturbation". And the source code of this module contains a lot of references towards formula and terms defined in this notebook. This notebook also contains a lot of the explanations on the implementation of that module. So if you make changes to this notebook, ALWAYS update the source code of that module to keep them consistent.
Here we discusses the theory of perturbation.
The theory of perturbation is very easy for non-degenerate states.
First of all, we often use the symbol like $\sum_i$. If the range of sum is omited (as in here), it means the sum is over all non-negative integers from $0$ to $+\infty$.
Assume the Hamiltonian of the system can be written as $$\hat{H} = \hat{H}_0+\lambda\hat{H}_1+\lambda^2\hat{H}_2+\cdots=\sum_i \lambda^i \hat{H}_i\tag{1.1.1.1}$$ Where $\hat{H}_0$ is the unperturbed Hamilton, $\hat{H}_i$is the $i^{th}$ order perturbation.
Let $E_n^{(0)}$ be the unperturbed eigen value of Hamiltonian with eigen state $\left|\Psi_n^{(0)}\right>$. The perturbed eigen value of Hamiltonian is called $E_n$, and the corresponding eigen state is $\left|\Psi_n\right>$, where $n=0,1,2,3,\ldots$.
We can write the Taylor expansion of $E_n$ as $$E_n = \sum_i\lambda^i E_n^{(i)}\tag{1.1.1.2}$$ and $\left|\Psi_n\right>$ as $$\left|\Psi_n\right>=\sum_i\lambda^i\left|\Psi_n^{(i)}\right>\tag{1.1.1.3}$$
We use the normalization condition $$\left<\Psi_k^{(0)}\middle|\Psi_l^{(0)}\right>=\delta_{kl}$$
From the difination above, we can write $$\hat{H}\left|\Psi_n\right>=E_n\left|\Psi_n\right>\tag{1.1.2.1}$$ and $$\hat{H}_0\left|\Psi_n^{(0)}\right>=E_n^{(0)}\left|\Psi_n^{(0)}\right>\tag{1.1.2.2}$$
Substitute $(1.1.1.1),(1.1.1.2),(1.1.1.3)$ into $(1.1.2.1)$, we get $$\left(\sum_i \lambda^i \hat{H}_i\right)\left(\sum_i\lambda^i\left|\Psi_n^{(i)}\right>\right)=\left(\sum_i\lambda^i E_n^{(i)}\right)\left(\sum_i\lambda^i\left|\Psi_n^{(i)}\right>\right)$$
Move everything to the left, we get $$\left[\sum_i\lambda^i\left(\hat{H}_i-E_n^{(i)}\right)\right]\left(\sum_i\lambda^i\left|\Psi_n^{(i)}\right>\right)=0$$
Expand this formula, we will get $$\left[\sum_{ij}\lambda^{i+j}\left(\hat{H}_i-E_n^{(i)}\right)\left|\Psi_n^{(j)}\right>\right]=0$$
Reorganize by setting $s=i+j$, we have $$\sum_s\left[\lambda^{s}\sum_{i=0}^s\left(\hat{H}_i-E_n^{(i)}\right)\left|\Psi_n^{(s-i)}\right>\right]=0$$
In order that the left side and the right side be the same series, we must have $$\sum_{i=0}^s\left(\hat{H}_i-E_n^{(i)}\right)\left|\Psi_n^{(s-i)}\right>=0\tag{1.1.2.3}$$
This equation contains the key information of $\hat{H}_i$, $E_n^{(i)}$ and $\left|\Psi_n^{(i)}\right>$. It is very important and we can retrieve everything we need from this equation.
$(1.1.2.3)$ can be written as $$\left(\hat{H}_0-E_n^{(0)}\right)\left|\Psi_n^{(s)}\right>+\left(\hat{H}_s-E_n^{(s)}\right)\left|\Psi_n^{(0)}\right>+\sum_{i=1}^{s-1}\left(\hat{H}_i-E_n^{(i)}\right)\left|\Psi_n^{(s-i)}\right>=0\tag{1.1.3.1}$$
Multiply $\left<\Psi_n^{(0)}\right|$ to the left, we get $$\left<\Psi_n^{(0)}\right|\hat{H}_s\left|\Psi_n^{(0)}\right>-E_n^{(s)}+\sum_{i=1}^{s-1}\left<\Psi_n^{(0)}\right|\left(\hat{H}_i-E_n^{(i)}\right)\left|\Psi_n^{(s-i)}\right>=0$$
So we get $E_n^{(s)}$ immediately $$E_n^{(s)}=\left<\Psi_n^{(0)}\right|\hat{H}_s\left|\Psi_n^{(0)}\right>+\sum_{i=1}^{s-1}\left<\Psi_n^{(0)}\right|\left(\hat{H}_i-E_n^{(i)}\right)\left|\Psi_n^{(s-i)}\right>\tag{1.1.3.2}$$
After getting $E_n^{(s)}$, we write $(1.1.3.1)$ as $$\left(\hat{H}_0-E_n^{(0)}\right)\left|\Psi_n^{(s)}\right>=-\sum_{i=1}^{s}\left(\hat{H}_i-E_n^{(i)}\right)\left|\Psi_n^{(s-i)}\right>\tag{1.1.3.3}$$
Our goal is to express $\left|\Psi_n^{(s)}\right>$ in the base $\left\{\left|\Psi_\alpha^{(0)}\right>\right\}$, to achieve this goal we write $\left|\Psi_n^{(s)}\right>$ as $$\left|\Psi_n^{(s)}\right>=\sum_\alpha C_{n\alpha}^{(s)}\left|\Psi_\alpha^{(0)}\right>\tag{1.1.3.4}$$
Substitute $(1.1.3.4)$ to $(1.1.3.3)$ and multiply $\left<\Psi_\beta^{(0)}\right|$ to the left, we get $$\left<\Psi_\beta^{(0)}\right|\left(\hat{H}_0-E_n^{(0)}\right)\left(\sum_\alpha C_{n\alpha}^{(s)}\left|\Psi_\alpha^{(0)}\right>\right)=-\sum_{i=1}^{s}\left<\Psi_\beta^{(0)}\right|\left(\hat{H}_i-E_n^{(i)}\right)\left|\Psi_n^{(s-i)}\right>\tag{1.1.3.5}$$
The left side of formula $(1.1.3.5)$ can be reorganized as $$\sum_{\alpha}\left<\Psi_\beta^{(0)}\right|\left(\hat{H}_0-E_n^{(0)}\right)\left|\Psi_\alpha^{(0)}\right>C_{n\alpha}^{(s)}=\sum_{\alpha}\left(E_\beta^{(0)}-E_n^{(0)}\right)C_{n\alpha}^{(s)}\delta_{\alpha\beta}=\left(E_\beta^{(0)}-E_n^{(0)}\right)C_{n\beta}^{(s)}$$
So we have $$\left(E_n^{(0)}-E_\beta^{(0)}\right)C_{n\beta}^{(s)}=\sum_{i=1}^{s}\left<\Psi_\beta^{(0)}\right|\left(\hat{H}_i-E_n^{(i)}\right)\left|\Psi_n^{(s-i)}\right>\tag{1.1.3.6}$$
For any $n\neq\beta$, we can get immediately $$C_{n\beta}^{(s)}=\frac{1}{E_n^{(0)}-E_\beta^{(0)}}\sum_{i=1}^{s}\left<\Psi_\beta^{(0)}\right|\left(\hat{H}_i-E_n^{(i)}\right)\left|\Psi_n^{(s-i)}\right>\tag{1.1.3.7}$$
In principle, the value of $C_{nn}^{(s)}$ is arbitrary, any value of $C_{nn}^{(s)}$ will give the correct answer. So we simply set $$C_{nn}^{(s)}=0\tag{1.1.3.8}$$
We have gotten the formula for non-degenerate states: we can use $(1.1.3.2)$ to calculate energy, and use $(1.1.3.7)$ and $(1.1.3.8)$ to calculate state. The alogrithm for degenerate states are similar but much more complicated compared with non-degenerate states.
The unperturbed Hamiltonian $\hat{H}_0$ have many eigenvalues, we call them $E^{(0)}_0, E^{(0)}_1, E^{(0)}_2, \ldots$. Each eigenvalue might be degenerate or not. If one eigenvalue is degenerate, then all the eigenvectors corresponding to this eigenvalue will form a subspace of the total Hilbert space of the states of the system, we call these spaces eigenspace of the corresponding eigenvalue. For a non-degenerate eigenvalue, the eigenspace is a 1 dimension subspace. The eigenspace of $E^{(0)}_i$ is denoted by $\mathscr{H^{(0)}_i}$. The dimension of $\mathscr{H^{(0)}_i}$ equals the degree of degeneracy of eigenvalue $E^{(0)}_i$.
After applying the first order perturbation, the eigenvalues and their eigenspaces might be splitted by the first order perturbation. In this case, the eigenvalue $E^{(0)}_i$, and its eigenspace $\mathscr{H^{(0)}_i}$ might split into several eigenvalues. We use $E^{(1)}_{i0}, E^{(1)}_{i1}, E^{(1)}_{i2}, \ldots$ to denote the change of eigenvalue (so the new eigenvalue will be $E^{(0)}_i+E^{(1)}_{i0}, E^{(0)}_i+E^{(1)}_{i1}, E^{(0)}_i+E^{(1)}_{i2}, \ldots$) due to first order perturbation. $\mathscr{H^{(0)}_i}$ will also be splitted, the new eigenspace will be denoted by $\mathscr{H^{(1)}_{i0}}, \mathscr{H^{(1)}_{i1}}, \mathscr{H^{(1)}_{i2}}, \ldots$. If there are still degeneracy after first order perturbation, second order perturbation might further split the eigenvalues and eigenspaces. This time, we use $E^{(2)}_{ij0}, E^{(2)}_{ij1}, E^{(2)}_{ij2}, \ldots$ and $\mathscr{H^{(2)}_{ij0}}, \mathscr{H^{(2)}_{ij1}}, \mathscr{H^{(2)}_{ij2}}, \ldots$ to denote the eigenvalue change of newly splitted eigenvalues and their newly splitted eigenspace. For higher order perturbation, this rule easily applies.
Assume the energy level is degerate up to the order $t$, the $t^{th}$ order degeneracy is $g_{t}(i_0i_1\cdots i_{t})$, that is: there is a $g_{t}(i_0i_1\cdots i_{t})$ dimension subspace $\mathscr{H^{(t)}_{i_0i_1\ldots i_t}}$ that all states in this space shares the same eigenvalues and their perturbation: $$E_{i_0}^{(0)}, E_{i_0i_1}^{(1)}, E_{i_0i_1i_2}^{(2)}, E_{i_0i_1i_2i_3}^{(3)}, \ldots, E_{i_0i_1\cdots i_{t}}^{(t)}$$ When $(t+1)^{th}$ order perturbation applies, $\mathscr{H^{(t)}_{i_0i_1\ldots i_t}}$ might split into several smaller subspaces $\mathscr{H^{(t+1)}_{i_0i_1\ldots i_ti_{t+1}}}$, where for any $i_{t+1}$, $$\mathscr{H^{(t+1)}_{i_0i_1\ldots i_ti_{t+1}}}\subseteq \mathscr{H^{(t)}_{i_0i_1\ldots i_t}}$$
Before calculating the $(t+1)^{th}$ order perturbation, we need first point out that how $\mathscr{H^{(t)}_{i_0i_1\ldots i_t}}$ is splitted into $\mathscr{H^{(t+1)}_{i_0i_1\ldots i_ti_{t+1}}}$.
For simplicity, we will now use a more compact notation by defining $n\equiv i_0i_1\cdots i_{t}$. We will also use $E_{n}^{(0)}, E_{n}^{(1)}, E_{n}^{(2)}, E_{n}^{(3)}, \ldots, E_{n}^{(t)}$ to denote $E_{i_0}^{(0)}, E_{i_0i_1}^{(1)}, E_{i_0i_1i_2}^{(2)}, E_{i_0i_1i_2i_3}^{(3)}, \ldots, E_{i_0i_1\cdots i_{t}}^{(t)}$. In the case of non-degenerate states, when we say $E_n^{(0)}$ we mean the $n^{th}$ eigenvalue of $\hat{H}_0$. $E_n^{(0)}$ is always different from $E_m^{(0)}$ unless $m=n$. This notation is unambiguous for the case of non-degenerate states. However, when there are degeneracy in energy levels, each of which might split after perturbation, the definition might be ambiguous. We may assign the each value of $E^{(0)}$ to several suffix $n$; in this case, different $n$ might correspond to the same value of $E^{(0)}$. Similarly for $E^{(1)}, E^{(2)},\ldots,E^{(t-1)}$.
Let $\left\{\left(\left|\Psi_{n\xi}^{(0)}\right>,\left|\Psi_{n\xi}^{(1)}\right>,\ldots,\left|\Psi_{n\xi}^{(t)}\right>\right)\middle| \xi=0,1,2,\ldots,g_t(n)\right\}$ be an arbitrary set of base of $\mathscr{H^{(t)}_{n}}$ and it's perturbation (like in the case of non-degenerate states). We are now trying to find the basis, which we call "suitable base" here, that will goes to different $\mathscr{H^{(t+1)}_{ni_{t+1}}}$. The suitable base will be linear combinations of the arbitrary picked base: $$\begin{cases} \left|\Phi_{n\eta}^{(0)}\right>=\sum_{\xi=0}^{g_{t-1}(n)}a_{\eta\xi}\left|\Psi_{n\xi}^{(0)}\right>\\ \left|\Phi_{n\eta}^{(1)}\right>=\sum_{\xi=0}^{g_{t-1}(n)}a_{\eta\xi}\left|\Psi_{n\xi}^{(1)}\right>\\ \ldots\\ \left|\Phi_{n\eta}^{(t-1)}\right>=\sum_{\xi=0}^{g_{t}(n)}a_{\eta\xi}\left|\Psi_{n\xi}^{(t)}\right> \end{cases}\tag{1.2.2.1}$$ $\left|\Phi_{n\eta}^{(0)}\right>$ must satisfy the normalization condition $\left<\Phi_{m\eta}^{(0)}\middle|\Phi_{n\xi}^{(0)}\right>=\delta_{mn}\delta_{\eta\xi}$
We write down $(1.1.2.3)$ for $\Phi$ with $s=t+1$ as $$\left(\hat{H}_0-E_n^{(0)}\right)\left|\Phi_{n\eta}^{(t+1)}\right>+\left(\hat{H}_{t+1}-E_{ni_{t+1}}^{(t+1)}\right)\left|\Phi_{n\eta}^{(0)}\right>+\sum_{i=1}^{t}\left(\hat{H}_i-E_n^{(i)}\right)\left|\Phi_{n\eta}^{(t+1-i)}\right>=0$$
Multiply $\left<\Psi_{n\chi}^{(0)}\right|$ to the left, we have $$\left<\Psi_{n\chi}^{(0)}\right|\left(\hat{H}_{t+1}-E_{ni_{t+1}}^{(t+1)}\right)\left|\Phi_{n\eta}^{(0)}\right>+\sum_{i=1}^{t}\left<\Psi_{n\chi}^{(0)}\right|\left(\hat{H}_i-E_n^{(i)}\right)\left|\Phi_{n\eta}^{(t+1-i)}\right>=0$$
Substitute $(1.2.2.1)$, we can get $$\sum_{\xi=0}^{g_{t}(n)}\left[\left<\Psi_{n\chi}^{(0)}\right|\hat{H}_{t+1}\left|\Psi_{n\xi}^{(0)}\right>+\sum_{i=1}^{t}\left<\Psi_{n\chi}^{(0)}\right|\left(\hat{H}_{i}-E_{n}^{(i)}\right)\left|\Psi_{n\xi}^{(t+1-i)}\right>-E_{ni_{t+1}}^{(t+1)}\delta_{\chi\xi}\right]a_{\eta\xi}=0\tag{1.2.2.2}$$
Let $$W_{\chi\xi}=\left<\Psi_{n\chi}^{(0)}\right|\hat{H}_{t+1}\left|\Psi_{n\xi}^{(0)}\right>+\sum_{i=1}^{t}\left<\Psi_{n\chi}^{(0)}\right|\left(\hat{H}_{i}-E_{n}^{(i)}\right)\left|\Psi_{n\xi}^{(t+1-i)}\right>\tag{1.2.2.3}$$ then $(1.2.2.2)$ can be written as $$\sum_{\xi=0}^{g_{t}(n)}\left[W_{\chi\xi}-E_{ni_{t+1}}^{(t+1)}\delta_{\chi\xi}\right]a_{\eta\xi}=0$$
This is a matrix equation $$\left(\left[\begin{array}{ccc} W_{00} & \cdots & W_{0g_{t}(n)}\\ \vdots & \ddots & \vdots\\ W_{g_{t}(n)0} & \cdots & W_{g_{t}(n)g_{t}(n)} \end{array}\right]-E_{ni_{t+1}}^{(t+1)}\left[\begin{array}{ccc} 1\\ & \ddots\\ & & 1 \end{array}\right]\right)\left[\begin{array}{c} a_{\eta0}\\ \vdots\\ a_{\eta g_{t}(n)} \end{array}\right]=0\tag{1.2.2.4}$$
Non-trival solution for this equation exists only when $$\det(W_{\chi\xi}-E_{ni_{t+1}}^{(t+1)}\delta_{\chi\xi})=0$$ $E_{ni_{t+1}}^{(t+1)}$ is the eigenvalues of matrix $W_{\chi\xi}$, and $a_{\eta\xi}$ is the eigenvectors.
We have just done two things: the first is to calculate the (splitted) $(t+1)^{th}$ order energy perturbation, the second is to point out the bases of $\mathscr{H^{(t+1)}_{ni_{t+1}}}$. Up to now, we know nothing about the $(t+1)^{th}$ order perturbation of states, and this is what we will do now.
Using our new notation for degenerate states, $(1.1.3.4)$ can be written as $$\left|\Phi_{n\xi}^{(s)}\right>=\sum_{m\alpha}C_{n\xi,m\alpha}^{(s)}\left|\Phi_{m\alpha}^{(0)}\right>$$
And $(1.1.3.6)$ can be written as $$\left(E_{n}^{(0)}-E_{m}^{(0)}\right)C_{n\xi,m\alpha}^{(s)}=\sum_{i=1}^{s}\left<\Phi_{m\alpha}^{(0)}\right|\left(\hat{H}_{i}-E_{n}^{(i)}\right)\left|\Phi_{n\xi}^{(s-i)}\right>$$
For $E_{n}^{(0)}\neq E_{m}^{(0)}$, we have $$C_{n\xi,m\alpha}^{(s)}=\frac{1}{E_{n}^{(0)}-E_{m}^{(0)}}\sum_{i=1}^{s}\left<\Phi_{m\alpha}^{(0)}\right|\left(\hat{H}_{i}-E_{n}^{(i)}\right)\left|\Phi_{n\xi}^{(s-i)}\right>\tag{1.2.3.1}$$
For $E_{n}^{(0)} = E_{m}^{(0)}$, we simply set $$C_{n\xi,m\alpha}^{(s)}=0\tag{1.2.3.2}$$
This section introduces how the perturbation theory is implemented in qutip.perturbation
In the perspective of data structure, we use forest, which is called 'energy level forest', to store energy level information. Trees in the energy level forest is called 'energy level tree'. The roots of energy level trees are $E_{i_0}^{(0)}$. Each $E_{i_0i_1\cdots i_{t}}^{(t)}$ correspond to one node in the energy level forest, the parent of $E_{i_0i_1\cdots i_{t}}^{(t)}$ is $E_{i_0i_1\cdots i_{t-1}}^{(t-1)}$.
We use list, which is called 'eigenspace list' to store information of the most splitted subspaces. When we say "most splitted subspaces" we means that if we have finished calculating $t^{th}$ order perturbation, in our datastructure we will store only $\mathscr{H^{(t)}_{i_0i_1\ldots i_t}}$ for all $i_0i_1\ldots i_t$. For any $k<t$, $\mathscr{H^{(k)}_{i_0i_1\ldots i_k}}$ can easily be calculated: $$\mathscr{H^{(k)}_{i_0i_1\ldots i_k}}=\bigcup_{i_{k+1}\ldots i_t}\mathscr{H^{(t)}_{i_0i_1\ldots i_t}}$$
Instances of this class are nodes in energy level forest. Each instance has the following attribute:
energy: the value of $E_{i_0i_1\cdots i_{s}}^{(s)}$.
degeneracy: the value of $g_s(i_0i_1\cdots i_{s})$.
parent: the reference to the parent node, $None$ if root.
children: a list of references to child nodes, empty if a leaf.
eigen_space: $None$ if a leaf, otherwise reference to $\mathscr{H^{(t)}_{i_0i_1\ldots i_t}}$.
The instance of this class stores one base vector and its perturbations, i.e. $\left[\left|\Psi_{i_0i_1\ldots i_t\xi}^{(0)}\right>, \left|\Psi_{i_0i_1\ldots i_t\xi}^{(1)}\right>,\ldots,\left|\Psi_{i_0i_1\ldots i_t\xi}^{(t)}\right>\right]$ for a specific $\xi$. The binary arithmetic operator $+-*/$ are overloaded as entrywise plus, entrywise minus, multiply by a number and divide by a number. Unitary operator $+-$ are also overloadded as entrywise.
Each instance has the following attribute:
vectors: the list $\left[\left|\Psi_{i_0i_1\ldots i_t\xi}^{(0)}\right>, \left|\Psi_{i_0i_1\ldots i_t\xi}^{(1)}\right>,\ldots,\left|\Psi_{i_0i_1\ldots i_t\xi}^{(t)}\right>\right]$, each element in the list is of type Qobj.
Instances of this class stores the information of $\mathscr{H^{(t)}_{i_0i_1\ldots i_t}}$. Each instance has the following attribute:
elp_node: the reference towards the corresponding leaf in energy level tree.
base_set: a list of bases(instance of qutip.perturbation.PerturbedBase).
There are also methods:
split(eigen_system_info): this method is often called during the calculation of a new order perturbation. in the case of degenerate state, new order perturbation might cause the eigenspace split into several parts. This method returns the newly splitted subspaces and update the sturcture of energy level forest.
energy_levels(): this method returns a list of energy corrections of different order perturbations corresponding to this subspace, i.e. the corresponding $\left[E_{i_0}^{(0)}, E_{i_0i_1}^{(1)}, E_{i_0i_1i_2}^{(2)}, E_{i_0i_1i_2i_3}^{(3)}, \ldots, E_{i_0i_1\cdots i_{t}}^{(t)}\right]$.
The instance of Perturbation serves as a calculator, so we need to create an instance of Perturbation for each job. When we create the instance of Perturbation we must give the Hamiltonian $\hat{H}_0$. We may also provide the eigenvalues and eigenvectors of $\hat{H}_0$. If the eigenvalues are not given but eigenstates are given, then eigenvalues will be calculated by $E = \frac{\left<\phi\middle|\hat{H}_0\middle|\phi\right>}{\left<\phi\middle|\phi\right>}$. If both the eigenvalues and eigenvectors are not given, the build-in eigenvalue calculator will automatically be called to calculate the eigenvalues and eigenvectors of $\hat{H}_0$.
Each instance of Perturbation has the following attributes:
hamiltonians: a list that stores the Hamiltonians of each order, i.e. $\left[\hat{H}_0, \hat{H}_1, \ldots\right]$.
order: the order of perturbation that we currently have calculated.
etol: the tolerance of energy, in practice, if $\left|E_m-E_n\right|<etol$, then we think the two energies are equal.
energy_level_trees: the energy level forest. it is a list of references to the root of each energy level tree.
eigen_spaces: is a list that store the information of the most splitted subspaces.
Perturbation(h0,zero_order_eigenstates=None,zero_order_energy_levels=None,etol=1E-8): This is the construction method. It initialize the instance of this class with H0 and its eigen values and vectors. It takes 4 parameters:
The first parameter is the unperturbed Hamiltonian $\hat{H}_0$.
The second parameter is the list of eigenstates of $\hat{H}_0$. This parameter has a default value $None$. If this value is $None$, then the build-in solver of qutip will be automatically called to solve for the eigen states and the third parameter will be ignored.
The third parameter is the list of energy levels of H0 corresponding to each value in the second parameter. This parameter has a default value $None$. If it's $None$ but the second parameter is not $None$, it will be automatically calculated by $E = \frac{\left<\phi\middle|\hat{H}_0\middle|\phi\right>}{\left<\phi\middle|\phi\right>}$.
The forth parameter $etol$ is the tolerance of energy, i.e. if $\left|E_m-E_n\right|<etol$, we regard the two energy are degenerate.
next_order(ht=(zero operator)): Calculate the $(self.order+1)^{th}$ order perturbation with the given ht. This method takes one parameter ht, it is the perturbation Hamiltonian of $t^{th}$ order, i.e. $\hat{H}_t$. It has a default value that is a zero operator.
categorize_eigenstates(eigenvalues,eigenstates): This is an inner method used by other methods of this class, it is not designed for user's usage. This method categorize eigenvalues and eigenvectors by the degeneracy. It returns such a data structure: $\begin{array}{ll} [ & [eigen\ value\ 1,[eigen\ vector\ 1,eigen\ vector\ 2,...]\ ],\\ & [eigen\ value\ 2,[eigen\ vector\ 1,eigen\ vector\ 2,...]\ ],\\ & \ldots\ ] \end{array}$
goto_converge(lamda,tol): Keep calculating higher order perturbation without any perturbation in Hamiltonian ($H_t=0$) until both (perturbed) eigenvalues and eigenstates converge. The parameter lamda is the $\lambda$ in the formula $H=H_0+\lambda H_1+\cdots$. The parameter tol is the tolerance used as criteria of converge, with default value the "etol" specified when the instance of the class Perturbation is created.
result(lamda): A method similiar to Qobj.eigenstates() except that the returned values are the values calculated by perturbation theory. The parameter lamda is a float number which is the $\lambda$ in the formula $H=H_0+\lambda H_1+\cdots$. This method have two return values: The first return value, eigvals, is an array of eigenvalues for operator $H=H_0+\lambda H_1+\cdots$ calculated by perturbation theory. The second return value, eigvecs, is an array of eigenkets for operator $H=H_0+\lambda H_1+\cdots$ calculated by perturbation theory.
Before we begin, we must import modules needed.
In [1]:
from qutip import *
from numpy import *
The Hamiltonian without perturbation is $$\hat{H}_{0}=\left[\begin{array}{ccc} 1\\ & 1\\ & & 1 \end{array}\right]$$ The first order perturbation is $$\hat{H}_{1}=\left[\begin{array}{ccc} 1\\ & 0\\ & & -1 \end{array}\right]$$ The total hamiltonian is $$\hat{H}=\hat{H}_0+\lambda\hat{H}_1$$ We chooose the base $$\left|\Psi_{1}^{(0)}\right\rangle =\frac{1}{\sqrt{2}}\left[\begin{array}{c} 1\\ 1\\ 0 \end{array}\right]$$ $$\left|\Psi_{2}^{(0)}\right\rangle =\frac{1}{\sqrt{2}}\left[\begin{array}{c} 1\\ -1\\ 0 \end{array}\right]$$ $$\left|\Psi_{3}^{(0)}\right\rangle =\left[\begin{array}{c} 0\\ 0\\ 1 \end{array}\right]$$ These three states are degenerate and have the same energy $E^{(0)}=1$
Now we do the calculation:
In [2]:
H0 = Qobj([[1,0,0],[0,1,0],[0,0,1]])
H1 = Qobj([[1,0,0],[0,0,0],[0,0,-1]])
base1 = Qobj([[1],[1],[0]])/sqrt(2)
base2 = Qobj([[1],[-1],[0]])/sqrt(2)
base3 = Qobj([[0],[0],[1]])
calculator = Perturbation(H0,[base1,base2,base3])
calculator.next_order(H1)
The first order corrections to energy level (with unit $\lambda$) are
In [3]:
for i in calculator.energy_level_trees[0].children:
print(i.energy)
The base states of each (splitted due to perturbation) subspace and their first order corrections are
In [4]:
for i in calculator.eigen_spaces:
for base in i.base_set:
j = base.vectors
print('zeroth order:\n',j[0].full(),'\n\nfirst order:\n',j[1].full(),'\n\n-------------------')
So the final result, considering only the zeroth and first order perturbation, with $\lambda=0.1$, is:
In [5]:
energylevels,states = calculator.result(0.1)
print("energy levels:",energylevels)
print("eigenstates:",states)
Now we let it keep calculating higher and higher order perturbation until converge:
In [6]:
calculator.goto_converge(0.1)
The final result is:
In [7]:
energylevels,states = calculator.result(0.1)
print("energy levels:",energylevels)
print("eigenstates:",states)
From the result we can see that, the result of first order perturbation is already the exact solution.
The Hamiltonian without perturbation is $$\hat{H}_{0}=\left[\begin{array}{ccc} 25739.860\\ & 50266.917 \end{array}\right]$$ The first order perturbation is $$\hat{H}_{1}=\left[\begin{array}{ccc} & 100\\ 100 & \end{array}\right]$$ The total hamiltonian is $$\hat{H}=\hat{H}_0+\lambda\hat{H}_1$$ We chooose the base $$\left|\Psi_{1}^{(0)}\right\rangle =\left[\begin{array}{c} 1\\ 0 \end{array}\right]$$ $$\left|\Psi_{2}^{(0)}\right\rangle =\left[\begin{array}{c} 0\\ 1 \end{array}\right]$$ These two states have energy $$E_1^{(0)}=25739.860$$ $$E_2^{(0)}=50266.917$$
Now we do the calculation:
In [8]:
H0 = Qobj([[25739.860,0],[0,50266.917]])
H1 = Qobj([[0,100],[100,0]])
calculator = Perturbation(H0)
calculator.next_order(H1)
calculator.next_order()
The first order corrections to energy level (with unit $\lambda$) are
In [9]:
for i in calculator.energy_level_trees:
for j in i.children:
print(j.energy)
The second order corrections to energy level (with unit $\lambda^2$) are
In [10]:
for i in calculator.energy_level_trees:
for j in i.children:
for k in j.children:
print(k.energy)
The base states of each (splitted due to perturbation) subspace and their first order corrections are
In [11]:
for i in calculator.eigen_spaces:
for base in i.base_set:
j = base.vectors
print('zeroth order:\n',j[0].full(),'\n\nfirst order:\n',j[1].full(),'\n\n-------------------')
So the final result, considering only the zeroth, first and second order perturbation, with $\lambda=1$, is:
In [12]:
energylevels,states = calculator.result(1)
print("energy levels:",energylevels)
print("eigenstates:",states)
Now we let it keep calculating higher and higher order perturbation until converge:
In [13]:
calculator.goto_converge(1)
The final result is:
In [14]:
energylevels,states = calculator.result(1)
print("energy levels:",energylevels)
print("eigenstates:",states)
Everytihing are the same as Question 2 excpet that now $$\hat{H}_{1}=\left[\begin{array}{ccc} & 5000\\ 5000 & \end{array}\right]$$
Now we do the calculation
In [15]:
H0 = Qobj([[25739.860,0],[0,50266.917]])
H1 = Qobj([[0,5000],[5000,0]])
calculator = Perturbation(H0)
calculator.next_order(H1)
calculator.next_order()
First order corrections (of unit $\lambda$) of energy are now:
In [16]:
for i in calculator.energy_level_trees:
for j in i.children:
print(j.energy)
Second order corrections (of unit $\lambda^2$) of energy are now:
In [17]:
for i in calculator.energy_level_trees:
for j in i.children:
for k in j.children:
print(k.energy)
The base states of each (splitted due to perturbation) subspace and their first order corrections are
In [18]:
for i in calculator.eigen_spaces:
for base in i.base_set:
j = base.vectors
print('zeroth order:\n',j[0].full(),'\n\nfirst order:\n',j[1].full(),'\n\n-------------------')