notebooks/uni_00_adiabatic_elimination.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adiabatic Elimination in the Simple Unidirectionally-coupled Configuration (Uni_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",
    "eta, g_L, g_R, gamma_L, gamma_R, kappa_L, kappa_R, omega_L, omega_R = symbols('\\\\eta, 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 = symbols('\\\\Delta_{L}, \\\\Delta_{R}', 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 = symbols('\\\\alpha_{L}, \\\\alpha_{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}}_{R}^{in} (t) = \\sqrt{\\eta} \\left(- \\sqrt{2} \\sqrt{\\kappa_{L}} \\hat{a}_{L}{\\left(t \\right)} + \\hat{a}_{L}^{in}{\\left(t \\right)}\\right) + \\sqrt{1 - \\eta} \\hat{a}_{R}^{in}{\\left(t \\right)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "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)} + \\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 \\Delta_{R} \\hat{a}_{R}{\\left(t \\right)} + i \\alpha_{R} g_{R} \\hat{b}_{R}{\\left(t \\right)} + i \\alpha_{R} g_{R} \\hat{b}_{R}{\\left(t \\right)}^{\\dagger} - 2 \\sqrt{\\eta} \\sqrt{\\kappa_{L}} \\sqrt{\\kappa_{R}} \\hat{a}_{L}{\\left(t \\right)} + \\sqrt{2} \\sqrt{\\eta} \\sqrt{\\kappa_{R}} \\hat{a}_{L}^{in}{\\left(t \\right)} + \\sqrt{2} \\sqrt{\\kappa_{R}} \\sqrt{1 - \\eta} \\hat{a}_{R}^{in}{\\left(t \\right)} - \\kappa_{R} \\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": [
    "# unidirectional transfer\n",
    "expr_a_L_out_t = a_L_in_t - sqrt(2 * kappa_L) * a_L_t\n",
    "expr_a_R_in_t = sqrt(eta) * expr_a_L_out_t + sqrt(1 - eta) * a_R_in_t \n",
    "display(Math('\\\\dot{\\\\hat{a}}_{R}^{in} (t) = ' + latex(expr_a_R_in_t)))\n",
    "\n",
    "# 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) + 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) + sqrt(2 * kappa_R) * a_R_in_t\n",
    "expr_da_R_dt_t = expr_da_R_dt_t.subs(a_R_in_t, expr_a_R_in_t).expand()\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)}$"
      ],
      "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)} - 2 \\sqrt{\\eta} \\sqrt{\\kappa_{L}} \\sqrt{\\kappa_{R}} \\hat{\\tilde{a}}_{L}{\\left(t \\right)} + \\sqrt{2} \\sqrt{\\eta} \\sqrt{\\kappa_{R}} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)} + \\sqrt{2} \\sqrt{\\kappa_{R}} \\sqrt{1 - \\eta} \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)} - \\kappa_{R} \\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}}}_{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).expand()\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).expand()\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).expand()\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).expand()\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)}$"
      ],
      "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} - 2 \\sqrt{\\eta} \\sqrt{\\kappa_{L}} \\sqrt{\\kappa_{R}} \\hat{\\tilde{a}}_{L}{\\left(t \\right)} + \\sqrt{2} \\sqrt{\\eta} \\sqrt{\\kappa_{R}} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)} + \\sqrt{2} \\sqrt{\\kappa_{R}} \\sqrt{1 - \\eta} \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)} - \\kappa_{R} \\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}}}_{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} \\hat{\\tilde{b}}_{L}{\\left(t \\right)}^{\\dagger}}{\\kappa_{L}} + \\frac{\\sqrt{2} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}}{\\sqrt{\\kappa_{L}}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\hat{\\tilde{a}}_{R} = - \\frac{2 i G_{L} \\sqrt{\\eta} \\hat{\\tilde{b}}_{L}{\\left(t \\right)}^{\\dagger}}{\\sqrt{\\kappa_{L}} \\sqrt{\\kappa_{R}}} + \\frac{i G_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)}^{\\dagger}}{\\kappa_{R}} - \\frac{\\sqrt{2} \\sqrt{\\eta} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}}{\\sqrt{\\kappa_{R}}} + \\frac{\\sqrt{2 - 2 \\eta} \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}}{\\sqrt{\\kappa_{R}}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\dot{\\hat{\\tilde{b}}}_{L} (t) = \\frac{G_{L} \\overline{G_{L}} \\hat{\\tilde{b}}_{L}{\\left(t \\right)}}{\\kappa_{L}} + \\frac{\\sqrt{2} i G_{L} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger}}{\\sqrt{\\kappa_{L}}} + \\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} i G_{R} \\sqrt{\\eta} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger}}{\\sqrt{\\kappa_{R}}} - \\frac{2 G_{R} \\sqrt{\\eta} \\overline{G_{L}} \\hat{\\tilde{b}}_{L}{\\left(t \\right)}}{\\sqrt{\\kappa_{L}} \\sqrt{\\kappa_{R}}} + \\frac{G_{R} \\overline{G_{R}} \\hat{\\tilde{b}}_{R}{\\left(t \\right)}}{\\kappa_{R}} + \\frac{\\sqrt{2} i G_{R} \\sqrt{1 - \\eta} \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}^{\\dagger}}{\\sqrt{\\kappa_{R}}} + \\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# 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",
    "    (conjugate(sqrt(2 - 2 * eta)), sqrt(2) * sqrt(1 - eta))\n",
    "]\n",
    "# mechanical mode equations\n",
    "expr_db_L_tilde_dt_adia_t = expr_db_L_tilde_dt_approx_t.subs(list_subs).expand().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().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{2 G_{R} \\sqrt{\\eta} \\overline{G_{L}}}{\\sqrt{\\kappa_{L}} \\sqrt{\\kappa_{R}}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\Gamma_{L} = \\frac{G_{L} \\overline{G_{L}}}{\\kappa_{L}}, \\quad \\Gamma_{R} = \\frac{G_{R} \\overline{G_{R}}}{\\kappa_{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}}}, \\quad \\eta_{R} = \\frac{\\sqrt{2} G_{R}}{\\sqrt{\\kappa_{R}}}$"
      ],
      "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)} + i \\eta_{L} \\hat{\\tilde{a}}_{L}^{in}{\\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{b}}}_{R} (t) = \\Gamma_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)} - \\chi \\hat{\\tilde{b}}_{L}{\\left(t \\right)} + \\eta_{R} \\left(- i \\sqrt{\\eta} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger} + i \\sqrt{1 - \\eta} \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}^{\\dagger}\\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 = 2 * sqrt(eta) * conjugate(G_L) * G_R / sqrt(kappa_L * kappa_R)\n",
    "expr_eta_L = sqrt(2) * G_L / sqrt(kappa_L)\n",
    "expr_eta_R = sqrt(2) * G_R / sqrt(kappa_R)\n",
    "expr_Gamma_L = G_L * conjugate(G_L) / kappa_L\n",
    "expr_Gamma_R = G_R * conjugate(G_R) / kappa_R\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, eta_L),\n",
    "    (expr_eta_R, eta_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)} + i \\eta_{L} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger} + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{b}}_{L}^{in}{\\left(t \\right)} - \\gamma \\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) = \\Gamma_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)} - \\chi \\hat{\\tilde{b}}_{L}{\\left(t \\right)} + \\eta_{R} \\left(- i \\sqrt{\\eta} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger} + i \\sqrt{1 - \\eta} \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}^{\\dagger}\\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) = \\left(2 \\Gamma_{L} - 2 \\gamma\\right) \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\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}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle}{d t} (t) = - \\chi \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\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)} - \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\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) = - \\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)} - \\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} = - \\frac{2 G_{L} \\sqrt{\\eta} \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle \\overline{G_{R}}}{\\kappa} - \\frac{2 G_{R} \\sqrt{\\eta} \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle \\overline{G_{L}}}{\\kappa} + \\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 G_{L} \\sqrt{\\eta} \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{R} \\rangle \\overline{G_{R}}}{\\kappa} - \\frac{2 G_{R} \\sqrt{\\eta} \\langle \\hat{\\tilde{b}}_{R}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle \\overline{G_{L}}}{\\kappa} + \\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{2 i G_{R} \\sqrt{\\eta} \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle \\overline{G_{L}}}{\\kappa \\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{2 i G_{L} \\sqrt{\\eta} \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle \\overline{G_{R}}}{\\kappa \\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{\\sqrt{\\eta} \\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)}{\\kappa \\left(\\Gamma_{R} - \\gamma\\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{4 \\Gamma_{L} \\Gamma_{R} \\eta \\langle \\hat{\\tilde{b}}_{L}^{\\dagger} \\hat{\\tilde{b}}_{L} \\rangle \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\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),\n",
    "    (conjugate(G_R) * G_R, Gamma_R * kappa)\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)} - \\eta_{L_{I}} \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)} + \\eta_{L_{R}} \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)} + \\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)} + \\eta_{L_{I}} \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)} + \\eta_{L_{R}} \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)} + \\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{P}}_{L}{\\left(t \\right)} - \\chi_{R} \\hat{\\tilde{Q}}_{L}{\\left(t \\right)} + \\sqrt{\\eta} \\eta_{R_{I}} \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)} - \\sqrt{\\eta} \\eta_{R_{R}} \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)} - \\eta_{R_{I}} \\sqrt{1 - \\eta} \\hat{\\tilde{X}}_{R}^{in}{\\left(t \\right)} + \\eta_{R_{R}} \\sqrt{1 - \\eta} \\hat{\\tilde{Y}}_{R}^{in}{\\left(t \\right)} + \\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{Q}}_{L}{\\left(t \\right)} - \\chi_{R} \\hat{\\tilde{P}}_{L}{\\left(t \\right)} - \\sqrt{\\eta} \\eta_{R_{I}} \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)} - \\sqrt{\\eta} \\eta_{R_{R}} \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)} + \\eta_{R_{I}} \\sqrt{1 - \\eta} \\hat{\\tilde{Y}}_{R}^{in}{\\left(t \\right)} + \\eta_{R_{R}} \\sqrt{1 - \\eta} \\hat{\\tilde{X}}_{R}^{in}{\\left(t \\right)} + \\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",
    "    (conjugate(sqrt(1 - eta)), sqrt(1 - eta))\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} \\overline{G_{L}}}{\\kappa_{L}}, \\quad \\Gamma_{R} = \\frac{G_{R} \\overline{G_{R}}}{\\kappa_{R}}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\chi = \\frac{2 G_{R} \\sqrt{\\eta} \\overline{G_{L}}}{\\sqrt{\\kappa_{L}} \\sqrt{\\kappa_{R}}} = 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}}} = i \\eta_{L_{I}} + \\eta_{L_{R}}, \\quad \\eta_{R} = \\frac{\\sqrt{2} G_{R}}{\\sqrt{\\kappa_{R}}} = 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} + \\sqrt{2} \\sqrt{\\gamma} \\hat{\\tilde{b}}_{L}^{in}{\\left(t \\right)} - \\gamma \\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) = \\Gamma_{R} \\hat{\\tilde{b}}_{R}{\\left(t \\right)} - \\chi \\hat{\\tilde{b}}_{L}{\\left(t \\right)} - i \\sqrt{\\eta} \\eta_{R} \\hat{\\tilde{a}}_{L}^{in}{\\left(t \\right)}^{\\dagger} + i \\eta_{R} \\sqrt{1 - \\eta} \\hat{\\tilde{a}}_{R}^{in}{\\left(t \\right)}^{\\dagger} + \\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)} - \\eta_{L_{I}} \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)} + \\eta_{L_{R}} \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)} + \\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)} + \\eta_{L_{I}} \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)} + \\eta_{L_{R}} \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)} + \\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{P}}_{L}{\\left(t \\right)} - \\chi_{R} \\hat{\\tilde{Q}}_{L}{\\left(t \\right)} + \\sqrt{\\eta} \\eta_{R_{I}} \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)} - \\sqrt{\\eta} \\eta_{R_{R}} \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)} - \\eta_{R_{I}} \\sqrt{1 - \\eta} \\hat{\\tilde{X}}_{R}^{in}{\\left(t \\right)} + \\eta_{R_{R}} \\sqrt{1 - \\eta} \\hat{\\tilde{Y}}_{R}^{in}{\\left(t \\right)} + \\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{Q}}_{L}{\\left(t \\right)} - \\chi_{R} \\hat{\\tilde{P}}_{L}{\\left(t \\right)} - \\sqrt{\\eta} \\eta_{R_{I}} \\hat{\\tilde{Y}}_{L}^{in}{\\left(t \\right)} - \\sqrt{\\eta} \\eta_{R_{R}} \\hat{\\tilde{X}}_{L}^{in}{\\left(t \\right)} + \\eta_{R_{I}} \\sqrt{1 - \\eta} \\hat{\\tilde{Y}}_{R}^{in}{\\left(t \\right)} + \\eta_{R_{R}} \\sqrt{1 - \\eta} \\hat{\\tilde{X}}_{R}^{in}{\\left(t \\right)} + \\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": "markdown",
   "metadata": {},
   "source": [
    "## Drift Matrix and Correlation Matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle A = \\left[\\begin{matrix}\\Gamma_{L} - \\gamma & 0 & 0 & 0\\\\0 & \\Gamma_{L} - \\gamma & 0 & 0\\\\- \\chi_{R} & \\chi_{I} & \\Gamma_{R} - \\gamma & \\delta\\\\- \\chi_{I} & - \\chi_{R} & - \\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}\\frac{\\eta_{L_{I}}^{2}}{2} + \\frac{\\eta_{L_{R}}^{2}}{2} + \\gamma \\left(2 n_{th} + 1\\right) & 0 & - \\frac{\\sqrt{\\eta} \\left(\\eta_{L_{I}} \\eta_{R_{I}} + \\eta_{L_{R}} \\eta_{R_{R}}\\right)}{2} & - \\frac{\\sqrt{\\eta} \\left(- \\eta_{L_{I}} \\eta_{R_{R}} + \\eta_{L_{R}} \\eta_{R_{I}}\\right)}{2}\\\\0 & \\frac{\\eta_{L_{I}}^{2}}{2} + \\frac{\\eta_{L_{R}}^{2}}{2} + \\gamma \\left(2 n_{th} + 1\\right) & \\frac{\\sqrt{\\eta} \\left(- \\eta_{L_{I}} \\eta_{R_{R}} + \\eta_{L_{R}} \\eta_{R_{I}}\\right)}{2} & - \\frac{\\sqrt{\\eta} \\left(\\eta_{L_{I}} \\eta_{R_{I}} + \\eta_{L_{R}} \\eta_{R_{R}}\\right)}{2}\\\\- \\frac{\\sqrt{\\eta} \\left(\\eta_{L_{I}} \\eta_{R_{I}} + \\eta_{L_{R}} \\eta_{R_{R}}\\right)}{2} & \\frac{\\sqrt{\\eta} \\left(- \\eta_{L_{I}} \\eta_{R_{R}} + \\eta_{L_{R}} \\eta_{R_{I}}\\right)}{2} & \\frac{\\eta_{R_{I}}^{2}}{2} + \\frac{\\eta_{R_{R}}^{2}}{2} + \\gamma \\left(2 n_{th} + 1\\right) & 0\\\\- \\frac{\\sqrt{\\eta} \\left(- \\eta_{L_{I}} \\eta_{R_{R}} + \\eta_{L_{R}} \\eta_{R_{I}}\\right)}{2} & - \\frac{\\sqrt{\\eta} \\left(\\eta_{L_{I}} \\eta_{R_{I}} + \\eta_{L_{R}} \\eta_{R_{R}}\\right)}{2} & 0 & \\frac{\\eta_{R_{I}}^{2}}{2} + \\frac{\\eta_{R_{R}}^{2}}{2} + \\gamma \\left(2 n_{th} + 1\\right)\\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, 0, 0], \n",
    "    [0, (Gamma_L - gamma), 0, 0], \n",
    "    [- chi_R, chi_I, (Gamma_R - gamma), delta], \n",
    "    [- chi_I, - chi_R, - delta, (Gamma_R - gamma)]\n",
    "])\n",
    "display(Math('A = ' + latex(A)))\n",
    "\n",
    "# noise matrix\n",
    "diag_L = eta_LR**2 / 2 + eta_LI**2 / 2 + gamma * (2 * n_th + 1)\n",
    "diag_R = eta_RR**2 / 2 + eta_RI**2 / 2 + gamma * (2 * n_th + 1)\n",
    "temp_same = - sqrt(eta) / 2 * (eta_LR * eta_RR + eta_LI * eta_RI)\n",
    "temp_diff = - sqrt(eta) / 2 * (eta_LR * eta_RI - eta_LI * eta_RR)\n",
    "D = Matrix([[diag_L, 0, temp_same, temp_diff], [0, diag_L, - temp_diff, temp_same], [temp_same, - temp_diff, diag_R, 0], [temp_diff, temp_same, 0, diag_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": "markdown",
   "metadata": {},
   "source": [
    "## Analytical Expressions for Quadrature Correlations"
   ]
  },
  {
   "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{\\eta_{L_{I}}^{2}}{4} + \\frac{\\eta_{L_{R}}^{2}}{4} + \\gamma n_{th} + \\frac{\\gamma}{2}}{\\Gamma_{L} - \\gamma}$"
      ],
      "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{\\eta_{L_{I}}^{2}}{4} + \\frac{\\eta_{L_{R}}^{2}}{4} + \\gamma n_{th} + \\frac{\\gamma}{2}}{\\Gamma_{L} - \\gamma}$"
      ],
      "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} \\eta_{R_{I}}^{2}}{4} - \\frac{\\Gamma_{L}^{3} \\eta_{R_{R}}^{2}}{4} - \\Gamma_{L}^{3} \\gamma n_{th} - \\frac{\\Gamma_{L}^{3} \\gamma}{2} - \\frac{\\Gamma_{L}^{2} \\Gamma_{R} \\eta_{R_{I}}^{2}}{2} - \\frac{\\Gamma_{L}^{2} \\Gamma_{R} \\eta_{R_{R}}^{2}}{2} - 2 \\Gamma_{L}^{2} \\Gamma_{R} \\gamma n_{th} - \\Gamma_{L}^{2} \\Gamma_{R} \\gamma - \\frac{\\Gamma_{L}^{2} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}}}{2} + \\frac{\\Gamma_{L}^{2} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}}}{2} + \\frac{\\Gamma_{L}^{2} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}}}{2} + \\frac{\\Gamma_{L}^{2} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}}}{2} + \\frac{5 \\Gamma_{L}^{2} \\eta_{R_{I}}^{2} \\gamma}{4} + \\frac{5 \\Gamma_{L}^{2} \\eta_{R_{R}}^{2} \\gamma}{4} + 5 \\Gamma_{L}^{2} \\gamma^{2} n_{th} + \\frac{5 \\Gamma_{L}^{2} \\gamma^{2}}{2} - \\frac{\\Gamma_{L} \\Gamma_{R}^{2} \\eta_{R_{I}}^{2}}{4} - \\frac{\\Gamma_{L} \\Gamma_{R}^{2} \\eta_{R_{R}}^{2}}{4} - \\Gamma_{L} \\Gamma_{R}^{2} \\gamma n_{th} - \\frac{\\Gamma_{L} \\Gamma_{R}^{2} \\gamma}{2} - \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}}}{2} + \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}}}{2} + \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}}}{2} + \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}}}{2} + \\frac{3 \\Gamma_{L} \\Gamma_{R} \\eta_{R_{I}}^{2} \\gamma}{2} + \\frac{3 \\Gamma_{L} \\Gamma_{R} \\eta_{R_{R}}^{2} \\gamma}{2} + 6 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} n_{th} + 3 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} - \\frac{\\Gamma_{L} \\chi_{I}^{2} \\eta_{L_{I}}^{2}}{4} - \\frac{\\Gamma_{L} \\chi_{I}^{2} \\eta_{L_{R}}^{2}}{4} - \\Gamma_{L} \\chi_{I}^{2} \\gamma n_{th} - \\frac{\\Gamma_{L} \\chi_{I}^{2} \\gamma}{2} + \\frac{\\Gamma_{L} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}} \\delta}{2} + \\frac{3 \\Gamma_{L} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}} \\gamma}{2} - \\frac{3 \\Gamma_{L} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}} \\gamma}{2} + \\frac{\\Gamma_{L} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}} \\delta}{2} - \\frac{\\Gamma_{L} \\chi_{R}^{2} \\eta_{L_{I}}^{2}}{4} - \\frac{\\Gamma_{L} \\chi_{R}^{2} \\eta_{L_{R}}^{2}}{4} - \\Gamma_{L} \\chi_{R}^{2} \\gamma n_{th} - \\frac{\\Gamma_{L} \\chi_{R}^{2} \\gamma}{2} - \\frac{3 \\Gamma_{L} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}} \\gamma}{2} + \\frac{\\Gamma_{L} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}} \\delta}{2} - \\frac{\\Gamma_{L} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}} \\delta}{2} - \\frac{3 \\Gamma_{L} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}} \\gamma}{2} - 2 \\Gamma_{L} \\eta_{R_{I}}^{2} \\gamma^{2} - \\frac{\\Gamma_{L} \\eta_{R_{I}}^{2} \\delta^{2}}{4} - 2 \\Gamma_{L} \\eta_{R_{R}}^{2} \\gamma^{2} - \\frac{\\Gamma_{L} \\eta_{R_{R}}^{2} \\delta^{2}}{4} - 8 \\Gamma_{L} \\gamma^{3} n_{th} - 4 \\Gamma_{L} \\gamma^{3} - \\Gamma_{L} \\gamma \\delta^{2} n_{th} - \\frac{\\Gamma_{L} \\gamma \\delta^{2}}{2} + \\frac{\\Gamma_{R}^{2} \\eta_{R_{I}}^{2} \\gamma}{4} + \\frac{\\Gamma_{R}^{2} \\eta_{R_{R}}^{2} \\gamma}{4} + \\Gamma_{R}^{2} \\gamma^{2} n_{th} + \\frac{\\Gamma_{R}^{2} \\gamma^{2}}{2} - \\frac{\\Gamma_{R} \\chi_{I}^{2} \\eta_{L_{I}}^{2}}{4} - \\frac{\\Gamma_{R} \\chi_{I}^{2} \\eta_{L_{R}}^{2}}{4} - \\Gamma_{R} \\chi_{I}^{2} \\gamma n_{th} - \\frac{\\Gamma_{R} \\chi_{I}^{2} \\gamma}{2} + \\frac{\\Gamma_{R} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}} \\gamma}{2} - \\frac{\\Gamma_{R} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}} \\gamma}{2} - \\frac{\\Gamma_{R} \\chi_{R}^{2} \\eta_{L_{I}}^{2}}{4} - \\frac{\\Gamma_{R} \\chi_{R}^{2} \\eta_{L_{R}}^{2}}{4} - \\Gamma_{R} \\chi_{R}^{2} \\gamma n_{th} - \\frac{\\Gamma_{R} \\chi_{R}^{2} \\gamma}{2} - \\frac{\\Gamma_{R} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}} \\gamma}{2} - \\frac{\\Gamma_{R} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}} \\gamma}{2} - \\Gamma_{R} \\eta_{R_{I}}^{2} \\gamma^{2} - \\Gamma_{R} \\eta_{R_{R}}^{2} \\gamma^{2} - 4 \\Gamma_{R} \\gamma^{3} n_{th} - 2 \\Gamma_{R} \\gamma^{3} + \\frac{\\chi_{I}^{2} \\eta_{L_{I}}^{2} \\gamma}{2} + \\frac{\\chi_{I}^{2} \\eta_{L_{R}}^{2} \\gamma}{2} + 2 \\chi_{I}^{2} \\gamma^{2} n_{th} + \\chi_{I}^{2} \\gamma^{2} - \\frac{\\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}} \\gamma \\delta}{2} - \\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}} \\gamma^{2} + \\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}} \\gamma^{2} - \\frac{\\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}} \\gamma \\delta}{2} + \\frac{\\chi_{R}^{2} \\eta_{L_{I}}^{2} \\gamma}{2} + \\frac{\\chi_{R}^{2} \\eta_{L_{R}}^{2} \\gamma}{2} + 2 \\chi_{R}^{2} \\gamma^{2} n_{th} + \\chi_{R}^{2} \\gamma^{2} + \\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}} \\gamma^{2} - \\frac{\\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}} \\gamma \\delta}{2} + \\frac{\\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}} \\gamma \\delta}{2} + \\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}} \\gamma^{2} + \\eta_{R_{I}}^{2} \\gamma^{3} + \\frac{\\eta_{R_{I}}^{2} \\gamma \\delta^{2}}{4} + \\eta_{R_{R}}^{2} \\gamma^{3} + \\frac{\\eta_{R_{R}}^{2} \\gamma \\delta^{2}}{4} + 4 \\gamma^{4} n_{th} + 2 \\gamma^{4} + \\gamma^{2} \\delta^{2} n_{th} + \\frac{\\gamma^{2} \\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 + 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": [
       "<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} \\eta_{R_{I}}^{2}}{4} - \\frac{\\Gamma_{L}^{3} \\eta_{R_{R}}^{2}}{4} - \\Gamma_{L}^{3} \\gamma n_{th} - \\frac{\\Gamma_{L}^{3} \\gamma}{2} - \\frac{\\Gamma_{L}^{2} \\Gamma_{R} \\eta_{R_{I}}^{2}}{2} - \\frac{\\Gamma_{L}^{2} \\Gamma_{R} \\eta_{R_{R}}^{2}}{2} - 2 \\Gamma_{L}^{2} \\Gamma_{R} \\gamma n_{th} - \\Gamma_{L}^{2} \\Gamma_{R} \\gamma - \\frac{\\Gamma_{L}^{2} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}}}{2} + \\frac{\\Gamma_{L}^{2} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}}}{2} + \\frac{\\Gamma_{L}^{2} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}}}{2} + \\frac{\\Gamma_{L}^{2} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}}}{2} + \\frac{5 \\Gamma_{L}^{2} \\eta_{R_{I}}^{2} \\gamma}{4} + \\frac{5 \\Gamma_{L}^{2} \\eta_{R_{R}}^{2} \\gamma}{4} + 5 \\Gamma_{L}^{2} \\gamma^{2} n_{th} + \\frac{5 \\Gamma_{L}^{2} \\gamma^{2}}{2} - \\frac{\\Gamma_{L} \\Gamma_{R}^{2} \\eta_{R_{I}}^{2}}{4} - \\frac{\\Gamma_{L} \\Gamma_{R}^{2} \\eta_{R_{R}}^{2}}{4} - \\Gamma_{L} \\Gamma_{R}^{2} \\gamma n_{th} - \\frac{\\Gamma_{L} \\Gamma_{R}^{2} \\gamma}{2} - \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}}}{2} + \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}}}{2} + \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}}}{2} + \\frac{\\Gamma_{L} \\Gamma_{R} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}}}{2} + \\frac{3 \\Gamma_{L} \\Gamma_{R} \\eta_{R_{I}}^{2} \\gamma}{2} + \\frac{3 \\Gamma_{L} \\Gamma_{R} \\eta_{R_{R}}^{2} \\gamma}{2} + 6 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} n_{th} + 3 \\Gamma_{L} \\Gamma_{R} \\gamma^{2} - \\frac{\\Gamma_{L} \\chi_{I}^{2} \\eta_{L_{I}}^{2}}{4} - \\frac{\\Gamma_{L} \\chi_{I}^{2} \\eta_{L_{R}}^{2}}{4} - \\Gamma_{L} \\chi_{I}^{2} \\gamma n_{th} - \\frac{\\Gamma_{L} \\chi_{I}^{2} \\gamma}{2} + \\frac{\\Gamma_{L} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}} \\delta}{2} + \\frac{3 \\Gamma_{L} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}} \\gamma}{2} - \\frac{3 \\Gamma_{L} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}} \\gamma}{2} + \\frac{\\Gamma_{L} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}} \\delta}{2} - \\frac{\\Gamma_{L} \\chi_{R}^{2} \\eta_{L_{I}}^{2}}{4} - \\frac{\\Gamma_{L} \\chi_{R}^{2} \\eta_{L_{R}}^{2}}{4} - \\Gamma_{L} \\chi_{R}^{2} \\gamma n_{th} - \\frac{\\Gamma_{L} \\chi_{R}^{2} \\gamma}{2} - \\frac{3 \\Gamma_{L} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}} \\gamma}{2} + \\frac{\\Gamma_{L} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}} \\delta}{2} - \\frac{\\Gamma_{L} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}} \\delta}{2} - \\frac{3 \\Gamma_{L} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}} \\gamma}{2} - 2 \\Gamma_{L} \\eta_{R_{I}}^{2} \\gamma^{2} - \\frac{\\Gamma_{L} \\eta_{R_{I}}^{2} \\delta^{2}}{4} - 2 \\Gamma_{L} \\eta_{R_{R}}^{2} \\gamma^{2} - \\frac{\\Gamma_{L} \\eta_{R_{R}}^{2} \\delta^{2}}{4} - 8 \\Gamma_{L} \\gamma^{3} n_{th} - 4 \\Gamma_{L} \\gamma^{3} - \\Gamma_{L} \\gamma \\delta^{2} n_{th} - \\frac{\\Gamma_{L} \\gamma \\delta^{2}}{2} + \\frac{\\Gamma_{R}^{2} \\eta_{R_{I}}^{2} \\gamma}{4} + \\frac{\\Gamma_{R}^{2} \\eta_{R_{R}}^{2} \\gamma}{4} + \\Gamma_{R}^{2} \\gamma^{2} n_{th} + \\frac{\\Gamma_{R}^{2} \\gamma^{2}}{2} - \\frac{\\Gamma_{R} \\chi_{I}^{2} \\eta_{L_{I}}^{2}}{4} - \\frac{\\Gamma_{R} \\chi_{I}^{2} \\eta_{L_{R}}^{2}}{4} - \\Gamma_{R} \\chi_{I}^{2} \\gamma n_{th} - \\frac{\\Gamma_{R} \\chi_{I}^{2} \\gamma}{2} + \\frac{\\Gamma_{R} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}} \\gamma}{2} - \\frac{\\Gamma_{R} \\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}} \\gamma}{2} - \\frac{\\Gamma_{R} \\chi_{R}^{2} \\eta_{L_{I}}^{2}}{4} - \\frac{\\Gamma_{R} \\chi_{R}^{2} \\eta_{L_{R}}^{2}}{4} - \\Gamma_{R} \\chi_{R}^{2} \\gamma n_{th} - \\frac{\\Gamma_{R} \\chi_{R}^{2} \\gamma}{2} - \\frac{\\Gamma_{R} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}} \\gamma}{2} - \\frac{\\Gamma_{R} \\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}} \\gamma}{2} - \\Gamma_{R} \\eta_{R_{I}}^{2} \\gamma^{2} - \\Gamma_{R} \\eta_{R_{R}}^{2} \\gamma^{2} - 4 \\Gamma_{R} \\gamma^{3} n_{th} - 2 \\Gamma_{R} \\gamma^{3} + \\frac{\\chi_{I}^{2} \\eta_{L_{I}}^{2} \\gamma}{2} + \\frac{\\chi_{I}^{2} \\eta_{L_{R}}^{2} \\gamma}{2} + 2 \\chi_{I}^{2} \\gamma^{2} n_{th} + \\chi_{I}^{2} \\gamma^{2} - \\frac{\\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}} \\gamma \\delta}{2} - \\chi_{I} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}} \\gamma^{2} + \\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}} \\gamma^{2} - \\frac{\\chi_{I} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}} \\gamma \\delta}{2} + \\frac{\\chi_{R}^{2} \\eta_{L_{I}}^{2} \\gamma}{2} + \\frac{\\chi_{R}^{2} \\eta_{L_{R}}^{2} \\gamma}{2} + 2 \\chi_{R}^{2} \\gamma^{2} n_{th} + \\chi_{R}^{2} \\gamma^{2} + \\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{I}} \\gamma^{2} - \\frac{\\chi_{R} \\sqrt{\\eta} \\eta_{L_{I}} \\eta_{R_{R}} \\gamma \\delta}{2} + \\frac{\\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{I}} \\gamma \\delta}{2} + \\chi_{R} \\sqrt{\\eta} \\eta_{L_{R}} \\eta_{R_{R}} \\gamma^{2} + \\eta_{R_{I}}^{2} \\gamma^{3} + \\frac{\\eta_{R_{I}}^{2} \\gamma \\delta^{2}}{4} + \\eta_{R_{R}}^{2} \\gamma^{3} + \\frac{\\eta_{R_{R}}^{2} \\gamma \\delta^{2}}{4} + 4 \\gamma^{4} n_{th} + 2 \\gamma^{4} + \\gamma^{2} \\delta^{2} n_{th} + \\frac{\\gamma^{2} \\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 + 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": [
       "<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{2 \\sqrt{\\eta} \\delta \\left(\\Gamma_{L} - \\gamma\\right) \\left(\\eta_{L_{I}} \\eta_{R_{R}} - \\eta_{L_{R}} \\eta_{R_{I}}\\right) \\left(\\delta^{4} + \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2}\\right) + 2 \\sqrt{\\eta} \\left(\\Gamma_{L} - \\gamma\\right) \\left(\\eta_{L_{I}} \\eta_{R_{I}} + \\eta_{L_{R}} \\eta_{R_{R}}\\right) \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right) \\left(\\delta^{4} + \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2}\\right) + \\left(\\eta_{L_{I}}^{2} + \\eta_{L_{R}}^{2} + 4 \\gamma n_{th} + 2 \\gamma\\right) \\left(\\chi_{I} \\delta^{5} + \\chi_{I} \\delta^{3} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + \\chi_{I} \\delta^{3} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\chi_{I} \\delta \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} - \\chi_{R} \\delta^{4} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right) - \\chi_{R} \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right) - \\chi_{R} \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{3} - \\chi_{R} \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{3}\\right)}{4 \\left(\\Gamma_{L} - \\gamma\\right) \\left(\\delta^{6} + \\delta^{4} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + 2 \\delta^{4} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + 2 \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{4} + \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{4}\\right)}$"
      ],
      "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{2 \\sqrt{\\eta} \\delta \\left(\\Gamma_{L} - \\gamma\\right) \\left(\\eta_{L_{I}} \\eta_{R_{I}} + \\eta_{L_{R}} \\eta_{R_{R}}\\right) \\left(\\delta^{4} + \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2}\\right) - 2 \\sqrt{\\eta} \\left(\\Gamma_{L} - \\gamma\\right) \\left(\\eta_{L_{I}} \\eta_{R_{R}} - \\eta_{L_{R}} \\eta_{R_{I}}\\right) \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right) \\left(\\delta^{4} + \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2}\\right) - \\left(\\eta_{L_{I}}^{2} + \\eta_{L_{R}}^{2} + 4 \\gamma n_{th} + 2 \\gamma\\right) \\left(\\chi_{I} \\delta^{4} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right) + \\chi_{I} \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right) + \\chi_{I} \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{3} + \\chi_{I} \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{3} + \\chi_{R} \\delta^{5} + \\chi_{R} \\delta^{3} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + \\chi_{R} \\delta^{3} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\chi_{R} \\delta \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2}\\right)}{4 \\left(\\Gamma_{L} - \\gamma\\right) \\left(\\delta^{6} + \\delta^{4} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + 2 \\delta^{4} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + 2 \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{4} + \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{4}\\right)}$"
      ],
      "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{- 2 \\sqrt{\\eta} \\delta \\left(\\Gamma_{L} - \\gamma\\right) \\left(\\eta_{L_{I}} \\eta_{R_{I}} + \\eta_{L_{R}} \\eta_{R_{R}}\\right) \\left(\\delta^{4} + \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2}\\right) + 2 \\sqrt{\\eta} \\left(\\Gamma_{L} - \\gamma\\right) \\left(\\eta_{L_{I}} \\eta_{R_{R}} - \\eta_{L_{R}} \\eta_{R_{I}}\\right) \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right) \\left(\\delta^{4} + \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2}\\right) + \\left(\\eta_{L_{I}}^{2} + \\eta_{L_{R}}^{2} + 4 \\gamma n_{th} + 2 \\gamma\\right) \\left(\\chi_{I} \\delta^{4} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right) + \\chi_{I} \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right) + \\chi_{I} \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{3} + \\chi_{I} \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{3} + \\chi_{R} \\delta^{5} + \\chi_{R} \\delta^{3} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + \\chi_{R} \\delta^{3} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\chi_{R} \\delta \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2}\\right)}{4 \\left(\\Gamma_{L} - \\gamma\\right) \\left(\\delta^{6} + \\delta^{4} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + 2 \\delta^{4} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + 2 \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{4} + \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{4}\\right)}$"
      ],
      "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{2 \\sqrt{\\eta} \\delta \\left(\\Gamma_{L} - \\gamma\\right) \\left(\\eta_{L_{I}} \\eta_{R_{R}} - \\eta_{L_{R}} \\eta_{R_{I}}\\right) \\left(\\delta^{4} + \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2}\\right) + 2 \\sqrt{\\eta} \\left(\\Gamma_{L} - \\gamma\\right) \\left(\\eta_{L_{I}} \\eta_{R_{I}} + \\eta_{L_{R}} \\eta_{R_{R}}\\right) \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right) \\left(\\delta^{4} + \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2}\\right) + \\left(\\eta_{L_{I}}^{2} + \\eta_{L_{R}}^{2} + 4 \\gamma n_{th} + 2 \\gamma\\right) \\left(\\chi_{I} \\delta^{5} + \\chi_{I} \\delta^{3} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + \\chi_{I} \\delta^{3} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\chi_{I} \\delta \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} - \\chi_{R} \\delta^{4} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right) - \\chi_{R} \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right) - \\chi_{R} \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{3} - \\chi_{R} \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{3}\\right)}{4 \\left(\\Gamma_{L} - \\gamma\\right) \\left(\\delta^{6} + \\delta^{4} \\left(\\Gamma_{R} - \\gamma\\right)^{2} + 2 \\delta^{4} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + 2 \\delta^{2} \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{2} + \\delta^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{4} + \\left(\\Gamma_{R} - \\gamma\\right)^{2} \\left(\\Gamma_{L} + \\Gamma_{R} - 2 \\gamma\\right)^{4}\\right)}$"
      ],
      "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": "markdown",
   "metadata": {},
   "source": [
    "## Difference in Phonon Numbers"
   ]
  },
  {
   "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 \\left(n_{th} + 1\\right) \\left(\\Gamma_{L}^{3} + 4 \\Gamma_{L}^{2} \\Gamma_{R} \\eta + \\Gamma_{L}^{2} \\Gamma_{R} - 4 \\Gamma_{L}^{2} \\gamma + 4 \\Gamma_{L} \\Gamma_{R}^{2} \\eta - \\Gamma_{L} \\Gamma_{R}^{2} - 8 \\Gamma_{L} \\Gamma_{R} \\eta \\gamma + 4 \\Gamma_{L} \\gamma^{2} + \\Gamma_{L} \\delta^{2} - \\Gamma_{R}^{3} + 4 \\Gamma_{R}^{2} \\gamma - 4 \\Gamma_{R} \\gamma^{2} - \\Gamma_{R} \\delta^{2}\\right)}{\\left(\\Gamma_{L} - \\gamma\\right) \\left(\\Gamma_{R} - \\gamma\\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)}$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# substitution list\n",
    "list_subs = [\n",
    "    (chi_R, sqrt(eta) / kappa * (G_L * conjugate(G_R) + conjugate(G_L) * G_R)),\n",
    "    (chi_I, I * sqrt(eta) / kappa * (G_L * conjugate(G_R) - conjugate(G_L) * G_R)),\n",
    "    (eta_LR, sqrt(2 / kappa) * (conjugate(G_L) + G_L) / 2),\n",
    "    (eta_LI, I * sqrt(2 / kappa) * (conjugate(G_L) - G_L) / 2),\n",
    "    (eta_RR, sqrt(2 / kappa) * (conjugate(G_R) + G_R) / 2),\n",
    "    (eta_RI, I * sqrt(2 / kappa) * (conjugate(G_R) - G_R) / 2),\n",
    "    (Gamma_L, conjugate(G_L) * G_L / kappa),\n",
    "    (Gamma_R, conjugate(G_R) * G_R / kappa)\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).expand().factor()\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), \n",
    "    (conjugate(G_R) * G_R, Gamma_R * kappa)\n",
    "]).expand().factor()\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} + 4 \\Gamma_{L}^{2} \\Gamma_{R} \\eta + \\Gamma_{L}^{2} \\Gamma_{R} - 4 \\Gamma_{L}^{2} \\gamma + 4 \\Gamma_{L} \\Gamma_{R}^{2} \\eta - \\Gamma_{L} \\Gamma_{R}^{2} - 8 \\Gamma_{L} \\Gamma_{R} \\eta \\gamma + 4 \\Gamma_{L} \\gamma^{2} + \\Gamma_{L} \\delta^{2} - \\Gamma_{R}^{3} + 4 \\Gamma_{R}^{2} \\gamma - 4 \\Gamma_{R} \\gamma^{2} - \\Gamma_{R} \\delta^{2}$"
      ],
      "text/plain": [
       "          3               2                             2                     \n",
       "\\Gamma_{L}  + 4⋅\\Gamma_{L} ⋅\\Gamma_{R}⋅\\eta + \\Gamma_{L} ⋅\\Gamma_{R} - 4⋅\\Gamm\n",
       "\n",
       "     2                                 2                             2        \n",
       "a_{L} ⋅\\gamma + 4⋅\\Gamma_{L}⋅\\Gamma_{R} ⋅\\eta - \\Gamma_{L}⋅\\Gamma_{R}  - 8⋅\\Ga\n",
       "\n",
       "                                                    2               2         \n",
       "mma_{L}⋅\\Gamma_{R}⋅\\eta⋅\\gamma + 4⋅\\Gamma_{L}⋅\\gamma  + \\Gamma_{L}⋅δ  - \\Gamma\n",
       "\n",
       "    3               2                             2               2\n",
       "_{R}  + 4⋅\\Gamma_{R} ⋅\\gamma - 4⋅\\Gamma_{R}⋅\\gamma  - \\Gamma_{R}⋅δ "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expr = (Gamma_L - Gamma_R) * ((Gamma_L + Gamma_R - 2 * gamma)**2 + delta**2) + 4 * Gamma_L * Gamma_R * eta * (Gamma_L + Gamma_R - 2 * gamma)\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{\\sqrt{\\eta} \\gamma \\left(n_{th} + 1\\right) \\left(G_{L} \\left(- \\Gamma_{L} - \\Gamma_{R} + 2 \\gamma + i \\delta\\right) \\overline{G_{R}} + G_{R} \\left(- \\Gamma_{L} - \\Gamma_{R} + 2 \\gamma - i \\delta\\right) \\overline{G_{L}}\\right)}{\\kappa \\left(\\Gamma_{L} - \\gamma\\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)}$"
      ],
      "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), \n",
    "    (conjugate(G_R) * G_R, Gamma_R * kappa)\n",
    "]).expand().factor().collect([\n",
    "    kappa**2, \n",
    "    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
}