text:
\section{Introduction}
Condensates of resonantly pumped exciton-polaritons in semiconductor
microcavities constitute a novel and exciting system for the study of
fundamental physical properties of superfluids out of
equilibrium~\cite{keeling09}, and for future device
applications~\cite{menon10,amo10}.  Being quantum superpositions of
light and matter they are privileged candidates for the realization of
the next generation of optical devices~\cite{menon10}, for example for
quantum information technologies.

From the point of view of the fundamental physical properties,
particularly interesting is the case of a coherently pumped polariton
superfluid in presence of defects. Here, in contrast to the
corresponding equilibrium case, a weak residual drag force is always
present even at extremely high polariton
densities~\cite{cancellieri10}. Nevertheless, properties
paradigmatical of an equilibrium superfluids, such as frictionless
flow of polariton bullets~\cite{amo09}, quantized vortices and
metastable persistent flow~\cite{sanvitto09}, and the appearance and
disappearance of {\u C}erenkov-like waves~\cite{amo09_b} have been
recently observed in coherently driven exciton-polaritons.

In view of the potential device applications, especially important is
the unique versatility of the polaritonic system, which, combined with
its high non-linear properties, have been already demonstrated to
produce parametric
scattering~\cite{carusotto04,savvidis00:prl,stevenson00} and
bistability~\cite{baas04}. Here, the implementation of logic
operations and gates comes in a natural way: By manipulating the
non-linear properties of the system using several lasers which
frequencies, angles of incidence, and intensities can be freely varied
externally.

In this paper we investigate a new realization of a two-component
polariton system coherently driven by two lasers with independently
tunable frequencies, angles of incidence and intensities.  Firstly, we
study the stability of the two polariton components when the two laser
intensities are varied.
Differently from the case of a single laser pump, where the system can
only be bistable, we disclose a rich phase diagram, where either one,
two or three stable states can coexist at given pumping conditions.
Then we suggest possible easily tunable multistable hysteresis cycles
when the two pumping lasers are varied up and down in intensities.

Alternatively, a multicomponent polariton fluid can be realized by
considering the polarization degrees of freedom. Multistability of
different polariton spin states has been recently proposed
theoretically~\cite{gippius07} and confirmed
experimentally~\cite{paraiso10} by the observation of three stable
spin states for a given excitation condition.  For the case of two
polarized components, multistability in space have also been
theoretically proposed~\cite{shelykh08,shelykh06} and experimentally
observed~\cite{adrados10}.  For the system with two pumping lasers
presented here, the same kind of spatial multistability is expected
but with much more complex features.  Additionally, superimposed to
the spatial multistability, interference fringes will appear due to
the difference in frequency and momentum of the two pumping lasers.
The analogy with the two-component polarized case, suggests that,
aside from the interest in investigating multistability, two-component
polariton condensates obtained with independent lasers can also be
used to realize switches~\cite{amo10} and memories.

The paper is organized as follows: in section II, we present the model
used to describe the steady state behavior of polaritons excited by
two continuous-wave lasers with different frequencies, wave-vectors
and intensities. The results obtained within a linear response
framework are shown in section III. In this section we study the
number of the possible solutions, their nature and possible Kerr or
parametric instabilities associated to them. As a consequence of the
presence of multiple stable solutions, different cycles of hysteresis
can be produced by varying, along different paths, the intensities of
the two pumping lasers. Finally, section IV contains the conclusions
drawn from our analysis.

\section{Model}
The dynamics of resonantly-driven microcavity
polaritons~\cite{ciuti03,carusotto04} can be described via a
Gross-Pitaevskii equation for coupled cavity ($\psi_C$) and exciton
($\psi_X$) fields generalized to include decay and resonant pumping
($\hbar=1$):
\begin{equation}
  i\partial_t \begin{pmatrix} \psi_X \\ \psi_C \end{pmatrix} =
  \begin{pmatrix} 0 \\ F \end{pmatrix} +
  \left[\hat{H}_0 + \begin{pmatrix} g_X|\psi_X|^2& 0 \\ 0 &
      0 \end{pmatrix}\right]
  \begin{pmatrix} \psi_X \\ \psi_C \end{pmatrix}\; .
\label{eq:model}
\end{equation}
The repulsive ($g_X>0$) exciton-exciton interaction induces a
non-linear dynamics of the eigenmodes of the single polariton
Hamiltonian (lower and upper polariton, $\omega_{LP,UP}(\vect{k})$):
\begin{equation}
  \hat{H}_0 = \begin{pmatrix} \omega_{X} (-i\nabla) - i \kappa_X &
    \Omega_R/2 \\ \Omega_R/2 & \omega_{C}(-i\nabla) - i
    \kappa_C \end{pmatrix} \; .
\end{equation}
Here, we assume the cavity dispersion to be quadratic,
$\omega_C(\vect{k})=\omega_C(0)+ k^2/(2m_C)$, with
$m_C=2\times10^{-5}m_0$ ($m_0$ is the bare electron mass), we will
neglect the exciton dispersion and consider the case of zero detuning
at normal incidence, $\omega_X(\vect{k}) = \omega_X
(0)=\omega_C(0)$. The Rabi frequency $\Omega_R=5.0$~$[$meV$]$ and the
excitonic and photonic decay rates, $\kappa_X=\kappa_C=0.05$~$[$meV$]$
are chosen in the range of experimental values.

Because of the continuous decay, a stationary state requires a
continuous injection of photons. Here, we consider two continuous-wave
laser fields,
\begin{equation}
  F(\vect{r},t)=F_1 e^{i (\vect{k}_1 \cdot \vect{r} - \omega_1 t)}
  +F_2 e^{i (\vect{k}_2 \cdot \vect{r} - \omega_2 t)}\; ,
\end{equation}
with independently tunable frequencies $\omega_{1,2}$ and momenta
$\vect{k}_{1,2}$, which can be experimentally changed by changing the
laser angle of incidence with respect to the growth direction.

We study the mean-field solutions of Eq.~\eqref{eq:model}
\begin{equation}
  \psi_{X,C}(\vect{r},t) =
  \psi^{ss}_{1_{X,C}}e^{i(\vect{k}_1\cdot\vect{r}-\omega_1 t)} +
  \psi^{ss}_{2_{X,C}}e^{i(\vect{k}_2\cdot\vect{r}-\omega_2 t)} \; ,
\label{eq:meanf}
\end{equation}
and their stability with respect to small fluctuations within a linear
response analysis. Substituting the expression~\eqref{eq:meanf}
into~\eqref{eq:model} we obtain 4 contributions, two of which
oscillate at the main frequencies $\omega_1$ and $\omega_2$ and the
additional two at the replica (or satellite state) frequencies
$\omega_1-\Delta\omega$ and $\omega_2+\Delta\omega$, where $\Delta
\omega=\omega_2 - \omega_1$. Similarly to what is done in the OPO
regime~\cite{whittaker05,wouters07} where replica states in addition
to the pump signal and idler states are neglected, here, we consider
only the terms oscillating at the main frequencies $\omega_1$ and
$\omega_2$. Later, see Eq.~\eqref{eq:fluct}, we analyse the dynamical
stability of the two-pump-frequency solution against the weak
population of satellite states $\omega_i \pm \omega$ via parametric
scattering processes. Through the paper, we will consider only
dynamically stable two-pump-frequency solutions.
In this approximation, we obtain the following mean-field equations
for $\psi^{ss}_{{1,2}_{X,C}}$:
\begin{equation}
  \begin{cases}
  [\omega_X-\omega_1-i\kappa_X+G_{12}]\psi_{1_X}^{ss}+
    \frac{\Omega_R}{2}\psi_{1_C}^{ss}=0 \\
  [\omega_C({\bf k}_1)-\omega_1-i\kappa_C]\psi_{1_C}^{ss}+
  \frac{\Omega_R}{2}\psi_{1_X}^{ss}+F_1 =0 \\
  [\omega_X-\omega_2-i\kappa_X+G_{21}]\psi_{2_X}^{ss}+
  \frac{\Omega_R}{2}\psi_{2_C}^{ss}=0 \\
  [\omega_C({\bf k}_2)-\omega_2-i\kappa_C]\psi_{2_C}^{ss}+
  \frac{\Omega_R}{2}\psi_{2_X}^{ss} +F_2 =0 \; ,
\end{cases}
\label{eq:stead}
\end{equation}
where $G_{ij}=g_X(|\psi_{i_X}^{ss}|^2+2|\psi_{j_X}^{ss}|^2)$ with
$i\neq j=1, 2$.
Note that the repulsive interaction term between excitons in different
states is two times larger the interaction term between excitons in
the same mode, resulting in a non-uniform blue-shift.  The mean-field
system of equations~\eqref{eq:stead} can have up to 9 solutions,
i.e. 6 solutions more than in the case of one pumping laser, but, as
discussed below, only a maximum of 3 solutions are stable.

The dynamical stability of the two-pump-frequency mean-field solution
can be established by adding small fluctuations,
\begin{multline}
  \psi_{X,C}(\vect{r},t) = e^{-i\omega_1 t}
  \left[e^{i\vect{k}_1\cdot\vect{r}} \psi^{ss}_{1_{X,C}} +
    \theta_{1_{X,C}}(\vect{r,t})\right] +\\
  e^{-i\omega_2 t}
  \left[e^{i\vect{k}_2\cdot\vect{r}} \psi^{ss}_{2_{X,C}} +
    \theta_{2_{X,C}}(\vect{r,t})\right] \; ,
\label{eq:fluct}
\end{multline}
where the fluctuation fields can be divided into particle-like and
hole-like excitations $\theta_{i_{X,C}} (\vect{r},t) = \sum_{\vect{k}}
[e^{-i\omega t + i\vect{k}\cdot \vect{r}} u_{i_{X,C}\vect{k}} +
  e^{i\omega t + i(2 \vect{k}_i -\vect{k})\cdot \vect{r}}
  v^*_{i_{X,C}\vect{k}}]$. Expanding Eq.~\eqref{eq:model} up to linear
terms in $\theta_{1,2_{X,C}}$, we obtain 4 terms oscillating at
frequencies $\omega_1-\Delta\omega \pm \omega$ and
$\omega_2+\Delta\omega \pm \omega$, which we neglect, and 4 terms
oscillating at $\omega_i \pm \omega$. In other words, we are checking
the stability of our solution, where only the two states with
frequencies $\omega_{1,2}$ are occupied, against the weak population
of the satellite states $\omega_i \pm \omega$ which can be populated
by parametric scattering processes. The fact that we consider only
linear terms in $u_{i_{X,C}\vect{k}}$ and $v^*_{i_{X,C}\vect{k}}$
implies that we can obtain only the threshold conditions for such
parametric processes, as well as the nature of the instability,
whether of Kerr-type or parametric-type --- see later. The equations
for $u_{i_{X,C}\vect{k}}$ and $v^*_{i_{X,C}\vect{k}}$ can be written
as an eigenvalue equation rearranging the excitations into an
8-component vector $\mathbb{U}^{\text{T}} =
(u_{1_X},u_{1_C},v_{1_X},v_{1_C},u_{2_X},u_{2_C},v_{2_X},v_{2_C})$:
\begin{equation}
  \left[\omega \mathbb{I} - \begin{pmatrix} \mathbb{L}_{11 \vect{k}}&
      \mathbb{L}_{12 \vect{k}} \\ \mathbb{L}_{21 \vect{k}} &
      \mathbb{L}_{22 \vect{k}} \end{pmatrix}\right]
  \mathbb{U}_{\vect{k}} = 0\; .
\label{eq:eigen}
\end{equation}
Here matrices $\mathbb{L}_{ij \vect{k}}$ with $i\ne j$ are given by
\begin{equation*}
  2g_X e^{i(\vect{k}_i - \vect{k}_j)\cdot \vect{r}} \begin{pmatrix}
    \psi_{i_X}^{ss}\psi_{j_X}^{ss\star} & 0 & \psi_{i_X}^{ss}
    \psi_{j_X}^{ss} & 0 \\ 0 & 0 & 0 & 0 \\ -
    \psi_{i_X}^{ss\star}\psi_{j_X}^{ss\star} & 0 & -
    \psi_{i_X}^{ss\star} \psi_{j_X}^{ss} & 0 \\ 0 & 0 & 0 & 0
  \end{pmatrix}
\end{equation*}
and $\mathbb{L}_{jj \vect{k}}$ are given by
\begin{widetext}
\begin{equation*}
  \begin{pmatrix}
  \omega_X - \omega_j - i\kappa_X + g_X|\psi_X^{ss}|^2 &
  \frac{\Omega_R}{2} & g_X\psi_{j_X}^{ss}\psi_{j_X}^{ss} & 0
  \\ \frac{\Omega_R}{2} & \omega_C(\vect{k}) - \omega_j - i\kappa_C &
  0 & 0 \\ -g_X \psi_{j_X}^{ss\star} \psi_{j_X}^{ss\star} & 0 &
  -\omega_X (2 \vect{k}_j- \vect{k}) + \omega_j - i\kappa_X -
  g_X|\psi_X^{ss}|^2 & -\frac{\Omega_R}{2} \\ 0 & 0 &
  -\frac{\Omega_R}{2} & - \omega_C (2 \vect{k}_j- \vect{k}) + \omega_j
  - i\kappa_C \end{pmatrix}\; ,
\end{equation*}
\end{widetext}
with $|\psi_X^{ss}|^2=2(|\psi_{1_X}^{ss}|^2+|\psi_{2_X}^{ss}|^2)$
being the total excitonic density. At given values of the pumping
strength $F_1$ and $F_2$, the solutions of the mean-field
equations~\eqref{eq:stead} are stable if all the eight eigenvalues
($\text{LP}_j^{\pm}(\vect{k})$, $\text{UP}_j^{\pm}(\vect{k})$) of
Eq.~\eqref{eq:eigen} have negative imaginary part for every value of
the momentum $\vect{k}$.

\section{Results}
For some choices of the system parameters, we find that the number of
stable solutions can be larger than one. In the case of one pumping
laser, the typical $S$-like shape dependence of the polariton field
intensity on the pump strength, also referred to as optical
bistability, can be explained in terms of the non-linear blue-shift
induced by the polariton-polariton
interaction~\cite{whittaker05,ciuti05}.  When the laser frequency is
well above the bare lower polariton dispersion, $\omega_p >
\omega_{LP} (\vect{k}_p)$, and the pump intensity increases from low
values, the polariton population remains small because it is hard for
the laser to inject polartions with a different energy. However,
increasing the pump power, the blue-shift pulls the polariton energy
towards resonance with the pump causing the population to grow
superlinearly and eventually to abruptly jump to a high value when the
pump intensity reaches a critical value $I_1$. In the opposite
situation, when the laser intensity is decreased from high values, the
polariton energy is blue-detuned close to the pumping laser frequency
and, therefore, the cavity is efficiently filled by the laser even at
low pumping intensities. In this case the polariton population jumps
down back at low polariton densities for a value $I_2$ of the pump
strength lower than $I_1$. The two jumps at different values of the
pump intensity cause therefore a hysteresis cycle. As explained below,
in the case of two-component fluids, the situation becomes even
richer.
\begin{figure}
\begin{center}
\includegraphics[width=1.0\linewidth,angle=0]{2D.pdf}
\includegraphics[width=1.0\linewidth,angle=0]{3D.pdf}
\end{center}
\caption{(Color online) 2D panels: phase diagram showing the number of
  stable solutions as a function of the rescaled pump intensities
  $F_{1,2}'=\sqrt{g_x} F_{1,2}$ [meV$^{3/2}$]. White, green, yellow
  and black regions correspond to respectively zero, one, two or three
  stable solutions. In panels I and II $\vect{k}_1=0.25$~$\mu$m$^{-1}$
  and $\vect{k}_2=0.7$~$\mu$m$^{-1}$ while in panels III and IV
  $\vect{k}_1=0.0$~$\mu$m$^{-1}$ and
  $\vect{k}_2=0.7$~$\mu$m$^{-1}$. In the left panels (I and III)
  $\omega_{1,2}=\omega_{LP} (\vect{k}_{1,2})+0.3$~$[$meV$]$, while in
  the right panels (II and IV) $\omega_{1,2}=\omega_{LP}
  (\vect{k}_{1,2})+0.4$~$[$meV$]$. The horizontal black lines lies at
  the three fixed values of $F_2'$ corresponding to the three panels
  of Fig.~\ref{fig:cuts}, while the blue diagonal line is the path
  used to plot Fig.~\ref{fig:hyst}. 3D panels: plots of
  $g_X|\psi_{X}^{ss}|^2$ [meV] as a function of $F_{1,2}'$ with
  parameters equal to panel I. Stable solutions with higher
  populations are shown in green, stable solution with the second
  higher population in yellow and third stable solution with lower
  population in black. All the solutions are shown in the left
  panel. Since the upper green branches hide a yellow upper branch,
  the right panel shows only the yellow and black solutions.}
\label{fig:2d3d}
\end{figure}

We fix both laser frequencies to be blue detuned with respect to the
bare polariton dispersion: $\omega_{1,2}=\omega_{LP}
(\vect{k}_{1,2})+0.3$~$[$meV$]$, with $\vect{k}_1=0.25$~$\mu$m$^{-1}$
and $\vect{k}_2=0.7$~$\mu$m$^{-1}$. We plot in panel I of
Fig.~\ref{fig:2d3d} the phase diagram showing the regions with a
different number of stable solutions (either one, two or three) as a
function of the two rescaled pumping intensities $F_{1,2}'=\sqrt{g_x}
F_{1,2}$ [meV$^{3/2}$]. In order to understand better the structure of
this phase diagram, we show in Fig.~\ref{fig:cuts} the total exciton
density $g_X |\psi_X^{ss}|^2$ when the pump intensities $F_2'$ is kept
constant at different values and $F_1'$ is varied. When the constant
pump $F'_2$ has a small value (Fig.~\ref{fig:cuts} top left panel),
the dependence of the population on the varying pump intensity $F_1'$
is similar to the one-fluid case showing bistability with an $S$-like
shape. For higher values of the constant pump $F_2'$ (Fig.~\ref{fig:cuts}
lower left panel), the number of possible solutions increases but just
two are found to be stable. Finally when $F_2'$ is further increased
the set of possible solutions further goes up but only a maximum of
$3$ are found to be stable.
\begin{figure}
\begin{center}
\includegraphics[width=1.0\linewidth,angle=0]{scurves1D.pdf}
\end{center}
\caption{(Color online) Stability curves of the total exciton density
  $g_X |\psi_X^{ss}|^2$ $[$meV$]$ (red dotted curves unstable
  solutions, black lines stable solutions) for fixed pump intensities
  as a function of $F_1'$, for $F_2'=0.00001$ $[$meV$^{3/2]}$ (top
  left), $F_2'=0.025$ $[meV^{3/2}]$ (bottom left), and $F_2'=0.08$
  $[$meV$^{3/2}]$ (right).  Points $A_i$ ($B_i$), $i=1,2,3$ correspond
  to the cases when the lower (upper) branch of the stability curve
  becomes unstable (see Fig.~\ref{fig:imaginary1} and
  ~\ref{fig:imaginary2}).  Point $C$ does not have a counterpart in
  the one-fluid case and correspond to the cases when the second high
  branch of the stability curve becomes unstable (see
  Fig.~\ref{fig:imaginary2}).}
\label{fig:cuts}
\end{figure}

The coexistence of three solutions, corresponding to the black regions
of Fig.~\ref{fig:2d3d}, can be understood as follows: when the two
pump intensities increase from low values, the polariton population is
small because it's energy is far below the laser frequencies,
resulting in one stable solution. In the opposite situation, when the
intensity of one of the lasers decreases from high values, the
polariton population is high and its dispersion is significantly blue
detuned with respect to the bare one. Such blue-shift can be sustained
by any of the two lasers, thus giving two additional stable solutions
for the same values of the pump intensities. Therefore a maximum of
three stable solutions can be expected. This is also evident while
considering the partial densities for particular values of the pump
strength at which three stable solutions are present (black region of
Fig.~\ref{fig:2d3d}), e.g. $F_1'\equiv \sqrt{g_X} F_1= 0.05$
$[$meV$^{3/2}]$ and $F_2'\equiv \sqrt{g_X} F_1= 0.08$
$[$meV$^{3/2}]$. Here, the solution with lower total polariton density
corresponds to partial densities $g_X|\psi_{1_X}^{ss}|^2=
0.009$~$[$meV$]$ and $g_X|\psi_{2_X}^{ss}|^2=0.023$~$[$meV$]$. The
other two solutions correspond to a high value of just one of the two
partial populations: $g_X|\psi_{1_X}^{ss}|^2= 0.008$~$[$meV$]$ and
$g_X|\psi_{2_X}^{ss}|^2=0.609$~$[$meV$]$ in one case and
$g_X|\psi_{1_X}^{ss}|^2= 0.646$~$[$meV$]$ and
$g_X|\psi_{2_X}^{ss}|^2=0.010$~$[$meV$]$ in the other. Note that this
situation is similar to the case of two-component condensates obtained
with two spins. However, while in the spin-dependent case the two
lasers pump the two spin populations with different intensities but at
the same angle and energy, here the two pumps are independent also in
angle and in energy.
This analogy is also visible in Eq.~\eqref{eq:stead}. However, the
difference between our system of equations and the spin-dependent case
is that here the interaction between different components is twice the
interaction between particles in the same component.

In panel II of Fig.~\ref{fig:2d3d} we plot the phase diagram for the
same parameters as in panel I but with the two pumping lasers $0.4$
$[$meV$]$ blue detuned with respect to the bare LP branch. We see that
the effect of the increased detuning is simply to stretch the phase
diagram. Since the two pumps are further apart from the LP branch, it
is more difficult to inject polaritons into the cavity, and thus the
need for higher pump intensities.
In panel IV of Fig.~\ref{fig:2d3d} we show that a similar phase
diagram can be obtained by changing $\vect{k}$ vector of pump 1 from
$0.25$ to $0.0$~$\mu$m$^{-1}$. We observe that the multistability is
quite robust with respect to the choice of the parameters and,
therefore, it should be within an experimental reach. An interesting
configuration is plotted in panel III of the same figure. Here a
region with no stable solutions appears in the central part of the
plot (white region).  The instability of this region can be understood
by noting that with a pump at $\vect{k}=0.0$ and just slightly blue
detuned from the LP and a pump at $\vect{k}=0.7$, close to the
inflection point of the LP, it is easy to satisfy phase matching
conditions for parametric scattering processes. For this set of
parameters the system is in a configuration unstable towards the
population of satellite states by scattering processes. For the other
three sets of parameters, shown in the remaining three panels of
Fig.~\ref{fig:2d3d}, it is also possible to find regions of the phase
diagram where no solutions are stable. These are the regions where the
proposed solution, given by Eq.~\ref{eq:meanf}, where only the two
frequency states, $\omega_1$ and $\omega_2$, are populated, is not a
stable solution because satellite states start also to be populated
--- our analysis giving the threshold for this to happen.

\begin{figure}
\begin{center}
\includegraphics[width=1.0\linewidth,angle=0]{spectraIMfourINST.pdf}
\end{center}
\caption{(Color online) Dispersion of the imaginary part of the
  excitation eigenfrequency $\omega=\text{LP}_j^{\pm}$. The three
  panels correspond to points $A_i$ with $i=1,2,3$ of
  Fig.~\ref{fig:cuts}, where the lower part of the stability curves
  become unstable. In blue (red) the parts corresponding to the
  scattering of two particles with $\vect{k}=\vect{k}_1$
  ($\vect{k}=\vect{k}_2$).}
\label{fig:imaginary1}
\end{figure}

To further discuss the stability of the system with respect to small
perturbations, we plot the dispersion of the imaginary part of the
excitation eigenfrequency $\omega=\text{LP}_j^{\pm}$ for several
points of the stability curves shown in Fig.~\ref{fig:cuts}.  We start
with the cases where the lower branch of the stability curve became
unstable at points $A_i$. For very small values of $F_2'$ the
imaginary part of the dispersion (top panel of
Fig.~\ref{fig:imaginary1}) shows two peaks for given values of
$\vect{k}$. One peak lies at higher value $\vect{k}_+=0.83$
$\mu$m$^{-1}$ and one peak at lower value $\vect{k}_-=-0.33$
$\mu$m$^{-1}$ with $\vect{k}_++\vect{k}_-=2\vect{k}_1$. This two-peaks
structure is a precursor of a parametric instability due to the
scattering between two particles in the component of the condensate
with momentum $\vect{k}_1$. This situation corresponds exactly to the
case of one component fluids. When the pump intensity $F_2'$ is
slightly increased (middle panel) we observe 6 other peaks appearing
in the imaginary part of the dispersion. Two of these new peaks are
such that $\vect{k}_++\vect{k}_-=2\vect{k}_1$ (blue lines) while the
other four can be combined to identify two different scattering
processes with $\vect{k}_++\vect{k}_- =2\vect{k}_2$ (red curves). This
more complicated structure of the imaginary parts of the eigenvalues
is consistent with the fact that with two components a richer
mechanisms of scattering might occur. When the intensity of $F_2'$ is
further increased (lower panel), still 4 different scattering may
occur but, in this case, it is the scattering between two particles
with $\vect{k}_2$ that induce the instability of the system.

\begin{figure}
\begin{center}
\includegraphics[width=1.0\linewidth,angle=0]{spectraIMkerrINST.pdf}
\end{center}
\caption{Dispersion of the imaginary part of the excitation
  eigenfrequency $\omega=\text{LP}_j^{\pm}$. The three panels
  correspond to points $B_i$ with $i=1,2,3$ and $C$ of
  Fig.~\ref{fig:cuts}, where the higher part of the stability curve
  becomes unstable.}
\label{fig:imaginary2}
\end{figure}

For the transition from stable to unstable regions of the higher
branch of the stability curve we plot the dispersion of the imaginary
part of the excitation eigenfrequency $\omega=\text{LP}_j^{\pm}$ for
points $B_i$ in Fig.~\ref{fig:imaginary2}. In analogy with the case of
fluids with one component, for low intensities of pump 2, the
imaginary part of the dispersion shows a peak at the wavevector of
pump 1, a clear precursor of a Kerr instability. When the intensity of
pump 2 is increased new peaks appear in the imaginary part. In the
case of point $B_2$ two new peaks are precursors of a parametric
instability for the state with $\vect{k}=0.25$ $[\mu m^{-1}]$ even if
the mechanism responsible for the instability of the solution is still
of Kerr type. It is only when the pump intensity $F_2'$ is further
increased (point $B_3$) that the two peaks at $k_-=0.15$ and
$k_+=0.35$ $[\mu m^{-1}]$ became more important and the mechanism of
instability is of the parametric type.  Finally in the lower right
panel (corresponding to point C) of Fig.~\ref{fig:imaginary2} a single
peak at $\vect{k}=\vect{k}_1$ is the precursor of a Kerr-type
instability that ends the region with three stable solutions in the
right panel of Fig.~\ref{fig:cuts}.

\begin{figure}
\begin{center}
\includegraphics[width=1.0\linewidth,angle=0]{histtot.pdf}
\end{center}
\caption{\label{fig:hyst}(Color online) Hysteresis cycles of
  $\alpha_1=|\psi_{1_X}^{ss}|^2/(|\psi_{1_X}^{ss}|^2+|\psi_{2_X}^{ss}|^2)$
  [dimensionless] as a function of $F_1'$ for different values of
  $F_2'$ --- stable solutions are in black, while the hysteresis cycle
  performed by the system is in brown. Panel A: $F_2'=0.14 - F_1'$ as
  in the blue line of Fig.~\ref{fig:2d3d}. Panel B: $F_2'= 0.08$
  $[$meV$^{3/2}]$ as in the horizontal black line of
  Fig.~\ref{fig:2d3d}.}
\end{figure}
Multistability also manifests itself in a hysteresis loop for the
populations and emission intensities obtained with a cycle of first
increasing and later decreasing the pumping intensities. Here, the
presence of three stable solutions gives more complicated loops than
the ones obtained for a bistable system in a single-component
polariton fluid. In order to study this aspect, we calculate the
exciton emission intensity at a given frequency $\omega_i$ normalized
to the total exciton emission intensity,
$\alpha_i=|\psi_{i_X}^{ss}|^2/(|\psi_{1_X}^{ss}|^2+|\psi_{2_X}^{ss}|^2)$,
along a closed path of varying pumping intensities. The panels A) and
B) of Fig.~\ref{fig:hyst} respectively show the hysteresis cycles of
$\alpha_1$ when the two pump intensities change along either the blue
or the higher horizontal black line of Fig.~\ref{fig:2d3d}. In panel
A), one starts from a low value of $F_1'$ taking
$F_2'=0.13$~$[$meV$^{3/2}]$ so that the population of state $2$ is
much higher than the population of state $1$, i.e. $\alpha_1 \ll
1$. Increasing $F_1'$ the two populations smoothly evolve until $F_1'
\approx 0.05$ $[$meV$^{3/2}]$ and $F_2' \approx 0.09$
$[$meV$^{3/2}]$. At this point $F_2'$ is too weak to sustain high
population densities in state 2 and, therefore, the system jumps to a
new stable configuration, in which the populations of both states are
low, i.e. $\alpha_1 \approx 0.5$. A further increase of $F_1'$
produces a smooth evolution of the two populations until $F_1'\approx
0.11$ $[$meV$^{3/2}]$ when the system jumps to a third configuration
with a population in state 1 much higher than in state 2,
i.e. $\alpha_1\approx 1$. When we revert the variation of the pumping
intensities along the same path, the jumps to states corresponding to
intermediate and low values of $\alpha_1$ are shifted to the left of
the ones just described for increasing $F_1'$. The multistable
hysteresis loop shown in Fig.~\ref{fig:hyst}A) is related to the fact
that the two pumping lasers are at different pumping angles,
$\vect{k}_i$, and pumping frequencies, $\omega_i$. Therefore the jumps
from low to high population for each component appear at different
values of the pumping intensities, producing the multistable behavior
of $\alpha_1$.

A similar situation occurs when the system evolves along a path on
which one of the pumping intensities remains constant, while the other
varies (black horizontal line at $F_2'=0.08$ of panel 1 in Fig~\ref{fig:2d3d}),
as shown in Fig.~\ref{fig:hyst}B). Starting the with $F_1=0$,
$\alpha_1$ increases smoothly from zero following the lower branch
up to $F_1'=0.06$
 $[$meV$^{3/2}]$. At this point $\alpha_1$ jumps from
values of the order of $0.01$ to $0.3$, corresponding to a population
of state 1 being smaller but non-negligible compared to state 2. As $F_1'$
further increases up to $0.08$ $[$meV$^{3/2}]$, $\alpha_1$ again jumps
abruptly to values of the order of $0.95$. In the reverse process,
$F_1'$ is decreased completing the loop. Also in this case, the jumps
from high to low values of $\alpha_1$ are shifted to the left because
of the different angles and energies at which the lasers are pumping the
two component of the fluid. It is worth noting that the length and the
height of the different plateaus of the hysteresis loops can be efficiently
tuned by carefully choosing the path followed by the intensities, angles
and frequencies of the pumping lasers.

\section{Conclusion}
To summarize, we have studied the stability of a two-component
exciton-polariton fluid under resonant excitation of two pumping
lasers with independently tunable frequencies, angles of incidence and
intensities. We have studied the effect of the detuning between the
laser pump and the bare LP branch, and discussed the different kinds
of instability that might occur for different values of the
parameters. We have shown that, even though the kind of instabilities
are the same as in the one component case (Kerr or parametric), here
the interplay between different instabilities in the two components of
the fluid, can produce a much richer picture.  Moreover, we have shown
that since each component of the fluid jumps between stable states of
its population at a different values of the pumping intensity, the
system sustain multistable hysteresis loops that can be easily
modulated by changing the parameters of the pumping lasers. Finally,
we believe that, due to the wide range of parameters for which the
system is multistable, and due to the increased number of degrees of
freedom with respect to the spin-dependent case, multiatability and
hysteresis loops should be within an experimental reach. Therefore,
the novel system of two-component fluid is a promising candidate for
the realization of optical switches and memories.

We thank I. Carusotto and D. Whittaker for discussions. This research
has been supported by the Spanish MEC (MAT2008-01555,
QOIT-CSD2006-00019) and CAM (S-2009/ESP-1503). F.M.M. acknowledges
financial support from the program Ram\'on y Cajal.

\newcommand\textdot{\.}


text:
\section{Introduction}

\noindent
To understand the nature of  singularities in theories of gravity is a challenging  problem. It  can be  considered that  a  Quantum Theory of Gravity  may  solve this problem.
As for now, we are very
far from the  Quantum Theory of Gravity, we can avoid the singularities with  regular black hole solutions.
The  regular black hole solution
 first were given   by Bardeen \cite{Bardeen}. Later, this solution of Bardeen was obtained from the field equations of the
 Einstein-Nonlinear electrodynamics
\cite{Beato}.
It is interesting to obtain
some  regular black hole solutions
of
the  theory which is $f(R)$ minimally coupled to the Non-linear electrodynamics.
In recent years,
various new regular black holes  were proposed and investigated increasingly in the literature \cite{Hawyard}-\cite{Rodrigues}
 (see for a review \cite{Ansoldi}).
  Since the non-minimally coupled $Y(R)$-Maxwell  models
  \cite{dereli1,dereli2,Liu,bamba1,bamba2,Sert,dereli3,dereli4,Sert2}
  have some solutions which
  can explain the rotation curves of galaxies
  and cosmic acceleration    of the universe,
then; it is natural to seek regular black hole solutions of  the non-minimally
coupled electromagnetic fields to gravity.
 We focus on this subject in this paper.

According to the
Penrose-Hawking singularity theorem \cite{Hawking},
 to arise
   a singularity inside the horizon of a     black hole,
  the strong energy condition (SEC)
  has to be satisfied. The regular black holes violate   the strong energy condition  in the central region inside the black hole.
   We find various   models with the non-minimally coupled
   $Y(R)$ function
   for some known regular metric functions. We calculate the energy conditions for the effective energy-momentum tensor of these models. Then we find that
 they lead to a negative tangential pressure in the central core,
and
 the effective equation of state with negative radial pressure $p_r= -\rho$ is everywhere, which is important for
the accelerated expansion  phase of the Universe. We see that at least SEC is violated  by these solutions in some central regions of  the black holes.

\section{ The Gravitational Model with $Y(R)F^2 $-Type Coupling } \label{model}

\noindent We start with  the action with the  $Y(R)F^2 $-type non-minimal coupling term \cite{dereli3,dereli4,bamba2}
\begin{equation}
        I[e^a,{\omega^a}_b,F] = \int_M{\{   \frac{1}{2\kappa^2} R*1 -\frac{1}{2}Y(R) F\wedge *F  + \lambda_a\wedge T^a \} }.
        \nonumber
\end{equation}
Here  $\{e^a\}$ is the co-frame 1-form,
  ${\{\omega^a}_b\}$ is  the  connection 1-form,  $F=dA=  \frac{1}{2} F_{ab} e^a  \wedge  e^b $ is the homogeneous electromagnetic
  field 2-form,
$\lambda_a $ is the  Lagrange multiplier 2-form whose variation leads to the torsion-free Levi-Civita connection.
 Then the connection can be found from $T^a =  de^a + \omega^{a}_{\;\;b} \wedge e^b=0$.   In this action,  $R$  is the curvature scalar which can be obtained by this operation  \   $    \iota_{ba} R^{ab}= R $  from the curvature tensor 2-forms
$ R^{a}_{\;\;b} = d\omega^{a}_{\;\;b} + \omega^{a}_{\;\;c} \wedge \omega^{c}_{\;\;b} $ via   the interior product $\imath_a $,  and
 $\kappa^2 = 8\pi G$ is   universal gravitational coupling constant.   We take the  space-time metric $g = \eta_{ab} e^a \otimes
  e^b$  with the signature  $(-+++)$.
  We set  the orientated volume element as
  $*1 = e^0 \wedge e^1 \wedge  e^2 \wedge e^3 $.

  We obtain  gravitational and electromagnetic field equations of the theory    by taking  infinitesimal variations of  the action according to  independent variations of
    $\{e^a\}$,
   ${\{\omega^a}_b\}$ and $\{ A\}$
    \cite{dereli3,dereli4}
 \begin{eqnarray}\label{einstein}
  - \frac{1}{2 \kappa^2}  R^{bc}
 \wedge *e_{abc} =    && \frac{1}{2} Y  (\iota_a F \wedge *F - F \wedge \iota_a *F)   + \frac{1}{2}Y_R F_{mn} F^{mn}*R_a
 \nonumber
 \\
 && + \frac{1}{2}  D [ \iota^b D(Y_R F_{mn} F^{mn} )]\wedge *e_{ab}
 \  ,
 \end{eqnarray}
 \begin{eqnarray}\label{maxwell1}
 d(*Y F) = 0 \; , \hskip 3 cm \  \   \   \   dF=0 \;
 \end{eqnarray}
 where  $Y_R = \frac{dY}{dR}$. The gravitational field equation (\ref{einstein}) can be  written as
 \begin{eqnarray}\label{Ga}
\frac{G^a}{\kappa^2} =  \tau^a
 \end{eqnarray}
 where $G_a =- \frac{1}{2 }  R^{bc}
 \wedge *e_{abc} = *R_a -\frac{1}{2}R*e_a $ is the Einstein tensor, and $\tau_a = \tau_{a,b} *e^b $ is the effective energy momentum tensor for this non-minimally coupled model,  which is equal to right hand side of (\ref{einstein}). The effective energy density, radial  pressure,  and tangential  pressures are found from   $  \rho =  \tau_{0,0} $, \
  $  p_r =\tau_{1,1} $, \
   $ p_t = \tau_{ 2,2} = \tau_{3,3} $
using the field equation (\ref{Ga}).
 \section{Regular Black Hole Solutions}

 \noindent We seek regular black hole solutions for the following   (1+3)-dimensional  spherically symmetric  static line element
 \begin{equation}\label{metric}
 g = -f^2(r)dt^2  +  f^{-2}(r)dr^2 + r^2d\theta^2 +r^2\sin(\theta)^2 d \phi^2
 \end{equation}
 and consider   the electromagnetic tensor $F$ which may have electric and  magnetic components
 \begin{eqnarray}\label{electromagnetic1}
 F   &=&  E(r) dr\wedge dt +     B(r)r^2 \sin\theta d \theta \wedge  d \phi  = E(r) e^1\wedge e^0 +   B(r)e^2 \wedge e^3.
 \end{eqnarray}
 The field equations of the model for these ansatz turn out to  five  equations (three gravitational, two electromagnetic) and four   unknown functions ($E, B, Y, f$). These equations can be found in \cite{sertd}. The homogeneous electromagnetic field condition $dF= 0$ determines the magnetic field as
 \begin{eqnarray}\label{Bq}
 B=\frac{q}{r^2}
 \end{eqnarray}
where  $q $ is a real integration constant representing  magnetic monopole charge. It is impossible  to solve these equations without any simplification. Then, we use the following constraint to simplify these equations
 \begin{eqnarray}\label{cond}
 Y_R( E^2-  B^2) =  \frac{1}{ \kappa^2}\ .
 \end{eqnarray}
 Under this constraint the number of  equations decreases to two
 \begin{eqnarray}\label{Yf}
 {f^2}'' -\frac{2}{r^2}(f^2 -1 )  = \kappa^2 Y (E^2 +B^2),
 \end{eqnarray}
 \begin{eqnarray}\label{YE}
 YE = \frac{q_e}{ r^2},
 \end{eqnarray}
 where $q_e$ is the electric charge.  We note that one can find the constraint  (\ref{cond})      by taking
differential of the equation (\ref{Yf}).
To show this,
 we rewrite the equation (\ref{Yf})  using the magnetic field   $B=\frac{q}{r^2}$ from  (\ref{Bq}) and the electric field  $E=\frac{q_e}{Yr^2}$  from (\ref{YE}) and
   find

\begin{eqnarray}\label{YeniY}
r^4\left( {f^2}'' -\frac{2}{r^2}(f^2 -1 )   \right) = \kappa^2( \frac{q_e^2}{Y}  + q^2 Y)\ .
\end{eqnarray}
After taking differential  of equation (\ref{YeniY})  we obtain
\begin{eqnarray}\label{difY}
\left(- {f^2}''' - \frac{4}{r} {f^2}'' + \frac{2}{r^2} {f^2}' + \frac{4}{r^3} (f^2 -1)\right) dr= \kappa^2(\frac{q_e^2}{Y^2r^4} -\frac{q^2}{r^4}    )dY \ .
\end{eqnarray}
We see that left hand side of the equation (\ref{difY}) is equal to differential of the curvature scalar  $R= - {f^2}'' - \frac{4}{r} {f^2}' - \frac{2}{r^2} (f^2 -1)$. Using the  equations (\ref{Bq}), (\ref{YE}) and $dR$ in (\ref{difY})  we obtain the following equation
\begin{eqnarray}
dR = \kappa^2(E^2 - B^2)dY
\end{eqnarray}
which is  differential form of the constraint (\ref{cond}).
 Thus,  it is obvious that we have two equations  (\ref{Yf}), (\ref{YE}) and three unknowns $(f, Y, E)$. Then, for a given model with a non-minimal function Y(R), we can  determine the metric function $f(r)$ from these equations. On the other hand, for a desiring metric we can reach a corresponding  model with a  non-minimal Y(R) function. In order to be  successful in this process,  we need  to solve  $r$  from $R(r)$ and  express the function Y depending   on $R$. Then we can determine the corresponding model.

\subsection{Regular Black Hole Solution-1}

 The   field  equations of this non-minimal $Y(R)F^2$ model  (\ref{Yf}), (\ref{YE}) accept
 the following  regular black hole   solution
 \begin{eqnarray}\label{metric1}
 f^2(r) = 1 - \frac{2 m }{r} \left( 1 -  \frac{1}{(1+ a^3 r^3)^{1/3}    } \right)
 \end{eqnarray}
with the magnetic field ($E=0$) and the non-minimal function
\begin{eqnarray}
B(r) &=& \frac{q}{r^2}\label{B} \\
Y(r) &=& \frac{
	8 ma^6r^7 }{\kappa^2 q^2(1+ a^3r^3)^{7/3} } \label{Y1}
\end{eqnarray}
where  we have  defined a new constant  $a= \frac{2m}{q^2}$. The metric function (\ref{metric1}) can be found  in \cite{Balart}  as an electrically charged solution of Einstein-Non-linear Electrodynamics.
 In this notation,
 we calculate  the curvature scalar $R$ and  the invariant 4-form $R_{ab}\wedge *R^{ab}$ for the metric function (\ref{metric1})
\begin{eqnarray}
R(r)  &=& \frac{8m a^3}{(1+ a^3r^3)^{7/3}} \label{R} \\
R_{ab}\wedge*R^{ab}
& =&
(  \frac{24 m^2}{r^6} - \frac{48m^2}{r^6p}- \frac{32m^2a^6}{p^7} +
\frac{24m^2}{r^6 p^2} +
\frac{40m^2 a^6}{p^8} +
\frac{32m^2a^{12} r^6}{ p^{14}}
\nonumber \\
&& -\frac{16m^2a^3}{r^3p^4}
+\frac{16m^2a^3}{r^3p^5}  )*1
 \end{eqnarray}
where we have defined $p= (1+ a^3r^3)^{1/3}$.
When we check the limits
\begin{eqnarray}
\lim_{r\rightarrow 0 } R & =&  8m a^3 =\frac{64 m^4}{q^6}   \\
 \lim_{r\rightarrow 0 } R_{ab}\wedge*R^{ab}& =& \frac{16}{3}m^2a^6*1 = \frac{2^{10}m^8}{3q^{12}}*1,
\end{eqnarray}
 they are regular at the center of black hole. It needs to solve $r$ from (\ref{R}) to rewrite the non-minimal function $Y$  in terms of $R$
\begin{eqnarray}\label{r}
r(R) =\frac{1}{a}\left( (  \frac{8m a^3 }{R} )^{3/7} - 1   \right)^{1/3}.
\end{eqnarray}
When we substitute the inverse function (\ref{r}) in the non-minimal function  (\ref{Y1}) we obtain
\begin{eqnarray}
Y(R)= \frac{\left[ (8ma^3)^{3/7} -R^{3/7} \right]^{7/3 }    } {\kappa^2 q^2 a^4}
\end{eqnarray}
and the
corresponding model is written as:
\begin{eqnarray}
L =     \frac{1}{2\kappa^2} R*1 - \frac{\left[ (8ma^3)^{3/7} -R^{3/7} \right]^{7/3 }    } {2\kappa^2 q^2 a^4} F\wedge *F  +\lambda_a\wedge T^a.
\end{eqnarray}

 After the duality transformation
 $ B\rightarrow - YE $, $q \rightarrow -q_e $ and $Y\rightarrow
 \frac{1}{Y}$;
 which is  given in  \cite{sertd},
 we   reach the field equations  of the model for  the electromagnetic tensor $F$ with only electric component ($B=0,\ E \neq  0$). As a consequence of this transformation, the same  metric function (\ref{metric1}) determines the electric field and the non-minimal function of this model as follows
 \begin{eqnarray}\label{Y11}
 Y(r) &= & \frac{\kappa^2 q_e^2 (1+ a^3r^3)^{7/3}
 	}{8ma^6 r^7}  \\
  E(r) &=& \frac{4 q_e}{\kappa^2  r^2 } \left( 1+ \frac{1}{a^3 r^3 }\right)^{-7/3} = \frac{q_e}{Y(r)r^2}  \label{E}\;
 \end{eqnarray}
with   $a= \frac{2m}{q_e^2}$. We see that the electric field is regular at the center of black hole, $\lim_{r\rightarrow 0 } E(r) = 0$.
 We can rewrite the non-minimal function  (\ref{Y11}) in terms of $R$ as
  \begin{eqnarray}
 Y(R)= \frac{ \kappa^2 q_e^2 a^4  } {\left[ (8ma^3)^{3/7} -R^{3/7} \right]^{7/3 }     }
 \end{eqnarray}
 and we write the corresponding Lagrangian of the model
 \begin{eqnarray}
L =     \frac{1}{2\kappa^2} R*1 - \frac{\kappa^2 q_e^2 a^4 }{2 \left[ (8ma^3)^{3/7} -R^{3/7} \right]^{7/3 }    }  F\wedge *F   +\lambda_a\wedge T^a
 \end{eqnarray}
via the duality transformation.

 The same metric function (\ref{metric1})   and an electric field different from (\ref{E}) only up to a scale factor was obtained from   a different theory with Einstein-nonlinear electrodynamics  in \cite{Balart}.

 It is important to check all the energy conditions for this solution.    The conditions which are calculated as below  have to be  equal or greater than zero for rising  a singularity  in  General Relativity.  But in the regular black holes at least the strong energy condition has to be violated. To show this
  we firstly calculate the energy density $\rho(r) $, the radial and tangential pressures $p_r, p_t$ for the metric function ($\ref{metric1}$) and the magnetic field (\ref{B}). We note that all the following results also can be obtained  from the solutions with  the electric field (\ref{E}) which has the electric charge $q_e=q$

  \begin{eqnarray}
  \rho(r) =  \frac{ 16 m^4 q^2 }{\kappa^2 ( q^6 + 8 m^3r^3)^{4/3}  }
  =
  -p_r(r)
  \hskip 2 cm
  p_t(r) = \frac{16 q^2m^4( 8 m^3 r^3 - q^6  ) }{ \kappa^2 ( q^6 + 8 m^3r^3)^{7/3}  } \; .
  \end{eqnarray}
  We find the following energy conditions   using the energy density $\rho(r)$, the radial and tangential  pressures  $p_r, p_t$
  \begin{eqnarray}
  DEC_1 & = & \rho   \geq 0,
  \\
  NEC_1 &  = & WEC_1 = \rho + p_{r} = 0\ ,
  \\
  NEC_2 &  = &  WEC_2 = \rho + p_{t} =
  \frac{ 2^8 m^7q^2 r^3}{
  	\kappa^2 ( q^6 + 8 m^3r^3)^{7/3}
  	} \ ,
  \\
  SEC \ \ & = &  \rho + p_r + 2p_t =  \frac{32m^4q^2(8m^3 r^3 -q^6)}{ 	\kappa^2 ( q^6 + 8 m^3r^3)^{7/3} }  \ ,
  \\
  DEC_2 & = & \rho -  p_r = 2\rho \ ,
  \\
  DEC_3 & = & \rho -  p_t = \frac{32 m^4 q^8 }{\kappa^2 ( q^6 + 8 m^3r^3)^{7/3}} \ .
  \end{eqnarray}
 Thus, we see that all the energy conditions are satisfied in the region $r\geq \frac{q^2}{2m}$ for the electrically or magnetically charged solutions.
 But, only the SEC is violated in the central  region $r < \frac{q^2}{2m} $.

\subsection{Regular Black Hole Solution-2}

 Secondly, we take the metric function from  \cite{Balart2} and   \cite{Rodrigues}  with the electric charge $q_e$
 \begin{eqnarray}\label{m2}
 f^2(r) = 1- \frac{2m}{r} e^{-\frac{q_e^2}{2mr}}\ .
 \end{eqnarray}

The Ricci scalar and $R_{ab}\wedge *R^{ab} $    are calculated using the metric function (\ref{m2}) as
 \begin{eqnarray}\label{R2}
 R(r) =  \frac{q_e^4}{2mr^5} e^{-\frac{q_e^2}{2mr}}
 \end{eqnarray}
 \begin{eqnarray}
R_{ab}\wedge*R^{ab} =  (  \frac{24 m^2}{r^6} - \frac{24m q_e^2}{r^7} + \frac{12 q_e^4 }{r^8} -
 \frac{2q_e^6}{mr^9 } +
 \frac{ q_e^8}{8m^2 r^{10} }  )e^{-\frac{q_e^2}{mr}}*1 \ .
 \end{eqnarray}

 When we check their limits
 \begin{eqnarray}
 \lim_{r\rightarrow 0 } R = 0, \hskip 1 cm \lim_{r\rightarrow 0 } R_{ab}\wedge*R^{ab} =0
 \end{eqnarray}
we see that  they are regular at the center of black hole.
 We find the solution of these differential equations (\ref{Yf}) and (\ref{YE}) for this regular metric function (\ref{m2}) as follows
 \begin{eqnarray}
 Y(r)& = & \frac{  2\kappa^2 mr}{ (8mr -q_e^2)e^{-\frac{q_e^2}{2mr}}  }\label{Y2} \\
 E(r) & = &\frac{q_e}{r^2Y(r)} =  \frac{ q_e(  8m r - q_e^2 ) e^{-\frac{q_e^2}{2mr }}   }{ 2m\kappa^2 r^3} \label{E2}\ .
 \end{eqnarray}
 While this electric field is regular at the center,  it has  the following  asymptotic  behavior
  \begin{eqnarray}
 E(r) =  \frac{4q_e}{\kappa^2 r^2} - \frac{5q_e^3}{2 m\kappa^2r^3} +\frac{3q_e^5}{4m^2\kappa^2 r^4} + O(\frac{1}{r^5})\ .
 \end{eqnarray}

  The inverse function of  $R(r)$ in (\ref{R2}) can be found  in terms of Lambert function \cite{Dence} as
  \begin{eqnarray}
  r(R) = -\frac{q_e^2}{10 m  } W^{-1} \left[- \frac{1}{10} \left(  \frac{2 q_e^6R}{ m^4} \right)^{1/5} \right]\ .
  \end{eqnarray}
  Then we rewrite the non-minimal function of this model as
   \begin{eqnarray}
   Y(R) = - \frac{}{}\frac{10^5 \kappa^2 m^4 W^5 \left[ -\frac{1}{10} \left(  \frac{2q_e^6 R}{m^4}\right)^{1/5}\right]}{ 2q_e^6 R \left( 4+ 5W\left[ - \frac{1}{10} \left(  \frac{2q_e^6 R}{m^4}\right)^{1/5}\right] \right) }\ .
   \end{eqnarray}
   From  the duality transformation $YE\rightarrow  - B $,\  \ $q_e\rightarrow -q$, \  $Y \rightarrow \frac{1}{Y}$ \ \cite{sertd}, we can find the  magnetic solution for this same metric (\ref{m2}), which corresponds to dual solution of the electrically charged solution (\ref{Y2}), (\ref{E2}). Then
the resulting magnetic field
  \begin{eqnarray}\label{B2}
  B(r) = \frac{q}{r^2}
  \end{eqnarray}
  and the corresponding non-minimal function
    \begin{eqnarray}
    Y(R) = -    \frac{ 2q^6 R \left( 4+ 5W\left[ - \frac{1}{10} \left(  \frac{2q^6 R}{m^4}\right)^{1/5}\right]\right)   }{ 10^5 \kappa^2 m^4 W^5 \left[ - \frac{1}{10} \left(  \frac{2q^6 R}{m^4}\right)^{1/5}\right] }
    \end{eqnarray}
    constructs the dual solution.
  We calculate the energy density and pressures for the metric function (\ref{m2}) with  the magnetic charge $q$ and the  the magnetic field  (\ref{B2}).
   We note that all the following results also can be obtained  from the solutions with  the electric field (\ref{E2}) which has the electric charge $q_e=q$
   \begin{eqnarray}
   \rho(r) =  \frac{q^2 e^{\frac{-q^2}{2 mr }}}{\kappa^2 r^4}
   =
      -p_r(r)
       \hskip 2 cm
       p_t(r) = \rho(r)  -   \frac{q^4 e^{\frac{-q^2}{2 mr }}}{4 \kappa^2 m r^5} \ .
   \end{eqnarray}

Now we calculate  all the energy conditions for the metric function (\ref{m2})
\begin{eqnarray}
DEC_1 & = & \rho   \geq 0,
\\
NEC_1 &  = & WEC_1 = \rho + p_{r} = 0\ ,
 \\
NEC_2 &  = &  WEC_2 = \rho + p_{t} =
\frac{q^2 e^{\frac{-q^2}{2 mr }}}{4 m \kappa^2 r^5} (8mr -q^2  ) \ ,
\\
SEC \ \ & = &  \rho + p_r + 2p_t =  \frac{q^2 e^{\frac{-q^2}{2 mr }}}{2 m \kappa^2 r^5} (4mr -q^2  )  \ ,
\\
DEC_2 & = & \rho -  p_r = 2\rho \ ,
\\
DEC_3 & = & \rho -  p_t =  \frac{q^4 e^{\frac{-q^2}{2 mr }}}{4 \kappa^2 m r^5} \ .
\end{eqnarray}

  We found that all the energy conditions are satisfied in the region $r \geq \frac{q^2}{4 m}$ for these electrically ($q_e=q$) or magnetically  charged solutions. But,   in the region $ \frac{q^2}{8m } \leq r  <  \frac{q^2}{4m }$ only the SEC is violated by this solution. Furthermore,
  in the central region $r < \frac{q^2}{8m }$  the conditions  $NEC_2, WEC_2$ together with  SEC  are violated.

\section{Conclusion}

\noindent
 We   have investigated  various  regular black hole solutions of the non-minimally coupled $Y(R)F^2$
 theory.  We found  electrically charged   or magnetically charged (after the duality transformation)   regular black hole solutions  which can be obtained  from the non-minimal model with some specific non-minimal functions $Y(R)$.
 We calculated all the energy conditions  for these solutions
 using the effective energy-momentum tensor that comes from  the non-minimally coupled $Y(R)F^2$ term.

The first regular black hole solution violates only the strong energy condition in a central region  $r < \frac{q^2}{2m} $,  inside the event horizon.
This solution is in agreement  with
the singularity theorem of General Relativity \cite{Zaslavski}.
But  the second regular black hole solution violates the weak energy condition  together with   the strong energy condition  in  the region $r < \frac{q^2}{8m} $, while it satisfies  all the energy  conditions in the outer  region  $r \geq \frac{q^2}{4m} $ for the electric or magnetic fields.
The same energy conditions  of the second regular black hole  with an electric field are also  found in \cite{Rodrigues} for a different theory
which is $f(R)$ minimally coupled to the Non-linear electrodynamics.

redpajama_set_name: RedPajamaArXiv