Solve the mean-field $\alpha-\Omega$ dynamo equations in the kinematic regime. That is, include the $\Omega$ effect term in the equation for $\dfrac{\partial \bar{B}_\phi}{\partial t}$ and the $\alpha$ effect term in the equation for $\dfrac{\partial \bar{B}_r}{\partial t}$. This requires specifying the overall magnitude and spatial dependence of $\Omega$ and $\alpha$.
where $q = − \dfrac{d \ln \Omega}{d \ln r}$ and $\alpha_0 > 0$ is the amplitude of the $\alpha$ effect. Note that $q > 0$ if $\Omega$ decreases with $r$, which is generally the case in galaxies, so $D < 0$.
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, its numerical solution using Crank-Nicholson method and some example seed fields, and finally calculated the decay rate.
In this chapter, we will begin with the kinematic regime of the solutions, beginning with the addition of the $\alpha-\Omega$ term in uor equations. Lets revisit the mean-field induction equation.
where $\mathcal{E} = \left( \alpha \bar{\mathbf{B}} \right) - \eta_t \left( \nabla \times \bar{\mathbf{B}} \right)$
We will solve the equations in the cylindrical coordinates (r, $\phi$, z) with the origin at the galactic centre and the z-axis parallel to the galactic angular velocity. However, to simplify things, lets make some approximations again.
But $\nabla \times \left( \nabla \times \bar{\mathbf{B}} \right) = \nabla \left( \nabla \cdot \bar{\mathbf{B}} \right) - \nabla^2 \bar{\mathbf{B}} $ and $\nabla \cdot \bar{\mathbf{B}} = 0$ (Gauss's Law). So we get,
$$ \boxed{ \dfrac{\partial \bar{\mathbf{B}}}{\partial t} = \nabla \times \left( \bar{\mathbf{V}} \times \bar{\mathbf{B}} \right) + \nabla \times \left(\alpha \bar{\mathbf{B}} \right) - \eta_T \left( \nabla \times \nabla \times \bar{\mathbf{B}} \right) } $$Here we have neglected all terms which involved $\dfrac{\partial}{\partial r}$ in $B$. We have also omitted the $\alpha^2$ term for now.
The dynamo number is defined as
$$ D = − \dfrac{\alpha_0 q \Omega_0 h^3}{\eta_t^2} $$Here $h$ is the thickness of the thin disk of the galaxy, $\alpha$ is the term responsible for the twisting of the toroidal fields into poloidal fields, $\Omega$ is the rotation rate fo the galaxy, causing the twisting of poloidal fields into toroidal fields, and $q = -\dfrac{\partial \ln \Omega}{\partial \ln r} = -\dfrac{r}{\Omega} \dfrac{\partial \Omega}{\partial r}$.
One can parameterize $\Omega$ as $$ \Omega = \dfrac{\Omega_0}{\sqrt{1 + \left( \dfrac{r}{r_{\omega}} \right)^2 }} $$
Thus $q$ solves as $$ q = \left( \dfrac{r^2}{r_{\omega}^2} \right) \left( 1 + \dfrac{r^2}{r_{\omega}^2} \right)^{-1} $$
For Milky Way galaxy, we have the $r_{\omega} = 2$ kpc, AND the peak rotational velocity $v$ AT $r_{\omega}$ 220 km/s. This gives us $\Omega_0 = \dfrac{220 \text{ km/s}}{2 \text{ kpc}} = 110$ km/s/kpc. Putting these, we get $q = 0.98$.
We will also get an estimate of the value of $\alpha$ $$ \begin{aligned} \alpha &= \dfrac{\tau^2 v^2 \Omega}{h} \\ &= \dfrac{(10 \text{ Myrs})^2 \times (10 \text{ km/s})^2 \times 110 \text{ km/s/kpc} }{ 100 \text{ pc}} \\ &= 11 \text{ km/s} \quad = \quad 0.1124 \text{ (100 pc/Myrs)} \end{aligned} $$
However, in this chapter, we will vary $\alpha_0$ and see how the magnetic fields behave.
With the values $\eta_T = 3.48 \times 10^{-2}\text{ (100 pc)}^2/\text{ Myr}$ (as chosen previously), $h = 100$ pc (typical thickness of thin disks), the critical dynamo number is found to be, $$ \begin{aligned} D &= − \dfrac{\alpha \times 0.98 \times 110 \text{ km/s/kpc} \times (100 \text{ pc})^3}{(3.48 \times 10^{-2}\text{ (100 pc)}^2/\text{ Myr})^2} \\ &= -0.948 \: \alpha \end{aligned} $$
where $\alpha$ is in km/s.
Now we begin to solve the equations numerically.
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.
Let $\bar{B}_r = P$ and $\bar{B}_\phi = Q$ for understanding purposes. Our equations modify to $$ \frac{\partial P}{\partial t} = - \frac{\partial (\alpha Q)}{\partial z} + \eta_T \frac{\partial^2 P}{\partial z^2} $$ $$ \frac{\partial Q}{\partial t} = -q \Omega P + \eta_T \frac{\partial^2 Q}{\partial z^2} $$
Lets first discretize equations, $$ \frac{P^{j+1}_{i} - P^{j}_{i}}{dt} = - \left. \dfrac{d \alpha}{dz}\right \vert_i \left( \frac{Q^{j+1}_{i} + Q^{j}_{i}}{2} \right) - \dfrac{\alpha_{i}}{2} \: \left( \frac{Q^{j+1}_{i+1} - Q^{j+1}_{i}}{dz} + \frac{Q^{j}_{i+1} - Q^{j}_{i}}{dz} \right) + \dfrac{\eta_T}{2} \: \left( \dfrac{P^{j+1}_{i+1} - 2P^{j+1}_{i} + P^{j+1}_{i-1}}{dz^2} + \dfrac{P^{j}_{i+1} - 2P^{j}_{i} + P^{j}_{i-1}}{dz^2} \right) $$ and $$ \frac{Q^{j+1}_{i} - Q^{j}_{i}}{dt} = -q\Omega \left( \frac{P^{j+1}_{i} - P^{j}_{i}}{2} \right) + \dfrac{\eta_T}{2} \: \left( \dfrac{Q^{j+1}_{i+1} - 2Q^{j+1}_{i} + Q^{j+1}_{i-1}}{dz^2} + \dfrac{Q^{j}_{i+1} - 2Q^{j}_{i} + Q^{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) P^{j+1}_{i} - \nu P^{j+1}_{i+1} - \nu P^{j+1}_{i-1} + \mu \alpha_i Q^{j+1}_{i+1} + \left( \dfrac{dt}{2} \left. \dfrac{d \alpha}{dz}\right \vert_i - \mu \alpha_i \right) Q^{j+1}_{i} = \left(1-2\nu \right) P^{j}_{i} + \nu P^{j}_{i+1} + \nu P^{j}_{i-1} - \mu \alpha_i Q^{j}_{i+1} + \left( - \dfrac{dt}{2} \left. \dfrac{d \alpha}{dz}\right \vert_i + \mu \alpha_i \right) Q^{j}_{i} $$ and $$ \left(1+2\nu \right) Q^{j+1}_{i} - \nu Q^{j+1}_{i+1} - \nu Q^{j+1}_{i-1} + \dfrac{dt \: q \Omega}{2} P^{j+1}_{i} = \left(1-2\nu \right) Q^{j}_{i} + \nu P^{j}_{i+1} + \nu P^{j}_{i-1} - \dfrac{dt \: q \Omega}{2} P^{j}_{i} $$
Let $U = \left[ P \:\: Q\right]^T$ This whole coupled equation can be simplified to single variable $U$ as $$ \left[\begin{array}{cc} (1+2 \nu) & \left( \dfrac{dt}{2} \left. \dfrac{d \alpha}{dz}\right \vert_i - \mu \alpha_i \right) \\ \dfrac{dt \: q \Omega}{2} & (1+2 \nu) \end{array}\right] U^{j+1}_{i} + \left[\begin{array}{cc}-\nu & \mu \alpha_i \\ 0 & -\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) & \left( - \dfrac{dt}{2} \left. \dfrac{d \alpha}{dz}\right \vert_i + \mu \alpha_i \right) \\ -\dfrac{dt \: q \Omega}{2} & (1-2 \nu) \end{array}\right] U^{j}_{i} + \left[\begin{array}{cc}\nu & -\mu \alpha_i \\ 0 & \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
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} \left( - \dfrac{dt}{2} \left. \dfrac{d \alpha}{dz}\right \vert_i + \mu \alpha_i \right) & -\mu \alpha_i & 0 & \cdots & 0 \\ 0 & \left( - \dfrac{dt}{2} \left. \dfrac{d \alpha}{dz}\right \vert_i + \mu \alpha_i \right) & -\mu \alpha_i & \ddots & \vdots \\ 0 & 0 & \left( - \dfrac{dt}{2} \left. \dfrac{d \alpha}{dz}\right \vert_i + \mu \alpha_i \right) & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & -\mu \alpha_i \\ 0 & \cdots & 0 & 0 & \left( - \dfrac{dt}{2} \left. \dfrac{d \alpha}{dz}\right \vert_i + \mu \alpha_i \right) \end{array}\right] \\ & \\ \left[\begin{array}{ccccc} -\dfrac{dt\: q \Omega}{2} & 0 & 0 & \cdots & 0 \\ 0 & -\dfrac{dt\: q \Omega}{2} & 0 & \ddots & \vdots \\ 0 & 0 & -\dfrac{dt\: q \Omega}{2} & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & 0 & -\dfrac{dt\: q \Omega}{2} \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} $$
We will use the same grid structure and boumndary conditions that we used in chapter 1.
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 *
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_Br(x):
return -(1-x**2)
def init_cond_Bphi(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 COMPONENTS'
plot_init_cond(z, init_cond_Br, init_cond_Bphi, title_1, title_2, global_title)
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 = matrix(len(z), 1+2*nu, mu, q*Omega*dt/2, 1+2*nu, -nu, mu, 0, -nu, -nu, 0, 0, -nu, alpha, dalpha_dt, dz, dt)
B = matrix(len(z), 1-2*nu, -mu, -q*Omega*dt/2, 1-2*nu, nu, -mu, 0, 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_Br(z), init_cond_Bphi(z), A, B)
B_r = solution[:len(z), :]
B_phi = 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(t, z, B_r, B_phi)
plt.show()
Here is an animation showing the evolution of $B_r$ 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.
Let us first run it again for a longer duration to get the effects clearly.
# 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 = 20 # 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 = matrix(len(z), 1+2*nu, mu, q*Omega*dt/2, 1+2*nu, -nu, mu, 0, -nu, -nu, 0, 0, -nu, alpha, dalpha_dt, dz, dt)
B = matrix(len(z), 1-2*nu, -mu, -q*Omega*dt/2, 1-2*nu, nu, -mu, 0, 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_Br(z), init_cond_Bphi(z), A, B)
B_r = solution[:len(z), :]
B_phi = solution[len(z):, :]
eta 0.7 alpha 0.4 omega 0.04087621516526248 Rw -1.430667530784187 Ra 2.857142857142857 Dynamo number -4.087621516526249
B_mid_r = B_r[int(len(z)/2), :]
B_mid_phi = B_phi[int(len(z)/2), :]
B_mid = np.sqrt(B_mid_r**2 + B_mid_phi**2)
plt.figure(figsize=(6, 4))
plt.plot(t, B_mid)
plt.plot(t[int(len(t)*3/4):], B_mid[int(len(t)*3/4):], 'r:', linewidth=4)
plt.xlabel('Time (in Myrs)')
plt.ylabel('Magnetic field (in $\mu G$)')
plt.yscale('log')
plt.title('Magnetic field at the midplane')
plt.show()
decay_rate = get_decay_rate(t[int(len(t)*3/4):], B_mid[int(len(t)*3/4):])
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: 6.550e-01 Dynamo number -4.087621516526249
For this example, we take an $\alpha$ for which the field is decaying.
def init_cond_Br(x):
return -(1-x**2)
def init_cond_Bphi(x):
return (1-x**2)#-(1-x**2)*np.cos(np.pi*x)
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 COMPONENTS'
plot_init_cond(z, init_cond_Br, init_cond_Bphi, title_1, title_2, global_title)
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 = 7/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 = matrix(len(z), 1+2*nu, mu, q*Omega*dt/2, 1+2*nu, -nu, mu, 0, -nu, -nu, 0, 0, -nu, alpha, dalpha_dt, dz, dt)
B = matrix(len(z), 1-2*nu, -mu, -q*Omega*dt/2, 1-2*nu, nu, -mu, 0, 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_Br(z), init_cond_Bphi(z), A, B)
B_r = solution[:len(z), :]
B_phi = solution[len(z):, :]
eta 0.7 alpha 1.4 omega 0.04087621516526248 Rw -1.430667530784187 Ra 10.0 Dynamo number -14.306675307841868
# Plot the solution in imshow
plot_diff(t, z, B_r, B_phi)
plt.show()
Here is an animation showing the evolution of $B_r$ and $B_{\phi}$ with time.
One can see that in this example, the fields are growing with time. We will calculate the growth rate at the midplane, i.e., $z=0$. To calculate the growth rate, we will find how the maxima of the field at the midplane varies with time.
Let us first run it again for a longer duration to get the effects clearly.
# 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 = 20 # 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
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 = matrix(len(z), 1+2*nu, mu, q*Omega*dt/2, 1+2*nu, -nu, mu, 0, -nu, -nu, 0, 0, -nu, alpha, dalpha_dt, dz, dt)
B = matrix(len(z), 1-2*nu, -mu, -q*Omega*dt/2, 1-2*nu, nu, -mu, 0, 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_Br(z), init_cond_Bphi(z), A, B)
B_r = solution[:len(z), :]
B_phi = solution[len(z):, :]
eta 0.7 alpha 1.4 omega 0.04087621516526248 Rw -1.430667530784187 Ra 10.0 Dynamo number -14.306675307841868
B_mid_r = B_r[int(len(z)/2), :]
B_mid_phi = B_phi[int(len(z)/2), :]
B_mid = np.sqrt(B_mid_r**2 + B_mid_phi**2)
plt.figure(figsize=(6, 4))
plt.plot(t, B_mid)
plt.plot(t[150:], B_mid[150:], 'r:', linewidth=4)
plt.xlabel('Time (in Myrs)')
plt.ylabel('Magnetic field (in $\mu G$)')
plt.yscale('log')
plt.title('Magnetic field at the midplane')
plt.show()
decay_rate = get_decay_rate(t[150:], B_mid[150:])
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: 2.664e-01 Dynamo number -14.306675307841868
# 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 = 20 # 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)
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
nu = eta_T*dt/(2*dz**2)
mu = dt/(2*dz)
A = matrix(len(z), 1+2*nu, mu, q*Omega*dt/2, 1+2*nu, -nu, mu, 0, -nu, -nu, 0, 0, -nu, alpha, dalpha_dt, dz, dt)
B = matrix(len(z), 1-2*nu, -mu, -q*Omega*dt/2, 1-2*nu, nu, -mu, 0, 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_Br(z), init_cond_Bphi(z), A, B)
B_r = solution[:len(z), :]
B_phi = solution[len(z):, :]
B_mid_r = B_r[int(len(z)/2), :]
B_mid_phi = B_phi[int(len(z)/2), :]
B_mid = np.sqrt(B_mid_r**2 + B_mid_phi**2)
decay_rate = get_decay_rate(t, B_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 0.325075 -7.575 -0.263914 -11.3125 0.055458 -9.44375 -0.096321 -10.378125 -0.018725 -10.845312 0.018766 -10.611719 0.000124 -10.494922 -0.009274 -10.55332 -0.004569 -10.58252 -0.002221 -10.597119 -0.001048 -10.604419 -0.000462 -10.608069 -0.000169 -10.609894 -2.2e-05 Critical Dynamo number Dc = -10.6099
In this task, we solved the coupled-diffusion equation with the $\alpha - \Omega$ terms involed. We specifically removed the $\alpha^2$ terms and all addtional terms for simplicity. The final solution was plotted for $B_r$ and $B_{\phi}$ for some simple arbitrary example seed fields.
We then varied the alpha term and solved to obtain the dynamo number for each of the cases.
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 $B_r$ and $B_{\phi}$ for $\alpha-\Omega$ dynamo.
Image showing the growth of magnetic fields $B_r$ and $B_{\phi}$ 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 = -10.61$.