# Relativistic calculation of non-dipole effects in high harmonic generation

I. A. Ivanov1<EMAIL_ADDRESS>Kyung Taec Kim1,2<EMAIL_ADDRESS>1
Center for Relativistic Laser Science, Institute for Basic Science (IBS),
Gwangju 61005, Republic of Korea  2Department of Physics and Photon Science,
GIST, Gwangju 61005, Korea

###### Abstract

We present results of relativistic calculations of even order harmonic
generation from various atomic targets. The even order harmonics appear due to
the relativistic non-dipole effects. We take these relativistic effects into
account by using an approach based on the solution of the time-dependent Dirac
equation. The spectra of the non-dipole even harmonics look qualitatively
similar to the spectra of the dipole harmonics obeying the same classical
cutoff rule. The temporal dynamics of the formation of the non-dipole
harmonics is, however, distinctly different from the process of dipole
harmonics formation. Even order harmonics emission is strongly suppressed at
the beginning of the laser pulse, and the emission times of the non-dipole
harmonics are shifted with respect to the bursts of the dipole emission. These
features are partly explained by a simple modification of the classical three-
step model which takes into account selection rules governing the emission of
harmonic photons.

###### pacs:

32.80.Rm 32.80.Fb 42.50.Hz

##  I Introduction

One can expect relativistic effects to play important role in the dynamics of
the processes of atomic or molecular interactions with strong laser pulses for
laser intensities over $10^{18}$ W/cm2 Reiss (1998), when, with increasing
ponderomotive energy, electron velocity can approach speed of light in the
vacuum. It has been realized since the pioneering paper by Reiss Reiss (1990),
however, that relativistic effects may reveal themselves even for moderately
intense ($10^{13}$-$10^{14}$ W/cm2) low frequency infrared (IR) laser fields.
For instance, even for the IR laser fields of intensity of the order of
$10^{13}$ W/cm2, the relativistic effects are visible in the photo-electron
spectra Ludwig et al. (2014a); Chelkowski et al. (2014, 2015); I.A.Ivanov et
al. (2016); Popov et al. (2006); Klaiber and Hatsagortsyan (2014); Yakaboylu
et al. (2015) in the tunneling regime of ionization, characterized by the
values $\gamma\lesssim 1$, where $\gamma=\omega\sqrt{2|I_{p}}/E_{0}$ is the
Keldysh parameter Keldysh (1965), and $\omega$, $E_{0}$ and $I_{p}$ are the
field frequency, field strength and ionization potential of the target system
expressed in atomic units. These relativistic non-dipole effects are due to
the influence of the magnetic field component of the laser pulse which induces
a non-negligible momentum transfer to the photoelectrons Smeenk et al.
(2011a); Ludwig et al. (2014b). Alternatively, if we prefer the photon picture
of light, one might say that an IR photon carries small momentum, but a large
number of the photons participating in the process of the tunneling ionization
Krausz and Ivanov (2009) deliver non-negligible momentum to the ionized
electron Chelkowski and Bandrauk (2018); Chelkowski et al. (2017).

The momentum delivered by the photons to the photo-electron was measured
experimentally under the typical parameters of the tunneling ionization regime
Smeenk et al. (2011b). This momentum manifests itself, on average, as a shift
of the photo-electron momentum distributions (PMD) in the pulse propagation
direction. More detailed picture, which emerges as a result of the complex
interplay of the magnetic and Coulomb forces, includes the so-called direct
electrons which never recollide with the parent ion and are driven in the
direction of the laser photon momentum, and the slow electrons which
experience recollisions and may acquire momentum opposite to the photon
momentum Chelkowski et al. (2015).

Theoretical study of these effects clearly necessitates methods which go
beyond the commonly used non-relativistic dipole approximation. A number of
theoretical procedures allowing to consider the relativistic non-dipole
effects have been described in the literature, including the relativistic
strong-field approximation Reiss (2013); Yakaboylu et al. (2013); Klaiber and
Hatsagortsyan (2014); Yakaboylu et al. (2015), time-dependent Schrödinger
equation (TDSE) with non-dipole corrections Chelkowski and Bandrauk (2018);
Chelkowski et al. (2017); Chelkowski and Bandrauk (2017); I.A.Ivanov et al.
(2016), an approach based on the non-dipole strong-field-approximation
Hamiltonian Jensen et al. (2020), and the time-dependent Dirac equation (TDDE)
Telnov and Chu (2020); Selstø et al. (2009); I.A.Ivanov (2015); Ivanov (2017).

The non-dipole effects manifest themselves as well in other processes
occurring when atoms or molecules interact with laser fields. The process
which will interest us in the present work is the process of the High Harmonic
Generation (HHG). The non-dipole effects are known to produce several
modifications in the HHG spectra. It was found Zhu and Wang (2016) that the
non-dipole interactions lead to decrease of harmonic intensity and shift of
odd order harmonics in the spectra. A detailed investigation of the effect of
the pulse magnetic field on harmonic spectra was reported in Potvliege et al.
(2000); Kylstra et al. (2001); Chirilă et al. (2002). It was found Potvliege
et al. (2000) that the non-dipole magnetic field effects result in the
emission of photons polarized along the propagation direction which, for the
laser pulse wavelength of $800$ nm and intensity of the order of $5\times
10^{15}$ W/cm2, is several orders of magnitude weaker than the photon emission
polarized parallel to the driving pulse polarization direction. For stronger
pulses with the intensities of the order of $10^{17}$ W/cm2, the magnetic
field effects start playing crucial role Chirilă et al. (2002). Electron drift
in the laser propagation direction due to the magnetic-field component of the
laser pulse prevents recollisions, and hence, as one could expect on the basis
of the picture provided by the celebrated three-step model of HHG Lewenstein
et al. (1994); Corkum (1993), leads to the decrease of the harmonic emission.

Perhaps one of the most striking manifestations of the non-dipole effects is
appearance of even order harmonics in the HHG spectra de Aldana and Roso
(2002); Bandrauk and Lu (2006); Mu-Xue et al. (2020), presenting an example of
a relatively small perturbation producing not only relatively minor
quantitative modifications of the spectra, but introducing a qualitative
change: harmonics with frequencies forbidden in the dipole approximation. The
appearance of the even order harmonics can be understood as a result of the
break-up of the well-known symmetry which the electron trajectories
responsible for the emission of the harmonic photons exhibit in the dipole
approximation Lewenstein et al. (1994). Magnetic field effects break this
symmetry, and thus make possible generation of even order harmonics. These
harmonics were studied theoretically in Bandrauk and Lu (2006), using
perturbative treatment of the non-dipole effects.

In the present paper we report a systematic theoretical study of the non-
dipole effects, in particular generation of the even order harmonics, from
various atomic targets. We use the TDDE as our main calculational tool, basing
on the previously developed procedure for the numerical solution of the time-
dependent Dirac equation I.A.Ivanov (2015); Ivanov (2017). The approach based
on the TDDE provides a complete non-perturbative description of the non-
dipole, as well as other relativistic effects.

Atomic units with $\hbar=1$, $e=1$, $m=1$, and $c\approx 137.036$ (here $e$
and $m$ are charge and mass of the electron, $c$\- speed of light) are used
throughout the paper.

##  II Theory

###  II.1 Numerical solution to the time-dependent Dirac equation

We solve the TDDE:

$i\frac{\partial\Psi(r,t)}{\partial t}=\hat{H}\Psi(r,t)$ (1)

following the procedure we described in I.A.Ivanov (2015); Ivanov (2017),
which we briefly recapitulate below for the readers convenience. In Eq. (1)
$\Psi(r,t)$ is a four-component bispinor and the Hamiltonian operator has a
form:

$\hat{H}=\hat{H}_{\rm atom}+\hat{H}_{\rm int}\ ,$ (2)

with:

$\hat{H}_{\rm atom}=c{\bm{\alpha}}\cdot{\hat{\bm{p}}}+c^{2}(\beta-I)+I\ V(r)\
,$ (3)

and

$\hat{H}_{\rm int}=c{\bm{\alpha}}\cdot{\hat{\bm{A}}}\ ,$ (4)

In Eq. (3):

$\displaystyle{\bm{\alpha}}=\left(\begin{array}[]{cc}{\bm{0}}&{\bm{\sigma}}\\\
{\bm{\sigma}}&{\bm{0}}\\\ \end{array}\right)$,
$\displaystyle\beta=\left(\begin{array}[]{cc}{\bm{I}}&{\bm{0}}\\\
{\bm{0}}&-{\bm{I}}\\\ \end{array}\right)$, $\displaystyle
I=\left(\begin{array}[]{cc}{\bm{I}}&{\bm{0}}\\\ {\bm{0}}&{\bm{I}}\\\
\end{array}\right)$, ${\bm{\sigma}}$ are Pauli matrices, ${\bm{0}}$ and
${\bm{I}}$ are $2\times 2$ null and identity matrices, $V(r)$ is the atomic
potential and $c=137.036$\- the speed of light. We subtracted from the field-
free atomic Hamiltonian (3) the constant term $Ic^{2}$ corresponding to the
rest mass energy of the electron.

We use a laser pulse linearly polarized in $z-$ and propagating in $x-$
directions. The vector potential of the pulse is defined in terms of the pulse
electric field:

$\bm{A}(x,t)=-\hat{{\bm{e}}_{z}}\int\limits_{0}^{u}E(\tau)\ d\tau\ ,$ (5)

where $u=t-x/c$. At any given point in space the pulse has a finite duration
$T_{1}$ so that $E(\tau)$ in Eq. (5) is non-zero only for $0<\tau<T_{1}$. As
targets, we will consider below a model atom with a short range (SR) Yukawa-
type potential $V(r)=-1.903e^{-r}/r$, hydrogen atom, and helium atom described
by means of an effective potential Sarsa et al. (2003). The target atom is
initially in the ground $s-$ state $|\phi_{0}\rangle$ with the ionization
potential (IP) of 0.5 a.u. for the hydrogen and Yukawa atoms and IP of 0.902
a.u. for the He atom.

The solution to Eq. (1) is expanded as a series in the basis bispinors:

$\Psi(\mathbf{r},t)=\sum\limits_{j\atop l=j\pm
1/2}\sum\limits_{M=-j}^{j}\Psi_{jlM}({\bm{r}},t),$ (6)

where:

$\Psi_{jlM}({\bm{r}},t)=\left(\begin{array}[]{c}g_{jlM}(r,t)\Omega_{jlM}({\bm{n}})\\\
f_{jlM}(r,t)\Omega_{jl^{\prime}M}({\bm{n}})\\\ \end{array}\right),$ (7)

and the two-component spherical spinors are defined as
$\displaystyle\Omega_{jlM}({\bm{n}})=\left(\begin{array}[]{c}C^{jM}_{l\
M-{1\over 2}{1\over 2}{1\over 2}}Y_{l,M-{1\over 2}}({\bm{n}})\\\ C^{jM}_{l\
M+{1\over 2}{1\over 2}-{1\over 2}}Y_{l,M+{1\over
2}}({\bm{n}})\end{array}\right)$, (here $C^{jM}_{lm{1\over 2}\mu}$ are the
Clebsch-Gordan coefficients, $Y_{lm}({\bm{n}})$\- spherical harmonics, and
${\bm{n}}={\bm{r}}/r$). Parameters $l$ and $l^{\prime}$ in Eq. (6) must
satisfy the relation $l+l^{\prime}=2j$.

To take into account the non-dipole effects due to the spatial dependence of
the laser fields, vector potential (5) is expanded in a series of spherical
harmonics at every time-step of the integration procedure. Substituting
expansion (6) and expansion for the vector potential in the TDDE (1), and
using well-known properties of spherical spinors Akhiezer and Berestetskii
(1965); Lifshitz and Berestetskii (1982), one obtains a system of coupled
differential equations for the radial functions $g_{jlM}(r,t)$ and
$f_{jlM}(r,t)$ in Eq. (7). This system has been solved using a relativistic
generalization of the well-known matrix iteration method (MIM) Nurhuda and
Faisal (1999), which we described in detail in I.A.Ivanov (2015).

Appropriate choice of the propagation technique is essential, as the Dirac
equation, as it is well-known, possesses some properties which are absent in
the case of the non-relativistic wave-equation. These properties are due to
the presence of the continuum of the negative energy states in the Dirac
Hamiltonian which makes the Dirac Hamiltonian unbounded from below. One
problem which this fact entails is the well-known problem of the collapse to
the negative energies continuum Hill and Krauthauser (1994), which may
manifest itself when basis set methods are used to construct approximations to
the bound states of the Dirac Hamiltonian Hill and Krauthauser (1994). We
avoid this problem, since we do not rely on the basis set methods. Initial
state of the system is prepared in our calculation by solving numerically the
eigenvalue equation for the field-free Dirac Hamiltonian employing shooting
method. A related problem is the so-called Zitterbewegung problem Dirac
(1964). Presence of a superposition of the states with positive and negative
energies implies that a solution to the TDDE should exhibit very fast
oscillations with characteristic frequencies of the order of $c^{2}$. Such
oscillations are indeed present and we can reproduce them in the framework of
our numerical procedure by using sufficiently small integration time-step
$\Delta$ I.A.Ivanov (2015). We had to use the time-step $\Delta$ of the order
of $10^{-6}$ a.u. in I.A.Ivanov (2015) to reproduce these oscillations. Use of
such small values for $\Delta$, if it were imperative, would make any
practical calculations impossible, of course. Fortunately, one can bypass this
problem by using an appropriate time-propagation technique. We discussed this
issue in greater detail in Ivanov and Kim (2015); I.A.Ivanov (2015). For
readers convenience we present a core of the argument below. From the purely
numerical point of view, presence of the fast oscillating terms in a system of
the ordinary differential equations (ODE) gives us an example of a numerically
stiff system of ODE, i.e. a system in which vastly different time-scales are
present. To solve such a system of ODE we must use a stable integration method
Shampine (1994), which ensures that while the numerical solution does not
reproduce very fast oscillations, it describes accurately the overall behavior
of the true solution. The integration procedure that we use provides such a
stability. We can illustrate this point using a simple example of a stiff
system of two ODE:

$i\dot{\bm{y}}={\bm{A}}\cdot{\bm{y}},$ (8)

with Hermitian matrix ${\bm{A}}=diag(\lambda_{1}(t),\lambda_{2}(t))$. To mimic
the problem at hand let us assume that $\lambda_{1}$ is of order of $1$, while
$\lambda_{2}$ has large negative value on the interval of time that we
consider. Short-time propagator in the MIM method is a unitary Crank-Nicholson
(CN) propagator Goldberg et al. (1967), which relates solution vectors
${\bm{y}}_{n+1}={\bm{y}}(t_{n+1})$ and ${\bm{y}}_{n}={\bm{y}}(t_{n})$ at times
$t_{n}$ and $t_{n+1}=t_{n}+\Delta$ as follows:

${{\bm{y}}}_{n+1}={1-{i\Delta\over 2}{\bm{A}}(t_{n+1/2})\over 1+{i\Delta\over
2}{\bm{A}}(t_{n+1/2})}{{\bm{y}}}_{n}\ ,$ (9)

where $t_{n+1/2}=t_{n}+\Delta/2$. One can see from Eq. (9) that if at the
$n-$th step of the propagation the second component of the vector ${\bm{y}}$
acquires a numerical error $\delta y^{(2)}_{n}$, the unitarity of the CN
propagation matrix in Eq. (9) makes this error remain bounded for $m>n$.

Spatial variables in the coupled differential equations for the radial
functions $g_{jlM}(r,t)$ and $f_{jlM}(r,t)$ were discretized on a grid with
the step size $\delta r=0.05$ a.u., the radial variable was restricted to an
interval $(0,R_{\rm max})$, with $R_{\rm max}=400$ a.u., and angular momenta
$j$ up to 70 were included in the expansion (6) in the calculations below. The
propagation time-step $\Delta$ was $0.05$ a.u. Before proceeding to the
description of the results of this calculation, it is instructive, however, to
discuss an alternative treatment of the non-dipole effects based on the
leading order perturbation theory (LOPT) expansion, as it provides a more
transparent physical picture of the non-dipole effects than the complete Dirac
equation. LOPT calculation described below was also used as an accuracy test
for our solution to the TDDE.

###  II.2 LOPT treatment of the non-dipole effects.

We are interested in a LOPT solution to the TDDE considering the non-dipole
effects as relativistic corrections.

The leading order relativistic corrections describing the non-dipole effects
in atom-field interaction can be obtained by expanding the minimal coupling
atom-field interaction Hamiltonian Lambropoulos and Petrosyan (2007); Sobelman
(1972); Chelkowski et al. (2015) in the velocity gauge:

$\hat{H}^{\rm min}_{\rm
int}(t)=\hat{\bm{p}}\cdot{\bm{A}}({\bm{r}},t)+{\hat{\bm{A}}^{2}({\bm{r}},t)\over
2}$ (10)

in powers of $c^{-1}$ I.A.Ivanov et al. (2016):

$\hat{H}_{\rm min}(t)=\hat{p}_{z}A(t)+{\hat{v}_{z}xE(t)\over c}+{A^{2}(t)\over
2}+O(c^{-2})\ ,$ (11)

where $\displaystyle E(t)=-{\partial A(t)\over\partial t}$ is the electric
field of the pulse, and the velocity operator
$\displaystyle\hat{\bm{v}}=\hat{\bm{p}}+{\bm{A}}(t)$ has been introduced. The
last term on the r.h.s. of Eq. (11) is a function of time only and can be
removed by a unitary transformation of the wave-function.

Including spin effects in the interaction Hamiltonian is not necessary, if we
are interested in the effects of the leading order in powers of $c^{-1}$
Chelkowski et al. (2015); Zhu and Wang (2016). The fact that the spin degrees
of freedom can be neglected in the leading order of the $c^{-1}$ expansion,
can be understood using the semi-classical picture of the spin effects, in
which additional force due to the presence of the spin degrees of freedom,
acting on the electron, is ${\bm{F}}=-\nabla U_{m}$, where
$U_{m}=-{\bm{\mu}}\cdot{\bm{H}}$, energy of the spin-magnetic field
interaction. Here ${\bm{H}}$ is the magnetic field and ${\bm{\mu}}$ is
electron’s magnetic moment related to the expectation value of electron’s spin
${\bm{\mu}}=-2{\bm{S}}/c$. Spatial gradient of ${\bm{H}}$ introduces an
additional factor of $c^{-1}$, making contribution of the force $\bm{F}$ an
effect of higher order in $c^{-1}$. As for the relativistic corrections to the
field-free atomic Hamiltonian, the so-called Breit-Pauli Hamiltonian Sobelman
(1972), it adds terms of the order of $c^{-2}$ to the non-relativistic atomic
Hamiltonian. We do not have, therefore, to include these corrections in the
LOPT treatment. To the leading order in powers of the $c^{-1}$-expansion, the
dynamics of the system can thus be described by the time-dependent Schrödinger
equation (TDSE):

$i{\partial\Psi({\bm{r}},t)\over\partial t}=\left(\hat{H}_{\rm
atom}+\hat{H}_{\rm d}(t)+\hat{H}_{\rm nd}(t)\right)\Psi({\bm{r}},t)\ ,$ (12)

where

$\hat{H}_{\rm atom}={\hat{\bm{p}}^{2}\over 2}+V(r)$ (13)

is atomic field-free Hamiltonian,

$\hat{H}_{\rm d}(t)=\hat{p}_{z}A(t)$ (14)

is the dipole part of the atom-field interaction and

$\hat{H}_{\rm nd}(t)={\hat{v}_{z}xE(t)\over c}$ (15)

is the non-dipole part of the atom-field interaction containing the effects of
the order of $c^{-1}$.

It is easy to check that the LOPT solution to the equation (12), with the non-
dipole term (15) considered as a perturbation, can be written as:

$\Psi^{\rm LOPT}({\bm{r}},t)=\Psi_{\rm d}({\bm{r}},t)+\Psi^{(1)}_{\rm
nd}({\bm{r}},t)\ ,$ (16)

where the LOPT non-dipole correction is given by the expression:

$\Psi^{(1)}_{\rm nd}({\bm{r}},t)=-i\int\limits_{0}^{t}\hat{U}_{\rm
d}(t,\tau)\hat{H}_{\rm nd}(\tau)\Psi_{\rm d}({\bm{r}},\tau)\ d\tau\ .$ (17)

As can be seen from Eq. (15) for the operator $\hat{H}_{\rm nd}$ this
correction is of the order of $c^{-1}$. In Eq. (16) and Eq. (17) $\Psi_{\rm
d}({\bm{r}},t)$ is the zero-order solution to the non-relativistic TDSE taking
into account only the dipole part of the atom-field interaction, $\hat{U}_{\rm
d}(t,\tau)$ is the evolution operator describing evolution of the system
driven by the non-relativistic dipole Hamiltonian. $\hat{U}_{\rm d}(t,\tau)$
satisfies the operator equation:

$i{\partial\hat{U}_{\rm d}(t,\tau)\over\partial t}=\left(\hat{H}_{\rm
atom}+\hat{H}_{\rm d}(t)\right)\hat{U}_{\rm d}(t,\tau)\ ,$ (18)

and the initial condition $\hat{U}_{\rm d}(\tau,\tau)=\hat{I}$. In practice,
we need not solve the operator equation (18). All we have to do to compute the
expression under the integral on the r.h.s of Eq. (16) for given $\tau$ and
$t$, is to propagate first the initial state wave-function on the interval
$(0,\tau)$ using the non-relativistic TDSE with the Hamiltonian (14),
obtaining thus a state vector $\Psi_{\rm d}(\tau)$. We act than on this vector
with the operator $\hat{H}_{\rm nd}(\tau)$ and propagate it further in time
till the moment $t$. The non-relativistic TDSE was solved using the well-
tested numerical procedure described in Ivanov (2014).

###  II.3 Calculation of electron velocity and HHG spectra

Once the solution to the TDDE (1) is obtained, expectations value of the
electron velocity can be obtained as Avetissian et al. (2011):

${\bm{v}}(t)=c\langle\Psi(t)|{\bm{\alpha}}|\Psi(t)\rangle\ .$ (19)

Harmonic spectra can then be calculated using the usual semi-classical
approach, in which the spectral intensity of the harmonic emission can be
expressed in terms of the Fourier transform of electron’s velocity:

$S_{a}(\Omega)\propto\left|\int\limits_{0}^{T_{1}}v_{a}(t)W(t)e^{i\Omega t}\
dt\right|^{2}\ .$ (20)

where $v_{a}(t)$ is either $x-$ or $z-$ component of the electron velocity for
the non-dipole and dipole harmonic intensities $S_{x}(\Omega)$ and
$S_{z}(\Omega)$, respectively. In the velocity form for the harmonics
intensity which we use here, we do not need to introduce additional powers of
harmonic frequency, which would be present had we used length or acceleration
forms Baggesen and Madsen (2011). The factor $W(t)$ in Eq. (20) is the window
function Reiff et al. (2020), for which we employ the Hann form:
$\displaystyle W(t)=\sin^{2}{\left(\pi t\over T_{1}\right)}$.

The most noticeable effects which the relativistic non-dipole corrections
produce are appearance of harmonic photons polarized in the laser propagation
direction Potvliege et al. (2000); Kylstra et al. (2001); Chirilă et al.
(2002) and appearance of even order harmonics in the HHG spectra Mishra et al.
(2012); Mu-Xue et al. (2020). The LOPT picture allows to explain these
features transparently. Substituting the expression Eq. (16) for the LOPT
wave-function into the matrix element:

$\langle\Psi^{\rm LOPT}(t)|\hat{\bm{v}}|\Psi^{\rm
LOPT}(t)\rangle\approx\hat{\bm{x}}v_{x}(t)+\hat{\bm{y}}v_{y}(t)+\hat{\bm{z}}v_{z}(t)\
,$ (21)

defining the leading order contributions to the expectation value of electron
velocity, one obtains:

$v_{z}(t)=\langle\Psi_{\rm d}(t)|\hat{v}_{z}|\Psi_{\rm d}(t)\rangle\ .$ (22)

For the geometry we use, the evolution operator $\hat{U}_{\rm d}(t,\tau)$
commutes with $\hat{l}_{z}$\- the $z-$ component of the angular momentum,
i.e., it is a conserved quantity for the quantum evolution driven by the
dipole Hamiltonian (13) and (14). ${\hat{l}}_{z}$, therefore, has a definite
value $l_{z}=0$ in the state described by the wave-function $\Psi_{\rm d}(t)$,
and the matrix element $\displaystyle\langle\Psi_{\rm
d}(t)|\hat{v}_{x}|\Psi_{\rm d}(t)\rangle$ vanishes because of the well-known
dipole selection rules Sobelman (1972). Leading order contribution to
$v_{x}(t)$, is, therefore, of the order of $c^{-1}$, and is given by the
expression:

$\displaystyle v_{x}(t)$ $\displaystyle=$ $\displaystyle\langle\Psi_{\rm
d}(t)|\hat{v}_{x}|\Psi^{(1)}_{\rm nd}\rangle+\langle\Psi^{(1)}_{\rm
nd}(t)|\hat{v}_{x}|\Psi_{\rm d}\rangle$ $\displaystyle=$ $\displaystyle 2{\rm
Re}\langle\Psi_{\rm d}(t)|\hat{v}_{x}|\Psi^{(1)}_{\rm nd}\rangle$
$\displaystyle=$ $\displaystyle 2{\rm
Im}\left(\int\limits_{0}^{t}\langle\Psi_{\rm d}(t)|\hat{p}_{x}\hat{U}_{\rm
d}(t,\tau)\hat{H}_{\rm nd}(\tau)|\Psi_{\rm d}(\tau)\rangle\ d\tau\right)\ .$

In the last line of Eq. (LABEL:velx) we used expression (17) for
$\Psi^{(1)}_{\rm nd}$. The same dipole selection rules Sobelman (1972) and the
structure of Eq. (17) ensure that the contribution of the order of $c^{-1}$ to
$v_{y}(t)$ is zero. The leading contribution of the non-dipole effects is,
therefore, non-zero only for the $x$-component of the electron velocity.
Orientation of the dipole velocity due to this relativistic contribution
results, thus, in the appearance of the harmonic photons polarized in the
propagation direction in accordance with the observations made in Potvliege et
al. (2000); Kylstra et al. (2001); Chirilă et al. (2002).

As we mentioned above, the appearance of the even order harmonics can be
understood as a result of violation of the symmetry of the electron
trajectories responsible for the emission of harmonic photons in the dipole
approximation Lewenstein et al. (1994). From the LOPT perspective this effect
can be explained as follows. As one can see from Eq. (14) and Eq. (15), the
dipole interaction operator (14) has odd parity, i.e. it couples states of
different parities, while the non-dipole operator (15) has even parity.
Employing a somewhat lousy language, we might say that the presence of these
two atom-field interaction Hamiltonians can be described as the presence of
two kinds of photons: the ”dipole” photons and the ”non-dipole” photons, whose
emission and absorption are governed by the operators (14) and (15),
respectively. Using these notions and the LOPT expression for $v_{x}(t)$ in
Eq. (LABEL:velx), contribution of the non-dipole interaction to the formation
of the $N-$th harmonic can be described as absorption of $N-1$ ”dipole”
photons and one ”non-dipole” photon, with subsequent recombination to atomic
ground state accompanied by emission of a harmonic photon with frequency
$N\omega$. Using the informal terminology which we adopted, one might say that
the emitted harmonic photon is of the ”dipole” nature since spontaneous
emission satisfies the dipole selection rules. Conservation of the total
parity for the combined system of atom and the ”dipole” and the ”non-dipole”
photons implies then that $N$ must necessarily be even.

Besides providing a simple physical picture of the appearance of even
harmonics, the LOPT approach which we described above, can be used as a test
of the accuracy of our solution to the TDDE. To do such a test we performed
calculations of the expectation values of electron velocity using TDDE and
LOPT approaches for the cosine-pulse form shown in Fig. 1, with the vector
potential in Eq. (5) given by the equation:
$\displaystyle{\bm{A}}(x,t)=-{\bm{e}}_{z}{E_{0}\over\omega}\sin^{2}{\left(\pi
u\over T_{1}\right)}\sin{\omega u}$ where $\omega=0.057$ a.u., $E_{0}=0.0534$
a.u., $u=t-x/c$. A comparison of the TDDE results obtained using Eq. (19) and
the LOPT results obtained using Eq. (LABEL:velx) for the $x-$ component of
electron velocity is shown in Fig. 2. The results of the LOPT treatment prove
to be virtually identical to the results of the TDDE calculation which is not
surprising given that the relativistic corrections could be expected to be
small for the field parameters we consider.

##  III Results

We report below results which we obtained from our TDDE calculations for
dipole $S_{z}(\Omega)$ and non-dipole $S_{x}(\Omega)$ harmonic intensities for
different targets. HHG spectra were obtained by computing electron velocity as
prescribed by Eq. (19) and using Eq. (20) to compute harmonic intensities.
Calculations were performed using the sine waveform shown in Fig. 1 with the
electric field given by the equation: $\displaystyle
E(u)=E_{0}\sin^{2}{\left(\pi u\over T_{1}\right)}\sin{\omega u}$. We report
below results for the base frequencies $\omega=0.114$ a.u. (wavelength of
$400$ nm) and $\omega=0.057$ a.u. (wavelength of $800$ nm).

In Fig. 3 we show HHG spectra that we obtained for the driving pulse
wavelength $\lambda=400$ nm and different field strengths for various targets.
Fig. 3 shows both dipole $S_{z}(\Omega)$ and non-dipole $S_{x}(\Omega)$
harmonic intensities. The vertical lines in the Figures show positions of the
classical cutoffs given by the well-known $3.17U_{p}+I_{p}$ (here
$U_{p}=E_{0}^{2}/4\omega^{2}$ and $I_{p}$ are ponderomotive and ionization
energies respectively) rule of the three-step model Lewenstein et al. (1994);
Corkum (1993). In Fig. 4 we zoom on the parts of the harmonic spectra more
closely to demonstrate the presence of odd and even harmonics in the dipole
and non-dipole spectra respectively.

Quite expectedly, behavior of the dipole intensity $S_{z}(\Omega)$ shown in
Fig. 3 agrees very well with the three-step model predictions, exhibiting a
sharp drop in magnitude after reaching the classical cutoff. The non-dipole
$S_{x}(\Omega)$ spectra mimic this behavior very closely. This may be not
surprising if we make use again of the LOPT picture of formation of the non-
dipole harmonics we presented above, which relied on the notions of ’dipole’
and ’non-dipole’ photons with operators describing their interactions with an
atom given by Eq. (14) and Eq. (15), respectively. We remind, that in the
framework of this picture the $N-$th non-dipole harmonic is produced as a
result of the absorption of $N-1$ ”dipole” photons and one ”non-dipole”
photon. As far as the harmonic spectra are concerned, the mechanism
responsible for the formation of the non-dipole harmonic emission differs thus
from the mechanism of the emission of the dipole harmonics only in the
replacement of one ’dipole’ photon with a ’non-dipole’ one. This replacement
leads to the replacement of the odd order harmonics in the spectra by the even
order ones and results in an overall drop in magnitude in the harmonic spectra
due to the presence of the additional factor of $c^{-1}$ in the non-dipole
interaction operator (15).

The energy and parity conservation considerations which lead us to the general
conclusions about the character of the non-dipole spectra do not tell us
anything about temporal dynamics of the formation of the non-dipole harmonics.
We can have a glimpse of this temporal dynamics by analyzing Gabor transforms
Gabor (1946) of dipole and non-dipole velocities:

$T_{a}(\Omega,t)=\int\limits_{0}^{T_{1}}v_{a}(\tau)\Phi^{*}(t,\tau,\Omega)d\tau\
,$ (24)

where
$\displaystyle\Phi(t,\tau,\Omega)=\exp{\left\\{i\Omega\tau-(t-\tau)^{2}/2(x_{0}T)^{2}\right\\}}$,
parameter $x_{0}$ determines resolution in the temporal domain, and $T$ is an
optical cycle of the laser field. Gabor transform, as well as closely related
wavelet transform, allows us to take a look simultaneously at both time and
frequency domains, and allows to determine, in particular, when different
harmonics are emitted Antoine et al. (1995); Tang et al. (2000); Wang et al.
(2005). We used $x_{0}=0.1$ in the calculations below. This value of $x_{0}$
gives us rather poor resolution in the frequency domain, but high resolution
in the time domain, which is of interest to us presently.

The absolute values $|T_{a}(\Omega,t)|$ for both dipole and non-dipole
velocities are shown in Fig. 5 and Fig. 6 for the SR Yukawa and hydrogen
atoms. One can see that, dynamically, formation of dipole and non-dipole
harmonics proceeds quite differently. For both Yukawa and hydrogen atoms
systems emission of the non-dipole harmonics is strongly suppressed at the
early stages of pulse development, and emission times for the non-dipole
harmonics are shifted with respect to the dipole radiation bursts. Such
behavior could be anticipated by looking at Fig. 2 which shows that $x-$
component of the velocity starts actually respond to the field only for times
approaching the midpoint of the pulse. The reason for this could be traced
back to the character of the fully quantum expression for the velocity
component $v_{x}$ in the second LOPT equation (LABEL:velx), with time
integration on the right-hand side of this equation smoothing out high
frequency oscillations. To elucidate this issue further we performed a simple
classical calculation of the emitted photon energy as a function of the
recombination time using the physical picture provided by the three-step
model. We assume that electron is ionized at the moment of time $t_{ion}$ and
returns to the parent ion at the moment of time $t_{ret}$, emitting a harmonic
photon with energy $E_{ret}+I_{p}$. As is usually assumed in the three-step
model calculations, we consider only the effect of the external field (5) on
the electron motion, neglecting completely ionic potential. The only
difference of our calculation and the traditional three-step model analysis of
the harmonic emission, is that we take into account effect of the Lorentz
force due to the magnetic field of the pulse. We simulate electron motion in a
plane (which is the $(x,z)$\- plane for the geometry we employ), solving the
set of the classical Newton equations, which for the fields configuration,
geometry and atomic units system we employ, can be written as:

$\displaystyle\ddot{x}$ $\displaystyle=$ $\displaystyle-{{v}_{z}\over c}E(t)$
$\displaystyle\ddot{z}$ $\displaystyle=$ $\displaystyle-E(t)+{{v}_{x}\over
c}E(t)$ $\displaystyle\ .$ (25)

Following the prescription of the traditional three-step model we solve
equations (25) with zero initial conditions imposed at the ionization time:
$v_{x}(t_{ion})=v_{z}(t_{ion})=0$ and $x(t_{ion})=z(t_{ion})=0$. We assume
that the electron trajectory returns to the origin, if at the moment of time
$t_{ret}$, $z-$coordinate of the electron trajectory changes sign.

Fig. 7(a) shows results of such a simulation, which qualitatively agree with
the dynamics of the dipole harmonics emission shown in Fig. 5 and Fig. 6, with
bursts of harmonics emission occurring every half cycle of the laser pulse. To
be able to apply this classical analysis to the emission of the non-dipole
harmonics we must, however, introduce one essentially quantum ingredient in
the model described by the classical equations(25). Emission of the non-dipole
radiation differs from the emission of the dipole harmonics in one important
aspect. For the geometry we employ, the dipole harmonics photon emission
process satisfies selection rule $\Delta M=0$, where $M$ is the $z-$
projection of the electron angular momentum. On the other hand, emission of
the non-dipole harmonic photon, as can be seen from the LOPT analysis we
presented above, must satisfy selection rule $\Delta M=\pm 1$. This means that
for the ground $s-$ state that we consider, non-dipole radiation can be
emitted only by electrons with non-zero angular momentum. We can incorporate
this fact in our classical model by introducing a filter parameter $f$ in the
simulations, and considering only those returning trajectories for which at
the moment of time $t_{ret}$ squared classical angular momentum value exceeds
the threshold value set by the filter parameter $f$. Results of such
calculations are shown in Fig. 7(b-d) for different values of the filter
parameter $f$. One can see that by increasing the value of the filter
parameter, we make the classical picture in Fig. 7 look more like the Gabor
transform results shown in Fig. 5 and Fig. 6. In particular, Fig. 7(c-d) show
the absence of the non-dipole harmonics emission during the first two cycles
of the laser pulse, the feature which is also demonstrated by the quantum
analysis based on the Gabor transform in Fig. 5 and Fig. 6. Applying non-zero
filter parameter does not change, however, the maximum energy $E_{ret}$ of the
returning electron, which explains why non-dipole harmonic emission spectra
exhibit essentially the same cutoffs as the dipole harmonic emission spectra.
This simple classical picture of the formation of the non-dipole harmonics,
which takes as quantum ingredient only the requirement that the electron
angular momentum on the returning trajectories should exceed certain threshold
value, agrees, thus, qualitatively with the fully quantum picture.

We also performed TDDE calculations for the pulse base frequency
$\omega=0.057$ a.u. (corresponding to the wavelength of $800$ nm). In Fig. 8
and Fig. 9 we show harmonic spectra we obtain from TDDE for the SR Yukawa and
hydrogen atoms. Fig. 10 shows results of the analysis of the temporal dynamics
of the harmonic formation based on the Gabor transform (24). These Figures
show essentially the same picture as the results we presented above for the
driving pulse wavelength of $400$ nm. The spectra of the non-dipole harmonics
follow closely the classical dipole cutoff rule, and differ in this respect
from the dipole emission spectra only in their intensity. Temporal pictures of
the harmonics formation in the dipole and the non-dipole cases are, however,
totally different. The main difference is, just as in the case of the driving
pulse wavelength of $400$ nm, the absence of the harmonic emission at the
early stages of the pulse development, the feature which we explained above
using the results of the classical calculations shown in Fig. 7(c,d).

The factor which is responsible for the difference in intensity between the
dipole and non-dipole harmonics is the additional factor of $c^{-1}$ which, as
one can see from Eq. (LABEL:velx) and Eq. (15), is present in the LOPT formula
for the $x-$component of the velocity. The presence of this factor in $v_{x}$
leads to a dampening factor of $c^{-2}$ in the expression for the non-dipole
harmonics intensity. It is rather difficult to obtain a more detailed insight
about relative magnitude of the dipole and non-dipole harmonic intensities
form the cumbersome LOPT expressions Eq. (22) and Eq. (LABEL:velx). One can,
however, obtain a simple estimate using the reasoning based not on the
Schrödinger picture that we have used so far, but on the equivalent Heisenberg
picture of the quantum mechanics (QM). In the latter, we remind, the operators
evolve in time, while the state vectors do not. We obtain, of course, the same
expectation values for all physical observables in both pictures.

In the Heisenberg picture time-evolution of the operators ${\hat{\bm{r}}}(t)$
and ${\hat{\bm{p}}}(t)$ is described by the equations Landau and Lifshitz
(1977): $\displaystyle i{\dot{\hat{\bm{r}}}}=[\hat{{\bm{r}}},\hat{H}]$,
$\displaystyle i{\dot{\hat{\bm{p}}}}=[{\hat{\bm{p}}},\hat{H}]$, where the
Hamiltonian operator in our problem is $\hat{H}=\hat{H}_{\rm
atom}+\hat{H}_{\rm d}(t)+\hat{H}_{\rm nd}(t)$, with $\hat{H}_{\rm atom}$,
$\hat{H}_{\rm d}(t)$ and $\hat{H}_{\rm nd}(t)$ given by Eq. (13), Eq. (14),
and Eq. (15), respectively. Calculating the commutators, one obtains the
following equations of motion:

$\displaystyle\dot{\hat{x}}$ $\displaystyle=$ $\displaystyle\hat{p}_{x}$
$\displaystyle\dot{\hat{p}}_{x}$ $\displaystyle=$
$\displaystyle-i[{\hat{p}}_{x},\hat{V}]-{\hat{v}_{z}\over c}E(t)$
$\displaystyle\dot{\hat{z}}$ $\displaystyle=$
$\displaystyle\hat{p}_{z}+A(t)+{{\hat{x}}\over c}E(t)$
$\displaystyle\dot{\hat{p}}_{z}$ $\displaystyle=$
$\displaystyle-i[{\hat{p}}_{z},\hat{V}]\ ,$

where $\hat{V}$ is atomic potential operator, $\hat{v}_{z}=\hat{p}_{z}+A(t)$,
$A(t)$ and $E(t)$ are the vector potential and the electric field of the
pulse. Eq. (LABEL:he) is the quantum-mechanical analogue of the classical
equations describing electron motion in the potential $V$ in presence of the
Lorenz force. It contains the same physical information and is, therefore,
equivalent to the LOPT equations Eq. (22) and Eq. (LABEL:velx), but it
provides a more clear physical picture and can be used as a starting point for
making simplifying assumptions.

From the first two equations (LABEL:he) one obtains:

$\ddot{\hat{x}}=-i[{\hat{p}}_{x},\hat{V}]-{\hat{v}_{z}\over c}E(t)\ ,$ (27)

We will make an assumption that one can omit the commutator
$[{\hat{p}}_{x},\hat{V}]$ in Eq. (27). Some justification for this operation
can be provided in the case of the SR Yukawa atom, when potential function
$V({\bm{r}})$ is effectively zero everywhere excepting a small neighborhood of
the atom. We obtain then from Eq. (27) a relation for the expectation values
of the electron acceleration $\displaystyle
a_{x}=\langle\phi_{0}|\ddot{\hat{x}}|\phi_{0}\rangle$ and velocity
$\displaystyle v_{z}=\langle\phi_{0}|{\hat{v}}_{z}|\phi_{0}\rangle$:

$a_{x}=-{v_{z}\over c}E(t)\ ,$ (28)

where $|\phi_{0}\rangle$ is the initial atomic state, which does not evolve in
time in the Heisenberg picture. Assuming further that $E(t)$ is a
monochromatic wave: $E(t)=E_{0}\cos{\omega t}$ and calculating Fourier
transforms of both sides of Eq. (28), we obtain a relation between the Fourier
transforms $\displaystyle\tilde{v}_{x}(\Omega)=\int v_{x}(t)e^{i\Omega t}\ dt$
and $\displaystyle\tilde{v}_{z}(\Omega)=\int v_{z}(t)e^{i\Omega t}\ dt$:

$-i\Omega\tilde{v}_{x}(\Omega)={E_{0}\over
2c}\left(\tilde{v}_{z}(\Omega+\omega)+\tilde{v}_{z}(\Omega-\omega)\right)\ ,$
(29)

from which, using the fact that for any complex numbers $z_{1}$, $z_{2}$:
$\displaystyle|z_{1}+z_{2}|^{2}\leq(|z_{1}|+|z_{2}|)^{2}$, we obtain an
inequality:

$\Omega^{2}S_{x}(\Omega)\leq{E_{0}^{2}\over
4c^{2}}\large(\sqrt{S_{z}(\Omega+\omega)}+\sqrt{S_{z}(\Omega-\omega)}\large)^{2}$
(30)

We see from Eq. (30) that for $\Omega>\omega$ we have:

$R(\Omega)={4c^{2}\omega^{2}\over
E_{0}^{2}}{S_{x}(\Omega)\over\large(\sqrt{S_{z}(\Omega+\omega)}+\sqrt{S_{z}(\Omega-\omega)}\large)^{2}}\leq
1\ .$ (31)

Introducing the magnitude $A_{0}=E_{0}/\omega$ of the pulse vector potential,
we can rewrite inequality (31) as:

${S_{x}(\Omega)\over\large(\sqrt{S_{z}(\Omega+\omega)}+\sqrt{S_{z}(\Omega-\omega)}\large)^{2}}\leq{A_{0}^{2}\over
4c^{2}}$ (32)

The ratio $R(\Omega)$ defined in Eq. (31) is shown in Fig. 11(a) for the SR
Yukawa potential and various pulse parameters. Of course, we cannot expect Eq.
(31) to provide a rigorous upper bound since deriving it we neglected atomic
potential in Eq. (27), which constitutes a rather drastic approximation. As
one can see from Fig. 11(a) inequality (31) can indeed be violated. One can
see, nevertheless, that Eq. (31), and consequently Eq. (32) provide reasonably
accurate estimates of the relative magnitude of the intensities of the dipole
and non-dipole harmonics.

While the non-zero expectation value $v_{x}$ and appearance of the non-dipole
harmonics is an entirely relativistic phenomena, the non-dipole effects also
modify slightly the velocity component $v_{z}$. The magnitude of this effect
is of the order of $c^{-2}$. This can be most easily seen from the Heisenberg
equations of motion (LABEL:he). The equation for $v_{z}(t)$ (the third of the
equations (LABEL:he)) contains the term ${\hat{x}}E(t)/c$ on the right-hand
side. Since the expectation value of $x$ is itself of the order of $c^{-1}$,
the resulting effect on $v_{z}(t)$ is of the order of $c^{-2}$, which will
produce a relativistic correction of the order of $c^{-2}$ for the dipole
harmonic intensity. We may expect, therefore, that the normalized difference:

${\Delta S_{z}(\Omega)\over S_{z}(\Omega)}={S_{z}(\Omega)-S^{\rm
nr}_{z}(\Omega)\over S^{\rm nr}_{z}(\Omega)}\ ,$ (33)

where $S_{z}(\Omega)$ is the dipole harmonics intensity obtained in the
present TDDE calculation and $S^{\rm nr}_{z}(\Omega)$ is the result of the
non-relativistic TDSE calculation, should be of the order of $c^{-2}$. i.e.,
we may expect $\Delta S_{z}(\Omega)/S_{z}(\Omega)\sim 10^{-4}$. That this is
indeed the case can be seen from Fig. 11(b), where we show results of the TDDE
and TDSE calculations performed for the same pulse parameters for the Yukawa
atom.

The analysis based on the Heisenberg equations of motion (LABEL:he) also
allows to give a simple explanation for the behavior of $v_{x}(t)$ shown in
Fig. 2, where the $x-$ component of the electron velocity starts responding to
the field only for the times approaching the midpoint of the pulse.
Integrating Eq. (28) we obtain for the expectation value $\displaystyle
v_{x}=\langle\phi_{0}|{\hat{v}}_{x}|\phi_{0}\rangle$ (assuming that it has
zero value at $t=0$):

$v_{x}(t)=-{1\over c}\int\limits_{0}^{t}v_{z}(\tau)E(\tau)\ d\tau\ .$ (34)

We could have obtained the same equation by integrating the first of the set
of the classical equations (25), which is not surprising given the great
formal similarity between the classical mechanics and the QM in the Heisenberg
picture. We show in Fig. 12(a) the expectation value $v_{z}(t)$ obtained in
the LOPT calculation for the cosine pulse with $E_{0}=0.0534$ a.u. and
$\omega=0.057$ a.u. We show only the LOPT result. Just as in the case of
$v_{x}(t)$, shown in Fig. 2, the TDDE and LOPT results for $v_{z}(t)$ differ
very slightly. In Fig. 12(b) we show the LOPT expectation value $v_{x}(t)$, as
well as the estimate for $v_{x}(t)$ that we obtain if we substitute the LOPT
value for $v_{z}(\tau)$ under the integral sign in Eq. (34). One can see that
the estimate thus obtained reproduces fairly well the general behavior of
$v_{x}(t)$. In particular, it reproduces the feature that we mentioned above:
the $x-$component of the velocity begins deviating from zero appreciably only
for the times approaching the midpoint of the pulse. We remind that effect of
the atomic potential on the motion in the $x-$direction was neglected in the
Heisenberg equation of motion (28) which we used to obtain the estimate (34).
The fact that the estimate (34) reproduces qualitative behavior of the
$x-$component of electron velocity shown in Fig. 2 tells us, therefore, that
this behavior is a result of the interplay of the motion in $x-$ and $z-$
directions which are mutually interconnected due to presence of the Lorentz
force.

##  IV Conclusion

We have presented results of the relativistic calculations of even harmonic
generation from various atomic targets. Our approach was based on the
numerical solution of the TDDE. The HHG spectra of the non-dipole even order
harmonics were found to look qualitatively similar to the spectra of the
dipole harmonics, obeying the same classical cutoff rules. The temporal
formation of the non-dipole harmonics, however, was found to be quite
different. The results of the Gabor transform analysis show that formation of
the non-dipole harmonics is strongly suppressed at the beginning of the laser
pulse, and bursts of the non-dipole radiation are shifted in time with respect
to the bursts of the dipole emission. These features are partly explained by a
simple generalization of the classical three-step model, which takes into
account the selection rules governing emission of harmonic photons. We modeled
the effect of these selection rules by using a filter parameter, which selects
the trajectories with angular momentum exceeding a certain threshold value at
the recollision time.

For the field parameters we considered the relativistic effects are still
relatively weak and could be described perturbatively. LOPT provides, as we
have seen, an adequate description of the non-dipole effects responsible for
the even order harmonics emission. Use of the TDDE, however, is technically
simpler than the calculations based on the LOPT, and opens the perspective of
making an excursion into the truly relativistic domain in the future. We
relied, therefore, on the TDDE-based approach in the present work. The present
approach can also be generalized relatively easily to include some quantum
electrodynamical (QED) effects, such as the vacuum polarization effects, or
the QED strong Coulomb field radiative corrections, which can be taken into
account by using effective potentials such as the Uehling potential Uehling
(1935a, b) or the radiative potential proposed in Flambaum and Ginges (2005).
The procedure we apply to solve the Dirac equation can also be used to study
the process of electron-positron pair production (PP) in strong
electromagnetic fields, which occurs when field strength reaches the
characteristic Schwinger field strength of $1.3\times 10^{16}$ V/cm. The
process of PP in both homogeneous and inhomogeneous electric fields has
received considerable interest in the literature Gies and Klingmüller (2005).
Theoretical treatment of PP in the semiclassical approximation relies on a
solution of the TDDE for a given field configuration Mocken et al. (2010). Our
procedure might prove useful for this purpose, especially in the case of the
spatially inhomogeneous field, which has been found to play an important role
in the PP Gies and Klingmüller (2005); Kohlfürst (2020).

The numerical procedure we employ can be relatively easily generalized for the
case of the many-electron relativistic Hamiltonians used in the quantum
chemistry calculations Nakai (2021); Liu (2020). Use of the representation of
the wave-function analogous to the expansion (6) would be, of course,
impractical for systems with more than one electron if we want to use such
expansions to represent the wave-function in the whole space. One may use,
however, the idea of the $R-$matrix approach, which separates the coordinate
space in the inner region, where a suitable basis set representation can be
used to represent many-electron wave-functions and the outer region, where one
has to concentrate on the description of a single electron motion, for which
the finite difference method might be better suited. Such a strategy has been
implemented with success in the framework of the so-called R-Matrix
incorporating Time method (RMT) Lysaght et al. (2012) which allows to solve
the non-relativistic TDSE for many-electron systems. One can use a similar
approach in the relativistic case, relying on the results of the stationary
quantum chemistry calculations Nakai (2021); Liu (2020) for the description of
the inner region, where many-electron effects are important, and using the
present procedure to solve the TDDE describing electron propagation in the
outer region.

## Acknowledgments

This work was supported by the Institute for Basic Science grant (IBS-R012-D1)
and the National Research Foundation of Korea (NRF), grant funded by the Korea
government (MIST) (No. 2022R1A2C3006025). Computational works for this
research were performed on the IBS Supercomputer Aleph in the IBS Research
Solution Center.

Figure 1: (Color online) Pulse shapes $E(t)$ employed in the calculations. Figure 2: (Color online) Expectation value of the $x-$ component of the electron velocity as a function of time obtained in TDDE and LOPT calculations Cosine pulse with $E_{0}=0.0534$ a.u., $\omega=0.057$ a.u. has been used in the calculation. Figure 3: (Color online) Dipole (dash green) and non-dipole (red solid) harmonic intensities for the pulse wavelength $\lambda=400$ nm for the SR Yukawa, hydrogen and He atoms. Vertical dash lines show cutoff positions.  Figure 4: (Color online) Dipole (dash green) and non-dipole (red solid) harmonic intensities for the pulse wavelength $\lambda=400$ nm for harmonics with orders $n\leq 20$ for the SR Yukawa and hydrogen atoms. Vertical solid and dash lines show positions of odd and even harmonics, respectively.  Figure 5: (Color online) Gabor transform $|T(\Omega,t)|$ for the pulse wavelength $\lambda=400$ nm and different field strengths for the SR Yukawa atom.  Figure 6: (Color online) Gabor transform $|T(\Omega,t)|$ for the pulse wavelength $\lambda=400$ nm and different field strengths for the hydrogen atom.  Figure 7: (Color online) Classical calculations of emitted photon energy as function of return time for the dipole harmonic radiation (a) and non-dipole harmonic radiation (b-d) with different filter parameters.  Figure 8: (Color online) Dipole (dash green) and non-dipole (red solid) harmonic intensities for the pulse wavelength $\lambda=800$ nm. for the SR Yukawa and hydrogen atoms. Vertical dash lines show cutoff positions.  Figure 9: (Color online) Dipole (dash green) and non-dipole (red solid) harmonic intensities for the pulse wavelength $\lambda=800$ nm for harmonics with orders $40\leq n\leq 60$ for the SR Yukawa and hydrogen atoms. Vertical solid and dash lines show positions of odd and even harmonics, respectively.  Figure 10: (Color online) Gabor transform $|T(\Omega,t)|$ for the pulse wavelength $\lambda=800$ nm for the SR Yukawa and hydrogen atoms.  a) |  b)
---|---
|
Figure 11: (Color online) (a) Estimate (31) for the ratio $R(\Omega)$. (b) Normalized difference $\left(S_{z}(\Omega)-S^{\rm nr}_{z}(\Omega)\right)/S^{\rm nr}_{z}(\Omega)$ of the TDDE and TDSE calculations for the dipole harmonic intensities. a) |  b)
---|---
|
Figure 12: (Color online) (a) LOPT expectation value $v_{z}(t)$. (b)
$v_{x}(t)$ obtained in the LOPT calculation and using Eq. (34). Cosine pulse
with $E_{0}=0.0534$ a.u., $\omega=0.057$ a.u. has been used in the
calculations.

## References

  * Reiss (1998) H. R. Reiss, Opt. Express 2, 261 (1998).
  * Reiss (1990) H. R. Reiss, J. Opt. Soc. Am. B 7, 574 (1990).
  * Ludwig et al. (2014a) A. Ludwig, J. Maurer, B. W. Mayer, C. R. Phillips, L. Gallmann, and U. Keller, Phys. Rev. Lett. 113, 243001 (2014a).
  * Chelkowski et al. (2014) S. Chelkowski, A. D. Bandrauk, and P. B. Corkum, Phys. Rev. Lett. 113, 263005 (2014).
  * Chelkowski et al. (2015) S. Chelkowski, A. D. Bandrauk, and P. B. Corkum, Phys. Rev. A 92, 051401(R) (2015).
  * I.A.Ivanov et al. (2016) I.A.Ivanov, J. Dubau, and K. T. Kim, Phys. Rev. A 94, 033405 (2016).
  * Popov et al. (2006) V. S. Popov, B. M. Karnakov, V. D. Mur, and S. G. Pozdnyakov, Sov. Phys. -JETP 102, 760 (2006).
  * Klaiber and Hatsagortsyan (2014) M. Klaiber and K. Z. Hatsagortsyan, Phys. Rev. A 90, 063416 (2014).
  * Yakaboylu et al. (2015) E. Yakaboylu, M. Klaiber, and K. Z. Hatsagortsyan, Phys. Rev. A 91, 063407 (2015).
  * Keldysh (1965) L. V. Keldysh, Sov. Phys. -JETP 20, 1307 (1965).
  * Smeenk et al. (2011a) C. T. L. Smeenk, L. Arissian, B. Zhou, A. Mysyrowicz, D. M. Villeneuve, A. Staudte, and P. B. Corkum, Phys. Rev. Lett. 106, 193002 (2011a).
  * Ludwig et al. (2014b) A. Ludwig, J. Maurer, B. W. Mayer, C. R. Phillips, L. Gallmann, and U. Keller, Phys. Rev. Lett. 113, 243001 (2014b).
  * Krausz and Ivanov (2009) F. Krausz and M. Ivanov, Rev. Mod. Phys. 81, 163 (2009).
  * Chelkowski and Bandrauk (2018) S. Chelkowski and A. D. Bandrauk, Phys. Rev. A 97, 053401 (2018).
  * Chelkowski et al. (2017) S. Chelkowski, A. D. Bandrauk, and P. B. Corkum, Phys. Rev. A 95, 053402 (2017).
  * Smeenk et al. (2011b) C. T. L. Smeenk, L. Arissian, B. Zhou, A. Mysyrowicz, D. M. Villeneuve, A. Staudte, and P. B. Corkum, Phys. Rev. Lett. 106, 193002 (2011b).
  * Reiss (2013) H. R. Reiss, Phys. Rev. A 87, 033421 (2013).
  * Yakaboylu et al. (2013) E. Yakaboylu, M. Klaiber, H. Bauke, K. Z. Hatsagortsyan, and C. H. Keitel, Phys. Rev. A 88, 063421 (2013).
  * Chelkowski and Bandrauk (2017) S. Chelkowski and A. D. Bandrauk, Molecular Physics 115, 1971 (2017).
  * Jensen et al. (2020) S. V. B. Jensen, M. M. Lund, and L. B. Madsen, Phys. Rev. A 101, 043408 (2020).
  * Telnov and Chu (2020) D. A. Telnov and S.-I. Chu, Phys. Rev. A 102, 063109 (2020).
  * Selstø et al. (2009) S. Selstø, E. Lindroth, and J. Bengtsson, Phys. Rev. A 79, 043418 (2009).
  * I.A.Ivanov (2015) I.A.Ivanov, Phys. Rev. A 91, 043410 (2015).
  * Ivanov (2017) I. A. Ivanov, Phys. Rev. A 96, 013419 (2017).
  * Zhu and Wang (2016) X. Zhu and Z. Wang, Optics Communications 365, 125 (2016).
  * Potvliege et al. (2000) R. M. Potvliege, N. J. Kylstra, and C. J. Joachain, J. Phys. B 33, L743 (2000).
  * Kylstra et al. (2001) N. J. Kylstra, R. M. Potvliege, and C. J. Joachain, J. Phys. B 34, L55 (2001).
  * Chirilă et al. (2002) C. C. Chirilă, N. J. Kylstra, R. M. Potvliege, and C. J. Joachain, Phys. Rev. A 66, 063411 (2002).
  * Lewenstein et al. (1994) M. Lewenstein, P. Balcou, M. Y. Ivanov, A. L’Huillier, and P. B. Corkum, Phys. Rev. A 49, 2117 (1994).
  * Corkum (1993) P. B. Corkum, Phys. Rev. Lett. 71, 1994 (1993).
  * de Aldana and Roso (2002) J. R. V. de Aldana and L. Roso, J. Phys. B 35, 1633 (2002).
  * Bandrauk and Lu (2006) A. D. Bandrauk and H. Z. Lu, Phys. Rev. A 73, 013412 (2006).
  * Mu-Xue et al. (2020) W. Mu-Xue, C. Si-Ge, L. Hao, and P. Liang-You, Chinese Physics B 29, 013302 (2020).
  * Sarsa et al. (2003) A. Sarsa, F. J. Gálvez, and E. Buendia, J. Phys. B 36, 4393 (2003).
  * Akhiezer and Berestetskii (1965) A. Akhiezer and V. Berestetskii, _Quantum Electrodynamics_ (John Wiley & Sons, 1965).
  * Lifshitz and Berestetskii (1982) E. M. Lifshitz and V. B. Berestetskii, _Quantum Electrodynamics_ (Pergamon Press, 1982).
  * Nurhuda and Faisal (1999) M. Nurhuda and F. H. M. Faisal, Phys. Rev. A 60, 3125 (1999).
  * Hill and Krauthauser (1994) R. N. Hill and C. Krauthauser, Phys. Rev. Lett. 72, 2151 (1994).
  * Dirac (1964) P. A. M. Dirac, _The Principles of Quantum Mechanics_ (McGrow-Hill, New York, 1964).
  * Ivanov and Kim (2015) I. A. Ivanov and K. T. Kim, Phys. Rev. A 92, 053418 (2015).
  * Shampine (1994) L. Shampine, _Numerical Solution of Ordinary Differential Equations_ (Chapman $\&$ Hall, New York, NY, 1994).
  * Goldberg et al. (1967) A. Goldberg, H. M. Schey, and J. L. Scwartz, Am. J. Phys. 35, 177 (1967).
  * Lambropoulos and Petrosyan (2007) P. Lambropoulos and D. Petrosyan, _Fundamentals of Quantum Optics and Quantum Information_ (Springer-Verlag, Berlin, 2007).
  * Sobelman (1972) I. I. Sobelman, _Introduction to the Theory of Atomic Spectra_ (Pergamon Press, New York, 1972).
  * Ivanov (2014) I. A. Ivanov, Phys. Rev. A 90, 013418 (2014).
  * Avetissian et al. (2011) H. K. Avetissian, A. G. Markossian, and G. F. Mkrtchian, arXiv: Atomic Physics (2011).
  * Baggesen and Madsen (2011) J. C. Baggesen and L. B. Madsen, J. Phys. B 44, 115601 (2011).
  * Reiff et al. (2020) R. Reiff, T. Joyce, A. Jaroń-Becker, and A. Becker, J. Phys. Commun 4, 065011 (2020).
  * Mishra et al. (2012) R. Mishra, D. Kalita, and A. Gupta, Eur. Phys. J. D 66, 169 (2012).
  * Gabor (1946) D. Gabor, J. Inst. Electr. Eng. 93, 429 (1946).
  * Antoine et al. (1995) P. Antoine, B. Piraux, and A. Maquet, Phys. Rev. A 51, R1750 (1995).
  * Tang et al. (2000) Y. Y. Tang, L. H. Yang, J. Liu, and H. Ma, _Wavelet Theory and Its Application to Pattern Recognition_ (World Scientific, 2000).
  * Wang et al. (2005) B. Wang, T. Cheng, X. Li, P. Fu, S. Chen, and J. Liu, Phys. Rev. A 72, 063412 (2005).
  * Landau and Lifshitz (1977) L. D. Landau and E. M. Lifshitz, _Quantum Mechanics_ (Pergamon Press, New York, 1977).
  * Uehling (1935a) E. A. Uehling, Phys. Rev. 48, 55 (1935a).
  * Uehling (1935b) E. A. Uehling, Physical Review 48, 55 (1935b).
  * Flambaum and Ginges (2005) V. V. Flambaum and J. S. M. Ginges, Phys. Rev. A 72, 052115 (2005).
  * Gies and Klingmüller (2005) H. Gies and K. Klingmüller, Phys. Rev. D 72, 065001 (2005).
  * Mocken et al. (2010) G. R. Mocken, M. Ruf, C. Müller, and C. H. Keitel, Phys. Rev. A 81, 022122 (2010).
  * Kohlfürst (2020) C. Kohlfürst, Phys. Rev. D 101, 096003 (2020).
  * Nakai (2021) H. Nakai, Bull. Chem. Soc. Jpn. 94, 1664 (2021).
  * Liu (2020) W. Liu, J. Chem. Phys. 152, 180901 (2020).
  * Lysaght et al. (2012) M. A. Lysaght, L. R. Moore, L. A. A. Nikolopoulos, J. S. Parker, H. W. van der Hart, and K. T. Taylor, Journal of Physics: Conference Series 388, 012027 (2012).


# Probing many-body systems near spectral degeneracies

K. Ziegler  Institut für Physik, Universität Augsburg, D-86135 Augsburg,
Germany

###### Abstract

The diagonal elements of the time correlation matrix are used to probe closed
quantum systems that are measured at random times. This enables us to extract
two distinct parts of the quantum evolution, a recurrent part and an
exponentially decaying part. This separation is strongly affected when
spectral degeneracies occur, for instance, in the presence of spontaneous
symmetry breaking. Moreover, the slowest decay rate is determined by the
smallest energy level spacing, and this decay rate diverges at the spectral
degeneracies. Probing the quantum evolution with the diagonal elements of the
time correlation matrix is discussed as a general concept and tested in the
case of a bosonic Josephson junction. It reveals for the latter characteristic
properties at the transition to Hilbert-space localization.

##  I Introduction

Symmetries play a central role in classical as well as in quantum many-body
systems. They determine the macroscopic behavior of these systems. Moreover,
symmetries of macroscopic states reflect symmetries and spontaneous symmetry
breaking of the underlying system. For instance, the groundstates of the Ising
model with ferromagnetic nearest-neighbor spin-spin coupling are
$|\uparrow,\dots,\uparrow\rangle$ and $|\downarrow,\dots,\downarrow\rangle$.
This two-fold degeneracy is a consequence of the symmetry of the Hamiltonian.
Then we experience spontaneous symmetry breaking when we add an arbitrarily
small perturbation (e.g., a small magnetic field or special boundary
conditions) that breaks the symmetry. This leads to a lifting of the
degeneracy of the two quantum states. The fact that already an arbitrarily
small symmetry breaking term can lift the macroscopic degeneracy indicates a
complex dynamics near this symmetry and the associated phase transition. In
more general terms, it is important to understand the evolution of many-body
systems near a symmetry.

In the following we will consider the unitary evolution of closed quantum
many-body systems. It is based on the idea that the extraction of information
about the quantum system in an experiment is limited. In other words, not all
properties or degrees of freedom of the quantum model are accessible by the
experiment. Typical exceptions are the return and transition probabilities for
quantum states. Moreover, quantum systems have a complex dynamics. Even for a
few particles the evolution can be quite erratic similar to a classical random
walk. Such a behavior suggests a statistical approach to extract generic
information about the quantum evolution, using averaged quantities. A
statistical approach is also supported by the fact that large sets of
experimental data are available whose properties can be treated statistically.
This idea is not new and found a very successful realization in the Random
Matrix Theory (RMT). It has been applied to many physical systems, such as
nuclei, atoms and mesoscopic systems [1, 2, 3, 4, 5, 6, 7]. The motivation for
the RMT is that there is no way of knowing the Hamiltonian of even a
relatively small many-body quantum system, such as an atomic nucleus. On the
other hand, the spectra of these systems, complex though, have some
characteristic features, such as level repulsion. Thus, instead of guessing a
specific Hamiltonian, a random ensemble of Hamiltonians is chosen, which
describes the generic features of a class of quantum systems. The class is
characterized by the invariance of the random ensemble with respect to
symmetry transformations. These are typically orthogonal, unitary or
symplectic transformations. Another application of the RMT has been recently
proposed for the description of random measurements. It is based on Dyson’s
circular matrix ensemble [3, 8, 9, 6]), which represents random unitary
matrices and has been used as a tool to determine the trace of powers of the
density matrix and the related Rényi entropy [10, 11, 12, 13, 14].

In contrast to these RMT approaches, we consider in the following a dynamical
approach in which only the time of a measurement is random, whereas the energy
levels $\\{E_{j}\\}$ of the Hamiltonian $H$ and the overlaps $\langle
E_{j}|\Psi_{0}\rangle$, $\langle E_{j}|\Psi\rangle$ of the energy eigenstates
$\\{|E_{j}\rangle\\}$ with a given initial state $|\Psi_{0}\rangle$ and a
measured state $|\Psi\rangle$ are not random. This leads to the time
correlation matrix (TCM) as the central tool for the definition of the
statistical model, instead of the random ensemble of Hamiltonians in the RMT.
We will employ this approach, which was previously described in Ref. [15], to
analyze the evolution of the return and transition probabilities. In more
concrete terms, for a given time $t_{k}$ we evaluate (in a calculation or in a
real experiment) the probability $p_{k}$ that the system is in a certain
state. Then we evaluate the probabilities $\\{p_{1},p_{2},...\\}$ at different
discrete and randomly chosen times $\\{t_{1},t_{2},...\\}$. This can be
translated into practical observations, where it is assumed that each
experiment is prepared in the same initial state and all measurements are
performed for the same final state of the evolution at different times. These
experiments provide an ensemble of probabilities $\\{p_{1},p_{2},...\\}$ with
the corresponding times $\\{t_{1},t_{2},...\\}$.

For given overlaps $\langle E_{j}|\Psi_{0}\rangle$, $\langle
E_{j}|\Psi\rangle$ we can immediately predict some restrictions for the
evolution of the probability $p_{k}$ in the $N$–dimensional Hilbert space.
When the overlaps vanish for some of the eigenstates $|E_{n}\rangle$, the
evolution cannot reach those states and the accessible Hilbert space is
restricted to the states $|E_{j}\rangle$ with $j\neq n$. This reduction of the
Hilbert space can be interpreted as Hilbert-space localization [16] or
Hilbert-space fragmentation [17]. This effect can be associated with
spontaneous symmetry breaking, induced by the choice of the initial and
measured states. In the case that the overlap to some states is not strictly
zero but very small, the access to those states may be negligible and can be
ignored. This corresponds to a complex dynamical behavior and requires a
careful analysis. It will be addressed briefly for the example of a bosonic
Josephson junction in Sect. III, where the mirror symmetry of the junction is
spontaneously broken.

This paper is organized as follows. After the definition of the TCM in Sect.
II we focus on the properties of its diagonal elements (Sect. II.1). Then the
effect of spectral degeneracies on the diagonal TCM elements are discussed in
Sect. II.2. In Sect. III we analyze the diagonal TCM elements in the specific
example of a bosonic Josephson junction.

##  II Time Correlation matrix

We consider the transition amplitude from $|\Psi_{0}\rangle$ to $|\Psi\rangle$

$u_{k}=\langle\Psi|e^{-iHt_{k}}|\Psi_{0}\rangle\ ,$ (1)

which is based on the unitary evolution with the Hamiltonian $H$ from the
initial state $|\Psi_{0}\rangle$. The probability to measure the state
$|\psi\rangle$ at time $t_{k}$ is given by
$|\langle\Psi|e^{-iHt_{k}}|\Psi_{0}\rangle|^{2}$. In other words,
$p_{k}=|u_{k}|^{2}$ is the probability to find the quantum system in the state
$|\Psi\rangle$ after the unitary evolution from the initial state
$|\Psi_{0}\rangle$ over the time $t_{k}$. Since the evolution is defined by
the Hamiltonian $H$, we consider its eigenstates
$\\{|E_{j}\rangle\\}_{j=1,...,N}$ and its corresponding eigenvalues
$\\{E_{j}\\}_{j=1,...,N}$ and write the amplitude in spectral representation
as

$u_{k}=\sum_{j=1}^{N}\langle\Psi|E_{j}\rangle\langle E_{j}|\Psi_{0}\rangle
e^{-iE_{j}t_{k}}\equiv\sum_{j=1}^{N}q_{j}e^{-iE_{j}t_{k}}\ .$ (2)

Although the phases are not directly accessible experimentally, their effect
is observable through the interference of quantum states. This enables us, for
instance, to relate the product of amplitudes at different times with
probabilities of interfering amplitudes:

$u_{k}^{*}u_{k^{\prime}}+u_{k^{\prime}}^{*}u_{k}=\frac{1}{2}\left(|u_{k}+u_{k^{\prime}}|^{2}-|u_{k}-u_{k^{\prime}}|^{2}\right)\
,\ \
i(u_{k}^{*}u_{k^{\prime}}-u_{k^{\prime}}^{*}u_{k})=\frac{1}{2}\left(|u_{k}+iu_{k^{\prime}}|^{2}-|u_{k}-iu_{k^{\prime}}|^{2}\right)\
,$ (3)

where the probabilities are observable in interferometric measurements. This
relation suggests to consider the correlation of the amplitudes $u_{k}$,
$u_{k^{\prime}}$ at different times through the TCM $\langle
u^{*}_{k}u_{k^{\prime}}\rangle_{\tau}$, where the average
$\langle...\rangle_{\tau}$ is taken with respect to the distribution of times
$\\{t_{k}\\}$ as a result of inaccurate clocks: The time is measured by a
clock in each laboratory, which counts time steps $\\{\tau_{n}\\}$. These
clocks have a limited accuracy, such that the time steps vary randomly. This
implies a sequence of measurements in each laboratory, where the clock
indicates $k$ time steps corresponding to the total evolution time
$t_{k}=\tau_{1}+\cdots+\tau_{k}$ for different values of $k$. Now we compare
the measurements of different laboratories. This provides a distribution of
results for $u^{*}_{k}u_{k^{\prime}}$ due to different inaccurate clocks,
where we assume that the fluctuations of the time steps $\\{\tau_{n}\\}$ are
independently and equally distributed. Then the TCM reads in spectral
representation

$\langle
u_{k}^{*}u_{k^{\prime}}\rangle_{\tau}=\sum_{j,j^{\prime}}q^{*}_{j}q_{j^{\prime}}\langle
e^{iE_{j}(\tau_{1}+\cdots+\tau_{k})}e^{-iE_{j^{\prime}}(\tau_{1}+\cdots+\tau_{k^{\prime}})}\rangle_{\tau}$
(4) $=\sum_{j,j^{\prime}}q^{*}_{j}q_{j^{\prime}}\cases{\langle
e^{i(E_{j}-E_{j^{\prime}})(\tau_{1}+\cdots+\tau_{k})}\rangle_{\tau}\langle
e^{-iE_{j^{\prime}}(\tau_{k+1}+\cdots+\tau_{k^{\prime}})}\rangle_{\tau}&$k^{\prime}>k$\cr\langle
e^{i(E_{j}-E_{j^{\prime}})(\tau_{1}+\cdots+\tau_{k^{\prime}})}\rangle_{\tau}\langle
e^{iE_{j^{\prime}}(\tau_{k^{\prime}+1}+\cdots+\tau_{k})}\rangle_{\tau}&$k^{\prime}<k$\cr\langle
e^{i(E_{j}-E_{j^{\prime}})(\tau_{1}+\cdots+\tau_{k})}\rangle_{\tau}&$k^{\prime}=k$\cr}\
.$ (5)

Defining $\lambda_{j}=\langle e^{iE_{j}\tau}\rangle_{\tau}$ and
$\lambda_{jj^{\prime}}=\langle e^{i(E_{j}-E_{j^{\prime}})\tau}\rangle_{\tau}$
the TCM elements become

$\langle
u_{k}^{*}u_{k^{\prime}}\rangle_{\tau}=\sum_{j,j^{\prime}}q^{*}_{j}q_{j^{\prime}}\cases{\lambda_{jj^{\prime}}^{k}\lambda_{j^{\prime}}^{k^{\prime}-k}&$k^{\prime}\geq
k$\cr\lambda_{jj^{\prime}}^{k^{\prime}}{\lambda^{*}_{j}}^{k-k^{\prime}}&$k^{\prime}<k$\cr}\
.$ (6)

The TCM decays exponentially with $|k-k^{\prime}|$, provided
$|\lambda_{j}|<1$. For fixed $|k-k^{\prime}|$ the TCM is constant for the
diagonal elements $\lambda_{jj}=1$, though. This reflects the fact that a
unitary evolution between the same energy eigenstates gives just a phase
factor $e^{-iE_{j}\tau}$ (cf. Eq. (2)). On the other hand, these phase factors
lead to a decay of different energy states due to interference effects after
the time average.

###  II.1 Diagonal elements of the TCM

The diagonal TCM element $\langle|u_{k}|^{2}\rangle_{\tau}$ is the probability
to measure the state $|\Psi\rangle$ at time $t_{k}$. Before time averaging,
the expression $|u_{k}|^{2}$ is a diagonal element of the density matrix
$\rho(t_{k})$ with respect to the state $|\Psi\rangle$. The trace of
$|u_{k}|^{2}$ with respect to all states $|\Psi\rangle$ of the underlying
Hilbert space is the spectral form factor, often used for the characterization
of many-body quantum chaos [18, 19]. We only mention this but will not study
it here.

According to Eq. (6) the average transition probability
$\langle|u_{k}|^{2}\rangle_{\tau}$ reads

$\langle|u_{k}|^{2}\rangle_{\tau}=P_{N}+\sum_{j,j^{\prime}=1;j^{\prime}\neq
j}^{N}q^{*}_{j}q_{j^{\prime}}\lambda_{jj^{\prime}}^{k}\ ,\ \
P_{N}=\sum_{j=1}^{N}|q_{j}|^{2}\ .$ (7)

The term $P_{N}$ describes the recurrent behavior, which does not depend on
time. It is the asymptotic transition probability for $k\to\infty$

$P_{N}=\lim_{k\to\infty}\langle|u_{k}|^{2}\rangle_{\tau}\ ,$ (8)

provided that the energy levels are not degenerated. The case of degenerate
energy levels is discussed in the next section. The second term in Eq. (7)
decays exponentially with time due to $|\lambda_{jj^{\prime}}|<1$, and only
this term describes a change of the transition probability during the
evolution of the quantum system. This result provides a separation of the
diagonal elements of the TCM into a static recurrent term $P_{N}$ and a
dynamic term that decays quickly.

The recurrent term $P_{N}$ stores important information regarding the
properties of the quantum system. Since
$|q_{j}|^{2}=|\langle\Psi|E_{j}\rangle|^{2}|\langle\Psi_{0}|E_{j}\rangle|^{2}$
is a product of the overlaps between the energy eigenstate $|E_{j}\rangle$
with the initial state and with the measured state, it provides a measure of
how much this energy eigenstate contributes to the transition
$|\Psi_{0}\rangle\to|\Psi\rangle$ during the unitary evolution. For instance,
the asymptotic behavior of the return probability to the initial state
$|\Psi_{0}\rangle\to|\Psi_{0}\rangle$ with the dimensionality $N$ of the
underlying Hilbert space describes Anderson localization when
$\lim_{N\to\infty}P_{N}>0$ and the absence of Anderson localization when
$\lim_{N\to\infty}P_{N}=0$ [16]. This can be understood by noting that the
normalization of quantum states implies
$\sum_{j=1}^{N}|\langle\Psi_{0}|E_{j}\rangle|^{2}=1$ and that for a localized
state only a few energy eigenstates have a nonzero overlap with
$|\Psi_{0}\rangle$. For a delocalized state, on the other hand, the overlap is
nonzero for a large number of energy eigenstates which is of the order of $N$.
An extreme case is given when these overlaps are equal. Then we have
$|\langle\Psi_{0}|E_{j}\rangle|^{2}=1/N$ due to the normalization, which
implies $P_{N}=1/N$. Anderson localization is associated with a random
Hamiltonian [20]. According to the above described picture, we can also
consider Hilbert-space localization for a deterministic Hamiltonian, which
depends strongly on the initial state. For an energy eigenstate, the system
will always remain in the latter under unitary evolution. More general, is the
initial state a superposition of $m$ energy eigenstates, the system will
always remain inside the $m$–dimensional Hilbert space, spanned by these
energy eigenstates. In the case, where the initial state is eigenstate of
$H_{0}$ of the Hamiltonian $H=H_{0}+\eta H_{1}$ and $\eta$ is a small
parameter, $\eta H_{1}$ is a small perturbation. In that case it is possible
that this perturbation provides an exponentially decaying evolution away the
initial state. This would be considered as exponential Hilbert-space
localization.

In the subsequent discussion we will focus on the diagonal elements of the
TCM, since the off-diagonal TCM elements decay exponentially with
$|k-k^{\prime}|$ according to Eq. (6).

###  II.2 Effect of spectral degeneracies

Assuming that there is a spectral degeneracy $E_{1}=E_{2}$, we have
$\lambda_{12}=\lambda_{21}=1$, and the diagonal TCM elements in Eq. (7) read
in this case

$\langle|u_{k}|^{2}\rangle_{\tau}=\sum_{j,j^{\prime}=1}^{N}q^{*}_{j}q_{j^{\prime}}\lambda_{jj^{\prime}}^{k}=P_{N}+q_{1}^{*}q_{2}+q_{2}^{*}q_{1}+\sum_{j,j^{\prime}=1;j^{\prime}\neq
j;(j,j^{\prime})\neq(1,2),(2,1)}^{N}q^{*}_{j}q_{j^{\prime}}\lambda_{jj^{\prime}}^{k}\
,$ (9)

such that the recurrent part of the transition probability becomes

$\lim_{k\to\infty}\langle|u_{k}|^{2}\rangle_{\tau}=P_{N}+q_{1}^{*}q_{2}+q_{2}^{*}q_{1}=|q_{1}+q_{2}|^{2}+\sum_{j=3}^{N}|q_{j}|^{2}\
.$ (10)

Thus, the effect of a spectral degeneracy is a change of the recurrent and the
decaying behavior, where the recurrent term changes by
$|q_{1}+q_{2}|^{2}-|q_{1}|^{2}-|q_{2}|^{2}$. This means that the diagonal
elements of the TCM are very sensitive in terms of spectral degeneracies.

After applying a discrete Fourier transformation to the decaying part of
$\langle|u_{k}|^{2}\rangle_{\tau}$ we obtain the function

${\tilde{U}}_{d}(e^{i\omega})=\sum_{k\geq 1}e^{i\omega
k}\sum_{j,j^{\prime}=1;k^{\prime}\neq
k;(j,j^{\prime})\neq(1,2),(2,1)}^{N}q^{*}_{j}q_{j^{\prime}}\lambda_{jj^{\prime}}^{k}=\sum_{j,j^{\prime}=1;j^{\prime}\neq
j;(j,j^{\prime})\neq(1,2),(2,1)}^{N}q^{*}_{j}q_{j^{\prime}}\frac{\lambda_{jj^{\prime}}}{e^{-i\omega}-\lambda_{jj^{\prime}}}\
,$ (11)

which is a function of $\omega$ on the interval $[0,2\pi)$. In other words,
${\tilde{U}}_{d}(z)$ is a sum of poles inside the unit circle due to
$|\lambda_{jj^{\prime}}|<1$. The poles $\lambda_{nn^{\prime}}$ and
$\lambda_{n^{\prime}n}$ approach the unit circle when we are getting closer to
a degeneracy of $E_{n}$ and $E_{n^{\prime}}$ This should be visible in
${\tilde{U}}_{d}(e^{i\omega})$. The corresponding decay time
$T_{d}=-1/\log|\lambda_{nn^{\prime}}|$ diverges due to
$|\lambda_{nn^{\prime}}|\sim 1$. Therefore, the decay time $T_{d}$ is a
measure for the distance from a spectral degeneracy; it diverges when we
approach the degeneracy. In general, we can define

$T_{m}=\max_{j,j^{\prime}=1,\ldots,N}-\frac{1}{\log|\lambda_{jj^{\prime}}|}$
(12)

as the largest decay time as a measure of level degeneracy.

##  III Example: bosonic Josephson junction

In this section we will study the diagonal TCM elements of a bosonic Josephson
junction (BJJ) with $N$ bosons as a closed quantum system. The motivation for
choosing this example is at least threefold: The model is (i) simple enough
but not trivial with interesting features based on tunneling and boson-boson
interaction, (ii) it can be solved exactly and (iii) it has been realized
experimentally [21, 22] with applications to commercial quantum computers
[23]. The BJJ consists of two identical wells filled with interacting bosons
and a tunneling junction between them. More formally, it is defined by the
Bose-Hubbard Hamiltonian [24]

$H=-\frac{J}{2}(a_{l}^{\dagger}a_{r}+a_{r}^{\dagger}a_{l})+U(n_{l}^{2}+n_{r}^{2}),\
\ \ n_{l,r}=a_{l,r}^{\dagger}a_{l,r}\ ,$ (13)

where $a^{\dagger}_{l,r}$ ($a_{l,r}$) are the creation (annihilation)
operators in the left and right well, respectively. The first term of $H$
describes tunneling of atoms between the wells, and for $U>0$ the second term
represents a repulsive particle-particle interaction that favors energetically
a symmetric distribution of bosons in the double well. Without tunneling
(i.e., for $J=0$) there are two-fold degenerate energy levels
$E_{k}=U[(N-k)^{2}+k^{2}]/2$ with eigenstates that are superpositions of the
product Fock state $|k,N-k\rangle$ ($\equiv|k\rangle\otimes|N-k\rangle$) and
its mirror image $|N-k,k\rangle$. This two-fold degeneracy is similar to the
two-fold degeneracy of the Ising model, mentioned in the Introduction. Thus,
tunneling between the wells plays the role of the symmetry-breaking term. In
contrast to the macroscopic Ising model though, an arbitrarily small tunneling
term may not be sufficient to cause symmetry breaking. Hence the following
study is more related to an Ising model with a finite number of spins than to
a macroscopic system. Nevertheless, the sensitive dynamics due to tunneling
between the two wells with degenerate energy levels will give us an insight
into the evolution near symmetry points. Another difference between the Ising
model and the BJJ is that without tunneling in the latter all energy levels
are two-fold degenerate. Therefore, the initial state can be prepared in any
of these energy levels to follow the evolution due to tunneling in the
vicinity of degenerate levels. This is important because the degenerate
groundstate may not be reached due to the energy conservation in the closed
quantum system.

For the following we will use $|\Psi_{0}\rangle=|0,N\rangle$ as the initial
state and $|\Psi\rangle=|N,0\rangle$ as the measured state. Then we define the
return probability (RP) and the transition probability (TP) as

$|u_{r,k}|^{2}=|\langle 0,N|e^{-iHt_{k}}|0,N\rangle|^{2}\ ,\ \ \
|u_{t,k}|^{2}=|\langle N,0|e^{-iHt_{k}}|0,N\rangle|^{2}\ .$ (14)

Since both states $|\Psi_{0}\rangle$, $|\Psi\rangle$ are eigenstates of $H$ in
the absence of tunneling ($J=0$), we get

$|u_{r,k}|^{2}=1\ ,\ \ \ |u_{t,k}|^{2}=0$ (15)

for any $k$ due to orthogonality. This reflects the fact that this pair of
states breaks the mirror symmetry of the double well. The opposite extreme is
the BJJ without boson-boson interaction ($U=0$), which is more complex and
will be discussed in the next subsection. As we will see, this case can be
described by simple functions for $|u_{r,k}|^{2}$, $|u_{t,k}|^{2}$. For the
interplay of tunneling and boson-boson interaction ($J,U\neq 0$) the behavior
is more complex and we will rely on the time averaged expressions
$\langle|u_{r,k}|^{2}\rangle_{\tau}$, $\langle|u_{t,k}|^{2}\rangle_{\tau}$
with $\tau_{k}={\bar{\tau}}+\tau^{\prime}_{k}$ and an exponential distribution
for $\tau^{\prime}_{k}$:

$\langle...\rangle_{\tau}=\int_{0}^{\infty}...\prod_{n\geq
1}e^{-\tau^{\prime}_{n}}d\tau^{\prime}_{n}\ ,$ (16)

where the time is measured in units of $\hbar/J$.

###  III.1 Non-interacting bosons

For $U=0$ the spectrum of $H$ consists of equidistant energy levels
$E_{j}=-J(N/2-j)$ ($j=0,1,...,N$) and eigenstates

$|E_{j}\rangle=\frac{2^{-N/2}}{\sqrt{j!(N-j)!}}(a_{l}^{\dagger}+a_{r}^{\dagger})^{j}(a_{l}^{\dagger}-a_{r}^{\dagger})^{N-j}|0,0\rangle\
.$

Then the RP and the TP at time $t_{k}$ for $N$ bosons read

$|u_{r,k}|^{2}=|\langle
0,N|e^{-iHt_{k}}|0,N\rangle|^{2}=|\cos^{N}(Jt_{k}/2)|^{2}\ ,\ \ \
|u_{t,k}|^{2}=|\langle
N,0|e^{-iHt_{k}}|0,N\rangle|^{2}=|\sin^{N}(Jt_{k}/2)|^{2}\ .$ (17)

For the average TP we get with
$t_{k}=k{\bar{\tau}}+\tau^{\prime}_{1}+\cdots+\tau^{\prime}_{k}$ and the
exponential distribution of Eq. (16)

$\langle|u_{t,k}|^{2}\rangle_{\tau}=2^{-2N}\sum_{l,l^{\prime}=0}^{N}{{N}\choose{l}}{{N}\choose{l^{\prime}}}(-1)^{l+l^{\prime}}\left[\frac{e^{-iJ{\bar{\tau}}(l-l^{\prime})}}{1-iJ(l-l^{\prime})}\right]^{k}$
$=2^{-2N}\sum_{l=0}^{N}{{N}\choose{l}}^{2}+2^{-2N}\sum_{l,l^{\prime}=0;l^{\prime}\neq
l}^{N}{{N}\choose{l}}{{N}\choose{l^{\prime}}}(-1)^{l+l^{\prime}}\left[\frac{e^{-iJ{\bar{\tau}}(l-l^{\prime})}}{1-iJ(l-l^{\prime})}\right]^{k}\
,$ (18)

while for the average RP is the same without the factor $(-1)^{l+l^{\prime}}$.
From these results we get for the asymptotic value at $k\sim\infty$ for both
probabilities

$P_{N}\sim\langle|u_{r,k}|^{2}\rangle_{\tau}\sim\langle|u_{t,k}|^{2}\rangle_{\tau}\sim
2^{-2N}\sum_{l=0}^{N}{{N}\choose{l}}^{2}\ ,$ (19)

which reflects the mirror symmetry of the BJJ.

Figure 1:  Signatures of a qualitative change of the evolution of a bosonic
Josephson junction upon an increasing interaction strength $u$. The plots
represent the dynamics of 20 bosons for $u=1,2$, where the top panel gives the
return probability $|u_{r,k}|^{2}$ and the bottom panel the transition
probability $|u_{t,k}|^{2}$.

###  III.2 Interacting bosons

Exact solutions exist for this model also for $J,U\neq 0$ but in contrast to
the non-interacting case they are complex and difficult to present in general.
For instance, the resolvent is a meromorphic function with polynomials of
order $N$ and $N+1$ [16, 25]. Therefore, we only plot the results for the RP
and the TP and their averaged counterparts here. For transparency, we choose
for all subsequent plots $N=20$ bosons.

First, the evolution of the RP and the TP are presented in Fig. 1 for two
values $u=1,2$ of the interaction parameter $u=NU/J$. This clearly reveals
that the RP dominates over the TP for increasing $u$, as we expect from the
results of the two limits $J=0$ and $U=0$. It is interesting to note that in a
mean-field (classical) approximation of the BJJ there is a sharp phase
transition in terms of the interaction parameter, where the mean-field TP is
completely suppressed when $u\geq u_{c}=2$ [26]. The strong interaction phase
is also called the self-trapping phase. The analogue of the latter in the
quantum BJJ is Hilbert-space localization, reflected by the scaling behavior
of the inverse participation ratio [16]. This also indicates the existence of
a critical $u_{c}$.

In Fig. 2 the effect of time averaging on $|u_{r,k}|^{2}$ and $|u_{t,k}|^{2}$
for $u=1$ is visualized. It reflects the smoothing of the strongly fluctuating
dynamics with a recurrent and a decaying contribution according to Eqs. (7),
(9). It is obvious that the separation of the recurrent and the decaying
behavior is not accessible without time averaging.

The existence of a critical interaction strength $u_{c}\approx 1.89...$ is
demonstrated in Fig. 3, where the $\langle|u_{r,k}|^{2}\rangle_{\tau}$ jumps
up upon increasing $u$ at $u_{c}$. Moreover,
$\langle|u_{t,k}|^{2}\rangle_{\tau}$ develops a characteristic peak at
$u_{c}$. This behavior reflects the appearance of nearly degenerate energy
levels, as described in Sect. II.2.

Finally, in Fig. 4 the change of the time scales for the decay of the average
TP $\langle|u_{t,k}|^{2}\rangle_{\tau}$ is visualized for $u=1.7,...,2.2$. The
decay is reduced by an increasing interaction strength $u$. This reflects the
fact that the splitting of the energy levels is reduced by the interaction, as
we would have expected.

Figure 2:  Comparison of the return probability $|u_{r,k}|^{2}$ and the
average return probability $\langle|u_{r,k}|^{2}\rangle_{\tau}$ (top panel)
and of the corresponding transition probabilities (bottom panel). The average
was taken with respect to the exponential distribution of Eq. (16). The
interaction parameter is $u=1$ and ${\bar{\tau}}=1/10$.  Figure 3:  The
critical regime of the Hilbert-space localization with $u_{c}\approx 1.89$ is
visualized with $\langle|u_{r,k}|^{2}\rangle_{\tau}$ and
$\langle|u_{t,k}|^{2}\rangle_{\tau}$ at $k=70$.  Figure 4:  Decay of
$\langle|u_{t,k}|^{2}\rangle_{\tau}$ for different interaction parameters
$u=1.7,\ldots,2.2$.

##  IV Discussion and conclusions

Our analysis of the quantum unitary evolution is strictly focused on the
result of a single measurement in each of many identical experiments, which
are subject of a unitary evolution. Averaging with respect to the statistical
outcome due to measurements at randomly distributed times leads to the TCM. We
have focused on the diagonal TCM elements to study the evolution of the
quantum system. The analysis of the off-diagonal TCM elements was the subject
of a previous work [15]. Similar to the off-diagonal TCM elements, the
diagonal TCM elements reveal a separation of the evolution into a static
recurrent part and a dynamic decaying part. We have found that the decay rate
of the latter is related to the spacing between energy levels, which diverges
when the spacing vanishes near a degeneracy. Thus, the decay rate is a
quantity that can be used to detect symmetry changes or the appearance of
spontaneous symmetry breaking. This has been observed in the example of the
BJJ: In this model the energy levels are two-fold degenerate in the limit
$u\to\infty$. This is reflected in Fig. 4, where the decay decreases with
increasing $u$.

Another interesting aspect of the BJJ is the transition to Hilbert-space
localization [16]. This transition was also detected with the helps of the
average RP and average TP in Fig. 3, where the average RP experience a jump to
a higher value for $u>u_{c}$. On the other hand, the average TP has only a
sharp peak near $u_{c}$ but has the same value away from $u_{c}$.

We can conclude that time averaging over an ensemble of measurements is
crucial for extracting the (static) recurrent behavior and the (dynamic)
decaying behavior. This can be formulated in terms of the TCM. The separation
of static and dynamic behavior is generic for the unitary quantum evolution.
Then the TCM provides a tool to analyze properties of the quantum system,
which are associated with spectral degeneracy in the case of phase
transitions. It can be applied to theoretical calculations as well as to
experimental data, collected from many experiments at identical quantum
systems. We have demonstrated in the case of the BJJ that the TCM approach
delivers interesting generic information. The BJJ can be considered as a
building block of the Bose-Hubbard model on a lattice. Therefore, the TCM
approach should be applicable to more complex quantum models, including
bosonic and fermionic Hubbard models as well as quantum spin systems.

###### Acknowledgements.

I am grateful to Eli Barkai for interesting discussions regarding the Hilbert-
space fragmentation. This work was supported by the Julian Schwinger
Foundation.

## References

  * Wigner [1955] E. P. Wigner, Characteristic vectors of bordered matrices with infinite dimensions, Annals of Mathematics 62, 548 (1955).
  * Porter and Thomas [1956] C. E. Porter and R. G. Thomas, Fluctuations of nuclear reaction widths, Phys. Rev. 104, 483 (1956).
  * Dyson [1962a] F. J. Dyson, Statistical theory of the energy levels of complex systems. i, Journal of Mathematical Physics 3, 140 (1962a), https://doi.org/10.1063/1.1703773 .
  * Porter [1966] C. E. Porter, Book Review: Statistical theories of spectra: fluctuations. C.E. PORTER (Academic Press, New York, 1965\. xv-576 p. 5.95 paper, 9.50 cloth), Nuclear Physics 78, 696 (1966).
  * Wigner [1967] E. P. Wigner, Random matrices in physics, SIAM Review 9, 1 (1967).
  * Mehta [2004] M. Mehta,  _Random Matrices_, Pure and applied mathematics No. v. 142 (Elsevier/Academic Press, 2004).
  * Beenakker [1997] C. W. J. Beenakker, Random-matrix theory of quantum transport, Rev. Mod. Phys. 69, 731 (1997).
  * Dyson [1962b] F. J. Dyson, Statistical theory of the energy levels of complex systems. ii, Journal of Mathematical Physics 3, 157 (1962b), https://doi.org/10.1063/1.1703774 .
  * Dyson [1962c] F. J. Dyson, Statistical theory of the energy levels of complex systems. iii, Journal of Mathematical Physics 3, 166 (1962c), https://doi.org/10.1063/1.1703775 .
  * van Enk and Beenakker [2012] S. J. van Enk and C. W. J. Beenakker, Measuring $\mathrm{Tr}{\rho}^{n}$ on single copies of $\rho$ using random measurements, Phys. Rev. Lett. 108, 110503 (2012).
  * Elben  _et al._ [2018] A. Elben, B. Vermersch, M. Dalmonte, J. I. Cirac, and P. Zoller, Rényi entropies from random quenches in atomic hubbard and spin models, Phys. Rev. Lett. 120, 050406 (2018).
  * Vermersch  _et al._ [2018] B. Vermersch, A. Elben, M. Dalmonte, J. I. Cirac, and P. Zoller, Unitary $n$-designs via random quenches in atomic hubbard and spin models: Application to the measurement of rényi entropies, Phys. Rev. A 97, 023604 (2018).
  * Li  _et al._ [2019] Y. Li, X. Chen, and M. P. A. Fisher, Measurement-driven entanglement transition in hybrid quantum circuits, Phys. Rev. B 100, 134306 (2019).
  * Skinner  _et al._ [2019] B. Skinner, J. Ruhman, and A. Nahum, Measurement-induced phase transitions in the dynamics of entanglement, Phys. Rev. X 9, 031009 (2019).
  * Ziegler [2021] K. Ziegler, Quantized dynamics in closed quantum systems, Journal of Physics A: Mathematical and Theoretical 54, 205303 (2021).
  * Cohen  _et al._ [2016] D. Cohen, V. I. Yukalov, and K. Ziegler, Hilbert-space localization in closed quantum systems, Phys. Rev. A 93, 042101 (2016).
  * Moudgalya and Motrunich [2021] S. Moudgalya and O. I. Motrunich, Hilbert space fragmentation and commutant algebras (2021), arXiv:2108.10324 [cond-mat.stat-mech] .
  * Kos  _et al._ [2018] P. Kos, M. Ljubotina, and T. c. v. Prosen, Many-body quantum chaos: Analytic connection to random matrix theory, Phys. Rev. X 8, 021062 (2018).
  * Chan  _et al._ [2018] A. Chan, A. De Luca, and J. T. Chalker, Solution of a minimal model for many-body quantum chaos, Phys. Rev. X 8, 041019 (2018).
  * Anderson [1958] P. W. Anderson, Absence of diffusion in certain random lattices, Phys. Rev. 109, 1492 (1958).
  * Bar-Gill  _et al._ [2009] N. Bar-Gill, G. Kurizki, M. Oberthaler, and N. Davidson, Dynamic control and probing of many-body decoherence in double-well bose-einstein condensates, Phys. Rev. A 80, 053613 (2009).
  * Juliá-Díaz  _et al._ [2012] B. Juliá-Díaz, T. Zibold, M. K. Oberthaler, M. Melé-Messeguer, J. Martorell, and A. Polls, Dynamic generation of spin-squeezed states in bosonic josephson junctions, Phys. Rev. A 86, 023615 (2012).
  * Barends  _et al._ [2013] R. Barends, J. Kelly, A. Megrant, D. Sank, E. Jeffrey, Y. Chen, Y. Yin, B. Chiaro, J. Mutus, C. Neill, P. O’Malley, P. Roushan, J. Wenner, T. C. White, A. N. Cleland, and J. M. Martinis, Coherent josephson qubit suitable for scalable quantum integrated circuits, Phys. Rev. Lett. 111, 080502 (2013).
  * Gati and Oberthaler [2007] R. Gati and M. K. Oberthaler, A bosonic josephson junction, Journal of Physics B: Atomic, Molecular and Optical Physics 40, R61 (2007).
  * Ziegler [2011] K. Ziegler, Dynamical creation of entangled bosonic states in a double well, Journal of Physics B: Atomic, Molecular and Optical Physics 44, 145302 (2011).
  * Milburn  _et al._ [1997] G. J. Milburn, J. Corney, E. M. Wright, and D. F. Walls, Quantum dynamics of an atomic bose-einstein condensate in a double-well potential, Phys. Rev. A 55, 4318 (1997).