text:
# Stadium Billiard with Moving Walls

## I Introduction

Consider the problem of particle in a box, where some piece of the wall is deformed periodically in time. As an example one may think of a particle in a cylinder with a moving piston. The particle has mass $`m`$ and kinetic energy $`E`$. The piston is pushed back and forth. The velocity in which the piston is displaced is $`\pm V`$, and the maximum displacement is $`A`$. At the end of each cycle the piston is back in its original location.

In the present Paper we are going to study what happens to a quantum mechanical particle in such a box, during one cycle of the driving. We assume that the particle is initially prepared in an energy eigenstate of the (un-deformed) box. We shall explain that it is important to specify whether the maximum displacement $`A`$ is less or larger compared with De-Broglie wavelength.

The problem of particle-in-a-box with a moving wall is a prototype example for study of driven systems lecnotes  that are described by a Hamiltonian $`(x(t))`$, where $`x(t)`$ is a time dependent parameter. In case of a piston, the parameter $`x`$ is the position of the piston, and we use the notation $`V=|\dot{x}|`$.

There is a lot of interest today in studies of quantum irreversibility. This is a relevant issue to the design of quantum computers, where the fidelity of a driving cycle is important. If at the end of a driving cycle the system is back in its initial state, then we say that the driving cycle has high fidelity. Obviously, the piston model can serve as a prototype example for studying of quantum irreversibility dsp .

The study of one-cycle driving also constitutes a bridge to the study of the response to multi-cycle (periodic) driving. It should be clear that the long time behavior of driven systems is determined by the short time dynamics. Therefore it is essential to have a good understanding of the latter.

At first sight one may think that it is simplest to study the one-dimensional (1D) box case, also known as ‘infinite square well potential with moving wall’ infwell\_a ; infwell\_b . The case of periodic driving is also known as the ‘Fermi acceleration problem’ jose . In a second sight one find out that the 1D-box case is actually the most complicated one. As in the case of the kicked rotator (standard map) qkr  there is a complicated route to chaos and stochasticity.

Things become much simpler if the motion in the box is chaotic to begin with. Driven chaotic systems exhibit, as a result of the driving, stochastic energy spreading of relatively simple nature lecnotes . This is the case that we want to consider in this Paper. Hence we consider the simplest ”chaotic box”, which is a two dimensional (2D) ”billiard” model. Specifically, we are going to introduce a detailed numerical study of the one-pulse response of a stadium billiard to a shape deformation. The stadium billiard is a recognized prototype system for ”quantum chaos” studies. Our simulations are feasible thanks to a new powerful technique for finding clusters of billiard eigenstates VS ; vergini\_thesis ; barnett\_thesis . Previous applications of this technique, to the study of restricted aspects of the present problem, have been reported in diego  and in prm .

The work which is presented in this Paper, is the first numerical study which is aimed in presenting a systematic analysis of non-perturbative features of a time dependent spreading process crs ; rsp ; frc . A procedure is developed for the determination of the various time stages in the evolution of the energy distribution. This allows the identification of the various $`V`$ regime (”adiabatic”, ”perturbative”, ”non-perturbative”, and ”sudden”) which were predicted in past theoretical studies.

## II Outline

In Sections 3 and Section 4 we define the model system, and briefly describe the classical picture. In particular we explain that the time dependent features that we study in this Paper are purely quantum-mechanical, and not of semiclassical origin.

In Section 5 we define the main object of our study, which is the energy spreading kernel $`P_t(n|m)`$. See Eq.(4). Regarded as a function of the level index $`n`$, it is the energy distribution after time $`t`$, while $`m`$ is the initial level. We also define the sqrt of the variance $`\delta E(t)`$, and the $`50\%`$ probability width $`\mathrm{\Gamma }(t)`$, that characterize the evolving energy distribution.

Our three phase strategy for analysis of energy spreading is presented in Section 6. The three phases are:     Study of the band profile.    Study of parametric evolution.    Study of the actual time evolution.    The relevant information regarding ”Phase I” is summarized in Section 7. Various approximations for $`P_t(n|m)`$, and in particular the notion of ”parametric evolution” are presented in Section 8. The relevant information regarding ”Phase II” is summarized in Section 9. The main concern of this Paper is with ”Phase III”.

For a given $`V`$ one should be able to characterize the nature of the dynamical scenario. The theoretical considerations regarding this issue have been discussed in crs ; rsp ; frc ; vrn , and for a concise review see dsp . Rather than duplicating these discussions, we are going to present in Section 10 a phenomenological definition, and a practical procedure, for the identification of the $`V`$ regimes. We would like to emphasize that this is the first time that the different dynamical scenarios (corresponding to the different $`V`$ regimes) are illustrated in a numerical simulation. The only other numerical studies (mainly diego ; rsp ) were too restricted in scope, and did not contain analysis of the stages in the evolution of the energy distribution.

Section 11 and Section 12 question the applicability of linear response theory (LRT) to the analysis of the energy spreading. Further discussion of ”non perturbative response”, and conclusions, are presented in Section 13.

## III The semiclassical picture

Consider a classical particle inside a box. Its kinetic energy is $`E`$ and the corresponding velocity is $`v_E=\sqrt{2E/m}`$. The shape of the box is externally controlled. The control parameter is denoted by $`x`$. We assume that $`x`$ has units of length, such that $`V=|\dot{x}|`$ is the typical wall velocity. Obviously different parts of the wall may have different velocities. In case of the ”piston model” only one piece of wall is moving (either inward or outward). However, we are not interested in the trivial conservative work which is being done, but only in the irreversible work. Therefore, rather than analyzing an actual ”piston model” configuration, it is wiser to consider a volume preserving deformation. In such case the conservative work is zero, and the major issue is the irreversible effect which is explained below.

As a results of the collisions of the particle with the deforming walls, there is a stochastic-like diffusion in energy space. The explanation is as follows: Each time that the particle collides with the moving wall it either gains or looses energy. To simplify the presentation let us assume head-on collisions. The change in energy is $`\pm 2mv_EV`$ depending on whether the wall is moving inward or outward at the point of the collision. For volume preserving deformation the ergodic average over $`\delta E`$ gives zero. Thus we have random-walk in energy space where the steps are $`\pm 2mv_EV`$. This leads to diffusion in energy space. As explained in frc  this diffusion is biased (the diffusion is stronger for larger $`E`$). This leads, in the long run, to a systematic (irreversible) increase of the average energy. For a more detailed presentation that takes the box geometry into account see wlf .

Do we have a corresponding (semiclassical) picture in the quantum-mechanical case? Again we remind the reader that we assume a volume preserving deformation. This implies that the energy levels of the system do not have a collective ”upward” or downward” change as a result of the deformation. The physics in which we are interested is related to the transitions between different levels. The question that we ask is whether these transitions are ”classical-like”.

Let us assume that we start with an eigenstate whose energy is $`E`$. Semiclassically it is as if we start with a microcanonical preparation. If the dynamics is of classical nature, then we expect that after a short time some of the probability will make a transition $`EE^{}`$ such that $`|E^{}E|2mv_EV`$. Naturally such description is meaningful only if the energy scale $`2mv_EV`$ is much larger than the mean level spacing. This leads for 1D-box to the condition

$`V>\mathrm{}/mL`$ (1)
In the general case (eg 2D box), having $`2mv_EV`$ much larger than the mean level spacing is not a sufficient condition for getting semiclassical behavior. Still, non-trivial analysis wld  reveals that the condition for getting semiclassical spreading in the general case is the same as in the 1D case, namely given by Eq.(1).

The above Eq.(1) is a necessary condition for having semiclassical transitions. In order to actually witness semiclassical transitions it is also necessary to have a time period much larger compared with the ballistic time ($`A/VL/v_E`$), and to have an amplitude that is much larger compared with De-Broglie wavelength ($`A\mathrm{}/(mv_E)`$). These additional conditions can be satisfied only in the semiclassical regime which is defined by Eq.(1), else they are not compatible.

## IV The numerical model

As a specific example for chaotic box, we consider the quarter stadium billiard. We define $`x`$ as the length of the straight edge, and adjust the radius parameter such that the total area is kept constant. For numerical reasons we do not analyze this model ”literally” but rather consider, as in previous study diego , a linearized version of the quarter stadium billiard Hamiltonian. Namely, we study a model Hamiltonian that has the matrix representation

$`(x(t))𝐄+\delta x(t)𝐅`$ (2)
Here $`𝐄`$ is an ordered diagonal matrix, that consists of the eigen-energies $`E_n`$ of the quarter stadium billiard with straight edge $`x=1`$. The eigen-energies were determined numerically. The perturbation due to $`\delta x`$ deformation, is represented by the matrix $`𝐅`$. Also this matrix has been determined numerically as explained in diego . We note that the fingerprints of the classical chaos are present in the statistical properties of the matrices $`𝐄`$ (level statistics) and $`𝐅`$ (band structure). The latter is discussed in Section 7.

The linearization Eq.(2) of the billiard Hamiltonian can be regarded as a valid approximation if the wall displacement parameter $`\delta x`$ is small compared with De-Broglie wavelength. This automatically excludes the possibility of addressing the semiclassical regime which has been discussed in the previous section. In our simulation the De-Broglie wavelength is roughly $`0.1`$. This implies that (at best) the maximum driving amplitude that can be allowed is $`A=0.2`$, so as to have $`|\delta x(t)|<0.1`$. \[Note that for presentation purpose we later re-define $`\delta x(t)`$ as $`\delta x(t)\delta x(0)`$.\] In the simulations we indeed have $`A=0.2`$, meaning that we allow a deformation which is comparable or possibly somewhat larger than that ”allowed” by the linearization.

However, inspite of the fact that Eq.(2) may be a ”bad” or even ”inadequate” approximation to real-world physics, it is still a totally ”legitimate” Hamiltonian from mathematical point of view. Moreover, the numerical model Eq.(2) contains all the physical ingredients that are relevant for the aim of the present study.

To summarize: In this paper we consider, as far as formulation is concerned, a chaotic billiard driven by volume preserving shape deformation. On the other hand, as far as numerics is concerned, we analyze a specific quantum mechanical model which is defined by Eq.(2). The numerical model is motivated by the quarter stadium billiard system Hamiltonian, but still it is not literally the same model. In particular, Eq.(2) does not possess the the semiclassical regime which has been discussed in the previous section.

## V The evolving energy distribution

In case of time independent Hamiltonian, the energy distribution does not change with time. In order to have an ”evolving” energy distribution, we have to make $`\delta x(t)`$ time dependent. One possibility is to assume linear driving. In such case we write $`\delta x(t)=Vt`$. But more generally, for a cycle, we write $`\delta x(t)=A\times f(t)`$, with the convention $`f(0)=f(T)=0`$. In some equations below, whenever a linear driving is concerned, $`A`$ can be replaced by $`V`$, which assumes the particular choice $`f(t)=t`$. For practical purpose, given $`f(t)`$, it is convenient to associate with it the a spectral function

$`\stackrel{~}{F}_t(\omega )=\left|{\displaystyle _0^t}\dot{f}(t^{})\text{e}^{i\omega t^{}}𝑑t^{}\right|^2`$ (3)
Some useful cases are summarized in Appendix A. In this Paper we are primarily interested in the case of triangular pulse Eq.(A).

Given the model Hamiltonian Eq.(2), and the driving scheme $`\delta x(t)`$, we can calculate the unitary evolution operator $`U(t)`$. This kernel propagates ”wavefunctions” in time. Equivalently, we can describe the quantum mechanical state using a probability density matrix. The propagator of the latter is denoted by $`\mathrm{\Lambda }(t)`$. In order to describe the evolution of the energy distribution we define the kernel

$`P_t(n|m)`$ $`=`$ $`\text{trace}(\rho _n\mathrm{\Lambda }(t)\rho _m)`$ (4)
$`=`$ $`|n(x(t))|U(t)|m(x(0))|^2`$
where $`\rho _n`$ can be interpreted as either the probability matrix or as the corresponding Wigner function that represents the eigenstate $`|n`$. In the latter case the trace operation should be interpreted as $`dPdQ/(2\pi \mathrm{})^d`$ integration over phase space, where $`(Q,P)`$ are the canonical coordinates of the particle, and $`d=2`$ is the dimensionality.

We shall use the notation $`P_t(r)=P_t(nm)=P_t(n|m)`$, with implicit average over the reference state $`m`$. We shall refer to $`P_t(r)`$ as the ”average spreading profile”. Whenever we have a wide distribution, we disregard the distinction between ”energy difference” and ”level difference”, and make the identification

$`E_nE_m\mathrm{}\omega _{nm}\mathrm{\Delta }\times r`$ (5)
where $`\mathrm{\Delta }`$ is the mean level spacing, and $`r=nm`$. (In our numerics $`\mathrm{\Delta }7.22`$).

Fig.1 displays the evolution as a function of time. Each column is a profile of the probability distribution $`P_t(r)`$. The first 7 time steps are log spaced, while the rest are linearly spaced. The $`V=100`$ evolution is predominantly parametric. As $`V`$ becomes smaller and smaller the deviation from parametric evolution becomes larger and larger. Some representative spreading profiles are presented in Fig.2.

There are various practical possibilities available for the characterization of the distribution $`P_t(r)`$. It turns out that the major features of this distribution are captured by the following three measures:

$`𝒫(t)`$ $`=`$ $`P_t(r=0)`$ (6)
$`\mathrm{\Gamma }(t)`$ $`=`$ 50% probability width (7)
$`\delta E(t)`$ $`=`$ $`\mathrm{\Delta }\times \left({\displaystyle \underset{r}{}}r^2P_t(r)\right)^{1/2}`$ (8)
The first measure is the survival probability $`𝒫(t)`$. It is mentioned here just for completeness of presentation. The second measure $`\mathrm{\Gamma }(t)`$ is the energy width of the central $`r`$ region that contains 50% of the probability. It is calculated as $`\mathrm{\Gamma }(t)=(r_{75\%}r_{25\%})\mathrm{\Delta }`$, where $`r_{25\%}`$ and $`r_{75\%}`$ are the values for which the cumulative energy distribution equals $`25\%`$ and $`75\%`$ repectively. Finally the energy spreading $`\delta E(t)`$ is defined above as the square-root of the variance.

Fig.3 shows how the width $`\mathrm{\Gamma }(t)`$ of the profile evolves. The lower panel of Fig.3 is a log-log plot. We see clearly that up to $`\delta x_c=0.006`$ the width is one level, which is an indication for the applicability of standard first-order perturbation theory. For larger $`\delta x`$ several levels are mixed non-perturbatively, and therefore the width becomes larger than one.

Fig.4 shows how the spreading $`\delta E(t)`$ of the profile evolves. The lower panel displays the relative spreading, which is defined as the ratio between the actual spreading and the parametric one. The notion of ”parametric evolution” is defined in the next paragraph. As $`V`$ becomes smaller, the departure from the parametric behavior happens earlier.

The ”sudden limit” ($`V\mathrm{}`$) of $`P_t(n|m)`$ will be denoted by $`P(n|m)`$. We shall refer to it as the ”parametric kernel”. It is formally obtained by the replacement $`\mathrm{\Lambda }(t)1`$, or equivalently $`U(t)1`$ in Eq.(4), namely,

$`P(n|m)`$ $`=`$ $`\text{trace}(\rho _n\rho _m)`$ (9)
$`=`$ $`|n(x(t))|m(x(0))|^2`$
Obviously its dependence on $`t`$ is exclusively via $`\delta x=x(t)x(0)`$, irrespective of $`V`$. The parametric evolution of $`P(r)`$ versus $`\delta x`$ for a deforming billiard has been studied in prm . The numerics in this Paper can be regarded as the extension of the numerics of prm  to the case of finite $`V`$.

## VI Three phase strategy

In the future we want to have a theory that allows the prediction of the (numerically simulated) time evolution. The minimum input which is required for such theory is the band profile of the $`𝐅`$ matrix. For precise definition of the bandprofile see Appendix B. Note the existence of a very efficient semiclassical recipe to find the band profile. This can save the need for a tedious quantum mechanical calculations.

Using the bandprofile one hopes to be able, in the ”second phase” lds , to calculated the parametric kernel $`P(n|m)`$. The bandprofile does not contain information about the correlations between the off-diagonal elements. Therefore one has to make the so-called RMT conjecture, namely to assume that the off-diagonal elements are effectively uncorrelated. It turns out that, upon using such conjecture, the non-universal features of the parametric evolution are lost. Still, the obtained results are qualitatively correct, and therefore such an approach is legitimate as an approximation.

However, finding the parametric kernel $`P(n|m)`$, using the bandprofile as an input, is not the subject of this paper. Therefore, as a matter of strategy, we would like to take the numerically determined parametric evolution as an input. Given the parametric $`P(n|m)`$, the question is whether we can calculate the actual evolution $`P_t(n|m)`$ for any finite $`V`$.

In previous works (see review in dsp ) we gave a negative answer to the above question. We have claimed that non-perturbative features of the dynamics cannot be deduced from the parametric analysis. To explain this observation let us assume that we have two model Hamiltonians, say $`_{\text{physical}}(x)`$ and $`_{\text{artificial}}(x)`$. Let us assume further that the two models have the same bandprofile and the same parametric kernel $`P(n|m)`$, Still we claim that for finite $`V`$ the two Hamiltonians may generate different temporal kernels $`P_t(n|m)`$. In particular, it has been argued that this is in fact the case if $`_{\text{artificial}}(x)`$ is an effective RMT model which is associated with $`_{\text{physical}}(x)`$.

In past publications the manifestation of non-perturbative features in case of driven physical (non-RMT) models has not been investigated numerically. On the theoretical side it is an open subject for further research rsp . The purpose of the following sections is to present a strategy for analysis of numerical simulations, that paves the way towards a theory for the non-perturbative aspects of the energy spreading.

## VII The bandprofile

The bandprofile of $`𝐅`$ is described by the spectral function $`\stackrel{~}{C}(\omega )`$, which is a Fourier transform of a correlation function $`C(\tau )`$. See Appendix B. If the collisions are uncorrelated, as in the case where the deformation involves only a small surface element, then $`C(\tau )`$ is a delta function, and the band profile is flat. This is not the case in the present model. The correlations between collisions cannot be neglected, and therefore $`C(\tau )`$ equals delta function plus a smooth component wlf . Due to the smooth component the band profile has a very pronounced non-universal structure. See Fig.5.

For large enough $`\omega `$ the contribution of the non-universal component vanishes. Thus the bandprofile should possess flat tails. These tails reflect the presence of the delta-function component in $`C(\tau )`$. For numerical reasons, due to truncation, the bandprofile of our system is multiplied by a Gaussian envelope. Therefore the tails of the effective bandprofile are not flat, but rather vanishingly small. The lack of flat tails in the numerical model can be loosely interpreted as having ’soft’ rather than ’hard’ walls.

Thus our system is characterized by a finite bandwidth. This is actually the generic case. Namely, for generic (”smooth”) Hamiltonian the correlation function $`C(\tau )`$ is non-singular, which implies finite bandwidth.

## VIII Approximations

A major issue in the studies of energy spreading is the knowledge how to combine tools or approximations in order to understand or calculate the kernel $`P_t(n|m)`$. In particular we have the following approximations:     The perturbative kernel $`P_t^{\text{prt}}(n|m)`$.    The semiclassical kernel $`P_t^{\text{scl}}(n|m)`$.    The Gaussian kernel $`P_t^{\text{sto}}(n|m)`$.    The parametric kernel $`P(n|m)`$.    The parametric kernel $`P(n|m)`$, that corresponds to the sudden limit ($`V\mathrm{}`$) has already been defined in Eq.(9). The other kernels will be defined in the present section. Obviously the kernels listed above are very different. Our aim is to clarify in what regime ($`V`$), in what time stage ($`t`$), and in what energy region ($`r`$), which of them is the valid approximation.

The semiclassical kernel $`P_t^{\text{scl}}(n|m)`$ is defined and obtained by assuming in Eq.(4) that the Wigner functions can be approximated by smeared microcanonical distributions. See Ref.vrn ; lds ; prm  for further details and numerical examples. Whenever $`P_t(n|m)P_t^{\text{scl}}(n|m)`$ we say that the spreading profile is ”semiclassical”. In order to have a valid semiclassical approximation in case of billiard systems the displacement $`\delta x`$ of the walls should be much larger than De-Broglie wavelength prm . Thus, for reasons that were explained in Section 3, the semiclassical approximation is not applicable for the numerical model that we consider in this Paper.

The perturbative kernel $`P_t^{\text{prt}}(n|m)`$ is obtained from perturbation theory frc ; vrn . The defining expression is:

$`P_t^{\text{prt}}(n|m)=`$ $`A^2\times \stackrel{~}{F}_t\left({\displaystyle \frac{E_nE_m}{\mathrm{}}}\right)\times `$
$`{\displaystyle \frac{\mathrm{\Delta }}{2\pi \mathrm{}}}\stackrel{~}{C}\left({\displaystyle \frac{E_nE_m}{\mathrm{}}}\right){\displaystyle \frac{1}{\mathrm{\Gamma }^2+(E_nE_m)^2}}`$ (10)
where $`\mathrm{\Gamma }`$ is determined by normalization. If we make the replacement $`\mathrm{\Gamma }0`$ we get the standard result of first order perturbation theory (see Eq.(45)). The presence of $`\mathrm{\Gamma }`$ in the denominator reflects corrections to infinite order. We note that $`P_t^{\text{prt}}(n|m)`$ constitutes a generalization of Wigner Lorentzian. We indeed would get from it a Lorentzian if the bandprofile were flat with no finite bandwidth.

Whenever $`P_t(n|m)P_t^{\text{prt}}(n|m)`$ we say that the spreading profile is ”perturbative”. The perturbative structure is characterized by having separation of scales:

$`\mathrm{\Gamma }(t)\delta E(t)`$ (11)
We say that $`P_t(n|m)`$ has a ”standard” perturbative structure Eq.(45), which is given by first order perturbation theory, if $`\mathrm{\Gamma }<\mathrm{\Delta }`$. This means that more than $`50\%`$ probability is concentrated in the initial level. We use the term ”core-tail” structure if we want to emphasize the existence of a finite non-perturbative ”core” region $`|r|<\mathrm{\Gamma }/\mathrm{\Delta }`$. The core width $`\mathrm{\Gamma }`$, as determined by ”perturbation theory” by imposing normalization on Eq.(VIII), constitutes a rough estimate for the width $`\mathrm{\Gamma }(t)`$ of the energy distribution (Eq.(7)).

The spreading $`\delta E(t)`$ for the core-tail structure is determined by the tail region, which is defined as $`|r|\mathrm{\Gamma }/\mathrm{\Delta }`$. If Eq.(VIII) were a Lorentzian, it would imply $`\delta E(t)=\mathrm{}`$. This is of course not the case, because the spectral functions provide a physical cutoff. Thus, the core-tail structure which is described by Eq.(VIII) is characterized by a ”tail” component that contains a vanishingly small probability but still dominates the variance.

If we do not have the separation of energy scales Eq.(11), then perturbation theory becomes useless vrn ; frc . In such case $`P_t(n|m)`$ becomes purely non-perturbative. In order to determine $`P_t(n|m)`$ we have to use tools that go beyond perturbation theory. In particular, for long times one can justify vrn ; frc  a stochastic approximation, leading to the Gaussian kernel

$`P_t^{\text{sto}}(n|m)={\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \frac{\mathrm{\Delta }}{\delta E(t)}}\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{E_nE_m}{\delta E(t)}}\right)^2\right]`$ (12)
Note that for very long times this Gaussian should be replaced by an appropriate solution of a diffusion equation. But this is not relevant to our simulations. For a critical discussion that incorporates estimates for the relevant time scales see vrn ; frc .

A final remark: For a Gaussian profile the $`50\%`$ width is $`\mathrm{\Gamma }(t)=1.35\times \delta E(t)`$. However, whenever a non-perturbative structure is concerned it is better, in order to avoid confusion, not to use the notation $`\mathrm{\Gamma }(t)`$. The notation $`\mathrm{\Gamma }`$ has been adopted in the common diagrammatic formulation of perturbation theory. In case of non-perturbative structure this formulation becomes useless, and therefore the significance of $`\mathrm{\Gamma }`$, as a distinct energy scale, is lost.

## IX Analysis of the Parametric evolution, $`\delta x`$ regimes

The parametric evolution of $`P(n|m)`$ is illustrated in Fig.1 (upper panel). For very small $`\delta x`$ we observe clearly a standard perturbative profile (Eq.(45)), whose structure is just a reflection of the bandprofile. By definition this holds in the standard perturbative (parametric) regime $`\delta x<\delta x_c`$. From the numerics (looking on the lower panel of either Fig.3 or Fig.6) we find that $`\delta x_c=0.006`$.

For $`\delta x>\delta x_c`$ the initial level starts to mix with neighboring levels. As a result a non-perturbative ‘core‘ component starts to develop. Thus we obtain a core-tail structure. The tail is the perturbative component. The main component of the tail is the ”first order tail” which can be calculated using first order perturbation theory. The tails in the vicinity of the core are growing slower, which can be regarded as a suppression of core-to-tail transitions due to the mixing. There are also higher order tails lds  that can be neglected.

Having $`\delta E\mathrm{\Gamma }`$ is an indication that $`\delta E`$ is dominated by the (perturbative) tail component of the core-tail structure Eq.(VIII). The condition $`\mathrm{\Gamma }\delta E`$, can be re-written as $`\delta x\delta x_{\text{prt}}`$, which constitutes a definition of the extended perturbative (parametric) regime. Unfortunately, the strong inequality $`\mathrm{\Gamma }\delta E`$ is nowhere satisfied in our numerics. Therefore the numerical definition of $`\delta x_{\text{prt}}`$ becomes ambiguous. One may naively define $`\delta x_{\text{prt}}`$ by transforming the weak inequality $`\mathrm{\Gamma }<\delta E`$ into a weak inequality $`\delta x<\delta x_{\text{prt}}`$. But this procedure is numerically meaningless: The quantities $`\mathrm{\Gamma }`$ and $`\delta E`$ are energy scales. As such their definition is arbitrary up to a prefactor of order unity. After some thinking one realizes that the only practical definition for $`\delta x_{\text{prt}}`$ is as the $`\delta x`$ where the $`P_t^{\text{prt}}`$ based prediction of $`\delta E`$ becomes significantly less than the actual spreading. From the numerics (Fig.6) we find that $`\delta x_{\text{prt}}=0.05`$. This definition is not ambiguous numerically because we compare ”variance based measure” ($`\delta E`$) to ”variance based measure” ($`\delta E^{\text{prt}}`$), rather than comparing ”variance based measure” ($`\delta E`$) to a ”width measure” ($`\mathrm{\Gamma }`$).

For $`\delta x\delta x_{\text{prt}}`$ the $`P_t^{\text{prt}}`$ based calculation of $`\delta E`$ gives saturation. This (non-physical) saturation is the consequence of having finite bandwidth. It should be regarded as an artifact that reflects the limited validity of perturbation theory.

In the non-perturbative (parametric) regime, namely for $`\delta x>\delta x_{\text{prt}}`$, the tails are no longer the dominant component. This is associated with a structural change in the spreading profile. Looking at the upper panel of Fig.1 we observe a core-tail structure up to the 13th time step. Then, the secondary lobe of the bandprofile is swallowed by the core. This happening is reflected in Fig.6 (upper panel) by a boost in the width $`\mathrm{\Gamma }`$. \[The notion of ”secondary lobe” should be clear by looking at the $`|r|20`$ region in the upper panel of Fig.2. Unfortunately its visibility in the corresponding image (upper panel of Fig.1) is quite poor. Still it is possible to follow the evolution of the $`r20`$ tail region, and to realize that it is swallowed by the core\].

To have a weak inequality $`\delta x>\delta x_{\text{prt}}`$ is not quite the same as to be in the (deep) non perturbative regime where $`\delta x\delta x_{\text{prt}}`$. For $`\delta x>\delta x_{\text{prt}}`$ we do not get a purely non-perturbative structure. Inspite of the failure of the core-tail picture Eq.(VIII), we still can make a phenomenological distinction between ”core” and ”tail” regions. One clearly observes that the tail region is ”pushed” outside because of the expanding core. This non-perturbative effect is not captured by Eq.(VIII). In a previous study prm  the crossover from a core-tail structure to a purely non-perturbative (semiclassical) structure was quite abrupt. In the present study the ”deep” non-perturbative regime, where all the tail components disappear, is not accessible due to numerical limitations. A ”purely non perturbative” structure would be obtained if all the tail components were swallowed by the expanding core.

## X Analysis of actual evolution, $`V`$ regimes

We start this section with a qualitative description of the actual evolution, which is based on looking in illustrations such as in Fig.1 and Fig.2. Later we present the quantitative analysis. In general one observes that for finite $`V`$ the evolution acquires stochastic-like features. At intermediate times ($`0<t<T`$) the spreading profile contains both parametric and stochastic ”components”. The parametric component is ”reversible”. In contrast to that the stochastic component of the spreading is not affected by the velocity reversal, and may lead to a final Gaussian distribution.

Let us describe what happens to the evolution, as a function of $`\delta x=Vt`$, as we make simulations with smaller and smaller $`V`$. On the basis of Eq.(VIII) we expect to observe a modulation of the tails of $`P_t(r)`$ in a way that is implied by the presence of envelope $`\stackrel{~}{F}_t(\omega )`$. This modulation is characterized by the energy scale $`\mathrm{}/t=\mathrm{}V/\delta x`$. If we make $`V`$ further and further smaller, the secondary lobes of the tails \[see previous section for definition\] do not have a chance to become visible, because they are suppressed by the narrower envelope $`\stackrel{~}{F}_t(\omega )`$. For small enough $`V`$ only the main core component of the bandprofile is left visible. For very small $`V`$ the spreading profile becomes very close to a Gaussian shape in a very early stage of the evolution.

The spreading at the end of the pulse period is very small both is the sudden limit (very large $`V`$, corresponding to multi-period driving with large frequency), and also in the adiabatic limit (very small $`V`$). Fig.2 illustrates representative spreading profiles, which are observed in the end of the pulse period. Also displayed are the first-order perturbative profile (Eq.(45)), and a Gaussian with the same width. Note that the first-order perturbative profile Eq.(45), unlike Eq.(VIII), does not have a proper normalization. We can define a ‘difference measure’ in order to quantify the deviation of the (actual) lineshape from Gaussian shape. One possible definition is

$`\sqrt{{\displaystyle \underset{r}{}}P^{\text{sto}}(r)[\mathrm{log}(P_t(r))\mathrm{log}(P_t^{\text{sto}}(r))]^2}`$ (13)
The lower panel of Fig.8 shows that for $`V<20`$ the spreading profile at the end of the pulse is close to Gaussian shape, while for $`V>20`$ this profile is predominantly of perturbative nature.

We should address now the issue of $`V`$ regimes. Namely, for a given $`V`$ we would like to characterize the nature of the dynamical scenario. Below we shall introduce a practical procedure for the identification of the various regimes. We are going to explain that in the numerical simulation we observe four different $`V`$ regimes:     Adiabatic regime ($`V<3`$).    Perturbative regime ($`3<V<7`$).    Non-perturbative regime ($`7<V<20`$).    Sudden regime ($`20<V`$).    Before we get into details, we would like to make a connection with the theoretical discussion of regimes in rsp . There we have considered a (sinusoidal) periodic driving that is characterized by amplitude $`A`$ and frequency $`\mathrm{\Omega }=2\pi /T`$. For sinusoidal driving the root-mean square ”velocity” is $`V=A\mathrm{\Omega }/\sqrt{2}`$. So fixing $`A`$ and changing $`V`$ in the present paper, is completely analogous to fixing $`A`$ and changing $`\mathrm{\Omega }`$ in Ref.rsp . Thus, the four $`V`$ regimes listed above correspond to a horizontal cut ($`A=\text{const}>A_{\text{prt}}`$) in the $`(\mathrm{\Omega },A)`$ regime diagram of rsp  (see Fig. 1 there). One should realize that the ”perturbative regime” is in fact the ”linear response (Kubo) regime” as defined in rsp , and that the ”sudden regime” corresponds to the regime of vanishing (high frequency) response.

In the adiabatic $`V`$ regime the spreading is due to near neighbor level transitions, and disregarding extremely short times (for which we can apply standard first order perturbation theory (Eq.(45)), it looks stochastic. This means that Eq.(12) is a quite satisfying approximation. For further details regarding the identification of the adiabatic regime see Section 12.

Outside of the adiabatic regime there is a time stage where perturbation theory (Eq.(VIII)) is valid. But after this time we have to use theoretical considerations that go beyond perturbation theory. In the following paragraph we divide the non-adiabatic regime into ”perturbative” $`V`$ regime, ”non-perturbative” $`V`$ regime, and ”sudden” regime. The basis for this distinction is the timing of the departure from parametric evolution. The departure from parametric evolution is best illustrated by the lower panel of Fig.4.

The perturbative $`V`$ regime is defined by the requirement of having the departure from parametric evolution happen before the breakdown of Eq.(VIII). Consequently, in this regime the departure time can be deduced from Eq.(VIII). For example, let us assume a simple bandprofile. In such case departure time is just the inverse of the bandwidth, which implies via Eq.(30) that it is simply the classical correlation time $`\tau _{cl}`$ of the chaotic motion. In the actual numerical analysis the bandprofile is not ”simple”, but has some structure. Rather than debating over the definition of $`\tau _{cl}`$, it is more practical to determine the departure time by inspection of Fig.4. For convenience we have indicated the location of $`\delta x_{\text{prt}}`$ by a vertical line. For $`V<7`$ the departure from parametric evolution happens before this line. This way one can determine the upper border ($`V=7`$) of the perturbative regime.

The non-perturbative $`V`$ regime is defined by the requirement of having the departure from parametric evolution happen after the breakdown of Eq.(VIII). In other words, it means that this departure cannot be captured by perturbation theory. A more careful definition of the non-perturbative regime should take into account the driving reversal at $`t=T/2`$. As explained in the next paragraph, the non-perturbative regime is further restricted by the requirement of having the departure from parametric evolution happens before the driving reversal.

The distinction between ”perturbative” and ”non-perturbative” $`V`$ regimes is related to the possibility to ”capture” the stochastic features of the energy spreading by perturbation theory. Non-perturbative features of the energy distribution during intermediate times are of no relevance if they are of parametric (non-stochastic) origin. Let us look for example, in the upper panel of Fig.1. We clearly see that the spreading profile at the end of the pulse (at $`t=T`$) is of standard perturbative nature. All the non-perturbative features that develop during the first half of the driving cycle, are completely reversed in the second half of the cycle.

Conceptually the simplest way to set a criterion for getting at the end of the pulse a perturbative structure, is to use ”fixed basis” perturbation theory (see Appendix C of frc ). If we adhere to the present (more physical) approach of using $`x`$-dependent basis, an equivalent method frc  is to determine the sudden time $`t_{\text{sdn}}`$. This is the time to resolve the expanding ”core”. It is determined as the time when the following inequality is violated:

$`\mathrm{\Gamma }(t)\mathrm{}/t`$ (14)
The regime where the condition $`t_{\text{sdn}}T`$ is violated is defined as the ”sudden regime”. In the sudden regime all the non-perturbative features are of parametric nature, and therefore the validity of perturbation theory survives at the end of the pulse. In the sudden regime the departure from parametric evolution becomes visible only after the driving reversal. One can determine the sudden regime simply by looking at the lower panel of Fig.8. The departure of the energy distribution (at $`t=T`$) from Gaussian shape is correlated with getting into the sudden regime.

The discussion above of $`V`$ regimes was based on looking for the breaktime of perturbation theory Eq.(VIII), on the one hand, and looking for the departure from parametric behavior on the other hand. The departure from parametric behavior is an indication for the appearance of a predominant stochastic component in the spreading profile. This departure does not imply that we get a Gaussian line shape. In the perturbative regime we get the Gaussian line shape after the breaktime of perturbation theory, which happens after the departure from parametric behavior. Therefore in the perturbative regime there are 3 stages in the evolution: a parametric stage, a perturbative stochastic stage, and a genuine stochastic stage. In contrast to that, in the non-perturbative regime we do not have an intermediate ”perturbative stochastic stage”, because the departure form parametric behavior happens after the breakdown of perturbation theory.

## XI LRT formula

The general LRT formula for the variance of the spreading is

$`\delta E(t)^2=A^2\times {\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi }}\stackrel{~}{F}_t(\omega )\stackrel{~}{C}(\omega )`$ (15)
The proportionality $`\delta EA`$ reflects having ”linear response”. Two spectral functions are involved: One is the the spectral content of the driving (Appendix A), and the other is the power spectrum of the fluctuations (Appendix B). The latter is the Fourier transform of a correlation function $`C(\tau )`$.

A special case is the sudden limit ($`V=\mathrm{}`$) for which $`\stackrel{~}{F}_t(\omega )=1`$ and accordingly

$`\delta E(t)=\sqrt{C(\tau =0)}\times A`$ (16)
Another special case is the response for persistent (either linear or periodic) driving where $`\stackrel{~}{F}_t(\omega )=t\times 2\pi \delta (\omega )`$ implies diffusive behavior:

$`\delta E=\sqrt{2D_Et}`$ (17)
In such case the expression for $`D_E`$ is known as Kubo formula, leading to a fluctuation-dissipation relation.

Finding the conditions for validity of LRT is of major importance. This should be regarded as the first step in the analysis of the response of a driven system. The classical derivation of Eq.(15) is quite simple, and for completeness we present it in Appendix C. For sake of the following discussion one can assume that the classical ”slowness” conditions for the validity of classical LRT are satisfied.

The quantum mechanical derivation of this formula is much more subtle crs ; frc ; vrn ; rsp  and leads to the distinction between ”adiabatic” and ”(extended) perturbative” and ”non-perturbative” regimes. The semiclassical limit is contained in the latter regime. It is known that the LRT formula does not hold in the adiabatic regime wilk , but it is valid in the (extended) perturbative regime crs ; frc . It is not necessarily valid in the non-perturbative regime crs ; rsp , but if the system has a classical limit, it must be valid again in the semiclassical regime. See further discussion in the next sections.

## XII Analysis: Comparison with LRT

Using the bandprofile as an input we can calculate the spreading using Eq.(15). In Fig.6 the calculation is done for the parametric evolution (dashed line given by Eq.(16)), while in Fig.7 the calculation is done with Eq.(15) for finite $`V`$ values. The latter figure should be compared with the upper panel of Fig.4. The agreement is very good unless the velocity $`V`$ is small. See later discussion of the QM-adiabatic regime. The spreading $`\delta E`$ at the end of the pulse is better illustrated in Fig.8.

There are three possible strategies for evaluating the bandprofile. The first is to use the semiclassical strategy with Eq.(30). The second is to make a careful numerical evaluation of the matrix elements, and to use the definition Eq.(28). This gives the thin line in Fig. 5. However, the most practical and economical procedure is simply to deduce the bandprofile from the spreading profile $`P(r)`$ that correspond to the smallest $`\delta x`$ value. This is the thick line in Fig. 5, which we regard as the most appropriate estimate.

We found out that the LRT formula Eq.(15) is not sensitive to the way in which the bandprofile is evaluated, unless the $`|nm|=1`$ elements of the $`𝐅_{nm}`$ matrix dominate the result. Let us denote by $`\sigma `$ the root-mean-square magnitude of these matrix elements. The sensitivity to $`\sigma `$ happens in the $`V`$ regime which is determined by the condition $`V<(\mathrm{\Delta })^2/(\mathrm{}\sigma )`$, which is the adiabaticity condition. This sensitivity can be used as a practical tool for the determination of the adiabatic regime. In the upper panel of Fig. 8 we display the result of the LRT calculation upon setting $`\sigma =0`$. We see that the LRT calculation implies QM-adiabatic behavior for $`V<3`$.

If $`V`$ is well away from the adiabatic regime, then near-neighbor level transitions have negligible contribution. In such case the LRT calculation is not sensitive to the exact value of $`\sigma `$. As $`V`$ becomes smaller, there is a larger relative weight to the near-neighbor matrix elements.

Deep in the QM-adiabatic regime the mean-level energy difference is resolved much before the levels are mixed, and therefore, as a result of recurrences, the probability stays concentrated in the initial level. This is of course a leading order description. In fact we cannot neglect higher order corrections.

In the adiabatic regime the spreading is dominated by transition between neighboring levels (whereas outside of the adiabatic regime the contribution of neighboring level transitions can be neglected). Therefore, it is only in the adiabatic regime where the quantum-mechanical calculation (using Eq.(15)) gives results that are different from the classical expectation ophir . Is it possible to witness a regime where the quantum mechanical (rather than the classical) LRT prediction is observed? Apparently the answer is positive, but not in a typical numerical experiment. The reason is that almost always the Landau-Zener (non LRT) mechanism for energy spreading takes over wilk .

We can verify the dominance of Landau-Zener mechanism as follows. The theoretical prediction is $`\delta E(t)\sqrt{D_Et}`$ with $`D_EV^{3/2}`$. In the upper panel of Fig. 8 we plot the spreading at the end of half pulse period ($`t=T/2=A/V`$) and at the end of full pulse period ($`t=T=2A/V`$). Hence we expect $`\delta EV^{1/4}`$. The agreement with this expectation is quite good if we consider the spreading after half pulse period. At the end of full pulse period we get that the spreading is larger than expected. This is apparently due to the non-adiabatic nature of the $`VV`$ switching.

## XIII Discussion and Conclusions

The numerical study that we have presented should be regarded as an application of a general procedure for the analysis of energy spreading. In the summary below we re-order the stages in this analysis in the way which is implied by the results of the previous sections.

The first step in the analysis is to find the band profile. This can be done using the semiclassical recipe (Appendix B) without any need to make heavy numerical simulations. Then it is possible to determine $`\delta E(t)`$ by using the LRT formula Eq.(15).

The second step in the analysis is to determine the adiabatic $`V`$ regime. This is done by checking whether the LRT calculation of $`\delta E(t)`$ is sensitive to $`\sigma `$. In order to get in the adiabatic regime LRT-based quantum corrections to the classical result, we have to take the level spacing statistics into account ophir ; robbins . We have pointed out the difficulty in observing such corrections. Rather we can improve over the LRT prediction by taking into account either higher order or Landau-Zener corrections to perturbation theory.

The third step in the analysis is to calculate the parametric perturbative profile $`P^{\text{prt}}(r)`$. This is associated with getting a rough estimate for the core width $`\mathrm{\Gamma }`$. Also the parametric scales $`\delta x_c`$ and $`\delta x_{\text{prt}}`$ can be determined on the basis of this analysis. The former is determined from inspecting $`\mathrm{\Gamma }`$, while the latter is determined by comparing $`\delta E`$ of the LRT calculation to the $`P^{\text{prt}}`$ based prediction.

The fourth step in the analysis is to distinguish between the perturbative and the non-perturbative regime. For this purpose we have to look for the timing of the departure from parametric behavior, as in the lower panel of Fig.4. For $`\delta E`$ one can use the LRT based calculation of Eq.(15).

The fifth step in the analysis is to identify the sudden regime. For this purpose we have to eliminate the ”sudden time” from Eq.(14). The problem is to estimate $`\mathrm{\Gamma }(t)`$ in the non-perturbative regime. If we have only the bandprofile as an input, then we can use the rough estimate $`\mathrm{\Gamma }(t)\delta E(t)`$. This is based on the assumption that in the non-perturbative regime the energy distribution is characterized by a single energy scale.

The sixth step in the analysis is to make a simulation of the parametric evolution. This is a relatively heavy task, but it is still much easier than making finite $`V`$ simulation. \[It requires merely diagonalizations, while a temporal simulation requires an iterative procedure\]. The main non-trivial effects that we have found in the analysis of the parametric evolution were:     Higher order tails grow up.    A non-perturbative core region develops.    The core-to-tail transitions are suppressed.    The tails are ”pushed out”.    Tail components are swallowed by the expanding core.    The last two items are in the spirit of the ”core-tail” theory, but go beyond the perturbative approximation Eq.(VIII). Knowledge of the parametric evolution allows accurate determination of $`\mathrm{\Gamma }`$, leading to a refined determination of the $`V`$ regimes.

The seventh step in the analysis is to make simulation of the actual evolution. We would like to emphasize that this is the first time that the different dynamical scenarios (corresponding to the different $`V`$ regimes) are illustrated in a numerical simulation. The only other numerical studies (mainly diego ; rsp ) were too restricted in scope, and did not contain analysis of the stages in the evolution of the energy distribution.

The only non-perturbative effect on the response that we have discussed so far is the Landau-Zener correction to the spreading in the adiabatic regime. Are there any other non-perturbative effects that affect the response? The immediate tendency is to regard LRT as the outcome of standard first order perturbation theory (Appendix D). Then the question that arises is what happens if Eq.(45) does not apply?

Let us recall the answer in case of the parametric evolution of $`P(n|m)`$. As $`\delta x`$ becomes larger we should be worried regarding the implications of having the effects that are listed at the end of Section 9. Having expanding ”core” that ”pushes out” the tails, and having growing higher-order tail components, may suggest that the first order calculation of the variance (Eq.(D)) should be ”corrected”, and should include ”higher order” terms. Does it mean that Eq.(16) underestimate the spreading? Or maybe we should assume that Eq.(VIII) provides the correct ”trend” of higher order corrections? The tails in the vicinity of the core are growing slower, which can be regarded as a suppression of core-to-tail transitions due to the mixing. Consequently we would conclude that Eq.(16) is an overestimation of spreading. Moreover, if we take Eq.(VIII) too seriously, beyond its regime of validity, we would conclude that the spreading has saturation for large $`\delta x`$.

For the parametric evolution these possible speculations turn out to be wrong. The above effects are exactly balanced, and the LRT formula Eq.(16) remains exact beyond any order of perturbation theory, which means that it is exact even in the non-perturbative regime where perturbation theory is not applicable. This claim has a simple derivation lds . Note also that if the system has a classical limit, then the validity of Eq.(16) can be established deep in the non-perturbative regime, where the semiclassical approximation becomes reliable.

The question is whether this delicate balance is violated in case of finite $`V`$. For example, it may be that due to incomplete core-tail recurrences we shall have enhanced spreading (compared with LRT). Unlike the parametric case, we do not have a theoretical proof that excludes such a possibility. In fact the contrary is the case. We have demonstrated that for an artificial (random matrix theory) model, the LRT formula cannot be trusted in the non-perturbative regime. Whether such effect is possible also for ”quantized” system that possess a good classical limit has been left an open question.

The possibility of having deviation from LRT in the non-perturbative regime was one important motivation for the present research. Clearly we did not witness such effect in our simulations (Fig.8). This reflects that there is a clash between the semiclassical limit and the RMT limit.

## Appendix A  <br>The spectral function $`\stackrel{~}{F}_t(\omega )`$

Given a function $`f(t^{})`$ that describes the shape of the driving pulse during the time interval $`0<t^{}<t`$, we define a spectral function $`\stackrel{~}{F}_t(\omega )`$ by Eq.(3). This spectral function describes the spectral content of the driving pulse. Below we list some useful driving schemes. We use the notation $`\mathrm{\Theta }()`$, which is defined by $`\mathrm{\Theta }(\text{False})=0`$ and $`\mathrm{\Theta }(\text{True})=1`$.

For step function we have

$`f(t^{})`$ $`=`$ $`\mathrm{\Theta }(0<t^{})`$ (18)
$`\stackrel{~}{F}_t(\omega )`$ $`=`$ $`1`$ (19)
For rectangular pulse we have

$`f(t^{})`$ $`=`$ $`\mathrm{\Theta }(0<t^{}<t)`$ (20)
$`\stackrel{~}{F}_t(\omega )`$ $`=`$ $`|1\text{e}^{i\omega t}|^2`$ (21)
If it is followed by a negative pulse pulse we get

$`f(t^{})`$ $`=`$ $`\mathrm{\Theta }(0<t^{}<t)\mathrm{\Theta }({\displaystyle \frac{T}{2}}<t^{}<t)`$ (22)
$`\stackrel{~}{F}_t(\omega )`$ $`=`$ $`|12\text{e}^{i\omega \frac{T}{2}}+\text{e}^{i\omega t}|^2\text{for}t>{\displaystyle \frac{T}{2}}`$ (23)
For linear driving we have

$`f(t^{})`$ $`=`$ $`t^{}`$ (24)
$`\stackrel{~}{F}_t(\omega )`$ $`=`$ $`t^2(\text{sinc}(\frac{1}{2}\omega t))^2`$ (25)
where $`\text{sinc}()=\mathrm{sin}()/()`$. Note that for very large $`t`$ we have

$`\stackrel{~}{F}_t(\omega )t\times 2\pi \delta (\omega )`$ (26)
Finally, for triangular pulse

$`f(t^{})`$ $`=`$ $`2{\displaystyle \frac{t^{}}{T}}\mathrm{\Theta }\left(0<t^{}<{\displaystyle \frac{T}{2}}\right)+`$
$`2(1{\displaystyle \frac{t^{}}{T}})\mathrm{\Theta }\left({\displaystyle \frac{T}{2}}<t^{}<T\right)`$ (27)
we get

$`\stackrel{~}{F}_t(\omega )`$ $`=`$ $`\left[{\displaystyle \frac{2}{T}}\right]^2t^2(\text{sinc}(\frac{1}{2}\omega t))^2)\text{for}0<t<{\displaystyle \frac{T}{2}}`$
$`\stackrel{~}{F}_t(\omega )`$ $`=`$ $`\left[{\displaystyle \frac{2}{T}}\right]^2\left|{\displaystyle \frac{12\text{e}^{i\omega \frac{T}{2}}+\text{e}^{i\omega t}}{\omega }}\right|^2\text{for}{\displaystyle \frac{T}{2}}<t<T`$
## Appendix B  <br>The spectral function $`\stackrel{~}{C}(\omega )`$

Let us denote by $`𝐅_{mn}`$ the matrix representation of some quantized observable $`F(Q,P)`$, in the basis which is determined by some quantized chaotic Hamiltonian $`(Q,P)`$. The bandprofile of $`𝐅_{mn}`$ is conveniently characterized by the spectral function

$`\stackrel{~}{C}(\omega )={\displaystyle \underset{n(m)}{}}\left|𝐅_{nm}\right|^22\pi \delta \left(\omega {\displaystyle \frac{E_nE_m}{\mathrm{}}}\right)`$ (28)
with implicit average over the reference state $`m`$. This is the power spectrum of the fluctuating quantity $`(t)`$, whose classical definition is $`(t)=F(Q(t),P(t))`$, with a corresponding quantum-mechanical definition within the Heisenberg picture. The power spectrum of a fluctuating quantity $`(t)`$ is defined as the Fourier transform of the corresponding correlation function $`C(\tau )`$.

Chaotic systems are characterized by fast decay of dynamical correlations. We assume a separation of time scales between the (short) classical correlation time, and the (long) quantum-mechanical Heisenberg time. Thus, for times during which we can ignore the recurrences, we expect the following quantal-classical correspondence:

$`C(\tau )C^{\text{cl}}(\tau )`$ (29)
This correspondence implies that the envelope of $`\stackrel{~}{C}(\omega )`$ is given by $`\stackrel{~}{C}^{\text{cl}}(\omega )`$, and the following semiclassical expression for the matrix elements follows mario :

$`\left|𝐅_{nm}\right|^2{\displaystyle \frac{\mathrm{\Delta }}{2\pi \mathrm{}}}\stackrel{~}{C}^{\text{cl}}\left({\displaystyle \frac{E_nE_m}{\mathrm{}}}\right)`$ (30)
Taking into account the level spacing statistics, we deduce the following relation ophir :

$`\stackrel{~}{C}(\omega )\widehat{R}(\omega )\stackrel{~}{C}^{\text{cl}}(\omega )`$ (31)
where $`\widehat{R}(\omega )_n\delta (\omega \omega _{nm})_m`$ is the two-point correlation function of the energy spectrum (Fourier transform of the spectral form factor). For the Gaussian orthogonal ensemble (GOE) of random matrix theory (RMT) the following result is well known:

$`\widehat{R}(\omega )|_{\text{GOE}}=1(\text{sinc}(\pi \mathrm{}\omega /\mathrm{\Delta }))^2`$ (32)
and more generally it is common to assume $`\widehat{R}(\omega )\omega ^\beta `$ for small $`\omega `$.

## Appendix C Classical derivation of the LRT formula

Consider $`(Q,P;x(t))`$, and define the following time dependent quantities:

$`(t)`$ $`=`$ $`{\displaystyle \frac{}{x}}(Q(t),P(t);x(t))`$ (33)
$`(t)`$ $`=`$ $`(Q(t),P(t);x(t))`$ (34)
$`^{}(t)`$ $`=`$ $`(Q(t),P(t);x(0))`$ (35)
Note that $`(t)`$ is the energy in the conventional sense, while $`^{}(t)`$ is the energy using a ”fixed basis”. We have the following relations

$`{\displaystyle \frac{d(t)}{dt}}={\displaystyle \frac{}{t}}`$ $`=\dot{x}(t)(t)`$ (36)
$`{\displaystyle \frac{d^{}(t)}{dt}}=[,_0]`$ $`\delta x(t)\dot{}(t)`$ (37)
The latter approximated equality strictly holds if the perturbation $`_0`$ is linear with respect to the perturbation parameter $`\delta x=xx_0`$. From the equations above it follows that

$`(t)(0)`$ $`=`$ $`{\displaystyle _0^t}\dot{x}(t^{})(t^{})𝑑t^{}`$ (38)
$`^{}(t)^{}(0)`$ $``$ $`{\displaystyle _0^t}\delta x(t^{})\dot{}(t^{})𝑑t^{}`$ (39)
Obviously the two latter expressions coincide for a cycle ($`x(t)=x(0)`$). Squaring Eq.(38), and performing a microcanonical average over the (implicit) initial conditions $`(Q(0),P(0))`$ one obtains

$`\delta E(t)^2=A^2{\displaystyle _0^t}{\displaystyle _0^t}\dot{f}(t^{})\dot{f}(t^{\prime \prime })C(t^{}t^{\prime \prime })𝑑t^{}𝑑t^{\prime \prime }`$ (40)
This can be written as Eq.(15).

## Appendix D First order perturbation theory and LRT

For convenience we assume that the perturbation is linear in $`\delta x=xx_0`$. In complete analogy with the classical analysis we can work either using ”fixed basis”, or else we can use the proper $`x`$-dependent basis (the so called adiabatic representation). The respective matrix representations of the Hamiltonian are

$`𝐄+\delta x(t)𝐅`$ (41)
$`𝐄+\dot{x}(t)𝐖`$ (42)
where $`𝐄`$ is a diagonal matrix. Note that in the first equation $`𝐄`$ and $`𝐅`$ are calculated for $`x=x_0`$, while in the second equation there is an implicit $`x(t)`$ dependence. The matrix elements of $`𝐅`$ are

$`𝐅_{nm}=n\left|{\displaystyle \frac{}{x}}\right|m`$ (43)
The off diagonal matrix elements of $`𝐖`$ are

$`𝐖_{nm}={\displaystyle \frac{\mathrm{}}{i}}n|{\displaystyle \frac{d}{dt}}m=i{\displaystyle \frac{\mathrm{}}{E_nE_m}}𝐅_{nm}`$ (44)
and we use the ‘gauge’ convention $`𝐖_{nm}=0`$ for $`n=m`$. (Only one parameter is being changed and therefore Berry’s phase is not an issue). The derivation of Eq.(42) is standard, and can be found in Section 11 of frc .

Using first order perturbation theory with Eq.(42) we get

$`P_t(n|m)=\delta _{nm}+A^2\stackrel{~}{F}_t\left({\displaystyle \frac{E_nE_m}{\mathrm{}}}\right)\times \left|{\displaystyle \frac{𝐖_{nm}}{\mathrm{}}}\right|^2`$ (45)
Note that for a cycle ($`f(t)=f(0)=0`$), the same result is obtained via first order perturbation theory with Eq.(41). As a global approximation Eq.(45) is valid only in the standard perturbative regime. In the extended perturbative regime it is valid only for the first order tail region (see Section 8). The expression for the first order tail region can be written in a concise way as

$`P_t(n|m)=A^2\stackrel{~}{F}_t(\omega _{nm})\times {\displaystyle \frac{\mathrm{\Delta }}{2\pi \mathrm{}}}\stackrel{~}{C}(\omega _{nm})\left[{\displaystyle \frac{1}{\mathrm{}\omega _{nm}}}\right]^2`$ (46)
The corresponding global approximation is given by Eq.(VIII).

Assuming that the variance is dominated by the first order tail component, we get the LRT result

$`\delta E(t)^2`$ $`=`$ $`{\displaystyle \underset{n}{}}P_t(n|m)\times (E_nE_m)^2`$
$`=`$ $`A^2\times {\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{d\omega }{2\pi \mathrm{}}}\stackrel{~}{F}_t(\omega )\stackrel{~}{C}(\omega )`$ (47)
An implicit average over $`m`$ is assumed. Note that the latter formula does not contain $`\mathrm{}`$. This ”restricted” quantal-classical correspondence holds only for the variance. Higher moments of the energy distribution are typically much smaller compared with the classical expectation, and scale like $`\mathrm{}`$ to the power of the moment order minus 2.

Acknowledgments: It is our pleasure to thank Tsampikos Kottos and Eduardo Vergini for useful discussions. The research was supported in part by the Israel Science Foundation (grant No.11/02). DAW gratefully acknowledges support from CONICET (Argentina). Research grants by CONICET and ECOS-SeTCIP.

Fig.1. Images of time evolution within one-pulse period for $`V=100`$ (upper panel) and for $`V=1`$ (lower panel). Each column is a profile of the probability distribution $`P_t(r)`$ for a different time step $`t`$. The first 7 time steps are log spaced, while the rest are linearly spaced. The $`V=100`$ evolution is predominantly parametric.

Fig.2. Spreading profile at the end of one pulse period for $`V=100`$ (upper panel) and for $`V=10`$ (lower panel). The thin solid lines are from standard first order perturbation theory, which means $`P_t^{\text{prt}}`$ of Eq.(VIII) with $`\mathrm{\Gamma }=0`$. The dashed lines are for Gaussian with the same variance.

Fig.3. The width $`\mathrm{\Gamma }(t)`$ as a function of time. The numerical determination of $`\mathrm{\Gamma }`$ is explained after Eq.(7). The horizontal axis is the scaled time $`Vt`$. The different curves correspond to the different velocities $`V=100,80,40,20,10,7,5,3,1,0.5,0.2,0.1,0.07,0.05`$. The solid lines highlight the velocities $`V=100`$ (upper most), and $`V=10`$, and $`V=1`$, and $`V=0.1`$. The lower panel is the same as in the upper figure, but in log-log scale, and only half period is displayed.

Fig.4. The spreading $`\delta E(t)`$ of Eq.(8) as a function of time, for different $`V`$ simulations. The upper panel is in one to one correspondence with the upper panel of the previous figure. In the lower panel the data of the upper panel are presented in a different way: the vertical axis is the relative spreading; The horizontal axis is log scale; And only half period is displayed. The location of $`\delta x_{\text{prt}}=0.05`$ is indicated by vertical line. See the text for more details.

Fig.5. The band profile, as defined in Appendix B. The two curves are the outcome of two different numerical procedures: The thin line is based on direct evaluation of matrix elements, while the thick line is deduced form the evolution over an infinitesimal time step (see Section 12 for further details). A third possible procedure (not displayed) is to use the semiclassical recipe Eq.(30)

Fig.6. Parametric evolution: The horizontal axis is $`\delta x`$ either in linear scale (upper panel) or in logarithmic scale (lower panel). The dotted line is the width $`\mathrm{\Gamma }`$ from (7). The thick line is $`\delta E`$ from (8), while the dashed line is the LRT estimate Eq.(16). The thin line is the perturbative estimate Eq.(8) using Eq.(VIII). From this plots we can determine $`\delta x_c=0.006`$ and $`\delta x_{\text{prt}}=0.05`$.

Fig.7. LRT calculation of the spreading versus time, using Eq.(15) with the bandprofile as an input. The horizontal axis is $`Vt`$. This Figure is in one to one correspondence with the upper panel of Fig.4. One observes that for small $`V`$ the calculation underestimates the observed spreading. See discussion in Section 12.

Fig.8. upper panel: Response for one-pulse versus $`V`$. The Thick line is $`\delta E`$ at the end of the pulse. The thin (solid and dotted) lines correspond to the LRT calculation Eq.(15). For the calculation of the lower thin lines we have set the near-level coupling $`\sigma =0`$. The dashed line is the spreading at the end of half pulse period. The dash-dot line is the slope that corresponds to Landau-Zener spreading mechanism. The lower panel is the difference Eq.(13) from Gaussian line shape.


format: text

text:
# 1 Preliminaries

## 1 Preliminaries

This work is devoted to the study of duality symmetry in the context of the $`BF`$ theory, viz., in models presenting a topological, first-order derivative coupling between forms of different ranks that is a dimensional extension of the duality between the self-dual (SD)  and Maxwell-Chern-Simons models (MCS) , shown by Deser and Jackiw  long time ago. To this end we investigate the existence of a constraint of self duality in the massive, non invariant model ($`SD_{BF}`$) and adopt a new dynamical embedding formalism , that is alternative to the master Lagrangian approach, to obtain the gauge invariant $`BF`$ model. Our study also includes the case of dynamical fermionic matter minimally coupled to the self-dual sector.

Duality is a fascinating symmetry concept allowing the connection of two opposite regimes for the same dynamics. It plays an important role in nowadays physics, both in the original contexts of condensed matter and Maxwell electromagnetism, as well as in the recent research of extended objects. The existence of such a symmetry within a model has important consequences - it can be used to derive (exact) non perturbative results since swapping opposite regimes allows a perturbative investigation of theories with large coupling constants.

The study of this symmetry has received renewed interest in recent research in diverse areas in field theory such as, supersymmetric gauge theories , sine-Gordon model , statistical systems  and, in the context of condensed matter models, applied for instance to planar high-TC materials, Josephson junction arrays  and Quantum Hall Effect . In particular the duality mapping has been of great significance in order to extend the bosonization program from two to three dimensions with important phenomenological consequences . It also plays preponderant role in the ADS/CFT correspondence  that illustrates the holographic principle .

The idea of duality has also been used in recent developments of string theory , where different vacua are shown to be related by duality. In this context a general procedure for constructing dual models was proposed by Busher  and generalized by Rocek and Verlind  that consists in lifting the global symmetry of the tensor fields with a new gauge field, whose field strength is then constrained to zero by the use of a Lagrange multiplier. Integrating, sequentially, the multiplier and the gauge field yields the original action while the dual action is obtained if one integrates the gauge field together with the original tensor field, keeping the Lagrange multiplier that then plays the role of dual field to the original tensor field. This line of research was used in the investigation of bosonization as duality by Burgess and Quevedo  and to discuss S-duality, the relation between strong and weak couplings in gauge theories . This procedure has also been shown to be related to canonical transformations . Recently, this line of research has been applied in the context of the topologically massive $`BF`$ theory, which is related to our interest here, to study its equivalence with the Stuckelberg construction of gauge invariant massive excitations .

The duality we are treating in this work deals with the equivalence between models describing the same physical phenomenon involving the presence of a topological term in a four dimensional spacetime. It is closely related to the odd-dimensional duality involving the Chern-Simons term (CST) , whose paradigm is the equivalence between SD  and MCS  theories in (2+1) dimensions. As shown in , in three dimensions there are two different ways to describe the dynamics of a single, freely propagating spin one massive mode, using either the SD theory  or the MCS theory. They also established the identification that relates the basic field of the SD model with the dual of the MCS field . This correspondence displays the way the gauge symmetry of the MCS representation, gets hidden in the SD representation . It is well understood by now that it is the presence of the topological and gauge invariant Chern-Simons term the responsible for the essential features manifested by the three dimensional field theories, while in the four dimensional context this role is played by the $`BF`$ term. To extend this duality symmetry relation and study its consequences in the context of four dimensional field theories with $`BF`$ term, is the main purpose of this work.

The study of gauge theories with a topological term, in (3+1) dimensions, has received considerable attention recently. Among other possibilities the $`BF`$ term is interesting for providing a gauge invariant mechanism to give mass to the gauge field and to produce statistical transmutation in (3+1) dimensions. Here $`B`$ is a Kalb-Ramond field, i.e., a totally antisymmetric tensor potential (a potential 2-form) while $`F=dA`$ is the field strength of the one-form potential $`A`$. An abelian antisymmetric tensor potential was probably first used in the context of the particle theory to describe a massless particle of zero-helicity . It reappeared later on in the context of fundamental strings , has been used to study cosmic strings  and to put topological charge (hair) on black holes. The free theory of a rank-2 antisymmetric tensor has also been intensively studied both classically  and quantically  and has been shown to be dynamically dual (under the Hodge mapping) to a massless scalar field (zero-form).

In this work we study duality in field theories involving the $`BF`$ term. To this end, in the next section, we investigate the gauge non invariant $`SD_{BF}`$ model, define a new, non-Hodge, (derivative) duality operation and show the existence of self-duality. Next, in Section III, we apply an iterative dynamical embedding procedure to construct an invariant theory out of the self-dual $`BF`$ model - the topologically massive $`BF`$ model ($`TM_{BF}`$). This is a gauge embedding procedure that is done with the inclusion of counter terms in the non invariant action, built with powers of the Euler vectors and tensors (whose kernels give the field equations for the potentials $`A`$ and $`B`$) to warrant the dynamical equivalence. Such construction discloses hidden gauge symmetries in such systems. One can then consider the non-invariant model as the gauge fixed version of a gauge theory. A deeper and more illuminating interpretation of these systems is then obtained. The advantage in having a gauge theory lies in the fact that the underlying gauge symmetry allows us to establish a chain of equivalence among different models by choosing different gauge fixing conditions. In Section IV we consider the minimal coupling with fermionic matter. Our results are discussed in the final section of the paper.

## 2 Self Dual $`BF`$ Theory

The model with a built-in SD constraint in (2+1) dimensions was proposed in  as an alternative to the concept of topologically massive modes proposed in . The former is a non gauge invariant, first order model, while the later is a second order gauge invariant formulation, both making use of the topological Chern-Simons term. In this section we want to formulate and study a first order, non gauge invariant model, making use of the topological $`BF`$ term and prove the existence of self duality property as a consequence of a built in SD constraint.

The model in question shows the coupling of a vector field potential $`A_\mu `$ with a tensor field potential $`B_{\mu \nu }`$  as,

$$_{SD}^{(0)}=\frac{1}{2}m^2A_\mu A^\mu \frac{1}{4}B_{\mu \nu }B^{\mu \nu }+\frac{\chi \theta }{4}ϵ^{\mu \nu \lambda \rho }B_{\mu \nu }F_{\lambda \rho },$$
(1)
where the superscript index in the Lagrangean is the counter of the iterative algorithm to be implemented in the sequel, $`\chi =\pm 1`$ will be shown to display the self or anti-self duality, $`\theta `$ is the coupling constant and the field strength of the basic potentials are,

$`F_{\mu \nu }`$ $`=`$ $`_\mu A_\nu _\nu A_\mu `$
$`H_{\mu \nu \lambda }`$ $`=`$ $`_\mu B_{\nu \lambda }+_\nu B_{\lambda \mu }+_\lambda B_{\mu \nu }.`$ (2)
The coefficients of the mass terms are so chosen to give mass dimension one and two, respectively, to the potentials $`A_\mu `$ and $`B_{\mu \nu }`$ which, consequently, keeps dimensionless the coupling constant $`\theta `$ in the $`BF`$ term. Here the potentials play an active role in the duality transformations. This shall be in contrast with the dynamical matter field, to be considered latter on that, although coupled to the potentials, are passive fields (spectators) in the duality mapping. It is immediate to work out the equations of motion of the basic potentials $`A_\mu `$ and $`B_{\mu \nu }`$ to obtain, respectively,

$`A^\mu `$ $`=`$ $`{\displaystyle \frac{\chi \theta }{2m^2}}ϵ^{\mu \nu \lambda \rho }_\nu B_{\lambda \rho }`$
$`B^{\mu \nu }`$ $`=`$ $`\chi \theta ϵ^{\mu \nu \lambda \rho }_\lambda A_\rho `$ (3)
satisfying the constraints

$`_\mu A^\mu `$ $`=`$ $`0`$ (4)
$`_\mu B^{\mu \nu }`$ $`=`$ $`0`$ (5)
identically. The equations (2) constitute a set of first-order coupled equations that can be combined into a decoupled second-order, massive, wave equations as

$`\left(_\alpha ^\alpha +{\displaystyle \frac{m^2}{\theta ^2}}\right)F`$ $`=`$ $`0;F=\{A_\mu ;B_{\mu \nu }\}`$ (6)
whose mass depends crucially on the value of the coupling constant.

Next, we discuss the self-duality inherent to the above theory. To this end we define a new derivative duality operation by means of a set of star-variables as

$`^{}A_\mu `$ $``$ $`{\displaystyle \frac{\theta }{2m^2}}ϵ_{\mu \nu \lambda \rho }^\nu B^{\lambda \rho }`$
$`^{}B_{\mu \nu }`$ $``$ $`\theta ϵ_{\mu \nu \lambda \rho }^\lambda A^\rho `$ (7)
With this definition we obtain, for the double duality operation, the relations

$`^{}\left(^{}F\right)=F;F=\{A_\mu ;B_{\mu \nu }\}`$ (8)
after use of the equations of motion (6). This is important because it validates the notion of self (or anti-self) duality

$`^{}F=\chi F;F=\{A_\mu ;B_{\mu \nu }\}`$ (9)
as a solution for the field equations, very much like the three-dimensional SD model. However, this conceptualization of duality operation and self-duality in four-dimensions is new.

Before we start the iterative procedure for the transformation of the $`SD_{BF}`$ model into a topological $`BF`$ model let us digress on the consequences of the self-duality relation (9). Notice first that under the usual gauge transformations of the potentials

$`A_\mu `$ $``$ $`A_\mu +_\mu \mathrm{\Lambda }`$
$`B_{\mu \nu }`$ $``$ $`B_{\mu \nu }+_\mu \mathrm{\Lambda }_\nu _\nu \mathrm{\Lambda }_\mu `$ (10)
the fields strengths $`F_{\mu \nu }`$ and $`H_{\mu \nu \lambda }`$ are left invariant. Therefore, although the basic potentials are gauge dependent their duals, defined in (2), are not. This situation parallels the three-dimensional case involving the Chern-Simons term which is the origin for the presence of a hidden (gauge) symmetry in the SD model of  while it is explicit in the topologically massive model of . Here too the $`SD_{BF}`$ model hides the gauge symmetry (2) that is explicit in the $`TM_{BF}`$ model. Let us next consider the direct automorphism

$$(A,B)^{}(^{}A,^{}B)$$
(11)
Since it is constructed with the dual fields it is automatically gauge invariant and reads

$$^{}(^{}A,^{}B)=\frac{1}{12m^2}H_{\alpha \beta \lambda }H^{\alpha \beta \lambda }\frac{1}{8}F_{\alpha \beta }F^{\alpha \beta }+\frac{1}{2}ϵ^{\alpha \beta \lambda \rho }^\sigma F_{\sigma \rho }H_{\alpha \beta \lambda }$$
(12)
It is a simple algebra to check that the equations of motion for the SD action (1) are also solutions for the field equation of (12). This exercise clearly shows an intimate connection between the $`SD_{BF}`$ with a gauge invariant version through a dual transformation. However, although establishing the dual connection, the result obtained in (12) produces a set of field equations involving higher derivatives that will produce more solutions than the original set. Besides they are not the usual $`TM_{BF}`$ model. In the next section we shall discuss a dynamical gauge embedding procedure that will clearly produce an equivalent gauge invariant model.

## 3 The Gauge Invariant $`BF`$ Theory

In previous works  we have used the dynamical gauge embedding formalism to study dual equivalence in (2+1) dimensions in diverse situations with models involving the presence of the topological Chern-Simons term. In this section we extend that technique to study duality symmetry among four dimensional models involving the presence of a topological $`BF`$ term.

Our basic goal is to transform the symmetry (2) that is hidden in the Lagrangian (1) into a local gauge symmetry by lifting the global parameter $`\mathrm{\Lambda }`$ into its local form, i.e., $`\mathrm{\Lambda }\mathrm{\Lambda }(x^\mu )`$. The method works by looking for an (weakly) equivalent description of the original theory which may be obtained by adding a function $`f(K_\mu ,M_{\mu \nu })`$ to the Lagrangian (1). Here $`K_\mu `$ and $`M_{\mu \nu }`$ are the Euler tensors, defined by the variation

$`\delta _{SD}^{(0)}`$ $`=`$ $`K_\mu \delta A^\mu +M_{\mu \nu }\delta B^{\mu \nu }`$ (13)
whose kernels give the equations of motion for the $`A_\mu `$ and $`B_{\mu \nu }`$ fields, respectively. The minimal requirement for $`f(K_\mu ,M_{\mu \nu })`$ is that it must be chosen such that it vanishes on the space of solutions of (1), viz. $`f(0,0)=0`$, so that the effective Lagrangian $`_{eff}`$

$$_{SD}^{(0)}_{eff}=_{SD}^{(0)}+f(K_\mu ,M_{\mu \nu })$$
(14)
is dynamically equivalent to $`_{SD}^{(0)}`$. To find the specific form of this function that also induces a gauge symmetry into $`_{SD}^{(0)}`$ we work iteratively. To this end we compute the variations (13) of $`_{SD}^{(0)}`$ to find the Euler tensors as

$`K_\mu `$ $`=`$ $`m^2A_\mu {\displaystyle \frac{\chi \theta }{2}}ϵ_{\mu \nu \rho \lambda }^\nu B^{\rho \lambda }`$
$`M_{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{1}{2}}B_{\mu \nu }+{\displaystyle \frac{\chi \theta }{2}}ϵ_{\mu \nu \lambda \rho }^\lambda A^\rho `$ (15)
and define the first-iterated Lagrangian as,

$$_{SD}^{(1)}=_{SD}^{(0)}a_\mu K^\mu b_{\mu \nu }M^{\mu \nu }$$
(16)
with the Euler tensors being imposed as constraints and the new fields, $`a_\mu `$ and $`b_{\mu \nu }`$, to be identified with ancillary gauge fields, acting as a Lagrange multipliers.

The transformation properties of the auxiliary fields $`a_\mu `$ and $`b_{\mu \nu }`$ accompanying the basic field transformations (2) is chosen so as to cancel the variation of $`_{SD}^{(0)}`$, which gives

$`\delta a_\mu `$ $`=`$ $`\delta A_\mu `$
$`\delta b_{\mu \nu }`$ $`=`$ $`\delta B_{\mu \nu }`$ (17)
A simple algebra then shows

$`\delta _{SD}^{(1)}`$ $`=`$ $`a_\mu \delta K^\mu b_{\mu \nu }\delta M^{\mu \nu }`$ (18)
$`=`$ $`\delta \left({\displaystyle \frac{m^2}{2}}a_\mu a^\mu +{\displaystyle \frac{1}{4}}b_{\mu \nu }b^{\mu \nu }\right)`$
where we have used (2) and (3). Because of (18), the second iterated Lagrangian is unambiguously defined as

$$_{SD}^{(2)}=_{SD}^{(1)}+\frac{m^2}{2}a_\mu a^\mu \frac{1}{4}b_{\mu \nu }b^{\mu \nu }$$
(19)
that is automatically gauge invariant under the combined local transformation of the original set of fields ($`A_\mu `$, $`B_{\mu \nu }`$) and the auxiliary fields ($`a_\mu `$, $`b_{\mu \nu }`$).

We have therefore succeed in transforming the global $`SD_{BF}`$ theory into a locally invariant gauge theory. We may now take advantage of the Gaussian character displayed by the auxiliary field to rewrite (19) as an effective action depending only on the original variables ($`A_\mu `$, $`B_{\mu \nu }`$). To this end we use (19) to solve for the fields $`a_\mu `$ and $`b_{\mu \nu }`$ (call the solutions $`\overline{a}_\mu `$ and $`\overline{b}_{\mu \nu }`$ collectively by $`\overline{h}_{\{\mu \}}`$), and replace it back into (19) to find

$`_{eff}`$ $`=`$ $`_{SD}^{(2)}_{h_{\{\mu \}}=\overline{h}_{\{\mu \}}}`$ (20)
$`=`$ $`_{SD}^{(0)}{\displaystyle \frac{1}{2m^2}}K_\mu K^\mu +M_{\mu \nu }M^{\mu \nu }`$
from which we identify the function $`f(K_\mu ,M_{\mu \nu })`$ in (14). This dynamically modified action can be rewritten to give the $`TM_{BF}`$ theory,

$$_{eff}=\frac{1}{12m^2}H_{\mu \nu \lambda }H^{\mu \nu \lambda }\frac{1}{4}F_{\mu \nu }F^{\mu \nu }\frac{\chi }{2\theta }ϵ^{\mu \nu \alpha \beta }A_\mu _\nu B_{\alpha \beta }.$$
(21)
after the scaling $`\theta A_\mu A_\mu `$ and $`\theta B_{\mu \nu }B_{\mu \nu }`$ is performed. Notice the inversion of the coupling constant $`\theta 1/\theta `$ resulting from the duality mapping. It becomes clear from the above derivation that the difference between these two models is given by a function of the Euler tensors of the $`SD_{BF}`$ model that vanishes over its space of solutions. This establishes the dynamical equivalence between the $`SD_{BF}`$ and the $`TM_{BF}`$ theory.

## 4 The minimal coupling with fermionic matter

Once the duality mapping between the free theories has been established one is ready to consider the requirements for the existence of duality when the coupling with dynamical matter is included. In this section we consider the case of U(1) charged fermionic matter. For clarity, we consider first the situation where only the vector field, in the self-dual representation, is minimally coupled to the fermionic current. This seems appropriate since it illustrates the main features of the duality via gauge embedding. The results for the full coupling are then quickly presented in the following subsection.

### 4.1 Single Coupling

Let us consider first a Dirac field minimally coupled to the vector field $`A_\mu `$ specified by the $`SD_{BF}`$ model. To this end let us introduce the fermionic current

$`J^\mu =\overline{\psi }\gamma ^\mu \psi `$ (22)
and define the rank-2 charge $`J^{\mu \nu }`$ as the (derivative) dual of the fermionic current $`J_\mu `$,

$$J^{\mu \nu }=\frac{\theta }{2m^2}ϵ^{\mu \nu \alpha \beta }_\alpha J_\beta $$
(23)
so that the interacting Lagrangian becomes

$$_{min}^{(0)}=_{SD}^{(0)}eA_\mu J^\mu +_D,$$
(24)
where $`M`$ is the fermion mass. Here the Dirac Lagrangian is,

$$_D=\overline{\psi }(i/M)\psi .$$
(25)
The fermionic field is treated as an spectator in the dual transformation between the gauge fields from the SD to the TM sectors. But to remain as a bystander field the coupling to the gauge fields and to itself has to be readjusted in the TM, as shown below.

As before, our basic strategy is to transform the hidden symmetry of the Lagrangian (24) into a local gauge symmetry. A variation of the Lagrangian (24) gives the Euler vectors as,

$`K_\mu K_\mu ^D`$ $`=`$ $`K_\mu eJ_\mu `$
$`M_{\mu \nu }M_{\mu \nu }^D`$ $`=`$ $`M_{\mu \nu }.`$ (26)
¿From now on we follow the same steps as in the free case just making the replacement $`K_\mu K_\mu ^D`$ and $`M_{\mu \nu }M_{\mu \nu }^D`$ to obtain an effective action as,

$`_{eff}=_{min}^{(0)}{\displaystyle \frac{1}{2m^2}}K_\mu ^DK^{D\mu }M_{\mu \nu }^DM^{D\mu \nu }`$ (27)
A further manipulation shows the presence of a new term, compared to the free case and of a non-minimal coupling

$`_{eff}`$ $`=`$ $`{\displaystyle \frac{1}{12m^2}}H_{\mu \nu \lambda }H^{\mu \nu \lambda }{\displaystyle \frac{1}{4}}F_{\mu \nu }F^{\mu \nu }{\displaystyle \frac{\chi }{2\theta }}ϵ^{\mu \nu \alpha \beta }A_\mu _\nu B_{\alpha \beta }.`$ (28)
$`+`$ $`\overline{\psi }(i/M)\psi {\displaystyle \frac{e^2}{2m^2}}J_\mu J^\mu {\displaystyle \frac{\chi }{\theta }}eJ^{\mu \nu }B_{\mu \nu }`$
The presence of these two terms involving the matter current are worth discussing since they play important role in preserving the structure of the fermionic sector upon application of the dual mapping upon the gauge sector. They are the Thirring-like self-interaction and the non-minimal magnetic-like interaction of the dual of the fermionic current with $`B_{\mu \nu }`$. It is interesting to observe that the minimal coupling involving the $`A_\mu `$ field became, through the dualization process, a non-minimal coupling for the $`B_{\mu \nu }`$ field. This swapping of the coupling is hardly a surprise. It is the manifestation, in the latent sector of duality, of the traditional duality inversion and, as far as we know, this phenomenon has not been reported before. Notice that both terms appear naturally as a consequence of the embedding algorithm. The presence of these terms, as we show in the sequel, are important to maintain the dynamical structure of the fermionic matter in both representations of the dual pair . This can be seen by carefully examining the dynamics of the fermionic sector for both theories. To see this we compute the fermionic equation for the $`SD_{BF}`$ sector, that reads

$`\left(i/M\right)\psi =eA^\mu \gamma _\mu \psi .`$ (29)
To eliminate the bosonic field $`A^\mu `$ in favor of the fermionic one we rewrite the equations of motion for the $`SD_{BF}`$ sector as,

$`A^\mu `$ $`=`$ $`{\displaystyle \frac{\chi \theta }{2m^2}}ϵ^{\mu \nu \lambda \rho }_\nu B_{\lambda \rho }+eJ^\mu `$ (30)
$`B^{\mu \nu }`$ $`=`$ $`\chi \theta ϵ^{\mu \nu \lambda \rho }_\lambda A_\rho `$ (31)
which can be combined to give, after using the constraint (4),

$$R^1A_\mu =eJ_\mu $$
(32)
where $`R`$ is a differential operator such that its inverse is the wave-operator, defined as

$`R^1=_\lambda ^\lambda +{\displaystyle \frac{m^2}{\theta ^2}}.`$ (33)
Substituting this equation back in the matter equation gives,

$`\left(i/M\right)\psi =e^2RJ^\mu \gamma _\mu \psi ,`$ (34)
that is a nonlinear, integro-differential equation, now written completely in terms of the fermionic fields.

Next we compare this result with the equations of motion for the fermionic matter from the $`TM_{BF}`$ sector, that reads

$`\left(i/M\right)\psi =\chi e^{}A^\mu \gamma _\mu \psi +{\displaystyle \frac{e^2}{m^2}}J^\mu \gamma _\mu \psi `$ (35)
where we have to obtain the bosonic functional $`^{}A^\mu ^{}A^\mu \left(\psi \right)`$ in terms of the fermions from the equations of motion for the gauge fields in the $`TM_{BF}`$ sector,

$`^{}A^\mu `$ $`=`$ $`{\displaystyle \frac{\chi \theta }{2m^2}}ϵ^{\mu \nu \lambda \rho }_\nu ^{}B_{\lambda \rho }`$
$`^{}B^{\mu \nu }`$ $`=`$ $`\chi \theta ϵ^{\mu \nu \lambda \rho }_\lambda ^{}A_\rho 2eG^{\mu \nu }`$ (36)
that can be combined to give,

$`^{}A^\mu \left(\psi \right)=\chi e\left({\displaystyle \frac{\theta ^2}{m^2}}R\right)J^\mu .`$ (37)
When inserted into (35) gives

$`\left(i/M\right)\psi =e^2RJ^\mu \gamma _\mu \psi ,`$ (38)
which agrees with the dynamical equation for the fermions in the SD side. We just stress the importance of the Thirring like self-interaction for the fermions and the non-minimal interaction with the tensor field to keep the dynamics of the latent, fermionic sector unaltered.

### 4.2 Double Coupling

It is now a simple task to consider the full coupling of the bosonic fields $`A_\mu `$ and $`B_{\mu \nu }`$ to fermionic matter. To this end we introduce the rank-2 current

$`G^{\mu \nu }=𝒞\overline{\psi }\gamma ^\mu \gamma ^\nu \psi `$ (39)
where $`𝒞`$ is a complex normalization constant and its dual

$`G^\mu =\theta ϵ^{\mu \lambda \rho \sigma }_\lambda G_{\rho \sigma }`$ (40)
The interacting Lagrangian now takes the form

$$_{min}^{(0)}=_{SD}^{(0)}eA_\mu J^\mu +gB_{\mu \nu }G^{\mu \nu }+_D,$$
(41)
with $`e`$ and $`g`$ being the strengths of the coupling with $`A_\mu `$ and $`B_{\mu \nu }`$, respectively. The effective, gauge invariant action is obtained directly from (20) just operating the replacement

$`K_\mu K_\mu ^D`$ $`=`$ $`K_\mu eJ_\mu `$
$`M_{\mu \nu }M_{\mu \nu }^D`$ $`=`$ $`M_{\mu \nu }+gG_{\mu \nu }`$ (42)
to produce

$`_{eff}`$ $`=`$ $`_{min}^{(0)}{\displaystyle \frac{1}{2m^2}}K_\mu ^DK^{D\mu }M_{\mu \nu }^DM^{D\mu \nu }`$ (43)
which, after some algebraic manipulation, gives

$`_{eff}`$ $`=`$ $`_D+{\displaystyle \frac{1}{12m^2}}H_{\mu \nu \lambda }H^{\mu \nu \lambda }{\displaystyle \frac{1}{4}}F_{\mu \nu }F^{\mu \nu }{\displaystyle \frac{\chi }{2\theta }}ϵ^{\mu \nu \alpha \beta }A_\mu _\nu B_{\alpha \beta }`$ (44)
$``$ $`{\displaystyle \frac{e^2}{2m^2}}J_\mu J^\mu +g^2G^{\mu \nu }G_{\mu \nu }{\displaystyle \frac{\chi e}{\theta }}J^{\mu \nu }B_{\mu \nu }+{\displaystyle \frac{\chi g}{\theta }}A_\mu G^\mu `$
From this result it becomes clear the full action of the dual mapping over the active and passive fields involved in the transformation. Notice the exchange of the minimal coupling adopted in the SD sector into a non minimal, magnetic like interaction in the TM sector, including a swapping between the fields and currents and the presence of the current-current interaction for the fermionic sector which is needed to maintain the dynamics of the spectator field unmodified. this is easily checked by just computing the equations of motion for the Dirac fields in both representations to obtain,

$`\left(i/M\right)\psi `$ $`=`$ $`{\displaystyle \frac{R}{\theta ^2}}\gamma _\mu \left[\left(e^2J^\mu \chi geG^\mu \right)+2m^2𝒞\left(\chi geJ^{\mu \nu }+g^2G^{\mu \nu }\right)\gamma _\nu \right]\psi `$ (45)
## 5 Conclusions

In this work we studied dual equivalence in four dimensional topological models, namely, between the $`BF`$ self-dual ($`SD_{BF}`$) and the $`BF`$ topologically massive ($`TM_{BF}`$) models using an iterative procedure of gauge embedding that produces the dual mapping. We defined a new derivative type of duality mapping, very much like the one adoted in the three-dimensional case and proved the self and antiself-duality property of the $`SD_{BF}`$ model, according to the relative sign of the topological term. Working out the free case firstly, where the A and B fields participate actively in the dual transformation we obseved, as expected, the traditional inversion in the coupling constant. The coupling to dynamical fermionic matter, which acts as a spectator field in the dual transformation, brought into the scene some new features. The apearance of a Thirring like self-interaction term in the dualized theory, that had already been observed in the (2+1) case, as well as the shift from minimal to non minimal coupling. However, in this case we observed a swapping of the couplings from a tensor to another. This is a new result due to the presence of tensors of distinct ranks participating actively in the dual transformation. We proved that the presence of these terms are demanded to maintain the equivalent dynamics in the fermionic sector in either representations of the duality. The cases where the active tensors appear non linearly and the coupling with bosonic matter are postponed for a forthcoming publication.

ACKNOWLEDGMENTS: This work is partially supported by CNPq, CAPES, FAPERJ and FUJB, Brazilian Research Agencies. CW thanks the Physics Department of UFPB for the kind hospitality during the course of this investigation.