notebooks/bi_00_adiabatic_elimination.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adiabatic Elimination in the Simple Bidirectionally-coupled Configuration (Bi_00)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dependencies\n",
    "from IPython.display import display, Math\n",
    "from sympy import init_printing, latex\n",
    "from sympy import conjugate, diff, E, Function, I, Matrix, solve, sqrt, Symbol, symbols\n",
    "from sympy.physics.quantum import Dagger\n",
    "# use mathjax renderer for faster load times\n",
    "init_printing(use_latex='mathjax')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# time\n",
    "t = symbols('t', real=True, positive=True)\n",
    "# positive parameters\n",
    "g_L, g_R, gamma_L, gamma_R, kappa_L, kappa_R, omega_L, omega_R = symbols('g_{L}, g_{R}, \\\\gamma_{L}, \\\\gamma_{R}, \\\\kappa_{L}, \\\\kappa_{R}, \\\\omega_{L}, \\\\omega_{R}', real=True, positive=True)\n",
    "# real-valued parameters\n",
    "Delta_L, Delta_R, lamb = symbols('\\\\Delta_{L}, \\\\Delta_{R}, \\\\lambda', real=True) \n",
    "# transformed parameters\n",
    "Gamma_L, Gamma_R, gamma, kappa, n_th = symbols('\\\\Gamma_{L}, \\\\Gamma_{R}, \\\\gamma, \\\\kappa, n_{th}', real=True, positive=True)\n",
    "delta = symbols('delta', real=True)\n",
    "chi, eta_L, eta_R, G_L, G_R = symbols('\\\\chi, \\\\eta_{L}, \\\\eta_{R}, G_{L}, G_{R}', complex=True)\n",
    "chi_R, chi_I, eta_LR, eta_LI, eta_RR, eta_RI = symbols('\\\\chi_{R}, \\\\chi_{I}, \\\\eta_{L_{R}}, \\\\eta_{L_{I}}, \\\\eta_{R_{R}}, \\\\eta_{R_{I}}', real=True)\n",
    "\n",
    "# classical amplitudes\n",
    "alpha_L, alpha_R, beta_L, beta_R = symbols('\\\\alpha_{L}, \\\\alpha_{R}, \\\\beta_{L}, \\\\beta_{R}', complex=True)\n",
    "# quantum fluctuations\n",
    "a_L_t = Function('\\\\hat{a}_{L}', commutative=False)(t)\n",
    "b_L_t = Function('\\\\hat{b}_{L}', commutative=False)(t)\n",
    "a_R_t = Function('\\\\hat{a}_{R}', commutative=False)(t)\n",
    "b_R_t = Function('\\\\hat{b}_{R}', commutative=False)(t)\n",
    "# input noises\n",
    "a_L_in_t = Function('\\\\hat{a}_{L}^{in}', commutative=False)(t)\n",
    "b_L_in_t = Function('\\\\hat{b}_{L}^{in}', commutative=False)(t)\n",
    "a_R_in_t = Function('\\\\hat{a}_{R}^{in}', commutative=False)(t)\n",
    "b_R_in_t = Function('\\\\hat{b}_{R}^{in}', commutative=False)(t)\n",
    "\n",
    "# transformed fluctuations\n",
    "a_L_tilde_t = Function('\\\\hat{\\\\tilde{a}}_{L}', commutative=False)(t)\n",
    "b_L_tilde_t = Function('\\\\hat{\\\\tilde{b}}_{L}', commutative=False)(t)\n",
    "a_R_tilde_t = Function('\\\\hat{\\\\tilde{a}}_{R}', commutative=False)(t)\n",
    "b_R_tilde_t = Function('\\\\hat{\\\\tilde{b}}_{R}', commutative=False)(t)\n",
    "# transformed noises# new bath operators\n",
    "a_L_in_tilde_t = Function('\\\\hat{\\\\tilde{a}}_{L}^{in}', commutative=False)(t)\n",
    "b_L_in_tilde_t = Function('\\\\hat{\\\\tilde{b}}_{L}^{in}', commutative=False)(t)\n",
    "a_R_in_tilde_t = Function('\\\\hat{\\\\tilde{a}}_{R}^{in}', commutative=False)(t)\n",
    "b_R_in_tilde_t = Function('\\\\hat{\\\\tilde{b}}_{R}^{in}', commutative=False)(t)\n",
    "\n",
    "# rotated mechanical quadratures\n",
    "Q_L_tilde_t = Function('\\\\hat{\\\\tilde{Q}}_{L}', commutative=False)(t)\n",
    "P_L_tilde_t = Function('\\\\hat{\\\\tilde{P}}_{L}', commutative=False)(t)\n",
    "Q_R_tilde_t = Function('\\\\hat{\\\\tilde{Q}}_{R}', commutative=False)(t)\n",
    "P_R_tilde_t = Function('\\\\hat{\\\\tilde{P}}_{R}', commutative=False)(t)\n",
    "# rotated noise quadratures\n",
    "X_L_in_tilde_t = Function('\\\\hat{\\\\tilde{X}}_{L}^{in}', commutative=False)(t)\n",
    "Y_L_in_tilde_t = Function('\\\\hat{\\\\tilde{Y}}_{L}^{in}', commutative=False)(t)\n",
    "Q_L_in_tilde_t = Function('\\\\hat{\\\\tilde{Q}}_{L}^{in}', commutative=False)(t)\n",
    "P_L_in_tilde_t = Function('\\\\hat{\\\\tilde{P}}_{L}^{in}', commutative=False)(t)\n",
    "X_R_in_tilde_t = Function('\\\\hat{\\\\tilde{X}}_{R}^{in}', commutative=False)(t)\n",
    "Y_R_in_tilde_t = Function('\\\\hat{\\\\tilde{Y}}_{R}^{in}', commutative=False)(t)\n",
    "Q_R_in_tilde_t = Function('\\\\hat{\\\\tilde{Q}}_{R}^{in}', commutative=False)(t)\n",
    "P_R_in_tilde_t = Function('\\\\hat{\\\\tilde{P}}_{R}^{in}', commutative=False)(t)\n",
    "\n",
    "# correlation expectation values\n",
    "b_L_tilde_dagger_b_L_tilde_expect_t = Function('\\\\langle \\\\hat{\\\\tilde{b}}_{L}^{\\\\dagger} \\\\hat{\\\\tilde{b}}_{L} \\\\rangle', complex=True)(t)\n",
    "b_L_tilde_dagger_b_R_tilde_expect_t = Function('\\\\langle \\\\hat{\\\\tilde{b}}_{L}^{\\\\dagger} \\\\hat{\\\\tilde{b}}_{R} \\\\rangle', complex=True)(t)\n",
    "b_R_tilde_dagger_b_L_tilde_expect_t = Function('\\\\langle \\\\hat{\\\\tilde{b}}_{R}^{\\\\dagger} \\\\hat{\\\\tilde{b}}_{L} \\\\rangle', complex=True)(t)\n",
    "b_R_tilde_dagger_b_R_tilde_expect_t = Function('\\\\langle \\\\hat{\\\\tilde{b}}_{R}^{\\\\dagger} \\\\hat{\\\\tilde{b}}_{R} \\\\rangle', complex=True)(t)\n",
    "b_L_tilde_dagger_b_L_tilde_expect = Symbol('\\\\langle \\\\hat{\\\\tilde{b}}_{L}^{\\\\dagger} \\\\hat{\\\\tilde{b}}_{L} \\\\rangle', complex=True)\n",
    "b_L_tilde_dagger_b_R_tilde_expect = Symbol('\\\\langle \\\\hat{\\\\tilde{b}}_{L}^{\\\\dagger} \\\\hat{\\\\tilde{b}}_{R} \\\\rangle', complex=True)\n",
    "b_R_tilde_dagger_b_L_tilde_expect = Symbol('\\\\langle \\\\hat{\\\\tilde{b}}_{R}^{\\\\dagger} \\\\hat{\\\\tilde{b}}_{L} \\\\rangle', complex=True)\n",
    "b_R_tilde_dagger_b_R_tilde_expect = Symbol('\\\\langle \\\\hat{\\\\tilde{b}}_{R}^{\\\\dagger} \\\\hat{\\\\tilde{b}}_{R} \\\\rangle', complex=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantum Langevin Equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{a}}_{L} (t) = i \\alpha_{L} g_{L} \\left(\\hat{b}_{L}{\\left(t \\right)} + \\hat{b}_{L}{\\left(t \\right)}^{\\dagger}\\right) + \\sqrt{2} \\sqrt{\\kappa_{L}} \\hat{a}_{L}^{in}{\\left(t \\right)} + i \\lambda \\hat{a}_{R}{\\left(t \\right)} + \\left(i \\Delta_{L} - \\kappa_{L}\\right) \\hat{a}_{L}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{b}}_{L} (t) = \\sqrt{2} \\sqrt{\\gamma_{L}} \\hat{b}_{L}^{in}{\\left(t \\right)} + i g_{L} \\left(\\alpha_{L} \\hat{a}_{L}{\\left(t \\right)}^{\\dagger} + \\overline{\\alpha_{L}} \\hat{a}_{L}{\\left(t \\right)}\\right) + \\left(- \\gamma_{L} - i \\omega_{L}\\right) \\hat{b}_{L}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{a}}_{R} (t) = i \\alpha_{R} g_{R} \\left(\\hat{b}_{R}{\\left(t \\right)} + \\hat{b}_{R}{\\left(t \\right)}^{\\dagger}\\right) + \\sqrt{2} \\sqrt{\\kappa_{R}} \\hat{a}_{R}^{in}{\\left(t \\right)} + i \\lambda \\hat{a}_{L}{\\left(t \\right)} + \\left(i \\Delta_{R} - \\kappa_{R}\\right) \\hat{a}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{b}}_{R} (t) = \\sqrt{2} \\sqrt{\\gamma_{R}} \\hat{b}_{R}^{in}{\\left(t \\right)} + i g_{R} \\left(\\alpha_{R} \\hat{a}_{R}{\\left(t \\right)}^{\\dagger} + \\overline{\\alpha_{R}} \\hat{a}_{R}{\\left(t \\right)}\\right) + \\left(- \\gamma_{R} - i \\omega_{R}\\right) \\hat{b}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# rate equations\n",
    "expr_da_L_dt_t = (- kappa_L + I * Delta_L) * a_L_t + I * g_L * alpha_L * (Dagger(b_L_t) + b_L_t) + I * lamb * a_R_t + sqrt(2 * kappa_L) * a_L_in_t\n",
    "display(Math('\\\\dot{\\\\hat{a}}_{L} (t) = ' + latex(expr_da_L_dt_t)))\n",
    "expr_db_L_dt_t = (- gamma_L - I * omega_L) * b_L_t + I * g_L * (conjugate(alpha_L) * a_L_t + Dagger(a_L_t) * alpha_L) + sqrt(2 * gamma_L) * b_L_in_t\n",
    "display(Math('\\\\dot{\\\\hat{b}}_{L} (t) = ' + latex(expr_db_L_dt_t)))\n",
    "expr_da_R_dt_t = (- kappa_R + I * Delta_R) * a_R_t + I * g_R * alpha_R * (Dagger(b_R_t) + b_R_t) + I * lamb * a_L_t + sqrt(2 * kappa_R) * a_R_in_t\n",
    "display(Math('\\\\dot{\\\\hat{a}}_{R} (t) = ' + latex(expr_da_R_dt_t)))\n",
    "expr_db_R_dt_t = (- gamma_R - I * omega_R) * b_R_t + I * g_R * (conjugate(alpha_R) * a_R_t + Dagger(a_R_t) * alpha_R) + sqrt(2 * gamma_R) * b_R_in_t\n",
    "display(Math('\\\\dot{\\\\hat{b}}_{R} (t) = ' + latex(expr_db_R_dt_t)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transformations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\hat{\\tilde{a}}_{L}{\\left(t \\right)} = e^{- i \\omega_{L} t} \\hat{a}_{L}{\\left(t \\right)}, \\quad \\hat{\\tilde{b}}_{L}{\\left(t \\right)} = e^{i \\omega_{L} t} \\hat{b}_{L}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\hat{\\tilde{a}}_{R}{\\left(t \\right)} = e^{- i \\omega_{L} t} \\hat{a}_{R}{\\left(t \\right)}, \\quad \\hat{\\tilde{b}}_{R}{\\left(t \\right)} = e^{i \\omega_{L} t} \\hat{b}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\Delta_{L} = \\omega_{L}, \\quad \\Delta_{R} = \\omega_{L}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)} = e^{- i \\omega_{L} t} \\hat{a}_{L}^{in}{\\left(t \\right)}, \\quad \\hat{\\tilde{b}}_{L}^{in}{\\left(t \\right)} = e^{i \\omega_{L} t} \\hat{b}_{L}^{in}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)} = e^{- i \\omega_{L} t} \\hat{a}_{R}^{in}{\\left(t \\right)}, \\quad \\hat{\\tilde{b}}_{R}^{in}{\\left(t \\right)} = e^{i \\omega_{L} t} \\hat{b}_{R}^{in}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{a}}}_{L} (t) = i \\alpha_{L} g_{L} \\hat{\\tilde{b}}_{L}{\\left(t \\right)}^{\\dagger} + i \\alpha_{L} g_{L} e^{- 2 i \\omega_{L} t} \\hat{\\tilde{b}}_{L}{\\left(t \\right)} + \\sqrt{2} \\sqrt{\\kappa_{L}} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)} - \\kappa_{L} \\hat{\\tilde{a}}_{L}{\\left(t \\right)} + i \\lambda \\hat{\\tilde{a}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{b}}}_{L} (t) = i \\alpha_{L} g_{L} \\hat{\\tilde{a}}_{L}{\\left(t \\right)}^{\\dagger} + \\sqrt{2} \\sqrt{\\gamma_{L}} \\hat{\\tilde{b}}_{L}^{in}{\\left(t \\right)} - \\gamma_{L} \\hat{\\tilde{b}}_{L}{\\left(t \\right)} + i g_{L} e^{2 i \\omega_{L} t} \\overline{\\alpha_{L}} \\hat{\\tilde{a}}_{L}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{a}}}_{R} (t) = i \\alpha_{R} g_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)}^{\\dagger} + i \\alpha_{R} g_{R} e^{- 2 i \\omega_{L} t} \\hat{\\tilde{b}}_{R}{\\left(t \\right)} + \\sqrt{2} \\sqrt{\\kappa_{R}} \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)} - \\kappa_{R} \\hat{\\tilde{a}}_{R}{\\left(t \\right)} + i \\lambda \\hat{\\tilde{a}}_{L}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{b}}}_{R} (t) = i \\alpha_{R} g_{R} \\hat{\\tilde{a}}_{R}{\\left(t \\right)}^{\\dagger} + \\sqrt{2} \\sqrt{\\gamma_{R}} \\hat{\\tilde{b}}_{R}^{in}{\\left(t \\right)} - \\gamma_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)} + i \\omega_{L} \\hat{\\tilde{b}}_{R}{\\left(t \\right)} - i \\omega_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)} + i g_{R} e^{2 i \\omega_{L} t} \\overline{\\alpha_{R}} \\hat{\\tilde{a}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# fluctuation transformations\n",
    "expr_a_L_tilde_t = a_L_t * E**(- I * omega_L * t)\n",
    "expr_b_L_tilde_t = b_L_t * E**(I * omega_L * t)\n",
    "expr_a_R_tilde_t = a_R_t * E**(- I * omega_L * t)\n",
    "expr_b_R_tilde_t = b_R_t * E**(I * omega_L * t)\n",
    "# noise transformations\n",
    "expr_a_L_in_tilde_t = a_L_in_t * E**(- I * omega_L * t)\n",
    "expr_b_L_in_tilde_t = b_L_in_t * E**(I * omega_L * t)\n",
    "expr_a_R_in_tilde_t = a_R_in_t * E**(- I * omega_L * t)\n",
    "expr_b_R_in_tilde_t = b_R_in_t * E**(I * omega_L * t)\n",
    "# substitution list\n",
    "list_subs = [\n",
    "    (a_L_t, a_L_tilde_t * a_L_t / expr_a_L_tilde_t),\n",
    "    (b_L_t, b_L_tilde_t * b_L_t / expr_b_L_tilde_t), \n",
    "    (a_R_t, a_R_tilde_t * a_R_t / expr_a_R_tilde_t),\n",
    "    (b_R_t, b_R_tilde_t * b_R_t / expr_b_R_tilde_t),\n",
    "    (expr_a_L_in_tilde_t, a_L_in_tilde_t),\n",
    "    (expr_b_L_in_tilde_t, b_L_in_tilde_t),\n",
    "    (expr_a_R_in_tilde_t, a_R_in_tilde_t),\n",
    "    (expr_b_R_in_tilde_t, b_R_in_tilde_t),\n",
    "    (Delta_L, omega_L), \n",
    "    (Delta_R, omega_L),\n",
    "]\n",
    "# display\n",
    "display(Math(latex(a_L_tilde_t) + ' = ' + latex(expr_a_L_tilde_t) + ', \\\\quad ' + latex(b_L_tilde_t) + ' = ' + latex(expr_b_L_tilde_t)))\n",
    "display(Math(latex(a_R_tilde_t) + ' = ' + latex(expr_a_R_tilde_t) + ', \\\\quad ' + latex(b_R_tilde_t) + ' = ' + latex(expr_b_R_tilde_t)))\n",
    "display(Math(latex(Delta_L) + ' = ' + latex(omega_L) + ', \\\\quad ' + latex(Delta_R) + ' = ' + latex(omega_L)))\n",
    "display(Math(latex(a_L_in_tilde_t) + ' = ' + latex(expr_a_L_in_tilde_t) + ', \\\\quad ' + latex(b_L_in_tilde_t) + ' = ' + latex(expr_b_L_in_tilde_t)))\n",
    "display(Math(latex(a_R_in_tilde_t) + ' = ' + latex(expr_a_R_in_tilde_t) + ', \\\\quad ' + latex(b_R_in_tilde_t) + ' = ' + latex(expr_b_R_in_tilde_t)))\n",
    "\n",
    "# transformed equations\n",
    "# differentiate and substitute rates in first step \n",
    "# substitute other variables in second step\n",
    "expr_da_L_tilde_dt_t = diff(expr_a_L_tilde_t, t).subs(diff(a_L_t), expr_da_L_dt_t).expand()\n",
    "expr_da_L_tilde_dt_t = expr_da_L_tilde_dt_t.subs(expr_a_L_tilde_t, a_L_tilde_t).subs(list_subs)\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{a}}}_{L} (t) = ' + latex(expr_da_L_tilde_dt_t)))\n",
    "expr_db_L_tilde_dt_t = diff(expr_b_L_tilde_t, t).subs(diff(b_L_t), expr_db_L_dt_t).expand()\n",
    "expr_db_L_tilde_dt_t = expr_db_L_tilde_dt_t.subs(expr_b_L_tilde_t, b_L_tilde_t).subs(list_subs)\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{b}}}_{L} (t) = ' + latex(expr_db_L_tilde_dt_t)))\n",
    "expr_da_R_tilde_dt_t = diff(expr_a_R_tilde_t, t).subs(diff(a_R_t), expr_da_R_dt_t).expand()\n",
    "expr_da_R_tilde_dt_t = expr_da_R_tilde_dt_t.subs(expr_a_R_tilde_t, a_R_tilde_t).subs(list_subs)\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{a}}}_{R} (t) = ' + latex(expr_da_R_tilde_dt_t)))\n",
    "expr_db_R_tilde_dt_t = diff(expr_b_R_tilde_t, t).subs(diff(b_R_t), expr_db_R_dt_t).expand()\n",
    "expr_db_R_tilde_dt_t = expr_db_R_tilde_dt_t.subs(expr_b_R_tilde_t, b_R_tilde_t).subs(list_subs)\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{b}}}_{R} (t) = ' + latex(expr_db_R_tilde_dt_t)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle G_{L} = \\alpha_{L} g_{L}, \\quad G_{R} = \\alpha_{R} g_{R}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\delta = - \\omega_{L} + \\omega_{R}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{a}}}_{L} (t) = i G_{L} \\hat{\\tilde{b}}_{L}{\\left(t \\right)}^{\\dagger} + \\sqrt{2} \\sqrt{\\kappa_{L}} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)} - \\kappa_{L} \\hat{\\tilde{a}}_{L}{\\left(t \\right)} + i \\lambda \\hat{\\tilde{a}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{b}}}_{L} (t) = i G_{L} \\hat{\\tilde{a}}_{L}{\\left(t \\right)}^{\\dagger} + \\sqrt{2} \\sqrt{\\gamma_{L}} \\hat{\\tilde{b}}_{L}^{in}{\\left(t \\right)} - \\gamma_{L} \\hat{\\tilde{b}}_{L}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{a}}}_{R} (t) = i G_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)}^{\\dagger} + \\sqrt{2} \\sqrt{\\kappa_{R}} \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)} - \\kappa_{R} \\hat{\\tilde{a}}_{R}{\\left(t \\right)} + i \\lambda \\hat{\\tilde{a}}_{L}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{b}}}_{R} (t) = i G_{R} \\hat{\\tilde{a}}_{R}{\\left(t \\right)}^{\\dagger} + \\sqrt{2} \\sqrt{\\gamma_{R}} \\hat{\\tilde{b}}_{R}^{in}{\\left(t \\right)} - \\gamma_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)} - i \\delta \\hat{\\tilde{b}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# effective values\n",
    "expr_delta = omega_R - omega_L\n",
    "expr_G_L = g_L * alpha_L\n",
    "expr_G_R = g_R * alpha_R\n",
    "# substitution list\n",
    "list_subs = [\n",
    "    (omega_R, delta + omega_L),\n",
    "    (expr_G_L, G_L), \n",
    "    (expr_G_R, G_R)\n",
    "]\n",
    "# display\n",
    "display(Math(latex(G_L) + ' = ' + latex(expr_G_L) + ', \\\\quad ' + latex(G_R) + ' = ' + latex(expr_G_R)))\n",
    "display(Math(latex(delta) + ' = ' + latex(expr_delta)))\n",
    "\n",
    "# approximated equations\n",
    "# substitute parameters and ignore fast rotating terms\n",
    "expr_da_L_tilde_dt_approx_t = expr_da_L_tilde_dt_t.subs(list_subs).expand().subs(E**(- 2 * I * omega_L * t), 0)\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{a}}}_{L} (t) = ' + latex(expr_da_L_tilde_dt_approx_t)))\n",
    "expr_db_L_tilde_dt_approx_t = expr_db_L_tilde_dt_t.subs(list_subs).expand().subs(E**(2 * I * omega_L * t), 0)\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{b}}}_{L} (t) = ' + latex(expr_db_L_tilde_dt_approx_t)))\n",
    "expr_da_R_tilde_dt_approx_t = expr_da_R_tilde_dt_t.subs(list_subs).expand().subs(E**(- 2 * I * omega_L * t), 0)\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{a}}}_{R} (t) = ' + latex(expr_da_R_tilde_dt_approx_t)))\n",
    "expr_db_R_tilde_dt_approx_t = expr_db_R_tilde_dt_t.subs(list_subs).expand().subs(E**(2 * I * omega_L * t), 0)\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{b}}}_{R} (t) = ' + latex(expr_db_R_tilde_dt_approx_t)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adiabatic Elimination"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\hat{\\tilde{a}}_{L} = \\frac{i G_{L} \\kappa_{R} \\hat{\\tilde{b}}_{L}{\\left(t \\right)}^{\\dagger}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} - \\frac{G_{R} \\lambda \\hat{\\tilde{b}}_{R}{\\left(t \\right)}^{\\dagger}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} + \\frac{\\sqrt{2} \\sqrt{\\kappa_{L}} \\kappa_{R} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} + \\frac{\\sqrt{2} i \\sqrt{\\kappa_{R}} \\lambda \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\hat{\\tilde{a}}_{R} = - \\frac{G_{L} \\lambda \\hat{\\tilde{b}}_{L}{\\left(t \\right)}^{\\dagger}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} + \\frac{i G_{R} \\kappa_{L} \\hat{\\tilde{b}}_{R}{\\left(t \\right)}^{\\dagger}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} + \\frac{\\sqrt{2} i \\sqrt{\\kappa_{L}} \\lambda \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} + \\frac{\\sqrt{2} \\kappa_{L} \\sqrt{\\kappa_{R}} \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{b}}}_{L} (t) = \\frac{\\sqrt{2} i G_{L} \\sqrt{\\kappa_{L}} \\kappa_{R} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} + \\frac{\\sqrt{2} G_{L} \\sqrt{\\kappa_{R}} \\lambda \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}^{\\dagger}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} + \\frac{G_{L} \\kappa_{R} \\overline{G_{L}} \\hat{\\tilde{b}}_{L}{\\left(t \\right)}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} - \\frac{i G_{L} \\lambda \\overline{G_{R}} \\hat{\\tilde{b}}_{R}{\\left(t \\right)}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} + \\sqrt{2} \\sqrt{\\gamma_{L}} \\hat{\\tilde{b}}_{L}^{in}{\\left(t \\right)} - \\gamma_{L} \\hat{\\tilde{b}}_{L}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{b}}}_{R} (t) = \\frac{\\sqrt{2} G_{R} \\sqrt{\\kappa_{L}} \\lambda \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} + \\frac{\\sqrt{2} i G_{R} \\kappa_{L} \\sqrt{\\kappa_{R}} \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}^{\\dagger}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} + \\frac{G_{R} \\kappa_{L} \\overline{G_{R}} \\hat{\\tilde{b}}_{R}{\\left(t \\right)}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} - \\frac{i G_{R} \\lambda \\overline{G_{L}} \\hat{\\tilde{b}}_{L}{\\left(t \\right)}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} + \\sqrt{2} \\sqrt{\\gamma_{R}} \\hat{\\tilde{b}}_{R}^{in}{\\left(t \\right)} - \\gamma_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)} - i \\delta \\hat{\\tilde{b}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# obtain solutions for optical modes\n",
    "sols = solve([expr_da_R_tilde_dt_approx_t, expr_da_L_tilde_dt_approx_t], [a_R_tilde_t, a_L_tilde_t])\n",
    "expr_a_L_tilde_adia_t = sols[0][1].expand()\n",
    "display(Math('\\\\hat{\\\\tilde{a}}_{L} = ' + latex(expr_a_L_tilde_adia_t)))\n",
    "expr_a_R_tilde_adia_t = sols[0][0].expand()\n",
    "display(Math('\\\\hat{\\\\tilde{a}}_{R} = ' + latex(expr_a_R_tilde_adia_t)))\n",
    "\n",
    "# substitution list\n",
    "list_subs = [\n",
    "    (a_L_tilde_t, expr_a_L_tilde_adia_t), \n",
    "    (a_R_tilde_t, expr_a_R_tilde_adia_t)\n",
    "]\n",
    "# mechanical mode equations\n",
    "expr_db_L_tilde_dt_adia_t = expr_db_L_tilde_dt_approx_t.subs(list_subs).expand()\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{b}}}_{L} (t) = ' + latex(expr_db_L_tilde_dt_adia_t)))\n",
    "expr_db_R_tilde_dt_adia_t = expr_db_R_tilde_dt_approx_t.subs(list_subs).expand()\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{b}}}_{R} (t) = ' + latex(expr_db_R_tilde_dt_adia_t)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\chi = \\frac{G_{R} \\lambda \\overline{G_{L}}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\Gamma_{L} = \\frac{G_{L} \\kappa_{R} \\overline{G_{L}}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}}, \\quad \\Gamma_{R} = \\frac{G_{R} \\kappa_{L} \\overline{G_{R}}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\eta_{L} = \\frac{\\sqrt{2} G_{L} \\sqrt{\\kappa_{L}} \\kappa_{R}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}}, \\quad \\eta_{R} = \\frac{\\sqrt{2} G_{R} \\kappa_{L} \\sqrt{\\kappa_{R}}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{b}}}_{1} (t) = \\Gamma_{L} \\hat{\\tilde{b}}_{L}{\\left(t \\right)} + \\eta_{L} \\left(i \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger} + \\frac{\\lambda \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}^{\\dagger}}{\\sqrt{\\kappa_{L}} \\sqrt{\\kappa_{R}}}\\right) + \\sqrt{2} \\sqrt{\\gamma_{L}} \\hat{\\tilde{b}}_{L}^{in}{\\left(t \\right)} - \\gamma_{L} \\hat{\\tilde{b}}_{L}{\\left(t \\right)} - i \\overline{\\chi} \\hat{\\tilde{b}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{b}}}_{R} (t) = \\Gamma_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)} - i \\chi \\hat{\\tilde{b}}_{L}{\\left(t \\right)} + \\eta_{R} \\left(i \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}^{\\dagger} + \\frac{\\lambda \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger}}{\\sqrt{\\kappa_{L}} \\sqrt{\\kappa_{R}}}\\right) + \\sqrt{2} \\sqrt{\\gamma_{R}} \\hat{\\tilde{b}}_{R}^{in}{\\left(t \\right)} - \\gamma_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)} - i \\delta \\hat{\\tilde{b}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# substituted values\n",
    "expr_chi = lamb * conjugate(G_L) * G_R / (kappa_L * kappa_R + lamb**2)\n",
    "expr_eta_L = sqrt(2) * G_L * sqrt(kappa_L) * kappa_R / (kappa_L * kappa_R + lamb**2)\n",
    "expr_eta_R = sqrt(2) * G_R * kappa_L * sqrt(kappa_R) / (kappa_L * kappa_R + lamb**2)\n",
    "expr_Gamma_L = G_L * kappa_R * conjugate(G_L) / (kappa_L * kappa_R + lamb**2)\n",
    "expr_Gamma_R = G_R * kappa_L * conjugate(G_R) / (kappa_L * kappa_R + lamb**2)\n",
    "# substitution list\n",
    "list_subs = [\n",
    "    (expr_chi, chi),\n",
    "    (conjugate(expr_chi), conjugate(chi)),\n",
    "    (expr_Gamma_L, Gamma_L),\n",
    "    (expr_Gamma_R, Gamma_R),\n",
    "    (expr_eta_L / sqrt(kappa_L) / kappa_R, eta_L / sqrt(kappa_L) / kappa_R),\n",
    "    (expr_eta_R / kappa_L / sqrt(kappa_R), eta_R / kappa_L / sqrt(kappa_R))\n",
    "]\n",
    "# display\n",
    "display(Math(latex(chi) + ' = ' + latex(expr_chi)))\n",
    "display(Math(latex(Gamma_L) + ' = ' + latex(expr_Gamma_L) + ', \\\\quad ' + latex(Gamma_R) + ' = ' + latex(expr_Gamma_R)))\n",
    "display(Math(latex(eta_L) + ' = ' + latex(expr_eta_L) + ', \\\\quad ' + latex(eta_R) + ' = ' + latex(expr_eta_R)))\n",
    "\n",
    "# substituted equations\n",
    "expr_db_L_tilde_dt_subs_t = expr_db_L_tilde_dt_adia_t.collect(E**(I * delta * t)).subs(list_subs).expand()\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{b}}}_{1} (t) = ' + latex(expr_db_L_tilde_dt_subs_t.collect(eta_L))))\n",
    "expr_db_R_tilde_dt_subs_t = expr_db_R_tilde_dt_adia_t.collect(E**(- I * delta * t)).subs(list_subs).expand()\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{b}}}_{R} (t) = ' + latex(expr_db_R_tilde_dt_subs_t.collect(eta_R))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Identical Cavities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\gamma_{L} = \\gamma_{R} = \\gamma, \\quad \\kappa_{L} = \\kappa_{R} = \\kappa$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{b}}}_{L} (t) = \\Gamma_{L} \\hat{\\tilde{b}}_{L}{\\left(t \\right)} + \\eta_{L} \\left(i \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger} + \\frac{\\lambda \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}^{\\dagger}}{\\kappa}\\right) + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{b}}_{L}^{in}{\\left(t \\right)} - \\gamma \\hat{\\tilde{b}}_{L}{\\left(t \\right)} - i \\overline{\\chi} \\hat{\\tilde{b}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{b}}}_{R} (t) = \\Gamma_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)} - i \\chi \\hat{\\tilde{b}}_{L}{\\left(t \\right)} + \\eta_{R} \\left(i \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}^{\\dagger} + \\frac{\\lambda \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger}}{\\kappa}\\right) + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{b}}_{R}^{in}{\\left(t \\right)} - \\gamma \\hat{\\tilde{b}}_{R}{\\left(t \\right)} - i \\delta \\hat{\\tilde{b}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# substitution list\n",
    "list_subs = [\n",
    "    (gamma_L, gamma),\n",
    "    (gamma_R, gamma),\n",
    "    (kappa_L, kappa),\n",
    "    (kappa_R, kappa)\n",
    "]\n",
    "# display\n",
    "display(Math(latex(gamma_L) + ' = ' + latex(gamma_R) + ' = ' + latex(gamma) + ', \\\\quad ' + latex(kappa_L) + ' = ' + latex(kappa_R) + ' = ' + latex(kappa)))\n",
    "\n",
    "# substituted equations\n",
    "expr_db_L_tilde_dt_identical_t = expr_db_L_tilde_dt_subs_t.subs(list_subs)\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{b}}}_{L} (t) = ' + latex(expr_db_L_tilde_dt_identical_t.collect(eta_L))))\n",
    "expr_db_R_tilde_dt_identical_t = expr_db_R_tilde_dt_subs_t.subs(list_subs)\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{b}}}_{R} (t) = ' + latex(expr_db_R_tilde_dt_identical_t.collect(eta_R))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Phonon Correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle}{d t} (t) = i \\chi \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle{\\left(t \\right)} + \\left(2 \\Gamma_{L} - 2 \\gamma\\right) \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle{\\left(t \\right)} - i \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle{\\left(t \\right)} \\overline{\\chi}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle}{d t} (t) = - i \\chi \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle{\\left(t \\right)} + i \\chi \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle{\\left(t \\right)} + \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma - i \\delta\\right) \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle}{d t} (t) = \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma + i \\delta\\right) \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle{\\left(t \\right)} + i \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle{\\left(t \\right)} \\overline{\\chi} - i \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle{\\left(t \\right)} \\overline{\\chi}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle}{d t} (t) = - i \\chi \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle{\\left(t \\right)} + \\left(2 \\Gamma_{R} - 2 \\gamma\\right) \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle{\\left(t \\right)} + i \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle{\\left(t \\right)} \\overline{\\chi}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# substitution list\n",
    "list_subs = [\n",
    "    (a_L_in_tilde_t, 0),\n",
    "    (b_L_in_tilde_t, 0),\n",
    "    (a_R_in_tilde_t, 0),\n",
    "    (b_R_in_tilde_t, 0),\n",
    "    (Dagger(b_L_tilde_t) * b_L_tilde_t, b_L_tilde_dagger_b_L_tilde_expect_t),\n",
    "    (Dagger(b_L_tilde_t) * b_R_tilde_t, b_L_tilde_dagger_b_R_tilde_expect_t),\n",
    "    (Dagger(b_R_tilde_t) * b_L_tilde_t, b_R_tilde_dagger_b_L_tilde_expect_t),\n",
    "    (Dagger(b_R_tilde_t) * b_R_tilde_t, b_R_tilde_dagger_b_R_tilde_expect_t)\n",
    "]\n",
    "\n",
    "# intermode correlations\n",
    "expr_db_L_tilde_dagger_b_L_tilde_expect_dt_t = (Dagger(expr_db_L_tilde_dt_identical_t).doit() * b_L_tilde_t + Dagger(b_L_tilde_t) * expr_db_L_tilde_dt_identical_t).expand().subs(list_subs)\n",
    "display(Math('\\\\frac{d \\\\langle \\\\hat{\\\\tilde{b}}_{L}^{\\\\dagger} \\\\hat{\\\\tilde{b}}_{L} \\\\rangle}{d t} (t) = ' + latex(expr_db_L_tilde_dagger_b_L_tilde_expect_dt_t.collect(b_L_tilde_dagger_b_L_tilde_expect_t))))\n",
    "expr_db_L_tilde_dagger_b_R_tilde_expect_dt_t = (Dagger(expr_db_L_tilde_dt_identical_t).doit() * b_R_tilde_t + Dagger(b_L_tilde_t) * expr_db_R_tilde_dt_identical_t).expand().subs(list_subs)\n",
    "display(Math('\\\\frac{d \\\\langle \\\\hat{\\\\tilde{b}}_{L}^{\\\\dagger} \\\\hat{\\\\tilde{b}}_{R} \\\\rangle}{d t} (t) = ' + latex(expr_db_L_tilde_dagger_b_R_tilde_expect_dt_t.collect(b_L_tilde_dagger_b_R_tilde_expect_t))))\n",
    "expr_db_R_tilde_dagger_b_L_tilde_expect_dt_t = (Dagger(expr_db_R_tilde_dt_identical_t).doit() * b_L_tilde_t + Dagger(b_R_tilde_t) * expr_db_L_tilde_dt_identical_t).expand().subs(list_subs)\n",
    "display(Math('\\\\frac{d \\\\langle \\\\hat{\\\\tilde{b}}_{R}^{\\\\dagger} \\\\hat{\\\\tilde{b}}_{L} \\\\rangle}{d t} (t) = ' + latex(expr_db_R_tilde_dagger_b_L_tilde_expect_dt_t.collect(b_R_tilde_dagger_b_L_tilde_expect_t))))\n",
    "expr_db_R_tilde_dagger_b_R_tilde_expect_dt_t = (Dagger(expr_db_R_tilde_dt_identical_t).doit() * b_R_tilde_t + Dagger(b_R_tilde_t) * expr_db_R_tilde_dt_identical_t).expand().subs(list_subs)\n",
    "display(Math('\\\\frac{d \\\\langle \\\\hat{\\\\tilde{b}}_{R}^{\\\\dagger} \\\\hat{\\\\tilde{b}}_{R} \\\\rangle}{d t} (t) = ' + latex(expr_db_R_tilde_dagger_b_R_tilde_expect_dt_t.collect(b_R_tilde_dagger_b_R_tilde_expect_t))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Difference in Phonon Numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d(\\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle + \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle)}{d t} = \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle \\left(2 \\Gamma_{L} - 2 \\gamma\\right) + \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle \\left(2 \\Gamma_{R} - 2 \\gamma\\right)$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d(\\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle - \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle)}{d t} = \\frac{2 i G_{L} \\lambda \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle \\overline{G_{R}}}{\\kappa^{2} + \\lambda^{2}} - \\frac{2 i G_{R} \\lambda \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle \\overline{G_{L}}}{\\kappa^{2} + \\lambda^{2}} + \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle \\left(- 2 \\Gamma_{L} + 2 \\gamma\\right) + \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle \\left(2 \\Gamma_{R} - 2 \\gamma\\right)$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d\\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle}{d t} = 0 \\Rightarrow \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle = - \\frac{G_{R} \\lambda \\left(\\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle - \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle\\right) \\overline{G_{L}}}{\\left(\\kappa^{2} + \\lambda^{2}\\right) \\left(\\delta + i \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)\\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d\\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle}{d t} = 0 \\Rightarrow \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle = - \\frac{G_{L} \\lambda \\left(\\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle - \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle\\right) \\overline{G_{R}}}{\\left(\\kappa^{2} + \\lambda^{2}\\right) \\left(\\delta - i \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)\\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{d\\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle}{d t} = 0 ~ \\mathrm{and} ~ \\frac{d\\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle}{d t} = 0 \\Rightarrow \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle - \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle = - \\frac{i \\lambda \\left(G_{L} \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle \\overline{G_{R}} - G_{R} \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle \\overline{G_{L}}\\right) \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)}{2 \\left(\\Gamma_{L} - \\gamma\\right) \\left(\\Gamma_{R} - \\gamma\\right) \\left(\\kappa^{2} + \\lambda^{2}\\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\Rightarrow \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle - \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle = - \\frac{\\Gamma_{L} \\Gamma_{R} \\lambda^{2} \\left(\\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle - \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle\\right) \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2}}{\\left(\\Gamma_{L} - \\gamma\\right) \\left(\\Gamma_{R} - \\gamma\\right) \\left(i \\Gamma_{L} + i \\Gamma_{R} - 2 i \\gamma - \\delta\\right) \\left(i \\Gamma_{L} + i \\Gamma_{R} - 2 i \\gamma + \\delta\\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# substitution list\n",
    "list_subs = [\n",
    "    (chi, expr_chi),\n",
    "    (gamma_L, gamma),\n",
    "    (gamma_R, gamma),\n",
    "    (kappa_L, kappa),\n",
    "    (kappa_R, kappa),\n",
    "    (b_L_tilde_dagger_b_L_tilde_expect_t, b_L_tilde_dagger_b_L_tilde_expect),\n",
    "    (b_L_tilde_dagger_b_R_tilde_expect_t, b_L_tilde_dagger_b_R_tilde_expect),\n",
    "    (b_R_tilde_dagger_b_L_tilde_expect_t, b_R_tilde_dagger_b_L_tilde_expect),\n",
    "    (b_R_tilde_dagger_b_R_tilde_expect_t, b_R_tilde_dagger_b_R_tilde_expect)\n",
    "]\n",
    "list_coll = [\n",
    "    b_R_tilde_dagger_b_R_tilde_expect,\n",
    "    - b_L_tilde_dagger_b_L_tilde_expect,\n",
    "    2 * I * chi_R,\n",
    "    2 * chi_I\n",
    "]\n",
    "display(Math('\\\\frac{d(' + latex(b_R_tilde_dagger_b_R_tilde_expect) + ' + ' + latex(b_L_tilde_dagger_b_L_tilde_expect) + ')}{d t} = ' + latex((expr_db_R_tilde_dagger_b_R_tilde_expect_dt_t + expr_db_L_tilde_dagger_b_L_tilde_expect_dt_t).subs(list_subs).expand().collect(list_coll))))\n",
    "display(Math('\\\\frac{d(' + latex(b_R_tilde_dagger_b_R_tilde_expect) + ' - ' + latex(b_L_tilde_dagger_b_L_tilde_expect) + ')}{d t} = ' + latex((expr_db_R_tilde_dagger_b_R_tilde_expect_dt_t - expr_db_L_tilde_dagger_b_L_tilde_expect_dt_t).subs(list_subs).expand().collect(list_coll))))\n",
    "\n",
    "# obtain solutions for phonon correlation\n",
    "sols = solve([\n",
    "    expr_db_L_tilde_dagger_b_R_tilde_expect_dt_t.subs(list_subs),\n",
    "    expr_db_R_tilde_dagger_b_L_tilde_expect_dt_t.subs(list_subs)\n",
    "], [\n",
    "    b_L_tilde_dagger_b_R_tilde_expect,\n",
    "    b_R_tilde_dagger_b_L_tilde_expect\n",
    "])\n",
    "expr_b_L_tilde_dagger_b_R_tilde_expect = sols[b_L_tilde_dagger_b_R_tilde_expect].expand().factor().subs(- Gamma_L - Gamma_R + 2 * gamma + I * delta, I * (delta + I * (Gamma_L + Gamma_R - 2 * gamma)))\n",
    "display(Math('\\\\frac{d' + latex(b_L_tilde_dagger_b_R_tilde_expect) + '}{d t} = 0 \\\\Rightarrow ' + latex(b_L_tilde_dagger_b_R_tilde_expect) + ' = ' + latex(expr_b_L_tilde_dagger_b_R_tilde_expect)))\n",
    "expr_b_R_tilde_dagger_b_L_tilde_expect = sols[b_R_tilde_dagger_b_L_tilde_expect].expand().factor().subs(Gamma_L + Gamma_R - 2 * gamma + I * delta, I * (delta - I * (Gamma_L + Gamma_R - 2 * gamma)))\n",
    "display(Math('\\\\frac{d' + latex(b_R_tilde_dagger_b_L_tilde_expect) + '}{d t} = 0 \\\\Rightarrow ' + latex(b_R_tilde_dagger_b_L_tilde_expect) + ' = ' + latex(expr_b_R_tilde_dagger_b_L_tilde_expect)))\n",
    "\n",
    "# obtain solutions for phonon correlation\n",
    "sols = solve([\n",
    "    expr_db_L_tilde_dagger_b_L_tilde_expect_dt_t.subs(list_subs),\n",
    "    expr_db_R_tilde_dagger_b_R_tilde_expect_dt_t.subs(list_subs)\n",
    "], [\n",
    "    b_L_tilde_dagger_b_L_tilde_expect,\n",
    "    b_R_tilde_dagger_b_R_tilde_expect\n",
    "])\n",
    "expr_n_b_diff = (sols[b_R_tilde_dagger_b_R_tilde_expect] - sols[b_L_tilde_dagger_b_L_tilde_expect]).expand().factor()\n",
    "display(Math('\\\\frac{d' + latex(b_L_tilde_dagger_b_L_tilde_expect) + '}{d t} = 0 ~ \\\\mathrm{and} ~ \\\\frac{d' + latex(b_R_tilde_dagger_b_R_tilde_expect) + '}{d t} = 0 \\\\Rightarrow ' + latex(b_R_tilde_dagger_b_R_tilde_expect) + ' - ' + latex(b_L_tilde_dagger_b_L_tilde_expect) + ' = ' + latex(expr_n_b_diff)))\n",
    "expr_n_b_diff_subs = expr_n_b_diff.subs([\n",
    "    (b_L_tilde_dagger_b_R_tilde_expect, expr_b_L_tilde_dagger_b_R_tilde_expect),\n",
    "    (b_R_tilde_dagger_b_L_tilde_expect, expr_b_R_tilde_dagger_b_L_tilde_expect),\n",
    "    (conjugate(G_L) * G_L, Gamma_L * (kappa**2 + lamb**2)),\n",
    "    (conjugate(G_R) * G_R, Gamma_R * (kappa**2 + lamb**2))\n",
    "]).expand().factor()\n",
    "display(Math('\\\\Rightarrow ' + latex(b_R_tilde_dagger_b_R_tilde_expect) + ' - ' + latex(b_L_tilde_dagger_b_L_tilde_expect) + ' = ' + latex(expr_n_b_diff_subs)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rotated Mechanical Quadrature Rates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\hat{\\tilde{b}}_{L}{\\left(t \\right)} = \\frac{\\sqrt{2} \\left(i \\hat{\\tilde{P}}_{L}{\\left(t \\right)} + \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}\\right)}{2}, \\quad \\hat{\\tilde{b}}_{R}{\\left(t \\right)} = \\frac{\\sqrt{2} \\left(i \\hat{\\tilde{P}}_{R}{\\left(t \\right)} + \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}\\right)}{2}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)} = \\frac{\\sqrt{2} \\left(\\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)} + i \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)}\\right)}{2}, \\quad \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)} = \\frac{\\sqrt{2} \\left(\\hat{\\tilde{X}}_{R}^{in}{\\left(t \\right)} + i \\hat{\\tilde{Y}}_{R}^{in}{\\left(t \\right)}\\right)}{2}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\hat{\\tilde{b}}_{L}^{in}{\\left(t \\right)} = \\frac{\\sqrt{2} \\left(i \\hat{\\tilde{P}}_{L}^{in}{\\left(t \\right)} + \\hat{\\tilde{Q}}_{L}^{in}{\\left(t \\right)}\\right)}{2}, \\quad \\hat{\\tilde{b}}_{R}^{in}{\\left(t \\right)} = \\frac{\\sqrt{2} \\left(i \\hat{\\tilde{P}}_{R}^{in}{\\left(t \\right)} + \\hat{\\tilde{Q}}_{R}^{in}{\\left(t \\right)}\\right)}{2}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{Q}}}_{L} (t) = \\Gamma_{L} \\hat{\\tilde{Q}}_{L}{\\left(t \\right)} - \\chi_{I} \\hat{\\tilde{Q}}_{R}{\\left(t \\right)} + \\chi_{R} \\hat{\\tilde{P}}_{R}{\\left(t \\right)} - \\eta_{L_{I}} \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)} + \\frac{\\eta_{L_{I}} \\lambda \\hat{\\tilde{Y}}_{R}^{in}{\\left(t \\right)}}{\\kappa} + \\eta_{L_{R}} \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)} + \\frac{\\eta_{L_{R}} \\lambda \\hat{\\tilde{X}}_{R}^{in}{\\left(t \\right)}}{\\kappa} + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{Q}}_{L}^{in}{\\left(t \\right)} - \\gamma \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{P}}}_{L} (t) = \\Gamma_{L} \\hat{\\tilde{P}}_{L}{\\left(t \\right)} - \\chi_{I} \\hat{\\tilde{P}}_{R}{\\left(t \\right)} - \\chi_{R} \\hat{\\tilde{Q}}_{R}{\\left(t \\right)} + \\eta_{L_{I}} \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)} + \\frac{\\eta_{L_{I}} \\lambda \\hat{\\tilde{X}}_{R}^{in}{\\left(t \\right)}}{\\kappa} + \\eta_{L_{R}} \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)} - \\frac{\\eta_{L_{R}} \\lambda \\hat{\\tilde{Y}}_{R}^{in}{\\left(t \\right)}}{\\kappa} + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{P}}_{L}^{in}{\\left(t \\right)} - \\gamma \\hat{\\tilde{P}}_{L}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{Q}}}_{R} (t) = \\Gamma_{R} \\hat{\\tilde{Q}}_{R}{\\left(t \\right)} + \\chi_{I} \\hat{\\tilde{Q}}_{L}{\\left(t \\right)} + \\chi_{R} \\hat{\\tilde{P}}_{L}{\\left(t \\right)} - \\eta_{R_{I}} \\hat{\\tilde{X}}_{R}^{in}{\\left(t \\right)} + \\frac{\\eta_{R_{I}} \\lambda \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)}}{\\kappa} + \\eta_{R_{R}} \\hat{\\tilde{Y}}_{R}^{in}{\\left(t \\right)} + \\frac{\\eta_{R_{R}} \\lambda \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)}}{\\kappa} + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{Q}}_{R}^{in}{\\left(t \\right)} - \\gamma \\hat{\\tilde{Q}}_{R}{\\left(t \\right)} + \\delta \\hat{\\tilde{P}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{P}}}_{R} (t) = \\Gamma_{R} \\hat{\\tilde{P}}_{R}{\\left(t \\right)} + \\chi_{I} \\hat{\\tilde{P}}_{L}{\\left(t \\right)} - \\chi_{R} \\hat{\\tilde{Q}}_{L}{\\left(t \\right)} + \\eta_{R_{I}} \\hat{\\tilde{Y}}_{R}^{in}{\\left(t \\right)} + \\frac{\\eta_{R_{I}} \\lambda \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)}}{\\kappa} + \\eta_{R_{R}} \\hat{\\tilde{X}}_{R}^{in}{\\left(t \\right)} - \\frac{\\eta_{R_{R}} \\lambda \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)}}{\\kappa} + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{P}}_{R}^{in}{\\left(t \\right)} - \\gamma \\hat{\\tilde{P}}_{R}{\\left(t \\right)} - \\delta \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# rotated quadratures\n",
    "expr_b_L_tilde_t = (Q_L_tilde_t + I * P_L_tilde_t) / sqrt(2)\n",
    "expr_b_R_tilde_t = (Q_R_tilde_t + I * P_R_tilde_t) / sqrt(2)\n",
    "expr_a_L_in_tilde_t = (X_L_in_tilde_t + I * Y_L_in_tilde_t) / sqrt(2)\n",
    "expr_b_L_in_tilde_t = (Q_L_in_tilde_t + I * P_L_in_tilde_t) / sqrt(2)\n",
    "expr_a_R_in_tilde_t = (X_R_in_tilde_t + I * Y_R_in_tilde_t) / sqrt(2)\n",
    "expr_b_R_in_tilde_t = (Q_R_in_tilde_t + I * P_R_in_tilde_t) / sqrt(2)\n",
    "# substitution list\n",
    "list_subs = [\n",
    "    (b_L_tilde_t, expr_b_L_tilde_t),\n",
    "    (Dagger(Q_L_tilde_t), Q_L_tilde_t),\n",
    "    (Dagger(P_L_tilde_t), P_L_tilde_t),\n",
    "    (b_R_tilde_t, expr_b_R_tilde_t),\n",
    "    (Dagger(Q_R_tilde_t), Q_R_tilde_t),\n",
    "    (Dagger(P_R_tilde_t), P_R_tilde_t),\n",
    "    (a_L_in_tilde_t, expr_a_L_in_tilde_t),\n",
    "    (Dagger(X_L_in_tilde_t), X_L_in_tilde_t),\n",
    "    (Dagger(Y_L_in_tilde_t), Y_L_in_tilde_t),\n",
    "    (b_L_in_tilde_t, expr_b_L_in_tilde_t),\n",
    "    (Dagger(Q_L_in_tilde_t), Q_L_in_tilde_t),\n",
    "    (Dagger(P_L_in_tilde_t), P_L_in_tilde_t),\n",
    "    (a_R_in_tilde_t, expr_a_R_in_tilde_t),\n",
    "    (Dagger(X_R_in_tilde_t), X_R_in_tilde_t),\n",
    "    (Dagger(Y_R_in_tilde_t), Y_R_in_tilde_t),\n",
    "    (b_R_in_tilde_t, expr_b_R_in_tilde_t),\n",
    "    (Dagger(Q_R_in_tilde_t), Q_R_in_tilde_t),\n",
    "    (Dagger(P_R_in_tilde_t), P_R_in_tilde_t),\n",
    "    (chi, chi_R + I * chi_I),\n",
    "    (eta_L, eta_LR + I * eta_LI),\n",
    "    (eta_R, eta_RR + I * eta_RI)\n",
    "]\n",
    "# display\n",
    "display(Math(latex(b_L_tilde_t) + ' = ' + latex(expr_b_L_tilde_t) + ', \\\\quad ' + latex(b_R_tilde_t) + ' = ' + latex(expr_b_R_tilde_t)))\n",
    "display(Math(latex(a_L_in_tilde_t) + ' = ' + latex(expr_a_L_in_tilde_t) + ', \\\\quad ' + latex(a_R_in_tilde_t) + ' = ' + latex(expr_a_R_in_tilde_t)))\n",
    "display(Math(latex(b_L_in_tilde_t) + ' = ' + latex(expr_b_L_in_tilde_t) + ', \\\\quad ' + latex(b_R_in_tilde_t) + ' = ' + latex(expr_b_R_in_tilde_t)))\n",
    "\n",
    "# express and substitute\n",
    "expr_dQ_L_tilde_dt_t = (Dagger(expr_db_L_tilde_dt_identical_t) + expr_db_L_tilde_dt_identical_t) / sqrt(2)\n",
    "expr_dQ_L_tilde_dt_t = expr_dQ_L_tilde_dt_t.expand().subs(list_subs).expand()\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{Q}}}_{L} (t) = ' + latex(expr_dQ_L_tilde_dt_t)))\n",
    "expr_dP_L_tilde_dt_t = I * (Dagger(expr_db_L_tilde_dt_identical_t) - expr_db_L_tilde_dt_identical_t) / sqrt(2)\n",
    "expr_dP_L_tilde_dt_t = expr_dP_L_tilde_dt_t.expand().subs(list_subs).expand()\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{P}}}_{L} (t) = ' + latex(expr_dP_L_tilde_dt_t)))\n",
    "expr_dQ_R_tilde_dt_t = (Dagger(expr_db_R_tilde_dt_identical_t) + expr_db_R_tilde_dt_identical_t) / sqrt(2)\n",
    "expr_dQ_R_tilde_dt_t = expr_dQ_R_tilde_dt_t.expand().subs(list_subs).expand()\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{Q}}}_{R} (t) = ' + latex(expr_dQ_R_tilde_dt_t)))\n",
    "expr_dP_R_tilde_dt_t = I * (Dagger(expr_db_R_tilde_dt_identical_t) - expr_db_R_tilde_dt_identical_t) / sqrt(2)\n",
    "expr_dP_R_tilde_dt_t = expr_dP_R_tilde_dt_t.expand().subs(list_subs).expand()\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{P}}}_{R} (t) = ' + latex(expr_dP_R_tilde_dt_t)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle G_{L} = \\alpha_{L} g_{L}, \\quad G_{R} = \\alpha_{R} g_{R}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\Delta_{L} = \\omega_{L}, \\quad \\Delta_{R} = \\omega_{R}, \\quad \\delta = - \\omega_{L} + \\omega_{R}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\Gamma_{L} = \\frac{G_{L} \\kappa_{R} \\overline{G_{L}}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}}, \\quad \\Gamma_{R} = \\frac{G_{R} \\kappa_{L} \\overline{G_{R}}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\chi = \\frac{G_{R} \\lambda \\overline{G_{L}}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} = i \\chi_{I} + \\chi_{R}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\eta_{L} = \\frac{\\sqrt{2} G_{L} \\sqrt{\\kappa_{L}} \\kappa_{R}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} = i \\eta_{L_{I}} + \\eta_{L_{R}}, \\quad \\eta_{R} = \\frac{\\sqrt{2} G_{R} \\kappa_{L} \\sqrt{\\kappa_{R}}}{\\kappa_{L} \\kappa_{R} + \\lambda^{2}} = i \\eta_{R_{I}} + \\eta_{R_{R}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\gamma_{L} = \\gamma_{R} = \\gamma, \\quad \\kappa_{L} = \\kappa_{R} = \\kappa$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{b}}}_{L} (t) = \\Gamma_{L} \\hat{\\tilde{b}}_{L}{\\left(t \\right)} + i \\eta_{L} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger} + \\frac{\\eta_{L} \\lambda \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}^{\\dagger}}{\\kappa} + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{b}}_{L}^{in}{\\left(t \\right)} - \\gamma \\hat{\\tilde{b}}_{L}{\\left(t \\right)} - i \\overline{\\chi} \\hat{\\tilde{b}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{b}}}_{R} (t) = \\Gamma_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)} - i \\chi \\hat{\\tilde{b}}_{L}{\\left(t \\right)} + i \\eta_{R} \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}^{\\dagger} + \\frac{\\eta_{R} \\lambda \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger}}{\\kappa} + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{b}}_{R}^{in}{\\left(t \\right)} - \\gamma \\hat{\\tilde{b}}_{R}{\\left(t \\right)} - i \\delta \\hat{\\tilde{b}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{Q}}}_{L} (t) = \\Gamma_{L} \\hat{\\tilde{Q}}_{L}{\\left(t \\right)} - \\chi_{I} \\hat{\\tilde{Q}}_{R}{\\left(t \\right)} + \\chi_{R} \\hat{\\tilde{P}}_{R}{\\left(t \\right)} - \\eta_{L_{I}} \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)} + \\frac{\\eta_{L_{I}} \\lambda \\hat{\\tilde{Y}}_{R}^{in}{\\left(t \\right)}}{\\kappa} + \\eta_{L_{R}} \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)} + \\frac{\\eta_{L_{R}} \\lambda \\hat{\\tilde{X}}_{R}^{in}{\\left(t \\right)}}{\\kappa} + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{Q}}_{L}^{in}{\\left(t \\right)} - \\gamma \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{P}}}_{L} (t) = \\Gamma_{L} \\hat{\\tilde{P}}_{L}{\\left(t \\right)} - \\chi_{I} \\hat{\\tilde{P}}_{R}{\\left(t \\right)} - \\chi_{R} \\hat{\\tilde{Q}}_{R}{\\left(t \\right)} + \\eta_{L_{I}} \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)} + \\frac{\\eta_{L_{I}} \\lambda \\hat{\\tilde{X}}_{R}^{in}{\\left(t \\right)}}{\\kappa} + \\eta_{L_{R}} \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)} - \\frac{\\eta_{L_{R}} \\lambda \\hat{\\tilde{Y}}_{R}^{in}{\\left(t \\right)}}{\\kappa} + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{P}}_{L}^{in}{\\left(t \\right)} - \\gamma \\hat{\\tilde{P}}_{L}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{Q}}}_{R} (t) = \\Gamma_{R} \\hat{\\tilde{Q}}_{R}{\\left(t \\right)} + \\chi_{I} \\hat{\\tilde{Q}}_{L}{\\left(t \\right)} + \\chi_{R} \\hat{\\tilde{P}}_{L}{\\left(t \\right)} - \\eta_{R_{I}} \\hat{\\tilde{X}}_{R}^{in}{\\left(t \\right)} + \\frac{\\eta_{R_{I}} \\lambda \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)}}{\\kappa} + \\eta_{R_{R}} \\hat{\\tilde{Y}}_{R}^{in}{\\left(t \\right)} + \\frac{\\eta_{R_{R}} \\lambda \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)}}{\\kappa} + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{Q}}_{R}^{in}{\\left(t \\right)} - \\gamma \\hat{\\tilde{Q}}_{R}{\\left(t \\right)} + \\delta \\hat{\\tilde{P}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{P}}}_{R} (t) = \\Gamma_{R} \\hat{\\tilde{P}}_{R}{\\left(t \\right)} + \\chi_{I} \\hat{\\tilde{P}}_{L}{\\left(t \\right)} - \\chi_{R} \\hat{\\tilde{Q}}_{L}{\\left(t \\right)} + \\eta_{R_{I}} \\hat{\\tilde{Y}}_{R}^{in}{\\left(t \\right)} + \\frac{\\eta_{R_{I}} \\lambda \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)}}{\\kappa} + \\eta_{R_{R}} \\hat{\\tilde{X}}_{R}^{in}{\\left(t \\right)} - \\frac{\\eta_{R_{R}} \\lambda \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)}}{\\kappa} + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{P}}_{R}^{in}{\\left(t \\right)} - \\gamma \\hat{\\tilde{P}}_{R}{\\left(t \\right)} - \\delta \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# final expressions\n",
    "display(Math(latex(G_L) + ' = ' + latex(expr_G_L) + ', \\\\quad ' + latex(G_R) + ' = ' + latex(expr_G_R)))\n",
    "display(Math(latex(Delta_L) + ' = ' + latex(omega_L) + ', \\\\quad ' + latex(Delta_R) + ' = ' + latex(omega_R) + ', \\\\quad ' + latex(delta) + ' = ' + latex(expr_delta)))\n",
    "display(Math(latex(Gamma_L) + ' = ' + latex(expr_Gamma_L) + ', \\\\quad ' + latex(Gamma_R) + ' = ' + latex(expr_Gamma_R)))\n",
    "display(Math(latex(chi) + ' = ' + latex(expr_chi) + ' = ' + latex(chi_R + I * chi_I)))\n",
    "display(Math(latex(eta_L) + ' = ' + latex(expr_eta_L) + ' = ' + latex(eta_LR + I * eta_LI) + ', \\\\quad ' + latex(eta_R) + ' = ' + latex(expr_eta_R) + ' = ' + latex(eta_RR + I * eta_RI)))\n",
    "display(Math(latex(gamma_L) + ' = ' + latex(gamma_R) + ' = ' + latex(gamma) + ', \\\\quad ' + latex(kappa_L) + ' = ' + latex(kappa_R) + ' = ' + latex(kappa)))\n",
    "# final mode rates\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{b}}}_{L} (t) = ' + latex(expr_db_L_tilde_dt_identical_t)))\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{b}}}_{R} (t) = ' + latex(expr_db_R_tilde_dt_identical_t)))\n",
    "# final quadrature rates\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{Q}}}_{L} (t) = ' + latex(expr_dQ_L_tilde_dt_t)))\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{P}}}_{L} (t) = ' + latex(expr_dP_L_tilde_dt_t)))\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{Q}}}_{R} (t) = ' + latex(expr_dQ_R_tilde_dt_t)))\n",
    "display(Math('\\\\dot{\\\\hat{\\\\tilde{P}}}_{R} (t) = ' + latex(expr_dP_R_tilde_dt_t)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle A = \\left[\\begin{matrix}\\Gamma_{L} - \\gamma & 0 & - \\chi_{I} & \\chi_{R}\\\\0 & \\Gamma_{L} - \\gamma & - \\chi_{R} & - \\chi_{I}\\\\\\chi_{I} & \\chi_{R} & \\Gamma_{R} - \\gamma & \\delta\\\\- \\chi_{R} & \\chi_{I} & - \\delta & \\Gamma_{R} - \\gamma\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle D = \\left[\\begin{matrix}d_{L} & 0 & 0 & 0\\\\0 & d_{L} & 0 & 0\\\\0 & 0 & d_{R} & 0\\\\0 & 0 & 0 & d_{R}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle V = \\left[\\begin{matrix}\\langle \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}\\hat{\\tilde{Q}}_{L}{\\left(t \\right)} \\rangle & \\frac{\\langle \\hat{\\tilde{P}}_{L}{\\left(t \\right)}\\hat{\\tilde{Q}}_{L}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}\\hat{\\tilde{P}}_{L}{\\left(t \\right)} \\rangle}{2} & \\frac{\\langle \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}\\hat{\\tilde{Q}}_{R}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}\\hat{\\tilde{Q}}_{L}{\\left(t \\right)} \\rangle}{2} & \\frac{\\langle \\hat{\\tilde{P}}_{R}{\\left(t \\right)}\\hat{\\tilde{Q}}_{L}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}\\hat{\\tilde{P}}_{R}{\\left(t \\right)} \\rangle}{2}\\\\\\frac{\\langle \\hat{\\tilde{P}}_{L}{\\left(t \\right)}\\hat{\\tilde{Q}}_{L}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}\\hat{\\tilde{P}}_{L}{\\left(t \\right)} \\rangle}{2} & \\langle \\hat{\\tilde{P}}_{L}{\\left(t \\right)}\\hat{\\tilde{P}}_{L}{\\left(t \\right)} \\rangle & \\frac{\\langle \\hat{\\tilde{P}}_{L}{\\left(t \\right)}\\hat{\\tilde{Q}}_{R}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}\\hat{\\tilde{P}}_{L}{\\left(t \\right)} \\rangle}{2} & \\frac{\\langle \\hat{\\tilde{P}}_{L}{\\left(t \\right)}\\hat{\\tilde{P}}_{R}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{P}}_{R}{\\left(t \\right)}\\hat{\\tilde{P}}_{L}{\\left(t \\right)} \\rangle}{2}\\\\\\frac{\\langle \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}\\hat{\\tilde{Q}}_{R}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}\\hat{\\tilde{Q}}_{L}{\\left(t \\right)} \\rangle}{2} & \\frac{\\langle \\hat{\\tilde{P}}_{L}{\\left(t \\right)}\\hat{\\tilde{Q}}_{R}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}\\hat{\\tilde{P}}_{L}{\\left(t \\right)} \\rangle}{2} & \\langle \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}\\hat{\\tilde{Q}}_{R}{\\left(t \\right)} \\rangle & \\frac{\\langle \\hat{\\tilde{P}}_{R}{\\left(t \\right)}\\hat{\\tilde{Q}}_{R}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}\\hat{\\tilde{P}}_{R}{\\left(t \\right)} \\rangle}{2}\\\\\\frac{\\langle \\hat{\\tilde{P}}_{R}{\\left(t \\right)}\\hat{\\tilde{Q}}_{L}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}\\hat{\\tilde{P}}_{R}{\\left(t \\right)} \\rangle}{2} & \\frac{\\langle \\hat{\\tilde{P}}_{L}{\\left(t \\right)}\\hat{\\tilde{P}}_{R}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{P}}_{R}{\\left(t \\right)}\\hat{\\tilde{P}}_{L}{\\left(t \\right)} \\rangle}{2} & \\frac{\\langle \\hat{\\tilde{P}}_{R}{\\left(t \\right)}\\hat{\\tilde{Q}}_{R}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}\\hat{\\tilde{P}}_{R}{\\left(t \\right)} \\rangle}{2} & \\langle \\hat{\\tilde{P}}_{R}{\\left(t \\right)}\\hat{\\tilde{P}}_{R}{\\left(t \\right)} \\rangle\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# drift matrix\n",
    "A = Matrix([\n",
    "    [(Gamma_L - gamma), 0, - chi_I, chi_R], \n",
    "    [0, (Gamma_L - gamma), - chi_R, - chi_I], \n",
    "    [chi_I, chi_R, (Gamma_R - gamma), delta], \n",
    "    [- chi_R, chi_I, - delta, (Gamma_R - gamma)]\n",
    "])\n",
    "display(Math('A = ' + latex(A)))\n",
    "\n",
    "# noise matrix\n",
    "d_L, d_R = symbols('d_{L}, d_{R}', real=True)\n",
    "diag_L = conjugate(eta_L) * eta_L * (1 + lamb**2 / kappa**2) / 2 + gamma * (2 * n_th + 1)\n",
    "diag_R = conjugate(eta_R) * eta_R * (1 + lamb**2 / kappa**2) / 2 + gamma * (2 * n_th + 1)\n",
    "D = Matrix([[d_L, 0, 0, 0], [0, d_L, 0, 0], [0, 0, d_R, 0], [0, 0, 0, d_R]])\n",
    "display(Math('D = ' + latex(D)))\n",
    "\n",
    "# correlation matrix\n",
    "u = [Q_L_tilde_t, P_L_tilde_t, Q_R_tilde_t, P_R_tilde_t]\n",
    "correlations_expect = Matrix([[Symbol('\\\\langle ' + latex(u_i) + latex(u_j) + ' \\\\rangle', real=True) for u_i in u] for u_j in u])\n",
    "V = (correlations_expect + correlations_expect.T) / 2\n",
    "display(Math('V = ' + latex(V)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\langle \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}\\hat{\\tilde{Q}}_{L}{\\left(t \\right)} \\rangle = \\frac{- \\frac{\\Gamma_{L}^{2} \\Gamma_{R} d_{L}}{2} + \\frac{\\Gamma_{L}^{2} \\gamma d_{L}}{2} - \\Gamma_{L} \\Gamma_{R}^{2} d_{L} + 3 \\Gamma_{L} \\Gamma_{R} \\gamma d_{L} - \\frac{\\Gamma_{L} \\chi_{I}^{2} d_{L}}{2} - \\frac{\\Gamma_{L} \\chi_{I}^{2} d_{R}}{2} - \\frac{\\Gamma_{L} \\chi_{R}^{2} d_{L}}{2} - \\frac{\\Gamma_{L} \\chi_{R}^{2} d_{R}}{2} - 2 \\Gamma_{L} \\gamma^{2} d_{L} - \\frac{\\Gamma_{R}^{3} d_{L}}{2} + \\frac{5 \\Gamma_{R}^{2} \\gamma d_{L}}{2} - \\frac{\\Gamma_{R} \\chi_{I}^{2} d_{L}}{2} - \\frac{\\Gamma_{R} \\chi_{I}^{2} d_{R}}{2} - \\frac{\\Gamma_{R} \\chi_{R}^{2} d_{L}}{2} - \\frac{\\Gamma_{R} \\chi_{R}^{2} d_{R}}{2} - 4 \\Gamma_{R} \\gamma^{2} d_{L} - \\frac{\\Gamma_{R} d_{L} \\delta^{2}}{2} + \\chi_{I}^{2} \\gamma d_{L} + \\chi_{I}^{2} \\gamma d_{R} + \\chi_{R}^{2} \\gamma d_{L} + \\chi_{R}^{2} \\gamma d_{R} + 2 \\gamma^{3} d_{L} + \\frac{\\gamma d_{L} \\delta^{2}}{2}}{\\Gamma_{L}^{3} \\Gamma_{R} - \\Gamma_{L}^{3} \\gamma + 2 \\Gamma_{L}^{2} \\Gamma_{R}^{2} - 7 \\Gamma_{L}^{2} \\Gamma_{R} \\gamma + \\Gamma_{L}^{2} \\chi_{I}^{2} + \\Gamma_{L}^{2} \\chi_{R}^{2} + 5 \\Gamma_{L}^{2} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R}^{3} - 7 \\Gamma_{L} \\Gamma_{R}^{2} \\gamma + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{I}^{2} + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{R}^{2} + 14 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R} \\delta^{2} - 4 \\Gamma_{L} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{L} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{L} \\gamma^{3} - \\Gamma_{L} \\gamma \\delta^{2} - \\Gamma_{R}^{3} \\gamma + \\Gamma_{R}^{2} \\chi_{I}^{2} + \\Gamma_{R}^{2} \\chi_{R}^{2} + 5 \\Gamma_{R}^{2} \\gamma^{2} - 4 \\Gamma_{R} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{R} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{R} \\gamma^{3} - \\Gamma_{R} \\gamma \\delta^{2} + 4 \\chi_{I}^{2} \\gamma^{2} + 4 \\chi_{R}^{2} \\gamma^{2} + 4 \\gamma^{4} + \\gamma^{2} \\delta^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\langle \\hat{\\tilde{P}}_{L}{\\left(t \\right)}\\hat{\\tilde{Q}}_{L}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}\\hat{\\tilde{P}}_{L}{\\left(t \\right)} \\rangle}{2} = 0$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\langle \\hat{\\tilde{P}}_{L}{\\left(t \\right)}\\hat{\\tilde{P}}_{L}{\\left(t \\right)} \\rangle = \\frac{- \\frac{\\Gamma_{L}^{2} \\Gamma_{R} d_{L}}{2} + \\frac{\\Gamma_{L}^{2} \\gamma d_{L}}{2} - \\Gamma_{L} \\Gamma_{R}^{2} d_{L} + 3 \\Gamma_{L} \\Gamma_{R} \\gamma d_{L} - \\frac{\\Gamma_{L} \\chi_{I}^{2} d_{L}}{2} - \\frac{\\Gamma_{L} \\chi_{I}^{2} d_{R}}{2} - \\frac{\\Gamma_{L} \\chi_{R}^{2} d_{L}}{2} - \\frac{\\Gamma_{L} \\chi_{R}^{2} d_{R}}{2} - 2 \\Gamma_{L} \\gamma^{2} d_{L} - \\frac{\\Gamma_{R}^{3} d_{L}}{2} + \\frac{5 \\Gamma_{R}^{2} \\gamma d_{L}}{2} - \\frac{\\Gamma_{R} \\chi_{I}^{2} d_{L}}{2} - \\frac{\\Gamma_{R} \\chi_{I}^{2} d_{R}}{2} - \\frac{\\Gamma_{R} \\chi_{R}^{2} d_{L}}{2} - \\frac{\\Gamma_{R} \\chi_{R}^{2} d_{R}}{2} - 4 \\Gamma_{R} \\gamma^{2} d_{L} - \\frac{\\Gamma_{R} d_{L} \\delta^{2}}{2} + \\chi_{I}^{2} \\gamma d_{L} + \\chi_{I}^{2} \\gamma d_{R} + \\chi_{R}^{2} \\gamma d_{L} + \\chi_{R}^{2} \\gamma d_{R} + 2 \\gamma^{3} d_{L} + \\frac{\\gamma d_{L} \\delta^{2}}{2}}{\\Gamma_{L}^{3} \\Gamma_{R} - \\Gamma_{L}^{3} \\gamma + 2 \\Gamma_{L}^{2} \\Gamma_{R}^{2} - 7 \\Gamma_{L}^{2} \\Gamma_{R} \\gamma + \\Gamma_{L}^{2} \\chi_{I}^{2} + \\Gamma_{L}^{2} \\chi_{R}^{2} + 5 \\Gamma_{L}^{2} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R}^{3} - 7 \\Gamma_{L} \\Gamma_{R}^{2} \\gamma + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{I}^{2} + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{R}^{2} + 14 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R} \\delta^{2} - 4 \\Gamma_{L} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{L} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{L} \\gamma^{3} - \\Gamma_{L} \\gamma \\delta^{2} - \\Gamma_{R}^{3} \\gamma + \\Gamma_{R}^{2} \\chi_{I}^{2} + \\Gamma_{R}^{2} \\chi_{R}^{2} + 5 \\Gamma_{R}^{2} \\gamma^{2} - 4 \\Gamma_{R} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{R} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{R} \\gamma^{3} - \\Gamma_{R} \\gamma \\delta^{2} + 4 \\chi_{I}^{2} \\gamma^{2} + 4 \\chi_{R}^{2} \\gamma^{2} + 4 \\gamma^{4} + \\gamma^{2} \\delta^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\langle \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}\\hat{\\tilde{Q}}_{R}{\\left(t \\right)} \\rangle = \\frac{- \\frac{\\Gamma_{L}^{3} d_{R}}{2} - \\Gamma_{L}^{2} \\Gamma_{R} d_{R} + \\frac{5 \\Gamma_{L}^{2} \\gamma d_{R}}{2} - \\frac{\\Gamma_{L} \\Gamma_{R}^{2} d_{R}}{2} + 3 \\Gamma_{L} \\Gamma_{R} \\gamma d_{R} - \\frac{\\Gamma_{L} \\chi_{I}^{2} d_{L}}{2} - \\frac{\\Gamma_{L} \\chi_{I}^{2} d_{R}}{2} - \\frac{\\Gamma_{L} \\chi_{R}^{2} d_{L}}{2} - \\frac{\\Gamma_{L} \\chi_{R}^{2} d_{R}}{2} - 4 \\Gamma_{L} \\gamma^{2} d_{R} - \\frac{\\Gamma_{L} d_{R} \\delta^{2}}{2} + \\frac{\\Gamma_{R}^{2} \\gamma d_{R}}{2} - \\frac{\\Gamma_{R} \\chi_{I}^{2} d_{L}}{2} - \\frac{\\Gamma_{R} \\chi_{I}^{2} d_{R}}{2} - \\frac{\\Gamma_{R} \\chi_{R}^{2} d_{L}}{2} - \\frac{\\Gamma_{R} \\chi_{R}^{2} d_{R}}{2} - 2 \\Gamma_{R} \\gamma^{2} d_{R} + \\chi_{I}^{2} \\gamma d_{L} + \\chi_{I}^{2} \\gamma d_{R} + \\chi_{R}^{2} \\gamma d_{L} + \\chi_{R}^{2} \\gamma d_{R} + 2 \\gamma^{3} d_{R} + \\frac{\\gamma d_{R} \\delta^{2}}{2}}{\\Gamma_{L}^{3} \\Gamma_{R} - \\Gamma_{L}^{3} \\gamma + 2 \\Gamma_{L}^{2} \\Gamma_{R}^{2} - 7 \\Gamma_{L}^{2} \\Gamma_{R} \\gamma + \\Gamma_{L}^{2} \\chi_{I}^{2} + \\Gamma_{L}^{2} \\chi_{R}^{2} + 5 \\Gamma_{L}^{2} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R}^{3} - 7 \\Gamma_{L} \\Gamma_{R}^{2} \\gamma + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{I}^{2} + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{R}^{2} + 14 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R} \\delta^{2} - 4 \\Gamma_{L} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{L} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{L} \\gamma^{3} - \\Gamma_{L} \\gamma \\delta^{2} - \\Gamma_{R}^{3} \\gamma + \\Gamma_{R}^{2} \\chi_{I}^{2} + \\Gamma_{R}^{2} \\chi_{R}^{2} + 5 \\Gamma_{R}^{2} \\gamma^{2} - 4 \\Gamma_{R} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{R} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{R} \\gamma^{3} - \\Gamma_{R} \\gamma \\delta^{2} + 4 \\chi_{I}^{2} \\gamma^{2} + 4 \\chi_{R}^{2} \\gamma^{2} + 4 \\gamma^{4} + \\gamma^{2} \\delta^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\langle \\hat{\\tilde{P}}_{R}{\\left(t \\right)}\\hat{\\tilde{Q}}_{R}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}\\hat{\\tilde{P}}_{R}{\\left(t \\right)} \\rangle}{2} = 0$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\langle \\hat{\\tilde{P}}_{R}{\\left(t \\right)}\\hat{\\tilde{P}}_{R}{\\left(t \\right)} \\rangle = \\frac{- \\frac{\\Gamma_{L}^{3} d_{R}}{2} - \\Gamma_{L}^{2} \\Gamma_{R} d_{R} + \\frac{5 \\Gamma_{L}^{2} \\gamma d_{R}}{2} - \\frac{\\Gamma_{L} \\Gamma_{R}^{2} d_{R}}{2} + 3 \\Gamma_{L} \\Gamma_{R} \\gamma d_{R} - \\frac{\\Gamma_{L} \\chi_{I}^{2} d_{L}}{2} - \\frac{\\Gamma_{L} \\chi_{I}^{2} d_{R}}{2} - \\frac{\\Gamma_{L} \\chi_{R}^{2} d_{L}}{2} - \\frac{\\Gamma_{L} \\chi_{R}^{2} d_{R}}{2} - 4 \\Gamma_{L} \\gamma^{2} d_{R} - \\frac{\\Gamma_{L} d_{R} \\delta^{2}}{2} + \\frac{\\Gamma_{R}^{2} \\gamma d_{R}}{2} - \\frac{\\Gamma_{R} \\chi_{I}^{2} d_{L}}{2} - \\frac{\\Gamma_{R} \\chi_{I}^{2} d_{R}}{2} - \\frac{\\Gamma_{R} \\chi_{R}^{2} d_{L}}{2} - \\frac{\\Gamma_{R} \\chi_{R}^{2} d_{R}}{2} - 2 \\Gamma_{R} \\gamma^{2} d_{R} + \\chi_{I}^{2} \\gamma d_{L} + \\chi_{I}^{2} \\gamma d_{R} + \\chi_{R}^{2} \\gamma d_{L} + \\chi_{R}^{2} \\gamma d_{R} + 2 \\gamma^{3} d_{R} + \\frac{\\gamma d_{R} \\delta^{2}}{2}}{\\Gamma_{L}^{3} \\Gamma_{R} - \\Gamma_{L}^{3} \\gamma + 2 \\Gamma_{L}^{2} \\Gamma_{R}^{2} - 7 \\Gamma_{L}^{2} \\Gamma_{R} \\gamma + \\Gamma_{L}^{2} \\chi_{I}^{2} + \\Gamma_{L}^{2} \\chi_{R}^{2} + 5 \\Gamma_{L}^{2} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R}^{3} - 7 \\Gamma_{L} \\Gamma_{R}^{2} \\gamma + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{I}^{2} + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{R}^{2} + 14 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R} \\delta^{2} - 4 \\Gamma_{L} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{L} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{L} \\gamma^{3} - \\Gamma_{L} \\gamma \\delta^{2} - \\Gamma_{R}^{3} \\gamma + \\Gamma_{R}^{2} \\chi_{I}^{2} + \\Gamma_{R}^{2} \\chi_{R}^{2} + 5 \\Gamma_{R}^{2} \\gamma^{2} - 4 \\Gamma_{R} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{R} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{R} \\gamma^{3} - \\Gamma_{R} \\gamma \\delta^{2} + 4 \\chi_{I}^{2} \\gamma^{2} + 4 \\chi_{R}^{2} \\gamma^{2} + 4 \\gamma^{4} + \\gamma^{2} \\delta^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\langle \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}\\hat{\\tilde{Q}}_{R}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}\\hat{\\tilde{Q}}_{L}{\\left(t \\right)} \\rangle}{2} = \\frac{- \\frac{\\Gamma_{L}^{2} \\chi_{I} d_{R}}{2} + \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{I} d_{L}}{2} - \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{I} d_{R}}{2} - \\frac{\\Gamma_{L} \\chi_{I} \\gamma d_{L}}{2} + \\frac{3 \\Gamma_{L} \\chi_{I} \\gamma d_{R}}{2} - \\frac{\\Gamma_{L} \\chi_{R} d_{R} \\delta}{2} + \\frac{\\Gamma_{R}^{2} \\chi_{I} d_{L}}{2} - \\frac{3 \\Gamma_{R} \\chi_{I} \\gamma d_{L}}{2} + \\frac{\\Gamma_{R} \\chi_{I} \\gamma d_{R}}{2} + \\frac{\\Gamma_{R} \\chi_{R} d_{L} \\delta}{2} + \\chi_{I} \\gamma^{2} d_{L} - \\chi_{I} \\gamma^{2} d_{R} - \\frac{\\chi_{R} \\gamma d_{L} \\delta}{2} + \\frac{\\chi_{R} \\gamma d_{R} \\delta}{2}}{\\Gamma_{L}^{3} \\Gamma_{R} - \\Gamma_{L}^{3} \\gamma + 2 \\Gamma_{L}^{2} \\Gamma_{R}^{2} - 7 \\Gamma_{L}^{2} \\Gamma_{R} \\gamma + \\Gamma_{L}^{2} \\chi_{I}^{2} + \\Gamma_{L}^{2} \\chi_{R}^{2} + 5 \\Gamma_{L}^{2} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R}^{3} - 7 \\Gamma_{L} \\Gamma_{R}^{2} \\gamma + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{I}^{2} + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{R}^{2} + 14 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R} \\delta^{2} - 4 \\Gamma_{L} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{L} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{L} \\gamma^{3} - \\Gamma_{L} \\gamma \\delta^{2} - \\Gamma_{R}^{3} \\gamma + \\Gamma_{R}^{2} \\chi_{I}^{2} + \\Gamma_{R}^{2} \\chi_{R}^{2} + 5 \\Gamma_{R}^{2} \\gamma^{2} - 4 \\Gamma_{R} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{R} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{R} \\gamma^{3} - \\Gamma_{R} \\gamma \\delta^{2} + 4 \\chi_{I}^{2} \\gamma^{2} + 4 \\chi_{R}^{2} \\gamma^{2} + 4 \\gamma^{4} + \\gamma^{2} \\delta^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\langle \\hat{\\tilde{P}}_{R}{\\left(t \\right)}\\hat{\\tilde{Q}}_{L}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}\\hat{\\tilde{P}}_{R}{\\left(t \\right)} \\rangle}{2} = \\frac{\\frac{\\Gamma_{L}^{2} \\chi_{R} d_{R}}{2} - \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{R} d_{L}}{2} + \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{R} d_{R}}{2} - \\frac{\\Gamma_{L} \\chi_{I} d_{R} \\delta}{2} + \\frac{\\Gamma_{L} \\chi_{R} \\gamma d_{L}}{2} - \\frac{3 \\Gamma_{L} \\chi_{R} \\gamma d_{R}}{2} - \\frac{\\Gamma_{R}^{2} \\chi_{R} d_{L}}{2} + \\frac{\\Gamma_{R} \\chi_{I} d_{L} \\delta}{2} + \\frac{3 \\Gamma_{R} \\chi_{R} \\gamma d_{L}}{2} - \\frac{\\Gamma_{R} \\chi_{R} \\gamma d_{R}}{2} - \\frac{\\chi_{I} \\gamma d_{L} \\delta}{2} + \\frac{\\chi_{I} \\gamma d_{R} \\delta}{2} - \\chi_{R} \\gamma^{2} d_{L} + \\chi_{R} \\gamma^{2} d_{R}}{\\Gamma_{L}^{3} \\Gamma_{R} - \\Gamma_{L}^{3} \\gamma + 2 \\Gamma_{L}^{2} \\Gamma_{R}^{2} - 7 \\Gamma_{L}^{2} \\Gamma_{R} \\gamma + \\Gamma_{L}^{2} \\chi_{I}^{2} + \\Gamma_{L}^{2} \\chi_{R}^{2} + 5 \\Gamma_{L}^{2} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R}^{3} - 7 \\Gamma_{L} \\Gamma_{R}^{2} \\gamma + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{I}^{2} + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{R}^{2} + 14 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R} \\delta^{2} - 4 \\Gamma_{L} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{L} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{L} \\gamma^{3} - \\Gamma_{L} \\gamma \\delta^{2} - \\Gamma_{R}^{3} \\gamma + \\Gamma_{R}^{2} \\chi_{I}^{2} + \\Gamma_{R}^{2} \\chi_{R}^{2} + 5 \\Gamma_{R}^{2} \\gamma^{2} - 4 \\Gamma_{R} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{R} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{R} \\gamma^{3} - \\Gamma_{R} \\gamma \\delta^{2} + 4 \\chi_{I}^{2} \\gamma^{2} + 4 \\chi_{R}^{2} \\gamma^{2} + 4 \\gamma^{4} + \\gamma^{2} \\delta^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\langle \\hat{\\tilde{P}}_{L}{\\left(t \\right)}\\hat{\\tilde{Q}}_{R}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}\\hat{\\tilde{P}}_{L}{\\left(t \\right)} \\rangle}{2} = \\frac{- \\frac{\\Gamma_{L}^{2} \\chi_{R} d_{R}}{2} + \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{R} d_{L}}{2} - \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{R} d_{R}}{2} + \\frac{\\Gamma_{L} \\chi_{I} d_{R} \\delta}{2} - \\frac{\\Gamma_{L} \\chi_{R} \\gamma d_{L}}{2} + \\frac{3 \\Gamma_{L} \\chi_{R} \\gamma d_{R}}{2} + \\frac{\\Gamma_{R}^{2} \\chi_{R} d_{L}}{2} - \\frac{\\Gamma_{R} \\chi_{I} d_{L} \\delta}{2} - \\frac{3 \\Gamma_{R} \\chi_{R} \\gamma d_{L}}{2} + \\frac{\\Gamma_{R} \\chi_{R} \\gamma d_{R}}{2} + \\frac{\\chi_{I} \\gamma d_{L} \\delta}{2} - \\frac{\\chi_{I} \\gamma d_{R} \\delta}{2} + \\chi_{R} \\gamma^{2} d_{L} - \\chi_{R} \\gamma^{2} d_{R}}{\\Gamma_{L}^{3} \\Gamma_{R} - \\Gamma_{L}^{3} \\gamma + 2 \\Gamma_{L}^{2} \\Gamma_{R}^{2} - 7 \\Gamma_{L}^{2} \\Gamma_{R} \\gamma + \\Gamma_{L}^{2} \\chi_{I}^{2} + \\Gamma_{L}^{2} \\chi_{R}^{2} + 5 \\Gamma_{L}^{2} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R}^{3} - 7 \\Gamma_{L} \\Gamma_{R}^{2} \\gamma + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{I}^{2} + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{R}^{2} + 14 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R} \\delta^{2} - 4 \\Gamma_{L} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{L} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{L} \\gamma^{3} - \\Gamma_{L} \\gamma \\delta^{2} - \\Gamma_{R}^{3} \\gamma + \\Gamma_{R}^{2} \\chi_{I}^{2} + \\Gamma_{R}^{2} \\chi_{R}^{2} + 5 \\Gamma_{R}^{2} \\gamma^{2} - 4 \\Gamma_{R} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{R} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{R} \\gamma^{3} - \\Gamma_{R} \\gamma \\delta^{2} + 4 \\chi_{I}^{2} \\gamma^{2} + 4 \\chi_{R}^{2} \\gamma^{2} + 4 \\gamma^{4} + \\gamma^{2} \\delta^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\langle \\hat{\\tilde{P}}_{L}{\\left(t \\right)}\\hat{\\tilde{P}}_{R}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{P}}_{R}{\\left(t \\right)}\\hat{\\tilde{P}}_{L}{\\left(t \\right)} \\rangle}{2} = \\frac{- \\frac{\\Gamma_{L}^{2} \\chi_{I} d_{R}}{2} + \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{I} d_{L}}{2} - \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{I} d_{R}}{2} - \\frac{\\Gamma_{L} \\chi_{I} \\gamma d_{L}}{2} + \\frac{3 \\Gamma_{L} \\chi_{I} \\gamma d_{R}}{2} - \\frac{\\Gamma_{L} \\chi_{R} d_{R} \\delta}{2} + \\frac{\\Gamma_{R}^{2} \\chi_{I} d_{L}}{2} - \\frac{3 \\Gamma_{R} \\chi_{I} \\gamma d_{L}}{2} + \\frac{\\Gamma_{R} \\chi_{I} \\gamma d_{R}}{2} + \\frac{\\Gamma_{R} \\chi_{R} d_{L} \\delta}{2} + \\chi_{I} \\gamma^{2} d_{L} - \\chi_{I} \\gamma^{2} d_{R} - \\frac{\\chi_{R} \\gamma d_{L} \\delta}{2} + \\frac{\\chi_{R} \\gamma d_{R} \\delta}{2}}{\\Gamma_{L}^{3} \\Gamma_{R} - \\Gamma_{L}^{3} \\gamma + 2 \\Gamma_{L}^{2} \\Gamma_{R}^{2} - 7 \\Gamma_{L}^{2} \\Gamma_{R} \\gamma + \\Gamma_{L}^{2} \\chi_{I}^{2} + \\Gamma_{L}^{2} \\chi_{R}^{2} + 5 \\Gamma_{L}^{2} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R}^{3} - 7 \\Gamma_{L} \\Gamma_{R}^{2} \\gamma + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{I}^{2} + 2 \\Gamma_{L} \\Gamma_{R} \\chi_{R}^{2} + 14 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R} \\delta^{2} - 4 \\Gamma_{L} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{L} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{L} \\gamma^{3} - \\Gamma_{L} \\gamma \\delta^{2} - \\Gamma_{R}^{3} \\gamma + \\Gamma_{R}^{2} \\chi_{I}^{2} + \\Gamma_{R}^{2} \\chi_{R}^{2} + 5 \\Gamma_{R}^{2} \\gamma^{2} - 4 \\Gamma_{R} \\chi_{I}^{2} \\gamma - 4 \\Gamma_{R} \\chi_{R}^{2} \\gamma - 8 \\Gamma_{R} \\gamma^{3} - \\Gamma_{R} \\gamma \\delta^{2} + 4 \\chi_{I}^{2} \\gamma^{2} + 4 \\chi_{R}^{2} \\gamma^{2} + 4 \\gamma^{4} + \\gamma^{2} \\delta^{2}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# all expressions\n",
    "expr_dV_dt = A * V + V * A.T + D\n",
    "# unique expressions\n",
    "idxs = [(0, 0), (0, 1), (1, 1), (2, 2), (2, 3), (3, 3), (0, 2), (0, 3), (1, 2), (1, 3)]\n",
    "expr_unique = [expr_dV_dt[idx[0], idx[1]] for idx in idxs]\n",
    "# for i in range(len(idxs)):\n",
    "#     display(Math(latex(expr_unique[i]) + ' = 0'))\n",
    "\n",
    "# solutions\n",
    "sols = solve(expr_unique, [V[idx[0], idx[1]] for idx in idxs])\n",
    "for idx in idxs:\n",
    "    display(Math(latex(V[idx[0], idx[1]]) + ' = ' + latex(sols[V[idx[0], idx[1]]])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{\\langle \\hat{\\tilde{P}}_{L}{\\left(t \\right)}\\hat{\\tilde{P}}_{L}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{P}}_{R}{\\left(t \\right)}\\hat{\\tilde{P}}_{R}{\\left(t \\right)} \\rangle}{2} - \\frac{\\langle \\hat{\\tilde{Q}}_{L}{\\left(t \\right)}\\hat{\\tilde{Q}}_{L}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{Q}}_{R}{\\left(t \\right)}\\hat{\\tilde{Q}}_{R}{\\left(t \\right)} \\rangle}{2} = - \\frac{\\gamma \\kappa^{2} \\left(\\Gamma_{L} - \\Gamma_{R}\\right) \\left(n_{th} + 1\\right) \\left(\\Gamma_{L}^{2} + 2 \\Gamma_{L} \\Gamma_{R} - 4 \\Gamma_{L} \\gamma + \\Gamma_{R}^{2} - 4 \\Gamma_{R} \\gamma + 4 \\gamma^{2} + \\delta^{2}\\right)}{\\Gamma_{L} \\Gamma_{R} \\lambda^{2} \\left(\\Gamma_{L}^{2} + 2 \\Gamma_{L} \\Gamma_{R} - 4 \\Gamma_{L} \\gamma + \\Gamma_{R}^{2} - 4 \\Gamma_{R} \\gamma + 4 \\gamma^{2}\\right) + \\kappa^{2} \\left(\\Gamma_{L}^{3} \\Gamma_{R} - \\Gamma_{L}^{3} \\gamma + 2 \\Gamma_{L}^{2} \\Gamma_{R}^{2} - 7 \\Gamma_{L}^{2} \\Gamma_{R} \\gamma + 5 \\Gamma_{L}^{2} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R}^{3} - 7 \\Gamma_{L} \\Gamma_{R}^{2} \\gamma + 14 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R} \\delta^{2} - 8 \\Gamma_{L} \\gamma^{3} - \\Gamma_{L} \\gamma \\delta^{2} - \\Gamma_{R}^{3} \\gamma + 5 \\Gamma_{R}^{2} \\gamma^{2} - 8 \\Gamma_{R} \\gamma^{3} - \\Gamma_{R} \\gamma \\delta^{2} + 4 \\gamma^{4} + \\gamma^{2} \\delta^{2}\\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# substitution list\n",
    "list_subs = [\n",
    "    (d_L, conjugate(eta_L) * eta_L * (1 + lamb**2 / kappa**2) / 2 + gamma * (2 * n_th + 1)),\n",
    "    (d_R, conjugate(eta_R) * eta_R * (1 + lamb**2 / kappa**2) / 2 + gamma * (2 * n_th + 1)),\n",
    "    (chi_R, lamb / (kappa**2 + lamb**2) * (G_L * conjugate(G_R) + conjugate(G_L) * G_R) / 2),\n",
    "    (chi_I, I * lamb / (kappa**2 + lamb**2) * (G_L * conjugate(G_R) - conjugate(G_L) * G_R) / 2),\n",
    "    (eta_L, sqrt(2 * kappa) * kappa * G_L / (kappa**2 + lamb**2)),\n",
    "    (eta_LR, sqrt(2 * kappa) * kappa / (kappa**2 + lamb**2) * (conjugate(G_L) + G_L) / 2),\n",
    "    (eta_LI, I * sqrt(2 * kappa) * kappa / (kappa**2 + lamb**2) * (conjugate(G_L) - G_L) / 2),\n",
    "    (eta_R, sqrt(2 * kappa) * kappa * G_R / (kappa**2 + lamb**2)),\n",
    "    (eta_RR, sqrt(2 * kappa) * kappa / (kappa**2 + lamb**2) * (conjugate(G_R) + G_R) / 2),\n",
    "    (eta_RI, I * sqrt(2 * kappa) * kappa / (kappa**2 + lamb**2) * (conjugate(G_R) - G_R) / 2),\n",
    "    (Gamma_L, conjugate(G_L) * kappa * G_L / (kappa**2 + lamb**2)),\n",
    "    (Gamma_R, conjugate(G_R) * kappa * G_R / (kappa**2 + lamb**2))\n",
    "]\n",
    "\n",
    "# phonon number difference\n",
    "expr_n_b_diff = (sols[V[3, 3]] + sols[V[2, 2]] - sols[V[1, 1]] - sols[V[0, 0]]) / 2\n",
    "expr_n_b_diff = expr_n_b_diff.subs(list_subs).cancel()\n",
    "# display(Math('n_{b_R} - n_{b_L} = ' + latex(expr_n_b_diff)))\n",
    "\n",
    "# resubstitute\n",
    "expr_n_b_diff = expr_n_b_diff.subs([\n",
    "    (conjugate(G_L) * G_L, Gamma_L * (kappa**2 + lamb**2) / kappa), \n",
    "    (conjugate(G_R) * G_R, Gamma_R * (kappa**2 + lamb**2) / kappa)\n",
    "]).expand().factor().collect([kappa**2, lamb**2 * Gamma_L * Gamma_R])\n",
    "display(Math(latex((V[3, 3] + V[2, 2]) / 2 - (V[1, 1] + V[0, 0]) / 2) + ' = ' + latex(expr_n_b_diff)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\Gamma_{L}^{3} \\Gamma_{R} - \\Gamma_{L}^{3} \\gamma + 2 \\Gamma_{L}^{2} \\Gamma_{R}^{2} - 7 \\Gamma_{L}^{2} \\Gamma_{R} \\gamma + 5 \\Gamma_{L}^{2} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R}^{3} - 7 \\Gamma_{L} \\Gamma_{R}^{2} \\gamma + 14 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R} \\delta^{2} - 8 \\Gamma_{L} \\gamma^{3} - \\Gamma_{L} \\gamma \\delta^{2} - \\Gamma_{R}^{3} \\gamma + 5 \\Gamma_{R}^{2} \\gamma^{2} - 8 \\Gamma_{R} \\gamma^{3} - \\Gamma_{R} \\gamma \\delta^{2} + 4 \\gamma^{4} + \\gamma^{2} \\delta^{2}$"
      ],
      "text/plain": [
       "          3                        3                      2           2       \n",
       "\\Gamma_{L} ⋅\\Gamma_{R} - \\Gamma_{L} ⋅\\gamma + 2⋅\\Gamma_{L} ⋅\\Gamma_{R}  - 7⋅\\G\n",
       "\n",
       "        2                                 2       2                        3  \n",
       "amma_{L} ⋅\\Gamma_{R}⋅\\gamma + 5⋅\\Gamma_{L} ⋅\\gamma  + \\Gamma_{L}⋅\\Gamma_{R}  -\n",
       "\n",
       "                        2                                         2           \n",
       " 7⋅\\Gamma_{L}⋅\\Gamma_{R} ⋅\\gamma + 14⋅\\Gamma_{L}⋅\\Gamma_{R}⋅\\gamma  + \\Gamma_{\n",
       "\n",
       "               2                      3                      2             3  \n",
       "L}⋅\\Gamma_{R}⋅δ  - 8⋅\\Gamma_{L}⋅\\gamma  - \\Gamma_{L}⋅\\gamma⋅δ  - \\Gamma_{R} ⋅\\\n",
       "\n",
       "                    2       2                      3                      2   \n",
       "gamma + 5⋅\\Gamma_{R} ⋅\\gamma  - 8⋅\\Gamma_{R}⋅\\gamma  - \\Gamma_{R}⋅\\gamma⋅δ  + \n",
       "\n",
       "        4         2  2\n",
       "4⋅\\gamma  + \\gamma ⋅δ "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr = (Gamma_L - gamma) * (Gamma_R - gamma) * ((Gamma_L + Gamma_R - 2 * gamma)**2 + delta**2)\n",
    "expr.expand()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\langle \\hat{\\tilde{P}}_{L}{\\left(t \\right)}\\hat{\\tilde{P}}_{R}{\\left(t \\right)} \\rangle}{2} + \\frac{\\langle \\hat{\\tilde{P}}_{R}{\\left(t \\right)}\\hat{\\tilde{P}}_{L}{\\left(t \\right)} \\rangle}{2} = - \\frac{\\gamma \\kappa^{2} \\lambda \\left(\\Gamma_{L} - \\Gamma_{R}\\right) \\left(n_{th} + 1\\right) \\left(G_{L} \\left(i \\Gamma_{L} + i \\Gamma_{R} - 2 i \\gamma + \\delta\\right) \\overline{G_{R}} + G_{R} \\left(- i \\Gamma_{L} - i \\Gamma_{R} + 2 i \\gamma + \\delta\\right) \\overline{G_{L}}\\right)}{2 \\left(\\kappa^{2} + \\lambda^{2}\\right) \\left(\\Gamma_{L} \\Gamma_{R} \\lambda^{2} \\left(\\Gamma_{L}^{2} + 2 \\Gamma_{L} \\Gamma_{R} - 4 \\Gamma_{L} \\gamma + \\Gamma_{R}^{2} - 4 \\Gamma_{R} \\gamma + 4 \\gamma^{2}\\right) + \\kappa^{2} \\left(\\Gamma_{L}^{3} \\Gamma_{R} - \\Gamma_{L}^{3} \\gamma + 2 \\Gamma_{L}^{2} \\Gamma_{R}^{2} - 7 \\Gamma_{L}^{2} \\Gamma_{R} \\gamma + 5 \\Gamma_{L}^{2} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R}^{3} - 7 \\Gamma_{L} \\Gamma_{R}^{2} \\gamma + 14 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} + \\Gamma_{L} \\Gamma_{R} \\delta^{2} - 8 \\Gamma_{L} \\gamma^{3} - \\Gamma_{L} \\gamma \\delta^{2} - \\Gamma_{R}^{3} \\gamma + 5 \\Gamma_{R}^{2} \\gamma^{2} - 8 \\Gamma_{R} \\gamma^{3} - \\Gamma_{R} \\gamma \\delta^{2} + 4 \\gamma^{4} + \\gamma^{2} \\delta^{2}\\right)\\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# phase correlations\n",
    "expr_phase_corr = sols[V[1, 3]].subs(list_subs).cancel().subs([\n",
    "    (conjugate(G_L) * G_L, Gamma_L * (kappa**2 + lamb**2) / kappa), \n",
    "    (conjugate(G_R) * G_R, Gamma_R * (kappa**2 + lamb**2) / kappa)\n",
    "]).expand().factor().collect([\n",
    "    kappa**2, \n",
    "    lamb**2 * Gamma_L * Gamma_R, \n",
    "    G_L * conjugate(G_R), \n",
    "    G_R * conjugate(G_L)\n",
    "])\n",
    "display(Math(latex(V[1, 3]) + ' = ' + latex(expr_phase_corr)))"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "b3ba2566441a7c06988d0923437866b63cedc61552a5af99d1f4fb67d367b25f"
  },
  "kernelspec": {
   "display_name": "Python 3.8.5 64-bit ('base': conda)",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}