\section{Motivation}
\label{sec:intro}
Einstein's theory of gravity (ETG) has already been investigated at low energy scales \cite{1}. However, it is also expected that this theory requires modification when it comes to high energy scales comparable to Planck scale. A strong predications related to phenomenon of higher dimensional geometries were arised in string theory and brane cosmology, so it is worthwhile to generalized gravity in higher dimensions \cite{2,2a,3,3a,3b}. In this context Lovelock \cite{4} introduced another gravity theory whose corresponding Lagrangian also contains high curvature terms along with the familiar Einstein-Hilbert Lagrangian terms in higher dimensions. It is worth-mentioning that in this theory equations of motion do not involve higher order derivatives of metric functions than the second order derivatives, due to this reason, it appears as a ghost-free theory \cite{5}. In four dimensions, ETG is recovered from Lovelock gravity; whereas in five dimensions, the quadratic curvature term commonly known as the Gauss–Bonnet term in string theory \cite{5,5a} appears in the gravitational action. Several black hole solutions in Lovelock gravities were derived in the literature \cite{5b,5c,5d,5e,5f,5g,5h,5i}.

 Assumption of taking the nonlinear terms of the invariants constructed from Riemann tensor in the gravitational Lagrangian, it seems natural to consider the nonlinear terms in the matter's sector as well. Therefore, for choosing an electromagnetic field as a source, it is worthwhile to use the nonlinear electromagnetic action instead of the linear Maxwell action. Therefore, in this work, we mainly focused on nonlinearly charged black hole solutions of Lovelock field equations, for doing this we will used the idea of non-linear electrodynamics (NLED). In cosmology the problem of the existence of big bang singularity and universe inflation are still not solved. It is widely believed that these problems can be solved when models of NLED are used. It should be noted that gravitational and electromagnetic interactions were very strong during creation of the universe. For this reason the nonlinear effects are not only playing a crucial role in early time universe but also are very important in the study of black hole's central singularities. From the beginning of the last century NLED got alot of attentions and several models have also been proposed for it. The problem of singularity which comes out in Maxwell's solution to the field of a point charge can be solved with the use of these models. The most famous proposal in this context is the formulation made by Born and Infeld \cite{81}. The action associated to the Born-Infeld (BI) formulation is constructed from low energy effective action corresponding to superstring theory \cite{82,83}. Due to this nonlinear action, electric field gives analytical expression and remains finite at the center of a charged particle. Like Born-Infeld theory of nonlinear electromagnetic fields, there exists also many other models of NLED introduced in literature \cite{84,85,86,87,88}. It was concluded that for these models too, the corresponding electric field at the center of charged particle appears to be finite. The important thing in BI model is its reduction towards familiar Maxwell's electrodynamics in the approximation of weak field. The properties of NLED are also appeared to be very comprehensible when coupled to gravity, this is due to the fact of strong electromagnetic fields domination in the early time of the universe. NLED also explains the physics at the center of charged particles and hence, it can be implemented to handle the difficulties of the essential singularity at the center of gravitating object like black hole. The exact black hole solution in ETG was given \cite{89} by use of a NLED Lagrangian which defines a large class of non-linear theories including BI and Euler-Heisenberg' electrodynamics \cite{90} whose weak field limit also approaches to Maxwell's electrodynamics. Other models of NLED in this regard are also considered as a sources of gravity in \cite{91,92,93,94,95}. The solution obtained with the use of these NLED models are singularity-free black hole solutions which asymptotically behave as the Reisnner-Nordstr\"{o}m (RN) black hole \cite{96,97}. There are alot of nonlinearly charged black hole solutions in ETG and other modified gravities coupled to NLED have been found in literature \cite{98,99,100,101,102,103,104,105,106,107,108,109,110,111,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,127a,127b,127c,127d,127e,127f,127g,127h,127i}. Moreover, instead of dark energy, non-linear electromagnetic fields are introduced to explain the inflation period in initial times where universe was born \cite{128,129}. Some models of NLED were also used to depict accelerated expansion of the universe \cite{130,131,132,133,134,135}. In $\Lambda$-Cold Dark Matter model, it has been shown that the universe's expansion is driven by the cosmological constant $\Lambda$, where the role of cosmological constant has been played by the non-zero trace of the energy-momentum tensor associated to NLED \cite{136,137}. Hence, like the above mentioned nonlinear charged black holes obtained in the framework of NLED, in this work we are also going to study charged black holes of Lovelock gravity coupled to recently proposed double-Logarithmic model of NLED \cite{88}.

The paper is constructed as follows. In the next section, we provide a brief introduction to Lovelock gravity and its coupling with double-Logarithmic electrodynamics has been performed. A polynomial equation which generates a family of magnetized Lovelock black hole solutions has been constructed where the source of Lovelock gravitational field is double-Logarithmic electrodynamics. In this section the gravitational field equations are solved and different thermodynamic quantities associated to such a Lovelock black holes are calculated in terms of event horizon. In Sections 3 and 4, black holes of Einstein and Gauss-Bonnet gravity have been studied, respectively. Finally, we summarize and conclude the paper in Section 5.

\section{Magnetized Lovelock black holes}

The action function defining Lovelock gravity with matter sources in diverse dimensions is written in the form
\begin{equation}
I=I_L+I_M,
\label{1}
\end{equation}
where $I_M$ corresponds to the action function of matter contents in the spacetime i.e. the double-Logarithmic electromagnetic field. The Lovelock action $I_L$ is given by
\begin{equation}
I_L=\frac{1}{2}\int d^Dx\sqrt{-g}\bigg[\sum_{p=0}^{s}\frac{\alpha_p}{2^p}\delta^{\mu_1...\mu_{2p}}_{\nu_1...\nu_{2p}} R^{\nu_1\nu_2}_{\mu_1\mu_2}...R^{\nu_{2p-1}\nu_{2p}}_{\mu_{2p-1}\mu_{2p}}\bigg],
\label{2}
\end{equation}
where $R^{\alpha\beta}_{\mu\nu}$ are the curvature tensor components and $\delta^{\mu_1...\mu_{2p}}_{\nu_1...\nu_{2p}}$ represents the generalized Kronecker delta having order $2p$ and $s=(d-1)/2$ where $s$ is a positive number. The coefficients $\alpha_p$ in (\ref{1}) are arbitrary constants in which $\alpha_0=-2\Lambda$, where $\Lambda$ is the cosmological constant.

Varying action (\ref{1}) with respect to the metric tensor, $g_{\mu\nu}$, the equations corresponding to the gravitational field can be obtained as
\begin{equation}
\sum_{p=0}^{s}\frac{\alpha_p}{2^{p+1}}\delta^{\nu\lambda_1...\lambda_{2p}}_{\mu\rho_1...\rho_{2p}} R^{\rho_1\rho_2}_{\lambda_1\lambda_2}...R^{\rho_{2p-1}\rho_{2p}}_{\lambda_{2p-1}\lambda_{2p}}=T^{\nu}_{\mu},
\label{4}
\end{equation}
where $T^{\nu}_{\mu}$ are the components of energy-momentum tensor given by
\begin{equation}
T_{\mu\nu}=-\frac{2}{\sqrt{-g}}\frac{\delta I_M}{\delta g^{\mu\nu}}.
\label{5}
\end{equation}
Now, the static and spherically symmetric line element in $D$-dimensions can be expressed in a general form as
\begin{equation}
ds^2=-f(r)dt^2+\frac{dr^2}{f(r)}+r^2 (h_{ab}dx^adx^b),
\label{6}
\end{equation}
where $h_{ab}dx^adx^b$ is the metric on $d-2$ dimensional hyper-surface whose scalar curvature is some constant $k$ such that the values $k=1,0,-1$, represent spherical, flat and hyperbolic spaces, respectively.

Now the Lagrangian density required for the definition of double-Logarithmic electromagnetic field is given by
 \begin{equation}\begin{split}
L_{EM}&=\frac{1}{2\beta}\bigg[\bigg(1-\sqrt{-2\beta P}\bigg)\log{\bigg(1-\sqrt{-2\beta P}\bigg)}+\bigg(1+\sqrt{-2\beta P}\bigg)\\&\times\log{\bigg(1+\sqrt{-2\beta P}\bigg)}\bigg],        \label{7}\end{split}
\end{equation}
 where $P=F_{\mu\nu}F^{\mu\nu}=2\bigg(\textbf{B}^2-\textbf{E}^2\bigg)$, is the Maxwell's invariant, $\textbf{E}$ represents the electric field, $\textbf{B}$ is the magnetic field and $F_{\mu\nu}$ is the electromagnetic field tensor which can be defined in terms of gauge potential $A_{\mu}$ as $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$. Variation of Eq. (\ref{1}) with respect to gauge potential $A_{\mu}$ yields the nonlinear electromagnetic field equations as
  \begin{equation}\begin{split}
 \partial^{\mu}\bigg[\frac{\sqrt{-g}}{\sqrt{2\beta P}}\log{\bigg(\frac{1-\sqrt{-2\beta P}}{1+\sqrt{-2\beta P}}\bigg)}F_{\mu\nu}\bigg]=0.  \label{8}\end{split}
 \end{equation}
 The corresponding energy-momentum tensor can be obtained in this case as
 \begin{equation}\begin{split}
 T_{\mu\nu}&=\frac{1}{2\beta}\bigg[\bigg(1-\sqrt{-2\beta P}\bigg)\log{\bigg(1-\sqrt{-2\beta P}\bigg)}+\bigg(1+\sqrt{-2\beta P}\bigg)\\&\times\log{\bigg(1+\sqrt{-2\beta P}\bigg)}\bigg]g_{\mu\nu}-\frac{2F_{\mu\lambda}F^{\lambda}_{\nu}}{\sqrt{-2\beta P}}\log{\bigg(\frac{1-\sqrt{-2\beta P}}{1+\sqrt{-2\beta P}}\bigg)}.        \label{9}\end{split}
 \end{equation}
 In order to construct a pure magnetized black hole solution, let us considered the case of only magnetic field contributions, i.e., $\textbf{B}\neq0$ and $\textbf{E}=0$, then $\textbf{B}=\frac{Q}{r^{d-2}}$ where $Q$ is the magnetic charge and the Maxwell's invariant can be expressed as $P=\frac{2Q^2}{r^{2d-4}}$.
   Let us introduce the metric function $f(r)$ in line element (\ref{6}) as
   \begin{equation}
   f(r)=k-r^2\Psi(r),\label{10}
   \end{equation}
   where the function $\Psi(r)$ is the polynomial defined as
   \begin{equation}
   P[\Psi(r)]=\sum_{p=0}^{s}\overline{\alpha}_p\Psi^p(r),\label{11}
   \end{equation}
   whose coefficients $\overline{\alpha}_p$ are defined as
   \begin{eqnarray}\begin{split}
   \overline{\alpha}_0&=\frac{\alpha_0}{(d-1)(d-2)},              \overline{\alpha}_1=1, \\
   \overline{\alpha}_p&=\prod_{i=3}^{2p}(d-i)\alpha_p.  \label{12}\end{split}
   \end{eqnarray}
   The general expression for $\overline{\alpha}_p$ in (\ref{12}) holds only for $p>3$. Therefore, using the above supposition of pure magnetic field in the expression of energy-momentum tensor, the gravitational field equations yields the Lovelock Polynomial for the magnetized Lovelock solution in the form as follows:
    \begin{eqnarray}\begin{split}
    P[\Psi(r)]&=\frac{2m}{(d-2)\Sigma_{d-2}r^{d-1}}-\frac{1}{\beta (d-2)(d-1)}\log{\bigg(1+\frac{4\beta Q^2}{r^{2d-4}}\bigg)}-\frac{8Q^2d}{(d-3)(d-1)r^{2d-4}}\\&F_1\bigg[1,\frac{d-3}{2d-4},\frac{3d-7}{2d-4},-\frac{4\beta Q^2}{r^{2d-4}}\bigg]+\frac{4Q r^{2-d}}{\sqrt{\beta}(d-2)(d-3)}\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r^{d-2}}\bigg)}, \label{13}\end{split}
    \end{eqnarray}
    which is satisfied by function $\Psi(r)$ and hence by the metric function $f(r)$ as well. Note that $C=2m/\Sigma_{d-2}(d-2)$ plays a role of integral constant associated with the black hole mass where
    \begin{equation}
    \Sigma_{d-2}=\frac{2\pi^{\frac{(d-1)}{2}}}{\Gamma{\bigg(\frac{d-1}{2}\bigg)}},\label{14}
    \end{equation}
      gives the volume of $d-2$-dimensional hypersurface. Next we want to calculate thermodynamical quantities associated to Lovelock black holes generated by polynomial equation (\ref{13}). These thermodynamical quantities can be computed in terms of the event horizon $r_h$ which satisfies the equation $f(r_h)=0$. So using Eq. (\ref{10}), we can write
    \begin{equation}
    r_h^2=\frac{k}{\Psi(r_h)}.\label{15}
    \end{equation}
    Thus, the total mass in terms of horizon's radius $r_h$ is given by
    \begin{eqnarray}\begin{split}
    m&=\frac{\Sigma_{d-2}}{2}\bigg[\sum_{p=0}^{s}\frac{\overline{\alpha}_pk^p(d-2)}{r_h^{-(d-2p-1)}}+\frac{8Q^2d(d-2)}{(d-3)(d-1)r_h^{d-3}}F_1\bigg[1,\frac{d-3}{2d-4},\frac{3d-7}{2d-4},-\frac{4\beta Q^2}{r_h^{2d-4}}\bigg]\\&+\frac{r_h^{d-1}}{\beta(d-1)}\log{\bigg(1+\frac{4\beta Q^2}{r_h^{2d-4}}\bigg)}-\frac{4Qr_h}{\sqrt{\beta}(d-3)}\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r_h^{d-2}}\bigg)}\bigg]. \label{16}\end{split}
    \end{eqnarray}
    The Hawking temperature \cite{138} of Lovelock black holes governed by the polynomial equation (\ref{13}) is defined in terms of surface gravity $\kappa_s$ as $T_H=\kappa_s/2\pi$. Thus, using the definition of surface gravity yields the Hawking temperature
     \begin{eqnarray}\begin{split}
   T_{H}&=\frac{1}{4\pi W(r_h)}\bigg[\sum_{p=0}^{s}\frac{\overline{\alpha}_pk^p(d-2p-1)}{r_h^{2p+1}}+\frac{1}{(d-2)r_h\beta}\log{\bigg(1+\frac{4\beta Q^2}{r_h^{2d-4}}\bigg)}\\&-\frac{8Q^2(d^2-4d+1)}{(d-1)(d-3)r_h(r_h^{2d-4}+4\beta Q^2)}+\frac{4Q\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r_h^{d-2}}\bigg)}}{\sqrt{\beta}(d-2)(d-3)r_h}\bigg], \label{17}\end{split}
    \end{eqnarray}
    where $W(r_h)$ is defined by
    \begin{equation}
    W(r_h)=\sum_{p=0}^{s}\frac{p\overline{\alpha}_pk^{p-1}}{r_h^{2p}}.\label{18}
    \end{equation}
    The entropy of the Lovelock black hole can be easily obtained with the use of Wald's method \cite{139,140}, this formulation says that
    \begin{equation}
    S=-2\pi\oint d^{d-2}x\sqrt{\gamma}\frac{\partial L}{\partial R_{abcd}}\epsilon_{abcd}, \label{19}
    \end{equation}
    in which $L$ is the Lagrangian density, $\gamma$ represents the determinant of the induced metric sector on the horizon and $\epsilon_{ab}$ defines the binormal to the horizon. Thus, using Lovelock action and the line element (\ref{6}), the entropy can be expressed in the following form
    \begin{equation}
    S=\frac{(d-2)\Sigma_{d-2}}{4}\sum_{p=1}^{s}\frac{\overline{\alpha}pk^{p-1}r_h^{d-2p}}{(d-2p)}, \label{20}.
    \end{equation}

    The heat capacity at constant charge $Q$ can be given by the relation
    \begin{equation}
    C_Q=T_{H}(r_h)\frac{dS}{dT_H}|_{Q}. \label{21}
    \end{equation}
    Differentiation of Eq. (\ref{17}) with respect to $r_h$ gives
     \begin{equation}\begin{split}
    \frac{\partial T_H}{\partial r_h}&=\frac{1}{4\pi W(r_h)}\bigg(\Xi_1(r_h)-\sum_{p=0}^{s}\frac{(2p+1)\overline{\alpha_p}(d-2p-1)k^{p}}{r_h^{2p+2}}\bigg)\\&-\frac{dW/dr_h}{4\pi W^2(r_h)}\bigg(\sum_{p=0}^{s}\frac{\overline{\alpha_p}k^p(d-2p-1)}{r_h^{2p+1}}+\Xi_2(r_h)\bigg), \label{22} \end{split}
    \end{equation}
    \begin{equation}\begin{split}
    \frac{dW}{dr_h}=-\sum_{p=1}^{s}\frac{2p^2\overline{\alpha_p}k^{p-1}}{r_h^{2p+1}}, \label{23},\end{split}
    \end{equation}
    \begin{equation}\begin{split}
    \Xi_1(r_h)&=-\frac{8Q^2r_h^{d-4}}{(d-3)\bigg(r_h^{2d-4}+4\beta Q^2\bigg)}-\frac{4Q}{\sqrt{\beta}(d-2)(d-3)r_h^2}\arctan{\bigg(\frac{2\sqrt{\beta}Q}{r_h^{d-2}}\bigg)}\\&+\frac{8Q^2(d^2-4d+1)\bigg(4\beta Q^2+(2d-3)r_h^{2d-4}\bigg)}{(d-1)(d-3)\bigg(r_h^{2d-3}+4\beta Q^2r_h\bigg)^2}-\frac{8Q^2}{r_h^2(r_h^{2d-4}+4\beta Q^2)}\\&-\frac{1}{\beta (d-2)r_h^2}\log{\bigg(1+\frac{4\beta Q^2}{r_h^{2d-4}}\bigg)}, \label{24} \end{split}
    \end{equation}
    \begin{equation}\begin{split}
    \Xi_2(r_h)&=\frac{1}{\beta (d-2)}\log{\bigg(1+\frac{4\beta Q^2}{r_h^{2d-4}}\bigg)}-\frac{8Q^2(d^2-4d+1)}{(d-1)(d-3)(r_h^{2d-3}+4\beta Q^2r_h)}\\&+\frac{4Q}{\sqrt{\beta}(d-2)(d-3)r_h}\arctan{\bigg(\frac{2\sqrt{\beta}Q}{r_h^{d-2}}\bigg)}. \label{25} \end{split}
    \end{equation}
    So, by using Eq. (\ref{20})-(\ref{25}), the heat capacity is expressed by the form
    \begin{eqnarray}\begin{split}
    C_Q=\frac{(d-2)\Sigma_{d-2}r_h^{d-1}W(r_h)\bigg(\Xi_2+\sum_{p=0}^{s}\frac{\overline{\alpha_p}k^p(d-2p-1)}{r_h^{2p+1}}\bigg)}{4\bigg(\Xi_1-\sum_{p=0}^{s}\frac{(2p+1)\overline{\alpha_p}(d-2p-1)k^p}{r_h^{2p+2}}-\frac{dW/dr_h}{W}\bigg(\sum_{p=0}^{s}\frac{\overline{\alpha_p}k^p(d-2p-1)}{r_h^{2p+1}}+\Xi_2\bigg)\bigg)}. \label{26}\end{split}
    \end{eqnarray}
    The expression of heat capacity is significant in the sense because of its role in thermodynamical stability of black holes. The region where it gives positive values indicates the stability while its negativity corresponds to the unstable thermodynamic system. The value of $r_h$ for which the sign of this quantity changes corresponds to the first order phase transition while for those values at which it becomes infinite and divergent defines the second order phase transition points of black hole.

\section{Magnetized Einsteinian black holes}

 The metric function associated to magnetized Einsteinian black holes in diverse dimensions is derived. This can be done by neglecting the higher curvature terms in the action function (\ref{1}), such that the coefficients $\alpha_p$ vanish for all $p\geq 2$. Note that, the source of gravitational field is again here the double-Logarithmic electromagnetic field. Therefore, by accepting this supposition, the metric function for $k=1$ from the polynomial equation (\ref{13}) can be determined as
      \begin{eqnarray}\begin{split}
      f(r)&=1-\frac{2m}{\Sigma_{d-2}(d-2)r^{d-3}}+\frac{\alpha_0r^2}{(d-1)(d-2)}+\frac{8dQ^2}{(d-1)(d-3)r^{2d-6}}\\&\times F_1\bigg[1,\frac{d-3}{2d-4},\frac{3d-7}{2d-4},\frac{-4\beta Q^2}{r^{2d-4}}\bigg]+\frac{r^2}{\beta (d-2)(d-1)}\log{\bigg(1+\frac{4\beta Q^2}{r^{2d-4}}\bigg)}\\&-\frac{4Q}{\sqrt{\beta}(d-2)(d-3)r^{d-4}}\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r^{d-2}}\bigg)}. \label{27}\end{split}
      \end{eqnarray}
      Choosing $d=4$, the above metric function becomes
      \begin{eqnarray}\begin{split}
      f(r)&=1-\frac{8\pi m}{r}-\frac{32Q^{2}}{3r^{2}}+\frac{\alpha_0r^2}{6}-\frac{4 Q^{\frac{3}{2}}}{3r\beta^{\frac{1}{4}}}\bigg[\arctan{\bigg(1+\frac{r}{\sqrt{Q}\beta^{\frac{1}{4}}}\bigg)}\\&-\arctan{\bigg(1-\frac{r}{\sqrt{Q}\beta^{\frac{1}{4}}}\bigg)}\bigg]+\frac{r^2}{6\beta}\log{\bigg(1+\frac{4\beta Q^2}{r^4}\bigg)}-\frac{2Q}{\sqrt{\beta}}\\&\times\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r^2}\bigg)}-\frac{2Q^{\frac{3}{2}}}{3\beta^{\frac{1}{4}}r}\log{\bigg[\frac{r^2-2\sqrt{Q}r\beta^{\frac{1}{4}}+2Q\sqrt{\beta}}{r^2+2\sqrt{Q}r\beta^{\frac{1}{4}}+2Q\sqrt{\beta}}\bigg]},\label{28}\end{split}
      \end{eqnarray}
which is now the metric function for non-asymptotically flat black hole of four-dimensional Einstein's gravity with double-Logarithmic electrodynamic source. If we choose $\beta\rightarrow 0$, then in this limit, the obtained metric function (\ref{27}) corresponds to the d-dimensional deSitter/antideSitter Reissner-Nordstrom like black hole with magnetic monopole charge for $\alpha_0$ positive or negative, respectively and is given by
\begin{eqnarray}\begin{split}
f(r)&=1-\frac{2m}{\Sigma_{d-2}(d-2)r^{d-3}}+\frac{\alpha_0r^2}{d^2-3d+2}+\frac{8Q^{2}d}{(d^2-4d+3)r^{2d-6}}+O(\beta).\label{29}\end{split}
\end{eqnarray}
Furthermore, Eq. (\ref{27}) shows that this solution has the behaviour of a non-asymptotically flat nonlinearly magnetically charged black hole , however, by choosing $\alpha_0=0$, the solution becomes an asymptotically flat.

In general, the Ricci and Kretschmann scalars for the line element (\ref{6}) are defined by the following expressions

\begin{eqnarray}\begin{split}
R(r)&=\bigg[(d-2)(d-3)\bigg(\frac{1-f(r)}{r^2}\bigg)-\frac{d^2f}{dr^2}-\frac{2(d-2)}{r}\frac{df}{dr}\bigg],\label{30}\end{split}
\end{eqnarray}
and
\begin{eqnarray}\begin{split}
K(r)&=\bigg[2(d-2)(d-3)\bigg(\frac{1-f(r)}{r^2}\bigg)^2-\bigg(\frac{d^2f}{dr^2}\bigg)^2+\frac{2(d-2)}{r^2}\bigg(\frac{df}{dr}\bigg)^2\bigg].\label{31}\end{split}
\end{eqnarray}
So, by using the metric function (\ref{27}), it can be clearly verified that both Ricci and Kretschmann scalars diverge at the center $r=0$, which confirms the existence of a true curvature singularity at this point and hence the resulting gravitating object is a black hole.
In $d$ dimensions, the mass of the black hole dependent on the event horizon $r_+$ can be constructed in the following form
    \begin{eqnarray}\begin{split}
     m&=\frac{\Sigma_{d-2}}{2}\bigg[(d-2)r_h^{d-3}+\frac{\alpha_0r_h^{d-1}}{(d-1)}+\frac{8Q^2d(d-2)}{(d-1)(d-3)r_h^{d-3}}F_1\bigg[1,\frac{d-3}{2d-4},\frac{3d-7}{2d-4},-\frac{4\beta Q^2}{r_h^{2d-4}}\bigg]\\&+\frac{r_h^{d-1}}{\beta (d-1)}\log{\bigg(1+\frac{4\beta Q^2}{r_h^{2d-4}}\bigg)}-\frac{4Qr_h}{\sqrt{\beta}(d-3)}\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r_h^{d-2}}\bigg)}\bigg].\label{32}\end{split}
          \end{eqnarray}

     The plots of mass in terms of horizon radius $r_h$ in different dimensions are shown in Fig. (\ref{skr1}) which shows that in each case there exists value $r_c$ such that for any value of $r_h$ greater than this critical horizon radius, no event horizon exists. This value $r_c$ corresponds to extremal value $m_e$ such that gravitating object possesses event horizon for all values of $m$ satisfying $m>m_e$. Fig. (\ref{Askar1a}) shows the plots of the metric function (\ref{27}) in different number of dimensions $d$. The point which corresponds to the intersection of associated curve with horizontal axis represents the location of the black hole's horizon for the chosen values of parameters involved in metric function.
   Using (\ref{27}), we obtain the Hawking temperature $T=\kappa/2\pi$ in this case as
    \begin{eqnarray}\begin{split}
    T_H(r_+)&=\frac{1}{4\pi}\bigg[\frac{d-3}{r_h}+\frac{\alpha_0r_h(d-1)}{(d^2-3d+2)}+\frac{r_h}{\beta (d-2)}\log{\bigg(1+\frac{4\beta Q^2}{r_h^{2d-4}}\bigg)}\\&-\frac{4Qr_h^{3-d}}{(d^2-5d+6)\sqrt{\beta}}\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r_h^{d-2}}\bigg)}-\frac{8Q^2r_h(d^2-3d-2)}{(d-1)(d-3)(r_h^{2d-4}+4\beta Q^2)}\bigg]. \label{33}\end{split}
    \end{eqnarray}

     \begin{figure}[h]
    	\centering
    	\includegraphics[width=0.8\textwidth]{mass1.png}
    	\caption{Plots of mass $M$ (Eq. (\ref{32})) vs $r_+$ for fixed values of $Q=10$, $\beta=1$, $\Sigma_{D-2}=100$ and $\alpha_0=-3$.}\label{skr1}
    \end{figure}
    \begin{figure}[h]
    	\centering
    	\includegraphics[width=0.8\textwidth]{metric1.png}
    	\caption{Plots of function $f(r)$ (Eq. (\ref{27})) for fixed values of $m=10$, $Q=10$, $\beta=1$, $\Sigma_{D-2}=100$ and $\alpha_0=-3$ .} \label{Askar1a}
    \end{figure}

    \begin{figure}[h]
    	\centering
    	\includegraphics[width=0.8\textwidth]{temperature1.png}
    	\caption{Plots of temperature $T_H$ (Eq. (\ref{33})) vs $r_h$ for fixed values of $Q=10$, $\beta=1$, $\Sigma_{D-2}=100$ and $\alpha_0=-3$.}\label{Saif1a}
    \end{figure}

Fig. (\ref{Saif1a}) shows the behaviour of Hawking temperature of Einsteinian black holes for different values of $D$ in the indicated domain of $r_+$. Now, the entropy associated tp Einsteinian black holes in this case yields
\begin{equation}
S=2\pi\Sigma_{d-2}r_h^{d-2},\label{34}
\end{equation}
which clearly says that in case of Einsteinian black holes Hawking area law is satisfied.
Differentiation of Hawking temperature (\ref{33}) yields
     \begin{eqnarray}\begin{split}
     \frac{dT_H}{dr_h}&=\frac{1}{4\pi}\bigg[\frac{\alpha_0}{(d-2)}-\frac{d-3}{r_h^2}+\frac{1}{\beta (d-2)}\log{\bigg(1+\frac{4\beta Q^2}{r_h^{2d-4}}\bigg)}\\&+\frac{4Qr_h^{2-d}}{(d-2)}\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r_h^{d-2}}\bigg)}-\frac{8Q^2(2d^2-8d+4)}{(d-1)(d-3)(r_h^{2d-4}+4\beta Q^2)}\\&+\frac{8Q^2r_h^{2d-4}(2d-4)(d^2-3d-2)}{(d-1)(d-3)(r_h^{2d-4}+4\beta Q^2)^2}\bigg].\label{35}\end{split}
     \end{eqnarray}
     The heat capacity at constant magnetic charge can be found from the general formula (\ref{21}), by using Eqs. (\ref{33})-(\ref{35}), as
     \begin{eqnarray}\begin{split}
    C_Q&=\frac{2\pi (d-2)\Sigma_{d-2}r_h^{d-2}\bigg((d^2-5d+6)+\alpha_0r_h^2+(d-2)r_h\zeta_1\bigg)}{\alpha_0r_h^2-(d-3)(d-2)+(d-2)r_h^2\zeta_2}, \label{36}\end{split}
     \end{eqnarray}
     where
     \begin{eqnarray}\begin{split}
     \zeta_1(r_h)&=\frac{r_h}{\beta(d-2)}\log{\bigg(1+\frac{4\beta Q^2}{r_h^{2d-4}}\bigg)}-\frac{4Q}{(d^2-5d+6)\sqrt{\beta}r_h^{d-3}}\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r_h^{d-2}}\bigg)}\\&-\frac{8Q^2(d^2-3d-2)r_h}{(d-1)(d-3)(r_h^{2d-4}+4\beta Q^2)}, \label{37}\end{split}
     \end{eqnarray}
  and
  \begin{eqnarray}\begin{split}
  \zeta_2(r_h)&=\frac{1}{\beta(d-2)}\log{\bigg(1+\frac{4\beta Q^2}{r_h^{2d-4}}\bigg)}+\frac{4Q}{(d-2)r_h^{d-2}}\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r_h^{d-2}}\bigg)}\\&-\frac{8Q^2(2d^2-8d+4)}{(d-1)(d-3)(r_h^{2d-4}+4\beta Q^2)}+\frac{8Q^2(2d-4)(d^2-3d-2)r_h^{2d-4}}{(d^2-4d+3)(r_h^{2d-4}+4\beta Q^2)^2}. \label{38}\end{split}
  \end{eqnarray}
\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{capacity1.png}
	\caption{Plot of function $C_Q$ (Eq. (\ref{36})) vs $r_h$ for fixed values of $Q=10$, $\beta=1$, $\Sigma_{D-2}=100$ and $\alpha_0=-3$.}\label{khan1}
\end{figure}
The graph of heat capacity is plotted in Fig. (\ref{khan1}). Note that, the stability of black hole is described by the region where heat capacity give positive values. The point on which sign of heat capacity changes corresponds to first order phase transition and the point on which heat capacity becomes indefinite implies the second order phase transition of black holes. The second order phase transition points are also determined from the condition $\frac{dT_H}{dr_h}=0$, because satisfaction of this condition implies heat capacity is divergent. The plot of $dT_H/dr_h$ is shown in Fig. (\ref{change1}). The points at which the curve associated to this quantity crossed the horizontal axis are the points corresponds the phase transitions of second order.
\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{change.png}
	\caption{Plot of function $\frac{dT_H}{dr_h}$ (Eq. (\ref{35})) vs $r_h$ for fixed values of $Q=10$, $\beta=1$, $\Sigma_{D-2}=100$ and $\alpha_0=-3$.}\label{change1}
\end{figure}

\section{Magnetized Gauss-Bonnet black holes}

  Now we want to discuss magnetized black holes in Gauss-Bonnet gravity coupled to double-Logarithmic electrodynamics. For this we consider $\alpha_2$ to be non-zero and $\alpha_p=0$ for $p\geq 3$ in the action function (\ref{1}) and considered $k=1$ in the polynomial equation (\ref{13}) to determine the expression for the metric function in two branches as
  \begin{eqnarray}\begin{split}
  f_{\pm}(r)&=1+\frac{r^2}{2\overline{\alpha}_2}\bigg(1\pm\sqrt{H(r)}\bigg), \label{39}\end{split}
  \end{eqnarray}
  where $H(r)$ is given by
  \begin{eqnarray}\begin{split}
  H(r)&=1-\frac{4\alpha_0\overline{\alpha}_2}{(d^2-3d+2)}+\frac{8\overline{\alpha}_2m}{\Sigma_{d-2}(d-2)r^{d-1}}-\frac{32\overline{\alpha_2}dQ^2}{(d-3)(d-1)r^{2d-4}}\\&\times F_1\bigg[1,\frac{d-3}{2d-4},\frac{3d-7}{2d-4},\frac{-4\beta Q^2}{r^{2d-4}}\bigg]-\frac{4\overline{\alpha_2}}{\beta(d^2-3d+2)}\log{\bigg(1+\frac{4\beta Q^2}{r^{2d-4}}\bigg)}\\&+\frac{16\overline{\alpha_2}Q}{\sqrt{\beta}(d-2)(d-3)r^{d-2}}\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r^{d-2}}\bigg)}, \label{40}\end{split}
  \end{eqnarray}
  where $\overline{\alpha_2}=(d-3)(d-4)\alpha_2$.
  The asymptotic value of $f_{\pm}(r)$ when $r\rightarrow\infty$ can be given as
  \begin{eqnarray}\begin{split}
  f_{\pm}(r)&=1+\frac{r^2}{2\overline{\alpha}_2}\bigg(1\pm\sqrt{1-\frac{4\overline{\alpha}_2\alpha_0}{(d-1)(d-2)}}\bigg). \label{41}\end{split}
  \end{eqnarray}
  This expression shows that for any positive value of $\overline{\alpha}_2$ the function $f_-(r)$ represents anti-de Sitter spacetime for $\alpha_0<0$ and de Sitter for $\alpha_0>0$. However, for any value of $\overline{\alpha_2}$, the positive branch $f_+(r)$ represents anti-deSitter when $\alpha_0>0$ whereas it represents deSitter for the values of $\alpha_0<0$. Like the previous case of Einsteinian black holes, here also both the curvature invariants have singularity at $r=0$. Thus, we conclude that our resulting solution (\ref{39}) of Gauss-Bonnet massive gravity also describes black hole. In order to study thermal stability, we will first compute mass of the black hole in terms of horizon radius by the condition $f(r_h)=0$, so that

  \begin{eqnarray}\begin{split}
  m&= \frac{\Sigma_{d-2}(d-2)}{2}\bigg[\alpha_2(d-3)(d-4)r_h^{d-5}+r_h^{d-3}+\frac{\alpha_0r_h^{d-1}}{(d-1)(d-2)}+\frac{r_h^{d-1}}{\beta (d-1)(d-2)}\\&\log{\bigg(1+\frac{4\beta Q^2}{r_h^{2d-4}}\bigg)}+\frac{8Q^2d}{(d-1)(d-3)r_h^{d-3}}F_1\bigg(1,\frac{d-3}{2d-4},\frac{3d-7}{2d-4},\frac{-4\beta Q^2}{r_h^{2d-4}}\bigg)\\&-\frac{4Qr_h}{\sqrt{\beta}(d-3)(d-4)}\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r_h^{d-2}}\bigg)}\bigg].\label{42}\end{split}
  \end{eqnarray}

  Figs. (\ref{MGB1}) shows the graphs of the mass function in terms of the horizon radius in different dimensions. Each of the plot indicate that there are negative values in the range of $m$ for different values of $r_h$ which donot corresponds to the horizons of the black hole. In other words, there exists value $r_1$ in each case such that all values of $r_h>r_1$ there doesnot correspond any horizon of black hole due to negativity of black hole's mass. This critical value $r_1$ corresponds to lower bound $m_l$ of the mass function.  Moreover, there also exists an upper bound $m_u$ on mass associated to value $r_2$ below which the associated value of $r_h$ represents the black hole event horizon. Thus, we have an interval $r_2\leq r_h\leq r_1$, in which all associated values corresponds to positive values of mass for which black hole horizons exist. The metric function $f(r)$ is plotted in Fig. (\ref{Askar2b}) for fixed values of parameters in different dimensions; the point at which the corresponding curve intersects the horizontal axis gives location of the event horizon. Similarly, the Hawking temperature $T_H$ of the obtained Gauss-Bonnet black hole is given by the following expression
  \begin{eqnarray}\begin{split}
  4\pi T_H(r_h)&=-\frac{2}{r_h}+\frac{r_h^4(r_h^2+2\alpha_2(d-3)(d-4)^{-1}}{4\alpha_2(d-3)(d-4)}\bigg[\frac{4\alpha_2^2(d-3)^2(d-4)^2(d-1)}{r_h^5}\\&+\frac{4(d-3)(d-4)\alpha_0\alpha_2}{(d-2)r_h}+\frac{4\alpha_2(d-3)(d-4)}{\beta (d-2)r_h}\log{\bigg(1+\frac{4\beta Q^2}{r_h^{2d-4}}\bigg)}\\&-\frac{16Q\alpha_2(d-4)}{\sqrt{\beta}(d-2)r_h^{d-1}}\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r_h^{d-2}}\bigg)}-\frac{32(d-4)(d^2-3d-2)Q^2\alpha_2}{(d-1)(r_h^{2d-3}+4\beta Q^2r_h)}\bigg].\label{43}\end{split}
  \end{eqnarray}

  \begin{figure}[h]
  	\centering
  	\includegraphics[width=0.8\textwidth]{mass2.png}
  	\caption{Plot of function $m(r_h)$ (Eq. (\ref{42})) for fixed values of $Q=10$, $\beta=1$, $\Sigma_{D-2}=100$, $\alpha_0=-3$ and $\alpha_2=2$.}\label{MGB1}
  \end{figure}

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{metric2.png}
	\caption{Plot of function $f_+(r)$ (Eq. (\ref{39})) for fixed values of $m=100$, $Q=10$, $\beta=1$, $\Sigma_{D-2}=100$, $\alpha_0=-3$ and $\alpha_2=2$.}\label{Askar2b}
\end{figure}
\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{temperature2.png}
	\caption{Plot of function $T_H(r_h)$ (Eq. (\ref{43})) for fixed values of $Q=10$, $\beta=0.1$, $\Sigma_{D-2}=100$, $\alpha_0=-3$ and $\alpha_2=2$.}\label{Saif2a}
\end{figure}
  The behaviors of Hawking temperature for such a black holes in different values of $d$ are shown in Fig. (\ref{Saif2a}). The entropy corresponding to the black hole solution (\ref{39}) can be obtained using Wald's method \cite{139,140} as
  \begin{eqnarray}
  S=(2d-4)\pi\Sigma_{d-2}r_h^{d-4}\bigg(\frac{r_h^2}{d-2}+2\alpha_2(d-3)\bigg),\label{44}
  \end{eqnarray}
  from which one can clearly sees that in this case the area law is not satisfied (as is the case in any other Lovelock black holes). Differentiation of (\ref{43}) gives
  \begin{eqnarray}\begin{split}
  4\pi\frac{dT_H}{dr_+}&=\frac{2}{r_h^2}+\frac{(d-4)r_h^4Z_1(r_h)}{r_h^2+2(d-4)(d-3)\alpha_2}-\frac{132(d-4)^2(d-3)Q\alpha_2r_h^{6-d}}{\sqrt{\beta}(d-2)(r_h^2+2(d-4)(d-3)\alpha_2)^2}\\&\times \arctan{\bigg(\frac{2\sqrt{\beta}Q}{r_h^{d-2}}\bigg)}+\frac{8(d-4)r_h^5(3r_h^2+4(d-3)(d-4)\alpha_2)Z_2(r_h)}{(r_h^2+2(d-3)(d-4)\alpha_2)^2},\label{45}\end{split}
  \end{eqnarray}
  where
  \begin{eqnarray}\begin{split}
  Z_2(r_h)&=\frac{8(d^2-3d-2)Q^2}{(d-1)r_h(r_h^{2d-4}+4Q^2\beta)}+\frac{(d-4)(d-3)^2(d-1)\alpha_2}{r_h^5}-\frac{(d-3)}{(d-2)r_h^3}\\&\times\big(\alpha_0r_h^2+d^2-3d+2\big)-\frac{(d-3)}{(d-4)\beta r_h}\log{\bigg(1+\frac{4\beta Q^2}{r^{2d-4}}\bigg)},\label{46}\end{split}
  \end{eqnarray}
  \begin{eqnarray}\begin{split}
  Z_1(r_h)&=\frac{4(d-1)Q}{(d-2)\sqrt{\beta}r_h^6}\arctan{\bigg(\frac{2\sqrt{\beta}Q}{r_h^{d-2}}\bigg)}+\frac{16Q^2(1+d^3-5d^2+5d)}{(d-1)r_h^2(r_h^{2d-4}+4\beta Q^2)}\\&-\frac{64(d-2)(d^2-3d-2)Q^4\beta}{(d-1)r_h^2(r_h^{2d-4}+4\beta Q^2)^2}-\frac{(d-3)(\alpha_0r_h^2+3d^2-9d+6)}{(d-2)r_h^4}\\&-\frac{5(d-4)(d-3)^2(d-1)\alpha_2}{r_h^6}-\frac{(d-3)}{(d-2)\beta r_h^2}\log{\bigg(1+\frac{4\beta Q^2}{r_h^{2d-4}}\bigg)}.\label{47}\end{split}
  \end{eqnarray}
 Finally, substituting the quantities calculated above into (\ref{21}), the heat capacity of our Gauss-Bonnet black hole can be expressed in the following form
 \begin{eqnarray}\begin{split}
 C_Q&=\frac{2\pi(d-2)\big(-8\alpha_2(d-3)(d-4)(r_h^2+2\alpha_2(d-3)(d-4))+r_h^5\Delta_1(r_h)\big)}{4r_h\alpha_2(d-4)(d-3)(r_h^2+2\alpha_2(d-3)(d-4))^{-1}}\\&\times\frac{\Sigma_{d-2}(r_h^{d-3}+2\alpha_2(d-3)(d-4)r_h^{d-5})}{\big((r_h^2+2\alpha_2(d-3)(d-4))^2\Delta_2(r_h)+8(d-4)r_h^5(3r_h^2+4(d-4)(d-3)\alpha_2)Z_2(r_h)\big)},\label{48}\end{split}
 \end{eqnarray}
 where
 \begin{eqnarray}\begin{split}
 \Delta_1(r_h)&=\frac{4\alpha_2^2(d-3)^2(d-4)^2(d-1)}{r_h^5}+\frac{4(d-3)(d-4)\alpha_0\alpha_2}{(d-2)r_h}\\&+\frac{4\alpha_2(d-3)(d-4)}{\beta (d-2)r_h}\log{\bigg(1+\frac{4\beta Q^2}{r_h^{2d-4}}\bigg)}-\frac{16Q\alpha_2(d-4)}{\sqrt{\beta}(d-2)r_h^{d-1}}\\&\times\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r_h^{d-2}}\bigg)}-\frac{32(d-4)(d^2-3d-2)Q^2\alpha_2}{(d-1)(r_h^{2d-3}+4\beta Q^2r_h)},\label{49}\end{split}
 \end{eqnarray}
 and
 \begin{eqnarray}\begin{split}
 \Delta_2(r_h)&=\frac{(d-4)r_h^4Z_1(r_h)}{(r_h^2+2\alpha_2(d-3)(d-4))}-\frac{132Q\alpha_2(d-4)^2(d-3)r_h^{6-d}}{\sqrt{\beta}(d-2)(r_h^2+2\alpha_2(d-3)(d-4))^2}\\&\times\arctan{\bigg(\frac{2Q\sqrt{\beta}}{r_h^{d-2}}\bigg)}+\frac{2}{r_h^2}.\label{50}\end{split}
 \end{eqnarray}
\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{capacity2.png}
	\caption{Plot of function $C_Q$ (Eq. (\ref{48})) vs $r_h$ for fixed values of $Q=10$, $\beta=1$, $\Sigma_{D-2}=100$, $\alpha_0=-3$ and $\overline{\alpha}_=2$.}\label{khan2}
\end{figure}
\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\textwidth]{change1.png}
	\caption{Plot of function $C_Q$ (Eq. (\ref{46})) vs $r_h$ for fixed values of $Q=10$, $\beta=1$, $\Sigma_{D-2}=100$, $\alpha_0=-3$ and $\overline{\alpha}_=2$.}\label{khan3}
\end{figure}
The graph of heat capacity is given in Fig. (\ref{khan2}). The region where the negativity of heat capacity occurs describes the black hole instability in it. The phase transition points are also clearly visible in each case, i.e., the points at which heat capacity changes sign represents the first order phase transitions and the points at which it diverges represents the second order phase transitions. It is worthwhile to note that, the second order phase transition points can also be understood from Fig. (\ref{khan3}), since this type phase transitions are described by the points where the associated curve in each case of dimensionality parameter intersect the horizontal axis.
\section{Summary and conclusion}

In this paper, the general static and spherically-symmetric line element is supposed and magnetized black holes are discussed in Lovelock, Gauss-Bonnet and Einstein gravities in the presence of double-Logarithmic electodynamics. After the coupling of Lovelock gravity with matter contents the resulting gravitational field equations are solved. In this process a polynomial equation (\ref{13}) is constructed which can generates the metric functions for magnetized black holes in Lovelock gravity of any order. We considered only pure magnetic field for magnetized solutions, because the problem with the electric type field doesnot yields the metric function in terms of elementary functions. Thermodynamics of these gravitating objects are then analyzed and in this context different thermodynamical quantities associated with the polynomial equation in terms of the event horizon are calculated. After this the two special cases described by polynomial equation (\ref{13}) and metric functions are constructed in terms of parameters $Q$ and $\beta$ corresponding to magnetized Einsteinian and Gauss-Bonnet black holes. It is shown that the obtained solutions are non-asymptotically flat and there exists a true curvature singularities at $r=0$ in each case. In addition to this, thermodynamic properties are studied and the mathematical expressions for the quantities like mass, Hawking temperature and heat capacity are computed. Each of these quantities are also plotted and their behaviors are discussed. The regions where the heat capacity and Hawking temperature are positive (negative) are pointed out which implies the black hole's thermal stability (instability). Furthermore, it is also shown that thermal phase transitions are also possible for our obtained black hole solutions. The first order phase transition is associated with those value of $r_h$ at which the heat capacity change signs, while the second order phase transitions of black holes corresponds to the zeros of $dT/dr_+=0$, or to the points at which the heat capacity diverges.

 It should be noted that when $\beta\rightarrow 0$ in each case, the resulting metric functions give non-asymptotically $d$-dimensional black holes of Einstein, Gauss-Bonnet and Lovelock gravities with Maxwel's electromagnetic source. However, for choosing $Q$ equal to zero, the calculated metric functions describe neutral black holes.

It would be very interesting to discuss Hawking radiations, quasinormal modes and grey body factors for the black holes in Einstein or Lovelock gravties coupled to this model of nonlinear electrodynamics. Further, the investigation of the causal structure and causality conditions will also give useful insights into the black hole solutions obtained in this paper. In addition to this, one can also used this double-Logarithmic electrodynamics for the analysis of accelerated expansion of the universe.


\section{Introduction}
In 2015,  the  LHCb collaboration  observed  two pentaquark candidates $P_c(4380)$ and $P_c(4450)$ in the $J/\psi p$ mass spectrum     with the statistical significances larger  than $9\sigma$ in the $\Lambda_b^0\to J/\psi K^- p$ decays \cite{LHCb-4380}.
In 2019, the LHCb collaboration  observed a  narrow pentaquark candidate $P_c(4312)$  with a statistical significance of  $7.3\sigma$ in the $\Lambda_b^0\to J/\psi K^- p$ decays,  and proved that the $P_c(4450)$ consists  of two narrow overlapping peaks $P_c(4440)$ and $P_c(4457)$
  with  the statistical significance of  $5.4\sigma$ \cite{LHCb-Pc4312}.
  The Breit-Wigner  masses and widths are
\begin{flalign}
 &P_c(4312) : M = 4311.9\pm0.7^{+6.8}_{-0.6} \mbox{ MeV}\, , \, \Gamma = 9.8\pm2.7^{+ 3.7}_{- 4.5} \mbox{ MeV} \, , \nonumber \\
 & P_c(4440) : M = 4440.3\pm1.3^{+4.1}_{-4.7} \mbox{ MeV}\, , \, \Gamma = 20.6\pm4.9_{-10.1}^{+ 8.7} \mbox{ MeV} \, , \nonumber \\
 &P_c(4457) : M = 4457.3\pm0.6^{+4.1}_{-1.7} \mbox{ MeV} \, ,\, \Gamma = 6.4\pm2.0_{- 1.9}^{+ 5.7} \mbox{ MeV} \,   .
\end{flalign}
 Very recently, the LHCb collaboration reported an evidence of a hidden-charm pentaquark candidate $P_{cs}(4459)$ with the strangeness $S=-1$ in the $J/\psi \Lambda$ invariant mass spectrum with  a statistical significance of  $3.1\sigma$ in the $\Xi_b^- \to J/\psi K^- \Lambda$ decays  \cite{LHCb-Pcs4459,LHCb-Pcs4459-2012},
the Breit-Wigner mass and width are
\begin{flalign}
 &P_{cs}(4459) : M = 4458.8 \pm 2.9 {}^{+4.7}_{-1.1} \mbox{ MeV}\, , \, \Gamma = 17.3 \pm 6.5 {}^{+8.0}_{-5.7} \mbox{ MeV} \, ,
\end{flalign}
but the spin and parity have not been determined yet.
 There have been  several  possible  assignments for  those hidden-charm  $P_c$ states, such as the diquark-diquark-antiquark type pentaquark states \cite{di-di-anti-penta-1,di-di-anti-penta-2,di-di-anti-penta-3,di-di-anti-penta-4,di-di-anti-penta-5,Wang1508-EPJC,WangHuang-EPJC-1508-12,
 WangZG-EPJC-1509-12,WangZG-NPB-1512-32,WangZhang-APPB,WZG-penta-IJMPA,Pc4312-penta-1,Pc4312-penta-2,Pc4312-penta-3},  the diquark-triquark  type  pentaquark states \cite{di-tri-penta-1,di-tri-penta-2,di-tri-penta-3},
 the  molecule-like  pentaquark  states \cite{mole-penta-before-Wu-PRL,mole-penta-before-Wu-PRC,mole-penta-before-Zhu,mole-penta-1,mole-penta-2,mole-penta-3,mole-penta-4,mole-penta-5,mole-penta-6,mole-penta-7,
 mole-penta-8,mole-penta-9,mole-penta-10,mole-penta-11,mole-penta-12,mole-penta-13,Pc4312-mole-penta-1,Pc4312-mole-penta-2,Pc4312-mole-penta-3,Pc4312-mole-penta-4,Pc4312-mole-penta-5,Pc4312-mole-penta-6,Pc4312-mole-penta-7,Pc4312-mole-penta-8,
Pc4312-mole-penta-9,Pc4312-mole-penta-10,Pcs4459-mole-penta-1,Pcs4459-mole-penta-2,Pcs4459-mole-penta-3}, the hadro-charmonium states \cite{Pc4312-hadrocharmonium-1,Pc4312-hadrocharmonium-2,Pc4312-hadrocharmonium-3},
 the re-scattering effects \cite{rescattering-penta-1,rescattering-penta-2}, etc.

 In Refs.\cite{Wang1508-EPJC,WangHuang-EPJC-1508-12,WangZG-EPJC-1509-12,WangZG-NPB-1512-32,WangZhang-APPB}, we construct the (light)diquark-(heavy)diquark-antiquark type pentaquark currents to study the spin-parity $J^P={\frac{1}{2}}^\pm$, ${\frac{3}{2}}^\pm$, ${\frac{5}{2}}^\pm$   hidden-charm pentaquark states with the strangeness  $S=0,\,-1,\,-2,\,-3$ via the QCD sum rules in an systematic way   by taking into account the contributions of the vacuum condensates up to  dimension 10, and explore the possible assignments of the $P_c(4380)$ and $P_c(4450)$. Now we confront the predicted masses of the pentaquark states with the strangeness $S=-1$ to the LHCb data \cite{LHCb-Pcs4459,LHCb-Pcs4459-2012}.

 In the case of the axialvector-diquark-scalar-diquark-antiquark type pentaquark states with the strangeness  $S=-1$, the predicted masses $M=4.51\pm0.11\,\rm{GeV}$ and $4.51\pm0.12\,\rm{GeV}$ for the pentaquark states with the spin-parity $J^P={\frac{1}{2}}^-$ and ${\frac{3}{2}}^-$ respectively are compatible with the LHCb data $4458.8 \pm 2.9 {}^{+4.7}_{-1.1} \mbox{ MeV}$ within uncertainties \cite{WangZG-EPJC-1509-12,WangZG-NPB-1512-32}.

   In the case of the axialvector-diquark-axialvector-diquark-antiquark type pentaquark states with the strangeness  $S=-1$, the predicted masses $M=4.47\pm0.15\,\rm{GeV}$, $4.51\pm 0.12\,\rm{GeV}$  and $4.52\pm0.12\,\rm{GeV}$ for the pentaquark states with the spin-parity $J^P={\frac{1}{2}}^-$, ${\frac{3}{2}}^-$ and ${\frac{3}{2}}^-$ respectively are also compatible with the LHCb data $4458.8 \pm 2.9 {}^{+4.7}_{-1.1} \mbox{ MeV}$ within uncertainties \cite{WangZG-EPJC-1509-12,WangZG-NPB-1512-32}.

  In the case of the scalar-diquark-axialvector-diquark-antiquark type pentaquark states with the strangeness  $S=-1$, the predicted mass $M=4.49\pm0.04\,\rm{GeV}$ for the pentaquark state with the spin-parity $J^P={\frac{3}{2}}^-$ is also compatible with the LHCb data $4458.8 \pm 2.9 {}^{+4.7}_{-1.1} \mbox{ MeV}$ within uncertainties \cite{WangZhang-APPB}.

  There are six possible pentaquark configurations with the spin-parity $J^P={\frac{1}{2}}^-$ or ${\frac{3}{2}}^-$ for the $P_{cs}(4459)$ in the old calculations, where
 we took into account the contributions of the vacuum condensates up to dimension  $10$, and used the energy scale formula $\mu=\sqrt{M_{P}^2-(2{\mathbb{M}}_c)^2}$  with the old value ${\mathbb{M}}_c=1.80\,\rm{GeV}$  to choose the optimal energy scales of the QCD spectral densities \cite{WangZG-EPJC-1509-12,WangZG-NPB-1512-32,WangZhang-APPB}.

After the discovery of the $P_c(4312)$,  we updated the old analysis by carrying out the operator product expansion up to the   vacuum condensates of $13$ in a consistent way and choosing the new value of the effective $c$-quark mass ${\mathbb{M}}_c=1.82\,\rm{GeV}$ \cite{WZG-tetraquark-Mc},  and  restudied the ground state mass spectrum of the diquark-diquark-antiquark type $uudc\bar{c}$ pentaquark states with the QCD sum rules  \cite{WZG-penta-IJMPA}.
The  predicted pentaquark  masses are improved slightly or remarkably  when the vacuum condensates of dimensions $11$ and $13$ are included, the lowest pentaquark state is $[ud][uc]\bar{c}$ ($0$, $0$, $0$, $\frac{1}{2}$)  with the spin-parity $J^P={\frac{1}{2}}^-$, which happens to  coincide with  the $P_c(4312)$.
Furthermore, the predicted pentaquark masses  support  assigning the $P_c(4440)$ to be the hidden-charm pentaquark state with  $J^{P}={\frac{1}{2}}^-$, ${\frac{3}{2}}^-$ or ${\frac{5}{2}}^-$ and assigning the $P_c(4457)$ to be the hidden-charm pentaquark state with $J^{P}={\frac{1}{2}}^-$ or ${\frac{3}{2}}^-$ \cite{WZG-penta-IJMPA}.

Now we take into account the flavor $SU(3)$ mass-breaking effects. For the traditional charmed baryon states, the  flavor $SU(3)$ mass-breaking effects due to the $s$-quark are $\Delta\approx 190\,\rm{MeV}$ and $120\,\rm{MeV}$ for the flavor antitriplets and sextets respectively from the Particle Data Group \cite{PDG}. We can take the average mass-breaking effects as  $\Delta=m_s=150\,\rm{MeV}$ in the present work, and make a crude estimation  for the masses of the pentaquark states with the strangeness  $S=-1$, which are presented in Table \ref{mass-1508-et al-Pcs}. From the table, we can see that the estimated  mass $4.46\pm0.11\,\rm{GeV}$ for the $[ud][sc]\bar{c}$ ($0$, $0$, $0$, $\frac{1}{2}$) pentaquark state is in excellent agreement with the experimental data $ 4458.8 \pm 2.9 {}^{+4.7}_{-1.1} \mbox{ MeV}$ from the  LHCb collaboration  \cite{LHCb-Pcs4459,LHCb-Pcs4459-2012}.

In this article, we extend our previous works \cite{Wang1508-EPJC,WangHuang-EPJC-1508-12,
 WangZG-EPJC-1509-12,WangZG-NPB-1512-32,WangZhang-APPB,WZG-penta-IJMPA} to study the ground state mass  of the diquark-diquark-antiquark type $udsc\bar{c}$ pentaquark state with the QCD sum rules, and explore the possible assignment of the $P_{cs}(4459)$ as the $[ud][sc]\bar{c}$ ($0$, $0$, $0$, $\frac{1}{2}$) pentaquark state with the spin-parity $J^P={\frac{1}{2}}^-$.

\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}\hline\hline
$[qq^\prime][q^{\prime\prime}c]\bar{c}$ ($S_L$, $S_H$, $J_{LH}$, $J$) &New analysis \cite{WZG-penta-IJMPA}    &$u\,{\rm or}\,d \to s$  &Assignments        \\ \hline

$[ud][uc]\bar{c}$ ($0$, $0$, $0$, $\frac{1}{2}$)                      &$4.31\pm0.11$   &$4.46\pm0.11$           &$?\,P_{cs}(4459)$         \\

$[ud][uc]\bar{c}$ ($0$, $1$, $1$, $\frac{1}{2}$)                      &$4.45\pm0.11$   &$4.60\pm0.11$           &     \\

$[uu][dc]\bar{c}+2[ud][uc]\bar{c}$ ($1$, $0$, $1$, $\frac{1}{2}$)     &$4.46\pm0.11$   &$4.61\pm0.11$           &      \\

$[uu][dc]\bar{c}+2[ud][uc]\bar{c}$ ($1$, $1$, $0$, $\frac{1}{2}$)     &$4.34\pm0.14$   &$4.49\pm0.14$           &      \\

$[ud][uc]\bar{c}$ ($0$, $1$, $1$, $\frac{3}{2}$)                      &$4.39\pm0.11$   &$4.54\pm0.11$           &       \\

$[uu][dc]\bar{c}+2[ud][uc]\bar{c}$ ($1$, $0$, $1$, $\frac{3}{2}$)     &$4.47\pm0.11$   &$4.62\pm0.11$           &    \\

$[uu][dc]\bar{c}+2[ud][uc]\bar{c}$ ($1$, $1$, $2$, $\frac{3}{2}$)     &$4.61\pm0.11$   &$4.76\pm0.11$           &                      \\

$[uu][dc]\bar{c}+2[ud][uc]\bar{c}$ ($1$, $1$, $2$, $\frac{3}{2}$)     &$4.52\pm0.11$   &$4.67\pm0.11$           &                     \\

$[uu][dc]\bar{c}+2[ud][uc]\bar{c}$ ($1$, $1$, $2$, $\frac{5}{2}$)     &$4.39\pm0.11$   &$4.54\pm0.11$           &      \\

$[ud][uc]\bar{c}$ ($0$, $1$, $1$, $\frac{5}{2}$)                      &$4.39\pm0.11$   &$4.54\pm0.11$           &     \\ \hline\hline
\end{tabular}
\end{center}
\caption{ The masses (in unit of GeV) of the pentaquark states with the strangeness $S=-1$, where the $S_L$ and $S_H$ denote the spins of the light diquarks and heavy diquarks respectively, $\vec{J}_{LH}=\vec{S}_L+\vec{S}_H$, $\vec{J}=\vec{J}_{LH}+\vec{J}_{\bar{c}}$, the $\vec{J}_{\bar{c}}$ is the angular momentum of the $\bar{c}$-quark.  }\label{mass-1508-et al-Pcs}
\end{table}

 The article is arranged as follows:   we derive the QCD sum rules for the mass and pole residue of  the
$[ud][sc]\bar{c}$ ($0$, $0$, $0$, $\frac{1}{2}$) pentaquark state in Sect.2;  in Sect.3, we present the numerical results and discussions; and Sect.4 is reserved for our
conclusion.

\section{QCD sum rules for  the  $udsc\bar{c}$ pentaquark state}
Firstly, let us write down  the two-point correlation function $\Pi(p)$ in the QCD sum rules,
\begin{eqnarray}\label{CF-Pi-Pi-Pi}
\Pi(p)&=&i\int d^4x e^{ip \cdot x} \langle0|T\left\{J(x)\bar{J}(0)\right\}|0\rangle \, ,
\end{eqnarray}
where
\begin{eqnarray}
 J(x)&=&\varepsilon^{ila} \varepsilon^{ijk}\varepsilon^{lmn}  u^T_j(x) C\gamma_5 d_k(x)\,s^T_m(x) C\gamma_5 c_n(x)\,  C\bar{c}^{T}_{a}(x) \, ,
\end{eqnarray}
where the $i$, $j$, $k$, $l$, $m$, $n$ and $a$ are color indices, the $C$ is the charge conjugation matrix \cite{Wang1508-EPJC,WangHuang-EPJC-1508-12,WangZG-EPJC-1509-12,WangZG-NPB-1512-32}. We choose the scalar-diquark-scalar-diquark-antiquark type current to interpolate the  $[ud][sc]\bar{c}$ ($0$, $0$, $0$, $\frac{1}{2}$) pentaquark state, see Table \ref{mass-1508-et al-Pcs}.
The current $J(x)$ couples potentially to the hidden-charm pentaquark state with the negative parity,
\begin{eqnarray}\label{Coupling-N}
\langle 0| J (0)|P^{-}(p)\rangle &=&\lambda_{-} U^{-}(p,s) \, ,
\end{eqnarray}
furthermore, it also couples potentially to the hidden-charm pentaquark state with the positive  parity,
\begin{eqnarray}\label{Coupling-P}
\langle 0| J (0)|P^{+}(p)\rangle &=&\lambda_{+} i\gamma_5 U^{+}(p,s) \, ,
\end{eqnarray}
as   multiplying $i \gamma_{5}$ to the current  $J(x)$ changes its parity \cite{Wang1508-EPJC,WangHuang-EPJC-1508-12,WangZG-EPJC-1509-12,WangZG-NPB-1512-32,WZG-penta-IJMPA,Chung82,Bagan93,Oka96,
WangHbaryon-1,WangHbaryon-2,WangHbaryon-3,WangHbaryon-4,WangHbaryon-5,WangHbaryon-6,Wang-cc-baryon-penta}, where the $U^\pm(p,s)$ are the Dirac spinors.

 At the hadron  side, we insert  a complete set  of intermediate pentaquark states with the same quantum numbers as the current operators  $J(x)$  and $i\gamma_5 J(x)$ into the correlation function
$\Pi(p)$ to obtain the hadronic representation
\cite{SVZ79-1,SVZ79-2,PRT85}, and  isolate the  lowest negative-parity and positive-parity hidden-charm pentaquark states, and obtain the
 results:
\begin{eqnarray}\label{CF-Hadron}
\Pi(p) & = & \lambda_{-}^2  {\!\not\!{p}+ M_{-} \over M_{-}^{2}-p^{2}  }+  \lambda_{+}^2  {\!\not\!{p}- M_{+} \over M_{+}^{2}-p^{2}  } +\cdots  \, ,\nonumber\\
&=&\Pi^1(p^2)\!\not\!{p}+\Pi^0(p^2)\, .
 \end{eqnarray}

 We obtain the spectral densities at the hadron  side through  dispersion relation,
\begin{eqnarray}
\frac{{\rm Im}\Pi^1(s)}{\pi}&=& \lambda_{-}^2 \delta\left(s-M_{-}^2\right)+\lambda_{+}^2 \delta\left(s-M_{+}^2\right) =\, \rho^1_{H}(s) \, , \\
\frac{{\rm Im}\Pi^0(s)}{\pi}&=&M_{-}\lambda_{-}^2 \delta\left(s-M_{-}^2\right)-M_{+}\lambda_{+}^2 \delta\left(s-M_{+}^2\right)
=\rho^0_{H}(s) \, ,
\end{eqnarray}
where we add the subscript $H$ to denote  the hadron side,
then we introduce the  weight functions $\sqrt{s}\exp\left(-\frac{s}{T^2}\right)$ and $\exp\left(-\frac{s}{T^2}\right)$ to obtain the QCD sum rules
at the hadron side,
\begin{eqnarray}
\int_{4m_c^2}^{s_0}ds \left[\sqrt{s}\,\rho^1_{H}(s)+\rho^0_{H}(s)\right]\exp\left( -\frac{s}{T^2}\right)
&=&2M_{-}\lambda_{-}^2\exp\left( -\frac{M_{-}^2}{T^2}\right) \, ,
\end{eqnarray}
where the $s_0$ is the continuum threshold parameter, and the $T^2$ is the Borel parameter.
We take the combinations  $\sqrt{s}\,\rho^1_{H}(s)\pm\rho^0_{H}(s)$ to separate the  contributions  of the pentaquark  states with the negative parity and positive parity unambiguously.

At the QCD side,  we contract the $u$, $d$, $s$ and $c$ quark fields in the correlation function $\Pi(p)$  with Wick theorem,
\begin{eqnarray}\label{QCD-Pi}
\Pi(p)&=&-i\,\varepsilon^{ila}\varepsilon^{ijk}\varepsilon^{lmn}\varepsilon^{i^{\prime}l^{\prime}a^{\prime}}\varepsilon^{i^{\prime}j^{\prime}k^{\prime}}
\varepsilon^{l^{\prime}m^{\prime}n^{\prime}}\int d^4x e^{ip\cdot x} \nonumber\\
&&   {\rm Tr}\Big[\gamma_5 D_{kk^\prime}(x) \gamma_5 C U^{T}_{jj^\prime}(x)C\Big] \,{\rm Tr}\Big[\gamma_5 C_{nn^\prime}(x) \gamma_5 C S^{T}_{mm^\prime}(x)C\Big] C C_{a^{\prime}a}^T(-x)C   \, ,
\end{eqnarray}
where
the $U_{ij}(x)$, $D_{ij}(x)$, $S_{ij}$ and $C_{ij}(x)$ are the full $u$, $d$, $s$ and $c$ quark propagators, respectively,
 \begin{eqnarray}
U/D_{ij}(x)&=& \frac{i\delta_{ij}\!\not\!{x}}{ 2\pi^2x^4}-\frac{\delta_{ij}\langle
\bar{q}q\rangle}{12} -\frac{\delta_{ij}x^2\langle \bar{q}g_s\sigma Gq\rangle}{192} -\frac{ig_sG^{a}_{\alpha\beta}t^a_{ij}(\!\not\!{x}
\sigma^{\alpha\beta}+\sigma^{\alpha\beta} \!\not\!{x})}{32\pi^2x^2} -\frac{\delta_{ij}x^4\langle \bar{q}q \rangle\langle g_s^2 GG\rangle}{27648} \nonumber\\
&&  -\frac{1}{8}\langle\bar{q}_j\sigma^{\mu\nu}q_i \rangle \sigma_{\mu\nu}+\cdots \, ,
\end{eqnarray}
\begin{eqnarray}
S_{ij}(x)&=& \frac{i\delta_{ij}\!\not\!{x}}{ 2\pi^2x^4}
-\frac{\delta_{ij}m_s}{4\pi^2x^2}-\frac{\delta_{ij}\langle
\bar{s}s\rangle}{12} +\frac{i\delta_{ij}\!\not\!{x}m_s
\langle\bar{s}s\rangle}{48}-\frac{\delta_{ij}x^2\langle \bar{s}g_s\sigma Gs\rangle}{192}+\frac{i\delta_{ij}x^2\!\not\!{x} m_s\langle \bar{s}g_s\sigma
 Gs\rangle }{1152}\nonumber\\
&& -\frac{ig_s G^{a}_{\alpha\beta}t^a_{ij}(\!\not\!{x}
\sigma^{\alpha\beta}+\sigma^{\alpha\beta} \!\not\!{x})}{32\pi^2x^2} -\frac{\delta_{ij}x^4\langle \bar{s}s \rangle\langle g_s^2 GG\rangle}{27648}-\frac{1}{8}\langle\bar{s}_j\sigma^{\mu\nu}s_i \rangle \sigma_{\mu\nu}  +\cdots \, ,
\end{eqnarray}
\begin{eqnarray}
C_{ij}(x)&=&\frac{i}{(2\pi)^4}\int d^4k e^{-ik \cdot x} \left\{
\frac{\delta_{ij}}{\!\not\!{k}-m_c}
-\frac{g_sG^n_{\alpha\beta}t^n_{ij}}{4}\frac{\sigma^{\alpha\beta}(\!\not\!{k}+m_c)+(\!\not\!{k}+m_c)
\sigma^{\alpha\beta}}{(k^2-m_c^2)^2}\right.\nonumber\\
&&\left. -\frac{g_s^2 (t^at^b)_{ij} G^a_{\alpha\beta}G^b_{\mu\nu}(f^{\alpha\beta\mu\nu}+f^{\alpha\mu\beta\nu}+f^{\alpha\mu\nu\beta}) }{4(k^2-m_c^2)^5}+\cdots\right\} \, ,\nonumber\\
f^{\alpha\beta\mu\nu}&=&(\!\not\!{k}+m_c)\gamma^\alpha(\!\not\!{k}+m_c)\gamma^\beta(\!\not\!{k}+m_c)\gamma^\mu(\!\not\!{k}+m_c)\gamma^\nu(\!\not\!{k}+m_c)\, ,
\end{eqnarray}
and  $t^n=\frac{\lambda^n}{2}$, the $\lambda^n$ is the Gell-Mann matrix
\cite{PRT85,Pascual-1984,WangHuang3900}.
We retain the terms $\langle\bar{q}_j\sigma_{\mu\nu}q_i \rangle$ and $\langle\bar{s}_j\sigma_{\mu\nu}s_i \rangle$ come from the Fierz re-arrangements   of the
$\langle q_i \bar{q}_j\rangle$ and $\langle s_i \bar{s}_j\rangle$  to  absorb the gluons  emitted from other quark lines to  extract the mixed condensate   $\langle\bar{q}g_s\sigma G q\rangle$ and $\langle\bar{s}g_s\sigma G s\rangle$, respectively \cite{WangHuang3900}.  Then we compute  all the integrals to obtain the correlation function $\Pi(p)$   at the quark-gluon  level, and finally obtain the QCD spectral densities through   dispersion relation,
\begin{eqnarray}\label{QCD-rho}
 \rho^1_{QCD}(s) &=&\frac{{\rm Im}\Pi^1(s)}{\pi}\, , \nonumber\\
\rho^0_{QCD}(s) &=&\frac{{\rm Im}\Pi^0(s)}{\pi}\, .
\end{eqnarray}
      There are two $c$-quark propagators and three light quark propagators  in Eq.\eqref{QCD-Pi}, if each $c$-quark line emits a gluon and each light quark line contributes  a quark-antiquark  pair, we obtain an operator $G_{\mu\nu}G_{\alpha\beta}\bar{u}u\bar{d}d\bar{s}s$ of dimension 13, so we should take into account the vacuum condensates at least up to dimension $13$ in a consistent way. Furthermore, we take into account the terms proportional to the $m_s$ to account for the flavor   $SU(3)$ mass-breaking effects.    The vacuum condensates
  $\langle\bar{q} q\rangle^2\langle\bar{q}g_s\sigma Gq\rangle $, $\langle\bar{q} q\rangle\langle\bar{s} s\rangle\langle\bar{q}g_s\sigma Gq\rangle $,  $\langle\bar{q} q\rangle^2\langle\bar{s}g_s\sigma Gs\rangle $, $\langle\bar{q} q\rangle \langle\bar{q}g_s\sigma Gq\rangle^2 $, $\langle\bar{s} s\rangle \langle\bar{q}g_s\sigma Gq\rangle^2 $, $\langle\bar{q} q\rangle \langle\bar{q}g_s\sigma Gq\rangle\langle\bar{s}g_s\sigma Gs\rangle  $,
  $\langle \bar{q}q\rangle^3\langle \frac{\alpha_s}{\pi}GG\rangle$, $\langle \bar{q}q\rangle^2\langle \bar{s}s\rangle\langle \frac{\alpha_s}{\pi}GG\rangle$ neglected in Refs. \cite{WangZG-EPJC-1509-12,WangZG-NPB-1512-32,WangZhang-APPB} are of dimension $11$ and $13$ respectively, and  are  associated with the inverse Borel parameters $\frac{1}{T^2}$, $\frac{1}{T^4}$, $\frac{1}{T^6}$ and $\frac{1}{T^8}$, which manifest themselves at  the small values of the Borel parameter $T^2$ and play an important role in determining the Borel windows.   We take the truncations $n\leq 13$ and $k\leq 1$ in a consistent way, just like in the new analysis in Ref.\cite{WZG-penta-IJMPA},
the quark-gluon operators of the orders $\mathcal{O}( \alpha_s^{k})$ with $k> 1$ and dimension $n>13$ are  discarded.

Now we  match the hadron side with the QCD side of the correlation function $\Pi(p)$, take the quark-hadron duality below the continuum threshold  $s_0$, and  obtain  the  QCD sum rules:
\begin{eqnarray}\label{QCDSR}
2M_{-}\lambda_{-}^2\exp\left( -\frac{M_{-}^2}{T^2}\right)&=& \int_{4m_c^2}^{s_0}ds \,\rho_{QCD}(s)\,\exp\left( -\frac{s}{T^2}\right)\,  ,
\end{eqnarray}
where $\rho_{QCD}(s)=\sqrt{s}\rho_{QCD}^1(s)+\rho_{QCD}^{0}(s)$, the explicit expressions of the QCD spectral densities are neglected for simplicity.

We differentiate     Eq.\eqref{QCDSR} with respect to  $\frac{1}{T^2}$, then eliminate the
 pole residue $\lambda_{-}$ and obtain the QCD sum rules for
 the mass of the hidden-charm  pentaquark state with the strangeness $S=-1$,
 \begin{eqnarray}
 M^2_{-} &=& \frac{-\int_{4m_c^2}^{s_0}ds \frac{d}{d(1/T^2)}\, \rho_{QCD}(s)\,\exp\left( -\frac{s}{T^2}\right)}{\int_{4m_c^2}^{s_0}ds \, \rho_{QCD}(s)\,\exp\left( -\frac{s}{T^2}\right)}\,  .
\end{eqnarray}

\section{Numerical results and discussions}
We take  the standard values of the  vacuum condensates
$\langle\bar{q}q \rangle=-(0.24\pm 0.01\, \rm{GeV})^3$,  $\langle\bar{s}s \rangle=(0.8\pm0.1)\langle\bar{q}q \rangle$,
 $\langle\bar{q}g_s\sigma G q \rangle=m_0^2\langle \bar{q}q \rangle$, $\langle\bar{s}g_s\sigma G s \rangle=m_0^2\langle \bar{s}s \rangle$,
$m_0^2=(0.8 \pm 0.1)\,\rm{GeV}^2$, $\langle \frac{\alpha_s
GG}{\pi}\rangle=0.012\pm0.004\,\rm{GeV}^4$    at the energy scale  $\mu=1\, \rm{GeV}$
\cite{SVZ79-1,SVZ79-2,PRT85,ColangeloReview}, and  take the $\overline{MS}$ masses $m_{c}(m_c)=(1.275\pm0.025)\,\rm{GeV}$
 and $m_s(\mu=2\,\rm{GeV})=(0.095\pm0.005)\,\rm{GeV}$
 from the Particle Data Group \cite{PDG}.
Furthermore,  we take into account
the energy-scale dependence of  the quark condensates, mixed quark condensates and $\overline{MS}$ masses according to  the renormalization group equation \cite{Narison-mix},
 \begin{eqnarray}
 \langle\bar{q}q \rangle(\mu)&=&\langle\bar{q}q\rangle({\rm 1 GeV})\left[\frac{\alpha_{s}({\rm 1 GeV})}{\alpha_{s}(\mu)}\right]^{\frac{12}{33-2n_f}}\, , \nonumber\\
 \langle\bar{s}s \rangle(\mu)&=&\langle\bar{s}s \rangle({\rm 1 GeV})\left[\frac{\alpha_{s}({\rm 1 GeV})}{\alpha_{s}(\mu)}\right]^{\frac{12}{33-2n_f}}\, , \nonumber\\
 \langle\bar{q}g_s \sigma Gq \rangle(\mu)&=&\langle\bar{q}g_s \sigma Gq \rangle({\rm 1 GeV})\left[\frac{\alpha_{s}({\rm 1 GeV})}{\alpha_{s}(\mu)}\right]^{\frac{2}{33-2n_f}}\, ,\nonumber\\
  \langle\bar{s}g_s \sigma Gs \rangle(\mu)&=&\langle\bar{s}g_s \sigma Gs \rangle({\rm 1 GeV})\left[\frac{\alpha_{s}({\rm 1 GeV})}{\alpha_{s}(\mu)}\right]^{\frac{2}{33-2n_f}}\, ,\nonumber\\
m_c(\mu)&=&m_c(m_c)\left[\frac{\alpha_{s}(\mu)}{\alpha_{s}(m_c)}\right]^{\frac{12}{33-2n_f}} \, ,\nonumber\\
m_s(\mu)&=&m_s({\rm 2GeV} )\left[\frac{\alpha_{s}(\mu)}{\alpha_{s}({\rm 2GeV})}\right]^{\frac{12}{33-2n_f}}\, ,\nonumber\\
\alpha_s(\mu)&=&\frac{1}{b_0t}\left[1-\frac{b_1}{b_0^2}\frac{\log t}{t} +\frac{b_1^2(\log^2{t}-\log{t}-1)+b_0b_2}{b_0^4t^2}\right]\, ,
\end{eqnarray}
  where $t=\log \frac{\mu^2}{\Lambda^2}$, $b_0=\frac{33-2n_f}{12\pi}$, $b_1=\frac{153-19n_f}{24\pi^2}$, $b_2=\frac{2857-\frac{5033}{9}n_f+\frac{325}{27}n_f^2}{128\pi^3}$,  $\Lambda=213\,\rm{MeV}$, $296\,\rm{MeV}$  and  $339\,\rm{MeV}$ for the flavors  $n_f=5$, $4$ and $3$, respectively  \cite{PDG}.
In the present work, we study  the hidden-charm pentaquark state $udsc\bar{c}$,  it is better to choose the flavor numbers $n_f=4$, and evolve those input parameters to a typical energy scale $\mu$, which satisfy the energy scale formula  $\mu=\sqrt{M_{P}-(2{\mathbb{M}}_c)^2}$ with the updated value ${\mathbb{M}}_c=1.82\,\rm{GeV}$ \cite{Wang1508-EPJC,WangHuang-EPJC-1508-12,WangZG-EPJC-1509-12,WangZG-NPB-1512-32,WZG-tetraquark-Mc,WangHuang3900,
Wang-tetra-formula-1,Wang-tetra-formula-2,Wang-tetra-IJMPA-1,Wang-tetra-IJMPA-2,Wang-tetra-IJMPA-3}. We tentatively  assign the $P_{cs}(4459)$ to be  the $[ud][sc]\bar{c}$ ($0$, $0$, $0$, $\frac{1}{2}$) pentaquark state and take into account the flavor $SU(3)$ mass-breaking effects, the optimal energy scale of the QCD spectral density is  $\mu=\sqrt{M_{P}-(2{\mathbb{M}}_c)^2}-m_s=2.4\,\rm{GeV}$.

In the QCD sum rules for  the  heavy, doubly-heavy, triply-heavy baryon states, hidden-charm pentaquark states and doubly-charm pentaquark states,  we usually choose the continuum threshold parameters as $\sqrt{s_0}=M_{gr}+ (0.5-0.8)\,\rm{GeV}$  \cite{Wang1508-EPJC,WangHuang-EPJC-1508-12,
 WangZG-EPJC-1509-12,WangZG-NPB-1512-32,WangZhang-APPB,WZG-penta-IJMPA,WangHbaryon-1,WangHbaryon-2,WangHbaryon-3,WangHbaryon-4,WangHbaryon-5,WangHbaryon-6,Wang-cc-baryon-penta},   where we use the subscript $gr$ to denote  the ground states.
 In the present work,  we choose the continuum threshold parameters to be  $\sqrt{s_0}= 5.15\pm0.10\,\rm{GeV}$, and verify the possible assignment of the $P_{cs}(4459)$ as the $[ud][sc]\bar{c}$ ($0$, $0$, $0$, $\frac{1}{2}$) pentaquark state.

 We obtain the  Borel  window $T^2=3.4-3.8\,\rm{GeV}^2$ via trial  and error, the pole contribution is about $(40-61)\%$, which is large enough to extract the pentaquark mass reliably, just as in our previous works on the hidden-charm pentaquark states and doubly-charm pentaquark states  \cite{Wang1508-EPJC,WangHuang-EPJC-1508-12,
 WangZG-EPJC-1509-12,WangZG-NPB-1512-32,WangZhang-APPB,WZG-penta-IJMPA,Wang-cc-baryon-penta}. In Fig.\ref{Pole-Borel}, we plot the pole contribution  with variation of the  Borel parameter $T^2$ at much larger range  than the Borel window.
  \begin{figure}
\centering
\includegraphics[totalheight=7cm,width=10cm]{Pole.EPS}
  \caption{ The pole contributions  with variation of the  Borel parameter $T^2$, where the $A$, $B$ and $C$ correspond to the continuum threshold parameter $\sqrt{s_0}=5.05\,\rm{GeV}$, $5.15\,\rm{GeV}$ and $5.25\,\rm{GeV}$, respectively,   the region between the two vertical lines is the Borel window. }\label{Pole-Borel}
\end{figure}

In Fig.\ref{OPE-Borel}, we plot the contributions of the higher dimensional  vacuum condensates with variation of the  Borel parameter $T^2$. From the figure, we can see that the higher dimensional  vacuum condensates  manifest themselves at the region $T^2< 2.5\,\rm{GeV}^2$, we should choose the value $T^2> 2.5\,\rm{GeV}^2$. Their values decrease monotonously and quickly with the increase of the Borel parameter, in the Borel window $T^2=3.4-3.8\,\rm{GeV}^2$, the contributions of the higher dimensional vacuum condensates are  $D(8)=-(21-26)\%$, $D(9)=(6-8)\%$, $D(10)=(2-3)\%$, $D(11)=-(1-2)\%$, $D(13)\ll1\%$, the convergent behavior is very good.

The higher dimensional vacuum condensates play an important role in obtaining the flat platform, especially those associated with the inverse Borel parameters  $\frac{1}{T^2}$, $\frac{1}{T^4}$, $\frac{1}{T^6}$ and $\frac{1}{T^8}$. In Fig.\ref{mass-D10-fig}, we plot the predicted pentaquark mass  with variation of the  Borel parameter $T^2$ with the  truncations of the operator product expansion up to the vacuum condensates of dimensions $n=10$ and $13$, respectively.  From the figure, we can see that the vacuum condensates of dimensions $11$ and $13$ play an important role to obtain the flat platform, we should take them into account in a consistent way.

\begin{figure}
\centering
\includegraphics[totalheight=7cm,width=10cm]{OPE.EPS}
  \caption{ The  contributions  of the higher dimensional  vacuum condensates $D(n)$ with variation of the  Borel parameter $T^2$, where the $n=8$, $9$, $10$, $11$ and $13$,   the region between the two vertical lines is the Borel window. }\label{OPE-Borel}
\end{figure}

\begin{figure}
\centering
\includegraphics[totalheight=7cm,width=10cm]{mass-D10.EPS}
  \caption{ The predicted  mass  with variation of the  Borel parameter $T^2$, where the $n=10$ and $13$ denote truncations of the operator product expansion, the expt denotes the experimental value. }\label{mass-D10-fig}
\end{figure}

Now we take into account  all uncertainties  of the input   parameters,
and obtain  the mass and pole residue of
 the   hidden-charm pentaquark state $udsc\bar{c}$ with the negative-parity, which are shown explicitly in Figs.\ref{mass}-\ref{pole-residue},
\begin{eqnarray}
M_P&=&4.47\pm0.11\,\rm{GeV}\, ,\nonumber\\
\lambda_P&=&\left(1.86\pm0.31\right)\times 10^{-3}\,\rm{GeV^6}\, .
\end{eqnarray}
The predicted mass  $M_P=4.47\pm0.11\,\rm{GeV}$ is in excellent agreement with the experimental data $ 4458.8 \pm 2.9 {}^{+4.7}_{-1.1} \mbox{ MeV}$ from the  LHCb collaboration  \cite{LHCb-Pcs4459,LHCb-Pcs4459-2012}, and supports assigning the $P_{cs}(4459)$   as the $[ud][sc]\bar{c}$ ($0$, $0$, $0$, $\frac{1}{2}$) pentaquark state
with the spin-parity $J^P={\frac{1}{2}}^-$. In Ref.\cite{WZG-penta-IJMPA}, we observe that the $P_c(4312)$ can be assigned to be the $[ud][uc]\bar{c}$ ($0$, $0$, $0$, $\frac{1}{2}$) pentaquark state with the spin-parity $J^P={\frac{1}{2}}^-$. The flavor  $SU(3)$ mass-breaking effect is about $\Delta=m_s=147\,\rm{MeV}$, the estimations presented in Table \ref{mass-1508-et al-Pcs} are reasonable and reliable. In Refs.\cite{WangZG-EPJC-1509-12,WangZG-NPB-1512-32}, we study the $J^P={\frac{1}{2}}^\pm$ and ${\frac{3}{2}}^\pm$ hidden-charm pentaquark states with the strangeness $S=0$, $-1$, $-2$ and $-3$ in a systematic way by carrying out the operator product expansion up to the vacuum condensates of dimension $10$ and choosing the old value ${\mathbb{M}}_c=1.80\,\rm{GeV}$, and observe that the flavor  $SU(3)$ mass-breaking effects  are $\Delta=(90-130)\,\rm{MeV}$ for the negative-parity pentaquark states. The new analysis supports a larger flavor  $SU(3)$ mass-breaking effect.

The constituent quark models have given many successful descriptions of the hadron spectroscopy, the simple constituent quark mass plus hyperfine spin-spin interaction model works well for the ground state mesons ($M$) and baryons ($B$),
\begin{eqnarray}
M_M=\sum\limits_{i=1,2}\, m_i+\frac{C_{12}}{m_1 m_2}\,\vec{s}_1 \cdot \vec{s}_2\, ,\nonumber\\
M_B=\sum\limits_{i=1,2,3}\,m_i +\sum\limits_{i<j}\,\frac{C_{ij}}{m_i m_j}\,\vec{s}_i \cdot \vec{s}_j\, ,
\end{eqnarray}
where the $m_i$  are the constituent quark masses, the $\vec{s}_i$ are the quark spins,  the coefficients $C_{ij}$ can be fitted phenomenologically \cite{mass-spin-M}.  Accordingly, we can study the ground state hidden-charm tetraquark ($T$) and pentaquark ($P$) states  with the simple constituent diquark mass plus hyperfine spin-spin interaction model,
\begin{eqnarray}
M_T= m_{\mathcal{ D}}+m_{\mathcal{ \bar{D}}}+\sum\limits_{i<j}\,\widetilde{C}_{ij}\,\vec{s}_i \cdot \vec{s}_j\, ,\nonumber\\
M_P= m_{\mathcal{ D}}+m_{\mathcal{ D}^\prime}+m_{\bar{c}}+\sum\limits_{i<j}\,\widetilde{C}_{ij}\,\vec{s}_i \cdot \vec{s}_j\, ,
\end{eqnarray}
where the coefficients $\widetilde{C}_{ij}=\frac{C_{ij}}{m_im_j}$, the $m_D$  are the constituent diquark masses,  and the $\vec{s}_i$ are the quark spins \cite{Karliner-B-cs,Ali-EPJC-2018,Karliner-S-junction}.
 If we want to take into account the light flavor $SU(3)$-breaking effects satisfactorily, we have to introduce additional parameters, the bounding terms or string-junctions beyond the quark masses, diquark masses and  coefficients $\widetilde{C}_{ij}$  \cite{Karliner-B-cs,Karliner-S-junction}. The parameters $\widetilde{C}_{ij}$ in the mesons,  baryons, tetraquark states and pentaquark states are not necessary to be the same, as they correspond to different quark structures and wave functions, the $SU(3)$ breaking effects in the mesons,  baryons, tetraquark states and pentaquark states can be quite different. Generally speaking, they can be fitted phenomenologically when the experimental data accumulate. For the traditional charmed baryon states, the  flavor $SU(3)$ mass-breaking effects due to the $s$-quark are $\Delta\approx 190\,\rm{MeV}$ and $120\,\rm{MeV}$ for the flavor antitriplets and sextets respectively from the Particle Data Group \cite{PDG}, which can be accounted  for satisfactorily by adding the bonding terms \cite{Karliner-B-cs}. For the hidden-charm pentaquark states, there are only five candidates $P_c(4312)$, $P_c(4380)$, $P_c(4440)$, $P_c(4457)$ and $P_{cs}(4459)$, their quantum numbers such as the spin, parity, $\cdots$, have not been determined yet, it is difficult to estimate the $SU(3)$ breaking effects in a simple and intuitive way. In the present work, we study the $SU(3)$ breaking effects with the QCD sum rules, and account for the $SU(3)$ breaking effects by the $s$-quark mass $m_s$, $s$-quark condensate $\langle\bar{s}s\rangle$ and $s$-quark mixed condensate $\langle\bar{s}g_s\sigma Gs\rangle$ , and obtain the mass gap about $\Delta=147\,\rm{MeV}$, and make crude estimations of the mass-spectrum of the hidden-charm pentaquark states with the strangeness $S=-1$,   direct but tedious calculations are still needed to make robust predictions.

\begin{figure}
\centering
\includegraphics[totalheight=7cm,width=10cm]{mass.EPS}
  \caption{ The mass  with variation of the  Borel parameter $T^2$. }\label{mass}
\end{figure}

\begin{figure}
\centering
\includegraphics[totalheight=7cm,width=10cm]{residue.EPS}
  \caption{ The pole residue with variation of the  Borel parameter $T^2$. }\label{pole-residue}
\end{figure}

\section{Conclusion}
  In this article, we study the scalar-diquark-scalar-diquark-antiquark type $udsc\bar{c}$ pentaquark state with the QCD sum rules in details. We take into account  all the    vacuum condensates up to dimension $13$ in a consistent way with the truncations $\mathcal{O}(m_s^k)$ for $k\leq1$, and use the energy scale formula $\mu=\sqrt{M_{P}-(2{\mathbb{M}}_c)^2}-m_s$ with the updated effective $c$-quark mass ${\mathbb{M}}_c=1.82\,\rm{GeV}$ and subtract the mass-breaking effect due to the small $s$-quark mass  to choose  the optimal  energy scales of the QCD spectral densities. The predicted mass $M_P=4.47\pm0.11\,\rm{GeV}$ is in excellent agreement with the experimental data $ 4458.8 \pm 2.9 {}^{+4.7}_{-1.1} \mbox{ MeV}$ from the  LHCb collaboration, and supports assigning the $P_{cs}(4459)$   to be the $[ud][sc]\bar{c}$ ($0$, $0$, $0$, $\frac{1}{2}$) pentaquark state
with the spin-parity $J^P={\frac{1}{2}}^-$. We take into account the flavor   $SU(3)$ mass-breaking effect about $150\,\rm{MeV}$, and estimate the mass spectrum of the diquark-diquark-antiquark type $udsc\bar{c}$ pentaquark states with the strangeness $S=-1$.

\section*{Acknowledgements}
This  work is supported by National Natural Science Foundation, Grant Number  11775079.