Galactic Dynamo

• Home
• Task 1
• Task 2
• Task 3
• References

3.1 - Problem Statement

To solve the mean field induction equation using vector potential $A$ in the Weyl gauge ($\phi=0$) (Brandenburg 2003, Sec. 5) and compare with the solutions for the equations in $B$.

What are the differences and which method is most efficient?

---

---

3.2 - Theory

In the previous chapter, we studied how the galactic magnetic field varies with time in the $z$-direction. For that we covered the diffusion equation, with the addition of the $\alpha-\Omega$ term in our equations. We studied about the $\alpha-\Omega$ dynamo, and how the galactic magnetic field varies with time in that case for dynamo numbers lower and greater than the critical dynamo number, using Crank-Nicolson numerical method and some example seed fields, and finally calculated the decay rate and then the critical dynamo number was calculated.

In this section, we will solve the induction equation in terms of the magnetic potential $A$ and try to obtain the same type of solutions.

Let's first derive the equations and set-up the floor in terms of the potential.

$$ \dfrac{\partial \mathbf{\bar{B}}}{\partial t} = \nabla \times \left[ \mathbf{\bar{V}} \times \mathbf{\bar{B}} + \mathcal{E} - \eta \left( \nabla \times \mathbf{\bar{B}} \right) \right] $$$$ \text{where} \qquad \mathcal{E} = \left( \alpha \mathbf{\bar{B}} \right) - \eta_t \left( \nabla \times \mathbf{\bar{B}} \right) $$

We will now implement the same equation using the vector potential $\mathbf{\bar{A}}$, using the relation $\mathbf{\bar{B}} = \nabla \times \mathbf{\bar{A}}$ with the Weyl gauge $\phi=0$. So we have

$$ \bar{A}= \bar{V} \times \bar{B}+\alpha \bar{B}-\eta \nabla \times \bar{B} $$$$ \bar{A}= \left(V_\phi B_z-B_\phi V_z+\eta \frac{\partial B_\phi}{\partial z}+\alpha B_r\right) \hat{r} +\left(V_z B_r-V_r B_z+\eta \frac{\partial B_z}{\partial r}-\eta \frac{\partial B_r}{\partial z}+\alpha B_\phi\right) \hat{\phi} +\left(V_r B_\phi-V_\phi B_r-\eta \frac{\partial B_\phi}{\partial r}-\eta \frac{B_\phi}{r}+\alpha B_z\right) \hat{z} $$

$$ \begin{aligned} \frac{\partial B_\phi}{\partial t} = & \: \frac{\partial A_r}{\partial z}-\frac{\partial A_z}{\partial r} \\ = & \: \frac{\partial}{\partial z}\left(V_\phi B_z-V_z B_\phi+\eta \frac{\partial B_\phi}{\partial z}+\alpha B_x\right) -\frac{\partial}{\partial r}\left(V_r B_\phi-V_\phi B_r-\eta \frac{\partial B_\phi}{\partial r}-\eta \frac{B_\phi}{r}+\alpha B_z\right) \end{aligned} $$

We used the change of variables $$ B_r = - \dfrac{1}{r} \dfrac{\partial \psi}{\partial z} $$ $$ B_\phi=\dfrac{T}{r} $$ $$ B_z= \dfrac{1}{r} \dfrac{\partial \psi}{\partial r} $$ $$ V_\phi=r \Omega $$ $$ S_z=r \dfrac{\partial \Omega}{\partial z} $$ $$ q=- \dfrac{\partial \ln \Omega}{\partial \ln r} $$ $$ \Lambda^{ \pm}=\dfrac{\partial^2}{\partial r^2} \pm \dfrac{1}{r} \dfrac{\partial}{\partial r}+\dfrac{\partial^2}{\partial z^2} $$

Multiplying both sides by $r$ and changing the variable, we get $$ \begin{aligned} \frac{\partial T}{\partial t}= & \: r \frac{\partial}{\partial z}\left(\Omega \frac{\partial \psi}{\partial r}\right)-\frac{\partial}{\partial z}\left(V_z T\right)+\frac{\partial}{\partial z}\left(\eta \frac{\partial T}{\partial z}\right)-\frac{\partial}{\partial z}\left(\alpha \frac{\partial \psi}{\partial z}\right) \\ & -r \frac{\partial}{\partial r}\left(\frac{V_r T}{r}\right)-r \frac{\partial}{\partial r}\left(\Omega \frac{\partial \Psi}{\partial z}\right)+r \frac{\partial}{\partial r}\left(\eta \frac{\partial}{\partial r}\left(\frac{T}{r}\right)\right. \\ & +r \frac{\partial}{\partial r}\left(\eta \frac{T}{r^2}\right)-r \frac{\partial}{\partial r}\left(\alpha \frac{1}{r} \frac{\partial \psi}{\partial r}\right) \end{aligned} $$

Upon simplification, we get $$ \begin{aligned} \frac{\partial T}{\partial t}= & -r \frac{\partial}{\partial r}\left(\frac{V_r T}{r}\right)-\frac{\partial}{\partial z}\left(V_z T\right)+q \Omega \frac{\partial \psi}{\partial z}+S_z \frac{\partial \psi}{\partial r}-\alpha \Lambda^{-} \psi \\ & -\frac{\partial \alpha}{\partial r} \frac{\partial \psi}{\partial r}-\frac{\partial \alpha}{\partial z} \frac{\partial \psi}{\partial z}+\eta \Lambda^{-} T \end{aligned} $$

After removing the outflow terms and $\alpha^2$-terms, we get $$ \frac{\partial T}{\partial t}=q \Omega \frac{\partial \psi}{\partial z}+\eta \Lambda^{-} T $$

Similarly we solve for $B_r$ $$ \frac{\partial B_r}{\partial t}=\frac{1}{r}\left(\frac{\partial A_z}{\partial \phi}-\frac{\partial}{\partial z}\left(r A_\phi\right)\right) $$

Taking azimuthal symmetry, $$ \frac{\partial B_r}{\partial t}=\frac{1}{r}\left(-r \frac{\partial}{\partial z}\left(A_\phi\right)\right)=-\frac{\partial A_\phi}{\partial z} $$

Now, putting $B_r=-\dfrac{1}{r} \dfrac{\partial \psi}{\partial z}$ gives $$ \frac{\partial \psi}{\partial t}=r A_\phi $$

$$ \begin{aligned} \frac{\partial \psi}{\partial t} & =r\left(V_z B_r-V_r B_z+\eta \frac{\partial B_z}{\partial r}-\eta \frac{\partial B_r}{\partial z}+\alpha B_\phi\right) \\ & =-r V_z \frac{1}{r} \frac{\partial \psi}{\partial z}-V_r \frac{\partial \psi}{\partial r}+r \eta \frac{\partial}{\partial z}\left(\frac{1}{r} \frac{\partial \psi}{\partial r}\right)+r \eta \frac{\partial}{\partial r}\left(\frac{1}{r} \frac{\partial \psi}{\partial z}\right)+\alpha T \\ & =-V_z \frac{\partial \psi}{\partial z}-V_r \frac{\partial \psi}{\partial r}+\alpha T+\eta\left(\frac{\partial^2}{\partial r^2}-\frac{1}{r} \frac{\partial}{\partial r}+\frac{\partial^2}{\partial z^2}\right) \psi \end{aligned} $$

Finally after simplification, we get $$ \frac{\partial \psi}{\partial t}=-V_z \frac{\partial \psi}{\partial z}-V_r \frac{\partial \psi}{\partial r}+\alpha T+\eta \Lambda^{-} \psi $$

Since we solving in $z$-coordinate, we will ignore all $r$-derivatives $\left(\dfrac{\partial}{\partial r}=0\right)$. This gives the final equations after expending $\Lambda^{-}$ as.

$$ \boxed{ \frac{\partial \psi}{\partial t} = \alpha T + \eta_T \frac{\partial^2 \psi}{\partial z^2} } \qquad \qquad \text{and} \qquad \qquad \boxed{ \frac{\partial T}{\partial t} = q \Omega \frac{\partial \psi}{\partial z} + \eta_T \frac{\partial^2 T}{\partial z^2} } $$

---

3.3 - Numerical Solutions

3.3.1 - Solver: Crank-Nicholson (modified for coupled equations)

We had used the Crank-Nicholson method for solving the 1D diffusion equation in the previous chapter, however, things get complicated when we mpove to coupled diffusion equations, like in this case. Let us understand the discretization scheme for coupled diffusion equations.

$$ \frac{\psi^{j+1}_{i} - \psi^{j}_{i}}{dt} = \alpha_i \left( \frac{T^{j+1}_{i} + T^{j}_{i}}{2} \right) + \dfrac{\eta_T}{2} \: \left( \dfrac{\psi^{j+1}_{i+1} - 2 \psi^{j+1}_{i} + \psi^{j+1}_{i-1}}{dz^2} + \dfrac{\psi^{j}_{i+1} - 2 \psi^{j}_{i} + \psi^{j}_{i-1}}{dz^2} \right) $$

and $$ \frac{T^{j+1}_{i} - T^{j}_{i}}{dt} = \dfrac{q \Omega}{2} \left( \frac{\psi^{j+1}_{i+1} - \psi^{j+1}_{i}}{dz} + \frac{\psi^{j}_{i+1} - \psi^{j}_{i}}{dz} \right) + \dfrac{\eta_T}{2} \: \left( \dfrac{T^{j+1}_{i+1} - 2 T^{j+1}_{i} + T^{j+1}_{i-1}}{dz^2} + \dfrac{T^{j}_{i+1} - 2 T^{j}_{i} + T^{j}_{i-1}}{dz^2} \right) $$

Putting $\mu = \dfrac{dt}{2 \: dz}$ and $\nu = \dfrac{\eta_T \: dt}{2 \: dz^2}$, we separate the present time-step $(j+1)$ and the past time-step $(j)$ as

$$ \left(1+2\nu \right) \psi^{j+1}_{i} - \nu \psi^{j+1}_{i+1} - \nu \psi^{j+1}_{i-1} - dz \mu \alpha_i T^{j+1}_{i} = \left(1-2\nu \right) \psi^{j}_{i} + \nu \psi^{j}_{i+1} + \nu \psi^{j}_{i-1} + dz \mu \alpha_i T^{j}_{i} $$

and $$ \left(1+2\nu \right) T^{j+1}_{i} - \nu T^{j+1}_{i+1} - \nu T^{j+1}_{i-1} - \mu q \Omega \psi^{j+1}_{i+1} + \mu q \Omega \psi^{j+1}_{i} = \left(1-2\nu \right) T^{j}_{i} + \nu \psi^{j}_{i+1} + \nu \psi^{j}_{i-1} + \mu q \Omega \psi^{j}_{i+1} - \mu q \Omega \psi^{j}_{i} $$

Let $U = \left[ \psi \:\: T\right]^T$ This whole coupled equation can be simplified to single variable $U$ as

$$ \left[\begin{array}{cc} (1+2 \nu) & - dz \mu \alpha_i \\ \mu q \Omega & (1+2 \nu) \end{array}\right] U^{j+1}_{i} + \left[\begin{array}{cc}-\nu & 0 \\ -\mu q \Omega & -\nu \end{array}\right] U^{j+1}_{i+1} + \left[\begin{array}{cc}-\nu & 0 \\ 0 & -\nu \end{array}\right] U^{j+1}_{i-1} = \left[\begin{array}{cc}(1-2 \nu) & dz \mu \alpha_i \\ -\mu q \Omega & (1-2 \nu) \end{array}\right] U^{j}_{i} + \left[\begin{array}{cc}\nu & 0 \\ \mu q \Omega & \nu \end{array}\right] U^{j}_{i+1} + \left[\begin{array}{cc}\nu & 0 \\ 0 & \nu \end{array}\right] U^{j}_{i-1} $$

Renaming these $2 \times 2$ matrices $E$, $F$, $G$, $H$, $I$ and $J$ from left to right respectively, we can write it as $$ E U^{j+1}_{i} + F U^{j+1}_{i+1} + G U^{j+1}_{i-1} = H U^{j}_{i} + I U^{j}_{i+1} + J U^{j}_{i-1} $$

Again we can now form the final matrix, however, this time, it is a bit different. $$ \widetilde{M}U^{j+1} = \widetilde{N}U^{j} $$ Here, the final matrices $\widetilde{M}$ and $\widetilde{N}$ are the tensor product of the above 6 matrices and the M, N matrices from the previous chapter. The final matrix would look something like this

$$ \widetilde{M} = \left[\begin{array}{cc} \left[\begin{array}{ccccc}1+2 \nu & -\nu & 0 & \cdots & 0 \\ -\nu & 1+2 \nu & -\nu & \ddots & \vdots \\ 0 & -\nu & 1+2 \nu & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -\nu \\ 0 & \cdots & 0 & -\nu & 1+2 \nu\end{array}\right] & \left[\begin{array}{ccccc} - dz \mu \alpha_i & 0 & 0 & \cdots & 0 \\ 0 & - dz \mu \alpha_i & 0 & \ddots & \vdots \\ 0 & 0 & - dz \mu \alpha_i & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & 0 & - dz \mu \alpha_i \end{array}\right] \\ & \\ \left[\begin{array}{ccccc} \mu q \Omega & -\mu q \Omega & 0 & \cdots & 0 \\ 0 & \mu q \Omega & -\mu q \Omega & \ddots & \vdots \\ 0 & 0 & \mu q \Omega & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -\mu q \Omega \\ 0 & \cdots & 0 & 0 & \mu q \Omega \end{array}\right] & \left[\begin{array}{ccccc}1+2 \nu & -\nu & 0 & \cdots & 0 \\ -\nu & 1+2 \nu & -\nu & \ddots & \vdots \\ 0 & -\nu & 1+2 \nu & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -\nu \\ 0 & \cdots & 0 & -\nu & 1+2 \nu\end{array}\right] \end{array}\right] $$

and

$$ \widetilde{N} = \left[\begin{array}{cc} \left[\begin{array}{ccccc}1-2 \nu & \nu & 0 & \cdots & 0 \\ \nu & 1-2 \nu & \nu & \ddots & \vdots \\ 0 & \nu & 1-2 \nu & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \nu \\ 0 & \cdots & 0 & \nu & 1-2 \nu\end{array}\right] & \left[\begin{array}{ccccc} dz \mu \alpha_i & 0 & 0 & \cdots & 0 \\ 0 & dz \mu \alpha_i & 0 & \ddots & \vdots \\ 0 & 0 & dz \mu \alpha_i & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & 0 & dz \mu \alpha_i \end{array}\right] \\ & \\ \left[\begin{array}{ccccc} -\mu q \Omega & \mu q \Omega & 0 & \cdots & 0 \\ 0 & -\mu q \Omega & \mu q \Omega & \ddots & \vdots \\ 0 & 0 & -\mu q \Omega & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \mu q \Omega \\ 0 & \cdots & 0 & 0 & -\mu q \Omega \end{array}\right] & \left[\begin{array}{ccccc}1-2 \nu & \nu & 0 & \cdots & 0 \\ \nu & 1-2 \nu & \nu & \ddots & \vdots \\ 0 & \nu & 1-2 \nu & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \nu \\ 0 & \cdots & 0 & \nu & 1-2 \nu\end{array}\right] \end{array}\right] $$

The final answer i.e., $U^{j+1}$ is obtained as $$ U^{j+1} = \widetilde{M}^{-1}\widetilde{N}U^{j} $$

3.3.2 - Grid and boundary conditions

We will use the same grid structure and boumndary conditions that we used in chapter 1.

---

---

THE CODE

For the full code with al functions and plotting scripts, please visit my github repository.

Now let us get into some examples to understand how the magnetic field varies with different dynamo numbers.

from my_code import *
from plotting import *

Example 1

We have already defined all the parameters in section 2.2.2.

For this example, we take an $\alpha$ for which the field is decaying.

def init_cond_S(x):
    return -(1-x**2)

def init_cond_T(x):
    return (1-x**2)

def source_term(x, t):
    return 0

z = np.linspace(-1, 1, 101)
title_1 = r'Parabolic seed field '+'\n'+r'$ y = x^2-1$'
title_2 = r'Parabolic seed field '+'\n'+r'$ y = 1-x^2$'
global_title = 'INITIAL CONDITIONS FOR MAGNETIC FIELD POTENTIAL COMPONENTS'

plt.figure(figsize=(11, 3.5))
plt.subplot(121)
plt.plot(z, init_cond_S(z))
plt.xlabel(r'$z$ (normalized to 100 pc)')
plt.ylabel(r'$B_r$ (in $\mu G$/100pc)')
plt.title(title_1)

plt.subplot(122)
plt.plot(z, init_cond_T(z))
plt.xlabel(r'$z$ (normalized to 100 pc)')
plt.ylabel(r'$B_{\phi}$ (in $\mu G$/100pc)')
plt.title(title_2)

plt.suptitle(global_title)
plt.tight_layout(pad=1)
plt.show()

# Constants and parameters
h = 5                               # 100pc
eta_T = 0.7#3.48e-2/h**2                     # magnetic diffusivity in (100pc)^2/Myr
Omega = 40*MYR*KM/(1000*PC)        # angular velocity, converted from km/s/kpc to 1/Myr
q = 0.98                            # shear parameter
t_max = 10                          # total simulation time
z_min = -1.0                        # minimum thickness of the disc
z_max = 1.0                         # thickness of the disc
dt = t_max/200                      # time step
dz = 0.01                           # spatial step in z direction
alpha_0 = 1/h#10*1e3*MYR/(100*PC)      # alpha effect, converted from km/s to 100pc/Myr

print('eta', eta_T, '\t',  'alpha',alpha_0, '\t', 'omega', Omega)
print('Rw', -q*Omega*h**2/eta_T)
print('Ra', alpha_0*h/eta_T)
print('Dynamo number ', -alpha_0*q*Omega*h**3/eta_T**2)

# Spatial grid
z = np.linspace(z_min, z_max, int((z_max - z_min) / dz) + 1)
t = np.linspace(0, t_max, int(t_max / dt) + 1)

# Coefficients for the matrix A and B
nu = eta_T*dt/(2*dz**2)
mu = dt/(2*dz)
alpha = alpha_0*np.sin(np.pi*z/2)
dalpha_dt = np.gradient(alpha, z)

A = mod_matrix(len(z), 1+2*nu, dz*mu, q*Omega*mu, 1+2*nu, -nu, 0, -q*Omega*mu, -nu, -nu, 0, 0, -nu, alpha, dalpha_dt, dz, dt)
B = mod_matrix(len(z), 1-2*nu, -dz*mu, -q*Omega*mu, 1-2*nu, nu, 0, q*Omega*mu, nu, nu, 0, 0, nu, alpha, dalpha_dt, dz, dt)

# Solve the diffusion equation in radial direction
solution = crank_nicolson_mod(len(z), len(t), init_cond_S(z), init_cond_T(z), A, B)

S = solution[:len(z), :]
T = solution[len(z):, :]

eta 0.7 	 alpha 0.2 	 omega 0.04087621516526248
Rw -1.430667530784187
Ra 1.4285714285714286
Dynamo number  -2.0438107582631244

"""
Plot the solution in both 1D and Heatmap format.
"""
def plot_diff_A(time_grid, spatial_grid, solution_r, solution_phi):

    # Create 2D plots
    plt.figure(figsize=(12, 8))
    plt.subplot(2, 2, 1)
    plt.suptitle(r'Diffusion of Magnetic potential $\psi$ and $T$')

    for i in (range(0, len(time_grid), int(len(time_grid)/5))):
        plt.plot(spatial_grid, solution_r[:, i], label=f'time = {time_grid[i]:.1f} Myrs')
    plt.xlabel(r'$z$ (normalized to 100 pc)')
    plt.ylabel(r'Magnetic potential $\psi$ (in $\mu G$/100pc)')
    plt.title(r'Diffusion of $\psi$')
    plt.grid()
    plt.legend(loc='upper right')

    # Create imshow plot
    plt.subplot(2, 2, 2)
    plt.contourf(*np.meshgrid(spatial_grid, time_grid), solution_r.T, 50, cmap='Spectral_r')
    plt.colorbar(label=r'Magnetic potential $\psi$ (in $\mu G$/100pc)')
    plt.title(r'Diffusion of $\psi$')
    plt.xlabel(r'$z$ (normalized to 100 pc)')
    plt.ylabel(r'Time (Myr)')


    # Create 2D plots
    plt.subplot(2, 2, 3)
    for i in (range(0, len(time_grid), int(len(time_grid)/5))):
        plt.plot(spatial_grid, solution_phi[:, i], label=f'time = {time_grid[i]:.1f} Myrs')
    plt.xlabel(r'$z$ (normalized to 100 pc)')
    plt.ylabel(r'Magnetic potential $T$ (in $\mu G$/100pc)')
    plt.title(r'Diffusion of $T$')
    plt.grid()
    plt.legend(loc='upper right')

    # Create imshow plot
    plt.subplot(2, 2, 4)
    plt.contourf(*np.meshgrid(spatial_grid, time_grid), solution_phi.T, 50, cmap='Spectral_r')
    plt.colorbar(label=r'Magnetic potential $T$ (in $\mu G$/100pc)')
    plt.title(r'Diffusion of $T$')
    plt.xlabel(r'$z$ (normalized to 100 pc)')
    plt.ylabel('Time (Myr)')
    plt.tight_layout(pad=3)

# Plot the solution in imshow
plot_diff_A(t, z, S, T)
plt.show()

Here is an animation showing the evolution of $\psi$ and $T$ with time.

One can see that in this example, the fields are decaying with time. We will calculate the decay rate at the midplane, i.e., $z=0$. To calculate the decay rate, we will find how the maxima of the field at the midplane varies with time.

# Constants and parameters
h = 5                               # 100pc
eta_T = 0.7#3.48e-2/h**2                     # magnetic diffusivity in (100pc)^2/Myr
Omega = 40*MYR*KM/(1000*PC)        # angular velocity, converted from km/s/kpc to 1/Myr
q = 0.98                            # shear parameter
t_max = 10                          # total simulation time
z_min = -1.0                        # minimum thickness of the disc
z_max = 1.0                         # thickness of the disc
dt = t_max/200                      # time step
dz = 0.01                           # spatial step in z direction
alpha_0 = 1/h#10*1e3*MYR/(100*PC)      # alpha effect, converted from km/s to 100pc/Myr

print('eta', eta_T, '\t',  'alpha',alpha_0, '\t', 'omega', Omega)
print('Rw', -q*Omega*h**2/eta_T)
print('Ra', alpha_0*h/eta_T)
print('Dynamo number ', -alpha_0*q*Omega*h**3/eta_T**2)

# Spatial grid
z = np.linspace(z_min, z_max, int((z_max - z_min) / dz) + 1)
t = np.linspace(0, t_max, int(t_max / dt) + 1)

# Coefficients for the matrix A and B
nu = eta_T*dt/(2*dz**2)
mu = dt/(2*dz)
alpha = alpha_0*np.sin(np.pi*z/2)
dalpha_dt = np.gradient(alpha, z)

A = mod_matrix(len(z), 1+2*nu, dz*mu, q*Omega*mu, 1+2*nu, -nu, 0, -q*Omega*mu, -nu, -nu, 0, 0, -nu, alpha, dalpha_dt, dz, dt)
B = mod_matrix(len(z), 1-2*nu, -dz*mu, -q*Omega*mu, 1-2*nu, nu, 0, q*Omega*mu, nu, nu, 0, 0, nu, alpha, dalpha_dt, dz, dt)

# Solve the diffusion equation in radial direction
solution = crank_nicolson_mod(len(z), len(t), init_cond_S(z), init_cond_T(z), A, B)

S = solution[:len(z), :]
T = solution[len(z):, :]

eta 0.7 	 alpha 0.2 	 omega 0.04087621516526248
Rw -1.430667530784187
Ra 1.4285714285714286
Dynamo number  -2.0438107582631244

S_mid = S[int(len(z)*3/4), :]
T_mid = T[int(len(z)*3/4), :]

A_mid = np.sqrt(S_mid**2 + T_mid**2)

plt.figure(figsize=(6, 4))
plt.plot(t, A_mid)
plt.plot(t[int(len(t)*3/4):], A_mid[int(len(t)*3/4):], 'r:', linewidth=4)
plt.xlabel('Time (in Myrs)')
plt.ylabel('Magnetic potential (in $\mu G$/100pc)')
plt.yscale('log')
plt.title('Magnetic potential at the maximum amplitude point')
plt.show()

decay_rate = get_decay_rate(t, A_mid)
if decay_rate < 0:
    print(r"decay rate at the mid plane:", format(-decay_rate, ".3e"))
else:
    print(r"growth rate at the mid plane:", format(decay_rate, ".3e"))

print('Dynamo number ', -alpha_0*q*Omega*h**3/eta_T**2)

decay rate at the mid plane: 3.458e-01
Dynamo number  -2.0438107582631244

Example 2

For this example, we take an $\alpha$ for which the field is decaying.

def init_cond_S(x):
    return -(1-x**2)

def init_cond_T(x):
    return (1-x**2)

def source_term(x, t):
    return 0

z = np.linspace(-1, 1, 101)
title_1 = r'Parabolic seed field '+'\n'+r'$ y = x^2-1$'
title_2 = r'Parabolic seed field '+'\n'+r'$ y = 1-x^2$'
global_title = 'INITIAL CONDITIONS FOR MAGNETIC FIELD POTENTIAL COMPONENTS'

plt.figure(figsize=(11, 3.5))
plt.subplot(121)
plt.plot(z, init_cond_S(z))
plt.xlabel(r'$z$ (normalized to 100 pc)')
plt.ylabel(r'$B_r$ (in $\mu G$/100pc)')
plt.title(title_1)

plt.subplot(122)
plt.plot(z, init_cond_T(z))
plt.xlabel(r'$z$ (normalized to 100 pc)')
plt.ylabel(r'$B_{\phi}$ (in $\mu G$/100pc)')
plt.title(title_2)

plt.suptitle(global_title)
plt.tight_layout(pad=1)
plt.show()

# Constants and parameters
h = 5                               # 100pc
eta_T = 0.7#3.48e-2/h**2                     # magnetic diffusivity in (100pc)^2/Myr
Omega = 40*MYR*KM/(1000*PC)        # angular velocity, converted from km/s/kpc to 1/Myr
q = 0.98                            # shear parameter
t_max = 10                          # total simulation time
z_min = -1.0                        # minimum thickness of the disc
z_max = 1.0                         # thickness of the disc
dt = t_max/200                      # time step
dz = 0.01                           # spatial step in z direction
alpha_0 = 2/h#10*1e3*MYR/(100*PC)      # alpha effect, converted from km/s to 100pc/Myr

print('eta', eta_T, '\t',  'alpha',alpha_0, '\t', 'omega', Omega)
print('Rw', -q*Omega*h**2/eta_T)
print('Ra', alpha_0*h/eta_T)
print('Dynamo number ', -alpha_0*q*Omega*h**3/eta_T**2)

# Spatial grid
z = np.linspace(z_min, z_max, int((z_max - z_min) / dz) + 1)
t = np.linspace(0, t_max, int(t_max / dt) + 1)

# Coefficients for the matrix A and B
nu = eta_T*dt/(2*dz**2)
mu = dt/(2*dz)
alpha = alpha_0*np.sin(np.pi*z/2)
dalpha_dt = np.gradient(alpha, z)

A = mod_matrix(len(z), 1+2*nu, dz*mu, q*Omega*mu, 1+2*nu, -nu, 0, -q*Omega*mu, -nu, -nu, 0, 0, -nu, alpha, dalpha_dt, dz, dt)
B = mod_matrix(len(z), 1-2*nu, -dz*mu, -q*Omega*mu, 1-2*nu, nu, 0, q*Omega*mu, nu, nu, 0, 0, nu, alpha, dalpha_dt, dz, dt)

# Solve the diffusion equation in radial direction
solution = crank_nicolson_mod(len(z), len(t), init_cond_S(z), init_cond_T(z), A, B)

S = solution[:len(z), :]
T = solution[len(z):, :]

eta 0.7 	 alpha 0.4 	 omega 0.04087621516526248
Rw -1.430667530784187
Ra 2.857142857142857
Dynamo number  -4.087621516526249

# Plot the solution in imshow
plot_diff_A(t, z, S, T)
plt.show()

Here is an animation showing the evolution of $S$ and $B_{\phi}$ with time.

One can see that in this example, the fields are decaying with time. We will calculate the decay rate at the midplane, i.e., $z=0$. To calculate the decay rate, we will find how the maxima of the field at the midplane varies with time.

# Constants and parameters
h = 5                               # 100pc
eta_T = 0.7#3.48e-2/h**2                     # magnetic diffusivity in (100pc)^2/Myr
Omega = 40*MYR*KM/(1000*PC)        # angular velocity, converted from km/s/kpc to 1/Myr
q = 0.98                            # shear parameter
t_max = 10                          # total simulation time
z_min = -1.0                        # minimum thickness of the disc
z_max = 1.0                         # thickness of the disc
dt = t_max/200                      # time step
dz = 0.01                           # spatial step in z direction
alpha_0 = 2/h#10*1e3*MYR/(100*PC)      # alpha effect, converted from km/s to 100pc/Myr

print('eta', eta_T, '\t',  'alpha',alpha_0, '\t', 'omega', Omega)
print('Rw', -q*Omega*h**2/eta_T)
print('Ra', alpha_0*h/eta_T)
print('Dynamo number ', -alpha_0*q*Omega*h**3/eta_T**2)

# Spatial grid
z = np.linspace(z_min, z_max, int((z_max - z_min) / dz) + 1)
t = np.linspace(0, t_max, int(t_max / dt) + 1)

# Coefficients for the matrix A and B
nu = eta_T*dt/(2*dz**2)
mu = dt/(2*dz)
alpha = alpha_0*np.sin(np.pi*z/2)
dalpha_dt = np.gradient(alpha, z)

A = mod_matrix(len(z), 1+2*nu, dz*mu, q*Omega*mu, 1+2*nu, -nu, 0, -q*Omega*mu, -nu, -nu, 0, 0, -nu, alpha, dalpha_dt, dz, dt)
B = mod_matrix(len(z), 1-2*nu, -dz*mu, -q*Omega*mu, 1-2*nu, nu, 0, q*Omega*mu, nu, nu, 0, 0, nu, alpha, dalpha_dt, dz, dt)

# Solve the diffusion equation in radial direction
solution = crank_nicolson_mod(len(z), len(t), init_cond_S(z), init_cond_T(z), A, B)

S = solution[:len(z), :]
T = solution[len(z):, :]

eta 0.7 	 alpha 0.4 	 omega 0.04087621516526248
Rw -1.430667530784187
Ra 2.857142857142857
Dynamo number  -4.087621516526249

S_mid = S[int(len(z)*3/4), :]
T_mid = T[int(len(z)*3/4), :]

A_mid = np.sqrt(S_mid**2 + T_mid**2)

plt.figure(figsize=(6, 4))
plt.plot(t, A_mid)
plt.plot(t[int(len(t)*3/4):], A_mid[int(len(t)*3/4):], 'r:', linewidth=4)
plt.xlabel('Time (in Myrs)')
plt.ylabel('Magnetic potential (in $\mu G$/100pc)')
plt.yscale('log')
plt.title('Magnetic potential at the maximum amplitude point')
plt.show()

decay_rate = get_decay_rate(t, A_mid)
if decay_rate < 0:
    print(r"decay rate at the mid plane:", format(-decay_rate, ".3e"))
else:
    print(r"growth rate at the mid plane:", format(decay_rate, ".3e"))

print('Dynamo number ', -alpha_0*q*Omega*h**3/eta_T**2)

growth rate at the mid plane: 3.601e-01
Dynamo number  -4.087621516526249

Calculation of Critical Dynamo number $D_c$

# Constants and parameters
h = 5                               # 100pc
eta_T = 0.7#3.48e-2/h**2                     # magnetic diffusivity in (100pc)^2/Myr
Omega = 40*MYR*KM/(1000*PC)        # angular velocity, converted from km/s/kpc to 1/Myr
q = 0.98                            # shear parameter
t_max = 10                          # total simulation time
z_min = -1.0                        # minimum thickness of the disc
z_max = 1.0                         # thickness of the disc
dt = t_max/200                      # time step
dz = 0.01                           # spatial step in z direction
# alpha_0 = 7/h#10*1e3*MYR/(100*PC)      # alpha effect, converted from km/s to 100pc/Myr
# alpha = alpha_0*np.sin(np.pi*z/2)
# dalpha_dt = np.gradient(alpha, z)


def f(D):
    # Spatial grid
    z = np.linspace(z_min, z_max, int((z_max - z_min) / dz) + 1)
    t = np.linspace(0, t_max, int(t_max / dt) + 1)

    # Coefficients for the matrix A and B
    nu = eta_T*dt/(2*dz**2)
    mu = dt/(2*dz)

    alpha_0 = -D/(q*Omega*h**3/eta_T**2)
    alpha = alpha_0*np.sin(np.pi*z/2)
    dalpha_dt = np.gradient(alpha, z)

    # Coefficients for the matrix A and B

    A = mod_matrix(len(z), 1+2*nu, dz*mu, q*Omega*mu, 1+2*nu, -nu, 0, -q*Omega*mu, -nu, -nu, 0, 0, -nu, alpha, dalpha_dt, dz, dt)
    B = mod_matrix(len(z), 1-2*nu, -dz*mu, -q*Omega*mu, 1-2*nu, nu, 0, q*Omega*mu, nu, nu, 0, 0, nu, alpha, dalpha_dt, dz, dt)

    # Solve the diffusion equation in radial direction
    solution = crank_nicolson_mod(len(z), len(t), init_cond_S(z), init_cond_T(z), A, B)

    S = solution[:len(z), :]
    T = solution[len(z):, :]
    S_mid = S[int(len(z)/2), :]
    T_mid = T[int(len(z)/2), :]

    A_mid = np.sqrt(S_mid**2 + T_mid**2)

    decay_rate = get_decay_rate(t, A_mid)

    return decay_rate

tol = 1e-4
print('Calculating critical Dynamo number...')
D_c, x_arr, y_arr = bisection(f, -0.1, -30, tol)
print('Critical Dynamo number Dc = ', np.round(D_c, 4))
# print('Value of alpha at D_c =', np.round(D_c/(q*Omega*h**3/eta_T**2)/(1e3*MYR/(100*PC)), 4), 'km/s')

plt.figure(figsize=(6, 4))
plt.plot(x_arr, y_arr)
plt.plot(x_arr, y_arr, 'ro')
plt.xlabel('Number of iterations')
plt.ylabel('Dynamo number')
plt.title('Convergence plot for Critical Dynamo number')
plt.grid()

plt.show()

Calculating critical Dynamo number...
-15.05 1.322168
-7.575 1.012576
-3.8375 0.278833
-1.96875 -0.396996
-2.903125 -0.020143
-3.370312 0.13748
-3.136719 0.060864
-3.019922 0.020932
-2.961523 0.000541
-2.932324 -0.009765
-2.946924 -0.004603
-2.954224 -0.002029
-2.957874 -0.000744
-2.959698 -0.000101
-2.960611 0.00022
-2.960155 5.9e-05
Critical Dynamo number Dc =  -2.9602

---

Results

In this task, we solved the coupled mean field induction equation with the $\alpha - \Omega$ terms involed in terms of the magnetic potential.

We specifically removed the $\alpha^2$ terms and all addtional terms for simplicity. The final solution was plotted for $\psi$ and $T$ for the same seed fields thatw e used in task 2.

We then varied the alpha term and solved to obtain the dynamo number for both the cases again. For the first example we took, we observed that the fields decayed with time.

We then increased the value of alpha in the second example and then observed that the fields grow exponentially now. This was because the $\alpha-\Omega$ effect powers the magnertic field and makes it grow. Animations for the plots have been showed below again.

---

Image showing the decay of magnetic fields $\psi$ and $T$ for $\alpha-\Omega$ dynamo.

Image showing the growth of magnetic fields $\psi$ and $T$ for $\alpha-\Omega$ dynamo.

---

These observations led us that the critical dynamo number, for which the decay rate is zero, shall have a dynamo number value in between these two values.

We iteratively calculated the value of the critical dynamo number and obtained it to be $D_c = -2.96$.

Comparison

We observe that upon changing the equations from $B_r$ and $B_{\phi}$ to the magnetic potential form $\psi$ and $T$, there is a change in the critical dynamo number. We kept everything else same as previous while varying just the alpha term for both the cases to obtain the critical dynamo number comparable.

---

---