# Analytical model for attosecond time delays and Fano’s propensity rules in
the continuum

D.I.R. Boll<EMAIL_ADDRESS>Instituto de Física Rosario, CONICET-
UNR, Blvd. 27 de Febrero 210 bis, 2000 Rosario, Argentina  L. Martini
Instituto de Física Rosario, CONICET-UNR, Blvd. 27 de Febrero 210 bis, 2000
Rosario, Argentina  O.A. Fojón  Instituto de Física Rosario, CONICET-UNR,
Blvd. 27 de Febrero 210 bis, 2000 Rosario, Argentina  Escuela de Ciencias
Exactas y Naturales, FCEIA, Universidad Nacional de Rosario, Argentina

###### Abstract

Extracting single photoionization time delays associated with atomic (or
molecular) species from attosecond time scale two-photon experiments usually
relies on the theoretical description of continuum-continuum transitions. The
available models for those processes predict a universal phase contribution,
independent of the angular quantum numbers of final states. However, a recent
experimental-theoretical study [Fuchs, _et al._ Optica 7, 154 (2020)]
determined a sizable time delay dependence on the angular momentum of near-
threshold photoelectrons. In this study, we present an analytical model for
the two-photon two-color transition matrix amplitudes that reproduces the
phase dependence on the angular quantum number of final states. Finally, we
show that our analytical model can also describe the generalized Fano’s
propensity rules [Busto, _et al._ Phys. Rev. Lett. 123, 133201 (2019)] for
laser-assisted photoionization.

##  I Introduction

The pioneering work by Eisenbud, Wigner, and Smith (EWS) on interpreting the
energy derivative of the phase of the wavefunction resulted in the concept of
time delays associated with quantum processes [1, 2, 3]. For many years, that
concept remained a theoretical curiosity because its experimental
determination was not feasible: no measurement protocol allowed to extract
time shifts on the sub-femtosecond scale. However, that barrier was removed in
2001 when attosecond science became a reality [4, 5]. Since then, techniques
such as streaking [4] or the reconstruction of attosecond beating by two-
photon transitions (RABBITT) [6, 5] are the methods of choice for triggering
and monitoring electron dynamics in its intrinsic attosecond time scale [7].
Indeed, in 2010 Schultze _et al._ [8] applied the streaking protocol to
determine the single photoionization (relative) time delay of electrons
emerging from the $2s$ and $2p$ shells in neon atoms. Later on, Klunder _et
al._ [9] obtained the relative time delays for electrons ionized from the $3s$
and $3p$ shells in argon atoms, employing the RABBITT scheme.

The experimental realization of a longstanding theoretical prediction prompted
a new and still active research area within the ultrafast community that
ultimately led to a deeper understanding of the role played by time in quantum
mechanics. In fact, as the streaking and RABBITT schemes belong to the set of
pump-and-probe measurement protocols, it was soon understood that information
about the EWS time delay is only accessible if the influence of the infrared
probe field is considered [10]. Further theoretical developments show that
observed time delays in RABBITT experiments comprise two contributions: an
EWS-like delay associated with the pump stage and the additional delay
originated from the probing process [11].

These theoretical developments were crucial for the successful obtention of
(relative) EWS time delay from RABBITT traces, corresponding to different
atomic [9] and molecular [12] species. However, the asymptotic approximations
on which they rely are equivalent to neglecting the effects of the centrifugal
potential on intermediate and final states. Consequently, the continuum-
continuum (cc) phase for different final states becomes a single contribution,
independent of the angular quantum numbers of those states. This _universal_
character of the measurement-induced phases should hold as long as the pump
stage generates rapidly escaping electrons. In such a scenario, the most
relevant contributions to one-photon cc transitions come from space regions
where short-range potentials, such as the centrifugal one, are masked by the
long tail of the Coulomb interaction between the photoelectron and the parent
ion.

The above qualitative argument will be no longer valid for near-threshold
photoelectron kinetic energies. Indeed, recent measurements for slow
photoelectrons [13] confirmed the breakdown of the universal character of the
delay induced during the probing stage. The theoretical analysis of these
near-threshold processes [11, 14, 15, 13] can be exclusively carried out
through numerical simulations that, despite being extremely useful for an
accurate evaluation of the observables, usually provide only partial
explanations for the physical origin behind the observed effects. Therefore,
the development of a complementary analytical method can be a valuable tool
for providing further insights into these phenomena [16].

In this study, we present an analytical model that quantitatively describes
the dependence of two-photon two-color transition amplitudes on the angular
quantum numbers of final states.

##  II Methodology

To establish a connection with some previous experimental and theoretical
results on time delays, we focus on the RABBITT scheme: atomic targets in the
gas phase irradiated with a combination of odd harmonics of a Ti:sapphire
laser plus a synchronized weak fraction of the fundamental infrared (IR)
laser. In this situation, the photoelectron spectra display two qualitatively
distinct sets of lines: a) dressed harmonic (DH) lines mainly populated by the
absorption of one photon from the source of odd harmonics and b) sideband (SB)
lines, populated by two-photon transitions, on which the time delay gets
imprinted (see Fig. 1 for details).

Following standard time-dependent second-order perturbation theory, it is
possible to show that sideband signal reads [5],

$\displaystyle S_{2q}$
$\displaystyle\propto\sum_{L,M}|M_{L,M}^{+}+M_{L,M}^{-}|^{2},$ (1)

where $M_{L,M}^{\pm}$ are the two-photon transition matrix amplitudes
connecting the initial bound state to a final state in the continuum with
azimuthal and magnetic quantum numbers $L$ and $M$, respectively. In Fig. 1 we
schematically show these reaction pathways for a generic sideband of order
$2q$ and an initial $s$ $(l_{i}=0)$ state. After solving the angular algebra
and the time integrals for infinitely long fields, the transition matrix
amplitudes above can be factorized as

$\displaystyle M_{L,M}^{\pm}=A_{L,M}e^{\pm
i\omega_{0}\tau}\left[T_{L,\lambda}^{l_{i}}(\Omega_{\mp})+T_{L,\lambda}^{l_{i}}(\pm\omega_{0})\right],$
(2)

where $A_{L,M}$ is the angular (algebra) factor and $\exp(\pm
i\omega_{0}\tau)$ is a phase factor allowing to control the sideband signal by
changing the delay $\tau$ between the harmonics and the infrared laser with
angular frequency $\omega_{0}$. The (pseudo) radial matrix elements in length
gauge are given by

$\displaystyle
T_{L,\lambda}^{l_{i}}(\Omega)=\frac{e^{i\sigma_{L}(k)}}{i^{L}}\operatorname*{\mathchoice{\ooalign{$\displaystyle\sum$\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.14pt}{\scalebox{0.7}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.2pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}\frac{\langle
R_{\varepsilon_{k},L}|r|R_{\varepsilon_{\kappa},\lambda}\rangle\langle
R_{\varepsilon_{\kappa},\lambda}|r|R_{\varepsilon_{i},l_{i}}\rangle}{\varepsilon_{i}+\Omega-\varepsilon_{\kappa}},$
(11)

where $\sigma_{L}(k)=\arg[\Gamma(L+1-i/k)]$ is the Coulomb phase-shift for
partial wave $L$ and momentum $k$, and we include the phase factor in front of
the r.h.s. of the above equation only for further convenience. The summation
(integration) runs over the entire spectrum of unperturbed states of the
target, with energy $\varepsilon_{\kappa}$, angular momentum $\lambda$, and
radial wavefunctions $R_{\varepsilon_{\kappa},\lambda}$. The initial and final
states of the system are also described by eigenstates of the field-free
Hamiltonian, with radial wavefunctions $R_{\varepsilon_{i},l_{i}}$ and
$R_{\varepsilon_{k},L}$, respectively.

Figure 1: (Color online) Schematic representation of pathways contributing to
consecutive dressed harmonic lines and the sideband in between. For low IR
intensities, dressed harmonic lines are mainly populated by transitions from
ground state to continuum states, triggered by the absorption of a single
photon with energy $\Omega_{\pm}=(2q\pm 1)\omega_{0}$ from the attosecond
pulse train. The interaction of these primary photoelectron distributions with
the IR, gives rise to the sidebands. Time-reversed contributions, where the IR
photon is absorbed or emitted first are not shown here.

The transition matrix amplitudes in Eq. (2) can be derived in closed-form only
for hydrogenic atoms [17, 18, 19]. Therefore, to obtain a simple analytical
model able to describe a broader range of targets, some approximations must be
made. Firstly, for typical infrared laser intensities in RABBITT experiments,
we can focus only on those processes where the extreme ultraviolet photon from
the source of high-order harmonics is absorbed first. This is equivalent to
neglect the radial matrix element $T_{L,\lambda}^{l_{i}}(\pm\omega_{0})$ in
Eq. (2), that accounts for the dressing of the initial state due to the
infrared laser field [20]. The remaining radial matrix element can be
evaluated by invoking the Dalgarno-Lewis method [21], to obtain

$\displaystyle T_{L,\lambda}^{l_{i}}(\Omega_{\mp})=\langle
R_{\varepsilon_{k},L}|r|\rho_{\varepsilon_{\pm},\lambda}\rangle,$ (12)

with $\rho_{\varepsilon_{\pm},\lambda}$ being the solution of the
inhomogeneous differential equation
$[\mathcal{H}_{\lambda}-\varepsilon_{\pm}]\rho_{\varepsilon_{\pm},\lambda}=-rR_{\varepsilon_{i},l_{i}}$,
where $\mathcal{H}_{\lambda}$ is the laser-free radial Hamiltonian for angular
momentum $\lambda$, and the energy parameter satisfy $\varepsilon_{\pm}=(2q\mp
1)\omega_{0}+\varepsilon_{i}$. From the graphical representation in Fig. 1, it
can be seen that $\varepsilon_{\pm}=\kappa_{\pm}^{2}/2$ is merely the energy
for the lower (upper) dressed harmonic line.

The theoretical developments in Ref. [11] follow from approximating
$\rho_{\varepsilon_{\pm},\lambda}$ and $R_{\varepsilon_{k},L}$ in Eq. (12) by
their asymptotic forms corresponding to large $r$. Such approximation allows
to obtain analytical expressions for the radial matrix elements that factorize
as
$T_{L,\lambda}^{l_{i}}(\Omega_{\mp})=|T_{L,\lambda}^{l_{i}}(\Omega_{\mp})|\exp(i\phi^{\pm})$,
with $\phi^{\pm}=\phi^{\pm}_{bc}+\phi^{\pm}_{cc}$. The first term
$(\phi^{\pm}_{bc})$ is the phase acquired during the ionization step,
triggered by the absorption of one photon with energy $\Omega_{\mp}$. The
second term $(\phi^{\pm}_{cc})$ is the measurement-induced phase due to one-
photon cc transitions stimulated by the infrared probe laser. One main
contribution in [11] is the derivation of an analytical expression for the
$\phi^{\pm}_{cc}$ phases, allowing thus access to the EWS time delay. Briefly,
by measuring the time offset in the oscillation of sideband signals, induced
by the controlled modification of the delay $\tau$, a finite differences
approximation to the EWS time delay can be obtained by subtracting the
measurement-induced time delay $(\phi^{+}_{cc}-\phi^{-}_{cc})/2\omega_{0}$ to
the total offset 111Usually, in experiments, an additional phase due to the
intrinsic chirp of harmonics in the train has to be compensated. One possible
solution is to consider the relative time delay for electrons emerging from
different atomic subshells. See Refs. [9, 12] for further details..

Our derivation closely follows Ref. [11] to model the radial matrix elements
from Eq. (12), but we use the exact radial wavefunction
$R_{\varepsilon_{k},L}$ for the final states instead of its asymptotic
approximation. In that case, we get

$\displaystyle\begin{split}T_{L,\lambda}^{l_{i},\pm}\simeq&-\pi(-i)^{L}e^{i\sigma_{L}(k)}\langle
R_{\varepsilon_{\pm},\lambda}|r|R_{\varepsilon_{i},l_{i}}\rangle
N_{\varepsilon_{k},L}N_{\varepsilon_{\pm}}\\\
&\times\int_{0}^{\infty}dre^{-ikr}r^{L+2}\exp\left[i\kappa_{\pm}r+i\Phi_{\kappa_{\pm},\lambda}(r)\right]\\\
&\times\,_{1}F_{1}(L+1+i/k;2L+2;2ikr)\end{split}$ (13) $\displaystyle\simeq$
$\displaystyle-\pi\langle
R_{\varepsilon_{\pm},\lambda}|r|R_{\varepsilon_{i},l_{i}}\rangle
N_{\varepsilon_{k},L}N_{\varepsilon_{\pm}}I_{k,\kappa_{\pm}}^{L,\lambda},$
(14)

where the expression for the final continuum states was taken from Ref. [23],
${}_{p}F_{q}(a_{1},\ldots,a_{p};b_{1},\ldots,b_{q};z)$ is the generalized
Hypergeometric function, and the normalization factors for intermediate and
final states in the continuum are given by

$\displaystyle N_{\varepsilon_{\pm}}$
$\displaystyle=\sqrt{\frac{2}{\pi\kappa_{\pm}}},$ (15) $\displaystyle
N_{\varepsilon_{k},L}$
$\displaystyle=\frac{1}{(2L+1)!}\frac{2(2k)^{L}}{\sqrt{1-\exp(-2\pi/k)}}\prod_{s=1}^{L}\sqrt{s^{2}+\frac{1}{k^{2}}},$
(16)

respectively. The position-dependent phase in Eq. (13) is given by
$\Phi_{\kappa,\lambda}(r)=\ln(2\kappa
r)/\kappa+\sigma_{\lambda}(\kappa)-\pi\lambda/2$ with $\sigma_{l}(\kappa)$ the
Coulomb phase-shift for partial wave $l$ and momentum $\kappa$ [11]. Operating
over the logarithmic phase and rearranging the factors,
$I_{k,\kappa_{\pm}}^{L,\lambda}$ can be obtained in closed form by taking a
limit for the integral of a closely related family of functions [24]

$\displaystyle I_{k,\kappa}^{L,\lambda}=$
$\displaystyle(2\kappa)^{i/\kappa}e^{-i\pi(\lambda+L)/2}e^{i\sigma_{L}(k)}e^{i\sigma_{\lambda}(\kappa)}\lim_{\epsilon\to
0^{+}}\int_{0}^{\infty}dre^{-(\epsilon+ik-i\kappa)r}r^{L+2+i/\kappa}\,_{1}F_{1}(L+1+i/k;2L+2;2ikr)$
(17)
$\displaystyle\begin{split}=&(2\kappa)^{i/\kappa}(-i)^{L}e^{-i\pi\lambda/2}e^{i\sigma_{L}(k)}e^{i\sigma_{\lambda}(\kappa)}\frac{\Gamma(L+3+i/\kappa)}{(ik-i\kappa)^{(L+3+i/\kappa)}}\,_{2}F_{1}\left(L+1+\frac{i}{k},L+3+\frac{i}{\kappa};2L+2;\frac{2ik}{ik-i\kappa+0^{+}}\right).\end{split}$
(18)

As the other factors entering Eq. (14) are strictly real [11], the bound-
continuum phase acquired during the ionization step and the measurement
induced phase, for two-photon processes as those depicted in Fig. 1, must be
encoded in the expression for $I_{k,\kappa}^{L,\lambda}$ above. As a matter of
fact, from Eq. (18) we can recognize the $\phi^{\pm}_{bc}$ phase in the
$\exp{[i\sigma_{\lambda}(\kappa)]}$ phase factor and, therefore, our
analytical model allows to disentangle the EWS phase contribution while the cc
phases depend on the angular momentum of final states.

It should be noted that our model reduces to the model of Dahlström _et al._
[11] when the exact radial wavefunction $R_{\varepsilon_{k},L}$ for the final
states is replaced by its asymptotic form corresponding to large $kr$ values.
Once that operation is performed, the obtention of the two-photon radial
matrix elements follows along the same line as in [11].

##  III Continuum-continuum phases

The model we present can be applied, in principle, to approximate the two-
photon transition matrix amplitudes from any initial state of atomic [9] or
molecular species [12, 25]. This goal can be achieved by replacing the bound-
continuum contributions in Eqs. (6) and (10) with the corresponding numerical
results obtained from single-center expansion calculations [12]. However, to
facilitate the comparison with previous theoretical and experimental data, we
will firstly concentrate on initial hydrogen $s$ states. In that case, the
dipole selection rules indicate that only final states with angular momentum
$L=0,2$ can be populated.

The results we get for $\phi_{cc,L}^{\pm}$ from our analytic continuum-
continuum radial matrix elements (ACC-RME) in Eq. (18), are shown in Fig. 2 by
the dashed $(L=0)$ and dash-dotted $(L=2)$ lines. From extensive testing, we
have found that our model works slightly better if momentum for intermediate
states is approximated as $\kappa_{\pm}\simeq k\mp\omega_{0}/k$. Following
this procedure, related to the soft-photon approximation [11] for asymptotic
intermediate states, does not change any of the conclusions reached below. To
assess the accuracy of our analytical model, we additionally perform two sets
of _ab-initio_ calculations. Firstly, we numerically solve the time-dependent
Schrodinger equation (TDSE) in the velocity gauge for hydrogen atoms,
employing the Qprop code [26]. Using a single harmonic and the IR, we obtain
the transition amplitudes $M_{L,M}^{\pm}$ for the absorption and emission
pathways separately. Besides, we calculate the same two-photon transition
matrix amplitudes in length gauge from second-order perturbation theory
(SOPT), based on the analytical results in Ref. [19]. For both approaches, we
obtain the cc phase associated to each transition amplitude by subtracting the
contributions $\phi^{\pm}_{bc}$ and $\pm\omega_{0}\tau$ to the total phase
222For comparison purposes, ACC-RME and TDSE results with the same kinetic
energy, were shifted by a global phase.. As can be clearly seen in Fig. 2, the
results we get for the _ab-initio_ calculations 333The data that support the
findings of this study is publicly available on Repositorio Institucional
CONICET Digital http://hdl.handle.net/11336/156629 are virtually identical to
each other, ruling out any gauge and/or convergence systematic errors.

Figure 2: (Color online) Dependence of cc phases on photoelectron kinetic
energy, as obtained from different methods for hydrogen atom. Empty circles
$(L=0)$ and squares $(L=2)$: TDSE. Filled circles and squares: SOPT. Long-
range amplitude-corrected results are given by the full thick lines. ACC-RME
results are represented by the dashed lines $(L=0)$ and the dashed-dotted
lines $(L=2)$.

Generally speaking, ACC-RME results in Fig. 2 for the emission pathways are
better than those for the absorption ones. This behavior is not surprising and
reflects a diminished descriptive power of our model when it comes to
representing absorption pathways. From Fig. (1), it is evident that
intermediate states for emission pathways will always have a larger kinetic
energy when compared to the sideband final state. In that case, description of
the intermediate states with asymptotic wavefunctions qualitatively captures
the kinematical conditions under which every model is expected to work better:
the lower the kinetic energy, the better the quality of wavefunctions
describing that state should be. The situation is reversed for the absorption
pathways. Now, intermediate states with lower kinetic energy are described
through approximate asymptotic states, which may be a poor description for
sufficiently low kinetic energies.

Nonetheless, ACC-RME results for the absorption pathway and $L=2$
quantitatively reproduce the accurate data for electron kinetic energies above
$\varepsilon_{k}\sim 8$ eV, being almost identical to those from the long-
range amplitude-corrected (LRAC) asymptotic approximation in Ref. [11]. On the
other hand, the $L=0$ results for the absorption pathway slightly overestimate
_ab-initio_ calculations, with the notable exception of results around
$E_{k}\sim 5$ eV, that are almost perfectly described. For the emission
pathway, the quantitative agreement between the _ab-initio_ and ACC-RME
results in Fig. 2 is prominent, considering the simplicity of our approach.
For low kinetic energy, our model faithfully reproduce the phase splitting
obtained through the exact calculations, whereas at high energy the results
for $L=0,2$ converge to the expected limit given by the LRAC approximation.

The above comparison of individual cc phases $\phi_{cc,L}^{\pm}$ obtained
through different methods is only possible by theoretical analysis. From the
experimental point of view, phase differences are the most detailed pieces of
information that angular distributions of RABBITT traces encode [14, 13].
Recently, these phase differences were experimentally obtained for helium
targets [13], and the simulations reported therein allow to conclude that only
single-active electron dynamics are triggered within this energy range, in
agreement with previous calculations for two-photon two-color processes [14,
29]. Therefore, we can directly compare the phase differences
$\phi_{cc,0}^{\pm}-\phi_{cc,2}^{\pm}$ obtained from our model with those from
accurate calculations and the available experimental data.

These comparisons are presented in Fig. 3, for the absorption and emission
processes separately. Our model globally captures the salient features
displayed by the phase differences obtained from _ab-initio_ and experimental
data: they decrease to eventually zero for larger photoelectron kinetic
energies, and, for the same kinetic energy, the results for the absorption
pathways are slightly larger than those for the emission ones. The convergence
to zero can be easily demonstrated by recalling that our model reduces to that
in Ref. [11] when final states can be accurately described by its asymptotic
forms corresponding to large $kr$ values. In that case, continuum-continuum
phases are independent of the angular momentum of final states and, therefore,
their difference is identical to zero. Alternatively, this limit can be
obtained by just replacing the Gauss hypergeometric function in Eq. (18) with
the asymptotic approximation ${}_{2}F_{1}(a,b;c;z)\sim\lambda_{1}(-z)^{-a}$,
for $|z|\to\infty$ [30]. For this second approach, we also obtain cc phases
identical to those in Eq. (22) from Ref. [11]. On the contrary, the
comparatively larger values obtained for the absorption pathway cannot be
easily demonstrated as this feature heavily relies on fine details of the
Gauss hypergeometric function.

As before, we observe a better agreement for the emission pathway where the
ACC-RME results quantitatively reproduces the experimental data and _ab-
initio_ simulations, for electron kinetic energies above $\sim 5$ eV. For
smaller kinetic energies, our model captures the trend set by the simulations
for hydrogen and the experimental data for helium. The results for the
absorption pathway are less accurate, as ACC-RME calculations overestimate
both the experimental data and the _ab-initio_ calculations on the entire
range of electron kinetic energies under analysis. Again, the model can
reproduce the trend showing an increase in the phase difference for smaller
electron kinetic energies. The results for electron kinetic energies below
$\sim 5$ eV are shown in dashed line to emphasize that for an 800 nm infrared
laser, the intermediate states of the absorption pathway will have a kinetic
energy below $\sim 3.5$ eV, a value perhaps too low to admit a description
with asymptotic states.

Figure 3: (Color online) Dependence of cc phase differences on photoelectron
kinetic energy. ACC-RME, SOPT and TDSE results correspond to atomic hydrogen
targets. The experimental results correspond to helium atoms [13]. Upper
panel: absorption process. Lower panel: emission process.

##  IV Analysis of generalized Fano’s propensity rules in the continuum

The primary aim of developing this model was to put forward a general
analytical expression able to reproduce phase differences observed for partial
waves populating the final channel in sidebands. The previous analysis
demonstrates that this goal is achieved. In addition, as a by-product, with
our model we can also study the branching ratio for partial waves upon cc
transitions. This concept can be traced back, at least, to a few experimental
and theoretical studies for the half-RABBITT sidebands obtained when atomic or
molecular targets are ionized with a combination of a free-electron laser
pulse and an infrared probing field [31, 32]. More recently [15], this concept
was applied to RABBITT experiments, generalized for different partial waves,
and rationalized as a superset of Fano’s propensity rules [33]. The
generalized Fano propensity rule states that upon absorption (emission) of one
IR photon for the continuum-continuum transition, the partial wave with larger
(smaller) angular momentum is favored [15].

Figure 4: (Color online) Moduli ratios for the absorption (red) and emission
(blue) pathways. ACC-RME results for different intermediate states are given
by full lines. Second-order perturbation theory calculations are given by
empty stars [31] and full squares (emission only). TDSE simulation results are
shown by empty squares (emission only). RPAE results [15] for different
intermediate states are shown by empty triangles and hexagons.

In Fig. 4, we compare the moduli ratios
$|T_{\lambda-1,\lambda}^{l_{i},\pm}/T_{\lambda+1,\lambda}^{l_{i},\pm}|$
obtained from our model with the benchmark data: random-phase approximation
with exchange (RPAE) [15], TDSE and SOPT calculations and the scarce
experimental data available in the literature [31, 32]. In any case, we
observe that results obtained from our model reproduce the behavior expected
from the generalized Fano propensity rules. For absorption pathways, the
radial matrix elements for partial waves with $L=\lambda+1$ are larger than
those for $L=\lambda-1$ and, therefore, their ratio is smaller than its
asymptotic value of one. For the emission pathways, the opposite trend is
observed and partial waves with smaller angular momentum dominates. Even more,
our analytical results show that this behavior is dictated exclusively by
continuum-continuum transitions taking place in a effective potential
$(Z_{net}-1)/r$, with $Z_{net}$ being the net charge of the target,
irrespective of the more subtle details that its electronic structure may
have.

For $\lambda=1$, the model almost perfectly describes _ab-initio_ and
experimental results for the emission pathway on the entire energy range.
Conversely, the results for the absorption pathway fail for electron kinetic
energies below $\sim 5$ eV. As before, this behavior may be ascribed to a poor
description of low kinetic energy intermediate states with asymptotic
wavefunctions. The approximately linear dependence on $\lambda$ of the kinetic
energy value below which model results deviate from the accurate calculations
emphasizes the previous conclusion. Further understanding of this trend would
require a direct comparison of the approximate intermediate states with those
obtained from the numerical resolution of the Dalgarno-Lewis method, which is
beyond the scope of the present study.

The final remark on these branching ratios relates to its universal character.
From Eq. (14), it is clear that two-photon radial matrix elements factorize as
the product of a radial matrix element for the pump stage and a single-photon
cc factor. Therefore, the agreement of our model results with those from more
accurate calculations and the experimental data for the branching ratios in
Fig. 4 is a direct demonstration that the observed variations are dictated by
the single-photon cc transitions. In turn, it means that correlation effects,
included in the RPAE calculations [15] and experimental data [31, 32], play
little to no role in these branching ratios.

The above analysis of branching ratios for the radial matrix elements can be
systematically performed from a theoretical perspective. On the contrary, the
empirical verification of these results has some nuances depending on the
initial state. On one hand, for atomic $s$ states there will be a single
intermediate state with $\lambda=1$, therefore, the experimental and
theoretical results can be directly compared. The agreement we observe in Fig.
4 indicates that the branching ratios are independent of the radial
distribution for the initial $s$ state. On the other hand, measurements for
initial states with angular momentum larger than zero will be hindered by the
fact that final states with $L=l_{i}$ are coherently populated by
contributions following two different quantum pathways. In such a scenario, a
branching ratio dependence on the radial distribution of the initial state
cannot be ruled out. This is in agreement with the findings in Ref. [25] (Fig.
5), where changing the radial distribution for a $1s$ initial state in
hydrogen does not modify an observable closely related to radial matrix
elements, whereas the same procedure for initial $2p$ orbitals induces a
dramatic change.

##  V Conclusions and outlook

To summarize, we develop an analytical model for the two-photon two-color
transition matrix amplitudes that quantitatively reproduce state-of-the-art
calculations and experimental data at low photoelectron kinetic energies,
without resorting to fitting free parameters. We demonstrate that the effects
induced by the centrifugal barrier onto the wavefunctions can be almost
completely recovered by just using accurate final states, even if intermediate
states are approximated by its asymptotic form corresponding to large $r$.
Besides, the functional form we obtain for the matrix elements is fully
compatible with the procedure to extract the EWS time delays from the sideband
signal oscillation in RABBITT experiments. Additionally, our results can be
readily adapted to upgrade models describing higher order cc processes [34],
as it only requires to change the description of final states. Finally, we
demonstrate that the partial waves branching ratio above a $\lambda$-dependent
photoelectron kinetic energy can also be accurately retrieved by considering
the centrifugal potential effects only on final states. In particular, our
model may be useful to disentangle the coherent contributions populating final
state partial waves with $L=l_{i}$ in two-color experiments from free-electron
laser facilities, where measurements can be performed away from threshold.

###### Acknowledgements.

The authors acknowledge the financial support from the Agencia Nacional de
Promoción Científica y Tecnológica, Project PICT No. 01912 (2015-3392), and
from the Consejo Nacional de Investigaciones Científicas y Técnicas de la
República Argentina, Project PIP No. 0784 (11220090101026). Part of the
results presented in this work have been obtained by using the facilities of
the CCT-Rosario Computational Center, member of the High Performance Computing
National System (SNCAD, MincyT-Argentina).

## References

  * Wigner [1955] E. P. Wigner, Lower limit for the energy derivative of the scattering phase shift, Phys. Rev. 98, 145 (1955).
  * Smith [1960] F. T. Smith, Lifetime matrix in collision theory, Phys. Rev. 118, 349 (1960).
  * de Carvalho and Nussenzveig [2002] C. de Carvalho and H. Nussenzveig, Time delay, Physics Reports 364, 83 (2002).
  * Hentschel  _et al._ [2001] M. Hentschel, R. Kienberger, C. Spielmann, G. A. Reider, N. Milosevic, T. Brabec, P. Corkum, U. Heinzmann, M. Drescher, and F. Krausz, Attosecond metrology, Nature 414, 509 (2001).
  * Paul  _et al._ [2001] P. M. Paul, E. S. Toma, P. Breger, G. Mullot, F. Augé, P. Balcou, H. G. Muller, and P. Agostini, Observation of a train of attosecond pulses from high harmonic generation, Science 292, 1689 (2001).
  * Véniard  _et al._ [1996] V. Véniard, R. Taïeb, and A. Maquet, Phase dependence of (N+1)-color (N$>$1) ir-xuv photoionization of atoms with higher harmonics, Phys. Rev. A 54, 721 (1996).
  * Goulielmakis  _et al._ [2010] E. Goulielmakis, Z.-H. Loh, A. Wirth, R. Santra, N. Rohringer, V. S. Yakovlev, S. Zherebtsov, T. Pfeifer, A. M. Azzeer, M. F. Kling, S. R. Leone, and F. Krausz, Real-time observation of valence electron motion, Nature 466, 739 (2010).
  * Schultze  _et al._ [2010] M. Schultze, M. Fieß, N. Karpowicz, J. Gagnon, M. Korbman, M. Hofstetter, S. Neppl, A. L. Cavalieri, Y. Komninos, T. Mercouris, C. A. Nicolaides, R. Pazourek, S. Nagele, J. Feist, J. Burgdörfer, A. M. Azzeer, R. Ernstorfer, R. Kienberger, U. Kleineberg, E. Goulielmakis, F. Krausz, and V. S. Yakovlev, Delay in photoemission, Science 328, 1658 (2010).
  * Klünder  _et al._ [2011] K. Klünder, J. M. Dahlström, M. Gisselbrecht, T. Fordell, M. Swoboda, D. Guénot, P. Johnsson, J. Caillat, J. Mauritsson, A. Maquet, R. Taïeb, and A. L’Huillier, Probing single-photon ionization on the attosecond time scale, Phys. Rev. Lett. 106, 143002 (2011).
  * Nagele  _et al._ [2011] S. Nagele, R. Pazourek, J. Feist, K. Doblhoff-Dier, C. Lemell, K. Tőkési, and J. Burgdörfer, Time-resolved photoemission by attosecond streaking: extraction of time information, Journal of Physics B: Atomic, Molecular and Optical Physics 44, 081001 (2011).
  * Dahlström  _et al._ [2013] J. Dahlström, D. Guénot, K. Klünder, M. Gisselbrecht, J. Mauritsson, A. L’Huillier, A. Maquet, and R. Taïeb, Theory of attosecond delays in laser-assisted photoionization, Chemical Physics 414, 53 (2013), attosecond spectroscopy.
  * Huppert  _et al._ [2016] M. Huppert, I. Jordan, D. Baykusheva, A. von Conta, and H. J. Wörner, Attosecond delays in molecular photoionization, Phys. Rev. Lett. 117, 093001 (2016).
  * Fuchs  _et al._ [2020] J. Fuchs, N. Douguet, S. Donsa, F. Martin, J. Burgdörfer, L. Argenti, L. Cattaneo, and U. Keller, Time delays from one-photon transitions in the continuum, Optica 7, 154 (2020).
  * Heuser  _et al._ [2016] S. Heuser, A. Jiménez Galán, C. Cirelli, C. Marante, M. Sabbar, R. Boge, M. Lucchini, L. Gallmann, I. Ivanov, A. S. Kheifets, J. M. Dahlström, E. Lindroth, L. Argenti, F. Martín, and U. Keller, Angular dependence of photoemission time delay in helium, Phys. Rev. A 94, 063409 (2016).
  * Busto  _et al._ [2019] D. Busto, J. Vinbladh, S. Zhong, M. Isinger, S. Nandi, S. Maclot, P. Johnsson, M. Gisselbrecht, A. L’Huillier, E. Lindroth, and J. M. Dahlström, Fano’s propensity rule in angle-resolved attosecond pump-probe photoionization, Phys. Rev. Lett. 123, 133201 (2019).
  * Armstrong  _et al._ [2021] G. S. J. Armstrong, M. A. Khokhlova, M. Labeye, A. S. Maxwell, E. Pisanty, and M. Ruberti, Dialogue on analytical and ab initio methods in attoscience, The European Physical Journal D 75, 209 (2021).
  * Klarsfeld [1970] S. Klarsfeld, Two-photon ionization of atomic hydrogen in the ground state, Lettere al Nuovo Cimento (1969-1970) 3, 395 (1970).
  * Arnous  _et al._ [1973] E. Arnous, S. Klarsfeld, and S. Wane, Angular distribution in the two-quantum atomic photoeffect, Phys. Rev. A 7, 1559 (1973).
  * Jayadevan and Thayyullathil [2001] A. P. Jayadevan and R. B. Thayyullathil, Two-photon ionization of atomic hydrogen above the one-photon ionization threshold, Journal of Physics B: Atomic, Molecular and Optical Physics 34, 699 (2001).
  * Cionga  _et al._ [1993] A. Cionga, V. Florescu, A. Maquet, and R. Taïeb, Target dressing effects in laser-assisted x-ray photoionization, Phys. Rev. A 47, 1830 (1993).
  * Dalgarno  _et al._ [1955] A. Dalgarno, J. T. Lewis, and D. R. Bates, The exact calculation of long-range forces between atoms by perturbation theory, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 233, 70 (1955).
  * Note [1] Usually, in experiments, an additional phase due to the intrinsic chirp of harmonics in the train has to be compensated. One possible solution is to consider the relative time delay for electrons emerging from different atomic subshells. See Refs. [9, 12] for further details.
  * Matsumoto [1991] A. Matsumoto, Multipole matrix elements for hydrogen atom, Physica Scripta 44, 154 (1991).
  * Landau and Lifshitz [1991] L. D. Landau and E. M. Lifshitz,  _Quantum mechanics: non-relativistic theory_ , 3rd ed., Course of theoretical physics No. 3 (Pergamon press, 1991).
  * Martini  _et al._ [2021] L. Martini, D. I. R. Boll, A. Palacios, and O. A. Fojón, Different approximations in the polarization control of two-color water molecule ionization, Journal of Physics B: Atomic, Molecular and Optical Physics 54, 135602 (2021).
  * Tulsky and Bauer [2020] V. Tulsky and D. Bauer, Qprop with faster calculation of photoelectron spectra, Computer Physics Communications 251, 107098 (2020).
  * Note [2] For comparison purposes, ACC-RME and TDSE results with the same kinetic energy, were shifted by a global phase.
  * Note [3] The data that support the findings of this study is publicly available on Repositorio Institucional CONICET Digital http://hdl.handle.net/11336/156629.
  * Boll  _et al._ [2020] D. I. R. Boll, L. Martini, O. A. Fojón, and A. Palacios, Off-resonance-enhanced polarization control in two-color atomic ionization, Phys. Rev. A 101, 013428 (2020).
  * Abramowitz and Stegun [1964] M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ , ninth dover printing, tenth gpo printing ed. (Dover, New York, 1964).
  * Meyer  _et al._ [2010] M. Meyer, D. Cubaynes, J. Dardis, P. Hayden, P. Hough, V. Richardson, E. Kennedy, J. Costello, S. Düsterer, W. Li, P. Radcliffe, H. Redlin, J. Feldhaus, S. Strakhova, E. Gryzlova, A. Grum-Grzhimailo, R. Taïeb, and A. Maquet, Two-color experiments in the gas phase at flash, Journal of Electron Spectroscopy and Related Phenomena 181, 111 (2010), electronic Spectroscopy and Structure (ICESS-11).
  * Fushitani  _et al._ [2012] M. Fushitani, A. Matsuda, and A. Hishikawa, Euv and soft x-ray photoelectron spectroscopy of isolated atoms and molecules using single-order laser high-harmonics at 42ev and 91ev, Journal of Electron Spectroscopy and Related Phenomena 184, 561 (2012).
  * Fano [1985] U. Fano, Propensity rules: An analytical approach, Phys. Rev. A 32, 617 (1985).
  * Bharti  _et al._ [2021] D. Bharti, D. Atri-Schuller, G. Menning, K. R. Hamilton, R. Moshammer, T. Pfeifer, N. Douguet, K. Bartschat, and A. Harth, Decomposition of the transition phase in multi-sideband schemes for reconstruction of attosecond beating by interference of two-photon transitions, Phys. Rev. A 103, 022834 (2021).


# A simplified algorithm for the topological entropy of multimodal maps

José M. Amigó
Centro de Investigación Operativa,
Universidad Miguel Hernández,
03202 Elche, Spain
<EMAIL_ADDRESS>&Angel Giménez
Centro de Investigación Operativa,
Universidad Miguel Hernández,
03202 Elche, Spain
<EMAIL_ADDRESS>
###### Abstract

A numerical algorithm to compute the topological entropy of multimodal maps is
proposed. This algorithm results from a closed formula containing the so-
called min-max symbols, which are closely related to the kneading symbols.
Furthermore, it simplifies a previous algorithm, also based on min-max
symbols, which was originally proposed for twice differentiable multimodal
maps. The new algorithm has been benchmarked against the old one with a number
of multimodal maps, the results being reported in the paper. In particular,
the comparison is favorable to the new algorithm, except in the unimodal case.

_K_ eywords Topological entropy; Multimodal maps; Min-max symbols.

##  1 Introduction

Let $f$ be a continuous selfmap of a compact interval $[a,b]\subset\mathbb{R}$
with a finite number of turning (or critical) points. Such maps are
generically called multimodal. Then, the topological entropy of $f$ [1, 2],
$h(f)$, can be calculated (along other possibilities) with the formula

$h(f)=\lim_{n\rightarrow\infty}\frac{1}{n}\log\ell_{n},$ (1)

where $\ell_{n}$ is shorthand for the lap number of $f^{n}$, i.e., the number
of maximal monotonicity segments of $f^{n}$, the $n$th iterate of $f$ [3, 4].

In [5, Sect. 7] a numerical algorithm to compute the topological entropy of
multimodal maps was proposed. Let us point out that this algorithm generalizes
and hence includes a previous one for unimodal maps published in [6]. The
algorithm builds on (1) by calculating $\ell_{n}$ with the help of the min-max
symbols of $f$ [5, 6, 7, 8], a generalization of the kneading symbols [9, 10].
The min-max symbols of a multimodal map not only locate the iterates of its
critical values up to the precision set by the partition defined by its
critical points, as the kneading symbols do, but they also display their
minimum/maximum (or “critical”) character. The interesting point is that such
an additional information supposes virtually no extra computational cost.
Indeed, it can be read recursively from a look-up table once the min-max
symbols of the critical values are known.

In this paper we propose a related algorithm which actually approximates the
value of $h(f)$ given by a closed formula involving also the min-max symbols
of $f$. The new algorithm eliminates a formal restriction that, as it turns
out, unnecessarily marred the applicability of the algorithm of [5]. At the
same time, it simplifies the computation scheme of the latter. We elaborate
upon these two points briefly.

With regard to the formal restriction, the theoretical results of [5] refer to
twice differentiable multimodal maps only. However, numerical simulations with
piecewise linear maps of constant slope (and alternating sign) suggested that
the algorithm of [5] could be applied as well to just continuous maps. In this
paper we justify the extension of results from smooth to just continuous maps.
Although the proof turns out to be straightforward, this generalization was
not explored in the previous papers [5, 6] just because these followed the
original approach in [7, 8], which only considered twice differentiable maps
for simplicity.

As for the simplification of the computation scheme, this has to do with the
boundary conditions (or the lack of them). Indeed, the algorithm of [5] keeps
track of the orbits of the boundary points, thus calculating the exact value
of the lap number $\ell_{n}$ in each computation loop. At the contrary, the
new algorithm dispenses with those orbits because they do not affect the limit
(1). In fact, Theorem 3 below shows that, as far as the computation of $h(f)$
is concerned, one may assume that $f$ is boundary-anchored, i.e.,
$f(\\{a,b\\})\subset\\{a,b\\}$. The result is a compact expression for the lap
number $\ell_{n}$ that makes possible a closed formula for $h(f)$.

In sum, we fill a theoretical gap in the application of the algorithm in [5]
by showing that continuity of the maps suffices. Moreover, we abridge the
numerical scheme by approximating $\ell_{n}$ in (1) with a formula which is
exact only for boundary-anchored maps but provides the right limit (1) for
$h(f)$.

This paper is organized as follows. To make the paper self-contained, we
review in Sect. 2 all the basic concepts, especially the concept of min-max
sequences, needed in the following sections. Most importantly, we extend in
Theorem 1 the transition rules for min-max symbols from twice differentiable
multimodal maps [5] to just continuous ones. In Sect. 3 we introduce some
instrumental results which lead in Sect. 4, together with Theorem 3, to a
closed formula for $h(f)$ containing the min-max symbols of $f$ (Theorem 4). A
formal proof of Theorem 3 has been shifted to the Appendix in order not to
interrupt the flow of ideas. Sect. 5 contains the main result of the paper,
namely, an algorithm for the topological entropy of (not necessarily smooth)
multimodal maps which approximates the value of $h(f)$ given in Theorem 4. As
way of illustration, this algorithm is put to test in Sect. 6. First, the new,
abridged algorithm is benchmarked in Sects. 6.1 to 6.3 against the full-
pledged one, Ref. [5], using smooth uni-, bi-, and trimodal maps,
respectively, borrowed from [5, 6]. Finally, in Sect. 6.4, we also compare
both algorithms via piecewise linear, four and five-modal maps of known
topological entropy. It turns out that, except in the unimodal case, the new
algorithm outperforms the old one.

##  2 Min-max sequences

For the reader’s convenience, we use the same notation as in [5] throughout.
Let $I$ be a compact interval $[a,b]\subset\mathbb{R}$ and $f:I\rightarrow I$
a piecewise monotone continuous map. Such a map is called $l$-modal if $f$ has
precisely $l$ turning points (i.e., points in $(a,b)$ where $f$ has a local
extremum). Sometimes we speak also of multimodal maps, in general, or of
unimodal maps in the particular case $l=1$. Furthermore, assume henceforth
that $f$ has local extrema at $c_{1}<...<c_{l}$ and is strictly monotone in
each of the $l+1$ intervals

$I_{1}=[a,c_{1}),I_{2}=(c_{1},c_{2}),...,I_{l}=(c_{l-1},c_{l}),I_{l+1}=(c_{l},b]\text{.}$

In this case we write $f\in\mathcal{M}_{l}(I).$ When the interval $I$ is clear
from the context or unimportant for the argument, we write just
$\mathcal{M}_{l}$.

Since the results we obtain below for the calculation of the topological
entropy do not depend on the shape of $f$, i.e., on whether $f(c_{1}$) is a
maximum (positive shape) or a minimum (negative shape), we assume, unless
otherwise stated, that $f$ has positive shape. This implies that
$f(c_{\text{odd}})$ are maxima, whereas $f(c_{\text{even}})$ are minima, where
here and hereafter “even” and “odd” stand for even and odd subindices,
respectively. Hence $f$ is strictly increasing on the intervals
$I_{\text{odd}}$, and strictly decreasing on the intervals $I_{\text{even}}$.
The points $f(c_{i})$, $1\leq i\leq l$, are called the critical values of $f$
although no differentiability of $f$ at $c_{i}$ is assumed when so doing.

Theorem 1. Let $f\in\mathcal{M}_{l}$ have a positive shape, and $n\geq 1$.
Then:

(a)

If $f^{n}(x)=c_{\text{odd}}$ then $f^{n+1}(x)$ is a maximum. If
$f^{n}(x)=c_{\text{even}}$ then $f^{n+1}(x)$ is a minimum.

(b)

If $f^{n}(x)$ is a minimum, then

$f^{n+1}(x)\text{ is a }\left\\{\begin{array}[]{cl}\text{minimum}&\text{if
}f^{n}(x)\in I_{\text{odd}},\\\ \text{maximum}&\text{if }f^{n}(x)\in
I_{\text{even}}.\end{array}\right.$

(c)

If $f^{n}(x)$ is a maximum, then

$f^{n+1}(x)\text{ is a }\left\\{\begin{array}[]{cl}\text{maximum}&\text{if
}f^{n}(x)\in I_{\text{odd}},\\\ \text{minimum}&\text{if }f^{n}(x)\in
I_{\text{even}}.\end{array}\right.$

Proof. (a) This is a trivial consequence of $f$ having a positive shape.

(b) Suppose that $f^{n}(x_{0})$ is a minimum with $f^{n}(x_{0})\in
I_{\text{odd}}$. Therefore, there exists a neighborhood of $x_{0}$,
$U(x_{0})$, such that $f^{n}(x_{0})\leq f^{n}(x)$ for all $x\in U(x_{0})$.
Without restriction we may assume that $U(x_{0})\subset
f^{-n}(I_{\text{odd}})$. It follows that

$f^{n+1}(x_{0})=f(f^{n}(x_{0}))\leq f(f^{n}(x))=f^{n+1}(x)$

for all $x\in U(x_{0})$ because $f^{n}(U(x_{0}))\subset I_{\text{odd}}$, an
interval where $f$ is increasing. We conclude that $f^{n+1}(x_{0})$ is a
minimum.

If $f^{n}(x_{0})\in I_{\text{even}}$ then we derive from $f^{n}(x_{0})\leq
f^{n}(x)$ that $f^{n+1}(x_{0})\geq f^{n+1}(x)$ for all $x\in U(x_{0})$ because
this time $f^{n}(U(x_{0}))\subset I_{\text{even}}$, an interval where $f$ is
decreasing.

(c) This case follows similarly to (b). $\square$

The itinerary of $x\in I$ under $f$ is a symbolic sequence

$\mathbf{i}(x)=(i_{0}(x),i_{1}(x),...,i_{n}(x),...)\in\\{I_{1},c_{1},I_{2},...,c_{l},I_{l+1}\\}^{\mathbb{N}_{0}}$

($\mathbb{N}_{0}\equiv\\{0\\}\cup\mathbb{N}$), defined as follows:

$i_{n}(x)=\left\\{\begin{array}[]{cl}I_{j}&\text{if }f^{n}(x)\in I_{j}\text{
}(1\leq j\leq l+1),\\\ c_{k}&\text{if }f^{n}(x)=c_{k}\text{ }(1\leq k\leq
l).\end{array}\right.$

The itineraries of the critical values,

$\gamma^{i}=(\gamma_{1}^{i},...,\gamma_{n}^{i},...)=\mathbf{i}(f(c_{i})),\;1\leq
i\leq l,$

are called the kneading sequences [9, 10] of $f$.

Definition 1 [5, 6, 7, 8]. The min-max sequences of an $l$-modal map $f$,

$\omega^{i}=(\omega_{1}^{i},\omega_{2}^{i},...,\omega_{n}^{i},...),\;1\leq
i\leq l,$

are defined as follows:

$\omega_{n}^{i}=\left\\{\begin{array}[]{ll}m^{\mathbf{\gamma}_{n}^{i}}&\text{if
}f^{n}(c_{i})\text{ is a minimum},\\\ M^{\mathbf{\gamma}_{n}^{i}}&\text{if
}f^{n}(c_{i})\text{ is a maximum.}\end{array}\right.$

where $\gamma_{n}^{i}$ are kneading symbols.

Thus, the min-max symbols $\omega_{n}^{i}$ have an exponential-like notation,
where the ‘base’ belongs to the alphabet $\\{m,M\\}$, and the ‘exponent’ is a
kneading symbol. The extra information of a min-max symbol $\omega_{n}^{i}$ as
compared to a kneading symbol $\gamma_{n}^{i}$ is contained, therefore, in the
base, which tell us whether $f^{n}(c_{i})$ is a minimum ($m$) or a maximum
($M$). Theorem 1 shows that once the symbol $\omega_{n}^{i}$ of a map with
positive shape is known, the symbol $\omega_{n+1}^{i}$ can be read from the
table

$\omega_{n}^{i}$ | $\rightarrow$ | $\omega_{n+1}^{i}$
---|---|---
$m^{c_{\text{even}}},M^{c_{\text{even}}}$ | $\rightarrow$ | $m^{\gamma_{n+1}^{i}}$
$m^{c_{\text{odd}}},M^{c_{\text{odd}}}$ | $\rightarrow$ | $M^{\gamma_{n+1}^{i}}$
$m^{I_{\text{odd}}},M^{I_{\text{even}}}$ | $\rightarrow$ | $m^{\gamma_{n+1}^{i}}$
$m^{I_{\text{even}}},M^{I_{\text{odd}}}$ | $\rightarrow$ | $M^{\gamma_{n+1}^{i}}$
(2)

Let us mention for completeness that if $f\in\mathcal{M}_{l}$ has negative
shape, then the transition rules from $\omega_{n}^{i}$ to $\omega_{n+1}^{i}$
read

$\omega_{n}^{i}$ | $\rightarrow$ | $\omega_{n+1}^{i}$
---|---|---
$m^{c_{\text{even}}},M^{c_{\text{even}}}$ | $\rightarrow$ | $M^{\gamma_{n+1}^{i}}$
$m^{c_{\text{odd}}},M^{c_{\text{odd}}}$ | $\rightarrow$ | $m^{\gamma_{n+1}^{i}}$
$m^{I_{\text{odd}}},M^{I_{\text{even}}}$ | $\rightarrow$ | $M^{\gamma_{n+1}^{i}}$
$m^{I_{\text{even}}},M^{I_{\text{odd}}}$ | $\rightarrow$ | $m^{\gamma_{n+1}^{i}}$
(3)

instead. This follows mutatis mutandis as in the proof of Theorem 1 .

The transition rules (2) and (3) substantiate our claim in the Introduction
that, from the point of view of the computational cost, min-max sequences and
kneading sequences are virtually equivalent.

Therefore, the kneading symbols of $f\in\mathcal{M}_{l}$, along with its
initial min-max symbols, i.e.

$\omega_{1}^{i}=\left\\{\begin{array}[]{ll}M^{\gamma_{1}^{i}}&\text{if
}i=1,3,...,2\left\lfloor\frac{l+1}{2}\right\rfloor-1,\\\
m^{\gamma_{1}^{i}}&\text{if
}i=2,4,...,2\left\lfloor\frac{l}{2}\right\rfloor,\end{array}\right.\;$ (4)

and the transition rules (2) allow to compute the min-max sequences of
$f\in\mathcal{M}_{l}$ in a recursive way.

A final ingredient (proper of min-max sequences) is the following. Let the
$i$th critical line, $1\leq i\leq l$, be the line $y=c_{i}$ in the Cartesian
product $I\times I$. Min-max symbols split into bad and good symbols with
respect to $i$th critical line. Geometrically, good symbols correspond to
local maxima strictly above the line $y=c_{i}$, or to local minima strictly
below the line $y=c_{i}$. All other min-max symbols are bad by definition with
respect to the $i$th critical line. We use the notation

$\mathcal{B}^{i}=\\{M^{I_{1}},M^{c_{1}},...,M^{I_{i}},M^{c_{i}},m^{c_{i}},m^{I_{i+1}},...,m^{c_{l}},m^{I_{l+1}}\\}$

for the set of bad symbols of $f\in\mathcal{M}_{l}$ with respect to the $i$th
critical line. There are $2(l+1)$ bad symbols and $2l$ good symbols with
respect to a given critical line.

Bad symbols appear in all results of [5, 6] concerning the computation of the
topological entropy of $f\in\mathcal{M}_{l}$ via min-max symbols. In this
sense we may say that bad symbols are the hallmark of this approach.

##  3 Auxiliary results

Let $s_{\nu}^{i}$, $1\leq i\leq l$, stand for the number of interior simple
zeros of $f^{\nu}(x)-c_{i}$, $\nu\geq 0$, i.e., solutions of $x-c_{i}=0$
($\nu=0$), or solutions of $f^{\nu}(x)=c_{i}$, $x\in(a,b)$, with
$f^{\mu}(x)\neq c_{i}$ for $0\leq\mu\leq\nu-1$, and $f^{\nu\prime}(x)\not=0$
($\nu\geq 1$). Geometrically $s_{\nu}^{i}$ is the number of transversal
intersections on the Cartesian plane $(x,y)$ of the curve $y=f^{\nu}(x)$ and
the straight line $y=c_{i}$, over the interval $(a,b)$. Note that
$s_{0}^{i}=1$ for all $i$.

To streamline the notation set

$s_{\nu}=\sum_{i=1}^{l}s_{\nu}^{i}$ (5)

for $\nu\geq 0$. In particular,

$s_{0}=\sum_{i=1}^{l}s_{0}^{i}=\sum_{i=1}^{l}1=l.$ (6)

According to [5, Eqn. (31)], the lap number of $f^{n}$, $\ell_{n}$, satisfies

$\ell_{n}=1+\sum\limits_{\nu=0}^{n-1}s_{\nu}=\ell_{n-1}+s_{n-1},$ (7)

for $n\geq 1$. In particular, $\ell_{1}=\ell_{0}+s_{0}=1+l$.

Furthermore, define

$K_{\nu}^{i}=\\{(k,\kappa),1\leq k\leq
l,1\leq\kappa\leq\nu:\omega_{\kappa}^{k}\in\mathcal{B}^{i}\\},$ (8)

($\nu\geq 1$, $1\leq i\leq l$), that is, $K_{\nu}^{i}$ collects the upper and
lower indices $(k,\kappa)$ of the bad symbols with respect to the $i$th
critical line in all the initial blocks

$\omega_{1}^{1},\omega_{2}^{1},...,\omega_{\nu}^{1};\;\;\omega_{1}^{2},\omega_{2}^{2},...,\omega_{\nu}^{2};\;\;...;\;\;\omega_{1}^{l},\omega_{2}^{l},...,\omega_{\nu}^{l};$

of the min-max sequences of $f$. We note for further reference that
$K_{\nu-1}^{i}\subset K_{\nu}^{i}$, the set-theoretical difference being

$K_{\nu}^{i}\backslash K_{\nu-1}^{i}=\\{(k,\nu),1\leq k\leq
l:\omega_{\nu}^{k}\in\mathcal{B}^{i}\\}.$ (9)

Finally, set

$S_{\nu}^{i}=2\sum_{(k,\kappa)\in K_{\nu}^{i}}s_{\nu-\kappa}^{k},$ (10)

where $S_{\nu}^{i}=0$ if $K_{\nu}^{i}=\emptyset$, and analogously to (5),

$S_{\nu}=\sum_{i=1}^{l}S_{\nu}^{i}.$ (11)

We say that $f\in\mathcal{M}_{l}$ is boundary-anchored if
$f\\{a,b\\}\subset\\{a,b\\}$. This boundary condition boils down to

$f(a)=a,\;\text{and }f(b)=\left\\{\begin{array}[]{ll}a&\text{if }l\text{ is
odd,}\\\ b&\text{if }l\text{ is even.}\end{array}\right.$ (12)

for multimodal maps with positive shape, and to

$f(a)=b,\;\text{and }f(b)=\left\\{\begin{array}[]{ll}b&\text{if }l\text{ is
odd,}\\\ a&\text{if }l\text{ is even.}\end{array}\right.$ (13)

for multimodal maps with negative shape. As we will see shortly, boundary-
anchored maps have some advantages when calculating the topological entropy.

Theorem 2. Let $f\in\mathcal{M}_{l}$ be boundary-anchored. Then

$s_{\nu}^{i}=1+\sum_{\mu=0}^{\nu-1}s_{\mu}-S_{\nu}^{i},$ (14)

Proof. Suppose for the time being that $f$ is twice differentiable on $[a,b]$
without any restriction at the boundaries. In this case, it was proved in [5,
Theorem 5.3] that

$s_{\nu}^{i}=1+\sum_{\mu=0}^{\nu-1}s_{\mu}-S_{\nu}^{i}-\alpha_{\nu}^{i}-\beta_{\nu}^{i},$
(15)

where $\alpha_{\nu}^{i}$,$\beta_{\nu}^{i}$ are binary variables that vanish if
$f$ is boundary-anchored. Moreover, the relation (15) follows from the
transition rules (2) (or (3) if $f$ has negative shape), which have been
proved to hold true also for continuous multimodal maps in Theorem 1. It
follows that (15) holds for continuous, multimodal maps as well. In
particular, (14) holds for the boundary-anchored ones. $\square$

Therefore, the boundary conditions (12) allow us to express $s_{\nu}^{i}$ with
the help of some $s_{0}^{k}$, $s_{1}^{k}$, …, $s_{\nu-1}^{k}$, $1\leq k\leq
l$, via (14) and (10). Sum (14) over $i$ from $1$ to $l$ to obtain the compact
relation

$s_{\nu}=l\left(1+\sum\limits_{\mu=0}^{\nu-1}s_{\mu}\right)-S_{\nu}$ (16)

between $s_{0}=l,s_{1},...,s_{\nu}$ and $S_{\nu}$, for all $\nu\geq 1$. By (7)
this equation can we rewritten as $s_{\nu}=l\ell_{\nu}-S_{\nu}$, hence

$\ell_{\nu}=\frac{1}{l}(s_{\nu}+S_{\nu}).$ (17)

##  4 A closed formula for the topological entropy of multimodal maps

According to [11, Lemma 4.4], the topological entropy of a multimodal map
depends only on the kneading sequences, i.e., on the itineraries of the
critical values, but not on the itineraries of the boundary points. This
entails that one may assume without restriction the boundary conditions (12)
or (13) when calculating the topological entropy of $l$-modal maps with
positive or negative shape, respectively. A formal justification is given by
the following theorem.

Theorem 3. Let $f\in\mathcal{M}_{l}(I)$. Then there exists
$F\in\mathcal{M}_{l}(J)$, where $J\supset I$, such that $h(F)=h(f)$ and $F$ is
boundary-anchored.

See [9, Lemma 7.7], and [11, proof of Lemma 4.4]. For the reader’s
convenience, a proof of Theorem 3 is given in the Appendix.

This being the case, (1) and (17) yield the following result.

Theorem 4. Let $f\in\mathcal{M}_{l}$. Then,

$\displaystyle h(f)$
$\displaystyle=\lim_{\nu\rightarrow\infty}\frac{1}{\nu}\log\frac{s_{\nu}+S_{\nu}}{l}$
(18)
$\displaystyle=\lim_{\nu\rightarrow\infty}\frac{1}{\nu}\log\frac{1}{l}\sum_{i=1}^{l}\left(s_{\nu}^{i}+2\sum_{(k,\kappa)\in
K_{\nu}^{i}}s_{\nu-\kappa}^{k}\right).$

Eqn. (18) provides a closed expression for $h(f)$ which includes the min-max
symbols of $f$.

##  5 A simplified algorithm for the topological entropy

An offshoot of the preceding section is that, when it comes to calculate the
topological entropy of a multimodal map, one can resort to the limit (18),
whether the map is boundary-anchored or not. Loosely speaking,

$h(f)\simeq\frac{1}{\nu}\log\frac{s_{\nu}+S_{\nu}}{l}$ (19)

for $\nu$ large enough.

As a matter of fact, the numerical algorithm below estimates $h(f)$ by
$\frac{1}{\nu}\log\frac{s_{\nu}+S_{\nu}}{l}$ to the desired precision. The
core of the algorithm consists of a loop over $\nu$. Each time the algorithm
enters the loop, the values of $s_{\nu-1}$ and $S_{\nu-1}$ are updated to
$s_{\nu}$ and $S_{\nu}$, and the current estimation of $h(f)$ is compared to
the previous one. Note that the computation of $S_{\nu}^{i}$, $1\leq i\leq l$,
requires $s_{0}^{i}=1,s_{1}^{i},...,s_{\nu-1}^{i}$, see (10), while the
computation of $s_{\nu}^{i}$, $1\leq i\leq l$, requires
$s_{0}^{i},s_{1}^{i},...,s_{\nu-1}^{i}$, and $S_{\nu}^{i}$, see (14).

We summarize next the algorithm resulting from (18) in the following scheme
(“$A\longrightarrow B$” stands for “$B$ is computed by means of $A$”).

(A1)

Parameters: $l\geq 1$ (number of critical points), $\varepsilon>0$ (dynamic
halt criterion), and $n_{\max}\geq 2$ (maximum number of loops).

(A2)

Initialization: $s_{0}^{i}=1$, and $K_{1}^{i}=\\{k,1\leq k\leq
l:\omega_{1}^{k}\in\mathcal{B}^{i}\\}$ ($1\leq i\leq l$).

(A3)

First iteration: For $1\leq i\leq l$,

$\begin{array}[]{rcll}s_{0}^{i},K_{1}^{i}&\longrightarrow&S_{1}^{i},S_{1}&\text{(use
(\ref{notation0}), (\ref{notation}))}\\\
s_{0}^{i},S_{1}^{i}&\longrightarrow&s_{1}^{i},s_{1}&\text{(use (\ref{sinu}),
(\ref{account}))}\end{array}$

(A4)

Computation loop. For $1\leq i\leq l$ and $\nu\geq 2$ keep calculating
$K_{\nu}^{i}$, $S_{\nu}^{i}$, and $s_{\nu}^{i}$ according to the recursions

$\begin{array}[]{rcll}K_{\nu-1}^{i}&\longrightarrow&K_{\nu}^{i}&\text{(use
(\ref{Kni2}), (\ref{transition}))}\\\
s_{0}^{i},s_{1}^{i},...,s_{\nu-1}^{i},K_{\nu}^{i}&\longrightarrow&S_{\nu}^{i},S_{\nu}&\text{(use
(\ref{notation0}), (\ref{notation}))}\\\
s_{0}^{i},s_{1}^{i},...,s_{\nu-1}^{i},S_{\nu}^{i}&\longrightarrow&s_{\nu}^{i},s_{\nu}&\text{(use
(\ref{sinu}), (\ref{account}))}\end{array}$ (20)

until (i)

$\left|\frac{1}{\nu}\log\frac{s_{\nu}+S_{\nu}}{l}-\frac{1}{\nu-1}\log\frac{s_{\nu}+S_{\nu}}{l}\right|\leq\varepsilon,$
(21)

or, else, (ii) $\nu=n_{\max}+1$.

(A5)

Output. In case (i) output

$h(f)=\frac{1}{\nu}\log\frac{s_{\nu}+S_{\nu}}{l}.$ (22)

In case (ii) output “Algorithm failed”.

As said above, the algorithm (A1)-(A5) simplifies the original algorithm [5],
which formally consists of the same five steps above but is based on the exact
value of the lap number $\ell_{\nu}$. This entails that the new algorithm
needs more loops to output $h(f)$ with the same parameter $\varepsilon$ in the
halt criterion (21), although this does not necessarily mean that the overall
execution time will be longer since now less computations are required. In
fact, we will find both situations in the numerical simulations of Sect. 6.

Furthermore, given a halt criterion $\varepsilon$, the execution time depends
as well on the units (i.e., on the base of the logarithm), whichever algorithm
is used. For instance, if logarithms to base $e$ are used (i.e., $h(f)$ in
nats) and $\nu=n_{\text{nat}}$ is the first time that the halt criterion,

$\left|\ln\frac{s_{\nu}+S_{\nu}}{l}-\ln\frac{s_{\nu-1}+S_{\nu-1}}{l}\right|\leq\varepsilon,$

happens to hold in the computation loop, then

$\displaystyle\left|\log_{2}\frac{s_{n_{\text{nat}}}+S_{n_{\text{nat}}}}{l}-\log_{2}\frac{s_{n_{\text{nat}}-1}+S_{n_{\text{nat}}-1}}{l}\right|$
$\displaystyle=\frac{1}{\ln
2}\left|\ln\frac{s_{n_{\text{nat}}}+S_{n_{\text{nat}}}}{l}-\ln\frac{s_{n_{\text{nat}}-1}+S_{n_{\text{nat}}-1}}{l}\right|$
$\displaystyle\leq\frac{\varepsilon}{\ln 2}=1.4427\varepsilon.$

Therefore, if the the halt criterion

$\left|\log_{2}\frac{s_{\nu}+S_{\nu}}{l}-\log_{2}\frac{s_{\nu-1}+S_{\nu-1}}{l}\right|\leq\varepsilon$

for the computation of $h(f)$ in bits does not hold when $\nu=n_{\text{nat}}$,
i.e.,

$\left|\ln\frac{s_{n_{\text{nat}}}+S_{n_{\text{nat}}}}{l}-\ln\frac{s_{n_{\text{nat}}-1}+S_{n_{\text{nat}}-1}}{l}\right|>(\ln
2)\varepsilon=0.6932\varepsilon\text{,}$

then the algorithm will not exit the computation loop. We conclude that
$n_{\text{bit}}\geq n_{\text{nat}}$ with both algorithms, where
$n_{\text{bit}}$ is the exit loop when logarithms to base 2 are employed.

Two final remarks:

R1.

The parameter $\varepsilon$ does not bound the error
$\left|h(f)-\frac{1}{\nu}\log\frac{s_{\nu}+S_{\nu}}{l}\right|$ but the
difference between two consecutive estimations, see (21). The number of exact
decimal positions of $h(f)$ can be found out by taking different
$\varepsilon$’s , as we will see in the next section. Equivalently, one can
control how successive decimal positions of
$\frac{1}{\nu}\log\frac{s_{\nu}+S_{\nu}}{l}$ stabilize with growing $\nu$.
Moreover, the smaller $h(f)$, the smaller $\varepsilon$ has to be chosen to
achieve a given approximation precision.

R2.

According to [4, Thm. 4.2.4], $\frac{1}{\nu}\log\ell_{\nu}\geq h(f)$ for any
$\nu$. We may expect therefore that the numerical approximations (22) converge
from above to the true value of the topological entropy with ever more
iterations, in spite of the relation $\ell_{\nu}=\frac{1}{l}(s_{\nu}+S_{\nu})$
holding in general for boundary-anchored maps only.

##  6 Numerical simulations

In this section we compute the topological entropy of a variety of multimodal
maps. To this end, a code for arbitrary $l$ was written with PYTHON, and run
on an Intel(R) Core(TM)2 Duo CPU. All the numerical results will be given with
six decimal positions for brevity.

Thus, in Sect. 6.1 to 6.3 we calculate the entropy of families of uni-, bi-,
and trimodal maps, respectively, taken from [6] (unimodal case) and [5]
(general case). Except for particular values of the parameters, these maps are
not boundary-anchored. The purpose of our choice is to compare our entropy
plots with the plots published in those references. To complete the picture,
we will consider non-smooth maps in Sect. 6.4. The natural choice are
piecewise linear maps of constant slope because, in this case, the exact value
of the topological entropy is known. In all sections, we are going to compare
numerically the performance of the algorithm presented in Sect. 5 with the
general algorithm presented in [5, Sect. 7] by means of single maps. For
brevity we shall refer to them as the new algorithm and the old one,
respectively.

As for the units, the nat is the usual choice in Applied Mathematics and
Physics, while the bit is the standard unit in Information Theory and
Communication Technologies. In the following subsections we are actually going
to use both of them despite the fact that, as shown in Sect. 5, computations
with Napierian logarithms are faster to a given precision. To be specific, we
use bits in Sect. 6.2, and 6.3 for the sake of comparison with the results
published in [5], which are given in that unit.

###  6.1 Simulation with 1-modal maps

Let $\alpha>0$, $-1<\beta\leq 0$, and
$f_{\alpha,\beta}:[-(1+\beta),(1+\beta)]\rightarrow[-(1+\beta),(1+\beta)]$ be
defined as [6, Eqn. (29)]

$f_{\alpha,\beta}(x)=e^{-\alpha^{2}x^{2}}+\beta.$

These maps have the peculiarity of showing direct and reverse period-doubling
bifurcations when the parameters are monotonically changed [6, Fig. 3(a)].

Fig. 1 shows the plot of $h(f_{2.8,\beta})$ vs $\beta$ calculated with the
algorithm of Sect. 7. Here $\varepsilon=10^{-4}$ and the parameter $\beta$ was
increased in steps of $\Delta\beta=0.001$ from $\beta=-0.999$ to $\beta=0$.
Upon comparing Fig. 1 with Fig. 3(b) of [6], we see that both plots coincide
visually, except for the two vanishing entropy tails. We conclude that
$\varepsilon=10^{-4}$ is not small enough to obtain reliable estimations of
the topological entropy for vanishing values of $h(f_{2.8,\beta})$. This fact
can also be ascertained numerically by taking different values of
$\varepsilon$, as we do in the table below.

Figure 1: Plot of $h(f_{2.8,\beta})$ in nats vs $\beta$, $-1<\beta\leq 0$
($\varepsilon=10^{-4},\Delta\beta=0.001$).

To compare the convergence speed and execution time of the old and the new
algorithm, we have computed $h(f_{2.8,-0.5})$ with both algorithms for
different $\varepsilon$’s. The number of loops $n$ needed to achieve the halt
condition $\varepsilon=10^{-d}$, $4\leq d\leq 7$, and the execution time $t$
(in seconds) are listed in Table 1. The columns $h_{old}$, $n_{old}$, and
$t_{old}$ were obtained with the old algorithm, while the columns $h_{new}$,
$n_{new}$, and $t_{new}$ were obtained with the new one. For
$\varepsilon=10^{-4}$ it exceptionally holds $t_{old}>t_{new}$, otherwise
$t_{old}<t_{new}$. Furthermore, we conclude from Table 1 that $h(f)=0.52...$
nats with either algorithm and $\varepsilon=10^{-6}$, both decimal digits
being exact. If $\varepsilon=10^{-7}$ the old algorithm fixes the third
decimal digit, $h(f)=0.524...$ nats, whereas the new algorithm does not.

| $h_{old}$ | $n_{old}$ | $t_{old}$ | $h_{new}$ | $n_{new}$ | $t_{new}$
---|---|---|---|---|---|---
$\varepsilon=10^{-4}$ | 0.531968 | 81 | 0.031076 | 0.534106 | 101 | 0.021248
$\varepsilon=10^{-5}$ | 0.526645 | 253 | 0.179558 | 0.527305 | 318 | 0.193149
$\varepsilon=10^{-6}$ | 0.524935 | 797 | 1.684213 | 0.525142 | 1004 | 1.912784
$\varepsilon=10^{-7}$ | 0.524391 | 2519 | 16.369158 | 0.524456 | 3174 | 18.900032
Table 1: Comparison of performances when computing $h(f_{2.8,-0.5})$ in nats.

Fig. 2 depicts the values of $h(f_{\alpha,\beta})$ for $2\leq\alpha\leq 3$,
$-1<\beta\leq 0$, $\varepsilon=10^{-4}$, and $\Delta\alpha,\Delta\beta=0.01$.

Figure 2: Level sets of $h(f_{\alpha,\beta})$ in nats vs $\alpha,\beta$,
$2\leq\alpha\leq 3$, and $-1<\beta\leq 0$
($\varepsilon=10^{-4},\Delta\alpha=\Delta\beta=0.01$).

###  6.2 Simulation with 2-modal maps

Let $0\leq v_{2}<v_{1}\leq 1$ and $f_{v_{1},v_{2}}:[0,1]\rightarrow[0,1]$ be
defined as [5, Sect. 8.1]

$f_{v_{1},v_{2}}(x)=(v_{1}-v_{2})(16x^{3}-24x^{2}+9x)+v_{2},$

These maps have convenient properties for numerical simulations as they share
the same fixed critical points,

$c_{1}=1/4,\;c_{2}=3/4,$

the critical values are precisely the parameters,

$f_{v_{1},v_{2}}(1/4)=v_{1},\;f_{v_{1},v_{2}}(3/4)=v_{2},$

and the values of $f$ at the endpoints are explicitly given by the parameters
as follows:

$f_{v_{1},v_{2}}(0)=v_{2},\;f_{v_{1},v_{2}}(1)=v_{1}.$

Fig. 3 shows the plot of $h(f_{1,v_{2}})$ vs $v_{2}$, $0\leq v_{2}<1$,
computed with the new algorithm, $\varepsilon=10^{-4}$, and $\Delta
v_{2}=0.001$. Again, this plot coincides visually with the same plot computed
with the old algorithm [5, Fig. 4] except for the vanishing entropy tail,
which indicates that $\varepsilon=10^{-4}$ is too large a value for obtaining
accurate estimates in that parametric region.

Figure 3: Plot of $h(f_{1,v_{2}})$ in bits vs $v_{2}$, $0\leq v_{2}\leq 1$
($\varepsilon=10^{-4},\Delta v_{2}=0.001$).

Table 2 displays the performance of the new algorithm as compared to the old
one when computing $h(f_{0.9,0.1})$. This time $t_{old}>t_{new}$ for
$\varepsilon=10^{-d}$, $4\leq d\leq 7$ (as in Table 1). Furthermore, we obtain
two correct decimal digits of the topological entropy,
$h(f_{0.9,0.1})=0.60...$ bits, with both algorithms and $\varepsilon=10^{-6}$.

| $h_{old}$ | $n_{old}$ | $t_{old}$ | $h_{new}$ | $n_{new}$ | $t_{new}$
---|---|---|---|---|---|---
$\varepsilon=10^{-4}$ | 0.619682 | 195 | 0.286922 | 0.622100 | 218 | 0.253133
$\varepsilon=10^{-5}$ | 0.606568 | 613 | 2.665108 | 0.607310 | 688 | 2.485049
$\varepsilon=10^{-6}$ | 0.602385 | 1938 | 26.238006 | 0.602622 | 2173 | 24.890648
$\varepsilon=10^{-7}$ | 0.601062 | 6125 | 271.074381 | 0.601137 | 6871 | 265.198039
Table 2: Comparison of performances when computing $h(f_{0.9,0.1})$ in bits.

Fig. 4 depicts the values of $h(f_{v_{1},v_{2}})$ for $0\leq v_{2}\leq
v_{1}-0.5$, $\varepsilon=10^{-4}$, and $\Delta v_{1},\Delta v_{2}=0.01$.

Figure 4: Level sets of $h(f_{v_{1},v_{2}})$ in bits vs $v_{1},v_{2}$, $0\leq
v_{2}\leq v_{1}-0.5$ ($\varepsilon=10^{-4},\Delta v_{1}=\Delta v_{2}=0.01$).

###  6.3 Simulation with 3-modal maps

Consider next the $3$-modal maps $f_{v_{2},v_{3}}:[0,1]\rightarrow[0,1]$
defined by the quartic polynomials [5, Sect. 8.2]

$\displaystyle f_{v_{2},v_{3}}(x)$
$\displaystyle=\frac{4\left[\left(2\sqrt{2}-1\right)v_{2}-2v_{3}\right]x}{2(2\sqrt{2}+1)v_{3}-7v_{2}}\left[4\left(1+2\sqrt{2}\right)(x-1)(1-2x)^{2}v_{3}\right.$
$\displaystyle+\left.\left(-56x^{3}+20\left(4+\sqrt{2}\right)x^{2}-\left(37+18\sqrt{2}\right)x+3\sqrt{2}+5\right)v_{2}\right],$

where $0\leq v_{2}<v_{3}\leq 1$. The critical points of $f_{v_{2},v_{3}}$ are

$c_{1}=\frac{-\sqrt{2}\text{$v_{2}$}-4\text{$v_{2}$}+12\sqrt{2}\text{$v_{3}$}-8\text{$v_{3}$}}{8\left(-7\text{$v_{2}$}+4\sqrt{2}\text{$v_{3}$}+2\text{$v_{3}$}\right)},\quad
c_{2}=1/2,\quad c_{3}=\frac{1}{4}(2+\sqrt{2}).$

Moreover this family verifies $f_{v_{2},v_{3}}(0)=0$,
$f_{v_{2},v_{3}}(c_{2})=v_{2}$, $f(c_{3})=v_{3}$, and

$f_{v_{2},v_{3}}(1)=\frac{4\left(5\sqrt{2}-8\right)\text{$v_{2}$}\left(\left(2\sqrt{2}-1\right)\text{$v_{2}$}-2\text{$v_{3}$}\right)}{-7\text{$v_{2}$}+4\sqrt{2}\text{$v_{3}$}+2\text{$v_{3}$}}.$

Fig. 5 shows the plot of $h(f_{v_{2},1})$ vs $v_{2}$, $0\leq v_{2}<1$,
computed with the new algorithm, $\varepsilon=10^{-4}$, and $\Delta
v_{2}=0.001$. Once more, this plot coincides visually with the same plot
computed with the old algorithm [5, Fig. 7 (left)] except for the vanishing
entropy tail, which again indicates that $\varepsilon=10^{-4}$ is too large a
value for obtaining accurate estimates in that parametric region.

Figure 5: Plot of $h(f_{v_{2},1})$ in bits vs $v_{2}$, $0\leq v_{2}<1$
($\varepsilon=10^{-4},\Delta v_{2}=0.001$).

Table 3 displays the performance of the new algorithm as compared to the old
one when computing $h(f_{0.7,1})$. Also this time $t_{old}>t_{new}$ for
$\varepsilon=10^{-d}$, $4\leq d\leq 7$ (as in Table 1 and 2). Furthermore, we
obtain two correct decimal digits of the topological entropy,
$h(f_{0.7,1})=0.69...$ bits, with both algorithms and $\varepsilon=10^{-6}$.

| $h_{old}$ | $n_{old}$ | $t_{old}$ | $h_{new}$ | $n_{new}$ | $t_{new}$
---|---|---|---|---|---|---
$\varepsilon=10^{-4}$ | 0.710321 | 162 | 0.422875 | 0.711709 | 177 | 0.351586
$\varepsilon=10^{-5}$ | 0.699339 | 511 | 4.028703 | 0.699793 | 557 | 3.448455
$\varepsilon=10^{-6}$ | 0.695855 | 1615 | 39.882275 | 0.696000 | 1759 | 34.853025
$\varepsilon=10^{-7}$ | 0.694752 | 5105 | 417.489291 | 0.694798 | 5561 | 368.900286
Table 3: Comparison of performances when computing $h(f_{0.7,1})$ in bits.

Fig. 6 depicts the values of $h(f_{v_{2},v_{3}})$ for $v_{2}+0.3\leq v_{3}\leq
1$, $\varepsilon=10^{-4}$, and $\Delta v_{2},\Delta v_{3}=0.01$.

Figure 6: Level sets of $h(f_{v_{2},v_{3}})$ in bits vs $v_{2},v_{3}$,
$v_{2}+0.3\leq v_{3}\leq 1$ ($\varepsilon=10^{-4},\Delta v_{1}=\Delta
v_{2}=0.01$).

###  6.4 Simulation with higher multimodal maps

Last but not least, we are going to compare the performance of the old and new
algorithms with the 4- and 5-modal maps of Fig. 7. These are piecewise linear
maps on $[0,1]$, with constant slope $s=\pm 1.5$, critical points

$c_{1}=\tfrac{3}{10},\;c_{2}=\tfrac{23}{60},\;c_{3}=\tfrac{7}{15},\;c_{4}=\tfrac{11}{20},$

and critical values

$f(c_{1})=f(c_{3})=0.450,\;f(c_{2})=f(c_{4})=0.325,$

in the $l=4$ case, while

$c_{1}=0.3,\;c_{2}=0.4,\;c_{3}=0.5,\;c_{4}=0.6,\;c_{5}=0.7,$

and

$f(c_{1})=f(c_{3})=f(c_{5})=0.45,\;f(c_{2})=f(c_{4})=0.30,$

in the $l=5$ case. By [4, Corollary 4.3.13],

$h(f)=\max\\{0,\ln\left|s\right|\\}=\ln 1.5=0.40547\text{ nats}$

in both cases. At variance with the previous examples in Sects. 6.1 to 6.3,
these two maps are non-smooth and boundary anchored.

Figure 7: A piecewise linear 4-modal map (left) and 5-modal map (right) with
constant slope $s=\pm 1.5$.

Table 4 summarizes the computational performance of both algorithms with the
$4$-modal map. As happened with the 2-, and 3-modal maps (Tables 2 and 3), the
new algorithm needs more computation loops but less execution time than the
old one for all $\varepsilon=10^{-d}$, $4\leq d\leq 7$.

| $h_{old}$ | $n_{old}$ | $t_{old}$ | $h_{new}$ | $n_{new}$ | $t_{new}$
---|---|---|---|---|---|---
$\varepsilon=10^{-4}$ | 0.421218 | 160 | 0.697776 | 0.422215 | 169 | 0.576800
$\varepsilon=10^{-5}$ | 0.410476 | 503 | 6.583444 | 0.410776 | 533 | 5.747668
$\varepsilon=10^{-6}$ | 0.407051 | 1589 | 65.236068 | 0.407147 | 1683 | 57.951979
$\varepsilon=10^{-7}$ | 0.405967 | 5021 | 678.706894 | 0.405997 | 5321 | 616.59469
Table 4: Comparison of performances when computing $h(f)$ in nats with the
$4$-modal map of Fig. 7 (left).

Likewise, Table 5 summarizes the computational performance of both algorithms
with the $5$-modal map. It is worth noting that now both algorithms need the
same number of loops for all halt criteria $\varepsilon$, and yet the new
algorithm is faster.

| $h_{old}$ | $n_{old}$ | $t_{old}$ | $h_{new}$ | $n_{new}$ | $t_{new}$
---|---|---|---|---|---|---
$\varepsilon=10^{-4}$ | 0.420542 | 152 | 0.848166 | 0.420542 | 152 | 0.644305
$\varepsilon=10^{-5}$ | 0.410239 | 480 | 8.231152 | 0.410239 | 480 | 6.501711
$\varepsilon=10^{-6}$ | 0.406978 | 1515 | 81.429872 | 0.406978 | 1515 | 65.307619
$\varepsilon=10^{-7}$ | 0.405944 | 4788 | 864.376277 | 0.405944 | 4788 | 695.24749
Table 5: Comparison of performances when computing $h(f)$ in nats with the
$5$-modal map of Fig. 7 (right).

As in the preceding simulations, we conclude from Table 4 and 5 that both
algorithms determine two correct decimal positions of the topological entropy
of the corresponding map, $h(f)=0.40...$ nats. But this time the halt
criterion $\varepsilon=10^{-6}$ does not suffice; here one has to set
$\varepsilon=10^{-7}$ to achieve the same precision.

A concluding observation. As anticipated in the remark R2 of Sect. 5 and
illustrated in the Tables 1-5, the values of $h_{new}$ converge from above
with ever more computation loops (or smaller values of the parameter
$\varepsilon$). This property follows for $h_{old}$ from [4, Thm. 4.2.4].

##  7 Conclusion

The main contributions of this paper are the following.

(i) In Theorem 1, we proved that the transition rules for min-max symbols (2)
and (3), which were derived in [5] for twice differentiable multimodal maps,
actually hold true for just continuous ones.

(ii) As a result of Theorem 1, we conclude that the validity of formula (15),
which was proved in [5, Theorem 5.3] for twice differentiable multimodal maps,
can be extended to continuous maps. For subsequent applications, only the
particularization of (15) to boundary-anchored maps (Theorem 2) is needed.

(iii) The results reviewed and proved in Sects. 2 and 3, leads to the closed
formula (18) for the topological entropy of multimodal maps. Previously we
proved in Theorem 3 that, although $\ell_{n}$ clearly depends on the boundary
conditions, the limit $h(f)=\lim_{n\rightarrow\infty}\frac{1}{n}\log\ell_{n}$
does not.

(iv) The numerical algorithm proposed in Sect. 5 for the computation of $h(f)$
amounts to a recursive scheme to approximate the limit in the closed formula
(18).

This algorithm is a simplification and, at the same time, a generalization of
the recursion scheme proposed in [5] for $h(f)$. Indeed, it is a
simplification because Eqn. (15) was used in [5] to compute the lap number
$\ell_{\nu}$, while the abridged expression (14) is used here. In other words,
the new algorithm does not track the orbits of the endpoints. And it is also a
generalization because we proved in Theorem 2 that (14) (and (15) for that
matter) holds not only for twice differentiable maps (as assumed in [5,
Theorem 5.3])) but also for just continuous ones. By the way, this point was
numerically checked in Sect. 6.4.

The performances of both algorithms, old and new, were compared in Sect. 6.1
to 6.4 using smooth and non-smooth $l$-modal maps with $1\leq l\leq 5$. In
view of the results summarized in Tables 1 to 5, the old algorithm performs
better in the unimodal case, while the opposite occurs in the other multimodal
cases.

##  8 Acknowledgements

We thank our referees for their constructive criticism. We are also grateful
to José S. Cánovas and María Muñoz Guillermo (Universidad Politécnica de
Cartagena, Spain) for clarifying discussions, and to Víctor Jiménez
(Universidad de Murcia, Spain) for the elegant proof in the Appendix. This
work was financially supported by the Spanish Ministerio de Economía y
Competitividad, grant MTM2012-31698.

## References

  * [1] R. Adler, A. Konheim and M. McAndrew.  Topological entropy.  Transactions of the American Mathematical Society, 114:309–319, 1965.
  * [2] Peter Walters.  An Introduction to Ergodic Theory, volume 79.  Springer Science and Business Media, 2000.
  * [3] M. Misiurewicz and W. Szlenk.  Entropy of piecewise monotone mappings.  Studia Mathematica, 67(1):45–63, 1980.
  * [4] L. Alsedà, J. Llibre and M. Misiurewicz.  Combinatorial Dynamics and Entropy in Dimension One.  World Scientific, Singapore, 2000.
  * [5] J.M. Amigó, R. Dilão and A. Giménez.  Computing the topological entropy of multimodal maps via min-max sequences.  Entropy, 14(4):742–768, 2012.
  * [6] R. Dilão and J.M. Amigó.  Computing the topological entropy of unimodal maps.  International Journal of Bifurcation and Chaos, 22(6):1250152, 2012\.
  * [7] J. Dias de Deus, R. Dilão and J. Taborda Duarte.  Topological entropy and approaches to chaos in dynamics of the interval.  Physics Letters A, 90:1–4, 1982.
  * [8] R. Dilão.  Maps of the interval, symbolic dynamics, topological entropy and periodic behavior (in Portuguese).  Ph.D. Thesis, Instituto Superior Técnico: Lisbon, 1985.
  * [9] J. Milnor and W. Thurston.  On iterated maps of the interval.  In James C. Alexander, editor, Dynamical Systems, pages 465–563, Berlin, Heidelberg, 1988. Springer Berlin Heidelberg.
  * [10] Welington De Melo and Sebastian Van Strien.  One-dimensional dynamics, volume 25.  Springer Science and Business Media, 2012.
  * [11] J.W. Milnor and C.P. Tresser.  On entropy and monotonicity for real cubic maps.  Communications in Mathematical Physics, 209:123–178, 1998.

##  9 APPENDIX

Let $g:X\rightarrow X$ be a continuous map of a compact Hausdorff space $X$
into itself. A point $x\in X$ is nonwandering with respect to the map $g$ if
for any neighborhood $U$ of $x$ there is an $n\geq 1$ (possibly depending on
$x$) such that $f^{n}(U)\cap U\neq\emptyset$. Fixed and periodic points are
examples of nonwandering points. The closed set of all nonwandering points of
$g$ is called its nonwandering set and denoted by $\Omega(g)$. According to
[4, Lemma 4.1.5],

$h(g)=h(\left.g\right|_{\Omega(g)}).$ (23)

Furthermore, if

$X=\bigcup\limits_{i=1}^{k}Y_{i}$

and all $Y_{i}$ are closed and $g$-invariant (i.e., $g(Y_{i})\subset Y_{i})$,
then [4, Lemma 4.1.10],

$h(g)=\max_{1\leq i\leq k}h(\left.g\right|_{Y_{i}}).$ (24)

To prove Theorem 3, suppose that $f$ is an $l$-modal selfmap of the compact
interval $I$ with positive shape (the proof for maps with negative shape is
analogous).

Set $I=[a,b]$, and $J=[a^{\prime},b^{\prime}]$ with $a^{\prime}\leq a<b\leq
b^{\prime}$. If $f(a)=a$, choose $a^{\prime}=a$; if $f(b)=a$ ($l$ odd) or
$f(b)=b$ ($l$ even), choose $b^{\prime}=b$. For definiteness, we suppose the
most general situation, namely, $a^{\prime}<a$ and $b<b^{\prime}$. Let
$F:J\rightarrow J$ be such that (i) $F$ is strictly increasing on
$[a^{\prime},a]$, (ii) $\left.F\right|_{[a,b]}=f$, and (iii) $F$ is strictly
decreasing ($l$ odd) or strictly increasing ($l$ even) on $[b,b^{\prime}]$. In
particular, $F$ may be taken piecewise linear on
$[a^{\prime},a]\cup[b,b^{\prime}]$. Thus, $F\in\mathcal{M}_{l}(J)$ has the
same critical points and values as $f$, has the same shape and is boundary-
anchored. Note that the shape enters in how $f$ is extended to $F$.

Moreover, it is easy to check that $\Omega(F)=\Omega(f)\cup C$, where $C$ is a
closed and $F$-invariant set that only contains fixed points. Thus,
$h(\left.F\right|_{C})=0$ and, according to (23) and (24),

$h(F)=h(\left.F\right|_{\Omega(F)})=\max\\{h(\left.F\right|_{\Omega(f)}),h(\left.F\right|_{C})\\}=h(\left.F\right|_{\Omega(f)})=h(\left.f\right|_{\Omega(f)})=h(f).\;\square$