Tutorial__GS__T03_GS_TDSE - fmauger1/QMol-grid GitHub Wiki
Tutorial 3: Time-dependent Schrödinger equation
This tutorial walks through the process of setting up and calculating the ionization yield of a one-dimensional hydrogen molecular ion model subjected to a few-cycle laser pulse using time-dependent Schrödinger equation.
Table of Contents
Introduction
In atomic units ($\hbar =1$), the time-dependent Schrödinger equation (TDSE) is given by
$$ i\partial_t \psi (x;t)=\hat{\mathcal{H}} \psi (x;t),~~~~~~(1) $$
where $\hat{\mathcal{H}}$ is the TDSE Hamiltonian operator
$$ \hat{\mathcal{H}} ={\hat{\mathcal{H}} }_0 +\hat{x} E(t)~~~~~~\textrm{with}~~~~~~{\hat{\mathcal{H}} }_0 =-\frac{\Delta }{2}+\hat{\mathcal{V}} ,~~~~~~(2) $$
decomposed between a time-independent atomic or molecular part ${\hat{\mathcal{H}} }_0$ and, for field-driven systems, a time-dependent external driving electric field $E(t)$ taken in the dipole approximation.
In this tutorial, we illustrate how to use the QMol-grid package to calculate the ionization yield of a one-dimensional hydrogen molecular ion model subjected to a 800-nm, 10-laser-cycle (total duration) laser pulse with a peak-field intensity of $5\times 10^{14}$ W/cm $^2$ . Specifically, it walks through (i) setting up all the components required for the TDSE and (ii) propagating and calculating the resulting ionization yield.
Setting up the TDSE model
Here we present a minimal example for a TDSE calculation. We discuss additional options and more advanced features in the next tutorial.
System model
We define a Schrödinger equation object corresponding to the hydrogen molecular ion we aim to model (see tutorial 1 and tutorial 2 for details regarding how to set up such systems).
% Molecular parameters
R = 2;
H1 = QMol_Va_softCoulomb('softeningParameter',sqrt(2),'name','H1','position',-R/2);
H2 = QMol_Va_softCoulomb('softeningParameter',sqrt(2),'name','H2','position', R/2);
% Schrödinger equation model
SE = QMol_SE( ...
'xspan', -35:.1:35, ...
'potential', QMol_SE_V('atom',{H1,H2}), ...
'numberWaveFunction', 1);
% Check model
SE.initialize;
SE.showDocumentation;
yielding (only showing the relevant portion of the run-time documentation)
(...)
=== Discretization =======================================================
* Domain discretization Cartesian grid
Grid = -35:0.1:35
Size = 701 (prime) points
V-01.21.000 (06/17/2024) F. Mauger
=== Potential ============================================================
* Potential
V-01.21.000 (06/17/2024) F. Mauger
* Atom-like center(s)
> H1, parameterized as (soft Coulomb)
Z = 1.00 | a = 1.41 | X0 = -1.00
> H2, parameterized as (soft Coulomb)
Z = 1.00 | a = 1.41 | X0 = 1.00
* Soft-Coulomb potential [Javanainen 1988] (soft Coulomb)
Parameterized as V(x) = -Z ./ sqrt( (x-X0).^2 + a^2 ).
V-01.21.000 (06/17/2024) F. Mauger
=== System model =========================================================
* Electronic structure wave functions
1 wave function(s)
(...)
Here we select a larger domain than in the previous tutorials to make room for the absorbing boundary conditions responsible for removing outgoing portions of the wave packet on the edges of the simulation box.
Initial condition
We want to start the TDSE simulation from the ground-state of the ${\textrm{H}}_2^+$ model, and then calculate the ground state associated with SE.
GSS = QMol_SE_eigs;
GSS.computeGroundState(SE);
yielding (again, only showing the relevant part of the printed output)
(...)
=== Wave-function energies ===============================================
Wave fcn Energy (-eV) Error(a.u.)
-------- ------------ -----------
1 26.866 3.898e-13
----------------------------------------------
=== Schrodinger-equation-component energies ==============================
Component Energy (a.u.) Energy (eV)
----------- ------------- -------------
Kinetic 0.086 2.332
Potential -1.073 -29.198
----------- ------------- -------------
Total -0.987 -26.866
----------------------------------------------
giving a ground-state energy reasonably close to the actual one of 30.005 [NIST] (one could match the ionization-potential energy by adjusting the softening parameter in the two H-atom pseudopotentials).
Note that the QMol-grid package decouples the initial condition definition from the TDSE propagation. Here we have elected to start from the molecular ground state but an arbitrary wave function can be defined in SE.waveFunction.waveFunction, provided it has a compatible shape.
Driving field
We now turn to the definition of the external driving field to be used in the TDSE simulations. Like all other input parameters in the QMol-grid package, the driving field must be specified in atomic units. The external components convertUnit and convertLaser can be used to facilitate this definition. For our case, we have
% Laser parameters
w = convertLaser.wavelength2frequency(800e-9);
E0 = convertLaser.intensity2amplitude(5e14);
T = 10*2*pi/w;
% Electric field
phi_CEP = w*T/2;
Env = @(t) sin(pi*t/T).^2.*(t>0).*(t<T);
EField = @(t) E0*cos(w*t-phi_CEP).*Env(t);
The Env function handle defines the overall envelope of the field, here a pulse with a total duration of 10 optical cycles shaped as a sine squared function. The .*(t>0).*(t<T) ensures that the field vanishes before $t=0$ and after $t=10$ o.c. From there, the EField function handle defines the pulse, as the product of a CW wave with frequency w corresponding to 800-nm wavelength, the peak amplitude corresponding to $5\times 10^{14}$ W/cm $^2$ , and the envelope. phi_CEP defines the carrier-envelope phase, here set to have the maximum of the envelope matching with a local maximum of the driving field. As a sanity check, we plot the electric field with
t = -20:.01:(T+20);
figure
plot(t,EField(t),'-','LineWidth',2)
hold off
set(gca,'box','on','FontSize',12,'LineWidth',2)
xlabel('time (a.u.)')
ylabel('E field (a.u.)')
xlim(t([1 end]))
producing
TDSE propagation
We now have all the elements to define the TDSE propagator. Here we choose a fourth-order Forest Ruth [Forest 1990] symplectic split-operator scheme [Mauger 2024] implemented in the QMol_TDSE_SSO_4FR class. From experience, this propagator provides a balanced computation efficacy (see TDSE documentation page for a list of available propagators).
ABC = QMol_TDSE_abs_CAP('length',20);
TDSE = QMol_TDSE_SSO_4FR( ...
'time', linspace(0,T+200,11), ...
'timeStep', 2e-2, ...
'electricField', EField, ...
'absorbingBoundary', ABC);
ABC defines the the complex-absorbing-potential (CAP) object to be used by the propagator TDSE, to remove outgoing parts of the wave packets from the simulations domain. The 'time' parameter defines the starting and final simulation times, as the first and last elements in the corresponding input vector, respectively. Intermediate times are used for the simulation-completion display but have no bearing on the propagation time step used in the simulation; this is independently set by the 'timeStep' input parameter (here 0.02 a.u.). For our ionization simulation, we let the propagation go on for 200 a.u. after the end of the driving pulse. This is to give enough time for all the ionized portions of the wave packet to reach the absorber region and be removed from the box. The driving external field, defined via its electric field, is registered with the 'electricField' parameter.
Like for Schrödinger equation objects, we can check the TDSE object configuration using the run-time documentation with
TDSE.initialize;
TDSE.showDocumentation;
yielding (as usual, only showing the relevant part of the printed output)
(...)
=== Time propagation =====================================================
* The TDSE model dynamics matches the canonical Hamiltonian formalism
discussed in [Mauger 2024].
* Symplectic split-operator propagation scheme [Mauger 2024],
with a VTV split motif.
V-01.21.000 (06/17/2024) F. Mauger
* Field driven propagation
dipole approximation, length gauge
* Forest-Ruth propagation scheme [Forest 1990] 4th order
V-01.21.000 (06/17/2024) F. Mauger
* Parameters
Time interval = 0.00 to 1303.19 a.u. (0.00 to 31.52 fs)
Time step = 0.020 a.u. (0.484 as)
* Absorbing boundaries complex absorbing potential
Length = 20 a.u. on both ends
Amplitude = 0.5 a.u. Shape = sin^1/8
V-01.21.000 (06/17/2024) F. Mauger
* Restart option is disabled
=== Output results =======================================================
=== References ===========================================================
[Forest 1990] E. Forest and R.D. Ruth, "Fourth-order symplectic
integration," Physica D: Nonlinear Phenomena 43, 105 (1990).
[Mauger 2024] F. Mauger, C. Chandre, M.B. Gaarde, K. Lopata, and K.J.
Schafer, "Hamiltonian formulation and symplectic split-operator
schemes for time-dependent density-functional-theory equations of
electron dynamics in molecules," Communications in Nonlinear Science
and Numerical Simulation 129, 107685 (2024).
[Mauger 2024b] F. Mauger and C. Chandre, "QMol-grid: A MATLAB package
for quantum-mechanical simulations in atomic and molecular systems,"
arXiv:2406.17938 (2024).
(...)
The selection of the length gauge for the field-driven propagation stems from our choice to input the electric field (this can be overwritten). The empty "Output results" section indicates that only the propagation of the wave function is performed and no on-the-fly observable is calculated and saved along the way.
TDSE propagation
To run the TDSE propagation simulation, we simply pass the Schrödinger equation object SE to the TDSE propagator.
TDSE.propagate(SE);
yielding (as usual, only showing the relevant part of the printed output)
(...)
=== Time propagation =====================================================
Iteration time (a.u.) Iteration Iteration progress (%)
start finish number 0 20 40 60 80 100
--------- --------- --------- -------------------------------
0.000 130.320 1 |||||||||||||||||||||||||||||||
130.320 260.640 6517 |||||||||||||||||||||||||||||||
260.640 390.960 13033 |||||||||||||||||||||||||||||||
390.960 521.280 19549 |||||||||||||||||||||||||||||||
521.280 651.600 26065 |||||||||||||||||||||||||||||||
651.600 781.920 32581 |||||||||||||||||||||||||||||||
781.920 912.240 39097 |||||||||||||||||||||||||||||||
912.240 1042.560 45613 |||||||||||||||||||||||||||||||
1042.560 1172.880 52129 |||||||||||||||||||||||||||||||
1172.880 1303.190 58645 |||||||||||||||||||||||||||||||
--------- --------- --------- -------------------------------
Time propagation finished without error. All results are saved in the
TDSE object and specified output MATLAB file(s).
##########################################################################
Like for the ground-state calculation of tutorial 1, the final wave function (at the end of the propagation) is stored in the SE object. To calculate the ionization yield, we look at how much density has left the final wave function.
dx = SE.xspan(2)-SE.xspan(1);
p = sum(abs(SE.waveFunction.waveFunction).^2)*dx;
fprintf('Ion yield = %#10.3e elec.\n',1-p);
yielding
Ion yield = 3.500e-01 elec.
i.e., a 35% ionization yield.
References
[Forest 1990] E. Forest and R.D. Ruth, "Fourth-order symplectic integration," Physica D: Nonlinear Phenomena 43, 105 (1990).
[Mauger 2024] F. Mauger, C. Chandre, M.B. Gaarde, K. Lopata, and K.J. Schafer, "Hamiltonian formulation and symplectic split-operator schemes for time-dependent density-functional-theory equations of electron dynamics in molecules," Communications in Nonlinear Science and Numerical Simulation 129, 107685 (2024).
[NIST] Electron-Impact Cross Sections for Ionization and Excitation Database, Hydrogen molecule ion ( ${\textrm{H}}_2^+$ ), https://physics.nist.gov/cgi-bin/Ionization/table.pl?ionization=H2%2B
Notes
The results displayed in this tutorial were generated using version 01.21 of the QMol-grid package.