\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. 



