Radiation-Reaction: Self-consistent Light-Matter interaction with real-time TDDFT

In this tutorial, we calculate the superradiant absorption-spectrum of a sodium-dimer chain.

The equations underlying the implementation are described in Ref. 1. The here used parameters are coarser than the published work. This implementation is so far only available for the LCAO mode.

We use real-time TDDFT LCAO mode and add the radiation-reaction potential for emission into a waveguide. The available implementation adds a local potential that accounts for the recoil that emerges as a consequence of the emission of light. This approach effectively embeds the electromagnetic environment into the KS equation and ensures a self-consistent interaction between light and matter. The keyword RRemission demands the cross-sectional area of the waveguide and the polarization of the propagating electric field. With decreasing cross-sectional area, the emission will be stronger. Thus far, 1d-waveguide emission is the only available electromagnetic environment but further additions are scheduled.

We will start by calculating the emission spectrum in z-direction for a single sodium dimer:

from ase import Atoms
from gpaw import GPAW
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.qed import RRemission
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
from gpaw.tddft.spectrum import photoabsorption_spectrum

# Sodium dimer chain
d = 1.104  # N2 bondlength
L = 8.0  # N2-N2 distance
atoms = Atoms('Na2',
              positions=[[L, 0, +d / 2], [L, 0, -d / 2]],
              cell=[L, 8.0, 8.0])
atoms.center()

# Ground-state calculation
calc = GPAW(mode='lcao', h=0.3, basis='dzp', xc='LDA',
            setups={'Na': '1'},
            convergence={'density': 1e-12},
            txt='gs_nad.out')
atoms.calc = calc
energy = atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')

# Time-propagation calculation
# Read converged ground-state file
td_calc = LCAOTDDFT('gs.gpw', rremission=RRemission(35.05, [0, 0, 1]))
# Attach any data recording or analysis tools
DipoleMomentWriter(td_calc, 'dm_nad.dat')
# Kick
td_calc.absorption_kick([0.0, 0.0, 1e-5])
# Propagate
td_calc.propagate(20, 5000)
# Please note that those parameter values are quite course
# and should be properly converged for subsequent applications.

# Calculate spectrum with small artificial broadening
photoabsorption_spectrum('dm_nad.dat', 'spec_nad.dat',
                         folding='Gauss', width=0.005,
                         e_min=0.0, e_max=10.0, delta_e=0.005)

The polarization of the electric field in the waveguide is aligned with the z-axis, exerting a recoil-force on the electronic coordinates that oscillate in this direction. As a result, the electronic motion is damped as energy is emitted into the waveguide, in line with Newtons third law.

Let us next add a second dimer parallel to the first one and orthorgonal to the chain-axis (H-aggregate configuration). For reasons that will become clear in the following, we can get away with half the propagation time.

from ase import Atoms
from gpaw import GPAW
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.qed import RRemission
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
from gpaw.tddft.spectrum import photoabsorption_spectrum

# Sodium dimer chain
d = 1.104  # N2 bondlength
L = 8.0  # N2-N2 distance
atoms = Atoms('Na4',
              positions=([[L, 0, +d / 2], [L, 0, -d / 2], [2 * L, 0, +d / 2],
                         [2 * L, 0, -d / 2]]), cell=[L + 2 * L, 8.0, 8.0])
atoms.center()

# Ground-state calculation
calc = GPAW(mode='lcao', h=0.3, basis='dzp', xc='LDA',
            setups={'Na': '1'},
            convergence={'density': 1e-12},
            txt='gs_nad2.out')
atoms.calc = calc
energy = atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')

# Time-propagation calculation
# Read converged ground-state file
td_calc = LCAOTDDFT('gs.gpw', rremission=RRemission(35.05, [0, 0, 1]))
# Attach any data recording or analysis tools
DipoleMomentWriter(td_calc, 'dm_nad2.dat')
# Kick
td_calc.absorption_kick([0.0, 0.0, 1e-5])
# Propagate
td_calc.propagate(20, 2500)
# Please note that those parameter values are quite course
# and should be properly converged for subsequent applications.

# Calculate spectrum with small artificial broadening
photoabsorption_spectrum('dm_nad2.dat', 'spec_nad2.dat',
                         folding='Gauss', width=0.005,
                         e_min=0.0, e_max=10.0, delta_e=0.005)

We smoothed the spectra in both cases with a sharp Lorentzian for visual appearance but the observed width corresponds to the correct lifetimes. In the perturbative limit, this corresponds to Wigner-Weisskop theory. Two features can be observed.

../../../_images/spectra_nad.png

1) The width of two sodium dimers is double the width of a single one. For a set of non-interacting matter-subsystems, their emission probability of a single photon scales linearly with the number of subsystem, i.e., the rate increases here linearly with the number of dimers. Intuitively, more charge oscillating provides a larger dipole which leads to a stronger emission. Clearly, the dimers are not entirely non-interacting for distance of 0.8 nm, which leads to the following observation.

2) The peak is shifted slightly to higher frequencies. The parallel configuration between the dimers results in a coupling between the dipolar excitations that shift the excitation energy to slightly higher values. This configuration is also known as H-aggregate. A head to end configuration on the other hand is known as J-aggregate, providing a red-shift.

It is important to realize here that the propagation time has to be long enough to allow the recoil-forces to substantially damp the oscillations, i.e., enough time is needed to capture the decay.

../../../_images/dipoles.png

Often this might be unfeasible or unsuitable such that an extrapolation from small cross-sectional areas and quick decay to larger areas and slower decay might be beneficial. A generalized implementation, including strong coupling and complex electromagnetic environments is currently in develop.

References

1