  This article deals with the study of the following nonlinear doubly nonlocal
equation:
  \begin{equation*}
  (-\Delta)^{s_1}_{p}u+\ba(-\Delta)^{s_2}_{q}u = \la a(x)|u|^{\delta-2}u+
b(x)|u|^{r-2} u,\; \text{ in }\; \Om, \; u=0 \text{ on } \mathbb{R}^n\setminus
\Om,
  \end{equation*}
  where $\Om$ is a bounded domain in $\mathbb{R}^n$ with smooth boundary, $1<
\de \le q\leq p<r \leq p^{*}_{s_1}$, with $p^{*}_{s_1}=\ds \frac{np}{n-ps_1}$,
$0<s_2 < s_1<1$, $n> p s_1$ and $\la, \ba>0$ are parameters. Here $a\in
L^{\frac{r}{r-\de}}(\Om)$ and $b\in L^{\infty}(\Om)$ are sign changing
functions. We prove the $L^\infty$ estimates, weak Harnack inequality and
Interior H\"older regularity of the weak solutions of the above problem in the
subcritical case $(r<p_{s_1}^*).$ Also, by analyzing the fibering maps and
minimizing the energy functional over suitable
  subsets of the Nehari manifold, we prove existence and multiplicity of weak
solutions to above
  convex-concave problem. In case of $\de=q$, we show the existence of
solution.
