Bfield_perturbation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "\n",
    "from astropy import units as u\n",
    "# from astropy.constants import c\n",
    "from astropy.constants import k_B, G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "MathJax.Hub.Config({\n",
       "    TeX: { equationNumbers: { autoNumber: \"AMS\" } }\n",
       "});\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%javascript\n",
    "MathJax.Hub.Config({\n",
    "    TeX: { equationNumbers: { autoNumber: \"AMS\" } }\n",
    "});"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using MHD wave relations and observables\n",
    "\n",
    "We can use equation 10.19 from Stahler & Palla's book\n",
    "\\begin{equation}\n",
    "\\delta u_n = \\delta u_i \\left(1 - \\frac{i\\omega}{n_i\\langle \\sigma_{in} u_i' \\rangle}\\right)^{-1}~,\n",
    "\\label{eq:1}\n",
    "\\end{equation}\n",
    "where \n",
    "$n_i\\langle \\sigma_{in} u_i' \\rangle$\n",
    "is \"the frequency with which a given natural atom or molecule is struck by ions\", $\\delta u_n$ is the perturbation on the neutral's velocity, and $\\delta u_i$ is the perturbation on the ion's velocity.\n",
    "\n",
    "The other relation to use is equation 10.21 from Stahler & Palla's book\n",
    "\\begin{equation}\n",
    "\\frac{\\omega^2}{k^2} = \\frac{B_o^2}{4\\pi \\rho_o} \n",
    "                       \\left(1 - \\frac{i\\omega}{n_i\\langle \\sigma_{in} u_i' \\rangle}\\right)~,\n",
    "\\label{eq:2}\n",
    "\\end{equation}\n",
    "which relates the $\\omega$ and $k$.\n",
    "\n",
    "Finally, equation 10.17 from Stahler & Palla's book relates the perturbation in the magnetic field as\n",
    "\\begin{equation}\n",
    "\\delta u_i = -\\frac{\\omega}{k} \\frac{\\delta B}{B_o}~.\n",
    "\\label{eq:3}\n",
    "\\end{equation}\n",
    "\n",
    "It is useful to define\n",
    "\\begin{equation}\n",
    "\\omega_o = n_i\\langle \\sigma_{in} u_i' \\rangle~,\n",
    "\\end{equation}\n",
    "and to assume that $\\omega$ is real. The dissipation timescale \n",
    "\n",
    "We estimate the velocity perturvations as the velocity dispersions, $\\delta u =\\sqrt{2}\\, \\sigma_{\\rm NT}$, derived from the spectral lines and therefore rewrite equation \\ref{eq:1} as \n",
    "\\begin{equation}\n",
    "\\left| \\frac{\\delta u_i}{\\delta u_n} \\right|^2 =\n",
    "1 + \\left(\\frac{\\omega}{\\omega_o}\\right)^{2}\n",
    "\\approx \\left(\\frac{\\sigma_{\\rm NT}({\\rm N_2H^+})}{\\sigma_{\\rm NT}({\\rm NH_3})}\\right)^2~.\n",
    "\\label{eq:4}\n",
    "\\end{equation}\n",
    "\n",
    "Now, equation \\ref{eq:2} can also be rewritten as\n",
    "\\begin{equation}\n",
    "\\left|\\frac{\\omega^2}{k^2}\\right| = \\left(\\frac{B_o^2}{4\\pi \\rho_o}\\right) \n",
    "                                    \\left[1 + \\left(\\frac{\\omega}{\\omega_o}\\right)^2\\right]^{1/2}\n",
    "                            \\approx \\frac{B_o^2}{4\\pi \\rho_o}\n",
    "                                    \\frac{\\sigma_{\\rm NT}({\\rm N_2H^+})}{\\sigma_{\\rm NT}({\\rm NH_3})}~.\n",
    "\\label{eq:5}\n",
    "\\end{equation}\n",
    "\n",
    "Combining equations \\ref{eq:3}, \\ref{eq:4} and \\ref{eq:5} we obtain\n",
    "$$\n",
    "\\frac{\\delta B}{B_o} = \\sqrt{\\sigma_{\\rm NT}({\\rm N_2H^+})\\, \\sigma_{\\rm NT}({\\rm NH_3})}\n",
    "\\frac{\\sqrt{8\\pi \\rho_o}}{B_o}\n",
    "$$\n",
    "or\n",
    "$$\n",
    "\\delta B= \\sqrt{8\\pi \\rho_o} \\sqrt{\\sigma_{\\rm NT}({\\rm N_2H^+})\\, \\sigma_{\\rm NT}({\\rm NH_3})}\n",
    "$$\n",
    "\n",
    "The term $n_i$ is estimated using the ionization degree from Caselli et al. (2002)\n",
    "$$\n",
    "x(e) = \\frac{n_i}{n({\\rm H_2})} = 5.2\\times 10^{-6} \\left( \\frac{n({\\rm H_2})}{{\\rm cm}^{-3}}\\right)^{-0.56} \n",
    "$$\n",
    "and the term\n",
    "$$\n",
    "\\langle \\sigma_{in} u_i' \\rangle = 1.69\\times 10^{-9} {\\rm cm^3 s^{-1}}\n",
    "$$\n",
    "is approximated by the Langevin term for HCO$^+$-H$_2$ collisions (McDaniels & Mason, 1973).\n",
    "\n",
    "### Wavelengths\n",
    "Using the relation \n",
    "$$\n",
    "\\lambda = \\frac{2\\pi}{k} = \\frac{2\\pi}{\\omega} \\frac{\\omega}{k}\n",
    "$$ \n",
    "and that \n",
    "$$\n",
    "\\frac{1}{\\omega} = \\frac{1}{\\omega_o}\n",
    "\\sqrt{\\frac{\\sigma_{\\rm NT}({\\rm NH_3})^2}{\\sigma_{\\rm NT}({\\rm N_2H^+})^2-\\sigma_{\\rm NT}({\\rm NH_3})^2}}\n",
    "$$\n",
    "then we can write\n",
    "\\begin{equation}\n",
    "\\lambda = \\sqrt{\\frac{\\pi}{\\rho_o}}  \\frac{B_o}{\\omega_o} \n",
    "\\sqrt{\\frac{\\sigma_{\\rm NT}({\\rm NH_3})\\, \\sigma_{\\rm NT}({\\rm N_2H^+})}{\\sigma_{\\rm NT}({\\rm N_2H^+})^2 - \\sigma_{\\rm NT}({\\rm NH_3})^2}}\n",
    "\\label{eq:lambda_in}\n",
    "\\end{equation}\n",
    "\n",
    "For a wave to propagate, the wavelength must be larger than a critical value (equation 10.23 from Stahler & Palla's book and equation 79 from Mouschovias et al. 2011)\n",
    "$$\n",
    "\\lambda_{min} = \\sqrt{\\frac{\\pi}{4\\rho_o}} \\frac{B_o}{\\omega_o}\n",
    "$$\n",
    "and therefore\n",
    "\\begin{equation}\n",
    "\\frac{\\lambda}{\\lambda_{min}} = 2\n",
    "\\sqrt{\\frac{\\sigma_{\\rm NT}({\\rm NH_3})\\, \\sigma_{\\rm NT}({\\rm N_2H^+})}{\\sigma_{\\rm NT}({\\rm N_2H^+})^2 - \\sigma_{\\rm NT}({\\rm NH_3})^2}}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "### Damping timescale\n",
    "In addition, the characteristic damping timescale is given by equation 80 in Mouschovias et al. (2011), \n",
    "\\begin{equation}\n",
    "    \\tau_d = \\frac{\\lambda^2\\,\\omega_o}{2\\pi^2\\,v_A^2}~,\n",
    "\\end{equation}\n",
    "where the Alfvén velocity is a function of density and magnetic field strength,\n",
    "\\begin{equation}\n",
    "    v_A^2 = \\frac{B_o^2}{4\\pi\\,\\rho_o}~. \\label{eq:Va}\n",
    "\\end{equation}\n",
    "Replacing relations \\ref{eq:lambda_in} and \\ref{eq:Va} we obtain\n",
    "\\begin{equation}\n",
    "    \\tau_d = \\frac{2}{\\omega_o}\\left(\\frac{\\sigma_{\\rm NT}({\\rm NH_3})\\, \\sigma_{\\rm NT}({\\rm N_2H^+})}{\\sigma_{\\rm NT}({\\rm N_2H^+})^2 - \\sigma_{\\rm NT}({\\rm NH_3})^2}\\right) = 0.2~\\textrm{Myr,}\n",
    "\\end{equation}\n",
    "which is much smaller than the crossing time of the coherent region\n",
    "\\begin{equation}\n",
    "    \\tau_{cross} = \n",
    "    \\frac{\\rm diameter}{c_s} \n",
    "    = 1.8~\\textrm{Myr,}\n",
    "\\end{equation}\n",
    "and therefore the wave must be continually injected if they are to persist for a crossing time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigma_thermal(mu_mol, tk=10*u.K):\n",
    "    \"\"\"\n",
    "    Returns the sound speed for temperature Tk and molecular weight mu.\n",
    "    This is also used to determine the thermal velocity dispersion of\n",
    "    a molecular transition.\n",
    "\n",
    "    \"\"\"\n",
    "    return np.sqrt(k_B * tk/(mu_mol * u.u)).to(u.km/u.s)\n",
    "\n",
    "\n",
    "def density_ion(dens_all, do_Caselli=True):\n",
    "    \"\"\"\n",
    "    1.3e-5 x n(H2)^{0.5}  (from McKee 1989) or 5.2e-6 x n*H2)^{0.44} \n",
    "    \"\"\"\n",
    "    if do_Caselli:\n",
    "        xe = 5.2e-6 * (dens_all/(u.cm**-3))**-0.56\n",
    "    else:\n",
    "        xe = 1.3e-5 * (dens_all/(u.cm**-3))**-0.5\n",
    "    return xe*dens_all\n",
    "\n",
    "\n",
    "def get_omega(sigma_ion=0.1*u.km/u.s, sigma_neutral=0.08*u.km/u.s, density=1e6/u.cm**3, do_Caselli=True):\n",
    "    \"\"\"\n",
    "    \"\"\"\n",
    "    n_i = density_ion(density, do_Caselli=do_Caselli)\n",
    "    return (sig_in_v_i * n_i * np.sqrt((sigma_ion/sigma_neutral)**2 - 1)).decompose()\n",
    "\n",
    "\n",
    "def get_wavelength(Bfield=100*u.uG, sigma_ion=0.1*u.km/u.s, sigma_neutral=0.08*u.km/u.s, \n",
    "                   n_H2=1e6/u.cm**3, do_Caselli=True):\n",
    "    \"\"\"\n",
    "    \"\"\"\n",
    "    rho0 = (n_H2 * u.u * 2.8).cgs\n",
    "    n_i = density_ion(n_H2, do_Caselli=do_Caselli).cgs\n",
    "    my_wave = np.sqrt(np.pi/rho0) * (Bfield.to(u.G)/(n_i*sig_in_v_i)) * np.sqrt(sigma_ion.cgs*sigma_neutral.cgs/(sigma_ion.cgs**2 - sigma_neutral.cgs**2))\n",
    "    return (my_wave.value) * u.cm\n",
    "\n",
    "\n",
    "def alven_speed(Bfield=100*u.G, n_H2=1e6/u.cm**3):\n",
    "    \"\"\"\n",
    "    \"\"\"\n",
    "    rho_o = (n_H2 * u.u * 2.8).cgs\n",
    "    V_a = Bfield.to(u.G)/np.sqrt(4*np.pi*rho_o)\n",
    "    return (V_a.value)*(u.cm/u.s)\n",
    "\n",
    "\n",
    "def delta_B(Bfield=100*u.uG, sigma_ion=0.1*u.km/u.s, sigma_neutral=0.08*u.km/u.s, \n",
    "                   n_H2=1e6/u.cm**3):\n",
    "    rho_o = (n_H2 * u.u * 2.8).cgs\n",
    "    delta_B_Bo = np.sqrt(sigma_ion * sigma_neutral) * np.sqrt(8 * np.pi * rho_o) / Bfield.to((u.g/u.cm)**(1/2)/u.s, equivalencies=equiv_B)\n",
    "    return delta_B_Bo.decompose()\n",
    "\n",
    "\n",
    "def wave_ratio(sigma_ion=0.1*u.km/u.s, sigma_neutral=0.08*u.km/u.s):\n",
    "    return 2*np.sqrt(sigma_ion * sigma_neutral / (sigma_ion**2 - sigma_neutral**2))\n",
    "\n",
    "\n",
    "def wave_damp(sigma_ion=0.1*u.km/u.s, sigma_neutral=0.08*u.km/u.s, density=1e6/u.cm**3, do_Caselli=True):\n",
    "    omega_o = get_omega(sigma_ion=sigma_ion, sigma_neutral=sigma_neutral, density=density, do_Caselli=do_Caselli)\n",
    "    return (2*(sigma_ion * sigma_neutral / (sigma_ion**2 - sigma_neutral**2)) / omega_o).to(u.Myr)\n",
    "\n",
    "\n",
    "sig_in_v_i = 1.69e-9*u.cm**3/u.s\n",
    "gauss_B = (u.g/u.cm)**(0.5)/u.s\n",
    "gauss2_B = (u.g/u.cm)/u.s**2\n",
    "equiv_B = [(u.G, gauss_B, lambda x: x, lambda x: x)]\n",
    "equiv_B2 = [(u.G**2, gauss2_B, lambda x: x, lambda x: x)]\n",
    "\n",
    "B_o = np.array([500.]) * u.uG\n",
    "c_sound = sigma_thermal(2.37, tk=9.7*u.K).cgs\n",
    "sigma_ion = 0.109 * u.km / u.s\n",
    "sigma_neutral = 0.088 * u.km / u.s\n",
    "n_H2 = 7e4 / u.cm**3\n",
    "rho_o = (n_H2 * u.u * 2.8).cgs\n",
    "n_H2_out = 5e3 / u.cm**3\n",
    "rho_o_out = (n_H2_out * u.u * 2.8).cgs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.05]\n"
     ]
    }
   ],
   "source": [
    "delta_B_Bo = np.round(delta_B(Bfield=B_o, sigma_ion=sigma_ion, n_H2=n_H2), decimals=2)\n",
    "print(delta_B_Bo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[25.] uG\n"
     ]
    }
   ],
   "source": [
    "delta_B = np.round(delta_B_Bo * B_o)\n",
    "print(delta_B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This corresponds to a variation of $25\\,\\mu$G, which corresponds to at least 5\\% of the field ($\\delta B/B_o$).\n",
    "\n",
    "### Wavelengths\n",
    "If we want to estimate the wavelength then we will need to estimate \n",
    "$n_i\\langle \\sigma_{in} u_i' \\rangle$ to derive $\\omega$ and then obtain $k$.\n",
    "\n",
    "We use the velocity dispersions for ions and neutrals, the mean density of the coherent core, and the magnetic field strengths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$[0.64] \\; \\mathrm{pc}$"
      ],
      "text/plain": [
       "<Quantity [0.64] pc>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.round((get_wavelength(Bfield=B_o, sigma_ion=sigma_ion, sigma_neutral=sigma_neutral, n_H2=n_H2)).to(u.pc), decimals=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Therefore, the wavelength is $<$0.64 pc, where the uncertainty is dominated by the magnetic field strength. The size of the coherent region is 0.34 pc (in diameter)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ratio between the wavelength and the Alfven cutoff wavelength is then:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$3.05 \\; \\mathrm{}$"
      ],
      "text/plain": [
       "<Quantity 3.05>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.round(wave_ratio(sigma_ion=sigma_ion, sigma_neutral=sigma_neutral), decimals=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculation of timescales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$0.2 \\; \\mathrm{Myr}$"
      ],
      "text/plain": [
       "<Quantity 0.2 Myr>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.round(wave_damp(sigma_ion=sigma_ion, sigma_neutral=sigma_neutral, density=n_H2, do_Caselli=True), decimals=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The crossing time is \n",
    "$$\\tau_{cross} = \\frac{diameter}{c_s}~,\n",
    "$$\n",
    "and for Barnard 5 the diameter is 0.34 pc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$1.8 \\; \\mathrm{Myr}$"
      ],
      "text/plain": [
       "<Quantity 1.8 Myr>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.round((0.34*u.pc/c_sound).to(u.Myr), decimals=1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}