\section{Introduction}
\qquad Chaos of the particle motion near the black hole has been studied for a long time \cite{dettmann1995chaos,Suzuki_199,letelier1997chaos,aless1999chaos}. For example, some authors have studied the stability of the particle geodesic motion near black holes \cite{cornish2003lyapunov,pradhan2013lyapunov,shahzadi2018circular, Dalui2018PresenceOH,dalui2019induction,dalui2019horizon,dalui2020near} and the quasinormal modes of black hole \cite{ching1994quasinormal,ching1995wave,nollert1996significance,horowitz1999quasinormal,cardoso2008geodesic,dolan2009expansion,pradhan2012stability,yang2012quasinormalmode,pradhan2012isco,pradhan2014circular} from the calculation of the Lyapunov exponent of particle motion orbits. In recent years, due to the proposal of AdS/CFT duality, a method for studying more complex quantum chaos by analyzing relatively simple black hole dynamics has been given, which has aroused more attention to chaos near black holes.
\\
\indent In 2015, Maldacena, Shenker and Stanford pointed out in quantum field theory that for a quantum system like a black hole, its Lyapunov exponent $\lambda$ (used to indicate the intensity of chaos) should have an upper bound \cite{Maldacena2016ABO}
\begin{eqnarray}
\lambda \leq 2\pi T,\label{chaos bound1}
\end{eqnarray}
where $T$ is the Hawking temperature of the black hole, and we take the natural unit $\hbar=1$. When we consider the black hole thermodynamics, the Hawking temperature $T$ and surface gravity $\kappa$ satisfy $2\pi T=\kappa$. So we know that for the black hole, the chaos bound \meqref{chaos bound1} has an equivalent form
\begin{eqnarray}
\lambda \leq \kappa. \label{chaos bound}
\end{eqnarray}

\indent To study the universality of this bound, Hashimoto and  Tanahashi calculated the Lyapunov exponent $\lambda$ of a test particle near the horizon to verify the establishment of chaos bound \cite{Hashimoto2017UniversalityIC}, where external forces were included. In \cite{Zhao2018StaticEO}, the authors specifically considered the situation when charged particle near a charged black hole maintained a static equilibrium on the horizon, and investigated the expansion near the horizon. They pointed out that when the higher-order expansion term is included, the Lyapunov exponent of the unstable particle at the equilibrium position $r_0$ near the outer horizon of the charged black hole $r_+$ will satisfy
\begin{eqnarray}
\lambda^2=\kappa^2+\gamma\left(r_0-r_+ \right)+\mathcal{O}\left(\left(r_0-r_+\right)^2\right),
\end{eqnarray}
where the $\gamma$ is a parameter which is decided by the black hole. As we show, when $\gamma$ is positive, the chaos bound \meqref{chaos bound} is violated. This means that by looking for an appropriate black hole solution, an example would show that the chaos bound can be violated near the horizon. Despite Schwarzschild black holes, RN black holes, and their corresponding AdS solutions have been studied extensively, and we still wish to deepen our understanding of chaos by studying some special black holes. A dyonic black hole solution with quasi-topological electromagnetism was obtained in \cite{Liu2019QuasitopologicalED}. We would like to study the geodesic behavior of this black hole solution and try to find some interesting phenomena.
\\
\indent Susskind proposed that when a particle radially falls into a black hole, the exponentially increasing of its momentum is consistent with the growth of the operator in the quantum chaotic system \cite{Susskind2018WhyDT}, then the connections between momentum and quantum complexity were constantly being explored \cite{susskind2019complexity,mohapatra2019sizemomentum,barbn2019momentumcomplexity,colangelo2020chaos}, where the quantum complexity is a very meaningful concept that can be related to many physics topics such as the black hole's action \cite{susskind2014addendum,brown2015complexity,brown2015},the shock wave geometry \cite{stanford2014complexity}, the accelerated expansion of the universe \cite{ge2017quantumy} and the partition function \cite{sun2019complexity}. Some authors conducted the research work about the increasing of momentum when a particle falls into black holes. For example, the momentum growth of particles when they fall into a black hole with the repulsive potential was discussed \cite{brown2018falling,ageev2018things,Zhou2018ParticleMA} and the case in an acoustic black hole was recently discussed in \cite{wang2019geometry}. Hereafter, we pay more attention to the behavior of particle's radial motion as it falls into a black hole in the radial direction and wish to connect it to the chaos bound.
\\
\indent The paper is organized as follows. In section \ref{A special dyonic black hole}, we provide some information about the black hole solution with quasi-topolodical electromagnetism \cite{Liu2019QuasitopologicalED}. In section \ref{Static equilibrium of charged particles}, we discuss the static equilibrium of charged particles near some special dyonic black holes to check whether the chaos bound $\lambda \leq \kappa$ will be violated. In section \ref{Geodesic motion of neutral particles}, we study the geodesic motion of neutral particles, including circular geodesic motion and radially falling. In section \ref{Conclusion and discussions}, we briefly summarize the paper. In appendix \ref{upper bound}, we consider some details about the ``maximal'' Lyapunov exponent $\lambda_s$ outside these black holes. In appendix \ref{Static equilibrium with only gravity}, we discussed a special case about static equilibrium with only gravity.
\section{A special dyonic black hole}\label{A special dyonic black hole}
\qquad Let us consider a black hole with quasi-topological electromagnetism proposed by \cite{Liu2019QuasitopologicalED}. In \cite{Liu2019QuasitopologicalED}, the authors considered a 4-dimensional gravitational theory including the pure cosmological constant $\Lambda_0$ and the minimum coupling electromagnetic interaction. Its Lagrangian is given \cite{Liu2019QuasitopologicalED}
\begin{eqnarray}
\mathcal{L}=\sqrt{-g}\left[R-2\Lambda_0 -\alpha_1 F^2-\alpha_2\left(\left(F^2\right)^2-2F^{(4)} \right) \right]\label{Dl},
\end{eqnarray}
where $\alpha_1$ and $\alpha_2$ are coupling constants and $F^2=F^{\mu \nu}F_{\mu \nu}$, $F^{(2)}=F_{\nu}^{mu}F_{\mu}^{nu}$, $F^{(4)}=F^{\mu}_{\nu}F^{\nu}_{\rho}F^{\rho}_{\sigma}F^{\sigma}_{\mu}$.
The corresponding Maxwell's field equation is
\begin{eqnarray}
&\nabla _{\mu}\tilde{F}_{\mu \nu}=0\label{DM},
\end{eqnarray}
where $\tilde{F}_{\mu \nu}=4\alpha_1 F^{\mu \nu}+8\alpha_2 \left( F^2 F_{\mu \nu}-2F^{\mu \rho}F^{\sigma}_{\ \rho}F_{\sigma}^{\ \nu}\right)$. The Einstein's field equation can be written as
\begin{eqnarray}
R_{\mu \nu}-\frac{1}{2}Rg_{\mu \nu}+\Lambda_0 g_{\mu \nu}=T_{\mu \nu}\label{DE}.
\end{eqnarray}
This  theory yields a static dyonic black hole solution
\begin{eqnarray}
ds^2=-f\left( r\right)dt^2+\frac{dr^2}{f\left( r\right)}+r^2d\Omega^2_{2,\epsilon}\label{metric},
\end{eqnarray}
where the metric $d\Omega^2_{2,\epsilon}$ corresponds to a two-dimensional hyperbolic, torus  and sphere respectively with $\epsilon$ takes values -1, 0, 1. The metric function can be expressed as \cite{Liu2019QuasitopologicalED}
\begin{eqnarray}
f\left(r \right)=-\frac{1}{3}\Lambda_0r^2+\epsilon -\frac{2M}{r}+\frac{\alpha_1p^2}{r^2}+\frac{q^2}{\alpha_1 r^2}\,_2F_1\left(\frac{1}{4},1;\frac{5}{4};-\frac{4p^2\alpha_2}{r^4\alpha_1} \right)\label{Dmf},
\end{eqnarray}
where $M$ is the mass of black holes, $\Lambda_0$ is the cosmological constant, $q$ is the electric charge, $p$ is the pressure from energy-momentum tensor and $\,_2F_1\left(\frac{1}{4},1;\frac{5}{4};-\frac{4p^2\alpha_2}{r^4\alpha_1} \right)$ is a hypergeometric function.
\\
\indent The electric and magnetic charges can be obtained from electromagnetic tensor $\tilde{F}_{\mu \nu}$ which is defined by Maxwell's field equation \meqref{DM}
\begin{eqnarray}
Q_e=\frac{1}{4\pi}\int \tilde{F}^{0r}=q,\qquad Q_m=\frac{1}{4\alpha_1 \pi}\int F=\frac{p}{\alpha_1}\label{charge}.
\end{eqnarray}
So the electric and magnetic potential functions are given by \cite{Liu2019QuasitopologicalED}
\begin{eqnarray}
\begin{aligned}
\Phi_e \left( r\right)&=\frac{q\,_2F_1\left(\frac{1}{4},1;\frac{5}{4};-\frac{4p^2\alpha_2}{r^4\alpha_1} \right)}{\alpha_1r},
\\
\Phi_m \left( r\right)&=-\frac{q^2\,_2F_1\left(\frac{1}{4},1;\frac{5}{4};-\frac{4p^2\alpha_2}{r^4\alpha_1}\right)}{4pr}+\frac{\alpha_1 q^2r^3}{4p\left(4\alpha_2p^2+\alpha_1r^4\right)}+\frac{\alpha_1^2p}{r}.
\end{aligned}\label{Dpf}
\end{eqnarray}
When the constant parameters ($M,p,q,\Lambda_0,\alpha_1 , \alpha_2$) take the proper value, different types of black hole solutions can be obtained.
\\
\indent Set $\alpha_1=1, \Lambda_0=0, \epsilon=1$, two types of dyonic black holes can be obtained as the constants $(\alpha_2, q, p)$  in the general solution \meqref{Dmf} take appropriate values. Note that in these black hole solutions, as $M$ takes different values, the black hole will also have different behaviors.
\subsubsection*{$\bullet$ Case 1: black holes with at most two horizons}
Consider a special condition of \meqref{Dmf} with the parameters
\begin{equation*}
(q^2,p^2,\alpha_2)=(\frac{5}{2},\frac{1}{2},2).
\end{equation*}
The corresponding metric function and electric potential function become
\begin{eqnarray}
\begin{aligned}
f\left(r\right)&=1-\frac{2M}{r}+\frac{1}{2r^2}\left(1+5_2F_1[\frac{1}{4},1;\frac{5}{4};-\frac{4}{r^4}]\right),
\\
\Phi_e \left(r\right)&=\frac{1}{2r}\left(1+5_2F_1[\frac{1}{4},1;\frac{5}{4};-\frac{4}{r^4}] \right).
\end{aligned}\label{case1}
\end{eqnarray}
It is one of the first type black hole solutions with a constant $M_0$ (its value is related to the value of $q$, $p$, $\alpha_2$ and the specific definition of $M_0$ can be found in \cite{Liu2019QuasitopologicalED}). When the black hole mass M takes different values, different situations can exist :
\begin{equation*}
\left\{
\begin{aligned}
M>M_0 \qquad &black\ holes\ with\ two\ horizons;
\\
M=M_0 \qquad &extreme\ black\ hole;
\\
M<M_0 \qquad &there\ is\ a\ naked\ singularity\ at\ r = 0.
\end{aligned}
\right.
\end{equation*}
For \meqref{case1}, the corresponding constant $M_0=1.6372384$. In the subsequent calculation process, for the first type of black hole solutions, we will consider an example with two horizons:
\\
\indent Specifically taking M = 1.996, the corresponding two horizons locate at
\begin{equation*}
r_-=0.28617, \qquad r_+=3.00002.
\end{equation*}
For this black hole, the metric function $f\left(r\right)$ is monotonically increasing outside the outer horizon and the increase rate $\frac{df\left(r\right)}{dr}$ is monotonically decreasing.
\subsubsection*{$\bullet$ Case 2: black holes with at most four horizons}
There is another special case for \meqref{Dmf} with parameters:
\begin{equation*}
(q^2,p^2,\alpha_2)=(\frac{20868}{443},\frac{396}{443},\frac{196249}{1584}).
\end{equation*}
The corresponding metric function and electric potential function can be given by
\begin{eqnarray}
\begin{aligned}
f(r)&=1-\frac{2M}{r}+\frac{12}{443r^2}\bigg(33+1739_2F_1[\frac{1}{4},1;\frac{5}{4};-\frac{443}{r^4}]\bigg),
\\
\Phi_e(r)&=\frac{12}{443r}\bigg(33+1739_2F_1[\frac{1}{4},1;\frac{5}{4};-\frac{443}{r^4}]\bigg).
\end{aligned}\label{case2}
\end{eqnarray}
The second type black hole solutions with three constants: $M_-, M_+, M_0$ (their values depend on the values of $q$, $p$, $\alpha_2$ and the specific definition of them can be found in \cite{Liu2019QuasitopologicalED}). Similarly, different situations are depending on the value of M:
\begin{equation*}
\left\{
\begin{aligned}
M<M_-  	  \quad &There\ is\ no\ horizon\ here\ but\ a\ naked\ singularity\ at\ r=0;
\\
M_-<M<M_+ \quad &black\ holes\ with\ two\ horizons\ and\ the \ equilibrium\ of\ Newton's\ potential;
\\
M_+<M<M_0 \quad &black\ holes\ with\ four\ horizons;
\\
M>M_0 \quad &black\ holes\ with\ two\ horizons,
\end{aligned}
\right.
\end{equation*}
where the corresponding $M _-, M _+$, and $M_0$ are given by
\begin{equation*}
M_-=6.6316,\quad M_+=6.7730,\quad M_0=6.9135.
\end{equation*}
In the subsequent calculations, we will take the values of M = 6.7, 6.8, 7.0 as concrete examples. The following are the corresponding horizon radii:
\begin{flushleft}
\textbf{Case 2-1: Black hole with Newtonian potential equilibrium}
\end{flushleft}
\begin{equation*}
M=6.7, \qquad \qquad r_-=0.66823,\quad r_+=1.54979.
\end{equation*}
 For this black hole, the metric function $f\left(r\right)$ is not monotonically increasing outside the outer horizon, and the Newton potential has equilibriums at $r=2.8480$ and $r=6.1421$. So it is possible for particles to maintain a static equilibrium at these positions without any external forces.
\begin{flushleft}
\textbf{Case 2-2: Black hole with cosmological horizons}
\end{flushleft}
\begin{equation*}
M=6.8,\qquad \qquad r_-=0.53736,\quad r_+=2.07702 ,\quad r_1=5.16615, \quad r_2=6.79477,
\end{equation*}
where $r_-$, $r_+$ are the two horizons of this black hole, and $r_1$, $r_2$ are the so-called cosmological horizons. The metric function $f\left(r\right)$ monotonically increases first and then monotonically decreases between the outer event horizon and the inner cosmology horizon. This black hole also has an equilibrium of the Newton potential at $r=3.06723$.
\begin{flushleft}
\textbf{Case 2-3: Black hole with two horizons and non-monotonic metric function increase rate}
\end{flushleft}
\begin{equation*}
M=7.0,\qquad \qquad r_-=0.40827,\quad r_+=8.34842.
\end{equation*}
 For this black hole, the metric function $f\left(r\right)$ is monotonically increasing outside the outer horizon, but the increase rate $\frac{df\left(r\right)}{dr}$ increases first and then decreases outside the outer horizon.
\section{Static equilibrium of charged particles outside the black hole}\label{Static equilibrium of charged particles}
\qquad In this section, we consider general discussion of charged particles outside the black hole and derive the Lyapunov exponent. We consider a 4-dimensional static spherically symmetric black hole
\begin{eqnarray}
ds^2=-f\left(r\right)dt^2+\frac{dr^2}{f\left(r\right)}+d\Omega^2 \label{smetric}.
\end{eqnarray}
At the horizon, the surface gravity $\kappa$ is
\begin{eqnarray}
\kappa=\left.-\frac{1}{\sqrt{g_{rr}}}\frac{d\sqrt{g_{tt}}}{dr}\right|_{horizon}=\left.\frac{1}{2}f^{'}\left(r \right)\right|_{horizon}\label{kappa},
\end{eqnarray}
where the prime `` $'$ '' means derivative with respect to $r$.
\\
\indent The Lagrangian of a particle which is near this black hole can be written as
\begin{eqnarray}
\mathcal{L}=-m\left(\sqrt{-g_{\mu \nu}\dot{x}^{\mu}\dot{x}^{\nu}}+V\left(r \right)\right)\label{sl1},
\end{eqnarray}
where $m$ is the mass of the particle  and $V\left(r \right)$ is the external potential in the radial direction. When $V^{'}\left(r \right)<0$, the potential $V\left(r \right)$ can provide the particle with a repulsive force away from the black hole, so that the particle may not fall into the black hole and maintain static equilibrium near the horizon. In this section, we take the static gauge $\tau =t$, so the dot `` . '' denotes derivative with respect to the proper time `` $t$ '' in this section.
\\
\indent When the particle maintains static equilibrium near the horizon, the particle's Lagrangian \meqref{sl1} describing its radial motion can be reduced to
\begin{eqnarray}
\mathcal{L}=-m\left(\sqrt{f\left(r \right)-\frac{\dot{r}}{f\left(r \right)}}+V\left(r \right)\right)\label{sl2}.
\end{eqnarray}
Note that when the particle maintains static equilibrium, we have $\dot{r}\ll 1$. After expanding the Lagrangian in terms of $\dot{r}$, we can obtain the effective Lagrangian
\begin{eqnarray}
\mathcal{L}_{eff}=\frac{\dot{r}^2}{2f\left(r \right)\sqrt{f\left(r \right)}} -V_{eff}\left(r \right),\qquad V_{eff}\left(r\right)=\sqrt{f\left(r \right)}+V\left(r \right)\label{effl1},
\end{eqnarray}
where $V_{eff}$ is effective potential and at the particle's equilibrium position $V_{eff}^{'}=0$. After expanding $V_{eff}$ at the particle's static equilibrium position $r=r_0$, we have the effective Lagrangian satisfying the relation
\begin{eqnarray}
\mathcal{L}_{eff}\sim \frac{\dot{r}^2}{2f\left(r \right)\sqrt{f\left(r \right)}}\left[\dot{r}^2+\lambda^2 \left(r-r_0 \right)^2 \right]\label{effl2},
\end{eqnarray}
where $\lambda$ is Lyapunov exponent, which obeys
\begin{eqnarray}
\lambda^2=-f\left(r \right)\sqrt{f\left(r \right)}\left[\left(\sqrt{f\left(r \right)} \right)^{''}+V^{''}\left(r \right) \right]\label{slambda}.
\end{eqnarray}
In the equilibrium position $r=r_0$, $V_{eff}^{'}$=0, if $V_{eff}$ has the maximum $V_{eff}^{''}<0$, the static equilibrium will be unstable and $\lambda^2>0$. If $V_{eff}$ has the minimum $V_{eff}^{''}>0$, the static equilibrium is stable and $\lambda^2<0$.
\\
\indent From the effective Lagrangian \meqref{effl2}, we can see that the particle should follow the equation of motion at the equilibrium position
\begin{eqnarray}
\ddot{r}=\lambda^2 \left(r-r_0 \right).
\end{eqnarray}
This equation of motion tells us the particle's trajectory in the radial direction, that is to say
\begin{eqnarray}
r=r_0+Ae^{\lambda t}+Be^{-\lambda t}.
\end{eqnarray}
If $\lambda$ is real, then this exponential increase in particle motion may imply the existence of chaos. Actually, when there is no perturbation in other directions, the position  $r=r_0$ (where $V_{eff}$ has the maximum) will be a separatrix of the phase space, which can provide a ``maximal'' Lyapunov exponent here \cite{levin2008homoclinic,perezgiz2008homoclinic,hackmann2009analytic,Hashimoto2017UniversalityIC}. When the external potential is strong enough, the position where $V_{eff}$ has the maximum can be found near the horizon.
\subsection{Chaos on the horizon}
\qquad Returning to the black hole solution \meqref{Dmf}, we calculate the Lyapunov exponent when a charged particle maintains its static equilibrium near the horizon of these black holes  and further verify the chaos bound $\lambda \leq \kappa$. Here, the external potential we are considering is provided by the electric field, and its form is
\begin{eqnarray}
V\left(r \right)=\frac{e}{m}\Phi_e\left(r \right),
\end{eqnarray}
where $e$ and $m$ are the charge and mass of the particle, $\Phi_e\left(r \right)$ is the electric potential function.
\\
\indent For the particle's static equilibrium at $r=r_0$, the Lyapulov exponent $\lambda$ satisfies the expression \cite{Zhao2018StaticEO}
\begin{eqnarray}
\lambda^2=\left.f\left(r \right)\sqrt{f\left(r \right)}\left(\frac{\Phi_e^{''}\left(r \right)}{\Phi_e^{'}\left(r \right)}\left(\sqrt{f\left(r \right)} \right)^{'}-\left(\sqrt{f\left(r \right)} \right)^{''} \right)\right|_{r=r_0}\label{selambda}.
\end{eqnarray}
When the static equilibrium position is infinitely close to the black hole outer horizon $r=r_+$, the Lyapulov exponent $\lambda$ should have an upper bound, that is, the chaos bound $\lambda \leq \kappa$. The content of this section is mainly to calculate the Lyapulov exponent $\lambda$ of the charged particle near the outer horizon of these black holes and to check whether they satisfy the chaos bound.
\\
\indent Some interesting phenomena would occur when we pay attention to the near-horizon behavior shown in \cite{Zhao2018StaticEO}. The metric function and potential function of the black hole can be expanded to the second-order on the black hole event horizon $r=r_+$:
\begin{eqnarray}
\begin{aligned}
f(r)&=f_1(r-r_+)+f_2(r-r_+)^2...
\\
\Phi_e (r)&=\Phi_{e0}+\Phi_{e1}(r-r_+)+\Phi_{e2}(r-r_+)^2...
\end{aligned}\label{efz}
\end{eqnarray}
Then, take \meqref{efz} in \meqref{kappa} and \meqref{selambda}
\begin{eqnarray}
\begin{aligned}
\kappa& =\frac{1}{2}f_1,
 \\
 \lambda^2&=\left.\frac{12f_1f_2(r-r_+)^2\Phi_{e2}+8f_2^2(r-r_+)^3\Phi_{e2}+f_1^2(\Phi_{e1}+6(r-r_+)\Phi_{e2})}{4(\Phi_{e1}+2(r-r_+)\Phi_{e2})}\right|_{r=r_0}.
\end{aligned}\label{nhlambda}
\end{eqnarray}
Expanding $\lambda^2$ in \meqref{nhlambda} at the horizon $r=r_+$,  we can obtain
\begin{eqnarray}
\lambda^2=\frac{f_1^2}{4}+\frac{f_1^2\Phi_{e2}}{\Phi_{e1}}(r_0-r_+)+\mathcal{O}((r_0-r_+)^2)\label{nhelambda1}.
\end{eqnarray}
So \meqref{nhelambda1} can be rewritten as \cite{Zhao2018StaticEO}
\begin{eqnarray}
 \lambda^2=\kappa^2 +\gamma(r_0-r_+)+\mathcal{O}((r_0-r_+)^2),
 \qquad
 \gamma =4\kappa^2\frac{\Phi_{e2}}{\Phi_{e1}}\label{nhelambda2}.
\end{eqnarray}
From \meqref{nhelambda2}, it can be seen that when $\gamma >0$, the chaos bound $\lambda \leq \kappa$ will be violated.
\\
\indent The main reason we consider this dyonic black hole as our research object is that it can provide some special black hole solutions. Since the calculation process is very complicated, here we only show the main results of the calculations in the tabular form.
\subsubsection*{$\bullet$ Case 1}
We consider the black hole solution given in \meqref{case1} and take it into  \meqref{kappa} and \meqref{selambda} to calculate the surface gravity of the particles at the black hole horizon and the Lyapulov exponent $\lambda$. Then, the parameter $\gamma$ was calculated through \meqref{nhelambda2}. The calculation results are shown in the following table.
\\
\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
Black hole & M&  Horizon &$\lambda^2$ &$\kappa^2$ &$\gamma$\\
\hline
Case 1 & 1.996&$r_-=0.28617 \quad r_+=3.00002$&0.01283&0.01283&-0.03156\\
\hline
\end{tabular}.
\caption{The parameters $\lambda^2$, $\kappa^2$ and $\gamma$ of \textbf{Case 1} at the outer horizon $r_+$. As we show, there are $\lambda = \kappa$ and $\gamma < 0$, which means the chaos bound is satisfied.}
\end{center}
\end{table}
\subsubsection*{$\bullet$ Case 2}
Similarly, we also calculate the black hole solution given in \meqref{case2}
\\
\\
\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
Black hole & M&Horizon&$\lambda^2$ &$\kappa^2$ &$\gamma$\\
\hline
\textbf{Case 2-1}&6.7&$r_-=0.66823,\quad  r_+=1.54979$&0.01469&0.01469&-0.01537\\
\hline
 \textbf{Case 2-2}&6.8&$r_-=0.53736,\quad r_+=2.07702$&0.00720&0.00720&0.00846\\
 & &$r_1=5.16615,\quad r_2=6.79477$& & & \\
\hline
\textbf{Case 2-3}&7.0&$r_-=0.40827,\quad  r_+=8.34842$&0.00049&0.00049&-0.00039
\\
\hline
\end{tabular}.
\caption{The parameters $\lambda^2$, $\kappa^2$ and $\gamma$ of three black hole solutions in \textbf{Case 2} at their outer horizon $r_+$. There is always $\lambda = \kappa$ at the outer horizon $r_+$ for ever cases. But as we see, $\gamma$ is positive for the black hole in \textbf{Case 2-2}, which means that there may be a violation of the chaos bound $\lambda \leq \kappa$.}
\end{center}
\end{table}

\indent The above is the result of all our calculations about $\lambda$, $\kappa$ and $\gamma$ for these dyonic black holes in \meqref{case1} and \meqref{case2}. These calculations show that $\lambda =\kappa$ is universally established on the black hole outside horizon, ignoring certain errors, which means that the chaos bound of $\lambda \leq \kappa$ is satisfied.
\\
\indent However, when considering the near-horizon behavior of the black hole and expanding the correlation function to the second order, we find the case of $\gamma>0$, which indicates that there is a violation of the chaos bound. In order to examine this situation, we need do more general study about the second-order expansion of these black holes.
\subsection{More discussion on the expansion}
\qquad To more clearly show the anomalous situation of $\gamma$ in the second-order expansion, we will discuss the two sets of metric functions and their corresponding potential functions above (\meqref{case1} and \meqref{case2}), and then plot $\gamma$ as a function of $M$.
\begin{figure}[h]
\centering
\subfigure[\scriptsize{\textbf{Case 1}}]{
\includegraphics[scale=0.25]{32case1.pdf}
}
\quad
\subfigure[\scriptsize{\textbf{Case 2}}]{
\includegraphics[scale=0.25]{32case2.pdf}
}
\caption{The parameter $\gamma$ as a function of the black hole mass $M$. (a) shows no causality violation for M in \textbf{Case 1}. In (b), $\gamma$ can be positive in the range $M \in (6.73442,6.9135)$, which implies that the chaos bound $\lambda \leq \kappa$ is violated in \textbf{Case 2}.}\label{gp}
\end{figure}
\\
\indent As shown in \mpref{gp}, in the range $M \in (6.73442,6.9135)$,   the chaos bound $\lambda \leq \kappa$ can be violated in the near horizon region.

\section{Geodesic motion of neutral particles}\label{Geodesic motion of neutral particles}
\qquad In this section, we study the Lyapunov exponent of the geodesic motion of neutral particles near the black hole,  including circular geodesic motion and radial falling. We extend the static equilibrium of charged particles in section \ref{Static equilibrium of charged particles} beyond the horizon, and take the ``maximal'' Lyapunov exponent as $\lambda_s$ which can be obtained from the static equilibrium. From \meqref {selambda}, we obtain the formula of $\lambda_s$ for charged black holes
\begin{eqnarray}
\lambda_s=\sqrt{f\left(r \right)\sqrt{f\left(r \right)}\left(\frac{\Phi_e^{''}\left(r \right)}{\Phi_e^{'}\left(r \right)}\left(\sqrt{f\left(r \right)} \right)^{'}-\left(\sqrt{f\left(r \right)} \right)^{''} \right)}.\label{ls}
\end{eqnarray}
\indent Arbitrary location outside the horizon where the effective potential $V_{eff}$ has a maximum, the Lyapunov exponent $\lambda_s$ obtained from static equilibrium is at its maximal value. When we consider $\lambda_s$ at the horizon, the value of $\lambda_s$ will return to the surface gravity $\kappa$. We propose a hypothesis here: the ``maximal'' Lyapunov exponent $\lambda_s$ is an inherent property determined by the nature of the black hole\footnote{Note that this is our assumption that \meqref{ls} is maximal at any position outside the black hole. In other words, we assume there is a renormalization group flow of $\lambda_s$ outside the black hole. It can flow to $\lambda_s = \kappa$ at the event horizon. We will compare \meqref{ls} with the Lyapunov exponents obtained from particle circular geodesic motion \cite{cardoso2008geodesic,pradhan2012stability,pradhan2014circular} and radial falling.}. It is worthwhile to explore the connection between the ``maximal'' Lyapunov exponent $\lambda_s$ with the Lyapunov exponent of particle geodesic motion.
\subsection{Circular geodesic motion of neutral particles}
\qquad A 4-dimensional spherically symmetric black hole is considered here. Its metric is
\begin{eqnarray}
ds^2=-f\left(r\right)dt^2+\frac{dr^2}{f\left(r\right)}+d\Omega^2_2\label{4dbh}.
\end{eqnarray}
When we focus on a neutral test particle moving on the equatorial plane of this black hole ($\theta =\frac{\pi}{2}, \dot{\theta}=0$), its Lagrangian can be written as
\begin{eqnarray}
\mathcal{L}=\frac{1}{2}\left[g_{tt}\dot{t}^2+g_{rr}\dot{r}^2+g_{\phi \phi}\dot{\phi}^2 \right],
\end{eqnarray}
where the dot `` . '' denotes derivative with respect to the proper time `` $\tau$ ''.
\\
\indent According to the generalized momentum expression $p_q=\frac{\partial L}{\partial \dot{q}}$, we can get its generalized momentum as
\begin{eqnarray}
\begin{aligned}
&p_t=-f\left(r\right)\dot{t}=-E=Const,
\\
&p_{\phi}=r^2\dot{\phi}=L=Const,
\\
&p_r=\frac{\dot{r}}{f\left(r\right)},
\end{aligned}\label{cgm}
\end{eqnarray}
where $E$ is the particle's energy which is relative to the infinity  and $L$ is the angular momentum of the particle. Using normalization with four velocity $g_{\mu \nu}u^{\mu}u^{\nu}=\eta$
\begin{eqnarray}
-E\dot{t}+L\dot{\phi}+\frac{\dot{r}^2}{f\left(r\right)}=\eta,
\end{eqnarray}
we can obtain
\begin{eqnarray}
\dot{r}^2=E^2-\left(\frac{L^2}{r^2}-\eta \right)f\left( r\right)\label{rd}.
\end{eqnarray}
Note that $\eta=\pm 1,0$ corresponds to space-like, time-like and null geodesics respectively.
\\
\indent For a trajectory $X_i(t=0)=X_{i0}$ in a $n$-dimensional autonomous system, its equation of motion satisfies
\begin{eqnarray}
\frac{dX_i}{dt}=H_i\left(X_j\right),
\end{eqnarray}
where $X_i$ is a vector of n-dimensional states, and $H_i\left(X_j\right)$ describes the non-linear evolution of the dynamical system. To analyze the stability of the trajectory $X_i(t=0)=X_{i0}$, we linearize the motion equation
\begin{eqnarray}
\frac{d\delta X_i\left(t\right)}{dt}=K_{ij}\delta X_j\left(t\right),
\end{eqnarray}
where $K_{ij}$ is the Jacobian matrix, and its formula is
\begin{eqnarray}
K_{ij}=\left.\frac{\partial{H_i}}{\partial{X_j}}\right|_{X_{i0}}.\label{Jm1}
\end{eqnarray}
The eigenvalue of $K_{ij}$ is the Lyapunov exponent which indicates the stability of the trajectory $X_i(t=0)=X_{i0}$.
\\
\indent For the circular geodesic motion of neutral particles, we have the Euler-Lagrange equation
\begin{eqnarray}
\frac{dp_q}{d\tau}=\frac{\partial \mathcal{L}}{\partial q}\label{elp}.
\end{eqnarray}
After setting the phase space variables $X_i\left(t\right)=\left(p_r,r\right)$ and considering the particles moving in a circular orbit with a radius of $r=r_0$, we have two equations
\begin{eqnarray}
\frac{dp_r}{d\tau}=\frac{\partial{\mathcal{L}}}{\partial r} \quad \textrm{and} \quad \frac{dr}{d\tau}=\frac{p_r}{g_{rr}}.
\end{eqnarray}
Then, \meqref{Jm1} can be rewritten as
\begin{eqnarray}
K_{ij}=\left.\left|
\begin{array}{cc}
0&\frac{d}{dr}\left(\frac{\partial \mathcal{L}}{\partial r}\right)\\
\frac{1}{g_{rr}}&0\\
\end{array}
\right|\right|_{r=r_0}.
\end{eqnarray}
The eigenvalues of this matrix can give the expression of the proper time Lyapunov exponent of circular geodesic motion $\lambda_p$, and we observe that $\lambda_p$ satisfies \cite{pradhan2012stability}
\begin{eqnarray}
\lambda_p^2=\left.\frac{1}{g_{rr}}\frac{d}{dr}\left(\frac{\partial \mathcal{L}}{\partial r}\right)\right|_{r=r_0}\label{cgp1}.
\end{eqnarray}
From the Lagrange's equation of geodesic motion
\begin{eqnarray}
\frac{d}{d\tau}\left(\frac{\partial \mathcal{L}}{\partial{\dot{r}}}\right)-\frac{\partial \mathcal{L}}{\partial{r}}=0,
\end{eqnarray}
and the formula
\begin{eqnarray}
\frac{d}{d\tau}\left(\frac{\partial \mathcal{L}}{\partial{\dot{r}}}\right)=\frac{d}{d\tau}\left(-g_{rr}\dot{r}\right)=-\dot{r}\frac{d}{dr}\left(g_{rr}\dot{r}\right)=-\frac{1}{2g_{rr}}\frac{d}{dr}\left(g_{rr}^2\dot{r}^2\right),
\end{eqnarray}
the expression of $\frac{\partial \mathcal{L}}{\partial{r}}$ can be written as
\begin{eqnarray}
\frac{\partial \mathcal{L}}{\partial{r}}=-\frac{1}{2g_{rr}}\frac{d}{dr}\left(g_{rr}^2\dot{r}^2\right).\label{pl}
\end{eqnarray}
 For the circular geodesic motion, we have the circular geodesic condition \cite{cgc}
\begin{eqnarray}
\dot{r}^2=\left(\dot{r}^2 \right)^{'}=0.
\end{eqnarray}
Under this condition, we take \meqref{pl} in \meqref{cgp1}, then the proper time Lyapunov exponent $\lambda_p$ in \meqref{cgp1} can be reduced to
\begin{eqnarray}
\lambda_p=\sqrt{\frac{\left(\dot{r}^2\right)^{''}}{2}}\label{cgp2}.
\end{eqnarray}
If we consider an alternative form of \meqref{elp}
\begin{eqnarray}
\frac{dP_q}{dt}=\frac{d\tau}{dt}\frac{\partial \mathcal{L}}{\partial q}\label{elc},
\end{eqnarray}
we can express the Lyapunov exponent $\lambda_c$ with coordinate time \cite{cardoso2008geodesic}
\begin{eqnarray}
\lambda_c=\sqrt{\frac{\left(\dot{r}^2\right)^{''}}{2\dot{t}^2}}\label{cgc}.
\end{eqnarray}
In \meqref{cgm}, $\dot{t}=\frac{E}{f\left(r\right)}$, so we obtain the relation from \meqref{cgp2} and \meqref{cgc}
\begin{eqnarray}
\lambda_c=\frac{f\left(r\right)}{E}\lambda_p.
\end{eqnarray}
Both $f\left(r\right)$ and $E$  are real, so the property of $\lambda_p$ and $\lambda_c$ is closely related. It is evident that the Lyapunov exponents $\lambda_p$ and $\lambda_c$ are the most important parameters directly verifying the stability of the motion. Only when $\lambda_p$ and $\lambda_c$ are both real, it  results in the unstable geodesic motion. In contrast, when one of the $\lambda_p$ and $\lambda_c$ is imaginary, the circular geodesic motion is stable; when $\lambda_p=0$ or $\lambda_c=0$, the circular geodesic motion is marginal, which means the circular geodesic motion can be easily broken.
\\
\indent The Lyapunov exponent is a measure of the deviation in time evolution of two adjacent trajectories in phase space, obviously it depends on the time coordinate. In the following, we will study the coordinate time Lyapunov exponent $\lambda_c$ of the neutral particle's geodesic motion.
\subsubsection{Circular time-like geodesic}
\qquad For time-like geodesic with $\eta=-1$
\begin{eqnarray}
\dot{r}^2=E^2-\left(\frac{L^2}{r^2}+1 \right)f\left( r\right)\label{tled}.
\end{eqnarray}
The circular geodesic condition $\dot{r}^2=\left(\dot{r}^2 \right)^{'}=0$ yields \cite{cardoso2008geodesic}
\begin{eqnarray}
E^2=\frac{2f^2\left(r\right)}{2f\left(r\right)-rf^{'} \left(r\right)}, \quad L^2=\frac{r^3f^{'}\left(r\right)}{2f\left(r\right)-rf^{'}\left(r\right)}.\label{tlEL}
\end{eqnarray}
After taking the metric functions of the four black holes given by \meqref{case1} and \meqref{case2} into \meqref{tlEL} and \meqref{tled}, the coordinate time Lyapulov exponent $\lambda_c$ of these circular time-like geodesic motion can be given by
\meqref{cgc}.
\\
\indent To have a clear picture, we compare the circular geodesic motion of RN black hole with the black hole solutions given by \meqref{case1} and \meqref{case2}. The metric function and the electrial function of RN black hole are
\begin{eqnarray}
\begin{aligned}
f\left(r\right)&=1-\frac{2M}{r}+\frac{Q^2}{r^2},
\\
\Phi_e\left(r\right)&=\frac{2Q}{r},
\end{aligned}\label{RN}
\end{eqnarray}
where $M$ is the mass of black hole, and $Q$ is the charge. For RN black hole, $\lambda_c^2$ and $\lambda_s^2$ can be calculated by \meqref{cgc} and \meqref{ls}
\begin{eqnarray}
\begin{aligned}
	\lambda_c^2&=\frac{-Mr^3+6M^2r^2-9MQ^2r+4Q^2}{r^6}
	\\
	\lambda_s^2&=\frac{M^2-Q^2}{r^4}
\end{aligned}
\end{eqnarray}
\indent We want to explore the relationship between the circular time-like geodesic motion's Lyapunov exponent $\lambda_c$ and ``maximal'' Lyapunov exponent $\lambda_s$ defined from a static equilibrium, so we plot figures of $\lambda_c^2$ and $\lambda_s^2$ in the region where the circular geodesic motion exists.
\\
\indent For the RN black hole, as we show in \mpref{cRN}, $\lambda_c^2$ decreases monotonically from the innermost circular orbit until it becomes 0 at the innermost stable circular orbit, the circular time-like geodesic motion can only exist beyond the innermost stable circular orbit. Compared with RN black hole, the black hole given by \textbf{Case 2-1} has the similar evolution trend of $\lambda_c^2$ in \mpref{c21c} with $r>12.4354$, but in the \mpref{c21b}, there is $\lambda_c^2 < 0$ corresponding to the existence of the stable circular time-like geodesic orbit in the range of $6.14209<r<6.44301$ where is closer to the horizon than in \mpref{c21c}. We guess that this special property of \textbf{Case 2-1} having a Newtonian potential equilibrium position allows it to have a stable circular time-like geodesic closer to horizon.
\begin{figure}[htp]
\centering
\subfigure[$2+\sqrt{3}<r<16$]{
\includegraphics[scale=0.24]{41rn1.pdf}
}
\quad
\subfigure[$8<r<36$]{
\includegraphics[scale=0.24]{41rn2.pdf}
}
\caption{$\lambda_c$ and $\lambda_s$ are defined by \meqref{cgc} and \meqref{ls}  as functions of radial coordinate $r$. For the RN black hole with $M=2$ and $Q=1$,  $\lambda_c^2$ decreases from the radius of the innermost circular orbit (the beginning of the coordinate axis). $\lambda_s^2$ decreases too, but $\lambda_s^2$ is always positive in contrast to $\lambda_c^2$. Near the innermost orbital radius, we have $\lambda_c^2 > \lambda_s^2$. Away from the black hole, $\lambda_c^2$ is negative, which means the circular geodesic motion is stable.}\label{cRN}
\end{figure}
\\
\indent In \cite{Hashimoto2017UniversalityIC}, the authors regarded the innermost circular geodesic motion as a equilibrium with a repulsive potential given by particle circular motion. But as shown in \mpref{cRN}, \mpref{cCase1} and \mpref{cCase2-3}, we can see that when the circular time-like geodesic motion approach the innermost circular orbit, $\lambda_s$ will not constraint $\lambda_c$ any more. The ``maximal'' Lyapunov exponent $\lambda_s$ seems to constraint $\lambda_c$ in \mpref{c21a} and \mpref{cCase2-2} with $\lambda_s=\lambda_c$ at the position where the Newton potential has a equilibrium. We think $\lambda_s$ restricts $\lambda_c$ in \mpref{c21a} and \mpref{cCase2-2} because the particle is close to the equilibrium position of the Newtonian potential where particle motion is slow. In the other word, we guess the main reason of $\lambda_c > \lambda_s$ around the innermost circular orbit in \mpref{cRN}, \mpref{cCase1} and \mpref{cCase2-3} may be when the particle approaches the innermost circular orbit its speed is close to the speed of light. There may be another reason is the movement in other directions makes Lyapunov exponent growth. To verify our conjecture, we discuss the particle's radial falling in Section \ref{RF}.
\begin{figure}[htp]
\centering
\includegraphics[scale=0.32]{41case1.pdf}
\caption{For the \textbf{Case 1}, the black hole defined by \meqref{case1} with $M=1.996$, $\lambda_c^2$ and $\lambda_s^2$ are analogous in the RN black hole.}\label{cCase1}
\end{figure}
\begin{figure}[htp]
\centering
\subfigure[\scriptsize{The range  $2.19387<r<2.84801$.}]{
\includegraphics[scale=0.16]{41case211.pdf}\label{c21a}
}
\quad
\subfigure[\scriptsize{The range  $6.14209<r<6.44301$. $\lambda_s$ is not shown, because of $\lambda_s^2<0$ in this range.}]{
\includegraphics[scale=0.16]{41case212.pdf}\label{c21b}
}
\quad
\subfigure[\scriptsize{The range  $r>12.4354$. $\lambda_s$ is not shown, because of $\lambda_s^2<0$ in this range.}]{
\includegraphics[scale=0.16]{41case213.pdf}\label{c21c}
}
\caption{For the \textbf{Case 2-1}, the black hole defined by \meqref{case2} with $M=6.7$, the circular time-like geodesic motion exists in three discontinuous range. As shown in (b), there is $\lambda_c^2 <0$ in the range $r \in (6.14209,6.44301)$ near the horizon, which means that in this range there can be a stable time-like circular geodesic motion. In (c), $\lambda_c^2$ has a RN-like behavior.}\label{cCase2-1}
\end{figure}
\begin{figure}[htp]
	\centering
	\includegraphics[scale=0.32]{41case22.pdf}
	\caption{For the \textbf{Case 2-2}, the black hole defined by \meqref{case2} with $M=6.8$, there is always $\lambda_s^2>\lambda_c^2$ in the range $2.68843<r<3.06723$.}\label{cCase2-2}
\includegraphics[scale=0.32]{41case23.pdf}
\caption{For the \textbf{Case 2-3}, the black hole defined by \meqref{case2} with $M=7$, the behavior of $\lambda_c^2$ and $\lambda_s^2$  is similar to that of RN black hole.}\label{cCase2-3}
\end{figure}

\subsubsection{Circular null geodesic and stable photon sphere}
\qquad To study the motion of the photon outside the black hole, we set $\eta =0$ in \meqref{rd}. Then, \meqref{rd} becomes
\begin{eqnarray}
\dot{r}^2=E^2-\frac{L^2}{r^2}f\left( r\right).\label{cnr}
\end{eqnarray}
From the circular geodesic condition $\dot{r}^2=\left(\dot{r}^2 \right)^{'}=0$, we have
\begin{subequations}
	\begin{align}
	2f\left(r_{cn}\right)&=r_{cn}f^{'}\left(r_{cn}\right),\label{cnc1}
	\\
	\frac{E^2}{L^2}&=\frac{f\left(r_{cn}\right)}{r_{cn}^2},\label{cnc2}
	\end{align}
\end{subequations}
where $r_{cn}$ is the radius of the circular null geodesic, which is defined by \meqref{cnc1}. After taking the black hole solutions we consider  \meqref{cnc1} to obtain $r_{cn}$, we substitute the radius $r_{cn}$, \meqref{cnc2} and \meqref{cnr} into \meqref{cgp2} and \meqref{cgc} to calculate the Lyapunov exponent. From the calculation results of the Lyapunov exponent as shown in Table.\ref{ps}, we can judge the stability of the circular null geodesics and verify the existence of stable photon spheres.
\\
\indent Using the properties of Lyapunov exponent $\lambda$, it can be tested that for \textbf{Case 2-1}, a stable photon sphere does exist at $r=6.44301$ which agrees with the result in \cite{Liu2019QuasitopologicalED}.
\begin{table}[h]
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
Black hole&Photon sphere radius&$\lambda^2$ &Stability\\
\hline
\textbf{Case 1}&r=4.72028&$\lambda^2 >0$&unstable\\
\hline
 &r=2.19387&$\lambda^2 >0$&unstable\\
\textbf{Case 2-1}&r=6.44301&$\lambda^2<0$&stable\\
 &r=12.4354&$\lambda^2 >0$&unstable\\
\hline
\textbf{Case 2-2}&r=2.68843&$\lambda^2 >0$&unstable
\\
\hline
\textbf{Case 2-3}&r=14.3266&$\lambda^2 >0$&unstable
\\
\hline
\end{tabular}.
\end{center}
\caption{Photon sphere radius and stability analysis. For the black hole of \textbf{Case 2-1}, there is $\lambda^2< 0$ at $r=6.44301$, which means the photon sphere is stable there.}\label{ps}
\end{table}

\subsection{Radial falling}\label{RF}
\qquad In this subsection, we study the geodesic motion of particles falling radially into the black hole. The metric of a 4-dimensional spherically black hole can be written as \meqref{4dbh}. When considering that the particle fall freely towards the black hole in the radial direction, \meqref{rd} can be reduced as
\begin{eqnarray}
\dot{r}^2=E^2+\eta f,
\end{eqnarray}
then we can obtain
\begin{eqnarray}
\frac{dr}{dt}=\sqrt{f^2-\eta \frac{f^3}{E^2}}.
\end{eqnarray}
When we expand around each point on the particle trajectory (denoted as $r_f$), there is
\begin{eqnarray}
\frac{dr}{dt}\sim \frac{2f^{'}+\frac{3\eta}{E^2}ff^{'}}{2\sqrt{1+\eta \frac{f}{E^2}}}\left(r-r_f \right)\label{fme}.
\end{eqnarray}
Similarly, we can get an exponential growth form of the coordinate $r$. Here, we discuss in two cases:
\subsubsection*{mass particle}
For the mass particle falling from infinity, there is time-like geodesic. Setting $\eta=-1, E=1$, we can obtain the Lyapunov exponent $\lambda_{mp}$
 from \meqref{fme}
 \begin{eqnarray}
 \lambda_{mp}= \frac{2f^{'}-3ff^{'}}{2\sqrt{1-f}}\label{lmp}.
 \end{eqnarray}
\subsubsection*{photon}
For the photon, there is null geodesic. With $\eta=0$, we have the Lyapunov exponent $\lambda_{ph}$ from \meqref{fme}
\begin{eqnarray}
\lambda_{ph}=f^{'}\label{lph}.
\end{eqnarray}
\indent We can see that there is $\lambda_{mp}\approx \lambda_{ph}=2\kappa$ near the horizon\footnote{Near the horizon $r_+$, we have $f=0$ for \meqref{fme}, then \meqref{fme} can be rewritten as $\frac{dr}{dt}\sim f^{'}\left(r-r_+\right)$. Because of the surface gravity $\kappa=\left.\frac{1}{2}f^{'}\right|_{r=r_+}$, we can obtain $r-r_+=e^{2\kappa t}$ at the horizon agreeing with \cite{Zhou2018ParticleMA,wang2019geometry}}. We study the relationship between these two exponents outside the horizon and the ``maximal'' Lyapunov exponent $\lambda_s$ ($\lambda_s$ can be calculated via \meqref{selambda}).
\\
\indent For simplicity, we first study the RN black hole in \meqref{RN}, then \meqref{lmp}, \meqref{lph} and \meqref{selambda} can be rewritten as
\begin{eqnarray}
\begin{aligned}
\lambda_{mp}&=\frac{\left(Q^2-Mr\right)\left(3Q^2-6Mr+r^2\right)}{r^4\sqrt{2Mr-Q^2}},
\\
\lambda_{ph}&=-\frac{2Q^2}{r^3}+\frac{2M}{r^2},
\\
\lambda_{s}&=\sqrt{\frac{M^2-Q^2}{r^4}},
\end{aligned}
\end{eqnarray}
\begin{figure}[htp]
\centering
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.25]{42rna.pdf}
}
\quad
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.25]{42rnb.pdf}
}
\caption{$\lambda_{mp}$, $\lambda_{ph}$ and $\lambda_s$ are the function of $r$. For the RN black hole with $M=2,Q=1$, we can see from (a) that there is $\lambda_{ph}>\lambda_{mp}>\lambda_s$ near the horizon, and $\lambda_{mp}^2$ approaches $\lambda_{ph}^2$. From (b), there is always $\lambda_{ph} \geq 2\lambda_s$.}\label{fRN}
\end{figure}
\\
\indent As shown in \mpref{fRN}, there are some interesting phenomena. For massive particles, in the region near the event horizon, the closer  to the horizon, the closer $\lambda_{mp}$ approaches  $\lambda_{ph}$. The reason maybe because  the speed of the massive particle approaches the speed of light. For  photons outside the event horizon, there is always
\begin{eqnarray}
\lambda_{ph} \geq 2\lambda_s\label{f2c}.
\end{eqnarray}
This concides with $\lambda_{OTOC} \geq 2\lambda_{chaos}$ \cite{Xu2020DoesSE}  similar as the formula about chaos. Considering the relationship between null geodesic and shock wave, we speculate that there may be some relationship between these two formulas, or even equivalent.
\\
\indent To study the property of \meqref{f2c}, we perform the same calculation for the black holes given by \meqref{case1} and \meqref{case2}. From \mpref{fcase1} and \mpref{fcase2-3}, we can see that $\lambda_{ph} \geq 2\lambda_s$. But we can see that $\lambda_{ph} \geq 2\lambda_s$ from \mpref{fcase2-1} and \mpref{fcase2-2}. We think the reason is that the metric function $f(r)$ is not monotonically increasing outside the horizon, and we could understand the potential relationship between $\lambda_{ph} \geq 2\lambda_s$ and $\lambda_{OTOC} \geq 2\lambda_{chaos}$ by studying more examples.
\begin{figure}[htp]
\centering
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.16]{42case1a.pdf}
}
\quad
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.16]{42case1b.pdf}
}
\quad
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.16]{42case1c.pdf}
}
\caption{For the black hole in \textbf{Case 1}, we can see from (a) that $\lambda_{mp}^2$ approaches $\lambda_{ph}^2$ near the horizon. From (b), there is always $\lambda_{ph} \geq 2\lambda_s$ outside the horizon. In(c), the metric function $f\left(r\right)$ increases monotonically.}\label{fcase1}
\end{figure}
\begin{figure}[htp]
\centering
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.16]{42case21a.pdf}
}
\quad
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.16]{42case21b.pdf}
}
\quad
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.16]{42case21c.pdf}
}
\caption{For the black hole in \textbf{Case 2-1}, $\lambda_{mp}^2$ approaches $\lambda_{ph}^2$ near the horizon in (a). But there is not always $\lambda_{ph} \geq 2\lambda_s$ as shown in (b). Then we observe from (c) that the metric function $f\left(r\right)$ doesn't increase monotonically outside the horizon.}\label{fcase2-1}
\end{figure}
\begin{figure}[htp]
\centering
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.16]{42case22a.pdf}
}
\quad
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.16]{42case22b.pdf}
}
\quad
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.16]{42case22c.pdf}
}
\caption{For the black hole in \textbf{Case 2-2}, there is not always $\lambda_{ph} \geq 2\lambda_s$ in (b). At the same time, the metric function $f\left(r\right)$ doesn't increase monotonically outside the horizon.}\label{fcase2-2}
\end{figure}
\begin{figure}[htp]
\centering
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.16]{42case23a.pdf}
}
\quad
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.16]{42case23b.pdf}
}
\quad
\subfigure[\scriptsize{}]{
\includegraphics[scale=0.16]{42case23c.pdf}
}
\caption{From the black hole in \textbf{Case 2-3}, we can see from (a) that $\lambda_{mp}^2$ approaches $\lambda_{ph}^2$ near the horizon. $\lambda_{ph} \geq 2\lambda_s$ always exists in (b). Finally, the metric function $f\left(r\right)$ increases monotonically outside the horizon.}\label{fcase2-3}
\end{figure}
\\
\\
\section{Conclusion and discussions}\label{Conclusion and discussions}
\qquad Chaos near black holes is an important subject in contemporary physics research. As far as we know, the nature of the chaos near the black hole will depend on the black hole. Therefore, studying as many different black holes as possible may help us understand chaos more clearly.
\\
\indent In this paper, we mainly perform relevant calculations on the black hole solutions given in \cite{Liu2019QuasitopologicalED}. First, we review the definition of $\lambda_s$ which is the ``maximal'' Lyapunov exponent obtained from static equilibrium \cite{Hashimoto2017UniversalityIC}. We calculate the Lyapunov exponent when the charged particle remains static equilibrium near the horizon of these black holes given in Section \ref{A special dyonic black hole}. The results show that the chaos bound does not seem to be universally satisfied. When considering the higher-order expansion term, for some black holes with specific parameters, the chaos bound is violated.
\\
\indent After studying the static equilibrium of the particles, we turn  to the geodesic motion of the particles. We study two types of geodesic motion: circular geodesic motion and radial falling. For circular geodesic motion, we obtain the coordinate time Lyapunov exponent $\lambda_c$ by the Jacobian matrix and there are some stable circular orbits. At the same time, we find that these Lyapunov exponents $\lambda_c$ of circular geodesic motion do not seem to be constrained by the ``maximal'' Lyapunov exponent $\lambda_s$ defined from static equilibrium. The reason is that the particle speed approaches the speed of light . For neutral particles, the circular geodesic motion on the horizon does not exist, because it will lead to superluminal prorogation of particles, so $\lambda_c$ cannot be defined on the horizon.
\\
\indent Then we study the particle's radial falling, and find that the Lyapunov exponent $\lambda_{mp}$ of the mass particle will approach the Lyapunov exponent $\lambda_{ph}$ of the photon in the region near the horizon. We speculate the reason may be that the falling speed of the mass particle is close to the speed of light. Note that the discussions are carried out in the geodesic dynamics. It would be interesting to include the backreaction of the particle motion to the background spacetime and study  many-body effect of the chaos behavior. We found that for RN-like black holes or black hole examples where the metric function $f(r)$ increases monotonically outside the horizon, the relation $\lambda_{ph} \geq 2\lambda_s$ is established.

\section*{Acknowledgement} We would like to thank Hong L$\rm\ddot{u}$, Shu Lin and Qing-Bing Wang for helpful discussions. The study was partially supported by NSFC,
China (grant No.11875184).\\

\section*{Appendix}
\begin{appendix}
\section{$V_{eff}$ and the ``maximal'' Lyapunov exponent $\lambda_s$}\label{upper bound}
\qquad In studying the static equilibrium of particles to obtain the ``maximal'' Lyapunov exponent $\lambda_s$, the most important point is to ensure that the static equilibrium position is just at the maximum value of the effective potential $V_{eff}$.
\\
\indent For the static equilibrium of charged particles, there is
\begin{eqnarray}
V_{eff}=\sqrt{f\left(r\right)}+\frac{e}{m}\Phi_e\left(r\right),
\end{eqnarray}
$V_{eff}^{'}=0$ can be used to determine the static equilibrium position, and on the the static equilibrium position,there is
\begin{eqnarray}
\frac{e}{m}=-\frac{\left(\sqrt{f\left(r\right)}\right)^{'}}{\Phi_e^{'}\left(r\right)}\label{hzb}.
\end{eqnarray}
When we consider \meqref{hzb}, we can obtain
\begin{eqnarray}
V_{eff}^{''}=\left(\sqrt{f\left(r\right)}\right)^{''}-\left(\sqrt{f\left(r\right)}\right)^{'}\frac{\Phi_e^{''}\left(r\right)}{\Phi_e^{'}\left(r\right)}\label{V2}.
\end{eqnarray}
Only when \meqref{V2} is negative, the effective potential $V_{eff}$ will have a maximum value at the static equilibrium position, which corresponds to the ``maximal'' Lyapunov exponent $\lambda_s$. As we shown in \mpref{V2}, for the black hole represented by \textbf{Case 2-1}, within a certain range, static equilibrium always only has a minimum value of effective potential, which will not be discussed in this article.
\begin{figure}[h]
\centering
\subfigure[\scriptsize{Case 1}]{
\includegraphics[scale=0.25]{acase1.pdf}
}
\quad
\subfigure[\scriptsize{Case 2-1}]{
\includegraphics[scale=0.25]{acase21.pdf}
}
\\
\subfigure[\scriptsize{Case 2-2}]{
\includegraphics[scale=0.25]{acase22.pdf}
}
\quad
\subfigure[\scriptsize{Case 2-3}]{
\includegraphics[scale=0.25]{acase23.pdf}
}
\caption{ $V_{eff}^{''}$ as a function of $r$ corresponds to \meqref{case1} and \meqref{case2}}\label{V2p}
\end{figure}
\section{Static equilibrium with only gravity}\label{Static equilibrium with only gravity}
\qquad Since the metric functions of most black holes are monotonically increasing outside the horizon, it is extremely difficult for us to study the situation where particles maintain static equilibrium only under the action of gravitational potential. Fortunately, \textbf{Case 2-1} has an equilibrium position that allows particle to maintain static equilibrium with only gravity.
\\
\indent When only gravitational potential is considered, the effective potential is
\begin{eqnarray}
V_{eff}\left(r \right)=\sqrt{f\left(r \right)}\label{gveff},
\end{eqnarray}
and \eqref{slambda} reduces to
\begin{eqnarray}
\lambda=\sqrt{-\frac{f\left(r \right) f^{''}\left(r\right)}{2}}\label{nl1}.
\end{eqnarray}
From the analysis of the effective potential curve, we can know that $V_{eff}$ has maximum and minimum values. We will study $\lambda$ at the maximum position of $V_{eff}$. Since this position is far from the black hole horizon, the ``maximal'' Lyapunov exponent $\lambda_s$ should not be given by surface gravity $\kappa$.
\\
\indent Near the black hole given by \textbf{Case 2-1}, $V_{eff}$ has a maximum at $r=2.84801$. At this position, from \meqref{nl1}, we can obtain the Lyapunov exponent of static equilibrium of neutral particles considering only radial direction
\begin{figure}[h]
\centering
\includegraphics[scale=0.35]{bveff.pdf}
\caption{When only the gravity is considered, for \textbf{Case 2-1}, $V_{eff}$ has two extreme points, namely the maximum point $r=2.84801$ and the minimum point $r=6.14209$.}
\end{figure}
\end{appendix}
\newpage