.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "api/gallery/plasma/waves/plot_04_debye_screening.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_api_gallery_plasma_waves_plot_04_debye_screening.py: Debye and Hückel: screening and the Debye length =================================================== Debye and Hückel (1923), studying electrolytes, asked what the potential around one ion looks like when the other mobile charges are free to rearrange. Charges of the opposite sign gather around it, those of the same sign are pushed away, and the Boltzmann-distributed cloud cancels the ion's field beyond a characteristic distance, the Debye length .. math:: \lambda_D = \sqrt{\frac{\varepsilon_0 k_BT}{n e^2}} = \frac{v_{th}}{\omega_p}. Linearizing the Poisson-Boltzmann equation gives the screened, Yukawa potential :math:`\phi = \frac{q}{4\pi\varepsilon_0 r}e^{-r/\lambda_D}`. Plasma physics adopted the result unchanged. This example solves the full nonlinear Poisson-Boltzmann equation around a test charge, compares it with the Debye-Hückel form, and tabulates :math:`\lambda_D` and the number of particles in a Debye sphere for real plasmas using :func:`~physicskit.plasma.waves.plasma_frequency`. .. GENERATED FROM PYTHON SOURCE LINES 25-32 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np from scipy.integrate import solve_bvp from physicskit.constants import ELECTRON_MASS, ELEMENTARY_CHARGE, K_B, VACUUM_PERMITTIVITY from physicskit.plasma.waves import plasma_frequency .. GENERATED FROM PYTHON SOURCE LINES 33-48 The screened potential -------------------------- In units :math:`\psi = e\phi/k_BT` and :math:`x = r/\lambda_D`, with Boltzmann electrons around a fixed positive test charge on a uniform ion background, .. math:: \psi'' + \frac{2}{x}\psi' = e^{\psi} - 1, with the bare Coulomb field :math:`\psi' = -Q/x^2` at the core and :math:`\psi\to0` far away. :math:`Q` measures the test charge's strength against the thermal energy. For small :math:`Q` the solution is Debye-Hückel's :math:`\psi = Q e^{-x}/x`; for large :math:`Q` the electrons pile up nonlinearly close in. .. GENERATED FROM PYTHON SOURCE LINES 48-86 .. code-block:: Python x_core, x_far = 0.05, 12.0 x = np.geomspace(x_core, x_far, 400) def solve(Q): def rhs(xx, y): return np.vstack([y[1], np.expm1(y[0]) - 2 * y[1] / xx]) def bc(ya, yb): return np.array([ya[1] + Q / x_core**2, yb[0]]) guess = np.vstack([Q * np.exp(-x) / x, -Q * np.exp(-x) * (1 + x) / x**2]) sol = solve_bvp(rhs, bc, x, guess, tol=1e-6, max_nodes=100000) return sol.sol(x)[0] fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4)) for Q, color in [(0.01, "steelblue"), (0.5, "darkorange"), (3.0, "firebrick")]: psi = solve(Q) debye = Q * np.exp(-x) / x ax1.loglog(x, psi / Q, color=color, label=f"Poisson-Boltzmann, Q = {Q}") ax2.semilogx(x, psi / debye, color=color, label=f"Q = {Q}") i1 = np.argmin(np.abs(x - 1)) print(f"Q = {Q:5}: psi / Debye-Huckel at r = lambda_D: {psi[i1] / debye[i1]:.3f}") ax1.loglog(x, 1 / x, "k:", label="bare Coulomb 1/x") ax1.loglog(x, np.exp(-x) / x, "k--", label=r"Debye-Hückel $e^{-x}/x$") ax1.set_ylim(1e-6, 30) ax1.set_xlabel(r"$r/\lambda_D$") ax1.set_ylabel(r"$\psi / Q$") ax1.set_title("A test charge, screened") ax1.legend(fontsize=7) ax2.axhline(1, color="k", lw=0.8) ax2.set_xlabel(r"$r/\lambda_D$") ax2.set_ylabel("Poisson-Boltzmann / Debye-Hückel") ax2.set_title("Linear theory holds for weak charges") ax2.legend(fontsize=8) fig1.tight_layout() .. image-sg:: /api/gallery/plasma/waves/images/sphx_glr_plot_04_debye_screening_001.png :alt: A test charge, screened, Linear theory holds for weak charges :srcset: /api/gallery/plasma/waves/images/sphx_glr_plot_04_debye_screening_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Q = 0.01: psi / Debye-Huckel at r = lambda_D: 0.999 Q = 0.5: psi / Debye-Huckel at r = lambda_D: 0.749 /Users/cpoli/Code/physicskit/examples/plasma/waves/plot_04_debye_screening.py:54: RuntimeWarning: overflow encountered in expm1 return np.vstack([y[1], np.expm1(y[0]) - 2 * y[1] / xx]) /Users/cpoli/Code/physicskit/.venv/lib/python3.12/site-packages/scipy/integrate/_bvp.py:41: RuntimeWarning: invalid value encountered in subtract df_dy[:, i, :] = (f_new - f0) / hi /Users/cpoli/Code/physicskit/.venv/lib/python3.12/site-packages/scipy/integrate/_bvp.py:567: RuntimeWarning: invalid value encountered in divide r_middle /= 1 + np.abs(f_middle) /Users/cpoli/Code/physicskit/.venv/lib/python3.12/site-packages/scipy/integrate/_bvp.py:568: RuntimeWarning: invalid value encountered in divide r1 /= 1 + np.abs(f1) /Users/cpoli/Code/physicskit/.venv/lib/python3.12/site-packages/scipy/integrate/_bvp.py:569: RuntimeWarning: invalid value encountered in divide r2 /= 1 + np.abs(f2) Q = 3.0: psi / Debye-Huckel at r = lambda_D: 0.468 .. GENERATED FROM PYTHON SOURCE LINES 87-95 Debye lengths of real plasmas --------------------------------- :math:`\lambda_D = v_{th}/\omega_p` with :math:`v_{th}=\sqrt{k_BT/m_e}`. A gas is a plasma only if :math:`\lambda_D` is much smaller than the system and a Debye sphere holds many particles, :math:`N_D = \tfrac43\pi n\lambda_D^3 \gg 1`, so the screening cloud is a smooth crowd rather than a few individuals. The dense solar core is the marginal case: only a handful of electrons per Debye sphere. .. GENERATED FROM PYTHON SOURCE LINES 95-123 .. code-block:: Python plasmas = { "ionosphere": (1e11, 1e3), "solar corona": (1e15, 1e6), "tokamak core": (1e20, 1e8), "solar core": (1e32, 1.5e7), "interstellar medium": (1e6, 1e4), } print(f"\n{'plasma':20s} {'n [m^-3]':>9s} {'T [K]':>8s} {'lambda_D':>12s} {'N_D':>9s}") lam, ND = {}, {} for name, (n, T) in plasmas.items(): v_th = np.sqrt(K_B * T / ELECTRON_MASS) lam[name] = v_th / plasma_frequency(n, ELEMENTARY_CHARGE, ELECTRON_MASS) ND[name] = 4 / 3 * np.pi * n * lam[name] ** 3 print(f"{name:20s} {n:9.0e} {T:8.0e} {lam[name]:10.2e} m {ND[name]:9.1e}") fig2, ax3 = plt.subplots(figsize=(6, 4)) for name in plasmas: ax3.loglog(lam[name], ND[name], "o", ms=8) ax3.annotate(name, (lam[name], ND[name]), textcoords="offset points", xytext=(6, 4), fontsize=8) ax3.axhline(1, color="k", ls=":", label=r"$N_D = 1$") ax3.set_xlabel(r"Debye length $\lambda_D$ [m]") ax3.set_ylabel(r"particles per Debye sphere $N_D$") ax3.set_title("Particles per Debye sphere") ax3.legend(fontsize=8) fig2.tight_layout() plt.show() .. image-sg:: /api/gallery/plasma/waves/images/sphx_glr_plot_04_debye_screening_002.png :alt: Particles per Debye sphere :srcset: /api/gallery/plasma/waves/images/sphx_glr_plot_04_debye_screening_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none plasma n [m^-3] T [K] lambda_D N_D ionosphere 1e+11 1e+03 6.90e-03 m 1.4e+05 solar corona 1e+15 1e+06 2.18e-03 m 4.4e+07 tokamak core 1e+20 1e+08 6.90e-05 m 1.4e+08 solar core 1e+32 2e+07 2.67e-11 m 8.0e+00 interstellar medium 1e+06 1e+04 6.90e+00 m 1.4e+09 .. GENERATED FROM PYTHON SOURCE LINES 124-128 Check ----- Weak charges are screened as Debye-Hückel's e^(-r/lambda_D)/r; strong ones beyond it (nonlinear Poisson-Boltzmann). lambda_D = sqrt(eps0 k T / n e^2), N_D = (4/3) pi n lambda_D^3. .. GENERATED FROM PYTHON SOURCE LINES 128-134 .. code-block:: Python i1 = np.argmin(np.abs(x - 1)) ratios = {Q: solve(Q)[i1] / (Q * np.exp(-x[i1]) / x[i1]) for Q in (0.01, 3.0)} assert abs(ratios[0.01] - 1) < 0.005 and ratios[3.0] < 0.6 for name, (n, T) in plasmas.items(): assert abs(lam[name] / np.sqrt(VACUUM_PERMITTIVITY * K_B * T / (n * ELEMENTARY_CHARGE**2)) - 1) < 1e-10 assert abs(ND[name] - 4 / 3 * np.pi * n * lam[name] ** 3) < 1e-9 * ND[name] .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.212 seconds) .. _sphx_glr_download_api_gallery_plasma_waves_plot_04_debye_screening.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_04_debye_screening.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_04_debye_screening.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_04_debye_screening.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_