PhysRevA-91-013834.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import *\n",
    "import numpy as np\n",
    "from qutip import *\n",
    "from scipy import linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#constants\n",
    "hb=6.63*1e-34/(2*np.pi)#Planck constant\n",
    "wl=2*np.pi*500*1e12 #frequency of driven laser\n",
    "Nc =4                   # Number of cavity states\n",
    "Nm1 =10                # Number of mech2 states\n",
    "wm1=2*np.pi*10**6        #mech1 frequency\n",
    "kappa =0.1             # Cavity damping rate   \n",
    "gamma1 =1e-6             # Mech1 damping rate\n",
    "n_th1 = 0    \n",
    "delta=1\n",
    "eta=1e-4\n",
    "g0=1e-4 #optomechanical coupling rate --cavity with mech1\n",
    "plist =np.linspace(0.0001,10,1500) #range of normalized detuning delta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a= tensor(destroy(Nc), qeye(Nm1))\n",
    "b1 = tensor(qeye(Nc), destroy(Nm1))\n",
    "num_a = a.dag()*a\n",
    "num_b1 = b1.dag()*b1\n",
    "    \n",
    "#quadratures\n",
    "x1 = (a + a.dag()) / np.sqrt(2)\n",
    "y1 = -1j * (a - a.dag()) / np.sqrt(2)\n",
    "q1 = (b1 + b1.dag()) / np.sqrt(2)\n",
    "p1 = -1j * (b1 - b1.dag()) / np.sqrt(2)\n",
    "R=[x1,y1,q1,p1]\n",
    "cc = np.sqrt(kappa)*a #cavity collapse\n",
    "cm1 = np.sqrt(gamma1*(1.0 + n_th1))*b1#mech1 \n",
    "cp1 = np.sqrt(gamma1*n_th1)*b1.dag()#mech1\n",
    "c_ops = [cc,cm1,cp1]\n",
    "En=[]#Entanglement\n",
    "var=[]#position variance for duffing oscillator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# definition of entanglement with one variable 'delta'\n",
    "for P in plist:\n",
    "    E=np.sqrt(2*0.1*wm1*0.001*P/(hb*wl))/wm1\n",
    "    H =delta*(num_a) + num_b1 - g0*(b1.dag()+b1)*num_a  + E*(a.dag()+a)+eta*(b1+b1.dag())**4\n",
    "    rho_ss= steadystate(H, c_ops)\n",
    "#covariance matrix\n",
    "    varxc=variance(x1,rho_ss)\n",
    "    varpc=variance(y1,rho_ss)\n",
    "    varq1=variance(q1,rho_ss)\n",
    "    varp1=variance(p1,rho_ss)\n",
    "    covxcpc=1/2* expect(x1*y1+y1*x1,rho_ss)-expect(x1,rho_ss)*expect(y1,rho_ss)\n",
    "    covxcq1=1/2* expect(x1*q1+q1*x1,rho_ss)-expect(x1,rho_ss)*expect(q1,rho_ss)\n",
    "    covxcp1=1/2* expect(x1*p1+x1*p1,rho_ss)-expect(x1,rho_ss)*expect(p1,rho_ss)\n",
    "    covpcp1=1/2* expect(p1*y1+y1*p1,rho_ss)-expect(p1,rho_ss)*expect(y1,rho_ss)\n",
    "    covpcq1=1/2* expect(q1*y1+y1*q1,rho_ss)-expect(q1,rho_ss)*expect(y1,rho_ss)\n",
    "    covq1p1=1/2* expect(q1*p1+p1*q1,rho_ss)-expect(q1,rho_ss)*expect(p1,rho_ss)\n",
    "    V=array([[varxc,covxcpc,covxcq1,covxcp1],[covxcpc,varpc,covpcq1,covpcp1],[covxcq1,covpcq1,varq1,covq1p1],[covxcp1,covpcp1,covq1p1,varp1]]).astype('float64')#  change the datatype\n",
    "    var.append(varq1)\n",
    "    def logarithmic_negativity(V):\n",
    "        A = array([[varxc,covxcpc],[covxcpc,varpc]])\n",
    "        B = array([[varq1,covq1p1],[covq1p1,varp1]])\n",
    "        C = array([[covxcq1,covxcp1],[covpcq1,covpcp1]])\n",
    "        sigma = linalg.det(A) +linalg.det(B) - 2*linalg.det(C)\n",
    "        eta = sigma/2 -np.sqrt(sigma**2 - 4 * linalg.det(V))/ 2\n",
    "        if eta < 0.0:\n",
    "            return 0.0\n",
    "        nu = np.sqrt(eta)\n",
    "        lognu = -np.log(2 * nu)\n",
    "        logneg = max(0, lognu)\n",
    "        return logneg\n",
    "    En.append(logarithmic_negativity(V))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAERCAYAAABy/XBZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xu4XeOdB/DvL6UXVdqa6BBDp+20DBpVpCi2qbrEM8UT\nZiRaWsY02hKU0WoqyXTUpVWEdCotqaqU1n1c4pDYSEQTBIlE0CAJEQm5CBEn5/zmj/dsZ5191t7r\n/r7v3uv7eZ7z7MvZa63fXnut9VvvZb1LVBVERFQ+A1wHQEREbjABEBGVFBMAEVFJMQEQEZUUEwAR\nUUkxARARldQmthcoIi8BWAOgG0Cnqu5lOwYiInKQAGAO/BVVXeVg2URE1MNFFZA4Wi4REQW4OBAr\ngPtEZLaInOxg+UREBDdVQPuq6jIRGQiTCBao6nQHcRARlZr1BKCqy3oeV4jIrQD2AvB+AhARDk5E\nRJSCqkqSz1utAhKRzURk857nHwVwMIB59Z9TVf6pYsyYMc5j8OWP64Lrguui+V8atksAnwJwa89Z\n/iYArlfVDssxEBERLCcAVX0RwG42l0lEROHYHdNjlUrFdQje4LroxXXRi+siG0lbd1QUEVHfYiIi\n8p2IQH1uBCYiIn8wARARlRQTABFRSTEBEBGVFBMAEVFJMQEQEZUUEwARUZSXXgJGjADm9Ru5pqXx\nOgAioih77QXMng1suSWwerXraELxOgAioiK8+KJ5XLPGbRw5YwIgIiopJgAiopJiAiAiKikmACKi\nkmICICIqKSYAIqKSYgIgIiopJgAiopJiAiAiKikmACKikmICICIqKSYAIqKSYgIgIiopJgAiopJi\nAiAiKikmACKikmICICIqKSYAIqKSYgIgIiopJgAiopJiAiAiKikmACKikmICICIqKSYAIqKSYgIg\nIiop6wlARAaIyBMicoftZRMRUS8XJYBRAOY7WC4REQVYTQAish2AoQB+Z3O5RETUn+0SwKUAzgag\nlpdLRER1NrG1IBE5HMByVX1SRCoApNFnx44d+/7zSqWCSqVSdHhERC2lWq2iWq1mmoeo2jkZF5Gf\nA/gmgI0APgLgYwBuUdXj6z6ntmIiIopl4EBg5Urz3NPjk4hAVRueWIdO4+JgKyIHAPihqn4j5H9M\nAETklzZNALwOgIiopJyUAJphCYCIvMMSABERtRMmACKikmICICIqKSYAIqKSYgIgIiopJgAiopJi\nAiAiKikmACKikip3AnjrLeDOO4HOTteREBFZV+4EcNRRwL/+KzBunOtIiIisK/dQENJz1fSOOwIL\nFthZJhG1Hg4FQURE7YQJgIiopJgAiIhKigmAiKikmACIiEqKCYCIqKSYAABvu3URERWJCYCIqKSY\nAIDeC8KIiEqECQBgFRARlRITABFRSTEBAKwCIqJSYgIAWAVERKXEBEBEVFJMAACrgIiolJgAiIhK\nigmAiKikmACIiEqKCYCIqKSYAIiISooJAOB1AEQuvPUWcPrpwOOPu46ktJgAiMiNsWOByy8H9tjD\ndSSlxQQA8DoAIhcWLXIdQekxAQCsAiKiUtrE5sJE5EMAHgLwwZ5l36Sq42zGQEREhtUEoKobRORA\nVX1HRD4AYIaI3KOqs2zG0Q+rgIiohKxXAanqOz1PPwSTgNzXv7AKiIhKyHoCEJEBIjIHwGsA7lPV\n2bZjII/dcANw1llMykQWWK0CAgBV7QbwJRHZAsBtIvLPqjrfdhx9sArIH8OHm8cjjgD2289tLERt\nznoCqFHVtSLyAIBDAfRJAGPHjn3/eaVSQaVSKTqYYudPya1d6zqC1rdgATBoELDFFq4joQJUq1VU\nq9VM8xC1ePATkb8D0Kmqa0TkIwDuBXChqt4d+Ixai6l25v+FLwDPPmtnmdRc7Te5807g8MPdxtLK\nnn4aGDwYGDgQeP1119GEO+oo4LbbzHPfT8IGDgRWrjTPPY1VRKCqiaozbJcAtgFwrYgMgGl/uDF4\n8CeinDz4oHlcscJtHOQ1291A5wLY3eYyiUrJ07NU8guvBCYiKikmAKJ2xJ5tFAMTABFRSTEBkJ9Y\nh01UOCYAIqKSYgIgIiopJgCADWZEVEpMAADrm4molJgAiNoRT2ooBiYAgFVAFG39eqCry3UU7W3+\nfOCaa5i8LHI2GqhXuMFRM2+/DWy+ObDLLsDcua6jaV8772wet9rKDAdOhWMJgNpfdzcwaxawYUO6\n6efN6/tIxeLIvNYwAVD7u+wyYMiQ3pvNEBGAFAlARD5URCBEhZk0yTzeeqvbOABTGiHyRJoSwBgR\n2S73SIja3c9+Bmy6qZ0qDnZsoBjSJIDfAegSka+IyLdE5KS8gyJqS+edZ0oAF1xQ/LLYsYFiSJMA\nFgPYD8CZAHYC0JFrREREZEWabqCLAPwWwAmquj7neIiIyJI0CeBQAJsB+KqIdALYU1V/kW9YRERU\ntMgEICJXA7gHwExVfQXAKwA2qOq7Pf9nKYCIqGbNGuCWW4Bhw4AttnAdTVNxSgCvqepNgdfvADhA\nRHYEsFhV7ygmNCKiFjRiBHD33ebvL39xHU1TcRLAiwAgIofDNPrOAvCAqt4vIlMBMAEQEdXcfbd5\nvMP/Q2PsXkCqeheAvwewEECtk7EHV9YQtRj20SdPxCkBXCAiBwKYDmA1gBWqWruc8Z3CIiPyBQ/Y\n1KbiJICfwFT7DAHwjwD+KiIbATwFYCsA1xQXHhERFSWyCkhVJ6rqk6p6laqepKp7Avg6gD8D+GDh\nERIR+aQF6vbjSnU/AFVdB2CaiLyVczxERH474gjg5ZeB7bd3HUlmmYaDVtXZeQVC1AfHsiGfLV/u\nOoJc8H4ARFQec+YAQ4fypjM9eEtIIiqP/fcH1q0D/vY3YOFC19E4528JYP584OabXUdBRO1k3Trz\n+NprbuPwhL8lgNoNoh991NzOj6hd8LoC8oS/JYCa555zHQFR8R55pHxnpWzod87/BEDtadUqc4vE\nxYtdR+Le7NnAvvsC22zjOhJ3Dj7YdQT+u/12YMyYXBOnv1VA1N5OOQW48Ubg9783DXJlNmuW6wjc\nCFaF3XefuzhaxZFHmseDDgL22y+XWbIEQG48+qh5XLTIbRzUWhYuBHbcsdwdRN58M7dZWU0AIrKd\niEwTkWdEZK6InGZz+UReKqIuvF3r108+2SSBo492HYk73d3Rn4nJdhXQRgBnquqTIrI5gMdFpENV\neVVG2bAnDKXR2ek6AvdyTABWSwCq+pqqPtnzfB2ABQAG2YyBPNGuZ6i+YIK1Z/Ro4Kyz7C2vqyu3\nWTlrAxCRTwPYDcBfXcVAFEu7HEwvuwzYdVfTA4vyc/75wCWXAO++a2d5rVoCqOmp/rkJwKiekgAR\nFe2MM4B584Dx411H0p5slWpbuA0AIrIJzMH/OlW9PewzYz/3ufefVxYsQMVOaERu2K4Oy/EAYt2G\nDa4jcK/n96tWq6hWq5lm5eI6gGsAzFfVyxt9YGywX/hOO1kIiahEWrn9Zc6c9NP+8Y/5xeFSz+9X\nqVRQqVTef3vcuHGJZ2W7G+i+AI4D8C8iMkdEnhCRQ23GQORcfZtCEW0MrXyQL8rPfuY6gnzk2Ahs\ntQSgqjMAfMDmMomsUgVGjQL23hsYPjz+NERxtXojMFHbHvSmTgWuuAIYMcJ1JI2167ovCyYAIk+t\nWeM6Amp3TABERB4IlqZslaxyXA4TAMAisWuTJwOTJrmOojy4vefHRQJo1UZgolDHHWceg42m7XaQ\n2rjR7vLa5epl6o9VQDnjzuKHHM9sGvrDH4Annih+OUHd3cA//IPdZTZLoO2WXF1yUQJo5SuBvcQd\nojxOOME82vzN3367+e0euf21rhZPACwBEBWNJUzKExuBc8Yd1L4ynfVGbV/c/loXSwBElEmZkiFl\nxwRA1EJsjP1DbrR4N1AmAJs6O81QAe+84zoSsqnoA8OUKeZGL3Pn+hFPmTRLAEUlepYAWtS4ccBB\nB/X2e6e+2vXAVPT3Ouwwc6MXn8cfsuXCC4Evf9n0vHKtqN89x8TCBGDTddeZx9tucxuHr8J2GFXT\nb583AonGdQT8+Mdme7E19j+HgrDo+utNdl+2zHUkZMtvf2t+82HDXEeSXpId9t570y8n7plhMJ6O\nDuCCC9qv9GXrrmcuEkCOWisBfPObJrv/9KeuI6GsGp3t17v2WvN4113FxlOkJAeGxx4rLo4whxwC\nnHsu8OCDdpdLXmitBFDDom57asEzqFiivlde3ztNCaCm2ZXK1BhLAG2gBX8461SB9evzm1/YwcrX\n3yEYq68xAtkaB22Mw2STz7+TR5gAKJ5jjwU22wx48cV85tdKO2jWWG191ywJwFadebvJUgJQdb7e\nmQAAXpgTx5//bB4nTy5uGa2QFNLE6Nv3Coun3RKArX06SwIYOhTYYQfgvffyjSmB1kwAqqbfc17F\nVt920LIK+x18SM55bx9FfacBGXbndqsCagVTpgBLlwILFyabrrTdQGuuv95c+fi977mOhMrG5xJA\nu1QBtdIJGRuBHZo4MZ/5+HCW2Y66uoAjjjD9zOPwdQdqlTaANWtMyTiNdisBuLgoy9ZQEDlq7QSQ\nl6QbS2cnMHOm/dv8tZqHHwbuuMP0M4/D1wQQlEcJYMkS4Omns80zzJIlpmT85JPJp/UpAbTAgbNd\nMAGkccYZwD77AGef7ToSv3V2Jvu8rwkg77g6OoDBg4HXX893vjUPP5x8Gp+qgPLQCo3AHmACAJJv\nLBMmmMcrr8w/lrJowZ0FQN+4X33VDPAXdSBv9F2XLMkvriTC4vGpBNBKmiWAorbx0jcC5y3LCj33\nXGC33fK9SKqVrFgBXHyxeczK16TQKK5vfAMYOzZ6dFdfv1dQu5UAKBYmgKwuuAB46ilT111GxxwD\nnHMOMHx4sula6UrgoGCMjz9uHmfOdBNLWmUoAfjQCNwCmAAom9ogYg88kGy6qMHgfNqZiuoFVHvf\nh+/KEoA7Dn9/JgCiJNJcrNZoB99zT+Dyy7PHlId2SwCt1Ag8eHB+XdoTYgIgf7g6E3rpJWDIkMbV\neFE7edoEAACnnx7vYJXnbUTLUAVUtLDSW5bt97vfNddx/PKXpnNBnGXngAmA/OEqAZx2GjBrlrlo\nLY2sZ5tR37ujA/joR02Dc1HarQTgg6Tbxfe+Z7qWH3JIMfGEYAIAgOeeA1audB0FuUoAUfePzXqW\nl/V7jR5tHseNyzafZvJKADNmALvsAjz6aD7zSyuPbemll8zd6G65pfH882y3eugh85j2au4UmABq\njjrKdQQU3IGOOcZdHHlLcmBI2jsqzY1c4t6NLY1KBXjmGeDgg/OZn0ujRpk7EDa7HamvHRdiYgKo\nmT7ddQQU3IGSXkVcpCLbAOJ8ttH0d90FbLNN/HnbUBsepR2ui2nW7tKCB/swTADkRtwDXdJ61BUr\ngIMOAm6/PV1cPmp0sBk/Pvk0ZZFHL6Bm67CIKqC4WrURWESuFpHlIvJ09KcdOPpoYNUq11GUVx71\n0KNHA1OnAkcemX1eNTZLAEVMn7fRo00Vj889h371K+Cxx7LNw7f1XgDbJYBJAOw1cSd1883AmDGu\noyivPHa4tWuzzyOpohOAbz10zj8fuO8+YM4c15E09sIL5jqLLGyXAOJOn+M1DlYTgKpOB+D3KfYb\nb7iOoLx8PeMqupgfNc+8l5llfsFpN9sseyw+i7OeWAVElBNfE0CURmdk8+aZkULbqQro3Xd7n3/4\nw+7isKFZyauI38TBfRA2sb7EGMYGnlcWLEDFURzUw/czm/nzTePvAQfkG09N2jaAXXc1jy++2Hz+\nRVch5WnNmt7nPsVVhKwlgIKHg65Wq6hWq5lm5X0CwE47OYqCAJhufYMH21lW2h1m553N49Kl+cWS\nRNYDuO2zySzLC94Fr8wJwGUvoB6VSgWVSuX91+NSXCjoogpIev7819lprmj0ubdD0Z56ypxh2xB1\nRW6UqDFU8lD0Tj5lCjByJPDee9mWmXSauJ9v8QufEsn6/Yqq0mnVRmARmQzgEQCfF5HFIvIdm8tP\nbORIYO+9ix2DpazCdq799rMfRxw2u4FOnw5cdRVw9dXppi9asF7cVu+kxYuB55+3s6wgX0sArdoI\nrKojVHVbVf2Qqm6vqpNsLj+xa64xj7/5jds4ymLduub/f+ghP0tjRdThB28zmfRAu349cP/9yeLx\nuQSwww7A5z/ftwHahjgJoMWxFxC1jgMOAC69NP/5Jqmjd7HjJ13mOeeYYSKKjsX2urB9jYev3UBz\n1BoJYOpUYOFC11FQPRfjJ02ebH+ZUVrtSuC8rgNowQNeIkVWAXmy7rzsBdTH888Dxx/vOgoK8/LL\nzf/vYiMvouGNQ0H0ctEG4EqR692T39T/EsCiRa4jsGP5cnOD+RUrXEdiR7A/eRKt1Ge+Jks30FWr\nzLbhizKVAOJcCOaiBJDjeve/BFAWRx0FzJwJVKvAvfe6jiYfjQ7Wa9fme4vDomUtAWTxyU/mP08b\nVUDB6wVaVZFtAHkdxFesAPbZBzjllFST+18CKIuZM81j7a5ArSbJQXDGDD/iyDJNkunb6Uw5eFbc\nTt8rTJG9gPJad7/6lRn47oc/TDU5EwDZl2SoZgfjo/Tjog2gyO+dZd5lqgJKUwK4+OL08467PseN\nAx5/3DzP2A7DBNCuurqAE09022um0Rg4watcy8DGlbxZ553mOgAbjcAuE06aweDOOSf98pJ8v5xu\nHM8E0K5uuw2YNAk47rhs88my0333u9mW7QuXbQC+sX1A9rXXURGNwEmGMsnpxlVMAO3KxY1RWk13\ndz4HsaVLmzd6tlNVSZ5tAM88E/0Z2yWOmjVr/GoELmikUSYA302enO7WdnkddPI6u73llt56yzz9\n6U/Rnwn7Dl1dwD/9E7D//tHT1+/k3d39z8CmTo03fZr/582HC8Fmz443vo+LKqCODuDjHweebnLn\nWteNwEwAJTBnjqnCyXpruyzy2umGDQP22COfedVs3AiMGJFu2mXLzDUmSa9mVgUOO6x/98y8xyiy\nXa1keyygZmMVBbmoArrggvifjbM+zj8f+MEPGk+XdDk1zz0Xv9G5ASYAny1enH7adq2XDu4IWQ66\nSRqi63e+jo5ky0p7oCzqjDfLgbTZAblRvM8+m/66j2ZVQMuWmYHv0shaRZqkDWD0aGDCBOCVV6I/\nG7W8oJR9/4OYAPJ0yimmdb6d6nx9TiT16zlJrGl7IjX6bZste9iwdMsq6ox30iTgiSfSTdvsgBcW\n74wZ5qZOe+5p7q1xxRVmuri/VbM2h223BT796Xjz8cH48cCDD5rneRwjchgdlVcC56k2bPQLL5j6\n5XbgczKLewYaJskNV+LMt9kB7bnn4sVUr6gEsH498OUvZ++eGqdh8s47zeP8+ebeGgCw447plhe2\nPoLDZtukarahyy/v+14zF19s/lSTrfuHHwYuuyxdnBGYAGzp6spWpUP9ZTlAbtiQbjqbCdHH5Nvs\njDzs9whLjEn2A5+vPB4/Hrjyyt7XRfUCitNRISVWAdly6ql2l+fbzlKELAkgSxtAmCKqynzq917T\n7Iw8brxJ1pWrbqBRVE0nDZdy2OaYAGz53/+N9znfDty+xROUpQooWALI48yttjPmub6KXvdpevQk\nrQJKc5B6+23Ta+aRR7LdGc02F4PBZcQqoCLY+HGjGtLyOiP1eacLOwON+73zbgPIotH8G637t94C\nLroIWLAgv+XG/Z2TNgKncfHFJkFPmAB86lP5zz/Mq6+aRtXPfCbe5z3px58VE4Bv4hzA3nvP9KrY\ne+/i71fsyYYaKqwEEDfetP3LbZYAGsU1erSpf84qTQkgjzYAkebbebB0FrwXQpH3gx40yDy+/Xa8\nz4dta0nGWMpjO2EVkKfyGm2xkenTzVWKV12Vfjl5xuNKliqg4MEkjzPLLAngjTfC3280r3nzki8j\nav4+VQE1YqM0umpVvHWxalX/q9DjJoAbbvBmv2qfBNDVBXz1q8B//IfrSIoXtiMsXmwSwvXXm9d5\nbGBvvtm+VUDBBNBox73xRjNcQZwDZZYEcNFF4e8Xve7jfK+5c4Hf/z78wqe0jcBpudwWV67s+/ra\na+NNFxbz8OHF3hMjAf+rgOLuUAsWmJU6Ywbwu98VG1OYvC6Rj3MAq9+o3nwT2GGH3tf//u/pl1/z\n0EPAAQe4v56h2frIUgKIqgKqVoFjjzXP/+//oufX6r2AGq27L37RPA4aBHz968mrgPKUdv4rVph9\npF6jdrQ4v2Xcs/1GMae9NiSIVUABrotUSeuU/+u/zAVjadR/17/9re/rPOpKa3XMcQbsciVLG0BU\nFdCSJX3nG/Y8TJrtMNjQGWdeeW3rSb5XrcE5jyqgtAeuRvvVM8+YhuNGXXu33tpcfPbaa33fz7Ie\nsyYAT0rW7VMCcC1pAvjFL0wxMs0Nv6Pm39VlZwgH179NWBVEnNFBgegEkPSgHPf/YXbfHbjnnv7v\n+1AFVFNbX6ed1niavBqBo2Kot8su5nHTTZtPXz/+T1cXMGBA9raQZhr9hnGmHzWq8Kqi9ikBxNHZ\n2b8uLy9pepWEXcYeZ8OIOvPKowTQaAedODH7vPOSpQ46KgEEv/+JJ/Y+nzQpfH5ZkmGj3yuqvSGr\nNAkgOERyliq4NKJ+3/nzk89v/Xrgjjt63+vqijdwXdwSQNLfNmj8+OZDqJeiCihN/+RGvvQlYODA\nvqPyZRFcZvCHTnIAXrgQ+M53wv/305+Gvx9cJ6tXA0OG9P3/xo3F9VMOXsafZgNMGtdDDzX+X/22\nkWRkyEYJ4Pvf7///4EnDmDHh8/vRj8xjmvXuqpogSQIIu+FN2hJAs/ebiVofSdv+uruBkSP73qO6\nqwuYNSt6WhslAAvaJwHEUbsDUW1EvqwaHfSTxHzwwaaXRU1tx3j9ddMDI0xw/mFnpFEJyOXGl2TZ\ny5Y1/3/9iJZJEkDwgBZcX7/+df/34qjdtCfNug1b1rXX+pUAwmJM2wic9sw17+sAuruBP/yh/3tx\nxCkBqDIBZFb0+CJ5SZsA6gfGqsVWXwzt6AD++Me+n6l/XtPs9oRz5wLbbNN/w6/XbH1eeSVw883N\np28kybrp7Gz+/2OO6fu6PgE0G/O92e+lmq6qsLs73fgwYQe2RYvyawT+9rej5xOnXSkqjqKvA0hy\no5Y4wr5T3CQTtr7qv/977/VtM2n22TRKUQWU5geJ2piTrvyVK80NzoM9Q+oFD7p5nLnVH/wOOQT4\n1rfMWXGanbVm5EjT8HzCCeljO/VU4Oij002b5Cwu6e9UnwDqE0RQs99r4MDG1XLNnHcecMklyacr\nugqoUZ91V1VAnZ3pumrfe2/yaZoJizfu9hm3x0/tpC3OZx3wPwHEXVHN6uDXrgXGju19nSb7TpzY\n/KAXtvw0P3Jth2k0XPHatdHzbVYCCHPppY3PVPKWZJ0kXX+vvtr3dbM7dzUrATS6MjfK+eenm67R\nQcdmv/o0VUBpG+EvuSSffvBZ/c//9H8vSwIo8uSmIK2dAIJnF40SgKo5cx43ru97aTRrHBo9uvd5\nd7c5w9lsM3MXpDQa1Wcff3zfRJR1QwSAM880d2p6/nnzeOCByW/h19UFzJwZ/bk8q4DqBX/jKHkP\nBZFF0nrivKpVgicKaRJA8GYoYfOYOjX8nrVZB7HLS1hpzVYCcL3N9WifBBDcmIPPx4/v282rKG+9\n1fv8gQeAk082Z/FXXJFufo0SQH0SylKPWW/DBlMSqFaBu++O/nz9TrDPPtHTJNnwwy7seecdYNq0\n5KWcesHGc9c7Y9ISQF5nj8F1+Otf992Gwz5bH8+0aX1f1///oIPiD65WlKTrKu8qIM+1dgLo7u5t\n6GtUAjjrrP7T/eQn6eOJs0EFSwNp5x/3DDxsg417cLzqKmDXXfsvO640B+GsCWDhQuBrX0tX1x70\n1FPpYipCkgN9VM+oRsLmVX+SceGFjafv6gqvlrznHtNg3WgZriU9GYpb/Zc1AUStq9qYXgXzPwFE\n9cbYcktTFG2UALbYov80S5emj8fGCJxA/C6NYb154m70I0f2HVky6YEw7AD97rtm4z7uuPBpsiaA\nmuuuiz+fKK4TQKPvGRbXttuaEmZSYcl69eq+r3/+88YHwEYJYOhQ4LOfNc9dr8cwSU9SDjkk3ueK\nTgBnnhl/XhlYTwAicqiIPCsiz4nIOZET1PpXN3P66UCl0vt6wgRTnF29GvjYx8KnifoBGvVSSDr+\n/uTJ4Y1NjdSqteImgGef7f9e/YVgcXt7jBwZb5k1YQeENWuAs8823ztMknHsk9y2MQvXB6533w1/\nP8+4wtpTwgZICxuSAjAl0qixq1yvxzBZqwobyVo6izr+xLnZfat1AxWRAQCuBHAIgJ0BDBeRHXNf\n0HnnmbFcPvGJxivppJOaz+Pkk8PfT7PSG13RG6a2YbzxBqrJl2R0dTXuFttsI03aYB2WANata149\nE+yNFSVwIVy1/n+1i/ry0Kx7rw2NDqzDhoW+XU2zjLCD/apVzacJHqSuvtrchCju5y2pRn2gqJOI\n+gEYAeC//zv+9J4kS9slgL0APK+qL6tqJ4AbABxRyJJqZ9AvvRT+/0Zjuqxc2bwBNMnVplmcdlr6\nBDB/ft8DZLBdodH6SKO+ERAwCSAvgfsoV7PMp6izwKI1SHLVNPM6ImQ3C0sKgOk++8lPmoHS4jrs\nsIYJq0jVqA9stVUxCw6rmk4ysGMe2+RHPpJ5FrZHAx0EIHi6tRQmKfijo6N5//GFC4td/oYNzQeA\niuP44/u+njrVdNFM0k0yjrDqgt12y2/+eZzln322P90OXaofNgMIL+Xef3/8m50ETZmSfJp2k+Tk\nJ4/rIHIoRfg/HHQZ7bFHvvM79NB859dKfvlL1xG0ljQHf0ru1luzzyOHK6NFLdbbichXAIxV1UN7\nXv8IgKrqRYHPeNiXjIjIf6qaqJHSdgL4AICFAL4GYBmAWQCGqyrL6EREllmtAlLVLhH5AYAOmAbo\nq3nwJyJyw2oJgIiI/OHVlcCJLxJrUyKynYhME5FnRGSuiFgaqtNPIjJARJ4QEQuDOvlNRLYUkb+I\nyIKe7WMizjV9AAADcElEQVRI9FTtSUTOEJF5IvK0iFwvIh90HZMtInK1iCwXkacD731CRDpEZKGI\n3CsiW0bNx5sEYO0isdawEcCZqrozgL0BfL/E6wIARgFIeMPXtnU5gLtVdScAgwGUsgpVRLYFcCqA\n3VX1izDV2ce6jcqqSTDHyqAfAbhfVb8AYBqAH0fNxJsEAJsXiXlOVV9T1Sd7nq+D2ckHuY3KDRHZ\nDsBQACnuINJeRGQLAPup6iQAUNWNqtrktmdt7wMAPioimwDYDMCrEZ9vG6o6HUD9pdxHAKj1470W\nwJGI4FMCCLtIrJQHvSAR+TSA3QD81W0kzlwK4GwAbKwC/hHAShGZ1FMlNlFEsl8O2oJU9VUAlwBY\nDOAVAKtV9X63UTm3taouB8xJJICtoybwKQFQHRHZHMBNAEb1lARKRUQOB7C8pzQkPX9ltgmA3QFM\nUNXdAbwDU+wvHRH5OMwZ7w4AtgWwuYiMcBuVdyJPmnxKAK8A2D7werue90qpp1h7E4DrVPV21/E4\nsi+Ab4jIIgB/AnCgiETczb6tLQWwRFVrQ+TeBJMQyuggAItU9U1V7QJwC4AYdyRqa8tF5FMAICJ/\nDyBySFGfEsBsAJ8TkR16WvOPBVDmXh/XAJivqpdHfrJNqeq5qrq9qn4GZnuYpqrHR03XrnqK90tE\n5PM9b30N5W0cXwzgKyLyYRERmHVRtgbx+lLxHQC+3fP8BACRJ47ejAXEi8R6ici+AI4DMFdE5sAU\n5c5VVY64RacBuF5ENgWwCMB3HMfjhKrOEpGbAMwB0NnzONFtVPaIyGQAFQBbichiAGMAXAjgLyJy\nIoCXAfxb5Hx4IRgRUTn5VAVEREQWMQEQEZUUEwARUUkxARARlRQTABFRSTEBEBGVFBMAEVFJMQEQ\nEZUUEwBRDCIyWkReEJGTROT7IvKbnsH6iFoWEwBRPI8BmKKqV6vqBABr0DvuClFLYgIgimcIgCcC\nr7eDGY6ZqGUxARDFsxd6EkDPULuDYO5aR9SyvBkNlMhzXwLwWRH5HICPARiqqiwBUEtjAiCKICKf\nAfCyqt7sOhaiPLEKiCjaEAAzXAdBlDcmAKImRGR/AP8JYGsRibzJNlEr4Q1hiIhKiiUAIqKSYgIg\nIiopJgAiopJiAiAiKikmACKikmICICIqKSYAIqKSYgIgIiqp/weDV7A81tc4CQAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ce28cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEQCAYAAABBQVgLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD2RJREFUeJzt3WuMbWV5B/D/AyhIQa0XRKUiSBQvWC8pGkRjbJsQUyWm\nRolFU6tSUeuFavrByzmHNNa0Rm0UoaipaKp4q00LSvWDhyo1SkMVBestKGoVtRhEOYjn+PbDnsnM\nmTMz71z2nrVnn98v2Zl1mzXPPNkz//2utfba1VoLAKzmkKELAGD6CQsAuoQFAF3CAoAuYQFA12FD\nF7BeVeXyLYANaK3VRr93W44sWmtje+zYsWPs26+0zXLL17Js8fxK09PYi9XW64Ve6MXW9mKzDt25\nc+fm/3tvoV27du0cd80PfOADx779Stsst3wtyxbPz0/v3r07T37yk7u1rMe4e7Haer1YfZ1erG+Z\nXiw/Pz+9a9eu7Ny5c9eqxaxmnKm7FY9RybTW2o4dO4YuYWroxQK9WKAXC+b+d274f++2PAzFyLhf\nMW1nerFALxboxfhUG8OxrK1UVW271QwwtKpKO9hOcAOwtYQFAF3CAoAuYQFAl7AAoEtYAAyhteQN\nb0g+9amhK1kTl84CDOGyy5KnPW00vQX/01w6C7Ad3XTT0BWsi7AAoEtYANAlLADoEhYAdAkLALqE\nBQBdwgKALmEBQJewAKBLWADQJSwA6BIWAHQJC4AhbLO7ZwsLALqEBcAQasMfLTEIYQFAl7AAoEtY\nANAlLADoEhYAdAkLALqEBQBdwgKALmEBQJewAKBLWAAMwY0EAZg1wgKALmEBQJewABiCW5QDMGuE\nBQBdwgKALmEBQJewAKBLWADQJSwA6BIWAHQJCwC6hAXAENx1FoBZIywA6BIWAHQJCwC6hAXAENyi\nHIBZIywA6BIWAHQJCwC6hAUAXcICgC5hAUCXsACgS1gADMFdZwGYNcICgC5hAUCXsACgS1gA0CUs\nAIbgFuUAzBphAUCXsACgS1gA0CUsAOgSFgB0CQsAuoQFwBDcdRaAWSMsAOgSFgB0CQsAuoQFAF3C\nAoAuYQEwBLcoB2DWCAsAuoQFAF3CAoAuYQFAl7AAoEtYAAzBXWcBmDXCAoAuYQFAl7AAoEtYANAl\nLADoEhYAQ3DXWQBmjbAAoEtYANC14bCoqidU1QnjLAaA6XTYejauqtclOSnJbUk+keRRSS6YQF0A\nTJF1hUWS61prf11Vd03y1CTfnUBNALNvm91IcL1h0arq91prVye5dBIFATB91hsWT06SqtqR0aGo\n/2itvWPcRQEwXdYbFh9JktbaVVV1lyQPH39JAEybdYVFa+2qRdN7qur68ZcEwLTZ7PsszhlLFQBM\nte7IoqrekuRJSX4+vyhJm/t6cpK3Taw6AKbCWg5D/WWSV7bW3rp0RVW9cvwlATBtuoehWmstyXuT\npEYeWVW/M7f6XROsDYApsaYT3K21n81N/k2SQ5M8rKr2JvnzJL+cUG0As2ub3aJ8Lecsnpfkv5N8\nNclnW2uXzy2/e5KXJzl/ohUCMLi1jCzOzOj9FKckOb6qTk9yZZJvJ7llgrUBMCXWEhZnt9b2JElV\nHZrkd5M8NqMA+ecJ1gbAlOiGxXxQzE3vS3LN3AOAg4QPPwIYwja766ywAKBLWADQJSwA6BIWAHQJ\nCwC6hAUAXcICgC5hAUCXsACgS1gADGGb3aJcWADQJSwA6BIWAHQJC4AhuOssALNGWADQJSwA6BIW\nAHQJCwC6hAXALNi7N/n61ye2e2EBMAue8Yzk5JOTSy+dyO6FBcAsuOyy0ddLLpnI7oUFAF3CAmCW\nTOid4cICYAhuUQ7AYIwsABiKsAAYgrvOAjBrhAUAXcICYJY4wQ3AUIQFAF3CAmCWOAwFwFCEBcAs\nMbIAYCjCAmCWTOgGhcICYJY4DAUwQ9yiHIDBGFkAzBB3nQVg1ggLgFniMBQAQxEWAHQJCwC6hAUA\nXcICYJY4wQ3AUIQFwCwxsgBgKMICgC5hAUCXsACgS1gATKs9e5Jdu5KvfnXt3zN/gvu1r01OOy25\n446xlCIsAKbFm9+cXHHFwvyb3pTs3Jmccsr69/XGNyaf/3zy6U+PpbTDxrIXADbn6quT17xmND0/\nOrj++s3vd+/eze8jRhbAtPjQh5Irrxy6iuH89KcHLtvIR68ufZ/Fb36zsXqWMLIAhvf97ydnnTWa\n3mafIDc2hyzz2n0cn9M9prAwsgCGt9yr6oPNpMJi377N7yPCApgGi0cTRhbjZWQBzIw9examb799\nuDqG5DAUQMcvf7n89Hb3la8kb33r2g4FjSssnOAGZtbSsLjXvVbf/te/Tp7znOQpT0nOPXeytW3G\nIx85+nr00ckLX7jydrfemvz4xwcuPxhHFlX1v1XV5h4HXPhbVZcvWt+q6hdbVRswsMVhcdtt/e0v\nvzz56EeTl7xkcjWN03XXrb7+HvdInvWs1bd52cuSL31pbf1ZbDud4K6qI5LcN8mFSR6Q5NCq+pcl\nm/0syb4k30vyf0nushW1AVNgvYehfrHNXkv23hi3ljfOXXBB8uhHJw996OrbLT0MdfPNyRve0N9/\nx1YdhvpCkrTWXpIkVbUvyR8t2eaUJL9M8uUkJyW5Z1Xdp7V20wF7u/DCiRYLjNn84ZSVrnT6zGcW\npj/4wdGo4QEPWPkwzBe+sDD9zncuv93evcmXv5w85jEr7+eznx0dKrrb3Rbqu/LK5IlPTL73veRr\nX0se+9jkmGNG87fckjziESv/np//fPLgByf3vOf+y9/xjuQhDxndeuO000aHpS6+eOX9vOpVyd3v\nnnzgAweuu/HG5PzzR/d82rMn+fa3R3XOu+qq/X/fV7965Z+zDtW24DK1qro5yW+31mpu/tYkv9Va\nO2TRNj9Kcsz8bJLbkzyhtXbNkn0drBfWAWxYJZn/H7wR03SC++gkNyS5KslZSY5IclySa5ZuuPjY\n2fFHHZUbzj57SwoENuimm5Ijjlh4Bb+cG24YveK+611H0yeckNz5zitvv3t3cuKJoxHIcm69NfnR\nj5L73z858sgD1998c/K5z41GDY9//GjZpz89eqV+6qnJwx6WXHLJ6ET60UcnX/zi6MT7iScu//Ou\nuCL5zndG0y9+8ejrRReNvj772aPbmcybPym/3FGSBz0o+eEPk2OPTQ47LPnGN5b/eWecsf9NB5c4\nIcl3V1y7Aa21iT8yOrTUFs3vTbJ3yTa3J7kto/MWv0rSkrx3mX01gInZu3dj33f++a2NDmQtv/72\n21t7//tb+8lP9l8+/z3zj5//fP/1r3/9wrpXvGL09X3v23+biy8+cD9LHov/B2/ksVVXQz0pSarq\nnVV1UpJDk3xiyTbvTnJ4kn9MMv/e/7/dovoARg49dGPf1ztAfvjhydln9y8LXnp+ZfFo7G1vG/2c\n5z53/21OO+3A/ZxxRrJjx+o/ax22JCxaa7ck+WGSc5N8M8m+1trTq+r2qpp/u+bT5+p5QZL7JflN\nkhu3oj6ATdvo2dRjj119/bnnJs98ZvKxj628zcMfPrqd+ctfvrDsk5/cWD0r2LJzFq21+y2z7IhF\n0ysceATYBjYaFtdeO7qK6swzR/NLRxZHHpl85CP9/Tz0oaNzK4sdccTy226A230AjMNGw+Le9x4d\nMpqEl740Of301S/TXSNhATAOm7mof/FoYjO3+LjPffafP/ro0XtJXvSije9zzjRdOguwfU3DO8DO\nOWf0Br1nPGPsuxYWAOOwmRv2jWtkcfjho3e0T4DDUADjMA0jiwkSFgDjMA3nLCZIWACMg5EFAF1G\nFgAc7IQFwDg4DAVA1zRcOjtBwgJgHI46augKJsqb8gDG4bzzRjcEfP7z1/+922BksSUfqzpOVdW2\nW80AXfMhcccdyZ3uNIHd16Y+VtVhKAC6hAUAXcICgC5hATBNNvoZ4BPmaiiAafCDHyT79iWHTOdr\neFdDARwEXA0FwMQJCwC6hMU2tnv37qFLmBp6sUAvFujF+AiLbcwfwgK9WKAXC/RifA76sFjvk2kt\n26+0zXLL17Js8fwkn/zj7sVq6/Vi9XV6sb5lerH8/Dh7ISyExYb3PUt/CGupZTPb68Xa1uvF6uuG\n7MW2vHR26BoAtqPNXDq77cICgK130B+GAqBPWADQJSwA6BIWAHQJCwC6ZiYsqurIqrq6qp46dC1D\nqqqTq+rCqvpwVb146HqGVFVnVtXFVfXBqvrDoesZUlWdUFXvrqoPD13LkOb+T7y3qv6hqp4zdD1D\nWu9zYmYuna2qXUluTXJ9a+0TQ9cztKqqJJe01p43dC1Dq6q7J/m71tqLhq5laFX14dbas4auYyhV\ndXaSn7XWLq+qS1trZw1d09DW+pyYqpFFVb2nqm6qqmuXLD+jqv6nqr5RVX+1zPf9QZLrk/wkyYbf\ndDJNNtqLuW2eluSyJDMRmpvpxZzXJblgslVujTH0YqZsoB/HJfne3PS+LSt0C0z8udFam5pHktOT\nPCrJtYuWHZLkW0mOT3KnJF9KcvLcuucmeWuS9yR5S5J/T/LxoX+PAXvxliT3XbT9ZUP/HgP34n5J\n3pTkKUP/DlPQi/vOzX9k6N9h4H78SZKnzk1/YOj6h+zFom3W9JyYqpFFa+1zSX62ZPGpSb7ZWvtu\na+3XSS5Ncubc9u9vrb2qtfaC1tp5Sf4pybu2tOgJ2WAvzkvy4Kr6+6q6KMnlW1r0hGyiF3+c5PeT\nPLOqztnKmidlE734VVVdmORRszTyWG8/knw8o+fDBUn+besqnbz19qKq7rGe58R2+Azu+2dh2Jgk\n38+oAQdorb1vSyoaTrcXrbUrk1y5lUUNZC29eHuSt29lUQNZSy9uTnLuVhY1oBX70Vq7LcmfDVHU\nQFbrxbqeE1M1sgBgOm2HsPhBkgcsmj9ubtnBSC8W6MUCvdiffiwYWy+mMSwq+1/RdHWSk6rq+Kq6\nc5KzkvzrIJVtPb1YoBcL9GJ/+rFgYr2YqrCoqg8k+c+MTtLeWFXPb63tS/IXST6V5Lokl7bWvjZk\nnVtBLxboxQK92J9+LJh0L2bmTXkATM5UjSwAmE7CAoAuYQFAl7AAoEtYANAlLADoEhYAdAkLALqE\nBQBdwgLGoKpeV1XfqqoXVNVLq+qiqjpq6LpgXIQFjMd/Jbmitfae1toFSW5J8qfDlgTjIyxgPB6X\n5JpF88cluW2gWmDshAWMx6mZC4uqOjajTyi7dNCKYIy2w8eqwnbw6CQPqqqTkhyd5KlzH+EJM0FY\nwCZV1YlJvtta+9jQtcCkOAwFm/e4JFcNXQRMkrCATaiqJyU5J8kxVXXM0PXApPikPAC6jCwA6BIW\nAHQJCwC6hAUAXcICgC5hAUCXsACg6/8B2sue+sBC50wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10ce87518>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig=figure(1)\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(plist,En,'r',lw=2)\n",
    "ylabel('$E_N$')\n",
    "xlabel('$P$')\n",
    "show()\n",
    "fig=figure(2)\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(plist,var,'r',lw=2)\n",
    "yticks = [0,0.1, 0.2, 0.3,0.4,0.5,0.6]\n",
    "ax.set_yticks(yticks)\n",
    "plt.xscale('log')\n",
    "ylabel('$b1n$')\n",
    "xlabel('$P$')\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [qutip]",
   "language": "python",
   "name": "Python [qutip]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}