{"id":769,"date":"2025-07-20T03:24:18","date_gmt":"2025-07-20T03:24:18","guid":{"rendered":"https:\/\/synchrotron-light.net\/wp\/?page_id=769"},"modified":"2025-07-21T05:26:10","modified_gmt":"2025-07-21T05:26:10","slug":"chapter-9-synchrotron-light-in-the-cosmos","status":"publish","type":"page","link":"https:\/\/synchrotron-light.net\/wp\/?page_id=769","title":{"rendered":"Chapter 9: Synchrotron light in the cosmos"},"content":{"rendered":"

Below is a set of python codes associated with Chapter 9 of Daniele Pelliccia and David M. Paganin, “Synchrotron Light: A Physics Journey from Laboratory to Cosmos” (Oxford University Press, 2025).<\/p>\n

In order to run any of these python codes, you will need to include the following header.<\/p>\n

import numpy as np \r\nimport pandas as pd\r\nimport matplotlib.pyplot as plt\r\nfrom matplotlib.colors import LogNorm\r\nfrom matplotlib.ticker import FormatStrFormatter\r\nfrom matplotlib import ticker\r\nimport matplotlib as mpl\r\nfrom sys import stdout\r\nfrom mpl_toolkits.mplot3d import Axes3D\r\nfrom scipy import special\r\nfrom skimage import exposure\r\nfrom scipy import special\r\n\r\nsynchrotron_style = {\r\n    'font.family': 'FreeSerif',\r\n    'font.size': 30,\r\n    'axes.linewidth': 2.0,\r\n    'axes.edgecolor': 'black',\r\n    'grid.color': '#5b5b5b',\r\n    'grid.linewidth': 1.2,\r\n}<\/pre>\n
Spectrum of a single particle as a function of the pitch angle<\/h2>\n
See Fig. 9.6.<\/p>\n
e = 1.602e-19  #C\r\nm = 9.10938356e-31 #kg\r\nc = 3.0e8  #m\/s\r\nB = 10  #T\r\nepsilon_0 = 8.854e-12  #F\/m\r\n\r\n\r\ngamma = 1000\r\nnsteps = 50000\r\nnc = 1.5*gamma**3\r\nnmax = int(8*c)\r\nn = np.logspace(0, 11,nsteps)\r\nalpha_label = [\"$0.5 \\\\times \\\\pi\/2$\",\"$0.75 \\\\times \\\\pi\/2$\",\"$0.9 \\\\times \\\\pi\/2$\", \\\r\n               \"$0.95 \\\\times \\\\pi\/2$\",\"$0.99 \\\\times \\\\pi\/2$\",\"$ \\\\pi\/2$\"]\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax= plt.subplots(figsize=(16,5.9))\r\n    for i, alpha in enumerate([np.pi\/4, 0.75*np.pi\/2,0.9*np.pi\/2,0.95*np.pi\/2,0.99*np.pi\/2, np.pi\/2]):\r\n        theta = 0\r\n        \r\n        beta_perp = np.sqrt(1.0-1.0\/gamma**2)*np.sin(alpha)\r\n        beta_par = np.sqrt(1.0-1.0\/gamma**2)*np.cos(alpha)\r\n        \r\n        gamma_perp = 1.0\/np.sqrt(1-beta_perp**2)\r\n        v_perp = beta_perp*c\r\n        \r\n        # Argument of the Bessel functions\r\n        arg_j = (n*beta_perp*np.cos(theta) \/ (1 - beta_par*np.sin(theta)) )\r\n        \r\n        power_perp = ((e**4 * B**2 *n**2 * beta_perp**2)\/(8*np.pi**2 * epsilon_0 * m**2 * c * gamma**2*(1-beta_par*np.sin(theta))**4 ))*\\\r\n                     (special.jvp(n,arg_j)**2 + ( (np.sin(theta) - beta_par) \/ (beta_perp*np.cos(theta))  )**2 *special.jv(n,arg_j)**2 \/ beta_perp**2 )\r\n        \r\n        ax.loglog(n, power_perp, lw=4, label=\"$\\\\alpha =$\"+alpha_label[i])\r\n        ax.set_ylim(1e-21,1e-13)\r\n        ax.set_xlim(0.5,1e11)\r\n        ax.tick_params(axis='x', labelsize=26, pad=5)\r\n        ax.tick_params(axis='y', labelsize=26)\r\n\r\n        ax.set_xlabel('$n$', fontsize=32, labelpad=1)\r\nplt.show()<\/pre>\n
\"\"<\/p>\n
Effect of electron energy distribution on synchrotron-emission spectrum<\/h2>\n
Simulation of a collection of Gaussian spectral lines<\/h3>\n
See Fig. 9.8.<\/p>\n
e = 1.602e-19  #C\r\nm = 9.10938356e-31 #kg\r\nc = 3.0e8  #m\/s\r\nB = 10  #T\r\n\r\nsigma=1e14\r\ngg = np.logspace(1,4,200)\r\ndg = gg - np.roll(gg, 1)\r\ndg[0] = dg[1]\r\nomega = np.logspace(14,20,500000)\r\n\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, axs = plt.subplots(1,2, figsize=(16,8))\r\n    \r\n    for gamma in gg:\r\n\r\n        omega_crit = 3*e*B*gamma**2 \/ (2*m)\r\n\r\n        em_power = (1\/gamma**2.5)*np.exp(-(omega - omega_crit)**2 \/ (2*sigma**2))\r\n        axs[0].loglog(omega, em_power, lw=2, c='k')\r\n    \r\n    axs[0].loglog(omega, 3.5e15*omega**(-1.25), lw=4, c='r')\r\n    axs[0].set_ylim(1e-11,0.1)\r\n    axs[0].set_xlim(1e14,1e20)\r\n    axs[0].tick_params(axis='x', labelsize=26, pad = 5)\r\n    axs[0].tick_params(axis='y', labelsize=26)\r\n    axs[0].set_xlabel('$\\\\omega$ (Hz)', fontsize=28, labelpad=3)\r\n    \r\n    em_power = np.zeros_like(em_power)\r\n    for ii, gamma in enumerate(gg):\r\n\r\n        omega_crit = 3*e*B*gamma**2 \/ (2*m)\r\n\r\n        em_power = em_power + (1\/gamma**2.5)*np.exp(-(omega - omega_crit)**2 \/ (2*sigma**2))*dg[ii]\r\n    \r\n    axs[1].loglog(omega, em_power, lw=2, c='k')\r\n    \r\n    axs[1].loglog(omega, 7.5e7*omega**(-0.75), lw=4, c='b')\r\n    \r\n    axs[1].set_ylim(1e-11,0.1)\r\n    axs[1].set_xlim(1e14,1e20)\r\n    axs[1].tick_params(axis='x', labelsize=26, pad = 5)\r\n    axs[1].tick_params(axis='y', labelsize=26)\r\n    axs[1].set_xlabel('$\\\\omega$ (Hz)', fontsize=28, labelpad=3)\r\n\r\nplt.show()<\/pre>\n
\"\"<\/p>\n
Spectrum emitted by an ensemble of electrons with power-law distribution of energies<\/h3>\n
See Fig. 9.9.<\/p>\n
Phi_pos = -0.1\r\nlambda_u = 0.01\r\nk_u = 2*np.pi\/lambda_u\r\nc = 3.0e8\r\nomega_u = 2*k_u*c\r\nprint(1e-11*omega_u)\r\nN = 10\r\nt = np.linspace(0,N*2*np.pi\/omega_u,1000000)\r\n\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax = plt.subplots(figsize=(12,5.8))\r\n    ax.plot(1e9*t, -np.ones_like(t)*np.sin(Phi_pos), lw = 4, label=\"$\\\\sin \\\\Phi^{+}$\")\r\n    ax.plot(1e9*t, -np.sin(Phi_pos-omega_u*t), lw = 4, label=\"$\\\\sin \\\\Phi^{-}$\")\r\n    ax.tick_params(axis='y', labelsize=26)\r\n    ax.tick_params(axis='x', labelsize=26)\r\n    ax.set_ylabel('Amplitude', fontsize=26)\r\n    ax.set_xlabel('$t$ (ns)', fontsize=26)\r\n    \r\n    ax.legend(fontsize=28, frameon=True, framealpha=1, handlelength=1)\r\nplt.show()<\/pre>\n
\"\"<\/p>\n
Qualitative spectral distribution of the synchrotron light emitted by an ensemble of charged particles<\/h2>\n
Qualitative synchrotron radiation spectrum, including self-absorption. See Fig. 9.12.<\/p>\n
omega1 = np.logspace(14,30,50000)\r\nomega = np.logspace(14,36,50000)\r\nom2 = 1.0e32\r\nkk = 0.2e-24\r\n\r\nfunc1 = 1e60\/(1e-23*omega1**3.252 + 1e60*omega1**(-0.002))\r\n\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax= plt.subplots(figsize=(16,7.8))\r\n    \r\n    ax.loglog(omega1, 3.5e-45*omega1**(2.5)*(func1) , lw=4, c='k')\r\n    ax.loglog(omega[35000:], tot_power[35000:], lw=4, c='k')\r\n\r\n    ax.loglog(omega, 3.5e38*omega**(-0.75), lw=2, c='b', ls='--')  \r\n    ax.loglog(omega, 3.5e-45*omega**(2.5), lw=2, c='g', ls='--')  \r\n    ax.loglog(omega, tot_power, lw=2, c='m', ls='--')\r\n    ax.set_ylim(1e1,1e25)\r\n    ax.set_xlim(1e19,1e34)\r\n    ax.tick_params(axis='x', labelsize=26, pad = 5)\r\n    ax.tick_params(axis='y', labelsize=26)\r\n    ax.set_ylim(5e1,1e25)\r\n    ax.set_xlabel('$\\\\omega$ (Hz)', fontsize=28, labelpad=3)\r\n\r\nplt.show()<\/pre>\n
\"\"<\/p>\n
Synchrotron cooling time<\/h2>\n
See Fig. 9.14.<\/p>\n
e = 1.602e-19  #C\r\nm = 9.10938356e-31 #kg\r\nc = 3.0e8  #m\/s\r\nB = 5e-8  #T\r\nepsilon_0 = 8.854e-12  #F\/m\r\nhbar = 6.5821e-16 # eV s\r\n\r\ngamma = np.logspace(1,13,2000)\r\nbeta = np.sqrt(1.0-1.0\/gamma**2)\r\nomega_c = 3*e*B*gamma**2 \/ (2*m)\r\nen_c = hbar*omega_c\r\n\r\nP_avg = (e**2 *B*gamma*beta)**2 \/ (9*np.pi*epsilon_0*m**2*c)\r\ntau_s = gamma*m*c**2 \/ P_avg\r\ntau_y = tau_s\/31556909\r\n\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax= plt.subplots(figsize=(8,7.8))\r\n    \r\n    ax.loglog(en_c, tau_s, lw=4)\r\n    \r\n    ax.tick_params(axis='x', labelsize=29, pad=5)\r\n    ax.tick_params(axis='y', labelsize=29)\r\n    ax.set_xlabel('Photon critical energy $\\\\mathcal{E}_{c}$ (eV)', fontsize=28, labelpad=3)\r\n    ax.set_ylabel('Cooling time (s)', fontsize=30, labelpad=3)\r\n    \r\n    ax2 = ax.twinx()\r\n    ax2.loglog(en_c, tau_y, lw=4)\r\n    ax2.tick_params(axis='y', labelsize=29)\r\n    ax2.set_ylabel('Cooling time (yr)', fontsize=30, labelpad=3)\r\n    ax.grid(False,axis='y')\r\nplt.show()<\/pre>\n
\"\"<\/p>\n
Synchrotron cooling: numerical calculation of the time-dependent Lorentz factor<\/h2>\n
See Fig. 9.15.<\/p>\n
e = 1.602e-19  #C\r\nm = 9.10938356e-31 #kg\r\nc = 3.0e8  #m\/s\r\nB = 5e-8  #T\r\nepsilon_0 = 8.854e-12  #F\/m\r\n\r\n# Define the constant\r\nA = (e**2 *B)**2 \/ (9*np.pi*epsilon_0*m**3*c**3)\r\n\r\nt = np.logspace(5,17.5,100000)\r\n\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax= plt.subplots(figsize=(12,7.8))\r\n    \r\n    gamma_label = ['$10^{2}$', '$10^{4}$', '$10^{6}$', '$10^{8}$']\r\n    for k, gamma_0 in enumerate([100,10000,1000000,100000000]):\r\n        gamma = np.zeros_like(t)\r\n        gamma[0] = gamma_0\r\n        for i,j in enumerate(t[:-1]):\r\n            gamma[i+1] = gamma[i] - A*(gamma[i]**2 - 1)*(t[i+1]-t[i])\r\n        \r\n        ax.loglog(t\/31556909, gamma, lw=4, label='$\\\\gamma(0) =$'+gamma_label[k])\r\n\r\n    ax.tick_params(axis='x', labelsize=28, pad=5)\r\n    ax.tick_params(axis='y', labelsize=28)\r\n    ax.set_ylabel('$\\\\gamma(t)$', fontsize=28, labelpad=3)\r\n    ax.set_xlabel('Time (yr)', fontsize=28, labelpad=3)\r\n    plt.yticks(np.logspace(0, 8, 9))\r\n        \r\n    plt.legend(fontsize = 28, frameon=True, framealpha = 1, handlelength = 1)    \r\n\r\nplt.show()<\/pre>\n
\"\"<\/p>\n
Synchrotron cooling: comparison of numerical and closed-form solution<\/h2>\n
See Fig. 9.16.<\/p>\n
e = 1.602e-19  #C\r\nm = 9.10938356e-31 #kg\r\nc = 3.0e8  #m\/s\r\nB = 5e-5  #T\r\nepsilon_0 = 8.854e-12  #F\/m\r\n\r\nA = (e**2 *B)**2 \/ (9*np.pi*epsilon_0*m**3*c**3)\r\n\r\nt = np.logspace(4,10.5,500000)\r\ngamma = np.zeros_like(t)\r\ngamma[0] = 1e4\r\nfor i,j in enumerate(t[:-1]):\r\n    gamma[i+1] = gamma[i] - A*(gamma[i]**2 - 1)*(t[i+1]-t[i])\r\n    \r\nA_cal = (gamma[0] - 1)\/(gamma[0] + 1)\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax= plt.subplots(figsize=(8,8))\r\n    ax.loglog(t, gamma, lw=6, label=\"Numerical solution\")\r\n    ax.loglog(t,gamma[0]\/(1+A*gamma[0]*t), lw=4, label=\"$\\\\gamma(0) \\gg 1$ solution\")\r\n    \r\n    ax.loglog(t, (1 + A_cal*np.exp(-2*A*t)) \/ (1 - A_cal*np.exp(-2*A*t)), 'r',lw=3, label=\"Analytical solution\" )\r\n    \r\n    ax.tick_params(axis='x', labelsize=30, pad=7)\r\n    ax.tick_params(axis='y', labelsize=30)\r\n    ax.set_ylabel('$\\\\gamma(t)$', fontsize=34, labelpad=3)\r\n    ax.set_xlabel('Time (yr)', fontsize=34, labelpad=3)\r\n\r\n    \r\n    plt.legend(fontsize = 28, frameon=True, framealpha = 1, handlelength = 1)  \r\n    plt.tight_layout()\r\nplt.show()<\/pre>\n
\"\"<\/p>\n
Fermi model for the acceleration of cosmic rays<\/h2>\n
See Fig. 9.21. Note: the code implements a Monte Carlo simulation, hence the plot will generally be different from the one reproduced in the figure, depending on the initialisation of the random number generator.<\/p>\n
def fermi_acceleration(n_particles,n_epochs, V, gamma_init, p_esc):\r\n    # Start with a monochromatic set of particles\r\n    gamma_particles = gamma_init*np.ones(n_particles)\r\n    \r\n    for i in range(n_epochs):\r\n        # Particles can escape the region with probability p_esc \r\n        prob1 = np.random.rand(gamma_particles.shape[0])\r\n        mask1 = prob1 > p_esc\r\n        gamma_particles=gamma_particles[mask1]\r\n\r\n        # The particles that are not lost, can either be accelerated or decelerated\r\n        prob2 = np.random.rand(gamma_particles.shape[0])\r\n        beta_particles = np.sqrt(1.0-1\/gamma_particles)\r\n        p_plus = (beta_particles + beta_mag)\/(2*beta_particles)\r\n        \r\n        # particles that gain energy\r\n        mask_plus = prob2 <= p_plus \r\n        gamma_particles[mask_plus] =  gamma_particles[mask_plus] + \\\r\n                                      2*gamma_particles[mask_plus]*beta_particles[mask_plus]*V\/c\r\n        # particles that lose energy    \r\n        mask_minus = prob2 > p_plus\r\n        gamma_particles[mask_minus] =  gamma_particles[mask_minus] - \\\r\n                                      2*gamma_particles[mask_minus]*beta_particles[mask_minus]*V\/c        \r\n    return gamma_particles \r\n\r\nc = 3e8\r\nbeta_mag = 1e-2\r\nV = beta_mag * c\r\nrest_energy = 0.511e6 \r\n\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax= plt.subplots(figsize=(8,8))\r\n    h, bins = np.histogram(rest_energy*gamma_final, bins=100)\r\n    ax.loglog(bins[1:],h, drawstyle='steps', lw=4)\r\n    ax.set_xlabel('Particle energy (eV)', rotation=0, fontsize=30)\r\n    ax.set_ylabel('Number of particles', fontsize=30, labelpad = 7)\r\n    ax.tick_params(axis='x', labelsize=26, pad=7)\r\n    ax.tick_params(axis='y', labelsize=26)\r\n\r\nplt.show()\r\n\r\ngamma_final = fermi_acceleration(200000,5000, V=V, gamma_init=5, p_esc=1e-6)\r\n<\/pre>\n
\"\"<\/p>\n
Energy equipartition<\/h2>\n
See Fig. 9.22.<\/p>\n
eta = 40 # As measured in cosmic rays\r\nB = np.linspace(0.1,100,1000)\r\nU = (1+eta)*B**(-1.5) + B**2\r\n\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax= plt.subplots(figsize=(8,8))\r\n    ax.loglog(B,U, 'k',lw=4)\r\n    ax.loglog(B, (1+eta)*B**(-1.5), 'r--', lw=3)\r\n    ax.loglog(B, B**2, 'b--', lw=3)\r\n    \r\n    ax.set_ylim(1,1e4)\r\n    ax.tick_params(axis='x', labelsize=26, pad = 7)\r\n    ax.tick_params(axis='y', labelsize=26)\r\n    ax.set_xlabel('B (T)', fontsize=28, labelpad=3)\r\nplt.show()<\/pre>\n
\"\"<\/p>\n
All-sky radio survey<\/h2>\n
The data plotted in the all-sky radio survey maps of Figs. 9.32 and 9.33 can be downloaded at the following pages:<\/p>\n