# Dielectric response of laser-excited silicon

S.A. Sato  Graduate School of Pure and Applied Sciences, University of
Tsukuba, Tsukuba 305-8571, Japan  K. Yabana  Center for Computational
Sciences, University of Tsukuba, Tsukuba 305-8577, Japan  Graduate School of
Pure and Applied Sciences, University of Tsukuba, Tsukuba 305-8571, Japan  Y.
Shinohara  Graduate School of Pure and Applied Sciences, University of
Tsukuba, Tsukuba 305-8571, Japan  T. Otobe  Advanced Photon Research Center,
JAEA, Kizugawa, Kyoto 619-0615, Japan  G.F. Bertsch  Department of Physics and
Institute for Nuclear Theory, University of Washington, Seattle 98195, U.S.A.

###### Abstract

We calculate the dielectric response of crystalline silicon following
irradiation by a high-intensity laser pulse, modeling the dynamics by time-
dependent density functional theory (TDDFT). The pump-probe measurements are
numerically simulated by solving the time-dependent Kohn-Sham equation with
the pump and probe fields included as external fields. As expected, the
excited silicon shows features of a particle-hole plasma in its response. We
compare the calculated response with a thermal model and with a simple Drude
model. The thermal model requires only a static DFT calculation to prepare
electronically excited matter and agrees rather well with the TDDFT for the
same particle-hole density. The Drude model with two fitted parameters
(electron effective mass and collision time) also shows fair agreement at the
lower excitation energies; the fitted effective masses are consistent with
carrier-band dispersions. The extracted Drude lifetimes range from 6 fs at
weak pumping fields to much lower values at high fields. However, we find that
the Drude model does not give a good fit to the imaginary dielectric function
at the highest fields. Comparing the thermal model with the Drude, we find
that the extracted lifetimes are in the same range, 1-13 fs depending on the
temperature. These short Drude lifetimes show that strong damping is possible
in the TDDFT, despite the absence of electron scattering. One significant
difference between the TDDFT response and the other models is that the
response to the probe pulse depends on the polarization of the pump pulse. We
also find that the imaginary part of the dielectric function can be negative,
particularly for the parallel polarization of pump and probe fields.

##  I introduction

The interaction of high-intensity and ultra-short electromagnetic fields with
condensed matter is an important subject from both fundamental and
technological points of view paris ; ps99 ; bk00 ; cm07 . To investigate
dynamics of electrons and phonons in real time, the pump-probe experimental
technique has been extensively employed. One example of its use is creating
coherent phonons and measuring their properties me97 . The vibration is
detected by measuring the change of reflectance of the probe pulse. However,
this requires a good understanding of the dielectric properties of the surface
excited by the pump pulse. Another example is the energy deposited by strong
laser pulses close to the damage threshold. They produce high-density
electron-hole pairs at the surface of dielectrics, causing strong reflection
for the probe pulse so00 . It is even now possible to measure the population
of high-density electron-hole pairs in the time resolution less than a
femtosecond sb13 ; sp13 ; mi11 . However, the existing theory describing these
effects is largely phenomenological. The dielectric properties of laser-
excited material are often modeled with the Drude model re04 ; me10 ; re10 ;
ka00 , assuming that excited electrons behave like free carriers.

We have been developing a first-principles theoretical approach to describe
electron dynamics in crystalline solids induced by the intense and ultrashort
laser pulses. We rely upon the time-dependent density functional theory
(TDDFT) rg84 , solving the time-dependent Kohn-Sham (TDKS) equation in real
time and real space yb96 . We have applied our approach to the optical
breakdown ot08 ; ot09 , coherent phonon generation sh10 ; sh12 , high harmonic
generation ot12 , and coupled dynamics of electrons and electromagnetic fields
in a multiscale description ya12 .

In the present paper, we apply the TDDFT to dielectric properties of a medium
excited by short, intense laser pulses. The method is to solve the TDKS
equation in the medium in the presence of an external electromagnetic field
having both pump and probe pulses. Thus, we simulate the pump-probe
experiments numerically. The theory describes the electron dynamics fully
quantum mechanically but assuming that the electrons only interact via a time-
dependent mean field. Thus, the theory is only expected to be justified before
the times when electron-electron collisions have substantially affected the
electronic structure. A separate issue is the creation of phonons. For the
excitation energies we consider here, the electron-electron collision time
sets a more stringent limit than the phonon interactions.

To interpret the results, we compare with several more simplified models. One
model that is often used is based on a Drude response of the excited
quasiparticles, which embedded in a dielectric medium me10 ; re10 . This
requires a number of parameters to be fitted. Another model ignores the
dynamics that created the excited electronic states, replacing it by a thermal
ensemble of electrons rc06 . One can carry out the TDDFT calculation of the
linear response of the thermal system and compare it directly with the
response to the pump-excited system. We will find that many features of
response can be understood even at a quantitative level with the simpler
treatments. However, there are also features that only appear in the full
pump-probe simulation.

The construction of the paper is as follows. In Sec. II, we describe a method
and results of numerical pump-probe simulation to extract dielectric
properties of excited silicon. In Sec. III, we present results of a thermal
model and compare them with the numerical pump-probe results. Our findings are
summarized in Sec. IV.

##  II Numerical Pump-probe experiments

In this section we carry out what we call numerical pump-probe experiments to
study the dielectric properties of the highly excited material. We examine the
electronic response in a unit cell of a crystalline solid irradiated by the
pump and probe laser pulses. Since the wavelength of the laser pulses is much
longer than a typical length scale of electron dynamics, we treat the laser
electric field as a spatially uniform field. The current induced by the probe
pulse will be used to investigate the dielectric properties of excited matter.

###  II.1 Calculation of electron dynamics

Our calculation method has been described in detail elsewhere be00 ; ya06 ;
be06 ; sh10 ; ya12 , so we only provide here the details germane to our study
here. The electrons dynamics is calculated using the TDKS equation,

$\displaystyle i\hbar\frac{\partial}{\partial
t}\psi_{i}(\vec{r},t)=h_{KS}(t)\psi_{i}(\vec{r},t),$ (1)

where $h_{KS}(t)$ is a time-dependent Kohn-Sham Hamiltonian. It differs from
the ordinary TDKS Hamiltonian by inclusion of the vector potential $\vec{A}$
in the kinetic term, $p^{2}/2m\rightarrow(\vec{p}+e\vec{A}/c)^{2}/2m$. There
is also a coupling to $\vec{A}$ in the pseudopotential for the ions; see Ref.
be00 ; ya06  for details. The electron-electron interaction in the TDKS
Hamiltonian is modeled in a simple a adiabatic local density approximation lda
. We calculate dynamics of valence electrons only, treating the electron-ion
interaction by the norm-conserving pseudopotential tm91 ; kb82 .

The current flowing within the unit cell is given by

$\displaystyle\vec{J}(t)=\sum_{i}\frac{1}{\Omega}\Re\left[\int_{\Omega}d\vec{r}\psi_{i}^{*}(\vec{r},t)\vec{j}(t)\psi_{i}(\vec{r},t)\right],$
(2)

where $\Omega$ is the volume of the unit cell and the current operator
$\vec{j}(t)$ is defined by

$\displaystyle\vec{j}(t)=-\frac{e}{m}\frac{1}{i\hbar}\left[\vec{r},h_{KS}(t)\right].$
(3)

The relation of the vector potential $A(t)$ in the unit cell to the external
electromagnetic field exciting the system depends on a number of factors
including possible macroscopic polarization fields. In the present analysis,
we assume a transverse geometry as discussed in Ref. ya12 . The sample is
treated as infinite in the direction of the polarization vector so there
appears no polarization field inside the solid. Of course the field is also
affected by the absorption and the reflection from the surface region, but we
don’t attempt here to express the results in terms of the incident laser
intensity. We take the following form for the vector potential of the pump
pulse in the medium $A_{P}(t)$,

$\displaystyle A_{P}(t)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}-c\frac{E_{0}}{\omega_{P}}\cos{(\omega_{P}t)}\sin^{2}(t/{\tau_{L}})&(0<t<{\tau_{L}})\\\
0&(otherwise),\end{array}\right.$ (6)

where $\omega_{P}$ and $\tau_{L}$ is the average frequency and the time length
of the laser pulse, respectively. $E_{0}$ is the maximum electric field
strength in the medium. This is related to the maximum intensity of the pulse
$I$ by $I_{v}=cE_{0}^{2}/8\pi$ in the vacuum and
$I_{m}=\epsilon^{1/2}cE_{0}^{2}/8\pi$ in the medium. Since the dielectric
function $\epsilon$ is not well-defined in the presence of a strong electric
field, we shall report our results using the field intensity corresponding to
the vacuum relationship.

Our computer code to solve the TDKS equation uses a three-dimensional grid
representation to represent orbital wave functions. The unit cell for the
silicon crystal treated has a length $a=10.26$ a.u. contains eight Si atoms.
The cubic unit cell is discretized into $16^{3}$ grid points. The four valence
electrons of Si atoms beyond the closed (1s2s1p) shells are treated
dynamically. The $k$-space is also discretized into $24^{3}$ grid points. The
time evolution is computed using a fourth-order Taylor expansion of the
operator $\exp(-ih_{KS}(t)\Delta t/\hbar)$ yb96 . We use a time step of
$\Delta t=0.08$ a.u. The number of time step is typically 24,000.

In Fig 1, we show an example of the calculated electron dynamics induced by
the intense pump pulse. Here the frequency of the pump pulse is set to
$\hbar\omega_{P}=1.55$ eV and the duration of the pump pulse is $\tau_{L}=18$
fs. These values will be used in all calculations of this paper. For this
figure, the electric field strength corresponds to an intensity of
$I=3.2\times 10^{12}$ W/cm2.

Panel a) of the figure shows the time profile of the electric field,
$E_{P}(t)=-\frac{\partial}{c\partial t}A_{P}(t)$. Panel b) shows the induced
current, calculated using the time-dependent orbitals in Eq. (2). The average
frequency $\hbar\omega=1.55$ eV is smaller than the direct band gap energy
(2.4 eV in the present calculation), so the initial current response is
nondissipative. This is seen by the phase difference between the current and
the electric field, which is shifted by $\pi/2$ at the beginning of the field
pulse ($t<5$ fs). As the intensity of the pulse increases, the system absorbs
energy by the excitation of electron-hole pairs. As a result, the phase
difference decreases. Note that the current shows a weak oscillation at high
frequency after the pulse has passed. Making a Fourier analysis, we find that
it is dominated by frequencies around 3.9 eV/$\hbar$. However, we have no
physical explanation of the oscillation.

Figure 1 (c) shows the excitation energy per Si atom. During the field pulse,
a rapid increase of the electronic excitation energy is seen. After the pulse
ends, the excitation energy is independent of time, showing that our
computational algorithm conserves the energy of the system.

Figure 1:  (a) The time profile of the applied electric field. (b) The current
induced by the applied electric field. (c) The excitation energy per Si atom.

We next discuss the number density of excited electron-hole pairs and the
electronic excitation energy when field pulses of different intensities
irradiate on the silicon crystal. To calculate the number density of excited
electron-hole pairs, we first define eigenstates of the Kohn-Sham Hamiltonian
of the excited system. We consider a system at a time $t_{f}$ sufficiently
after the applied field pulse ends, and denote the Kohn-Sham Hamiltonian of
Eq. (1) at the time $t_{f}$ as $h_{KS}(t_{f})$. We introduce the Kohn-Sham
orbitals which satisfy the following eigenvalue equations.

$\displaystyle\hat{h}_{KS}(t_{f})\phi_{i}^{t_{f}}=\epsilon_{i}^{t_{f}}\phi_{i}^{t_{f}}$
(8)

Note that the Kohn-Sham Hamiltonian $\hat{h}_{KS}(t_{f})$ is different from
that of the initial state, due to the change in electron density in the
excited system. Using eigenfunctions $\phi_{i}^{t_{f}}$, we may define the
number density of electron-hole pairs by

$\displaystyle
n_{eh}=\sum_{i}\left\\{1-\sum_{j}\left|\langle\phi_{j}^{t_{f}}|\psi_{i}(t_{f})\rangle\right|^{2}\right\\},$
(9)

where the sum over $i,j$ is taken over occupied orbitals, and
$|\psi_{i}(t_{f})\rangle$ is the orbital of the TDKS equation at the time
$t_{f}$. We use the final-state definition of the orbitals because it
facilitates the comparison to a thermal model that will be discussed later.

In Fig. 2, we show the number density of electron-hole pairs (top panel) and
the electronic excitation energy per atom (middle panel). The bottom panel
shows the electronic excitation energy divided by the number of electron-hole
pairs. As seen from the figure, both excitation energy and the number of
excited electrons increase with increasing the applied field intensity. At low
intensity region, they scale with the square of the field intensity. This is
because two photons are required for electrons to be excited across the direct
band gap. As seen from the bottom panel, the electronic excitation energy per
excited electron-hole pairs is given by 3.1 eV, which coincides with the two-
photon energy of the applied field pulse. As the field intensity increases
above $10^{11}$ W/cm2, both the number density of electron-hole pairs and the
electronic excitation energy deviated from the two-photon curve. As will be
seen later, the dielectric property of excited matter also shows a large
change from that in the ground state at field intensities above this value.

Figure 2:  Electronic properties of the crystalline silicon in the final state
following the pulsed excitation. Shown are: electronic excitation energy per
Si atom (top panel); the number density of electron-hole pairs (middle panel);
the excitation energy per electron-hole pair (bottom panel). Quantities are
shown as a function of the maximum pump intensity determined as
$I=cE_{0}^{2}/8\pi$.

###  II.2 Dielectric function from numerical pump-probe calculation

To extract dielectric properties of excited matter, we compare two
calculations, one solving the TDKS equation of Eq. (1) with the vector
potential containing pump and probe pulses and the other containing the pump
pulse only. We denote the electric field of the pump pulse as $E_{P}(t)$ and
that of the probe pulse as $E_{p}(t)$. We denote the current in the
calculations containing the pump pulse only as $J_{P}(t)$ and that in the
calculations containing both pump and probe pulses as $J_{Pp}(t)$. We define
the current induced by the probe pulse as the difference,

$\displaystyle J_{p}(t)=J_{Pp}(t)-J_{P}(t).$ (10)

From the difference of the induced currents, we may extract the conductivity
and the dielectric function of excited matter.

From the probe current $J_{p}(t)$, we may extract the electric conductivity
$\sigma(\omega)$ and the dielectric function $\epsilon(\omega)$ of excited
matter by the following equations:

$\displaystyle\sigma(\omega)=\frac{\int dtJ_{p}(t)e^{i\omega t}}{\int
dtE_{p}(t)e^{i\omega t}}$ (11) $\displaystyle\epsilon(\omega)=1+\frac{4\pi
i\sigma(\omega)}{\omega},$ (12)

In principle, the above-defined conductivity and dielectric function depend
also on the time delay $\tau_{Pp}$ between the pump and probe pulses. We will
later show that the dependence on delay-time is rather small in the TDDFT
calculations.

In practice we employ the vector potential of the form of Eq. (6) as the pump
pulse. As for the probe pulse, we use the same functional form as Eq. (6)
delayed by an amount $\tau_{Pp}$ from the pump pulse,

$\displaystyle A_{p}(t)$ $\displaystyle=$
$\displaystyle-c\frac{e_{0}}{\omega_{p}}\cos{(\omega_{p}(t-\tau_{Pp}))}$ (13)
$\displaystyle\times\sin^{2}((t-\tau_{Pp})/\tau_{L})$

for $\tau_{Pp}<t<\tau_{L}+\tau_{Pp}$ and zero otherwise.

In Fig. 3, we show typical time profiles of the electric fields and the
induced currents for a delay time of $\tau_{Pp}=19$ fs. The pump pulse is the
same as in Fig. 1, with a maximum intensity of $3.2\times 10^{11}$ W/cm2. The
probe intensity is a factor of 200 smaller, which we deem to be sufficiently
weak to extract the linear response. In the left panels of Fig. 3, we show
electric fields of pump and probe pulses, $E_{P}(t)+E_{p}(t)$, in (a), pump
pulse, $E_{P}(t)$, in (b), and probe pulse, $E_{p}(t)$, in (c), as functions
of the time. The right panels show currents induced by the pump and probe
pulses, $J_{Pp}(t)$, in (d), by the pump pulse only, $J_{P}(t)$, in (e), and
the difference of the currents, $J_{p}(t)$ of Eq. (10), in (f).

Figure 3:  Left panels show the electric field of pump and probe pulses in
(a), pump pulse in (b), and the probe pulse in (c). Right panels show the
current induced by the pump plus probe pulse in (d), the current by the pump
pulse only in (e), and the difference of the currents shown in (d) and (e) in
(f).

The next step is to calculate the dielectric function from the probe current
using Eqs. (11) and (12). The pump-probe calculation using the probe pulse of
Eq. (13) and the probe current of Eq. (10) gives us dielectric properties
around the average frequency $\hbar\omega_{p}$. To investigate dielectric
properties for a wide frequency region, we repeat the pump-probe calculations
for different frequencies of the probe pulses.

In Fig. 4, we show typical calculations using a number of probe pulses of
differing frequencies. Panels (a) and (b) show the absolute values of the
Fourier transforms of $E_{p}(t)$ and $J_{p}(t)$, respectively. Panel (c) shows
the real part of the dielectric function which is calculated using Eqs. (11)
and (12). The curve is composed of a number of curves with different colors
for each probe frequency. One can see that the overlap is very good for the
different average probe frequencies, validating our method to extract the
dielectric function.

We next ask how sensitive is the extracted dielectric function to the time
delay of the probe pulse? Since there are no dissipative processes in the
TDDFT under ALDA, the properties of the system should not change after some
initial period when the phases of the excited orbitals become incoherent.
Figure 5 shows how the extracted dielectric function depends on delay time for
one of the cases. We have selected delay times over a range that corresponds
to a full cycle of the pump pulse, since that frequency could be imprinted on
the phases of the particle at later times. The range of the delay times is
19.00 fs, 19.67 fs, 20.33fs, and 21.67fs. The latter three delay times
correspond to a quarter, a half, and one period of the pump pulse
$2\pi/\omega_{P}$ added to the first time. One can see that real part is
practically independent of the delay. The imaginary part, however, shows
variation although qualitatively the functions are similar. We found the same
behavior for other cases as well. Namely, the real part is independent of
delay, even extending the delay to very large times. The imaginary part is
only qualitatively similar for different delay times. In the sequel, we will
analyze all the results using the dielectric function at $\tau_{Pp}=19$ fs,
and one should remember that the imaginary part is less well defined than the
real part.

Figure 4:  The top two panels show the Fourier transformations of the probe
electric field $E_{p}(\omega)$ and the probe current $J_{p}(\omega)$. The
bottom panel shows the real part of the dielectric function extracted from
$E_{p}(\omega)$ and $J_{p}(\omega)$ through Eqs. (11) and (12). The pump pulse
has an intensity $I=3.2\times 10^{12}$ W/cm2 and an average frequency
$\omega_{P}=1.55$ eV/$\hbar$. The polarization directions of the pump and
probe pulses are taken to be parallel.

Figure 5:  Extracted dielectric functions as a function of the delay time
$\tau_{Pp}$. The pump pulse has an intensity $I=3.2\times 10^{12}$ W/cm2 and
an average frequency $\omega_{P}=1.55$ eV/$\hbar$. The probe frequency is
$\omega_{p}=1.55$ eV/$\hbar$ and its delay times for the four graphed lines
are $\tau_{Pp}=19,0,19.67,20.33$ and $21.67$ fs. The polarization directions
of the pump and probe pulses are taken to be parallel.  Figure 6:  Real (a)
and imaginary (b) parts of the dielectric functions of Si excited by field
pulses of three intensities. The dielectric functions are deduced using pump-
probe calculations. The polarization direction of the probe pulse is taken to
be parallel to that of pump pulse. The dielectric function of silicon in the
ground state is also shown.

We have carried out the pump-probe simulations for several intensities of the
pump pulse. The results for the dielectric functions are shown in Fig 6. The
real and the imaginary parts are presented in panels (a) and (b),
respectively.

The distinguishing feature in the response is the negative divergence at small
frequencies. This arises from the quasiparticles in the excited system, as we
will see more quantitatively later. The imaginary part of the response is not
quite as simple to analyze. The plasmon in a free electron gas at zero
temperature is undamped and the imaginary part of the dielectric function
vanishes at low frequency. Figure 6 shows, however, that ${\rm
Im}[\epsilon(\omega)]$ becomes large at low frequency. The quasiparticle
response is thus far from that of a free electron gas.

An interesting feature of the TDDFT response is that the dielectric tensor is
not isotropic in the excited crystal, even though the crystal symmetry is
cubic. This may be seen in Fig. 7, comparing the dielectric functions for the
probe polarization either parallel or perpendicular to the pump polarization.

Figure 7:  Comparison of dielectric functions of excited silicon probed with
either parallel or perpendicular direction to the pump fields.

The real part of the dielectric function shows the low-frequency plasmon
response more strongly for the parallel component. One may notice that ${\rm
Im}[\epsilon(\omega)]$ is negative at some frequencies. This might indicate a
population inversion that could sustain a growth of intensity at those
frequencies. However, one should carry the full calculation of the pulse
propagating in space as well as time to assert that the excited medium can
amply the pulses.

###  II.3 Comparison with free-carrier models

Dielectric properties of solids excited by intense and ultrashort laser pulses
are often modeled employing a simplified dielectric function, adding a Drude-
like component to the dielectric function in the ground state me10 ; re10 . In
this subsection, we will examine how well the dielectric function of excited
matter in the TDDFT calculation may be described by a simplified dielectric
model.

First we consider an embedded Drude model, the dielectric function given as a
sum of the ground state response and the Drude response of free carriers

$\displaystyle\epsilon_{ED}(\omega)=\epsilon_{0}(\omega)-4\pi
i\frac{e^{2}n_{eh}}{m^{*}\omega(\omega+i/\tau)}.$ (14)

Here $\epsilon_{0}(\omega)$ is the dielectric function in the ground state,
$n_{eh}$ is the electron-hole density, $m^{*}$ is the reduced mass of
electron-holes, and $\tau$ is the Drude damping time. For the dielectric
function in the ground state, $\epsilon_{0}(\omega)$, we will use the values
obtain from the TDDFT calculation. The number density of electron-hole pairs,
$n_{eh}$, is extracted from the calculation using Eq. (9). We treat $m^{*}$
and $\tau$ as parameters, fitting to the calculated $\epsilon(\omega)$.

Sokolowski-Tinten and von der Lind proposed a more complicated model for the
dielectric function excited by strong laser fields so00 , which we shall call
the SL model. They consider three physical effects for the dielectric response
of laser-excited semiconductor: (i) state and band filling, (ii)
renormalization of the band structure, and (iii) the free-carrier response.
The SL dielectric function is parameterized as

$\displaystyle\epsilon_{SL}(\omega)$ $\displaystyle=$ $\displaystyle
1+\left[\epsilon_{0}(\omega+\Delta
E_{gap})-1\right]\frac{n_{0}-n_{eh}}{n_{0}}$ (15) $\displaystyle-4\pi
i\frac{e^{2}n_{eh}}{m^{*}\omega(\omega+i/\tau)},$

where $\Delta E_{gap}$ is the change of the band gap by the laser irradiation
and $n_{0}$ is the density of electrons which contribute to the dielectric
response. For $\Delta E_{gap}$, we use a change of single particle energies,
$\epsilon_{i}^{t_{f}}$ of Eq. (8), after the laser pulse ended. We treat the
active number of valence electrons, $n_{0}$, as a fitting parameter.

Figure 8:  Dielectric function of silicon excited by two field pulses,
$3.2\times 10^{11}{\rm W/cm^{2}}$ (left panels) and $3.2\times 10^{12}{\rm
W/cm^{2}}$ (right panels), are fitted by the embedded Drude model and by the
SL model.

Figure 8 shows the fits obtained in the embedded Drude model and the SL model
for two intensities of the pump field. The polarization directions of the pump
and probe pulses are taken to be parallel. At the lower intensity, the
dielectric functions are well fitted by both models. At the higher intensity,
the real part is well described by both models at lower frequencies, but the
SL model fits better above the direct band gap. However, neither model does
well for the imaginary part in the high intensity case.

In the fitting procedure, we found the effective mass is sensitive to the real
part of the dielectric function and can be determined without ambiguity. The
effective mass for the pump pulse of $3.2\times 10^{11}$ W/cm2 is given by
$m^{*}=0.18$, while the effective mass for the pump pulse of $3.2\times
10^{12}$ W/cm2 is given by $m^{*}=0.32$. The collision time for the weaker
pump pulse case of $3.2\times 10^{11}$ W/cm2 is determined to be about 6 fs.
However, the Drude damping time $\tau$ cannot be determined well for high
intensity case, because the frequency dependence of the imaginary-part
dielectric function is very different from the Drude behavior.

The effective mass and its change with excitation energy may be understood
from the band structure. A weak pump pulse excites electrons at specific
$k$-points by two-photon absorption, while a strong pump pulse excites
electrons at various $k$-points by tunnel and multi-photon excitations. The
effective mass of electrons depends very much on the their positions in the
bands; only the lowest bands have the very small effective masses.

##  III Thermal model

The numerical pump-probe experiments reported in the previous section are
applicable to the excited matter immediately after the pump irradiation,
perhaps for a time period of a few tens of femtoseconds. In this section, we
will investigate dielectric properties of thermally excited matter, which
should be more appropriate at later times. We assume that the electronic
states are described by a thermal ensemble, but the ions have not yet had time
to respond.

###  III.1 Dielectric function at finite temperature

We describe the thermally excited matter by static density functional theory
at finite electron temperature. In the calculations, the electrons population
is described by the Fermi-Dirac distribution of a given temperature, while
atomic positions are frozen at the ground-state positions. The electron
density at temperature $T$, $\rho^{T}(\vec{r})$, is given by

$\displaystyle\rho^{T}(\vec{r})=\sum_{i}n_{i}^{T}|\phi_{i}(\vec{r})|^{2},$
(16)

where $n_{i}^{T}$ is the temperature-dependent occupation number of Fermi-
Dirac distribution,

$\displaystyle n_{i}^{T}=\frac{1}{1+{\rm e}^{(\epsilon_{i}-\mu)/k_{B}T}},$
(17)

Here $\epsilon_{i}$ is the single particle energy, $\mu$ is the chemical
potential, and $k_{B}T$ is the temperature in energy units. We note that all
the quantities related to the orbitals, $\phi_{i}$, $\epsilon_{i}$, and $\mu$
depend on the temperature $T$ due to the self-consistency requirement.

In Fig. 9, we show calculated results at several electron temperatures. The
top panel (a) shows the excitation energy per atom, the middle panel (b) shows
the number of excited electrons, and the bottom panel (c) shows the excitation
energy per excited electron. As seen from the figure, both excitation energy
and the number of excited electrons monotonically increase as the electron
temperature increases. At very low temperature, we find the excitation energy
per excited electrons is rather small, 1.2 eV for $k_{B}T=0.2$ eV. This value
should approach to the energy of the indirect gap, 0.52 eV in our calculation,
in the low temperature limit.

Figure 9:  Basic electron properties of the thermal model of crystalline
silicon. The top panel shows the excitation energy per Si atom. The middle
panel shows the number of excited electron per Si atom. The bottom panel shows
the excitation energy per excited electron.

We investigate dielectric properties of crystalline silicon at finite
temperature by the real-time method, applying a distorting vector potential of
step function in time be00 . In Fig. 10, we show dielectric functions of
silicon at several electron temperatures. As seen from the real-part of the
dielectric function, all responses at finite temperatures show a Drude-like
like behavior at low frequencies. This behavior is more or less similar to
those in our numerical pump-probe calculations shown in Fig. 6. The low energy
component of the imaginary part shows absorptive contributions at low
frequencies, increasing monotonically as the temperature increases.

Figure 10:  The dielectric functions of the finite temperature model at
several temperatures. Top panel shows the real part of the dielectric
function, and the bottom panel shows the imaginary part.

A good way to display a plasmon contribution to the response is to plot the
imaginary part of the inverse dielectric function, ${\rm Im}\epsilon^{-1}$.
This is shown in Fig. 11 for several temperatures up to $k_{B}T=1.2$ eV. At
the lowest temperature one sees a very sharp plasmon peak, located at an
energy of $\sim 0.4$ eV. The plasmon excitation energy increases with
temperature, due to the increased density of electron-hole pairs. We note that
the width of the plasmon also increases with temperature, up to about
$k_{B}T\approx 0.6$ eV. Beyond that, the width does not change very much, up
to the maximum temperature considered.

Figure 11:  The imaginary part of the inverse dielectric function for various
electronic temperatures.

The thermal dielectric function presented here can be well fitted by the SL
model. For the basic parameters of the SL model, we employ the calculated
dielectric function at zero temperature for $\epsilon_{0}$, the calculated
electron-hole pair density for $n_{eh}$, and the calculated shift of the gap
energy for $\Delta E_{gap}$. Other three parameters, $m^{*}$, $\tau$, and
$n_{0}$ are treated as fitting parameters. The fit is carried out by
minimizing the mean square error as given by

$\displaystyle
I_{error}=\int_{\omega_{i}}^{\omega_{f}}d\omega\left|\epsilon_{T}^{-1}(\omega)-\epsilon_{SL}^{-1}(\omega)\right|^{2},$
(18)

in the interval $\hbar\omega_{i}=0.35$ eV and $\hbar\omega_{f}=5.0$ eV. The
$\epsilon_{T}(\omega)$ is the dielectric function in the thermal model. The
quality of the fit is shown in Fig. 12 for a temperature of $k_{B}T=1.2$ eV in
the thermal model. The fit is very good except for the ${\rm Im}\epsilon$ at
the lowest frequencies. In particular, the plasmon feature in the inverse
dielectric function is very well reproduced.

Figure 12:  Comparison of the thermal model and a fit with the SL model. The
electronic temperature in the thermal model is $k_{B}T=1.2$ eV.

In Fig. 13, we show the fitted effective mass $m^{*}$ and the collision time
$\tau$ as functions of the temperature in the thermal model. The top panel
shows that the effective mass $m^{*}$ gets heavier as the temperature
increases. This feature was also found for the dielectric function extracted
from the numerical pump-probe experiments, see Fig. 8. As we discussed in Sec.
II, the change of effective mass may be understood by the change of the
distribution of the electron-hole quasiparticles in $k$-space.

The bottom panel of Fig. 13 shows that the damping time $\tau$ becomes very
small as the electron temperature increases. The value of $\tau$ reaches a
saturated value of 1.2 fs at $k_{B}T\approx 0.8$ eV. At first sight this is
puzzling, because there are no collision effects in either the TDKS equation
or in the thermal model. In spite of this fact, we obtain a plasmon feature
with large damping, corresponding to collision times as short as 1.2 fs in the
thermal model. We believe that the damping arises from the elastic scattering
of high-energy quasiparticles from ionic core potentials, but we have no
quantitative understanding of its magnitude or dependence on the quasiparticle
distribution. We note that TDDFT treatment of linear response describes the
dielectric function of metals fairly well, including the width of plasmon seen
in the inverse dielectric function be00 .

Figure 13:  Parameters of the SL model determined by the fitting procedure to
the thermal model. Top panel shows the effective mass $m^{*}$ and the bottom
panel shows the Drude damping time $\tau$.

###  III.2 Comparison with numerical pump-probe experiments

The difference between the numerical pump-probe calculations presented in Sec.
II and the thermal model comes entirely from the different electron-hole
distributions in the excited system to be probed. In this subsection, we
compare their predicted dielectric functions.

We first need to assume a correspondence between the excited systems that we
wish to compare. Two possibilities come to mind, namely consider systems of
equal excitation energy or of equal densities of particle-hole excitations. In
general, the laser-excited system will have a higher excitation energy for the
same number of particle-hole pairs. Since the plasmon characteristics are
closely tied to the number of quasiparticles, we shall use that measure to
make the comparison. One comparison will be with $n_{ph}=0.015$ /Atom; this is
obtained by a pump pulse of $3.2\times 10^{11}$ W/cm2 or by a thermal system
with electron temperature $k_{B}T=0.3$ eV. Another comparison will be with
$n_{ph}=0.3$ /Atom, requiring a pump pulse of $3.2\times 10^{12}$ W/cm2 or a
temperature of $k_{B}T=1.3$ eV. The two dielectric functions are shown in Fig.
14.

Figure 14:  Comparisons of the dielectric function of the numerical pump-probe
calculation and the thermal model. Left-hand panels: $n_{eh}=0.015$ /Atom;
right-hand panels: $n_{eh}=0.3$ /Atom.

As seen from the figure, the dielectric function of silicon excited by the
low-intensity pump-pulse ($I=3.2\times 10^{11}$ W/cm2) is close to the thermal
model. At the stronger pump intensity of $I=3.2\times 10^{12}$ W/cm2, the real
part of the dielectric function in the thermal model lies between the
perpendicular and parallel dielectric functions of the numerically pumped
system. In fact, the finite temperature function is close to the parallel case
at the higher frequencies. On the other hand, the imaginary part of the
dielectric function agrees with the perpendicular case at high frequencies,
but is much larger than either polarizations at low frequency. We thus
conclude that, making correspondence between the numerical pump-probe
calculations and the thermal electron model using the number density of
electron-hole pairs, the dielectric functions show a reasonable
correspondence, especially when the excitations are not very violent. The
difference between two calculations comes from the different $k$-space
distributions of electron-hole pairs. It seems that the difference is more
evident for the imaginary part.

##  IV SUMMARY

To investigate a change of dielectric properties of bulk silicon immediately
after ultrashort laser pulses, we made numerical pump-probe experiments
solving the time-dependent Kohn-Sham equation in real time including electric
fields of both pump and probe pulses. The simulation makes it possible to
investigate dielectric properties of excited matter before any dissipation or
dephasing effects start to become significant.

For a comparison, we have also constructed a thermal model by solving the
static Kohn-Sham equation with finite-temperature Fermi function occupation
factors. Its dielectric response was then computed by applying the linear
response theory using the usual real-time method. We found that the thermal
model works very well at lower excitation energies, but becomes unreliable at
the higher excitation energy where one finds differences in the parallel and
perpendicular dielectric functions.

An even more simple model can be constructed using ingredients of the Drude
model of free-electron dynamics. In general, the real part of the dielectric
function was found to be well described by a Drude-like contribution of the
excited quasiparticles embedded in the dielectric medium corresponding to the
ground state. As for the imaginary part, the dielectric function in the
thermal model is reasonably described by the embedded Drude model. The
dielectric function in the numerical pump-probe experiments shows rather
different behavior, even negative values for the imaginary part. The
difference comes from the nonequilibrium distributions of electrons and holes
in $k$-space in the numerical pump-probe experiments.

In the embedded Drude model, there are three parameters determining the
quasiparticle plasmon contribution, namely the density of quasiparticles,
their effective mass $m^{*}$, and the collision time $\tau$. The density of
quasiparticles is known from the TDDFT or thermal calculation, but the other
quantities are fit. From the real part of the dielectric function, we find
increase of the effective mass as the pump field intensity increases as
expected from the band structure. We find the calculated dielectric functions
show substantial imaginary part. In the thermal electron model, the collision
time of as short as 1.2 fs gives reasonable fit. This short value for the
collision time is unexpected, since there are no explicit collision terms in
the time-dependent Kohn-Sham equation that we solve. We believe the short
collision time comes from the elastic scattering of electrons from atoms
rather than from electron-electron or electron-phonon interactions, but we
lack a simple model to exhibit this aspect of the response.

## Acknowledgments

This work is supported by the Grants-in-Aid for Scientific Research No.
23340113, No. 23104503, No. 21340073, and No. 21740303. The numerical
calculations were performed on the supercomputer at the Institute of Solid
State Physics, University of Tokyo, and T2K-Tsukuba at the Center for
Computational Sciences, University of Tsukuba. G.F.B. acknowledges support by
the National Science Foundation under Grant No. PHY-0835543 and by the US
Department of Energy under Grant No. DE-FG02-00ER41132.

## References

  * (1) See, e.g. the Conference proceedings of http://lap2012.sciencesconf.org.
  * (2) M.D. Perry, B.C. Stuart, P.S. Banks, M.D. Feit, V. Yanovsky, A.M. Rubenchik, J. Appl. Phys. 85, 6803 (1999).
  * (3) T. Brabec, F. Krausz, Rev. Mod. Phys. 72, 545 (2000).
  * (4) A. Couairon and A. Mysyrowicz, Phys. Rep. 441, 47 (2007).
  * (5) R. Merlin, Solid State Commun. 102, 207 (1997).
  * (6) K. Sokolowski-Tinten and D. von der Linde, Phys. Rev. B 61, 2643 (2000).
  * (7) M. Schultze, E.M. Bothschafter, A. Sommer, S. Holzner, W. Schweinberger, M. Fiess, M. Hofstetter, R. Kienberger, V. Apalkov, V.S. Yakovlev, M.I. Stockman, and F. Krausz, Nature 493, 75 (2013).
  * (8) A. Schiffrin, T. Paasch-Colberg, N. Karpowicz, V. Apalkov, D. Gerster, S. Mühlbrandt, M. Korbman, J. Reichert, M. Schultze, S. Holzner, J. Barth, R. Kienberger, R. Ernstorfer, V. Yakovlev, M. I. Stockman, and F. Krausz, Nature 493, 70 (2013).
  * (9) A.V. Mitrofanov, A.J. Verhoef, E.E. Serebryannikov, J. Lumeau, L. Glebov, A.M. Zheltikov, A. Baltuska, Phys. Rev. Lett. 106, 147401 (2011).
  * (10) A. Kaiser, B. Rethfeld, M. Vicanek, and G. Simon, Phys. Rev. B61, 11437 (2000).
  * (11) B. Rethfeld, Phys. Rev. Lett. 92, 187401 (2004).
  * (12) N. Medvedev and B. Rethfeld, J. Appl. Phys. 108, 103112 (2010)
  * (13) B. Rethfeld, O. Brenk, N. Medvedev, H. Krutsch, D.H.H. Hoffmann, Appl. Phys. A101 19, (2010).
  * (14) E. Runge and E. K. U. Gross, Phys. Rev. Lett. 52, 997 (1984).
  * (15) K. Yabana and G.F. Bertsch, Phys. Rev. B54, 4484 (1996).
  * (16) T. Otobe, M. Yamagiwa, J.-I. Iwata, K. Yabana, T. Nakatsukasa, and G.F. Bertsch, Phys. Rev. B77, 165104 (2008).
  * (17) T. Otobe, K. Yabana, and J.-I. Iwata, J. Phys. Cond. Matter 21, 064224 (2009).
  * (18) Y. Shinohara, K. Yabana, Y. Kawashita, J.-I. Iwata, T. Otobe, G.F. Bertsch, Phys. Rev. B82 155110 (2010).
  * (19) Y. Shinohara, S.A. Sato, K. Yabana, T. Otobe, J.-I. Iwata, G.F. Bertsch, J. Chem. Phys. 137, 22A527 (2012).
  * (20) T. Otobe, J. Appl. Phys. 111, 093112 (2012).
  * (21) K. Yabana, T. Sugiyama, Y. Shinohara, T. Otobe, G.F. Bertsch, Phys. Rev. B 85 045134 (2012).
  * (22) V. Recoules, J. Clérouin, G. Zérah, P.M. Anglade, S. Mazevet, Phys. Rev. Lett. 96, 055503 (2006).
  * (23) G.F. Bertsch, J.-I. Iwata, A. Rubio, and K. Yabana, Phys. Rev. B 62 7998 (2000).
  * (24) K. Yabana, T. Nakatsukasa, J.-I. Iwata, and G.F. Bertsch, phys. stat. sol. (b)243, 1121 (2006).
  * (25) G.F. Bertsch and K. Yabana, Introduction to Computational Methods in Many Body Physics, Chap. 3, eds. M. Bonitz and D. Semkat, Rinton Press 2006.
  * (26) J.P. Perdew and A. Zunger, Phys. Rev. B23, 5048 (1981).
  * (27) N. Troullier and J.L. Martins, Phys. Rev. B43, 1993 (1991).
  * (28) L. Kleinman and D. M. Bylander, Phys. Rev. Lett. 48, 1425 (1982).


††thanks: They contribute equally to this work. ††thanks: They contribute
equally to this work.

# Rényi Entropy Dynamics and Lindblad Spectrum for Open Quantum System

Yi-Neng Zhou  Institute for Advanced Study, Tsinghua University, Beijing
100084, China  Liang Mao  Institute for Advanced Study, Tsinghua University,
Beijing 100084, China  Department of Physics, Tsinghua University, Beijing
100084, China  Hui Zhai<EMAIL_ADDRESS>Institute for Advanced Study,
Tsinghua University, Beijing 100084, China

###### Abstract

In this letter we point out that the Lindblad spectrum of a quantum many-body
system displays a segment structure and exhibits two different energy scales
in the strong dissipation regime. One energy scale determines the separation
between different segments, being proportional to the dissipation strength,
and the other energy scale determines the broadening of each segment, being
inversely proportional to the dissipation strength. Ultilizing a relation
between the dynamics of the second Rényi entropy and the Lindblad spectrum, we
show that these two energy scales respectively determine the short- and the
long-time dynamics of the second Rényi entropy starting from a generic initial
state. This gives rise to opposite behaviors, that is, as the dissipation
strength increases, the short-time dynamics becomes faster and the long-time
dynamics becomes slower. We also interpret the quantum Zeno effect as specific
initial states that only occupy the Lindblad spectrum around zero, for which
only the broadening energy scale of the Lindblad spectrum matters and gives
rise to suppressed dynamics with stronger dissipation. We illustrate our
theory with two concrete models that can be experimentally verified.

For a closed quantum system, the energy spectrums of Hamiltonian fully
determine the time scales of its dynamics. For an open quantum system, when
the environment is treated by the Markovian approximation, the couplings
between system and environment are controlled by a set of dissipation
operators. In this case, the dynamics of the system is governed by the
Lindblad equation which contains the contributions from both the Hamiltonian
and the dissipation operators open . Obviously, the spectrum of the
Hamiltonian alone can no longer determine the time scales of the entire
dynamics, and a natural question is then what energy scales set the time
scales of dynamics of an open quantum system.

There are various directions to approach this issue, and the answer also
relies on what type of dynamics that we are concerned with. Here let us focus
on the dissipation driven dynamics. There are still different physical
intuitions from different perspectives. One intuition is from the perturbation
theory when the dissipation strength is weaker compared with the typical
energy scales of the Hamiltonian Pan . In this regime, by treating the
dissipation perturbatively, it leads to a scenario that the dissipation
dynamics becomes faster when the dissipation strength is stronger. Another
intuition is from the studies of the quantum Zeno effect Zeno_paradox ;
Zeno1990 ; quantum_Zeno ; Zeno_review , which states that frequent
measurements can slow down the dynamics, provided that the typical time
interval between two successive measurements are shorter than the intrinsic
time scale of the system. Since the measurement can also be understood in term
of dissipations in the Lindblad master equation, it provides another scenario
that the dissipation dynamics is suppressed when the dissipation becomes
stronger, in the regime that the dissipation strength is stronger compared
with the typical energy scales of the Hamiltonian. It seems that these two
scenarios respectively work on different parameter regimes and the results are
also opposite to each other. It will be interesting to see that there actually
exists a framework that can unify these two scenarios.

When a system is coupled to a Markovian environment, the entropy of the system
will increase in time. The entropy dynamics of an open quantum many-body
system is a subject that attracts lots of interests recently entang_review ;
Qi_Cricuit ; Chen ; Kitaev ; YuChen ; Zhai . In this letter, we address the
issue of typical time scales of the entropy increasing dynamics of a quantum
many-body system coupled to a Markovian invironment, and especially, we should
focus on the second Rényi entropy, for the reason that will be clear below,
and answer the question whether the entropy dynamics is faster or slower when
the dissipation strength increases.

Figure 1: Schematic of the mapping between the Lindblad equation (left) and
the Schödinger like equation in a doubled system (right). Here
$\hat{L}\hat{\rho}$ denotes the r.h.s. of Eq. 1.

Our studies are based on a mapping between the Lindblad master equation and a
non-unitary evolution of wave function in a doubled space, as shown in Fig. 1.
Let us first review this mapping Operator_Schmidt ; Mixed_state . Considering
a density matrix $\hat{\rho}$, and given a set of complete bases
$\\{|n\rangle\\},(n=1,\dots,\mathcal{D}_{\text{H}})$ of the Hilbert space with
dimension $\mathcal{D}_{\text{H}}$ (say, the eigenstates of the Hamiltonian
$\hat{H}$ with eigenenergies $E_{n}$), the density matrix $\hat{\rho}$ can be
expressed as $\hat{\rho}=\sum_{mn}\rho_{mn}|m\rangle\langle n|$. By the
operator-to-state mapping, we can construct a wave function
$\Psi_{\rho}=\sum_{mn}\rho_{mn}|m\rangle\otimes|n\rangle$, which contains
exact the same amount information as $\hat{\rho}$. Here $\Psi_{\rho}$ is a
wave function on a system whose size is doubled compared to the original
system, and we will refer these two copies of original system as the “left”
(L) and the “right” (R) systems. Under this mapping, for instance, a density
matrix of a pure state $\hat{\rho}=|\psi\rangle\langle\psi|$ is mapped to a
product state $\Psi_{\rho}=|\psi\rangle\otimes|\psi\rangle$ in the double
system, and a thermal density matrix at temperature $T$ as
$\hat{\rho}=\sum_{n}e^{-E_{n}/(k_{\text{b}}T)}|n\rangle\langle n|$ is mapped
to a thermofield double state at temperature $T/2$ as $\Psi_{\rho}=\sum
e^{-E_{n}/(k_{\text{b}}T)}|n\rangle\otimes|n\rangle$ in the double system.

For an open system coupled to a Markovian environment, the density matrix
obeys the Lindblad master equation given by

$\hbar\frac{d\hat{\rho}}{dt}=-i[\hat{H},\hat{\rho}]+\sum\limits_{\mu}\gamma_{\mu}\left(2\hat{L}_{\mu}\hat{\rho}\hat{L}_{\mu}^{\dagger}-\\{\hat{L}^{\dagger}_{\mu}\hat{L}_{\mu},\hat{\rho}\\}\right),$
(1)

where $\hat{L}_{\mu}$ stand for a set of dissipation operators, and
$\gamma_{\mu}$ are their corresponding dissipation strengths. After the
mapping, the wave function $\Psi_{\rho}$ in the double system satisfies a
Schrödinger-like equation

$i\hbar\frac{d\Psi_{\rho}}{dt}=\left(\hat{H}_{\text{s}}-i\hat{H}_{\text{d}}\right)\Psi_{\rho}.$
(2)

Here $\hat{H}_{\text{s}}$ is the Hermitian part of the Hamiltonian determined
by system itself, and it is given by

$\hat{H}_{\text{s}}=\hat{H}_{\text{L}}\otimes\hat{I}_{\text{R}}-\hat{I}_{\text{L}}\otimes\hat{H}^{\text{T}}_{\text{R}},$
(3)

where operators with subscript “L” and “R” respectively stand for operators
acting on the left and the right systems, and “T” stands for the transpose,
and $\hat{I}$ represents the identity operator. $-i\hat{H}_{\text{d}}$ is the
non-Hermitian part of the Hamiltonian determined by the dissipation operators,
which is given by

$\displaystyle\hat{H}_{\text{d}}=\sum_{\mu}\gamma_{\mu}$
$\displaystyle\left[-2\hat{L}_{\mu,\text{L}}\otimes\hat{L}^{\text{*}}_{\mu,\text{R}}\right.$
$\displaystyle\left.+(\hat{L}^{\dagger}_{\mu}\hat{L}_{\mu})_{\text{L}}\otimes\hat{I}_{\text{R}}+\hat{I}_{\text{L}}\otimes(\hat{L}^{\dagger}_{\mu}\hat{L}_{\mu})^{\text{*}}_{\text{R}}\right],$
(4)

where the superscript * stands for taking complex conjugation. We can
diagnolize this non-Hermitian Hamiltonian
$\hat{H}_{\text{s}}-i\hat{H}_{\text{d}}$, which leads to a set of eigenstates
as

$(\hat{H}_{\text{s}}-i\hat{H}_{\text{d}})|\Psi^{l}_{\rho}\rangle=\epsilon_{l}|\Psi^{l}_{\rho}\rangle,$
(5)

where $\epsilon_{l}$ is in general a complex number, and we denote them as
$\epsilon_{l}=\alpha_{l}-i\beta_{l}$. This spectrum, originated from the
Lindblad equation, is referred to as the Lindblad spectrum. The full Lindblad
spectrum has been studied for a number of models before Prosen1 ; Prosen2 ;
Universal_spectra ; tenfold ; local_random_Liouvillians ; Wang . Here we would
like to make several useful comments on the Lindblad spectrum. i) $\alpha_{l}$
and $-\alpha_{l}$ always appear in pairs in the spectrum; ii) $\beta_{l}$ is
always non-negative; iii) If $\hat{L}_{\mu}$ are all hermitian, there always
exists a zero-energy eigenstate with $\epsilon_{l}=0$, and this eigenstate is
labelled as $l=0$ and is given by
$|\Psi^{l=0}_{\rho}\rangle=\frac{1}{\sqrt{\mathcal{D}_{\text{H}}}}\sum_{n}|n\rangle\otimes|n\rangle$.

Figure 2: The dynamics of the second Rényi entropy $S^{(2)}$ as a function of
$t\gamma$. $\gamma$ is the dissipation strength. Different curves have
different $\gamma$ in unit of $J$. The inset show the long-time behavior of
$S^{(2)}$ as functions of $tJ$ and $tJ^{2}/\gamma$. The dashed line is a
fitting of initial slop based on Eq. 13. (a) is for the Bose-Hubbard model
with $U=J$ and the number of sites $L=6$, and the number of bosons $N=3$. (b)
is for hard core bosons model with $V=J$, $L=8$ and $N=4$. The initial state
is taken as the ground state of $\hat{H}$.

Rényi Entropy and Lindblad Spectrum. Here we bring out a close relation
between the dynamics of the second Rényi entropy and the Lindblad spectrum.
For any density matrix $\hat{\rho}(t)$, the second Rényi entropy $S^{(2)}(t)$
is given by

$e^{-S^{(2)}}=\text{Tr}(\hat{\rho}^{2})=\sum\limits_{mn}\rho_{mn}(t)\rho_{nm}(t).$
(6)

On the other hand, in the double system, the total amplitude of the wave
function is given by

$|\Psi_{\rho}|^{2}=\sum\limits_{mn}\rho_{mn}(t)\rho^{*}_{mn}(t).$ (7)

Since the density matrix is always Hermitian, it gives
$\rho_{nm}(t)=\rho^{*}_{mn}(t)$, and therefore, we have

$e^{-S^{(2)}}=|\Psi_{\rho}|^{2}.$ (8)

An initial state $\Psi_{\rho}(0)$ in the double space can be expanded as
$\Psi_{\rho}(0)=\sum_{l}c_{l}|\Psi^{l}_{\rho}\rangle$, the subsequent
evolution is given by

$\Psi(t)=e^{-i\hat{H}_{\text{s}}t-\hat{H}_{\text{d}}t}|\Psi_{\rho}(0)\rangle=\sum_{l}c_{l}e^{-i\alpha_{l}t-\beta_{l}t}|\Psi^{n}_{\rho}\rangle$
(9)

and therefore

$e^{-S^{(2)}}=|\Psi_{\rho}|^{2}=\sum\limits_{n}|c_{l}|^{2}e^{-2\beta_{l}t}.$
(10)

Since the evolution in double system is non-unitary and all $\beta_{l}$ are
non-negative, the total amplitude of the wave function always decays in time.
Hence, by this entropy-amplitude relation Eq. 8, the decaying of
$|\Psi_{\rho}|^{2}$ gives rise to the increasing of $S^{(2)}$. Note that for
any initial density matrix with trace unity and for hermitian $\hat{L}_{\mu}$,
$c_{l=0}$ always equals $1/\sqrt{\mathcal{D}_{\text{H}}}$. This mode always
does not decay in time because $\beta_{l=0}=0$. If there is no other
eigenmodes with $\beta_{l}=0$, $l=0$ mode is the only remaining mode at
infinite long time, which gives a maximum second Rényi entropy
$\log\mathcal{D}_{\text{H}}$. Before reaching that limit, the imaginary parts
of the Lindblad spectrum of occupied states determine the time scales of the
Rényi entropy dynamics. Our discussion below will be based on this connection.

Models. Although our discussion below is quite general for quantum many-body
systems, we illustrate the results with two concrete models. The first model
is the Bose-Hubbard model, which reads

$\hat{H}=-J\sum\limits_{\langle
ij\rangle}(\hat{b}^{\dagger}_{i}\hat{b}_{j}+\text{h.c.})+\frac{U}{2}\sum\limits_{i}\hat{n}_{i}(\hat{n}_{i}-1),$
(11)

where $\hat{b}_{i}$ is the boson annihilation operator at site-$i$, and
$\hat{n}_{i}=\hat{b}^{\dagger}_{i}\hat{b}_{i}$ is the boson number operator at
site-$i$. $\langle ij\rangle$ denotes nearest neighbor sites. $J$ and $U$ are
respectively the hopping and the on-site interaction strengths. For the second
model, we consider hard-core bosons, which prevent two bosons to occupy the
same site. In addition, we introduce the nearest-neighbor repulsion, and the
model reads

$\hat{H}=-J\sum\limits_{\langle
ij\rangle}(\hat{b}^{\dagger}_{i}\hat{b}_{j}+\text{h.c.})+V\sum\limits_{\langle
ij\rangle}\hat{n}_{i}\hat{n}_{j}.$ (12)

In one-dimension, these two models are quite different, because the second
model can be mapped to a spinless fermion model with nearest neighbor
repulsion, and can also be mapped to a spin model with nearest neighbor
couplings, but the first model cannot. In both cases, we take all
$\hat{n}_{i}$ as the dissipation operators and we set the dissipation
strengthes uniformly as $\gamma$. In the numerical results shown below, we
have choose $J\sim U$ or $J\sim V$ such that $J$ sets the typical energy scale
of the Hamiltonian part, and therefore, strong and weak dissipations
respectively mean $\gamma/J>1$ or $\gamma/J<1$. Below we will show that both
models exhibit similar features, which supports that our results are quite
universal.

Figure 3: The Lindblad spectrum for strong dissipation case (a1,a2,b1,b2) with
$\gamma=5J$ and for weak dissipation case (c1,c2) with $\gamma=0.2J$. The red
points mark the eigenstates with significant occupation ($|c_{l}|^{2}\geqslant
1/\mathcal{D}_{\text{H}}$) by the initial state. For (a1) and (a2) in the
first raw, the initial state is taken as
$\Psi_{\rho}=|\psi_{\text{g}}\rangle\otimes|\psi_{\text{g}}\rangle$, where
$|\psi_{\text{g}}\rangle$ is the ground state of $\hat{H}$. For (b1) and (b2)
in the second raw, the initial states are taken as the zero-energy eigenstate
of $\hat{H}_{\text{d}}$, that are $|111000\rangle$ for (b1) and
$|11110000\rangle$ for (b2) in Fock bases. The left column (a1,b1,c1) are for
the Bose-Hubbard model with $U=J$ and the number of sites $L=6$, and the
number of bosons $N=3$. The right column (a2,b2,c2) are for hard core bosons
model with $V=J$, $L=8$ and $N=4$.

Dynamics of the Rényi Entropy. We first consider the short-time behavior of
the Rényi entropy dynamics. We apply the short-time expansion to Eq. 9 and
ultilize the relation Eq. 8, and to the leading order of entropy change, we
obtain

$\lim\limits_{t\rightarrow
0}\frac{dS^{(2)}}{dt}=2\frac{\langle\Psi_{\rho}(0)|\hat{H}_{\text{d}}|\Psi_{\rho}(0)\rangle}{\langle\Psi_{\rho}(0)|\Psi_{\rho}(0)\rangle}.$
(13)

The physical meaning of the r.h.s. of Eq. 13 in original system is the
fluctuation of the dissipation operators. For instance, if the initial state
is a pure state and $\hat{\rho}(0)=|\psi(0)\rangle\langle\psi(0)|$, then
$|\Psi_{\rho}(0)\rangle=|\psi(0)\rangle\otimes|\psi(0)\rangle$, and Eq. 13 can
be rewritten as

$\displaystyle\lim\limits_{t\rightarrow 0}\frac{dS^{(2)}}{dt}=$ $\displaystyle
4\sum\limits_{\mu}\gamma_{\mu}\left(\langle\psi(0)|\hat{L}^{\dagger}_{\mu}\hat{L}_{\mu}|\psi(0)\rangle-|\langle\psi(0)|\hat{L}_{\mu}|\psi(0)\rangle|^{2}\right).$
(14)

Suppose all $\gamma_{\mu}$ are taken as the same $\gamma$, this result shows
that the time-dependence of $S^{(2)}$ is governed by a dimensionless time
$\gamma t$. In other word, the larger $\gamma$ is, the faster the Rényi
entropy dynamics increases. This $\gamma t$ scaling is shown in Fig. 2 for two
different models, where one can see that the short-time parts of $S^{(2)}$
curves with different $\gamma$ collapse into a single line when plotted in
term of $\gamma t$. The dashed lines compare the short-time behavior with the
slope given by Eq. 13 and Eq. 14.

In Fig. 2, one also finds that $S^{(2)}$ no longer obeys the $\gamma t$
scaling when $\gamma t>1$. Moreover, in the strong dissipation regime, the
insets plotted in term of $tJ$ show an opposite trend at long-time, that is,
the larger $\gamma$ is, the slower the Rényi entropy increases. In fact, the
long-time behavior of $S^{(2)}$ exhibits a $t/\gamma$ scaling. As shown in the
insets of Fig. 2, when the long-time part of $S^{(2)}$ curves with different
$\gamma$ are ploted in term of $tJ^{2}/\gamma$, they all collapse into a
single curve.

Lindblad Spectrum with Strong Dissipation. This opposite behavior between
short- and long-time can be understood very well in term of the Lindblad
spectrum. As one can see from Fig. 3(a,b), for strong dissipation, the main
feature of the Lindblad spectrum is that it separates into segments along the
imaginary axes of the spectrum, and the separation between segments are
approximately $2\gamma$. For each segment, the width along the imaginary axes
is approximately given by $J^{2}/\gamma$. This feature can be understand by
perturbation treatment of $\hat{H}_{\text{s}}-i\hat{H}_{\text{d}}$. Since the
dissipation strength is stronger than the typical energy scales of the
Hamiltonian, we can treat $\hat{H}_{\text{s}}$ as a perturbation to
$\hat{H}_{\text{d}}$. To the zeroth order of $\hat{H}_{\text{d}}$, the
spectrum is purely imaginary and different segments are separated by
$2\gamma$. More importantly, it worth emphasizing that the eigenstates of
$\hat{H}_{\text{d}}$ are usually highly degenerate, for instance, when
different $\hat{L}_{\mu}$ commute with each other and are related by a
symmetry, such as $\hat{L}_{\mu}$ being $\hat{n}_{i}$ in our examples.
Usually, $\hat{H}_{\text{s}}$ and $\hat{H}_{\text{d}}$ do not commute with
each other, and the perturbation in $\hat{H}_{\text{s}}$ lifts the degeneracy
of the imaginary parts and gives rise to a broadening of the order of
$J^{2}/\gamma$, due to the nature of the second order perturbation.

We call these eigenstates with imaginary energies of the order of a few times
of $\gamma$ as “high imaginary energy states”, and these eigenstates with
imaginary energies of the order of a few times of $J^{2}/\gamma$ as “low-lying
imaginary energy states”. For a generic initial state, both two types of
eigenstates are occupied. Quite generally, the occupations of the “high
imaginary energy states” are significant, for instance, when the initial state
is taken as the eigenstates of $\hat{H}_{\text{s}}$. With the relation between
the Rényi entropy dynamics and the Lindblad spectrum discussed above, it is
clear that the short-time dynamics is dominated by these “high imaginary
energy states” that gives a dynamics scaled by $t\gamma$. Nevertheless, when
$\gamma t>1$, the weights on these “high imaginary energy states” mostly decay
out and the long-time dynamics is therefore dominated by the “low-lying
imaginary energy states” that gives a dynamics scaled by $tJ^{2}/\gamma$.

Figure 4: The dynamics of the second Rényi entropy $S^{(2)}$ as a function of
$tJ^{2}/\gamma$ for specific initial state. $\gamma$ is the dissipation
strength. Different curves have different $\gamma$ in unit of $J$. The inset
show the short-time behavior of $S^{(2)}$ as functions of $tJ$ and $t\gamma$.
(a) is for the Bose-Hubbard model with $U=J$ and the number of sites $L=6$,
and the number of bosons $N=3$. (b) is for hard core bosons model with $V=J$,
$L=8$ and $N=4$. The initial states are taken as the zero-energy eigenstate of
$\hat{H}_{\text{d}}$, that are $|111000\rangle$ for (a) and $|11110000\rangle$
for (b) in Fock bases.

Quantum Zeno Effect Revisited. Here we consider a specific initial state that
satisfies $\hat{H}_{\text{d}}|\Psi(0)\rangle=0$. In other word, such initial
states do not exhibit fluctuation of dissipation operators. Thus, according to
Eq. 13 and Eq. 14, the initial slop of $S^{(2)}$ is zero. Moreover, in the
strong dissipation regime, the populations of the “high imaginary energy
states” are strongly suppressed by the “gap” between different segments and
their contribution becomes negligible, and such initial states mainly populate
the “low-lying imaginary energy states”, as we shown in Fig. 3(b). Therefore,
the entire dynamics of the second Rényi entropy is set by the energy scale
$J^{2}/\gamma$ and it obeys the $t/\gamma$ scaling. This is shown in Fig. 4
for two models. To contrast such specific initial states with generic states
discussed above, we plot in the inset of Fig. 4 the short-time behavior of
$S^{(2)}$ as a function of $t\gamma$ and $tJ$. Unlike the results shown in
Fig. 2, the short-time dynamics with $t\gamma<1$ are quite different, because
it does not exhibit linear behavior and different curves do not collapse into
a single line in term of $t\gamma$.

For these initial states, that the dynamics is slower with stronger
dissipation is reminiscent of the quantum Zeno effect. In fact, the quantum
Zeno effect can indeed be understood in this way. Introducing
$\\{|M\rangle\\},(M=1,\dots,\mathcal{D}_{\text{H}})$ as a set of complete and
orthogonal measurement bases, we define the projection operators as
$\hat{P}_{M}=|M\rangle\langle M|$, and the frequent measurement process can
also be described by the Lindblad equation Eq. 1 with dissipation operator
$\hat{L}_{\mu}$ given by all $\hat{P}_{M}$. With such dissipation operators,
the Lindblad spectrum exhibits a set of “low-lying imaginary energy states”
with energy scale given by $J^{2}/\gamma$. It can be shown that, as long as
the initial state density matrix is diagonal in the measurement bases, the
initial states satisfy $\hat{H}_{\text{d}}|\Psi(0)\rangle=0$.

From Strong to Weak Dissipation. Finally we show that when $\gamma$ decreases
and eventually becomes weaker compared with the typical energy scales in the
Hamiltonian, the segments structure in the Lindblad spectrum disappears, as we
shown in Fig. 3(c). Thus, the entropy dynamics for generic states no longer
display the feature of two time scales. The quantum Zeno effect also
disappears even for the specific initial states, and this is understandable
because in this regime, the typical time interval between two measurements is
already longer than the intrinsic evolution time of the system.

Summary. In this work, we establish a relation between the Rényi entropy
dynamics and the Lindblad spectrum in double space. At the strong dissipation
regime, the Lindblad spectrum exhibits a segment structure, in which we can
introduce the “high imaginary energy eigenstates” and the “low-lying imaginary
energy eigenstates”. For a generic initial state with significantly occupied
“high imaginary energy eigenstates”, the former dominates the short-time
dynamics and the latter dominates the long-time dynamics, which respectively
give rise to $t\gamma$ scaling and $t/\gamma$ scaling. For a specific initial
state with only “low-lying imaginary energy eigenstates” significantly
occupied, the dynamics is dominated by $t/\gamma$ scaling, and we show the
quantum Zeno effect belongs to this class. We illustrate our results with two
concrete models. The second Rényi entropy can now been measured in ultracold
atomic gases in optical lattices, and in fact, it has been measured in the
Bose-Hubbard model with or without disorder Greiner1 ; Greiner2 ; Greiner3 .
The dissipation operators and their strenghes can also now be controlled in
ultracold atomic gases BHMExp , our predictions can therefore be verified
directly in the experimental setup.

Acknowledgment. We thank Lei Pan, Tian-Shu Deng, Tian-Gang Zhou and Pengfei
Zhang for helpful discussions. This work is supported by Beijing Outstanding
Young Scientist Program, NSFC Grant No. 11734010, MOST under Grant No.
2016YFA0301600.

Note Added. When finishing this work, we become aware of a work in which
similar behaviors of the Lindblad spectrum in strong dissipation regime are
also discussed full_spectrum .

## References

  * (1) H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, Oxford, 2007).
  * (2) L. Pan, X. Chen, Y. Chen, and H. Zhai, Nat. Phys. 16, 767(2020).
  * (3) B. Misra and E. C. G. Sudarshan, Journal of Mathematical Physics 18, 756 (1977).
  * (4) W. M. Itano, D. J. Heinzen, J. J. Bollinger, and D. J. Wineland, Phys. Rev. A 41, 2295 (1990).
  * (5) A. G. Kofman and G. Kurizki, Nature 405, 546 (2000).
  * (6) K. Koshino and A. Shimizu, Phys. Rep. 412, 191 (2005).
  * (7) L. Aolita, F. de Melo and L. Davidovich, Rep. Prog. Phys. 78, 042001 (2015)
  * (8) P. Lorenzo, S. Christoph, and X.-L., Qi, J. High Energ. Phys. 2020, 63 (2020).
  * (9) Y. Chen, X.-L. Qi and P. Zhang, J. High Energ. Phys. 2020, 121 (2020).
  * (10) P. Dadras, A. Kitaev, arXiv: 2011.09622
  * (11) Y. Chen, arXiv: 2012.00223
  * (12) K. Su, P. Zhang and H. Zhai, arXiv: 2101.*****
  * (13) J. E. Tyson, J. Phys. A: Math. Gen. 36, 10101 (2003).
  * (14) M. Zwolak and G. Vidal, Phys. Rev. Lett. 93, 207205 (2004).
  * (15) T. Prosen, Phys. Rev. Lett. 109, 090404 (2012)
  * (16) M. V. Medvedyeva, F. H. L. Essler, T. Prosen, Phys. Rev. Lett. 117, 137202 (2016)
  * (17) S. Denisov, T. Laptyeva, W. Tarnowski, D. Chru?ci?ski, and K. Zyczkowski, Phys. Rev. Lett. 123, 140403 (2019)
  * (18) S. Lieu, M. McGinley, and N. R. Cooper, Phys. Rev. Lett. 124, 040401 (2020).
  * (19) K. Wang, F. Piazza, and D. J. Luitz, Phys. Rev. Lett. 124, 100604 (2020).
  * (20) D. Yuan, H. Wang, Z. Wang, D. L. Deng, arXiv: 2009.00019
  * (21) R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Rispoli, M. Greiner, Nature 528, 77 (2015).
  * (22) A. M. Kaufman, M. Eric Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, M. Greiner, Science 353, 794 (2016)
  * (23) A. Lukin, M. Rispoli, R. Schittko, M. Eric Tai, A. M. Kaufman, S. Choi, V. Khemani, J. Léonard, M. Greiner, Science 364, 256 (2019)
  * (24) R. Bouganne, M. B. Aguilera, A. Ghermaoui, J. Beugnon, F. Gerbier, Nat. Phys. 16, 2125 (2020).
  * (25) V. Popkov and C. Presilla, arXiv:2101.05708