solving rate equations for the multi-photon multiple ionization dynamics due to intense x-ray radiation¶
CFEL 2022, Ludger Inhester
In [ ]:
import numpy as np
from scipy.integrate import solve_ivp
from matplotlib import pyplot as plt
AU2SEC = 2.418884326500e-17
FS2AU = 1e-15 / AU2SEC
AU2FS = 1 / FS2AU
AU2M = 0.5291772083e-10
M2AU = 1. / AU2M
In [ ]:
def solveREQ(pIni,
processData,
configurations,
fluence,
pulseEnv,
tStart=-150 * FS2AU,
tEnd=150 * FS2AU,
dt=0.1 * FS2AU,
method='LSODA'):
"""
solves the rate equation and returns results for the population
Mandatory arguments:
pIni: initial population, array of shape(nconfigurations)
processData: process data, dictionary obtained from readProcessData
configurations: array of configurations, array of shape (nconfigurations, number of occupations)
fluence: fluence value, scalar (number of photons per area in atomic units)
pulseEnv: pulseEnvelop, function f(t), t in atomic units
optional arguments:
tStart: time to start integration (i.e., where p=pIni) in atomic units, default -100*FSAU
tEnd: time to end integration in atomic units, default 100*FSAU
dt: timestep in atomic units for which to provide results, default 0.1*FSAU.
The timestep is also used as maximal timestep for the integrator.
It must be equal or smaller than the pulse duration.
method: one of 'RK23', 'RK45', 'BDF', 'LSODA', 'Radau', 'DOP853'
returns:
res: return results object of scipy.integrate.solve_ivp
in particular
res['t'] time values
res['y'] population values at different times
"""
confDict = {str(c): i for i, c in enumerate(configurations)}
gammas = [{
'i': confDict[str(proc['iconfig'])],
'f': confDict[str(proc['fconfig'])],
'rate': proc['rate']
} for proc in processData if proc['type'] == 'A' or proc['type'] == 'F']
sigmas = [{
'i': confDict[str(proc['iconfig'])],
'f': confDict[str(proc['fconfig'])],
'pcs': proc['pcs']
} for proc in processData if proc['type'] == 'P' or proc['type'] == 'V']
def REQ(t, p):
"""
returns dp/dt (the right-hand-side of the diff. eq.)
"""
pDot = np.zeros(p.shape)
Ipulse = pulseEnv(t) * fluence
for s in sigmas:
pDot[s['i']] += -s['pcs'] * Ipulse * p[s['i']]
pDot[s['f']] += +s['pcs'] * Ipulse * p[s['i']]
for g in gammas:
pDot[g['i']] += -g['rate'] * p[g['i']]
pDot[g['f']] += +g['rate'] * p[g['i']]
return (pDot)
def REQjac(t, p):
"""
returns d/dpj dp_i/dt (jacobian of the right-hand-side of the diff. eq.)
"""
pDot = REQ(t, p)
pDotjac = np.zeros((p.shape[0], p.shape[0]))
Ipulse = pulseEnv(t) * fluence
for s in sigmas:
pDotjac[s['i'], s['i']] += -s['pcs'] * Ipulse
pDotjac[s['f'], s['i']] += +s['pcs'] * Ipulse
for g in gammas:
pDotjac[g['i'], g['i']] += -g['rate']
pDotjac[g['f'], g['i']] += +g['rate']
return (pDotjac)
times = np.arange(tStart, tEnd, dt)
if method == 'RK45' or method == 'RK23' or method == 'Radau' or method == 'DOP853':
res = solve_ivp(REQ, (tStart - dt, tEnd + dt),
pIni,
t_eval=times,
method=method,
max_step=dt)
elif method == 'LSODA' or method == 'BDF':
res = solve_ivp(REQ, (tStart - dt, tEnd + dt),
pIni,
t_eval=times,
method='LSODA',
jac=REQjac,
max_step=dt)
else:
raise ValueError(f"method={method} unknown")
if not res.success:
print(res.message)
return res
In [ ]:
def readProcessData(filename):
"""
reads in and returns process data generated by xatom from file filename
"""
processData = []
with open(filename, 'r') as f:
lines = f.readlines()
for line in lines:
line = line.rstrip()
els = line.split()
if len(els) == 0:
continue
configs = line.split('#')[1].split('->')
if (els[0] == 'P' or els[0] == 'V'):
processData.append({
"type":
els[0],
"i":
els[1],
"f":
els[2],
"pcs":
float(els[3]),
"energy":
els[4],
"iconfig":
np.array([int(n) for n in configs[0].split()]),
"fconfig":
np.array([int(n) for n in configs[1].split()]),
})
else:
processData.append({
"type":
els[0],
"i":
els[1],
"f":
els[2],
"rate":
float(els[3]),
"energy":
els[4],
"iconfig":
np.array([int(n) for n in configs[0].split()]),
"fconfig":
np.array([int(n) for n in configs[1].split()]),
})
return processData
In [ ]:
def gaussian(x, sigma=10 * FS2AU, x0=0.):
"""
Gaussian function:
x0 is the center
sigma is the stdev
"""
return 1. / np.sqrt(2 * np.pi) * 1 / sigma * np.exp(-0.5 *
((x - x0) / sigma)**2)
Initial 3-configuration example¶
In [ ]:
pIni = np.array([1., 0., 0., 0.])
sigma = 1.
Gamma = 3 / 2
processData = [
{
"iconfig": np.array([1, 1]),
"fconfig": np.array([0, 1]),
"pcs": sigma,
"type": "P"
},
{
"iconfig": np.array([0, 1]),
"fconfig": np.array([0, 0]),
"pcs": sigma,
"type": "P"
},
{
"iconfig": np.array([0, 1]),
"fconfig": np.array([1, 0]),
"rate": Gamma,
"type": "A"
},
]
configurations = [
np.array([1, 1]),
np.array([0, 1]),
np.array([0, 0]),
np.array([1, 0])
]
fluence = 1.
pulseEnv = lambda t: 1. if t > 0 else 0.
res = solveREQ(pIni,
processData,
configurations,
fluence,
pulseEnv,
tStart=0,
tEnd=10,
dt=0.1,
method='LSODA')
In [ ]:
times = res['t']
p = res['y']
# plot populations and compare to analytical solution
fig, ax = plt.subplots()
ax.plot(times, p[0, :], label=r'$p_1$')
ax.plot(times, p[1, :], label=r'$p_2$')
ax.plot(times, p[2, :], label=r'$p_3$')
ax.plot(times, p[3, :], label=r'$p_4$')
ax.set_xlabel('time')
ax.set_ylabel('population')
ax.legend()
plt.show()
Exercise 3a put your code in the cells below¶
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
Rate Equation Solving with Calculated Atomic Data¶
read an preprocess the process data (rates and cross sections calculated by xatom)¶
In [ ]:
# read in process data: neon, 2000eV
processData = readProcessData('results_Ne_2000.dat')
# get a list of involved electronic configurations
configurations = np.unique(np.array([p['iconfig'] for p in processData] +
[p['fconfig'] for p in processData]),
axis=0)
Nc = configurations.shape[0]
# compute number of electrons for each configuration
Nelectron = np.array([np.sum(c) for c in configurations])
# initial population we start with the configuration that has the maximum number of electrons
startCfgIdx = np.argmax(Nelectron)
pIni = np.zeros(Nc)
pIni[startCfgIdx] = 1.
# Assume that nuclear charge is identical to maximum number of electron,
# i.e., the configuration with the most electrons is neutral
nuclearCharge = np.max(Nelectron)
solve rate equations¶
In [ ]:
%%time
# set the pulse parameters:
# the fluence number of photons per micrometer square
nPhotonsPerSqmuM = 1e12
# the pulse width (FWHM)
fwhmPulse = 25 * FS2AU
# convert fluence into atomic units and pulse width FWHM into stdev (atomic units)
fluence = nPhotonsPerSqmuM / (1e-6 * M2AU)**2
sigmaPulse = fwhmPulse / (2. * np.sqrt(2. * np.log(2.)))
# pulse envelope: simple Gaussian with stddev 25fs and t0 = 0fs
pulseEnvelope = lambda t: gaussian(t, sigma=sigmaPulse, x0=0.0 * FS2AU)
res = solveREQ(pIni,
processData,
configurations,
fluence,
pulseEnvelope,
method='LSODA')
# p is population, times is timepoints
p = res['y']
times = res['t']
# compute the average charge per time
avgCharge = nuclearCharge - np.dot(Nelectron, p)
# compute the population in a given charge state per time
chargePops = np.zeros((nuclearCharge + 1, p.shape[1]))
for q in range(chargePops.shape[0]):
idx = np.where(nuclearCharge - Nelectron == q)[0]
chargePops[q, :] = np.sum(p[idx, :], axis=0)
plot/analyze the results¶
In [ ]:
# plot average charge
fig, ax = plt.subplots()
ax2 = ax.twinx()
ax.plot(times / FS2AU, avgCharge, label='avg. charge')
ax2.plot(times / FS2AU,
fluence * pulseEnvelope(times) * 1 / AU2FS * 1 / ((1e6 * AU2M)**2),
label='pulse',
c='black',
alpha=0.3)
ax.set_xlabel('time (fs)')
ax.set_ylabel('avg. charge')
ax2.set_ylabel(r'intensity (n. photons / $\mathrm{\mu m}^2$ fs)')
ax.legend()
ax2.legend()
plt.show()
In [ ]:
# plot charge population
fig, ax = plt.subplots()
ax2 = ax.twinx()
for q in range(chargePops.shape[0]):
ax.plot(times / FS2AU, chargePops[q, :], label=f'q={q}')
ax2.plot(times / FS2AU,
fluence * pulseEnvelope(times) * 1 / AU2FS * 1 / ((1e6 * AU2M)**2),
label='pulse',
c='black',
alpha=0.3)
ax.legend()
ax.set_xlabel('time (fs)')
ax.set_ylabel('population')
ax2.set_ylabel(r'intensity (n. photons / $\mathrm{\mu m}^2$ fs)')
plt.show()
In [ ]:
# ratio of PP vs PA processes:
# i.e. either A or P starting from configuration 1,2,6
# determine the singly core ionized configuration
cfgP = [
p['fconfig'] for p in processData
if np.all(p['iconfig'] == np.array(configurations[startCfgIdx]))
and p['type'] == 'P'
][0]
print(
f"The singly core ionized configuration has electronic occupation numbers {cfgP}"
)
# index of singly core ionized configuration
idxP = np.argmax(np.all(configurations == cfgP, axis=1))
# find the total photoionization c.s. and the Auger-Meitner decay rate for configuration idxP
sigmaPP = np.sum([
p['pcs'] for p in processData
if np.all(p['iconfig'] == configurations[idxP]) and p['type'] == 'P'
])
GammaPA = np.sum([
p['rate'] for p in processData
if np.all(p['iconfig'] == configurations[idxP]) and p['type'] == 'A'
])
# probability for PP is
# int dt sigma * F(t) p_P(t)
# probability for PA is
# int dt Gamma * p_P(t)
dt = times[1] - times[0]
probPP = np.sum(sigmaPP * fluence * pulseEnvelope(times) * p[idxP, :]) * dt
probPA = np.sum(GammaPA * p[idxP, :]) * dt
print(
f"The probability for PP is {probPP}.\nThe probability for PA is {probPA}.\nRatio PP/PA is {probPP/probPA}"
)
Exercise 3b, 3c, 3d in the cells below¶
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: