Functions.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import sympy\n",
    "from scipy import constants, stats, linalg"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple Python tools for modelling few-level atom-light interactions\n",
    "\n",
    "This notebook outlines all the functions used by the `OBE_Tools` module to find solutions of the optical Bloch equations for an arbitrary number of levels both within and beyond the weak-probe regime. Different solution methods are available for different situations. We also provide some functions for looking at the time evolution of the system, and performing Doppler averaging in thermal ensembles. Examples of using the functions outlined here can be found in the \"Examples\" notebook. More information on the physics can be found in the paper.\n",
    "\n",
    "\n",
    "### Notes\n",
    "\n",
    "- $\\hbar = 1$\n",
    "- Parameter/variable names are chosen to match the notation in the paper. Parameters referring to fields that couple states $i$ and $j$ have the subscript $ij$, for example $\\Omega_{ij}$, `Omega_ij`. Parameters that are related to an individual atomic state $|n\\rangle$ have a single subscript $n$ such as $\\Gamma_n$, `Gamma_n`. Variables with upper/lower-case letters refer to parameters denoted by upper/lower-case Greek letters respectively, for example `Gamma` denotes $\\Gamma$ whereas `gamma` denotes $\\gamma$.\n",
    "\n",
    "The code makes a number of assumptions:\n",
    "- That the number of atomic levels is equal to 1 plus the number of fields. This is consistent with a ladder system in which each level is coupled to the ones above and below it (except for the top and bottom levels).\n",
    "- That each level only decays to the level below it, and there is no decay out of the bottom level. \n",
    "- That laser linewidths ($\\gamma$) do not matter (unless they are explicitly specified).\n",
    "\n",
    "\n",
    "# Introduction\n",
    "\n",
    "Often in modelling dynamics of atom-light systems it is useful to be able to numerically solve the optical Bloch equations (OBEs), especially since for systems with more than 2 levels the analytic solutions are difficult to derive and require significant approximations. The OBEs arise from using the Master equation to find the time dependence of the density matrix through\n",
    "\n",
    "$$\\frac{\\partial \\hat{\\rho}}{\\partial t} = -\\frac{i}{\\hbar}\\left[\\hat{H}, \\hat{\\rho}\\right] + \\hat{\\mathcal{L}}(\\hat{\\rho})$$\n",
    "\n",
    "where $\\hat{H}$ is the Hamiltonian of the atom-light system, $\\hat{\\rho}$ is the density matrix and $\\hat{\\mathcal{L}}$ describes the decay/dephasing in the system.\n",
    "\n",
    "\n",
    "\n",
    "# Steady-state solutions\n",
    "\n",
    " \n",
    "\n",
    "## Matrix Solution\n",
    "\n",
    "This solution method uses `Sympy` to derive the OBEs from an interaction Hamiltonian and decay/dephasing operator, then expresses these as matrix equations and uses a singular value decomposition to solve them. It makes no assumptions about the relative strengths of the field so is valid both within and beyond the weak-probe regime.\n",
    "\n",
    "First we set up two functions to create a list of the components of the density matrix $\\rho_{i,j}$ and then set up the density matrix $\\hat{\\rho}$ itself. \n",
    "These functions return a list/matrix populated with `Sympy` symbolic objects which will allow us to set up the equations and manipulate the terms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_rho_list(levels = 3):\n",
    "    rho_list = []\n",
    "    for i in range(levels):\n",
    "        for j in range(levels):\n",
    "            globals()['rho_'+str(i+1)+str(j+1)] = sympy.Symbol(\\\n",
    "            'rho_'+str(i+1)+str(j+1))\n",
    "            rho_list.append(globals()['rho_'+str(i+1)+str(j+1)])\n",
    "    return rho_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[rho_11, rho_12, rho_13, rho_21, rho_22, rho_23, rho_31, rho_32, rho_33]\n",
      "Type: <class 'list'>\n",
      "Length: 9\n"
     ]
    }
   ],
   "source": [
    "print(create_rho_list())\n",
    "print('Type: {}'.format(type(create_rho_list())))\n",
    "print('Length: {}'.format(len(create_rho_list())))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For an n-level system the density matrix will have $n^2$ components. For a 2-level system with ground state $|1\\rangle$ and excited state $|2\\rangle$ these components would be $\\rho_{11}, \\rho_{12}, \\rho_{21}$ and $\\rho_{22}$. The usefulness of creating a list instead of the properly structured $n\\times n$ matrix is purely for ease of referring to particular elements later in the code.\n",
    "\n",
    "Similarly if we put these components into the usual density matrix structure we get $$\\hat{\\rho} = \\begin{pmatrix} \\rho_{11} & \\rho_{12}\\\\ \\rho_{21} & \\rho_{22} \\end{pmatrix}$$\n",
    "\n",
    "Again this matrix is populated with Sympy symbolic objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_rho_matrix(levels = 3):\n",
    "    rho_matrix = numpy.empty((levels, levels), dtype = 'object')\n",
    "    for i in range(levels):\n",
    "        for j in range(levels):\n",
    "            globals()['rho_'+str(i+1)+str(j+1)] = \\\n",
    "            sympy.Symbol('rho_'+str(i+1)+str(j+1))\n",
    "            rho_matrix[i,j] = globals()['rho_'+str(i+1)+str(j+1)]\n",
    "    return numpy.matrix(rho_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[rho_11 rho_12 rho_13]\n",
      " [rho_21 rho_22 rho_23]\n",
      " [rho_31 rho_32 rho_33]]\n",
      "Type: <class 'numpy.matrix'>\n",
      "Shape: (3, 3)\n"
     ]
    }
   ],
   "source": [
    "print(create_rho_matrix())\n",
    "print('Type: {}'.format(type(create_rho_matrix())))\n",
    "print('Shape: {}'.format(create_rho_matrix().shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we set up the Hamiltonian of the interaction between the atom and the light. If we have a simple 2-level system where the levels are coupled by a field with Rabi frequency $\\Omega_{12}$ which is detuned from resonance by $\\Delta_{12}$ then (ignoring $\\hbar$) our Hamiltonian will look like \n",
    "$$\\hat{H} = \\frac{1}{2} \\begin{pmatrix} 0 & \\Omega_{12} \\\\ \\Omega_{12} & -2\\Delta_{12} \\end{pmatrix}$$\n",
    "\n",
    "Similarly for a 3-level system where levels $|1\\rangle$ and $|2\\rangle$ are coupled by $\\Omega_{12}$ with detuning $\\Delta_{12}$ and levels $|2\\rangle$ and $|3\\rangle$ are coupled by $\\Omega_{23}$ with detuning $\\Delta_{23}$, our Hamiltonian will be\n",
    "$$ \\hat{H} = \\frac{1}{2} \\begin{pmatrix} 0 & \\Omega_{12} & 0 \\\\ \\Omega_{12} & -2\\Delta_{12} & \\Omega_{23} \\\\ 0 & \\Omega_{23} & -2(\\Delta_{12} + \\Delta_{23}) \\end{pmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Hamiltonian(Omegas, Deltas):\n",
    "    levels = len(Omegas)+1\n",
    "    H = numpy.zeros((levels,levels))\n",
    "    for i in range(levels):\n",
    "        for j in range(levels):\n",
    "            if numpy.logical_and(i==j, i!=0):\n",
    "                H[i,j] = -2*(numpy.sum(Deltas[:i]))\n",
    "            elif numpy.abs(i-j) == 1:\n",
    "                H[i,j] = Omegas[numpy.min([i,j])]\n",
    "    return numpy.matrix(H/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function will set up the Hamiltonian for any n-level ladder system where each level is only coupled to the levels directly above and below it (ladder system). The values of the Rabi frequencies $\\Omega_{ij}$ and detunings $\\Delta_{ij}$ of each field must be supplied as lists/arrays of the same length (length $n - 1$ for an n-level system). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.   0.5  0. ]\n",
      " [ 0.5 -2.   1.5]\n",
      " [ 0.   1.5 -6. ]]\n",
      "Type: <class 'numpy.matrix'>\n",
      "Shape: (3, 3)\n"
     ]
    }
   ],
   "source": [
    "Omegas = [1,3]\n",
    "Deltas = [2,4]\n",
    "H = Hamiltonian(Omegas, Deltas)\n",
    "print(H)\n",
    "print('Type: {}'.format(type(H)))\n",
    "print('Shape: {}'.format(H.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Decay and Dephasing\n",
    "\n",
    "Next we set up the matrix form of the decay operator $\\hat{\\mathcal{L}}$. We will split this into two separate parts for simplicity: one part describing the atomic decay $\\hat{L}_{\\rm{atom}}$ and another describing dephasing due to the finite linewidths of the driving fields $\\hat{L}_{\\rm{dephasing}}$. The total operator will then be the sum of these parts $\\hat{\\mathcal{L}} = \\hat{L}_{\\rm{atom}} + \\hat{L}_{\\rm{dephasing}}$.\n",
    "\n",
    "First we focus on the Linblad superoperator $\\hat{L}_{\\rm{atom}}$ describing spontaneous decay. For this simple approach we assume that each level can only decay to the level below it at a rate $\\Gamma_n$, and that there is no decay out of the bottom level. This means that the coherences (off-diagonal elements) will be as $L_{ij,\\, i\\neq j} = -\\frac{\\Gamma_i + \\Gamma_j}{2}\\rho_{ij}$ while the populations (diagonal elements) will be given by $L_{ij, \\, i=j} = \\Gamma_{i+1}\\rho_{i+1, j+1} - \\Gamma_{i}\\rho_{ij}$. This means that for a 3-level system we have\n",
    "\n",
    "$$\\hat{L}_{\\rm{atom}} = \\begin{pmatrix} \n",
    "\\Gamma_2 \\rho_{22} - \\Gamma_1\\rho_{11} & - \\frac{\\Gamma_1 + \\Gamma_2}{2}\\rho_{12} & -\\frac{\\Gamma_1 + \\Gamma_3}{2}\\rho_{13} \\\\ \n",
    "-\\frac{\\Gamma_2 + \\Gamma_1}{2}\\rho_{21} & \\Gamma_3 \\rho_{33} - \\Gamma_2\\rho_{22} & -\\frac{\\Gamma_2 + \\Gamma_3}{2}\\rho_{23} \\\\\n",
    "-\\frac{\\Gamma_3 + \\Gamma_1}{2}\\rho_{31} & -\\frac{\\Gamma_3 + \\Gamma_2}{2}\\rho_{32} &  - \\Gamma_3\\rho_{33} \\end{pmatrix}$$\n",
    "Making the assumption that there is no decay out of the lower level ($\\Gamma_1 = 0$) simplifies the expression to \n",
    "$$\\hat{L}_{\\rm{atom}} = \\begin{pmatrix} \n",
    "\\Gamma_2 \\rho_{22} & - \\frac{\\Gamma_2}{2}\\rho_{12} & -\\frac{\\Gamma_3}{2}\\rho_{13} \\\\ \n",
    "-\\frac{\\Gamma_2}{2}\\rho_{21} & \\Gamma_3 \\rho_{33} - \\Gamma_2\\rho_{22} & -\\frac{\\Gamma_2 + \\Gamma_3}{2}\\rho_{23} \\\\\n",
    "-\\frac{\\Gamma_3}{2}\\rho_{31} & -\\frac{\\Gamma_3 + \\Gamma_2}{2}\\rho_{32} &  - \\Gamma_3\\rho_{33} \\end{pmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def L_decay(Gammas):\n",
    "    levels = len(Gammas)+1\n",
    "    rhos = create_rho_matrix(levels = levels)\n",
    "    Gammas_all = [0] + Gammas\n",
    "    decay_matrix = numpy.zeros((levels, levels), dtype = 'object')\n",
    "    for i in range(levels):\n",
    "        for j in range(levels):\n",
    "            if i != j:\n",
    "                decay_matrix[i,j] = -0.5*(Gammas_all[i]\\\n",
    "                                          +Gammas_all[j])*rhos[i,j]\n",
    "            elif i != levels - 1:\n",
    "                into = Gammas_all[i+1]*rhos[1+i, j+1]\n",
    "                outof = Gammas_all[i]*rhos[i, j]\n",
    "                decay_matrix[i,j] = into - outof\n",
    "            else:\n",
    "                outof = Gammas_all[i]*rhos[i, j]\n",
    "                decay_matrix[i,j] = - outof\n",
    "    return numpy.matrix(decay_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2*rho_22 -1.0*rho_12 -0.5*rho_13]\n",
      " [-1.0*rho_21 -2*rho_22 + rho_33 -1.5*rho_23]\n",
      " [-0.5*rho_31 -1.5*rho_32 -rho_33]]\n",
      "Type: <class 'numpy.matrix'>\n",
      "Shape: (3, 3)\n"
     ]
    }
   ],
   "source": [
    "Gammas = [2,1]\n",
    "L_atom = L_decay(Gammas)\n",
    "print(L_atom)\n",
    "print('Type: {}'.format(type(L_atom)))\n",
    "print('Shape: {}'.format(L_atom.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again this function returns an $n\\times n$ matrix populated by multiples of Sympy symbolic objects. \n",
    "\n",
    "The dephasing due to the laser linewidths is generally a smaller effect, so can often be neglected. If we do include it, it only affects the coherences (off-diagonal density matrix elements). Following the treatment in Gea-Banacloche et al, the elements of the dephasing matrix can be expressed as\n",
    "$$L_{ij, i\\neq j} = -\\left( \\gamma_{i,i+1} + \\dots + \\gamma_{j-1, j}\\right) \\rho_{ij}$$ where $\\gamma_{n,n+1}$ is the linewidth of the field coupling the $n$ and $n+1$ levels."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def L_dephasing(gammas):\n",
    "    levels = len(gammas)+1\n",
    "    rhos = create_rho_matrix(levels = levels)\n",
    "    deph_matrix = numpy.zeros((levels, levels), dtype = 'object')\n",
    "    for i in range(levels):\n",
    "        for j in range(levels):\n",
    "            if i != j:\n",
    "                deph_matrix[i,j] = -(numpy.sum(\\\n",
    "                gammas[numpy.min([i,j]):numpy.max([i,j])]))*rhos[i,j]\n",
    "    return numpy.matrix(deph_matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 -rho_12 -2*rho_13]\n",
      " [-rho_21 0 -rho_23]\n",
      " [-2*rho_31 -rho_32 0]]\n",
      "Type: <class 'numpy.matrix'>\n",
      "Shape: (3, 3)\n"
     ]
    }
   ],
   "source": [
    "gammas = [1,1]\n",
    "L_laser = L_dephasing(gammas)\n",
    "print(L_laser)\n",
    "print('Type: {}'.format(type(L_laser)))\n",
    "print('Shape: {}'.format(L_laser.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have expressions for the Hamiltonian and Linbladian, we can put them into the Master equation to get an expression for the time evolution of the density matrix. We have that \n",
    "$$ \\frac{\\partial \\hat{\\rho}}{\\partial t} = -i(\\hat{H}\\hat{\\rho} - \\hat{\\rho}\\hat{H}) + \\hat{L}_{\\rm{atom}} + \\hat{L}_{\\rm{dephasing}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Master_eqn(H_tot, L):\n",
    "    levels = H_tot.shape[0]\n",
    "    rho = create_rho_matrix(levels = levels)\n",
    "    return -1j*(H_tot*rho - rho*H_tot) + L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2*rho_22 - 1.0*I*(-0.5*rho_12 + 0.5*rho_21)\n",
      "  -2.0*rho_12 - 1.0*I*(-0.5*rho_11 + 2.0*rho_12 - 1.5*rho_13 + 0.5*rho_22)\n",
      "  -2.5*rho_13 - 1.0*I*(-1.5*rho_12 + 6.0*rho_13 + 0.5*rho_23)]\n",
      " [-2.0*rho_21 - 1.0*I*(0.5*rho_11 - 2.0*rho_21 - 0.5*rho_22 + 1.5*rho_31)\n",
      "  -2*rho_22 + rho_33 - 1.0*I*(0.5*rho_12 - 0.5*rho_21 - 1.5*rho_23 + 1.5*rho_32)\n",
      "  -2.5*rho_23 - 1.0*I*(0.5*rho_13 - 1.5*rho_22 + 4.0*rho_23 + 1.5*rho_33)]\n",
      " [-2.5*rho_31 - 1.0*I*(1.5*rho_21 - 6.0*rho_31 - 0.5*rho_32)\n",
      "  -2.5*rho_32 - 1.0*I*(1.5*rho_22 - 0.5*rho_31 - 4.0*rho_32 - 1.5*rho_33)\n",
      "  -rho_33 - 1.0*I*(1.5*rho_23 - 1.5*rho_32)]]\n",
      "Type: <class 'numpy.matrix'>\n",
      "Shape: (3, 3)\n"
     ]
    }
   ],
   "source": [
    "L = L_atom + L_laser\n",
    "master_rhs = Master_eqn(H,L)\n",
    "print(master_rhs)\n",
    "print('Type: {}'.format(type(master_rhs)))\n",
    "print('Shape: {}'.format(master_rhs.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have an expression for the time evolution of each of the components of the density matrix in terms of other components. \n",
    "To make this complex system of equations more convenient to solve we want to write them in the form \n",
    "$$\\frac{\\partial \\hat{\\rho}_{\\rm{vect}}}{\\partial t} = \\hat{M}\\hat{\\rho}_{\\rm{vect}}$$\n",
    "where $\\hat{\\rho}_{\\rm{vect}}$ is a column vector containing all of the elements of $\\hat{\\rho}$ in the order returned by `create_rho_list`, and $\\hat{M}$ is a matrix of coeficients. \n",
    "\n",
    "We can create this matrix of coefficients $\\hat{M}$ by using the `Sympy` functions `expand` and `coeff`. The `expand` function collects together terms containing the same symbolic object. The `coeff` function allows us to extract the multiplying coefficient for a specific symbolic object. (This works since the only terms we have are linear with respect to the density matrix elements, e.g. we never have terms in $\\rho^2$ or higher)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def OBE_matrix(Master_matrix):\n",
    "    levels = Master_matrix.shape[0]\n",
    "    rho_vector = create_rho_list(levels = levels)\n",
    "    coeff_matrix = numpy.zeros((levels**2, levels**2), dtype = 'complex')\n",
    "    count = 0\n",
    "    for i in range(levels):\n",
    "        for j in range(levels):\n",
    "            entry = Master_matrix[i,j]\n",
    "            expanded = sympy.expand(entry)\n",
    "            for n,r in enumerate(rho_vector):\n",
    "                coeff_matrix[count, n] = complex(expanded.coeff(r))\n",
    "            count += 1\n",
    "    return coeff_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0. +0.j   0. +0.5j  0. +0.j   0. -0.5j  2. +0.j   0. +0.j   0. +0.j\n",
      "   0. +0.j   0. +0.j ]\n",
      " [ 0. +0.5j -2. -2.j   0. +1.5j  0. +0.j   0. -0.5j  0. +0.j   0. +0.j\n",
      "   0. +0.j   0. +0.j ]\n",
      " [ 0. +0.j   0. +1.5j -2.5-6.j   0. +0.j   0. +0.j   0. -0.5j  0. +0.j\n",
      "   0. +0.j   0. +0.j ]\n",
      " [ 0. -0.5j  0. +0.j   0. +0.j  -2. +2.j   0. +0.5j  0. +0.j   0. -1.5j\n",
      "   0. +0.j   0. +0.j ]\n",
      " [ 0. +0.j   0. -0.5j  0. +0.j   0. +0.5j -2. +0.j   0. +1.5j  0. +0.j\n",
      "   0. -1.5j  1. +0.j ]\n",
      " [ 0. +0.j   0. +0.j   0. -0.5j  0. +0.j   0. +1.5j -2.5-4.j   0. +0.j\n",
      "   0. +0.j   0. -1.5j]\n",
      " [ 0. +0.j   0. +0.j   0. +0.j   0. -1.5j  0. +0.j   0. +0.j  -2.5+6.j\n",
      "   0. +0.5j  0. +0.j ]\n",
      " [ 0. +0.j   0. +0.j   0. +0.j   0. +0.j   0. -1.5j  0. +0.j   0. +0.5j\n",
      "  -2.5+4.j   0. +1.5j]\n",
      " [ 0. +0.j   0. +0.j   0. +0.j   0. +0.j   0. +0.j   0. -1.5j  0. +0.j\n",
      "   0. +1.5j -1. +0.j ]]\n",
      "Type: <class 'numpy.ndarray'>\n",
      "Shape: (9, 9)\n"
     ]
    }
   ],
   "source": [
    "coeffs = OBE_matrix(master_rhs)\n",
    "print(coeffs)\n",
    "print('Type: {}'.format(type(coeffs)))\n",
    "print('Shape: {}'.format(coeffs.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Singular Value Decomposition\n",
    "\n",
    "We now have an expression for $\\frac{\\partial \\hat{\\rho}_{\\rm{vect}}}{\\partial t}$ in terms of the matrix $\\hat{M}$ multiplied by $\\hat{\\rho}_{\\rm{vect}}$. In the steady state, $\\frac{\\partial \\hat{\\rho}_{\\rm{vect}}}{\\partial t} = 0$ so we need to solve the expression \n",
    "$$\\hat{M}\\hat{\\rho}_{\\rm{vect}} = 0.$$\n",
    "\n",
    "This can be solved numerically by performing a singular value decomposition (SVD) on the matrix $\\hat{M}$. The full details of the theory of the SVD are beyond this discussion, it is suffice to say that the SVD transforms $\\hat{M}$ into three matrices such that $$\\hat{M} = U\\Sigma V^{*},$$ where $V^{*}$ denotes the conjugate transpose of $V$. The matrix $\\Sigma$ is diagonal, the elements corresponding to the singular values of $\\hat{M}$. The solution to the system of equations is then the column of $V$ corresponding to the zero singular value. Since we have the same number of equations as we have unknowns we expect only one non-trivial solution. If none of the singular values are zero ($\\hat{M}$ is non-singular) then there is no non-trivial solution.\n",
    "\n",
    "We perform the SVD using the built-in `numpy.linalg.svd` routine which returns the matrices $U$ and $V^*$ along with an array of singular values $\\Sigma$. Before returning the solution we normalise the sum of the populations to be 1. The function returns an array of values for $\\hat{\\rho}_{\\rm{vect}}$, the order of which is the same as the elements of the output of the `create_rho_list` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def SVD(coeff_matrix):\n",
    "    levels = int(numpy.sqrt(coeff_matrix.shape[0]))\n",
    "    u,sig,v = numpy.linalg.svd(coeff_matrix) # perform the singular value decomposition\n",
    "    abs_sig = numpy.abs(sig)\n",
    "    minval = numpy.min(abs_sig)\n",
    "    if minval>1e-12:\n",
    "        # if there is no zero singular value (within floating point tolerance)\n",
    "        # return an empty array\n",
    "        print('ERROR - Matrix is non-singular')\n",
    "        return numpy.zeros((levels**2))\n",
    "    index = abs_sig.argmin() # find the position of the zero solution\n",
    "    rho = numpy.conjugate(v[index,:]) # extract the solution, which is the column \n",
    "    # of V at the position of the zero singular value\n",
    "    pops = numpy.zeros((levels))\n",
    "    for l in range(levels):\n",
    "        pops[l] = numpy.real(rho[(l*levels)+l]) # extract the terms of \\rho \n",
    "        # corresponding to populations (e.g. rho_11, rho_22 etc.)\n",
    "    t = 1/(numpy.sum(pops))\n",
    "    rho_norm = rho*t.item() #normalise such that the populations sum to 1\n",
    "    return rho_norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.92430159+0.00000000e+00j 0.09827196+1.24080432e-01j\n",
      " 0.0090494 +3.44112962e-02j 0.09827196-1.24080432e-01j\n",
      " 0.06204022+1.73894604e-17j 0.01416663+4.55273048e-03j\n",
      " 0.0090494 -3.44112962e-02j 0.01416663-4.55273048e-03j\n",
      " 0.01365819+7.47203377e-18j]\n",
      "Type: <class 'numpy.ndarray'>\n",
      "Shape: (9,)\n"
     ]
    }
   ],
   "source": [
    "result = SVD(coeffs)\n",
    "print(result)\n",
    "print('Type: {}'.format(type(result)))\n",
    "print('Shape: {}'.format(result.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For ease of calculation we combine all of these functions into one function that takes lists of the different parameters ($\\Omega$, $\\Delta$, $\\Gamma$ etc.) and outputs the steady state solution for the density matrix in the order that `rho_list` is created. Note that it outputs the solution of the entire density matrix, so we need to extract the element of interest (usually $\\rho_{21}$). All we need to do now is pass the system parameters to this function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def steady_state_soln(Omegas, Deltas, Gammas, gammas = []):\n",
    "    L_atom = L_decay(Gammas)\n",
    "    if len(gammas) != 0: \n",
    "        L_laser = L_dephasing(gammas)\n",
    "        L_tot = L_atom + L_laser\n",
    "    else:\n",
    "        L_tot = L_atom\n",
    "    H = Hamiltonian(Omegas, Deltas)\n",
    "    Master = Master_eqn(H, L_tot)\n",
    "    rho_coeffs = OBE_matrix(Master)\n",
    "    soln = SVD(rho_coeffs)\n",
    "    return soln"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.92430159+0.00000000e+00j 0.09827196+1.24080432e-01j\n",
      " 0.0090494 +3.44112962e-02j 0.09827196-1.24080432e-01j\n",
      " 0.06204022+1.73894604e-17j 0.01416663+4.55273048e-03j\n",
      " 0.0090494 -3.44112962e-02j 0.01416663-4.55273048e-03j\n",
      " 0.01365819+7.47203377e-18j]\n",
      "Type: <class 'numpy.ndarray'>\n",
      "Shape: (9,)\n"
     ]
    }
   ],
   "source": [
    "also_result = steady_state_soln(Omegas, Deltas, Gammas, gammas)\n",
    "print(also_result)\n",
    "print('Type: {}'.format(type(also_result)))\n",
    "print('Shape: {}'.format(also_result.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Measurable Quantities\n",
    "\n",
    "We now have the values of the steady state density matrix elements, from which we need to extract useful measurable quantities to compare to experiment. Often the experimental quantity of interest is the probe transmission, where the probe is the laser coupling the lowest two levels. The absorption coefficient of the system $\\alpha$ can be found from the imaginary part of the complex susceptibility $\\chi$ through $\\alpha = k\\textrm{Im}{[\\chi]}$. The susceptibility of the medium to the probe is related to the steady state coherence term $\\rho_{21}$ through $$\\chi = -\\frac{2N|d_{12}|^2}{\\hbar \\epsilon_0 \\Omega_{p}}\\rho_{21}$$ where $N$ is the atomic number density and $d_{12}$ is the dipole matrix element of the $|1\\rangle \\rightarrow |2\\rangle$ transition. From these expressions it can be seen that the absorption of the probe $\\alpha_{12} \\propto -\\textrm{Im}[\\rho_{21}]$ and hence the probe transmission can be approximated as $1 + \\textrm{Im}[\\rho_{21}]$.\n",
    "\n",
    "Since $\\rho_{21} = \\rho_{12}^*$, we will usually extract $\\rho_{12}$ from the vector of the solution as this is always at the same index (`[1]`) whereas the position of $\\rho_{21}$ will vary depending on the number of levels. The probe absorption is then $\\alpha_{12} \\propto \\textrm{Im}[\\rho_{12}]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytic Solutions\n",
    "\n",
    "Obviously because every parameter step requires solving an eigenvalue problem, this method is NOT fast. However it does allow you to compute results for parameters outside of the weak probe regime (defined as $\\Omega_{12}^2/\\Gamma_3(\\Gamma_2 + \\Gamma_3)\\ll 1$). If your system is in the weak probe regime, there are explicit analytic solutions, especially for low numbers of levels (usually <5).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fast_3_level(Omegas, Deltas, Gammas, gammas = []):\n",
    "    Delta_12, Delta_23 = Deltas[:]\n",
    "    Omega_12, Omega_23 = Omegas[:]\n",
    "    Gamma_2, Gamma_3 = Gammas[:]\n",
    "    if len(gammas) != 0:\n",
    "        gamma_12, gamma_23 = gammas[:]\n",
    "    else:\n",
    "        gamma_12, gamma_23 = 0, 0\n",
    "    bracket_1 = -1j*(Delta_12 + Delta_23) + (Gamma_3/2) + gamma_12 + gamma_23\n",
    "    bracket_2 = -1j*Delta_12 + (Gamma_2/2) + gamma_12 + (Omega_23**2)/(4*bracket_1)\n",
    "    return -1j*Omega_12/(2*bracket_2)\n",
    "\n",
    "def fast_4_level(Omegas, Deltas, Gammas, gammas = []):\n",
    "    Omega_12, Omega_23, Omega_34 = Omegas[:]\n",
    "    Delta_12, Delta_23, Delta_34 = Deltas[:]\n",
    "    Gamma_2, Gamma_3, Gamma_4 = Gammas[:]\n",
    "    if len(gammas) != 0:\n",
    "        gamma_12, gamma_23, gamma_34 = gammas[:]\n",
    "    else:\n",
    "        gamma_12, gamma_23, gamma_34 = 0,0,0\n",
    "    bracket_1 = 1j*(Delta_12 + Delta_23 + Delta_34) - gamma_12 - gamma_23 - gamma_34 - (Gamma_4/2)\n",
    "    bracket_2 = 1j*(Delta_12 + Delta_23) - (Gamma_3/2) - gamma_12 - gamma_23 + (Omega_34**2)/(4*bracket_1)\n",
    "    bracket_3 = 1j*Delta_12 - (Gamma_2/2) - gamma_12 + (Omega_23**2)/(4*bracket_2)\n",
    "    return (1j*Omega_12)/(2*bracket_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analytic solution for an arbitrary number of levels\n",
    "\n",
    "It is possible to formulate an expression for the analytic solution in the weak probe regime for a ladder scheme with any number of levels. This has to be done iteratively. \n",
    "\n",
    "For an system with $n$ levels (hence $n-1$ driving fields) we have an expression for the probe coherence which needs to be evaluated for all $m=1, 2 \\dots n-2, n-1$\n",
    "\n",
    "$$\\rho_{21} = \\frac{2iK_1}{\\Omega_{12}} $$\n",
    "where\n",
    "$$K_m = \\left(\\frac{\\Omega_{m, m+1}}{2}\\right)^2\\left[Z_m + K_{m+1}\\right]^{-1}$$\n",
    "and $$Z_m = \\sum_{k=1}^m \\left(i\\Delta_{k, k+1} - \\gamma_{k, k+1}\\right) -\\frac{\\Gamma_{m+1}}{2}.$$\n",
    "\n",
    "For example, if we have a 3-level scheme with two driving fields 1 and 2, then the expression will look like this\n",
    "\n",
    "$$\\rho_{21} = \\frac{2i}{\\Omega_{12}}\\left(\\frac{\\Omega_{12}}{2}\\right)^2\\left[Z_1 + \\left(\\frac{\\Omega_{23}}{2}\\right)^2\\left[Z_2\\right]^{-1}\\right]^{-1}$$\n",
    "where\n",
    "$$ Z_1 = i\\Delta_{12} - \\gamma_{12} - \\Gamma_2/2$$\n",
    "and\n",
    "$$ Z_2 = i(\\Delta_{12} + \\Delta_{23}) - (\\gamma_{12} + \\gamma_{23}) - \\Gamma_3/2$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This has been encapsulated into the function `fast_n_level` which takes the same sets of parameters as the `steady_state_soln` matrix method and returns the value of $\\rho_{21}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def term_n(n, Deltas, Gammas, gammas = []):\n",
    "    if len(gammas) == 0:\n",
    "        gammas = numpy.zeros((len(Deltas)))\n",
    "    # n>0\n",
    "    return 1j*(numpy.sum(Deltas[:n+1])) - (Gammas[n]/2) - numpy.sum(gammas[:n+1])\n",
    "\n",
    "def fast_n_level(Omegas, Deltas, Gammas, gammas = []):\n",
    "    n_terms = len(Omegas)\n",
    "    term_counter = numpy.arange(n_terms)[::-1]\n",
    "    func = 0\n",
    "    for n in term_counter:\n",
    "        prefact = (Omegas[n]/2)**2\n",
    "        term = term_n(n, Deltas, Gammas, gammas)\n",
    "        func = prefact/(term+func)\n",
    "    return (2j/Omegas[0])*func"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.008606577250114471-0.006320645282155274j)\n",
      "Type: <class 'complex'>\n",
      "(0.008606577250114473-0.0063206452821552754j)\n",
      "Type: <class 'numpy.complex128'>\n"
     ]
    }
   ],
   "source": [
    "Omegas = [0.1,5]\n",
    "Deltas = [5,0]\n",
    "Gammas = [5,1]\n",
    "gammas = [0.1,0.1]\n",
    "\n",
    "explicit = fast_3_level(Omegas, Deltas, Gammas, gammas)\n",
    "iterative = fast_n_level(Omegas, Deltas, Gammas, gammas)\n",
    "\n",
    "print(explicit)\n",
    "print('Type: {}'.format(type(explicit)))\n",
    "print(iterative)\n",
    "print('Type: {}'.format(type(iterative)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thermal Ensembles - Doppler Averaging\n",
    "\n",
    "Atoms in a thermal ensemble see frequency-shifted fields due to the Doppler effect. For now, we will only consider the 1D case, where all beams are either co- or counter-propagating. An atom moving in the same direction as the propagation of a laser field will see a red-shifted frequency, whereas an atom moving in the opposite direction would see a blue-shifted frequency. This will modify the detunings in the Hamiltonian as $\\Delta_{\\rm{eff}} = \\Delta - \\vec{k}\\cdot \\vec{v}$ where $\\vec{k}$ is the wavevector of the laser light and $\\vec{v}$ is the atomic velocity. Hence for an atom moving in the propagation direction we have $\\Delta_{\\rm{eff}} = \\Delta- kv$ and for an atom moving opposite to the direction of propagation we get $\\Delta_{\\rm{eff}} = \\Delta + kv$. The actual value of the sign does not matter (as the velocity distribution will be symmetric about 0), as long as we are consistent with the signs of the wavevectors relative to each other.\n",
    "\n",
    "### Distribution of Velocities\n",
    "\n",
    "The atoms will have a distribution of velocities as given by the Maxwell-Boltzmann distribution\n",
    "\n",
    "$$f_{v}(v_i) = \\sqrt\\frac{m}{2\\pi k_B T} \\exp \\left(-\\frac{mv_{i}^2}{2k_BT}\\right) $$\n",
    "\n",
    "where $m$ is the mass of an atom, $k_B$ is the Boltzmann contant and $T$ is the temperature of the vapour in Kelvin. We can either define this function explicitly, or recognise that this is a Gaussian distribution centred on 0 with width $\\sqrt{\\frac{k_BT}{m}}$ and use the built-in `scipy.stats.norm` functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "def MBdist(velocities, temp = 20, mass = 132.90545):\n",
    "    m = constants.m_u*mass #kg\n",
    "    kb = constants.k #m**2 kg s**-2 K**-1\n",
    "    temp_k = temp + 273.15 #K\n",
    "    sigma = numpy.sqrt((kb*temp_k)/m)\n",
    "    return stats.norm(0, sigma).pdf(velocities)\n",
    "\n",
    "def MBdist2(velocities, temp = 20, mass = 132.90545):\n",
    "    kb = 1.38e-23 #m**2 kg s**-2 K**-1\n",
    "    m = mass*1.66e-27 #kg\n",
    "    T = temp+273.15 #Kelvin\n",
    "    f = numpy.sqrt(m/(2*numpy.pi*kb*T))*numpy.exp((-m*velocities**2)/(2*kb*T))\n",
    "    return f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.22904507e-06 4.76995444e-05 3.59413448e-04 1.38138388e-03\n",
      " 2.70816060e-03 2.70816060e-03 1.38138388e-03 3.59413448e-04\n",
      " 4.76995444e-05 3.22904507e-06]\n",
      "Type: <class 'numpy.ndarray'>\n",
      "Shape: (10,)\n"
     ]
    }
   ],
   "source": [
    "vs = numpy.linspace(-500,500,10)\n",
    "\n",
    "MB1 = MBdist(vs)\n",
    "\n",
    "print(MB1)\n",
    "print('Type: {}'.format(type(MB1)))\n",
    "print('Shape: {}'.format(MB1.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Averaging over velocity classes\n",
    "\n",
    "Once we have calculated the response of the vapour (usually the probe coherence) for each velocity class, we need to integrate. We can either do the integration manually (using the rectangle rule for example) or using a built-in routine (such as `scipy.integrate.simpson`). In the case of manual integration, for each value of detuning we multiply the probe coherence for each velocity class by its relative abundance as given by the Maxwell-Boltzmann distribution. This abundance can be calculated using the `MBdist` (or `MBdist2`) functions we defined above. Then we multiply by the width of the velocity class, and sum over all velocities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Time-Dependent Solutions\n",
    "\n",
    "If we are not considering the steady state, we need to consider the time evolution of either the state vector $|\\psi\\rangle$ or the density matrix components $\\rho_{\\rm{vect}}$.\n",
    "\n",
    "First let us consider the state vector for a two-level system $|\\psi(t)\\rangle = c_1(t)|1\\rangle + c_2(t)|2\\rangle$ which can be written as\n",
    "$$|\\psi(t)\\rangle = \\begin{pmatrix} c_1(t) \\\\ c_2(t) \\end{pmatrix}.$$\n",
    "\n",
    "The time evolution of the state vector is of course governed by the Schrödinger equation\n",
    "$$i\\hbar\\frac{\\partial |\\psi\\rangle}{\\partial t} = \\hat{H}|\\psi\\rangle.$$\n",
    "\n",
    "The state at a time $t$ is then given by\n",
    "$$|\\psi(t)\\rangle = \\exp\\left(-i\\hat{H}t\\right) |\\psi(0)\\rangle$$\n",
    "where $\\hat{H}$ is the time-independent Hamiltonian describing the evolution of the state vector, and $|\\psi(0)\\rangle$ is the state vector at time $t=0$. (We have set $\\hbar=1$ for simplicity).\n",
    "\n",
    "For a two-level system, the interaction Hamiltonian can be written\n",
    "$$ \\hat{H} = \\frac{\\hbar}{2}\\begin{pmatrix} 0 & \\Omega \\\\ \\Omega & -2\\Delta \\end{pmatrix}.$$\n",
    "\n",
    "We can use the `Hamiltonian` function to construct the Hamiltonian for an n-level system by supplying lists of the Rabi frequencies and detunings.\n",
    "\n",
    "We can then define a function for the time-dependent part of the equation. In the above example the term `operator` is equal to $-i\\hat{H}$ but we will see later that this is not always the case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_evolve(operator, t, p_0):\n",
    "    exponent = operator*t\n",
    "    term = linalg.expm(exponent) #linalg.expm does matrix exponentiation\n",
    "    return numpy.matmul(term, p_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Often, at $t=0$ all population is in the ground state, so \n",
    "$$|\\psi(t)\\rangle = \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.+0.j 0.-1.j]\n",
      "Type: <class 'numpy.ndarray'>\n",
      "Shape: (2,)\n"
     ]
    }
   ],
   "source": [
    "Omegas = [1] # 2pi MHz\n",
    "Deltas = [0] # 2pi MHz\n",
    "H =  Hamiltonian(Omegas, Deltas) #Create Hamiltonian\n",
    "psi_0 = numpy.transpose(numpy.asarray([1,0])) #state vector at t=0\n",
    "t_pi = Omegas[0]*numpy.pi\n",
    "\n",
    "sol_pi = time_evolve(-1j*H, t_pi, psi_0)\n",
    "\n",
    "print(sol_pi)\n",
    "print('Type: {}'.format(type(sol_pi)))\n",
    "print('Shape: {}'.format(sol_pi.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considering only the state vector, it is not possible to consider the effects of decay. To do this we need to consider the density matrix $\\rho$ and its time evolution. \n",
    "\n",
    "The time evolution of the density matrix is governed by the Lindblad Master equation which can be written\n",
    "\n",
    "$$\\frac{\\partial \\rho}{\\partial t} = -\\frac{i}{\\hbar}\\left[\\hat{H}, \\rho\\right] + \\mathcal{L}(\\rho)$$ where $\\hat{H}$ is the total Hamiltonian in the interaction picture and $\\mathcal{L}(\\rho)$ is an operator describing decay and dephasing within the system.\n",
    "We can rewrite this in terms of a vector $\\rho_{\\rm{vect}}$ that contains all the elements of the density matrix, and a matrix of time-independent coefficients $\\hat{M}$ found from the Master equation as\n",
    "$$\\frac{\\partial \\rho_{\\rm{vect}}}{\\partial t} = \\hat{M} \\rho_{\\rm{vect}}.$$\n",
    "\n",
    "We can use some of the previously defined functions to construct the operator $\\hat{M}$.\n",
    "\n",
    "Then we can find the values of $\\rho_{\\rm{vect}}$ at any time $t$ through\n",
    "$$\\rho_{\\rm{vect}}(t) = \\exp\\left(\\hat{M}t\\right)\\rho_{\\rm{vect}}(0)$$\n",
    "where $\\rho_{\\rm{vect}}(0)$ is the value of the density matrix in vector form at $t=0$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_dep_matrix(Omegas, Deltas, Gammas, gammas = []):\n",
    "    # first create decay/dephasing operators\n",
    "    L_atom = L_decay(Gammas)\n",
    "    if len(gammas) != 0: \n",
    "        L_laser = L_dephasing(gammas)\n",
    "        L_tot = L_atom + L_laser\n",
    "    else:\n",
    "        L_tot = L_atom\n",
    "    H = Hamiltonian(Omegas, Deltas) #create the total Hamiltonian\n",
    "    Master = Master_eqn(H, L_tot) #put Hamiltonian and decay/dephasing operators into Master equation\n",
    "    rho_coeffs = OBE_matrix(Master) #create matrix of coefficients\n",
    "    return rho_coeffs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0. +0.j    0. +0.05j  0. +0.j    0. -0.05j  5. +0.j    0. +0.j\n",
      "   0. +0.j    0. +0.j    0. +0.j  ]\n",
      " [ 0. +0.05j -2.6-5.j    0. +2.5j   0. +0.j    0. -0.05j  0. +0.j\n",
      "   0. +0.j    0. +0.j    0. +0.j  ]\n",
      " [ 0. +0.j    0. +2.5j  -0.7-5.j    0. +0.j    0. +0.j    0. -0.05j\n",
      "   0. +0.j    0. +0.j    0. +0.j  ]\n",
      " [ 0. -0.05j  0. +0.j    0. +0.j   -2.6+5.j    0. +0.05j  0. +0.j\n",
      "   0. -2.5j   0. +0.j    0. +0.j  ]\n",
      " [ 0. +0.j    0. -0.05j  0. +0.j    0. +0.05j -5. +0.j    0. +2.5j\n",
      "   0. +0.j    0. -2.5j   1. +0.j  ]\n",
      " [ 0. +0.j    0. +0.j    0. -0.05j  0. +0.j    0. +2.5j  -3.1+0.j\n",
      "   0. +0.j    0. +0.j    0. -2.5j ]\n",
      " [ 0. +0.j    0. +0.j    0. +0.j    0. -2.5j   0. +0.j    0. +0.j\n",
      "  -0.7+5.j    0. +0.05j  0. +0.j  ]\n",
      " [ 0. +0.j    0. +0.j    0. +0.j    0. +0.j    0. -2.5j   0. +0.j\n",
      "   0. +0.05j -3.1+0.j    0. +2.5j ]\n",
      " [ 0. +0.j    0. +0.j    0. +0.j    0. +0.j    0. +0.j    0. -2.5j\n",
      "   0. +0.j    0. +2.5j  -1. +0.j  ]]\n",
      "Type: <class 'numpy.ndarray'>\n",
      "Shape: (9, 9)\n"
     ]
    }
   ],
   "source": [
    "Omegas = [0.1,5]\n",
    "Deltas = [5,0]\n",
    "Gammas = [5,1]\n",
    "gammas = [0.1,0.1]\n",
    "M = time_dep_matrix(Omegas, Deltas, Gammas, gammas)\n",
    "print(M)\n",
    "print('Type: {}'.format(type(M)))\n",
    "print('Shape: {}'.format(M.shape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}