text:
# Gravitational collapse of Type II fluid in higher dimensional space-times

## I Introduction

The gravitational collapse of sufficiently massive star, under fairly general conditions, will result into a singularity is the fact established by the singularity theorems of Hawking and Ellis he . However these theorems do not indicate whether the singularity will be hidden behind the event horizon or whether it will be visible to an outside observer. The singularity will not be visible is indicated in the relativistic literature by the phrase cosmic censorship conjecture (CCC). It will not be visible to an observer outside the event horizon is the weak form while it will not be visible to any observer, even the one who is sitting on the collapsing star, is the strong form of the conjecture. This conjecture was articulated by Penrose rp  (see r1  for reviews). The weak form essentially states that gravitational collapse from a regular initial data never creates the space-time singularity visible to distant observers, i.e., any singularity that forms must be hidden within a black hole. For the strong form on the other hand, nothing is supposed to emanate out of the singularity and hence it is not visible to anyone even if one is infinitesimally close to it.

There is no general theory of nature and visibility of singularities. There do exist a number of exact solutions of the Einstein equation which admit depending upon the initial data black holes (BH) or naked singularities (NS) (see, e.g., Joshi r1 ). In particular the Vaidya solution pc  is extensively used to show that the end state of collapse for a regular initial data results into a naked singularity. In the absence of a general result, the study of various examples of collapse becomes pertinent to examine the validity of the phenomenon of occurrence of naked singularities. The two most investigated cases are that of inhomogeneous dust and the Vaidya null fluid. A more general matter distribution should be considered. Undoubtedly, as singularity is approached, matter would be in highly dense and exotic state. The strange quark matter is the densest state of matter known which is produced in quark-hadron phase transition in the early Universe or at ultra high energy neutron stars converting into strange quark matter stars ew . It would therefore be interesting to study the collapse of the mixture of the Vaidya fluid with the strange matter (SQM), which is a Type II fluid he . This has been studied in the usual 4-dimensional space-time gd1  and here we wish to study it in higher dimensions. The presence of SQM tends to shrink the initial data window leading to NS. As is expected, the increase in dimensions would also tend to enhance the shrinkage because of gravity getting stronger in higher dimensions gd .

Recent developments in string theory indicate that gravity may be truly higher-dimensional(HD) interaction, becoming effectively $`4D`$ at lower energies. The main question is how to bring this theoretical development on the anvil of observations? One of the ways is to study its effects in astrophysical settings. That is where the phenomenon of gravitational collapse in HD attains relevance. In this study we would be finding the effects of HD as well as SQM on the end product of collapse of null SQM. For this, we first obtain the exact solution describing collapse of null SQM fluid sphere, and then investigate their effect in the context of CCC.

In section II, we obtain the generalized Vaidya solution for null quark fluid in higher dimensional spherically symmetric space-time. For definiteness we shall refer this metric as HD generalized Vaidya metric. We adapt this solution to study end state of collapse in section III . Finally we shall conclude with a discussion on the effect of SQM on HD Vaidya collapse.

## II Type II fluid collapse in higher dimensions

Let us begin with the general $`(n+2)`$-dimensional spherically symmetric space-time, in advanced Eddington time coordinate $`v`$, described by the metric pd ; jr :

$$ds^2=e^{\psi (v,r)}dv\left[f(v,r)e^{\psi (v,r)}dv+2dr\right]+r^2d\mathrm{\Omega }_n^2$$
(1)
where $`0r\mathrm{}`$ is the proper radial coordinate, $`\mathrm{}v\mathrm{}`$ is an avanced time coordiante, and

$`d\mathrm{\Omega }^2={\displaystyle \underset{i=1}{\overset{n}{}}}\left[{\displaystyle \underset{j=1}{\overset{i1}{}}}\mathrm{sin}^2\theta _j\right]d\theta _i^2=d\theta _1^2+sin^2\theta _1d\theta _2^2+`$
$`sin^2\theta _1sin^2\theta _2d\theta _3^2+...`$
$`+sin^2\theta _1sin^2\theta _2...sin^2\theta _{n1}d\theta _n^2`$ (2)
is the metric on $`n`$-sphere in spherical polar coordinates. It proves useful to introduce a local mass function $`m(v,r)`$ defined by $`f=12m/r`$ bi . Here $`m(v,r)`$ is an arbitrary function of advanced time $`v`$ and radial coordinate $`r`$. When $`m=m(v)`$, it is the Vaidya solution in higher dimensions pd ; gd . The usual Vaidya solution in 4-space-time follows for $`m=m(v)`$ and $`n=2`$.

We wish to find the general solution of the Einstein equation for the matter field given by Eq.(3) for the metric (1), which contains two arbitrary functions. Hence it is general in the retarded coordinates used. It is the field equation $`G_1^0=0`$ that leads as in the $`4D`$ case as well as in the present HD case to $`e^{\psi (v,r)}=g(v)`$ gd1 . This could be absorbed by writing $`d\stackrel{~}{v}=g(v)dv`$. Hence, without loss of generality, the metric (2) takes the form,

$$ds^2=\left[1\frac{2m(v,r)}{(n1)r^{n1}}\right]dv^2+2dvdr+r^2d\mathrm{\Omega }_n^2$$
(3)
The energy-momentum tensor for the Type II (null SQM) fluid is given by ww ; vh ,

$$T_{ab}=\mu l_al_b+(\rho +p)(l_a\eta _b+l_b\eta _a)+pg_{ab}$$
(4)
$`l_a=\delta _a^0,n_a={\displaystyle \frac{1}{2}}\left[1{\displaystyle \frac{2m(v,r)}{r}}\right]\delta _a^0\delta _a^1`$
$`l^a=\delta _1^a,n^a=\delta _0^a{\displaystyle \frac{1}{2}}\left[1{\displaystyle \frac{2m(v,r)}{r}}\right]\delta _1^a`$
$`l_al^a=n_an^a=0l_an^a=1.`$ (5)
Here $`\rho ,p`$ are the strange quark matter energy density and thermodynamic pressure while $`\mu `$ is the energy density of the Vaidya null radiation. The null vector $`l_a`$ is a double null eigenvector of $`T_{ab}`$. Physically occurring distribution is null radiation flowing in the radial direction corresponding to $`\rho =p=0`$, the Vaidya space-time of radiating star. When $`\mu =0`$, $`T_{ab}`$ reduces to degenerate Type I fluid pd , and it represents string dust for $`\mu =0=p`$. The energy condition for such a distribution are as follows he :
(i) the weak and strong energy conditions,

$$\mu >0,\rho 0,p0$$
(6)
(ii) the dominant energy condition,

$$\mu >0\rho p0.$$
(7)
By choosing properly the mass function $`m(v,r)`$, it is possible to satisfy the energy conditions. In the particular case of $`m=m(v)`$, the energy condition reduces to $`\mu 0`$.

The field equations are pd ,

$`8\pi \mu ={\displaystyle \frac{n\dot{m}}{(n1)r^n}},`$ (8a)
$`8\pi \rho ={\displaystyle \frac{nm^{}}{(n1)r^n}},`$ (8b)
$`8\pi p={\displaystyle \frac{m^{\prime \prime }}{(n1)r^{n1}}}.`$ (8c)
Here and in what follows dash and dot denote derivative respectively $`/r`$ and $`/v`$. The part $`\mu l_al_b`$ of $`T_{ab}`$ in (4) is the component of matter field that moves along the null hypersurface $`v=const.`$. In particular when $`p=\rho =0`$ we have the Vaidya solution in higher dimensions. Thus the distribution in (4) represents Vaidya radiating star in Type II fluid universe in higher dimensions. Note that when $`\mu =0`$, it is static perfect fluid space-time, which will yield regular space-time only when $`mr^{n+1}`$ which is the de Sitter space. As in the $`4D`$ case, the equation of state for the zero temperature quark matter is taken hc :

$$p=\frac{1}{k}(\rho 4B)$$
(9)
where $`k>0`$ is a constant. Inserting the above equation in the field equations, we find that

$$m^{\prime \prime }=\frac{nm^{}}{kr}+\frac{32\pi B}{k}r^{n1}.$$
(10)
If we choose the function $`m(v,r)`$ such that

$$m(v,r)=m_0(v,r)+\frac{\mathrm{\Lambda }}{(n+1)}r^{n+1},$$
(11)
where $`m_0(v,r)`$ is an unknown function, then we find that

$$m_0^{\prime \prime }=\frac{n}{k}m_0^{},B=\frac{32\pi \mathrm{\Lambda }}{n(k+1)}.$$
(12)
It admits the general solution

$$m_0=S(v)r^{n/k+1}+M(v).$$
(13)
Now, one can easily calculate $`m`$, $`\mu `$, $`\rho `$ and $`p`$, which are

$`m=M(v)+S(v)r^{n/k+1}+{\displaystyle \frac{\mathrm{\Lambda }}{n+1}}r^{n+1}`$ (14)
$`\mu ={\displaystyle \frac{n}{8\pi (n1)r^n}}\left[\dot{M}(v)+\dot{S}(v)r^{n/k+1}\right]`$ (15)
$`\rho ={\displaystyle \frac{n}{8\pi (n1)r^n}}\left[{\displaystyle \frac{kn}{k}}S(v)r^{n/k}+\mathrm{\Lambda }r^n\right]`$ (16)
$`p={\displaystyle \frac{n}{8\pi (n1)r^n}}\left[{\displaystyle \frac{(kn)}{k^2}}S(v)r^{n/k+1}\mathrm{\Lambda }r^n\right].`$ (17)
This is the general solution of collapsing null strange quark fluid in higher dimensions. The metric explicitly reads as

$`ds^2=[1{\displaystyle \frac{2M(v)}{(n1)r^{(n1)}}}{\displaystyle \frac{2S(v)}{(n1)r^{[(n2)k+n]/k}}}`$
$`{\displaystyle \frac{\mathrm{\Lambda }}{n^21}}r^2]dv^2+2dvdr+r^2d\mathrm{\Omega }^2.`$ (18)
Clearly, all the energy conditions would be satisfied for $`n2`$ because it would ensure $`\rho 0`$ and $`p0`$, while $`\mu 0`$ would be taken care of when we choose the mass functions for both null radiation and SQM. The case, $`S(v)=0`$, corresponds to the HD Vaidya de Sitter solution pd ; gd , whereas the Vaidya solution is obtained by taking $`S(v)=B=0`$. The solutions in Refs. gd1 ; hc  can be recovered by setting $`n=2`$. In the next section, we shall adapt the solution obtained to study the end state of the collapse for $`k=3`$ in Eq. (9), for which the matter is quark plasma.

## III Existence of a naked singularity

The physical scenario here is that of a radial influx of null fluid in an initially empty region of the higher dimensional non-flat but empty space. The first shell arrives at $`r=0`$ at time $`v=0`$ and the final at $`v=T`$. A central singularity of growing mass is developed at $`r=0`$. For $`v<0`$ we have $`m(v)=\mathrm{\hspace{0.33em}0}`$, i.e., the higher dimensional de Sitter like space, and for $`v>T`$, $`\dot{m}(v)=\mathrm{\hspace{0.33em}0}`$, $`m(v)`$ is positive definite. The metric for $`v=0`$ to $`v=T`$ is the higher dimensional generalized Vaidya metric, and for $`v>T`$ we have the higher dimensional Schwarzschild solution. For the HD Vaidya region, we choose,

$`2M(v)=\alpha (n1)v^{n1}(\alpha >0),`$ (19a)
$`2S(v)=\beta (n1)v^{[(n2)k+n]/k}(\beta >0),`$ (19b)
Radial null geodesics of the metric (3) must satisfy

$$\frac{dr}{dv}=\frac{1}{2}\left[1\frac{2M(v)}{(n1)r^{(n1)}}\frac{2S(v)}{(n1)r^{[(n2)k+n]/k}}\frac{\mathrm{\Lambda }}{n^21}r^2\right].$$
(20)
Clearly the above differential equation has a singularity at $`r=0`$, $`v=0`$. We note that a singularity forms when a shell hits the origin at $`r=0`$, where the density diverges. If the singularity is at least locally naked (for brevity we have addressed it as naked in this paper), there must be a light ray coming out from it. The critical ray is the Cauchy horizon defined as the first outgoing radial null geodesic emerging from the singularity. Therefore, by investigating the behavior of radial null geodesics near the singularity, one can find out if outgoing null geodesics meet the singularity in the past. In order to determine the nature of the limiting value of $`y`$ at $`r=0`$, $`v=0`$ on a singular geodesic, we let $`y_0=lim_{r0v0}y=lim_{r0v0}v/r`$. Using Eq. (20) and L’Hôpital’s rule we get

$`y_0=\underset{r0v0}{lim}y=\underset{r0v0}{lim}{\displaystyle \frac{v}{r}}=\underset{r0v0}{lim}{\displaystyle \frac{dv}{dr}}=`$
$`\underset{r0v0}{lim}{\displaystyle \frac{2}{1\alpha y^{n1}\beta y^{n(k+1)/k}\frac{\mathrm{\Lambda }}{n^21}r^2}}`$ (21)
which can be rearranged as

$$\beta y_0^{[(n2)k+n]/k+1}+\alpha y_0^ny_0+2=0$$
(22)
This algebraic equation is the key equation which governs the behavior of the tangent vector near the singular point.

If the singularity is naked, Eq. (22) must have one or more positive roots $`y_0`$, i.e., at least one outgoing geodesic that will terminate in the past at the singularity. Hence in the absence of positive roots, the collapse will always lead to a black hole. Thus, the occurrence of positive roots implies that the strong CCC is violated, though not necessarily the weak one.

The limit $`\beta 0`$ corresponds to HD Vaidya collapse gd , and in that case, Eq. (22) admits positive roots for $`\alpha \alpha _C^v`$. Hence singularities are naked for $`\alpha (0,\alpha _C^v]`$, and black holes form otherwise. Thus $`\alpha =\alpha _C^v`$ is the critical value at which the transition occurs, the end state of collapse switches from naked singularities to black holes. This is in agreement with earlier work (see, e.g., gd ).

We want to investigate the changes, at least qualitatively, in the above picture of HD Vaidya collapse in the presence of a SQM. It is verified numerically that as $`n`$ increases, the threshold value $`\beta _T`$ decreases (see Table 2). For $`\beta >\beta _T`$, it is BH for all $`\alpha `$. Further, the critical parameter $`\alpha _C^s`$ for HD SQM is smaller than corresponding critical parameter $`\alpha _C^v`$ for HD Vaidya , i.e, the naked singularity spectrum of the collapsing HD Vaidya region gets covered and hence the initial data space $`(0,\alpha _C^v]`$ for a naked singularity in HD Vaidya collapse contracts due to the presence of SQM. Thus, at least at a qualitative level, the presence of SQM facilitates the formation of black holes. However, it is interesting to see that for each $`\beta `$, there exists a $`\alpha _C^s`$ such that singularities are always naked for all $`\alpha (0,\alpha _C^s]`$. That is for each $`\beta `$ there exists a non zero measure set of $`\alpha `$ values giving rise to NS and consequently violating CCC (see Tables 1, 2).

### III.1 Other Cases

####

Case $`k=n`$.
This corresponds to

$$\mu =\frac{n}{8\pi (n1)r^n}\left[\dot{M}(v)+\dot{S}(v)\right],\rho =p=\frac{n\mathrm{\Lambda }}{8\pi (n1)}$$
(23)
and the algebraic equation takes the form

$$(\alpha +\beta )y_0^ny_0+2=0$$
(24)
The metric in this case takes the form of the HD Vaidya - de Sitter metric, i.e., the HD collapse of null fluid in an expanding de Sitter back ground where $`\mathrm{\Lambda }`$ is generated by the bag constant. The singularity is naked for

$$(\alpha +\beta )\lambda _C=\frac{1}{n}\left(\frac{n1}{2n}\right)^{n1}$$

The equal roots at $`\lambda _C`$ is $`y_o=2n/(n1)`$. It is remarkable to note that both critical parameters and tangents to outgoing geodesics are dependent on dimension of space-time. Further, we note that $`y_0`$ is bounded below by value 2, $`y_02`$ as $`\lambda 0`$ or $`D\mathrm{}`$.

####

Case $`k\mathrm{}`$.
We have

$`\mu ={\displaystyle \frac{n}{16\pi (n1)r^2}}[\alpha y^{n2}+\beta y^{n3}],p=0,`$
$`\rho ={\displaystyle \frac{n\beta }{16\pi r^2}}y^{n2}`$ (25)
and the algebraic equation becomes:

$$\alpha y_0^n+\beta y_0^{n1}y_0+2=0,$$

For $`n=2`$, this would admit a positive root for $`\alpha 1/8(\beta 1)^2`$, giving the range for naked singularity as obtained in jdj . This is simply the null fluid collapse in the background of constant potential which is characterized by $`T_0^0=T_1^1=const./r^2`$, as is the case in Eq.(25) above. The condition for naked singularity modifies to $`\beta ^2<4\alpha +108\alpha ^2+36\alpha \beta +\beta ^2`$ when $`n=3`$.

####

Other $`k`$
We also note that as $`k`$ increases, so does the threshold value $`\beta _T`$. However, to conserve the space we shall not present the details.

## IV Discussion

Motivated by the development of superstring and other field theories, there is considerable interest in the models with extra dimensions from viewpoint of both cosmology rs  and gravitational collapse hd . In this paper, We have shown that the 4D spherically symmetric solution describing null SQM fluid go over to $`(n+2)`$-dimensional spherically symmetric solution and essentially retaining its physical behavior. The solution obtained includes many previously known solutions to Einstein equations, such as, HD Vaidya de Sitter solution. We have also used this solution to study the end state of collapsing star and showed that there exists a regular initial data which leads to naked singularity. Higher dimensions tend to shrink the naked singularity initial data space alternatively enlargement of black hole initial data space, of corresponding $`4D`$ collapse. We have also studied the effect of SQM on Vaidya collapse (both in $`4D`$ and $`HD`$). The naked singularity spectrum of Vaidya collapse shrinks further with introduction of SQM and thus atleast qualitatively SQM favors black hole in comparison to naked singularities. It is also interesting to see that there exists threshold value $`\beta _T`$, which decreases with increase in $`n`$, such that for $`\beta >\beta _T`$, the end state of the collapse is always black hole for all $`\alpha `$. However the threshold touches zero only as $`n\mathrm{}`$. That means for finite dimensions, there would always be a window howsoever narrow in the initial data leading to NS. However for sufficiently large $`n`$, it would practically be black hole and CCH would be obeyed. This case includes several previous cases of spherical collapse in $`4D`$ and $`HD`$.

###### Acknowledgements.

One of the authors (SGG) would like to thank IUCAA, Pune for hospitality while this work was done.


format: text

text:
# Generalized Competing Glauber-type Dynamics and Kawasaki-type Dynamics

## I Introduction

In recent years, there have been continuing efforts towards a clear picture of the self-organization phenomena in the phase transitions of magnetic systems. Most of the works have been concentrated on Ising and Ising-like systems, governed by two competing dynamics: Glauber’s single-spin flipping mechanism and Kawasaki’s spin-pair exchange mechanism, both with a probability. The system is coupled to a heat bath at a given temperature while being subject to an external energy flux. Glauber’s flipping mechanism is to simulate the contact of the system with the heat bath. Changing the order parameter, it favors lower system energy. On the other hand, Kawasaki’s exchange mechanism is to simulate the influence of the energy flux. Keeping the order parameter conserved, it favors higher system energy. With these two competing mechanisms and the corresponding master equation, people expect to obtain the evolution of the system. As exact treatment is not possible for $`2D`$ Ising model, consequently Monte Carlo simulation and methods as Dynamic Pair Approximation have been employed. The results obtained helped to determine the interesting phase diagrams. (However, people are surprised to find there contradictions in the predictions of MC simulations and the dynamical mean-field theory, since they are both proved successful to yield good qualitative results in other studies. Though revisions of MC simulations are made and more accurate versions of mean-field theory are presented, the puzzle still remains.)

In our earlier studies, we have presented single-spin transition mechanism and spin-pair redistribution mechanism. These two dynamics are natural generalizations of Glauber’s single-spin flipping mechanism and Kawasaki’s spin-pair exchange mechanism respectively. They have similar mathematical expressions, and become counterparts of each other in the non-conserved and conserved processes respectively. As an example of the applications, we studied kinetic Gaussian model with both of them respectively. Our study shows that, in translational-invariant lattices, the dynamic critical exponent $`z=1/v=2`$ is independent of space dimensionality and the governing dynamical mechanism. Its dynamic properties are summarized in Sec. IV.

In this article, we formulate the competing dynamics combining single-spin transition and spin-pair redistribution. As these two mechanisms themselves are universal, the combined one is also applicable to arbitrary systems, and it can be deemed as a generalization. In Sec. II, we first briefly review the two mechanisms respectively and then give the formulation of the competing mechanism. In Sec. III, we explain the differences between our mechanism and that adopted conventionally, taking $`1D`$ Ising model as an example. In Sec. IV, as the chief task of this article, we apply it to kinetic Gaussian model and report the findings. In Sec. V we summarize our study with some discussions.

## II The competing mechanism

First we briefly review the single-spin transition mechanism and spin-pair redistribution mechanism.

### A Single-spin transition mechanism

Glauber’s single-spin flipping mechanism allows an Ising system to evolve with its spins flipping to their opposite. In single-spin transition mechanism, a single spin $`\sigma _i`$ may change itself to any possible values, $`\widehat{\sigma }_i`$, and the master equation is

$$\frac{d}{dt}P(\{\sigma \},t)=\underset{i}{}\underset{\widehat{\sigma }_i}{}\left[W_i(\sigma _i\widehat{\sigma }_i)P(\{\sigma \},t)W_i(\widehat{\sigma }_i\sigma _i)P(\{\sigma _{ji}\},\widehat{\sigma }_i,t)\right],$$
(1)
The transition probability is in a normalized form determined by a heat Boltzmann factor,

$$W_i(\sigma _i\widehat{\sigma }_i)=\frac{1}{Q_i}\mathrm{exp}\left[\beta _i(\widehat{\sigma }_i,\underset{<ij>}{}\sigma _j)\right],Q_i=\underset{\widehat{\sigma }_i}{}\mathrm{exp}\left[\beta _i(\widehat{\sigma }_i,\underset{<ij>}{}\sigma _j)\right].$$
(2)
One can clearly see that this mechanism favors a lower Hamiltonian of the system. Based on the master equation Eqs.(1), one can prove that the time expectations of single-spin and correlation functions are

$`{\displaystyle \frac{d}{dt}}\sigma _{i_1}(t)\sigma _{i2}(t)\mathrm{}\sigma _{i_n}(t)`$ (3)
$`=`$ $`n\sigma _{i_1}\sigma _{i_2\mathrm{}}\sigma _{i_n}+{\displaystyle \underset{\{\sigma \}}{}}\left\{{\displaystyle \underset{k=1}{\overset{n}{}}}\left[\left({\displaystyle \underset{j(k)=1}{\overset{n}{}}}\sigma _{i_j}\right)\left({\displaystyle \underset{\widehat{\sigma }_{i_k}}{}}\widehat{\sigma }_{i_k}W_{i_k}(\sigma _{i_k}\widehat{\sigma }_{i_k})\right)\right]\right\}P(\{\sigma \},t).`$ (4)
When $`n=1`$, it is

$$\frac{d}{dt}q_k\left(t\right)=q_k\left(t\right)+\underset{\left\{\sigma \right\}}{}\left[\underset{\widehat{\sigma }_k}{}\widehat{\sigma }_kW_k\left(\sigma _k\widehat{\sigma }_k\right)\right]P(\left\{\sigma \right\};t).$$
(5)
### B Spin-pair redistribution mechanism

Kawasaki’s spin-pair exchange mechanism allows an Ising system to evolve with its nearest neighbors exchanging their spin values. In spin-pair redistribution mechanism, two neighboring spins, $`\sigma _j\sigma _l`$, may change to any possible values, $`\widehat{\sigma }_j\widehat{\sigma }_l`$, as long as their sum are conserved, and the master equation is

$`{\displaystyle \frac{d}{dt}}P(\{\sigma \},t)`$ $`=`$ $`{\displaystyle \underset{jl}{}}{\displaystyle \underset{\widehat{\sigma }_j,\widehat{\sigma }_l}{}}[W_{jl}(\sigma _j\sigma _l\widehat{\sigma }_j\widehat{\sigma }_l)P(\left\{\sigma \right\};t)`$ (7)
$`+W_{jl}(\widehat{\sigma }_j\widehat{\sigma }_l\sigma _j\sigma _l)P(\{\sigma _{ij},\sigma _{lk}\},\widehat{\sigma }_j,\widehat{\sigma }_l;t)].`$
The redistribution probability is also in a normalized form determined by a heat Boltzmann factor,

$$W_{jl}\left(\sigma _j\sigma _l\widehat{\sigma }_j\widehat{\sigma }_l\right)=\frac{1}{Q_{jl}}\delta _{\sigma _j+\sigma _l,\widehat{\sigma }_j+\widehat{\sigma }_l}\mathrm{exp}\left[\beta _{jl}(\widehat{\sigma }_j,\widehat{\sigma }_l,\left\{\sigma _m\right\}_{mj,l})\right],$$
(8)
where the normalization factor $`Q_{jl}`$ is

$`Q_{jl}={\displaystyle \underset{\widehat{\sigma }_j,\widehat{\sigma }_l}{}}\delta _{\sigma _j+\sigma _l,\widehat{\sigma }_j+\widehat{\sigma }_l}\mathrm{exp}\left[\beta _{jl}(\widehat{\sigma }_j,\widehat{\sigma }_l,\left\{\sigma _m\right\}_{mj,l})\right]..`$

(Here it clearly favors a lower system Hamiltonian, but in the combined mechanism we shall change the sign before $`\beta _{jl}`$ and make it turn to the opposite.) The time expectation of single-spin is

$$\frac{d}{dt}q_k\left(t\right)=2dq_k\left(t\right)+\underset{w}{}\underset{\left\{\sigma \right\}}{}\left[\underset{\widehat{\sigma }_k,\widehat{\sigma }_{k+w}}{}\widehat{\sigma }_kW_{k,k+w}\left(\sigma _k\sigma _{k+w}\widehat{\sigma }_k\widehat{\sigma }_{k+w}\right)\right]P(\left\{\sigma \right\};t),$$
(9)
where $`d`$ is the space dimensionality and $`_w`$ means a summation taken over the nearest neighbors.

### C The competing mechanism

With the competing mechanisms, single-spin transition with probability $`p`$ and spin-pair redistribution with probability $`1p`$, the master equation can be written as,

$$\frac{d}{dt}P(\left\{\sigma \right\},t)=pG_{me}+\left(1p\right)K_{me},$$
(10)
where the Glauber-type

$$G_{me}=\underset{i}{}\underset{\widehat{\sigma }_i}{}\left[W_i\left(\sigma _i\widehat{\sigma }_i\right)P(\left\{\sigma \right\},t)+W_i\left(\widehat{\sigma }_i\sigma _i\right)P(\left\{\sigma _{ji}\right\},\widehat{\sigma }_i;t)\right],$$
(11)
and the Kawasaki-type

$`K_{me}`$ $`=`$ $`{\displaystyle \underset{jl}{}}{\displaystyle \underset{\widehat{\sigma }_j,\widehat{\sigma }_l}{}}[W_{jl}(\sigma _j\sigma _l\widehat{\sigma }_j\widehat{\sigma }_l)P(\left\{\sigma \right\};t)`$ (13)
$`+W_{jl}(\widehat{\sigma }_j\widehat{\sigma }_l\sigma _j\sigma _l)P(\{\sigma _{ij},\sigma _{lk}\},\widehat{\sigma }_j,\widehat{\sigma }_l;t)].`$
The Glauber-type mechanism is used to simulate the contact of the system with the external heat bath, and the transition probability is of the form given by Eqs.(2). This mechanism favors lower energy of the system. The spin-pair redistribution mechanism is used to simulate the energy flux. The redistribution probability given above favors a lower system Hamiltonian, but what we need here is to the contrary. We can turn this tendency to its opposite if we change the sign before $`\beta _{jl}`$ in the redistribution probability. It has the following form,

$$W_{jl}\left(\sigma _j\sigma _l\widehat{\sigma }_j\widehat{\sigma }_l\right)=\frac{1}{Q_{jl}}\delta _{\sigma _j+\sigma _l,\widehat{\sigma }_j+\widehat{\sigma }_l}\mathrm{exp}\left[\beta _{jl}(\widehat{\sigma }_j,\widehat{\sigma }_l,\left\{\sigma _m\right\}_{mj,l})\right],$$
(14)
$`Q_{jl}={\displaystyle \underset{\widehat{\sigma }_j,\widehat{\sigma }_l}{}}\delta _{\sigma _j+\sigma _l,\widehat{\sigma }_j+\widehat{\sigma }_l}\mathrm{exp}\left[\beta _{jl}(\widehat{\sigma }_j,\widehat{\sigma }_l,\left\{\sigma _m\right\}_{mj,l})\right].`$

This normalized form implies that the tendency is toward a higher system Hamiltonian. The word, competition, is in fact that between the two opposite directions either favored by one mechanism.

As given above, we have already obtained the time expectation of single-spin with either mechanism respectively. One can prove that with competing mechanisms it will be,

$$\frac{d}{dt}q_k\left(t\right)=pQ_k^G+\left(1p\right)Q_k^K,$$
(15)
where the Glauber-type,

$$Q_k^G=q_k\left(t\right)+\underset{\left\{\sigma \right\}}{}\left[\underset{\widehat{\sigma }_k}{}\widehat{\sigma }_kW_k\left(\sigma _k\widehat{\sigma }_k\right)\right]P(\left\{\sigma \right\};t),$$
(16)
and the Kawasaki-type,

$$Q_k^K=2dq_k\left(t\right)+\underset{w}{}\underset{\left\{\sigma \right\}}{}\left[\underset{\widehat{\sigma }_k,\widehat{\sigma }_{k+w}}{}\widehat{\sigma }_kW_{k,k+w}\left(\sigma _k\sigma _{k+w}\widehat{\sigma }_k\widehat{\sigma }_{k+w}\right)\right]P(\left\{\sigma \right\};t).$$
(17)
There are equations available for correlation functions with competing mechanisms, but in later studies we find this single-spin equation is enough to yield satisfying results.

## III On Ising Model

As mentioned above, these years great efforts have been contributed to Ising model with competing dynamics. Because exact analytical treatment is too hard, most of the studies have been either approximation or Monte Carlo simulation. There are some differences between the conventional method and ours. In the following paragraph we present our considerations.

The results one may expect directly depend on the expression of the transition (flipping, exchange, redistribution) probabilities. We think that there are two requirements: First, this probability should contain the Hamiltonian, and thus naturally favors either higher energy (Kawasaki-type) or lower energy (Glauber-type). Second, introducing temperature into it, we require that the transition be influenced by the heat noise. It is the first requirement that makes the two mechanisms compete and in all the studies it has been well adopted. However, due to the difficulties of actual practice, most of them used the simplified versions. In most of the studies, the temperature factor has not been introduced into the exchange probability, while in the flipping probability it has been only partly combined. Typically for ferromagnetic Ising model it has been set as: Glauber-type

$`W_i=\mathrm{min}\{1,\mathrm{exp}\left(\mathrm{}E_i/KT\right)\},`$

Kawasaki-type

$`W_{ij}=\{\begin{array}{c}1,\mathrm{}E_{ij}>0\\ 0,\mathrm{}E_{ij}0\end{array}.`$

In the new mechanism the transition and redistribution probabilities do not take the simplified versions. Besides some mathematical aspects such as normalization, their difference lies in the role of temperature. With the new mechanism there is greater influence of heat on the system. The fact, order for lower temperature and disorder for higher temperature, is not universally true in the phase diagrams obtained in earlier studies. However, we believe that this expectation will be unshakable if the system is governed by the new mechanism. In order to further study it, we applied our method to $`1D`$ Ising model. It is well known that, due to the heat noise, there is only paramagnetic phase in $`1D`$ Ising model. The analytical results we obtained confirms this conclusion, however one increases the energy flux, for all temperatures. (The details are in Appendix A, but we suggest you read it later for your easier understanding of our method.) $`1D`$ Ising model governed by the conventional mechanism has been studied in Ref. with approximation method and MC simulation. Further application on $`2D`$ Ising model is beyond the scope of this article. There is no better or worse, since the two mechanisms have different characteristics, but we think the comparison will be interesting and also feasible in practice.

## IV On the kinetic Gaussian model

In this section we apply the new mechanism to the $`3D`$ kinetic Gaussian model and report our findings with the phase diagram. One and two-dimensional cases are quite similar. First we briefly review some basic properties of this model.

Gaussian model, proposed by T. H. Berlin and M. Kac, at first in order to make Ising model more tractable, is an continuous-spin model. It has the same Hamiltonian as the Ising model (three dimensional),

$$\beta =K\underset{i,j,k=1}{\overset{N}{}}\underset{w}{}\sigma _{ijk}\left(\sigma _{i+w,jk}+\sigma _{ij+w,k}+\sigma _{ij,k+w}\right),$$
(18)
where $`K=J/k_BT`$ and $`_w`$ means summation over nearest neighbors. Compared with the Ising model, it has two extensions: First, the spins $`\sigma _{ijk}`$ can take any real value between $`(\mathrm{},+\mathrm{})`$. Second, to prevent the spins from tending to infinity, the probability of finding a given spin between $`\sigma _{ijk}`$ and $`\sigma _{ijk}+d\sigma _{ijk}`$ is assumed to be the Gaussian-type distribution

$$f\left(\sigma _{ijk}\right)d\sigma _{ijk}=\sqrt{\frac{b}{2\pi }}\mathrm{exp}\left(\frac{b}{2}\sigma _{ijk}^2\right)d\sigma _{ijk},$$
(19)
where $`b`$ is a distribution constant independent of temperature. Being an extension of Ising model, Gaussian model is quite different however. In the equilibrium case, on translational invariant lattices it is exactly solvable, and later as a starting point to study the unsolvable models it has also been investigated with mean field theory and the momentum-space renormalization-group method.

As an example of the applications of single-spin transition and spin-pair redistribution mechanism, we have studied kinetic Gaussian model with both of them respectively. We summarize its dynamic properties as follows. The inherent dynamical competition of this model is that: the system tries to lower its Hamiltonian with the spins tending to infinity, while the Gaussian-type probability serves to restrict this tendency. Above the critical temperature, the prevailing heat noise permits only a disordered state, whereas below the critical point some kind of order will appear. Our study shows that, on translational-invariant lattices, the dynamic critical exponent $`z=1/v=2`$ is independent of space dimensionality and the governing dynamical mechanism.

Now we turn to treat the $`3D`$ kinetic Gaussian model with the competing dynamics. $`1D`$ and $`2D`$ systems can be treated in a similar way and they have qualitatively the same properties. In earlier studies we obtained the time expectation of single-spin. With the competing mechanism we can borrow these equations from Ref. and .

(1) With Glauber dynamics:

$$\frac{d}{dt}q_{ijk}\left(t\right)Q_{ijk}^G=q_{ijk}+\frac{K}{b}\left(q_{i1,j,k}+q_{i+1,j,k}+q_{i,j1,k}+q_{i,j+1,k}+q_{i,j,k1}+q_{i,j,k+1}\right).$$
(20)
(2) With Kawasaki dynamics: since we have changed the sign before system Hamiltonian in the redistribution probability,

$`\beta =K{\displaystyle \underset{i,j,k=1}{\overset{N}{}}}{\displaystyle \underset{w}{}}\sigma _{ijk}\left(\sigma _{i+w,jk}+\sigma _{ij+w,k}+\sigma _{ij,k+w}\right),`$

$`\beta =\left(K\right){\displaystyle \underset{i,j,k=1}{\overset{N}{}}}{\displaystyle \underset{w}{}}\sigma _{ijk}\left(\sigma _{i+w,jk}+\sigma _{ij+w,k}+\sigma _{ij,k+w}\right),`$

all the expressions will remain the same if we switch $`K`$ to $`K`$. Here we will have to do the same,

$`{\displaystyle \frac{d}{dt}}q_{ijk}\left(t\right)`$ $``$ $`Q_{ijk}^K={\displaystyle \frac{1}{2\left[b+\left(K\right)\right]}}b\{[(q_{i+1,j,k}q_{ijk})(q_{ijk}q_{i1,j,k})]`$ (28)
$`+[(q_{i,j+1,k}q_{ijk})(q_{ijk}q_{i,j1,k})]+[(q_{i,j,k+1}q_{ijk})(q_{ijk}q_{i,j,k1})]\}`$
$`+{\displaystyle \frac{\left(K\right)}{2\left[b+\left(K\right)\right]}}[2(2q_{i1,j,k}q_{i1,j+1,k}q_{i1,j1,k})+(2q_{i1,j,k}q_{ijk}q_{i2,j,k})`$
$`+2\left(2q_{i+1,j,k}q_{i+1,j+1,k}q_{i+1,j1,k}\right)+\left(2q_{i+1,j,k}q_{ijk}q_{i+2,j,k}\right)`$
$`+2\left(2q_{i,j1,k}q_{i,j1,k+1}q_{i,j1,k1}\right)+\left(2q_{i,j1,k}q_{ijk}q_{i,j2,k}\right)`$
$`+2\left(2q_{i,j+1,k}q_{i,j+1,k+1}q_{i,j+1,k1}\right)+\left(2q_{i,j+1,k}q_{ijk}q_{i,j+2,k}\right)`$
$`+2\left(2q_{i,j,k1}q_{i1,j,k1}q_{i+1,j,k1}\right)+\left(2q_{i,j,k1}q_{ijk}q_{i,j,k2}\right)`$
$`+2(2q_{i,j,k+1}q_{i+1,j,k+1}q_{i1,j,k+1})+(2q_{i,j,k+1}q_{ijk}q_{i,j,k+2})].`$
This expression is in fact not as complex as it seems. Each term in a bracket is a second order derivative of $`q`$, and they will cancel each other if a summation is taken over all the spins. Thus $`_{ijk}Q_{ijk}^K=0`$.

Then with the competing mechanisms,

$$\frac{d}{dt}q_{ijk}\left(t\right)=pQ_{ijk}^G+\left(1p\right)Q_{ijk}^K.$$
(29)
It is natural to first consider using the system Hamiltonian,

$`=J{\displaystyle \underset{i,j}{}}\sigma _i\sigma _j,`$

to characterize its behavior. If we find the Hamiltonian is decreasing, then the system is believed to be evolving to Ferromagnetic phase; if the Hamiltonian is stable around zero, the system must be in Paramagnetic phase; if the Hamiltonian is increasing, then the system is evolving towards Antiferromagnetic phase. Unfortunately, in Gaussian model although we can conveniently obtain the exact result for $`\sigma _i\sigma _j`$ and the average value, $`_{ij}\sigma _i\sigma _j/N^2`$, we find it difficult to obtain an analytical result of the system Hamiltonian. In order to differentiate between these phases, instead we study the following two aspects. First, magnetization,

$$M\left(t\right)=\frac{1}{N}\underset{k}{}q_k\left(t\right),$$
(30)
and second,

$$M^{}\left(t\right)=\frac{1}{N}\underset{ijk}{}q_{ijk}^{}\left(t\right)=\frac{1}{N}\underset{ijk}{}()^{i+j+k}q_{ijk}\left(t\right).$$
(31)
For pure ferromagnetic phase, if without any extra conditions, we expect nonzero magnetization, that is to say, $`M\left(t\right)0`$. At the same time, if the system is divided into two penetrating subsets, in one of which the spins have $`i+j+k`$ being an odd number, and in the other one the spins have $`i+j+k`$ being even, then one will find these two subsets are almost identical, and that leads to $`M^{}\left(t\right)=0`$. So $`M\left(t\right)0`$ and $`M^{}\left(t\right)=0`$ is the characteristic of ferromagnetic phase. For pure antiferromagnetic phase, we expect such a situation that the system can be divided into two penetrating opposing sublattices in the way mentioned above, one of positive spin and one of negative. Except the direction of the spins, these two sublattices are identical. This leads to $`M^{}\left(t\right)0`$ but $`M\left(t\right)=0`$. In disordered paramagnetic phase, we expect disorder of the whole system and both $`M\left(t\right)`$ and $`M^{}\left(t\right)`$ to be zero. Thus, if we get the evolution of $`M\left(t\right)`$ and $`M^{}\left(t\right)`$, we can decide in which phase the system is.

With Eqs.(29), (30) and (31) we can write the time derivative of $`M\left(t\right)`$ and $`M^{}\left(t\right)`$ as,

$`{\displaystyle \frac{dM\left(t\right)}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{k}{}}{\displaystyle \frac{d}{dt}}q_k\left(t\right)`$ (32)
$`=`$ $`p\left(1+{\displaystyle \frac{6K}{b}}\right)M\left(t\right).`$ (33)
and

$`{\displaystyle \frac{dM^{}\left(t\right)}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{k}{}}{\displaystyle \frac{d}{dt}}q_k^{}\left(t\right)`$ (34)
$`=`$ $`\left[p\left(1+{\displaystyle \frac{6K}{b}}\right)+6\left(1p\right){\displaystyle \frac{6Kb}{bK}}\right]M^{}\left(t\right)`$ (35)
The details of Eqs.(35) are in Appendix B. Because we want $`K/b`$ to be positive, the inequality

$`p\left(1+{\displaystyle \frac{6K}{b}}\right)+6\left(1p\right){\displaystyle \frac{6Kb}{bK}}>0,`$

becomes

$$\frac{1}{12p}\left(36+41p+\sqrt{\left(12962808p+1561p^2\right)}\right)<\frac{K}{b}<1,$$
(36)
for $`0<p<1`$.

Fig.1 shows the phase diagram of $`3D`$ kinetic Gaussian model. It is divided into several regions by the following three curves,

$$\mathrm{exp}\left(K/b\right)=e^{1/6},$$
(37)
$$\mathrm{exp}\left(K/b\right)=1/e,$$
(38)
and

$$\mathrm{exp}\left(K/b\right)=\mathrm{exp}\left[\frac{1}{12p}\left(36+41p+\sqrt{\left(12962808p+1561p^2\right)}\right)\right].$$
(39)
Phase Diagrams for $`1D`$ and $`2D`$ models are quite similar.

We know that the critical point of $`3D`$ Gaussian model is $`K_c=b/2d`$, where $`d`$ is the space dimensionality. In the region above line (37), the temperature is higher than the critical value. Both $`M\left(t\right)`$ and $`M^{}\left(t\right)`$ are approaching zero exponentially, and this is identified as paramagnetic phase. The overwhelming heat noise permits no observable magnetization, and the Hamiltonian is static in equilibrium. When the temperature is below $`T_c`$, the heat noise becomes secondary and some kind of order appears.

For the region below line (37) and above curve (39), as well as that below line (38), we have exponentially increasing $`M\left(t\right)`$ while $`M^{}\left(t\right)`$ is approaching zero, and this corresponds to ferromagnetic phase. In this phase, Glauber-type mechanism prevails and the energy is decreasing. There is observable homogeneous magnetization, the direction of which depends on that of the initial magnetization.

In the region between curve (38) and line (39) we have both exponentially increasing $`M\left(t\right)`$ and $`M^{}\left(t\right)`$, and it can not be simply identified as ferromagnetic or antiferromagnetic. We give it a name, heterophase. In this region, the Kawasaki mechanism and the energy flux control the system, and the energy has a tendency to increase. The spin values are interesting. Our analytical results show that $`M\left(t\right)=M\left(0\right)e^{At}`$ and $`M^{}\left(t\right)=M^{}\left(0\right)e^{Bt}`$. This means that, if the initial phase is ferromagnetic, $`M^{}\left(0\right)=0`$, but $`M\left(0\right)0`$, later we will observe $`\left|M\left(t\right)\right|`$ increasing but $`M^{}\left(t\right)`$ staying at zero. If initially the system is antiferromagnetic, $`M\left(0\right)=0`$, but $`M^{}\left(0\right)0`$, later we will have $`\left|M^{}\left(t\right)\right|`$ increasing but $`M\left(t\right)`$ staying at zero. If initially a disordered paramagnetic phase is given, $`M\left(0\right)=0`$, and also $`M^{}\left(0\right)=0`$, we will get both zero $`M\left(t\right)`$ and $`M^{}\left(t\right)`$. From this one may get confused,—what is the real picture? How is the energy sure to increase if the spin values are dependent on other conditions? Direct computer simulation reveals the key. For example, we apply the periodic boundary condition and initially set the system as,

$`\mathrm{}+1,+1,1,+1,+1,1,+1,+1,1,+1,`$
$`+1,1,+1,+1,1,+1,+1,1,+1,+1,\mathrm{}`$
This leads to $`M\left(0\right)0`$, but $`M^{}\left(0\right)=0`$. After $`1`$ second, it becomes (we only give approximate values),

$`\mathrm{}0.7,1.1,0.2,0.8,0.6,0.3,0.6,0.6,0.2,0.6,`$
$`0.6,0.2,0.6,0.6,0.3,0.6,0.8,0.2,1.1,0.7,\mathrm{}`$
After $`10`$ seconds,

$`\mathrm{}9,13,5,17,3,17,3,15,6,11,`$
$`11,6,15,3,17,3,17,5,13,9,\mathrm{}`$
After $`24`$ seconds,

$`\mathrm{}249,6941,6735,12089,10041,13228,8902,9922,3755,3435,`$
$`3435,3755,9922,8902,13228,10041,12089,6735,6941,249,\mathrm{}`$
This shows the routine of self-organization the system chooses in this specific case, and it is a self-explanatory picture of how the energy manages to increase while $`M^{}\left(t\right)=0`$.

However, there is one exception. If at the beginning the system is set in a pure homogeneous ferromagnetic phase, we will find that the system remains in this homogeneous state. How can one explain this? Actually the system does have a tendency to increase its energy in this region, but it needs a hint to know exactly how. Here we cite the words of Tome and Oliveira: It is ”the result, in the case considered here, of a far from equilibrium process, namely, the continuous flux of energy into the system. Thus an instability of the usual (equilibrium) solutions leads the system toward states with spatial self-organized structure.” The self-organization routine is closely related with a basic phenomenon of symmetry loss. A system under conditions that lead to paramagnetic phase is like a ball placed on the ground, for which all directions are identical. A system set in such a heterophase is like a ball placed on the top of another sphere. The upper ball has a tendency to fall down, but the direction depends on some kind of disturbance. If there is no disturbance at all, it will stay still. It will even have its energy increasing if the lower ball is being lifted up. Surely in reality we have never seen a ball stay still on the top of another, because disturbance can not be avoided. So discussion here is only to help qualitatively explain this strange phase and the mechanism hidden behind.

## V Results and some Discussion

In this article, we have given a systematic formulation of the new competing mechanism. The Glauber-type single-spin transition mechanism with probability $`p`$ simulates the contact of the system with the heat bath and the Kawasaki-type spin-pair redistribution mechanism with probability $`1p`$ simulates an external energy flux. These two mechanisms themselves are natural generalizations of Glauber’s single-spin flipping mechanism and Kawasaki’s spin-pair exchange mechanism. Thus, on the one hand, this mechanism is in principle applicable to arbitrary systems, while on the other hand, our formulation is able to contain a mechanism that just directly combines single-spin flipping and spin-pair exchange in their original forms (not simplified). Compared with the conventional one, the new mechanism does not assume the simplified versions. Their difference lies in the different role the system temperature plays. (As we have emphasized before, there is no better or worse, since they are different.) We applied the new mechanism to $`1D`$ Ising model and used the analytical results to exemplify this difference. With the new mechanism, there is greater influence of temperature and the fact, order for lower temperature and disorder for higher temperature, will be universally true.

In Sec. IV, we applied this mechanism to $`3D`$ kinetic Gaussian model. $`1D`$ and $`2D`$ models can be treated following a similar way and have qualitatively the same properties. We used $`M\left(t\right)`$ and $`M^{}\left(t\right)`$ (their definition given in that section) to characterize the system. In the phase diagram of $`3D`$ model, we confirm the expectation, order for lower temperature and disorder for higher temperature. For temperature above $`T_c`$, the system evolves to paramagnetic phase. For temperature lower than $`T_c`$, we observe ferromagnetic phase and another heterophase (instead of the antiferromagnetic phase as guessed). This interesting heterophase is the result of the energy flux and self-organization. We have analyzed it in details in that section, and hope it will help make clear the self-organization process in phase transitions. With regard to Ising model, we also hope it will be a good reference when people are trying to solve the puzzling differences between the results yielded by approximation method and MC simulation.

## ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China under Grant No. 10075025.

## A Calculational details on the kinetic Ising model

We apply the new mechanism to $`1D`$ Ising model. We derive the time expectation of single-spin using Eqs.(15),

$`{\displaystyle \frac{d}{dt}}q_k\left(t\right)=pQ_k^G+\left(1p\right)Q_k^K.`$

In this combined form, the Glauber-type term is as Eqs.(16),

$$Q_k^G=q_k\left(t\right)+\underset{\left\{\sigma \right\}}{}\left[\underset{\widehat{\sigma }_k}{}\widehat{\sigma }_kW_k\left(\sigma _k\widehat{\sigma }_k\right)\right]P(\left\{\sigma \right\};t),$$
(A1)
and the Kawasaki-type is as Eqs.(17),

$`Q_k^K`$ $`=`$ $`2q_k\left(t\right)+{\displaystyle \underset{\left\{\sigma \right\}}{}}[{\displaystyle \underset{\widehat{\sigma }_k,\widehat{\sigma }_{k+1}}{}}\widehat{\sigma }_kW_{k,k+1}(\sigma _k\sigma _{k+1}\widehat{\sigma }_k\widehat{\sigma }_{k+1})`$ (A3)
$`+{\displaystyle \underset{\widehat{\sigma }_k,\widehat{\sigma }_{k1}}{}}\widehat{\sigma }_kW_{k,k1}(\sigma _k\sigma _{k1}\widehat{\sigma }_k\widehat{\sigma }_{k1})\left]P\right(\left\{\sigma \right\};t).`$
First we calculate (A1).

$`{\displaystyle \underset{\widehat{\sigma }_k}{}}\widehat{\sigma }_kW_k\left(\sigma _k\widehat{\sigma }_k\right)`$ (A4)
$`=`$ $`{\displaystyle \frac{\sigma _k\mathrm{exp}\left[K\sigma _k\left(\sigma _{k1}+\sigma _{k+1}\right)\right]\sigma _k\mathrm{exp}\left[K\sigma _k\left(\sigma _{k1}+\sigma _{k+1}\right)\right]}{\mathrm{exp}\left[K\sigma _k\left(\sigma _{k1}+\sigma _{k+1}\right)\right]+\mathrm{exp}\left[K\sigma _k\left(\sigma _{k1}+\sigma _{k+1}\right)\right]}}`$ (A5)
$`=`$ $`\{\begin{array}{c}0,\sigma _{k1}=\sigma _{k+1}\\ \mathrm{tanh}2K,\sigma _{k1}=\sigma _{k+1}=1\\ \mathrm{tanh}2K,\sigma _{k1}=\sigma _{k+1}=1\end{array}`$ (A9)
$`=`$ $`{\displaystyle \frac{1}{2}}\left(\sigma _{k1}+\sigma _{k+1}\right)\mathrm{tanh}2K.`$ (A10)
So

$`Q_G^k`$ $`=`$ $`q_k\left(t\right)+{\displaystyle \underset{\left\{\sigma \right\}}{}}\left[{\displaystyle \underset{\widehat{\sigma }_k}{}}\widehat{\sigma }_kW_k\left(\sigma _k\widehat{\sigma }_k\right)\right]P(\left\{\sigma \right\};t)`$ (A11)
$`=`$ $`q_k\left(t\right)+{\displaystyle \frac{1}{2}}\left(q_{k1}+q_{k+1}\right)\mathrm{tanh}2K.`$ (A12)
Second, we calculate (A3),

$`{\displaystyle \underset{\widehat{\sigma }_k,\widehat{\sigma }_{k+1}}{}}\widehat{\sigma }_kW_{k,k+1}\left(\sigma _k\sigma _{k+1}\widehat{\sigma }_k\widehat{\sigma }_{k+1}\right)`$
$`=`$ $`{\displaystyle \frac{\sigma _k\mathrm{exp}\left[K\left(\sigma _{k1}\sigma _k+\sigma _{k+1}\sigma _{k+2}\right)\right]+\sigma _{k+1}\mathrm{exp}\left[K\left(\sigma _{k1}\sigma _{k+1}+\sigma _k\sigma _{k+2}\right)\right]}{\mathrm{exp}\left[K\left(\sigma _{k1}\sigma _k+\sigma _{k+1}\sigma _{k+2}\right)\right]+\mathrm{exp}\left[K\left(\sigma _{k1}\sigma _{k+1}+\sigma _k\sigma _{k+2}\right)\right]}}`$
$`=`$ $`\{\begin{array}{c}\left(\sigma _k+\sigma _{k+1}\right)/2,\sigma _k=\sigma _{k+1}\\ \mathrm{tanh}\left[K\left(\sigma _{k1}\sigma _{k+2}\right)\right],\sigma _k=\sigma _{k+1}\end{array}`$
$`=`$ $`{\displaystyle \frac{\sigma _k+\sigma _{k+1}}{2}}+{\displaystyle \frac{\left(\sigma _k\sigma _{k+1}\right)^2}{4}}`$
$`\times \left[{\displaystyle \frac{\left(\sigma _{k1}1\right)\left(\sigma _{k+2}+1\right)}{4}}+{\displaystyle \frac{\left(\sigma _{k1}+1\right)\left(\sigma _{k+2}1\right)}{4}}\right]\mathrm{tanh}\left(2K\right).`$
$`=`$ $`{\displaystyle \frac{\sigma _k+\sigma _{k+1}}{2}}+{\displaystyle \frac{1}{4}}\left(\sigma _{k1}\sigma _{k+2}\sigma _{k1}\sigma _k\sigma _{k+1}+\sigma _k\sigma _{k+1}\sigma _{k+2}\right)\mathrm{tanh}\left(2K\right)`$
And similarly,

$`{\displaystyle \underset{\widehat{\sigma }_k,\widehat{\sigma }_{k1}}{}}\widehat{\sigma }_kW_{k,k1}\left(\sigma _k\sigma _{k1}\widehat{\sigma }_k\widehat{\sigma }_{k1}\right)`$
$`=`$ $`{\displaystyle \frac{\sigma _k+\sigma _{k1}}{2}}+{\displaystyle \frac{1}{4}}\left(\sigma _{k+1}\sigma _{k2}\sigma _{k+1}\sigma _k\sigma _{k1}+\sigma _k\sigma _{k1}\sigma _{k2}\right)\mathrm{tanh}\left(2K\right).`$
So

$`Q_k^K`$ $`=`$ $`2q_k\left(t\right)+{\displaystyle \underset{\left\{\sigma \right\}}{}}[{\displaystyle \underset{\widehat{\sigma }_k,\widehat{\sigma }_{k+1}}{}}\widehat{\sigma }_kW_{k,k+1}(\sigma _k\sigma _{k+1}\widehat{\sigma }_k\widehat{\sigma }_{k+1})`$ (A15)
$`+{\displaystyle \underset{\widehat{\sigma }_k,\widehat{\sigma }_{k1}}{}}\widehat{\sigma }_kW_{k,k1}(\sigma _k\sigma _{k1}\widehat{\sigma }_k\widehat{\sigma }_{k1})\left]P\right(\left\{\sigma \right\};t)`$
$`=`$ $`{\displaystyle \frac{1}{2}}\left(q_{k1}+q_{k+1}2q_k\right)+{\displaystyle \frac{1}{4}}\mathrm{tanh}\left(2K\right)\left(q_{k+1}q_{k2}+q_{k1}q_{k+2}\right)`$ (A17)
$`+{\displaystyle \frac{1}{4}}\mathrm{tanh}\left(2K\right)\left(\sigma _k\sigma _{k+1}\sigma _{k+2}2\sigma _{k+1}\sigma _k\sigma _{k1}+\sigma _k\sigma _{k1}\sigma _{k2}\right).`$
Thus, combining (A12) and (A17) one will get,

$`{\displaystyle \frac{d}{dt}}q_k\left(t\right)`$ $`=`$ $`p\left[q_k\left(t\right)+{\displaystyle \frac{1}{2}}\left(q_{k1}+q_{k+1}\right)\mathrm{tanh}2K\right]`$ (A20)
$`+(1p)[{\displaystyle \frac{1}{2}}(q_{k1}+q_{k+1}2q_k)+{\displaystyle \frac{1}{4}}\mathrm{tanh}(2K)(q_{k+1}q_{k2}+q_{k1}q_{k+2})`$
$`+{\displaystyle \frac{1}{4}}\mathrm{tanh}(2K)(\sigma _k\sigma _{k+1}\sigma _{k+2}2\sigma _{k+1}\sigma _k\sigma _{k1}+\sigma _k\sigma _{k1}\sigma _{k2})].`$
In order to decide in which phase the system is, we have suggested two quantities in Sec. IV,

$`M\left(t\right)={\displaystyle \underset{k}{}}q_k\left(t\right),`$

and

$`M^{}\left(t\right)={\displaystyle \underset{k}{}}()^kq_k\left(t\right){\displaystyle \underset{k}{}}q_k^{}\left(t\right).`$

We use their evolving tendency to characterize the system behavior. Obviously, (1) in a homogeneous ferromagnetic phase, we expect $`\left|M\left(t\right)\right|/N1`$ and $`\left|M^{}\left(t\right)\right|/N0`$; (2) in an antiferromagnetic phase consisting of two penetrating and opposing sublattices, we will have $`\left|M^{}\left(t\right)\right|/N1`$ and $`\left|M\left(t\right)\right|/N0`$; (3) in disordered paramagnetic phase, both $`M^{}\left(t\right)/N`$ and $`M\left(t\right)/N`$ will approach zero. (A detailed analysis can be found in that section.) With Eqs.(A20) we can get,

$`{\displaystyle \frac{d}{dt}}M\left(t\right){\displaystyle \underset{k}{}}{\displaystyle \frac{d}{dt}}q_k\left(t\right)=p\left(1\mathrm{tanh}2K\right)M\left(t\right).`$

Thus

$$M\left(t\right)=M\left(0\right)\mathrm{exp}\left[p\left(1\mathrm{tanh}2K\right)t\right].$$
(A21)
At the same time,

$`{\displaystyle \frac{d}{dt}}()^kq_k\left(t\right)`$
$``$ $`{\displaystyle \frac{d}{dt}}q_k^{}\left(t\right)=p\left[q_k^{}\left(t\right)+{\displaystyle \frac{1}{2}}\left(q_{k1}^{}q_{k+1}^{}\right)\mathrm{tanh}2K\right]`$
$`+(1p)[{\displaystyle \frac{1}{2}}(q_{k1}^{}q_{k+1}^{}2q_k^{})+{\displaystyle \frac{1}{4}}\mathrm{tanh}(2K)(q_{k+1}^{}q_{k2}^{}q_{k1}^{}q_{k+2}^{})`$
$`+{\displaystyle \frac{1}{4}}\mathrm{tanh}(2K)(\sigma _k^{}\sigma _{k+1}^{}\sigma _{k+2}^{}+2\sigma _{k+1}^{}\sigma _k^{}\sigma _{k1}^{}\sigma _k^{}\sigma _{k1}^{}\sigma _{k2}^{})],`$
and

$`M^{}\left(t\right)`$ $``$ $`{\displaystyle \underset{k}{}}{\displaystyle \frac{d}{dt}}q_k^{}\left(t\right)=\left[p\left(1\mathrm{tanh}2K\right)+\left(1p\right)\left(2+\mathrm{tanh}2K\right)\right]M^{}\left(t\right)`$
$`=`$ $`\left[2+p+\left(12p\right)\mathrm{tanh}2K\right]M^{}\left(t\right).`$
Thus

$$M^{}\left(t\right)=M^{}\left(0\right)\mathrm{exp}\left\{\left[2+p+\left(12p\right)\mathrm{tanh}2K\right]t\right\}.$$
(A22)
In $`1D`$ Ising model, as given in Eqs.(A21) and (A22), $`M\left(t\right)`$ and $`M^{}\left(t\right)`$ are both approaching zero exponentially, and we can make the conclusion that, however one increases the energy flux, the system will stay in paramagnetic phase at arbitrary temperature.

## B Calculational details of Eqs.(35)

Kawasaki-type:

$`{\displaystyle \frac{d}{dt}}q_{ijk}^{}\left(t\right)`$ $`=`$ $`{\displaystyle \frac{1}{2\left(bK\right)}}b\{[(q_{i,j,k+1}^{}q_{ijk}^{})(q_{ijk}^{}+q_{i,j,k1}^{})]`$ (B8)
$`+[(q_{i+1,j,k}^{}q_{ijk}^{})(q_{ijk}^{}+q_{i1,j,k}^{})]+[(q_{i,j+1,k}^{}q_{ijk}^{})(q_{ijk}^{}+q_{i,j1,k}^{})]\}`$
$`+{\displaystyle \frac{K}{2\left(bK\right)}}[2(2q_{i1,j,k}^{}q_{i1,j+1,k}^{}q_{i1,j1,k}^{})+(2q_{i1,j,k}^{}q_{ijk}^{}q_{i2,j,k}^{})`$
$`+2\left(2q_{i+1,j,k}^{}q_{i+1,j+1,k}^{}q_{i+1,j1,k}^{}\right)+\left(2q_{i+1,j,k}^{}q_{ijk}^{}q_{i+2,j,k}^{}\right)`$
$`+2\left(2q_{i,j1,k}^{}q_{i,j1,k+1}^{}q_{i,j1,k1}^{}\right)+\left(2q_{i,j1,k}^{}q_{ijk}^{}q_{i,j2,k}^{}\right)`$
$`+2\left(2q_{i,j+1,k}^{}q_{i,j+1,k+1}^{}q_{i,j+1,k1}^{}\right)+\left(2q_{i,j+1,k}^{}q_{ijk}^{}q_{i,j+2,k}^{}\right)`$
$`+2\left(2q_{i,j,k1}^{}q_{i1,j,k1}^{}q_{i+1,j,k+1}^{}\right)+\left(2q_{i,j,k1}^{}q_{ijk}^{}q_{i,j,k2}^{}\right)`$
$`+2(2q_{i,j,k+1}^{}q_{i+1,j,k+1}^{}q_{i1,j,k+1}^{})+(2q_{i,j,k+1}^{}q_{ijk}^{}q_{i,j,k+2}^{})].`$
Glauber-type:

$$\frac{d}{dt}q_{ijk}^{}\left(t\right)=q_{ijk}^{}+\frac{K}{b}\left(q_{i1,j,k}^{}q_{i+1,j,k}^{}q_{i,j+1,k}^{}q_{i,j1,k}^{}q_{i,j,k+1}^{}q_{i,j,k+1}^{}\right).$$
(B9)
And then

$`{\displaystyle \frac{dM^{}\left(t\right)}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{k}{}}{\displaystyle \frac{d}{dt}}q_k^{}\left(t\right)`$
$`=`$ $`\left[p\left(1+{\displaystyle \frac{6K}{b}}\right)+\left(1p\right)6{\displaystyle \frac{6Kb}{bk}}\right]M^{}\left(t\right).`$
Caption of figures

Fig.1. The phase diagram of $`3D`$ Gaussian model with competing dynamics, Glauber-type with probability $`p`$, and Kawasaki-type with probability $`1p`$. The system exhibits paramagnetic phase (Para), Ferromagnetic phase (Ferro), and heterophase (Hetero).


text:
# Vacuum polarization in the spacetime of charged nonlinear black hole

## I Introduction

One of the most intriguing open questions in modern theoretical physics is the issue of a final stage of the black hole evaporation. Although a definite answer would be possible only with the full machinery of quantum theory of gravity, some important preliminary results may be, in principle, obtained within the framework of semiclassical theory. However, as the semiclassical approach could not be used to describe the evolution of the system in the Planck regime, it can, at best, tell us about a tendency rather than the limit itself. Unfortunately, even this simplified programme is hard to execute as the semiclassical Einstein field equations comprise a rather complicated set of nonlinear partial differential equations, and, moreover, the source term — the renormalized stress-energy tensor — should be known for a wide class of nonstatic metrics. It is natural, therefore, that in order to make the back reaction problem tractable, one should refer to some approximations.

It seems that the most promising approach consists in constructing the approximate stress-energy tensor in the geometry of a static black hole and subsequent computation of the semiclassical corrections to the classical metric. Although evaluated in the static background, such corrections, as has been pointed out in Ref. , are relevant because they give direct information about the influence of the quantum effects on the temperature of a black hole. Thus far this programme, initiated in Ref. , has been carried out for massless fields in Schwarzschild spacetime  and for massive scalar fields with arbitrary curvature coupling in Reissner-Nordström (RN) geometry , where, among other things, quantum corrections to the geometry, entropy, and trace anomaly were computed. The most important ingredient of the approach is therefore the renormalized stress-energy tensor of the quantized field propagating in a spacetime of the static black hole constructed in a physically interesting state. However, direct evaluation of this object is rather complicated.

It is believed that the physical content of the theory of quantum fields in curved spacetime is contained in the (renormalized) effective action, $`W_R,`$ a useful quantity allowing evaluation of the stress-energy tensor by means of the standard formula

$$\frac{2}{g^{1/2}}\frac{\delta }{\delta g_{ab}}W_{ren}=T^{ab}_{ren}.$$
(1)
Unfortunately, the effective action is a nonlocal functional of the metric and its exact form is unknown. In the attempts to construct the renormalized stress-energy tensor one is forced, therefore, to employ numerical methods or to accept some approximations.

For the quantized massive fields in the large mass limit, i.e. when the Compton length is much smaller than the characteristic radius of a curvature, the nonlocal contribution to the effective action can be neglected, and the series expansion in $`m^2`$ of the renormalized effective action, $`W_R,`$ may be constructed with the aid of the DeWitt-Schwinger method . Thus constructed renormalized effective action has the form

$$W_{ren}=\frac{1}{32\pi ^2m^2}d^4xg^{1/2}\underset{n=3}{\overset{\mathrm{}}{}}\frac{(n3)!}{(m^2)^{n2}}[a_n(x,x^{})],$$
(2)
where $`[a_n(x,x^{})]`$ is a coincidence limit of the n-th Hadamard-Minakshisundaram-DeWitt-Seely (HMDS) coefficient. The coefficients are local, geometrical objects constructed from curvature tensor and its covariant derivatives of rapidly growing complexity. So far, only coefficients for $`n4`$ are known, but it seems the effective action constructed from higher order terms would be practically intractable, especially in the attempts to calculate the stress-energy tensor.

The first nonvanishing term of the effective action constructed for the massive scalar field with arbitrary curvature coupling satisfying

$$\left(\mathrm{}+\xi R+m^2\right)\varphi =\mathrm{\hspace{0.17em}0},$$
(3)
where $`\xi `$ is the coupling constant and $`m`$ is the mass of the field, form the coincidence limit of the HMDS coefficient $`\left[a_3\right],`$ is given by :

$`W_{ren}^{(1)}`$ $`=`$ $`{\displaystyle \frac{1}{192\pi ^2m^2}}{\displaystyle }d^4xg^{1/2}[{\displaystyle \frac{1}{2}}(\eta ^2{\displaystyle \frac{\eta }{15}}{\displaystyle \frac{1}{315}})R\mathrm{}R+{\displaystyle \frac{1}{140}}R_{pq}\mathrm{}R^{pq}`$ (7)
$`\eta ^3R^3+{\displaystyle \frac{1}{30}}\eta RR_{pq}R^{pq}{\displaystyle \frac{1}{30}}\eta RR_{pqab}R^{pqab}`$
$`{\displaystyle \frac{8}{945}}R_q^pR_a^qR_p^a+{\displaystyle \frac{2}{315}}R^{pq}R_{ab}R_{pq}^{ab}+{\displaystyle \frac{1}{1260}}R_{pq}R_{cab}^pR^{qcab}`$
$`+{\displaystyle \frac{17}{7560}}R_{ab}^{}{}_{}{}^{pq}R_{pq}^{}{}_{}{}^{cd}R_{cd}^{}{}_{}{}^{ab}{\displaystyle \frac{1}{270}}R_{pq}^{ab}R_{cd}^{pq}R_{ab}^{cd}]`$
$`=`$ $`{\displaystyle \frac{1}{192\pi ^2m^2}}{\displaystyle \underset{i=1}{\overset{10}{}}}\alpha _iW_{\left(i\right)},`$ (8)
where $`\eta =\xi 1/6`$ and $`\alpha _{i\text{ }}`$ are numerical coefficients that stand in front of geometric terms in $`W_{ren}^{(1)}.`$ Differentiating functionally $`W_{ren}^{\left(1\right)}`$ with respect to a metric tensor one obtains

$`T^{ab}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{10}{}}}\alpha _i\stackrel{~}{T}^{\left(i\right)ab}={\displaystyle \frac{1}{96\pi ^2m^2g^{1/2}}}{\displaystyle \underset{i=1}{\overset{10}{}}}\alpha _i{\displaystyle \frac{\delta W_{\left(i\right)}}{\delta g_{ab}}}`$ (9)
$`=`$ $`T^{\left(0\right)ab}+\eta T^{\left(1\right)ab}+\eta ^2T^{\left(2\right)ab}+\eta ^3T^{\left(3\right)ab},`$ (10)
where each $`\stackrel{~}{T}^{\left(i\right)ab}`$ is rather complicated expression constructed from the curvature tensor, its covariant derivatives, and contractions. Such calculations have been undertaken in Refs. , where generic expressions for the first nonvanishing order of the renormalized stress-energy tensor were obtained. They generalize earlier results of Frolov and Zel’nikov for vacuum type-D geometries . One can easily extend this result to the case of spinor and vector field as the analogous expressions differ only by numerical coefficients $`\alpha _i`$.

It should be emphasised, however, that the above assumptions place severe limitations on the domain of validity of the obtained approximation. Especially, it would be meaningless, at least in this formulation, to consider the massless limit of the approach. Therefore, the result, which consists of approximately one hundred local geometrical terms, may be used in any spacetime provided the temporal changes of the background geometry are slow and the mass of the field is sufficiently large. Because of the complexity of the thus obtained tensor it will be not presented here, and for its full form and technical details the reader is referred to .

For quantized massive scalar fields with arbitrary curvature coupling in static and spherically-symmetric geometries there exists a different method invented by Anderson, Hiscock, and Samuel . Their calculations were based on the WKB approximation of the solutions of the radial equation and summation of thus obtained mode functions. Both methods are equivalent: to obtain the term proportional to $`m^2`$ one has to use 6-th order WKB, and explicit calculations carried out for the Reissner-Nordström spacetime as well as in wormhole geometries considered in Ref.  yielded identical results. Moreover, detailed numerical analyses carried out in Ref.  and briefly reported in  show that for $`mM2`$ ($`M`$ being the black hole mass), the accuracy of the Schwinger-DeWitt approximation in the Reissner-Nordström geometry is quite good (1% or better).

Geometry of the Reissner-Nordström black hole is singular as $`r0`$ and there is a natural desire to construct its regular generalizations. It is expected that a good candidate for the source term of the Einstein equations is the (classical) stress-energy tensor of nonlinear electrodynamics. Moreover, renewal of interests in nonlinear electrodynamics that is seen recently is motivated by the observation that such theories arose as limiting cases of certain formulations of string theory. Unfortunately, the no go theorem proved in Ref.  clearly shows that for Lagrangians $`(F),`$ $`(F=F_{ab}F^{ab}),`$ with the Maxwell weak-field limit there are no spherically-symmetric static black hole solutions with a regular center.

Recently, employing the Schwinger-DeWitt approximation we have constructed the renormalized stress-energy tensor of the quantized conformally coupled massive scalar field in the spacetime of the electrically charged regular black hole being an exact solution of the equations of the nonlinear electrodynamics and the Einstein field equations. Such a solution has been proposed by Ayón-Beato and García (ABG) in Ref. . It should be noted, however, the ABG line element is not a solution of the standard nonlinear electrodynamics and the effective geometry (i.e. the geometry seen by photons of the nonlinear theory) is singular. Fortunately, the ABG solution may be reinterpreted as describing a magnetically charged regular solution of the coupled equations of standard nonlinear electrodynamics and gravitation with much more regular behaviour of the effective geometry . Moreover, it has been shown recently that it is possible to combine electric and magnetic line elements to obtain a regular electric solution with a magnetic core .

For small and intermediate charges as well as at large distances, the geometry of the ABG black hole resembles that of RN; noticeable differences appear near the extremality limits. It is interesting therefore to analyse how the similarities in the metric structure of the ABG and RN black holes are reflected in the structure of the stress-energy tensors. As for the conformally coupled massive scalar field the linear, quadratic, and cubic terms in $`\eta `$ in the stress-energy tensor (9) are absent, it is anticipated that such similarities do occur. Explicit calculations carried out in  confirmed this hypothesis and showed that for small values of $`q`$ the appropriate tensors are practically indistinguishable, and, as expected, important differences appear near and at the extremality limit.

In this paper we shall extend the results of Ref.  and investigate a much more complicated case of arbitrary $`\eta `$ with special emphasis put on the minimal coupling. Although complexity of the approximate stress-energy tensor constructed for the ABG line element practically prevents its direct examination, it is possible to extract interesting information expanding $`T_a^b`$ into a power series and retaining a few leading terms. We shall show that such analyses could be carried out for small charges, large distances, and in the vicinity of the event horizon of the extremal ABG black hole. Moreover, on general grounds one may easily estimate the role played by $`\stackrel{~}{T}^{\left(1\right)ab}`$ and $`\stackrel{~}{T}^{\left(3\right)ab}.`$ To gain a deeper understanding of the problem, however, we employ numerical calculations.

The paper is organized as follows. In Sec. II the essentials of the ABG black hole geometry that are neccessary in further development are briefly described. In Sec. III certain features of the approximate stress-energy tensor of the massive scalar field in the ABG geometry are discussed and compared with the appropriate tensors constructed in the RN geometry. The results of our numerical analyses are presented in Sec. IV in which we discuss behaviour of the components of $`T_a^b`$ in some detail and present it graphically.

## II Geometry

The general Reissner-Nordström solution of the Einstein-Maxwell equations describing a static and spherically symmetric black hole of mass $`M,`$ electric charge $`e,`$ and magnetic monopole charge $`e_m`$ has a remarkably simple form

$$ds^2=A\left(r\right)dt^2+A^1\left(r\right)dr^2+r^2\left(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2\right),$$
(11)
with

$$A\left(r\right)=1\frac{2M}{r}+\frac{e^2+e_m^2}{r^2}.$$
(12)
Since the charges enter (12) as a sum of squares the metric structure remains unchanged under the changes of charges leaving $`Q^2=e^2+e_m^2`$ constant. For $`Q^2/M^2<1,`$ the equation $`A\left(r\right)=0`$ has two positive roots

$$r_\pm =M\pm \sqrt{M^2Q^2},$$
(13)
and the larger one, $`r_+,`$ determines the location of the event horizon, whereas the smaller, $`r_{}`$, gives the position of the inner horizon. In the limit $`Q^2=M^{2\text{ }}`$ horizons merge at $`r=M`$ and the black hole degenerates to the extremal one. At $`r=0`$ the Reissner-Nordström solution has an irremovable curvature singularity, which, although hidden to an external observer for $`Q^2/M^21,`$ remains a somewhat unwanted feature of the solution. The solution for $`Q^2>M^{2\text{ }}`$ is clearly unphysical.

At least for purely electrical solutions the nonlinear electrodynamics does not remedy the situation. Indeed, consider ralization of the nonlinear electrodynamics with the action functional of the form

$$S=\frac{1}{16\pi }d^4x\sqrt{g}\left[R\left(F\right)\right],$$
(14)
where $`R`$ is a curvature scalar, $`F=F_{ab}F^{ab},`$ and $`\left(F\right)`$ is an arbitrary function with Maxwell asymptotic in a weak field limit, i. e. $`\left(F\right)F`$ and $`\frac{d}{dF}\left(F\right)1`$ as $`F0`$. According to a well known theorem , there are no regular, static, and spherically symmetric solutions of general relativity coupled to nonlinear electrodynamics describing a black hole with nonzero electric charge. However, as has been explicitly demonstrated recently by Bronnikov, adopted hypotheses leave room for appropriate regular solutions with a nonzero magnetic charge. In this regard it is interesting to note that within a different formulation of the nonlinear electrodynamics proposed in  (a $`𝒫`$ framework according to the nomenclature of Refs. ) obtained from the standard one (the $``$ formulation) (14) by means of a Legendre transformation, Ayón-Beato and García constructed a regular black hole solution with a nonzero electric charge and mass . Their solution has the simple form (11), where

$$A\left(r\right)=1\frac{2M}{r}\left(1\mathrm{tanh}\frac{e^2}{2Mr}\right).$$
(15)
Bronnikov also demonstrated that any spherically symmetric solution constructed within the  $``$ framework has its counterpart with the same metric tensor constructed within the $`𝒫`$ framework, and therefore the electric solution (11) has the magnetic companion with $`e`$ replaced by the magnetic charge $`e_m.`$ Moreover, Burinskii and Hildebrandt have proposed recently a regular hybrid model in which electrically and magnetically charged solutions were combined in such a way that the electric field does not extend to the center of the black hole .

Although the magnetic and electric ABG solutions are precisely of the type (11) with (15), the geometries seen by photons of the nonlinear theory are different. It is because photons of the nonlinear theory move along null geodesics of the effective metric, and the latter is singular for the electric solution . It should be noted, however, that otherwise the physical geometry is regular and well behaving.

Since geometries outside the event horizon are described by the same line element and since we are going to consider neutral scalar field only, our results will hold for any particular realization of the ABG black hole as the only concern here is the metric structure of the spacetime. Therefore in what follows we denote both electric and magnetic charge by $`e.`$

Inspection of the metric potentials reveals interesting features of the ABG solution: its curvature invariants are finite as $`r`$ $`0`$ and at large $`r`$ it approaches the Reissner-Nordström solution. Moreover, for small and intermediate values of the ratio $`e^2/M^2,`$ the radial coordinate of the event horizon, $`r_+,`$ is close to the event horizon of the RN black hole; significant differences occur near the extremality limit.

It has been shown in  that the location of the horizons of the ABG black hole may be expressed in terms of the Lambert function $`W`$ . Indeed, making use of the substitution $`r=Mx`$, $`e^2=q^2M^2,`$ it could be demonstrated that the location of the horizons is given by the real branches of the Lambert functions:

$$x_\pm =\frac{4q^2}{4W(\epsilon ,q^2/4\mathrm{exp}(q^2/4))q^2},$$
(16)
where $`\epsilon `$ is $`0`$ for the event horizon, $`x_{+\text{ }},`$ and $`1`$ for the inner one, $`x_{}.`$ The functions $`W(0,s)`$ and $`W(1,s)`$ are the only real branches of the Lambert function with the branch point at $`s=\mathrm{e}^1,`$ where e is the base of natural logarithms and

$$W(0,\mathrm{e}^1)=W(1,\mathrm{e}^1)=1.$$
(17)
Consequently, for

$$q=\mathrm{\hspace{0.17em}2}\sqrt{w}$$
(18)
the horizons merge at

$$x_{extr}=\frac{4w}{1+w},$$
(19)
where $`w=W(0,\mathrm{e}^1)=0.2785.`$ Numerically, one has $`x_{extr}=\mathrm{\hspace{0.17em}0.871}`$ for $`\left|e\right|/M=\mathrm{\hspace{0.17em}1.055}.`$

For small values of $`q,`$ the ABG line element resembles that of RN. Indeed, expanding the function $`A\left(x\right)`$ one has

$$A\left(x\right)=1\frac{2}{x}+\frac{q^2}{x^2}+𝒪\left(q^4\right).$$
(20)
In spite of similarities between ABG and RN geometries for $`q1,`$ there are substantial differences in the extremality limit: the extremal RN line element is described by (11) with

$$A(x)=\left(1\frac{1}{x}\right)^2,$$
(21)
whereas near the event horizon of the extremal ABG black hole the function $`A(x)`$ has the following expansion

$$A(x)=(xx_{extr})^4+\frac{(1+w)^3}{32w^2}(xx_{extr})^5+O(xx_{extr})^6.$$
(22)
The form of the above expansion indicates that the proper length between the event horizon and any point located at $`r>r_{extr}`$ is infinite.

## III The renormalized stress-energy tensor

Structure of terms in the effective potential (8) and in Eq. (9) indicates that for the conformally coupled massive scalar field, i.e. for $`\eta =0,`$ there are substantial simplifications of the renormalized stress-energy tensor as in this very case the third, fourth, and fifth term in (8) do not contribute to the final result. Moreover, from (8) it is evident that in $`R=0`$ geometries and for arbitrary curvature coupling the functional derivative of the first and the third terms in (8) with respect to the metric tensor vanishes and, therefore, the approximate renormalized stress-energy tensor of the massive field in RN spacetime has a general form

$$T_a^b=C_a^b+\eta D_a^b.$$
(23)
On the other hand, however, for the ABG metric one has

$$R=\frac{q^4}{M^2x^5}\mathrm{tanh}\left(\frac{q^2}{2x}\right)\left[1\mathrm{tanh}^2\left(\frac{q^2}{2x}\right)\right],$$
(24)
which for small $`q`$ is $`𝒪\left(q^6\right).`$ Detailed calculations carried out for $`q1`$ show that for the ABG solution neither $`\stackrel{~}{T}_a^{\left(1\right)b}`$ nor $`\stackrel{~}{T}_a^{\left(3\right)b}`$ contribute importantly to the result.

Because of the similarities of the metric structure of the RN and ABG solutions, the overall behaviour of the renormalized stress-tensors for conformal coupling should be qualitatively similar and comparable in magnitude at least for small and intermediate values of $`q.`$ However, the differences between the line element (21) and the expansion (22), give strong evidence for differences in the appropriate components of the approximate stress-energy tensors. For the conformally coupled massive scalar field this statement has been confirmed by extensive numerical calculations reported in Ref. . Now we shall analyse the case of an arbitrary coupling with curvature. First, observe that contribution of $`\eta ^2\stackrel{~}{T}^{\left(1\right)ab}`$ and $`\eta ^3\stackrel{~}{T}^{\left(3\right)ab}`$ to the stress-energy tensor could be made arbitrarily great by a suitable choice of the conformal coupling. It should be noted, however, that such great values of $`\eta `$ are clearly unphysical and should be rejected. Therefore, as the particular case of the conformal coupling has been considered earlier, here we shall confine ourselves mostly to the important and physically interesting case of minimal coupling $`\eta =1/6.`$

In order to construct $`C_a^b`$ and $`D_a^b`$ one has to either calculate the curvature tensor and its covariant derivatives to required order for the line element (11) with (12) and employ the general formulas presented in  or to make use of the method proposed by Anderson, Hiscock and Samuel . Both approaches yield the same result, which reads:

$$C_t^t=\frac{810x^4q^2855x^4202x^2q^4+1878x^31152x^3q^22307x^2q^2+3084xq^41248q^6}{30240M^6\pi ^2m^2x^{12}},$$
(25)
$$D_t^t=\frac{1008x^3q^2+819q^6+2604x^2q^2+728x^2q^42712xq^4+360x^4792x^3}{720M^6\pi ^2m^2x^{12}},$$
(26)
$$C_t^t=\frac{842x^2q^4+444q^6+162x^4q^2462x^31488x^3q^21932xq^4+315x^4+2127x^2q^2}{30240M^6\pi ^2m^2x^{12}},$$
(27)
$$D_r^r=\frac{504xq^4208x^2q^4588x^2q^2+336x^3q^2117q^6+216x^3144x^4}{720M^6\pi ^2m^2x^{12}},$$
(28)
$$C_\theta ^\theta =\frac{2202x^3486x^4q^2945x^43044x^2q^4+4884x^3q^29909x^2q^2+10356xq^43066q^6}{30240M^6\pi ^2m^2x^{12}}$$
(29)
and

$$D_\theta ^\theta =\frac{1176x^3q^2+1053q^6+3276x^2q^2+832x^2q^43408xq^4+432x^41008x^3}{720M^6\pi ^2m^2x^{12}}.$$
(30)
The analogous tensor in the ABG geometry is more complicated, and besides the linear terms it contains also terms which are quadratic and cubic in $`\eta .`$ Each component of the stress-energy tensor consists of more than three hundred terms and has a general form

$$\frac{1}{96\pi ^2m^2M^6}\left(1\mathrm{tanh}\frac{q^2}{2x}\right)\underset{i,j,k}{}\alpha _{ijk}\frac{q^{2i}}{x^j}\mathrm{tanh}^k\frac{q^2}{2x},$$
(31)
where $`\alpha _{ijk}`$ for a given $`\eta `$ are numerical coefficients $`(0i6,`$ $`8j15`$, $`0k8),`$ and for obvious reasons will not be presented here.

Since the form of the Eq. (31) differs considerably from the stress-energy tensor given by (23) with (2530), it could be expected that the runs of both tensors have nothing in common. On the other hand similarities in the metric structure of the RN and ABG black holes discussed earlier strongly suggest the opposite. Unfortunately, the complexity of the renormalized stress-energy tensor of the massive scalar field in the ABG geometry practically invalidates direct analytical treatment. One can, however, obtain interesting and important information analysing certain limiting cases. Indeed, expanding the stress-energy tensor for $`q1`$ one has

$$T_a^b^{ABG}=T_a^b+\mathrm{\Delta }_a^{\left(1\right)b}+𝒪\left(q^8\right),$$
(32)
where $`T_a^b`$ is given by (2330) and

$`96\pi ^2m^2\mathrm{\Delta }_t^{(1)}^t`$ $`=`$ $`{\displaystyle \frac{q^6\left(1204911416x+2660x^2\right)}{420M^6x^{12}}}+{\displaystyle \frac{2q^6\left(617596x+140x^2\right)\eta }{5M^6x^{12}}}`$ (34)
$`{\displaystyle \frac{12q^6\left(459366x+70x^2\right)\eta ^2}{M^6x^{12}}},`$
$`96\pi ^2m^2\mathrm{\Delta }_r^{(1)}^r`$ $`=`$ $`{\displaystyle \frac{q^6\left(15551064x+140x^2\right)}{420M^6x^{12}}}{\displaystyle \frac{2q^6\left(121140x+40x^2\right)\eta }{5M^6x^{12}}}`$ (36)
$`+{\displaystyle \frac{12q^6\left(8184x+20x^2\right)\eta ^2}{M^6x^{12}}},`$
$`96\pi ^2m^2\mathrm{\Delta }_\theta ^{(1)}^\theta `$ $`=`$ $`{\displaystyle \frac{q^6\left(52493528x+560x^2\right)}{420M^6x^{12}}}+{\displaystyle \frac{2q^6\left(779720x+160x^2\right)\eta }{5M^6x^{12}}}`$ (38)
$`{\displaystyle \frac{12q^6\left(540423x+80x^2\right)\eta ^2}{M^6x^{12}}}.`$
Inspection of expansions (3438), which are valid for any $`\eta `$ and $`x,`$ shows that in this order the term $`\stackrel{~}{T}_{\left(3\right)}^{ab}`$ is absent in the final result. It is because $`\stackrel{~}{T}_a^{\left(3\right)b}𝒪\left(q^{12}\right).`$ The second term $`\stackrel{~}{T}_a^{\left(1\right)b},`$ which vanishes in the RN geometry is now $`𝒪\left(q^6\right).`$

A similar expansion valid for large values of radial coordinate has the following form

$$T_a^b^{ABG}=T_a^b+\mathrm{\Delta }_a^{\left(1\right)b}+\mathrm{\Delta }_a^{\left(2\right)b}+𝒪\left(x^{13}\right),$$
(39)
where $`T_a^b`$ is given by (2330), $`\mathrm{\Delta }_a^{\left(1\right)b}`$ by (3438), and

$$96\pi ^2m^2\mathrm{\Delta }_t^{(2)}{}_{}{}^{t}=\frac{q^8\left(2574+1015q^2\right)}{210M^6x^{12}}\frac{14q^8\left(39+20q^2\right)\eta }{5M^6x^{12}}+\frac{105q^8\left(27+8q^2\right)\eta ^2}{M^6x^{12}},$$
(40)
$$96\pi ^2m^2\mathrm{\Delta }_r^{(2)}{}_{}{}^{r}=\frac{q^8\left(738+35q^2\right)}{378M^6x^{12}}+\frac{4q^8\left(51+28q^2\right)\eta }{9M^6x^{12}}\frac{5q^8\left(333+112q^2\right)\eta ^2}{3M^6x^{12}},$$
(41)
$$96\pi ^2m^2\mathrm{\Delta }_\theta ^{(2)}{}_{}{}^{\theta }=\frac{q^8\left(12456875q^2\right)}{1890M^6x^{12}}\frac{112q^8\left(12+5q^2\right)\eta }{9M^6x^{12}}+\frac{5q^8\left(1989+560q^2\right)\eta ^2}{3M^6x^{12}}.$$
(42)
Note that in Eqs. (4042) there is no contribution from the term $`\stackrel{~}{T}_a^{\left(3\right)b}`$ and, consequently, the result which is valid for any combination of couplings and charges is independent of $`\eta ^3.`$

In the vicinity of the event horizon of the extremal black hole the renormalized stress-energy tensor has the following expansion

$$96\pi ^2m^2T_a^b^{ABG}=t_a^{\left(1\right)b}+t_a^{\left(2\right)b}\left(xx_{extr}\right)+𝒪\left(xx_{extr}\right)^2,$$
(43)
$`t_t^{\left(1\right)t}=t_r^{\left(1\right)r}=`$ $`{\displaystyle \frac{1}{4096}}{\displaystyle \frac{\left(1+w\right)^6\left(2+w\right)\left(w1\right)^2\eta ^3}{w^6M^6}}{\displaystyle \frac{1}{245760}}{\displaystyle \frac{\left(1+w\right)^6\left(4+w+w^3+2w^2\right)\eta }{w^6M^6}}`$ (45)
$`+{\displaystyle \frac{1}{1935360}}{\displaystyle \frac{\left(1+w\right)^6\left(w^3+3w^2+3w+5\right)}{w^6M^6}},`$
$`t_t^{\left(2\right)t}=t_r^{\left(2\right)r}=`$ $`{\displaystyle \frac{3}{16384}}{\displaystyle \frac{\left(1+w\right)^7\left(w+3\right)\left(w^22w+1\right)\eta ^3}{w^7M^6}}{\displaystyle \frac{1}{2580480}}{\displaystyle \frac{\left(1+w\right)^7\left(w+3\right)\left(w^2+3\right)}{w^7M^6}}`$ (47)
$`+{\displaystyle \frac{1}{983040}}{\displaystyle \frac{\left(1+w\right)^7\left(w+3\right)\left(3w^22w+7\right)\eta }{w^7M^6}},`$
$`t_\theta ^{\left(1\right)\theta }`$ $`=`$ $`t_\varphi ^{\left(1\right)\varphi }={\displaystyle \frac{1}{8192}}{\displaystyle \frac{\left(1+w\right)^6\left(w+5\right)\left(w1\right)^2\eta ^3}{w^6M^6}}+{\displaystyle \frac{1}{491520}}{\displaystyle \frac{\left(1+w\right)^6\left(w+13+3w^2+w^3\right)\eta }{w^6M^6}}`$ (49)
$`{\displaystyle \frac{1}{3870720}}{\displaystyle \frac{\left(1+w\right)^6\left(3w+17+3w^2+w^3\right)}{w^6M^6}}`$
and

$`t_\theta ^{\left(2\right)\theta }`$ $`=`$ $`t_\varphi ^{\left(2\right)\varphi }={\displaystyle \frac{3}{32768}}{\displaystyle \frac{\left(1+w\right)^7\left(25w+3423w^2+9w^3+5w^4\right)\eta ^3}{M^6w^7}}`$ (52)
$`{\displaystyle \frac{1}{16384}}{\displaystyle \frac{\left(1+w\right)^7\left(10w8+w^2+4w^3+w^4\right)\eta ^2}{M^6w^7}}`$
$`+{\displaystyle \frac{1}{1966080}}{\displaystyle \frac{\left(1+w\right)^7\left(185w190111w^2+5w^415w^3\right)\eta }{M^6w^7}}.`$
Since $`\stackrel{~}{T}_a^{\left(1\right)b}`$ vanishes in the limit $`xx_{extr}`$ there are no terms proportional to $`\eta ^2`$ in (4549).

The components of the stress-energy tensor of the massive scalar field are regular functions of the radial coordinate and are finite on the event horizon. Moreover, it could be demonstrated that the difference between radial and time components of $`T_a^b^{ABG}`$ factorizes:

$$T_t^t^{ABG}T_r^r^{ABG}=\left[1\frac{2M}{r}\left(1\mathrm{tanh}\frac{e^2}{2Mr}\right)\right]F\left(r\right),$$
(53)
where $`F\left(r\right)`$ is a regular function and, consequently, the stress-energy tensor is finite in a freely falling frame. To demonstrate this let us consider a slightly more general line element

$$ds^2=f\left(x^1\right)\left(dx^0\right)^2+g\left(x^1\right)\left(dx^1\right)^2+\left(x^1\right)^2d\mathrm{\Omega }^2\text{.}$$
(54)
For a radial motion the vectors of the frame are the four-velocity $`e_0^\alpha =u^\alpha `$ and a unit length spacelike vector $`e_1^\alpha =n^\alpha `$. Then, using the geodesic equation, one finds

$$e_{\left(0\right)}^a=u^a=(\frac{\gamma }{f},\sqrt{\left(\frac{\gamma ^2}{f}1\right)\frac{1}{g}},0,0)$$
(55)
and

$$e_{\left(1\right)}^a=n^a=(\frac{1}{f}\sqrt{\gamma ^2f},\frac{\gamma }{\sqrt{fg}},0,0),$$
(56)
where $`\gamma `$ is the energy per unit mass along the geodesic. A simple calculation shows that the components $`T_{(0)(0)}`$, $`T_{(0)(1)}`$ and $`T_{(1)(1)\text{ }}`$in a freely falling frame (independent of the function $`g\left(x^1\right)`$) are

$$T_{(0)(0)}=\frac{\gamma ^2(T_1^1T_0^0)}{f}T_1^1\text{,}$$
(57)
$$T_{(1)(1)}=\frac{\gamma ^2(T_1^1T_0^0)}{f}+T_1^1\text{,}$$
(58)
$$T_{(0)(1)}=\frac{\gamma \sqrt{\gamma ^2f}(T_1^1T_0^0)}{f}\text{}$$
(59)
One concludes, therefore, that if all components of $`T_a^b`$ and

$$\frac{(T_1^1T_0^0)}{f}$$
(60)
are finite on the horizon, the stress-energy tensor in a freely falling frame is finite as well. It is expected that constructed formulas satisfactorily approximate the exact stress-energy tensor. On the other hand however, to establish the regularity of the exact $`T_a^b^{ABG}`$ in a freely falling frame one has to explicitly demonstrate that (60) is satisfied.

## IV Numerical results

Considerations of the previous section are limited to analytically tractable special cases. However, to gain insight into the overall behavior of the stress-energy tensor for any combination of couplings and charges one has to refer to numerical calculations — our complete but rather complicated analytical formulas are, unfortunately, not of much help in this regard. Below we describe the main features of the constructed tensors and present them graphically, fixing our attention on the physically interesting case of minimal coupling. Related plots showing the radial dependence of $`C_a^b`$ and $`T_a^{(0)b}`$ can be found in Ref. , where the discussion of the stress-energy tensor of the massive scalar field with conformal coupling with curvature has been carried out.

The case of arbitrary $`\eta `$ is much more complicated as the tensor $`T_a^b^{ABG},`$ given by (10), is modified by the presence of additional terms. However, numerical analysis performed for $`\eta =1/6`$ reveals that the contribution of $`\eta ^3T_a^{(3)b}`$ is neglegible for $`q`$ even as large as $`1.0`$. Moreover, a closer examination indicates that as one approaches the extremality limit, the magnitude of this very term becomes comparable with $`\eta ^2T_a^{(2)b}`$ and $`\eta T_a^{(1)b}`$ only in the closest vicinity of the event horizon. Our calculations also suggest that irrespective of the exact value of $`q,`$ each component of $`T_a^{(1)b}`$ is a monotonic function of the radial coordinate in a large neighbourhood of the event horizon, if not in the whole range $`r>r_+.`$ For the extremal case such a behaviour is illustrated in Figs. 1–3, where the time, radial, and angular components of $`T_a^{(1)b}`$ are displayed. For comparison and completeness we also present the run of $`T_a^{(0)b},`$ $`C_a^b,`$ and $`D_a^b,`$ where the latter two tensors have been calculated for the RN line element with the aid of Eqs. (2530). It should be noted that even in the extremal case, when the differences between two considered geometries are most prominent, the tensors $`T_a^{(1)b}`$ and $`D_a^b`$ behave in a similar manner though differing noticeably in magnitude. For small values of $`q,`$ the curves constructed for both types of black hole are almost indistinguishable, which could be easily established from Eq. (32).

On the other hand however, the components of $`T_a^{(2)b}`$ exhibit quite different but very regular behavior, which will be described in some detail. For $`q1`$ it can be, of course, inferred from the approximate formulas (3238); greater values of $`q`$ require numerical examination. Specifically, for small charges each component of the considered tensor is negative at $`r=r_+,`$ and before approaching zero as $`r\mathrm{},`$ it changes sign once for the radial component and twice for time and angular components. With the increase of $`q,`$ all components become strongly negative at $`r_+,`$ with their zeros shifted towards larger values of the radial coordinate. Near the extremality limit, the tensor $`T_a^{(2)b}`$ is very sensitive to changes of $`q,`$ which is illustrated in Figs. 4–6. If the value of ratio $`\left|e\right|/M`$ slightly exceeds 1.0, for all components there occurs a negative minimum not far away from the event horizon. For the extremal ABG black hole $`T_a^{(2)b}`$ vanishes at $`r_+,`$ as is clearly seen from Eqs. (4349), and exhibits oscillatory-like behavior with a rapidly decreased amplitude and increased intervals between zeros. The angular component attains a very distinct maximum near the event horizon whereas for time and radial component there are minima close to $`r_+`$.

The competition of the terms described above results in the overall behaviour of the stress-energy tensor of minimally coupled massive fields in the geometry of the ABG black hole. Numerical calculations indicate that for small and intermediate values of $`q,`$ up to about $`0.8,`$ $`T_a^b^{ABG}`$ still resembles that evaluated for conformal coupling. For larger values of $`q,`$ however, the change of curvature coupling leads to a considerable modification of the results as both magnitudes and radial variations become significantly different. This is illustrated in Figs. 7–9, where values of $`q`$ at and near the extremality limit were chosen. The oscillatory-like behaviour of $`T_a^{(2)b}`$ for $`q`$ close to $`q_{extr}`$ is refected by the presence of local extrema visible in the stress-energy tensor, and, for the angular component, an inflection point. It is in a sharp contrast with the almost monotonic behavior of the stress-energy tensor of conformally coupled massive fields depicted in Fig. 7–9 by the dashed lines. Moreover, it should be noted that for $`\eta =1/6,`$ when the extremality limit is approached, there are substantial changes of values of $`T_a^b^{ABG}`$ at the event horizon as well as in the behavior of the stress-energy tensor in a narrow strip near $`r_{extr}`$. There appear also certain new features for the time and radial components which are not present in the case of the conformal coupling. Indeed, for $`0.987q1.032,`$ the energy density $`\rho =T_t^t^{ABG}`$ is positive at the event horizon whereas for the same values of $`q,`$ the horizon value of the radial pressure $`p_r=T_r^r^{ABG}`$ is negative. It should be noted in this regard that for $`\eta =0`$ the radial component of the stress-energy tensor is always positive there. On the other hand, the angular pressure is positive on the event horizon for $`\left|e\right|/M0.922`$ and negative for larger values that is very similar to the previously studied case of conformal coupling.

## V Concluding remarks

In this work our goal was to construct the renormalized stress-energy tensor of the quantized massive field in the spacetime of a nonlinear black hole and to investigate how the choice of the curvature coupling affects the results. A regular electrically charged solution of this type has been recently proposed by Ayón-Beato and García in $`𝒫`$ formulation of nonlinear electrodynamics and reinterpreted by Bronnikov as a regular magnetically charged solution of the standard $``$ formulation. For small and intermediate values of the ratio $`\left|e\right|/M`$ the metric structure of the nonlinear solution closely resembles that of Reissner and Nordström and the similarities in the line elements are reflected in the behavior of the stress-energy tensor of the conformally coupled massive scalar fields; notable differences appear near and at the extremality limit.

As the general $`T_a^b^{ABG}`$ contains terms that are quadratic and cubic in $`\eta `$, the case of arbitrary coupling is more complicated. Again, for small and intermediate values of $`q`$ there is a similarity between the $`T_a^b`$ evaluated for the minimally coupled scalar field in the ABG geometry and its RN counterpart. Modifications for $`q`$ slightly exceeding 0.8, which are mainly due to the term $`T_a^{(2)b},`$ are noticeable although there are still some similarities. Comparison of $`T_a^b^{ABG}`$ for charges close to $`q_{extr}`$ indicates that the radial dependences for both couplings are completely different. Indeed, the behaviour of the stress-energy tensor for the minimal coupling is far more complicated as compared with its almost monotonic radial dependence for $`\eta =0`$. Moreover, the former is even more sensitive to changes of q, especially near the extremality limit. In addition, new features occur as, for example, positivity of the energy density at the event horizon for certain values of $`q.`$

It should be emphasized, that being local the Schwinger-DeWitt approximation does not describe particle creation, and, therefore, it must not be applied in strong or rapidly varying gravitational fields. However, it is expected that for sufficiently massive fields the method provides a good approximation of the exact renormalized stress-energy tensor.

Finally, we remark that although complicated, the derived stress-energy tensor may be employed as a source term of the semi-classical Einstein field equations. Indeed, preliminary calculations indicate that it is possible to find an analytical perturbative solution to the back reaction, and what makes this issue even more interesting and worth further studies is the regularity of the geometry of the ABG black hole. We intend to return to this group of problems elsewhere.