\section{Introduction and statement of  results}\label{sec1}}
\setcounter{equation}{0}
\par
In this paper,  we study the existence and asymptotic behavior of positive vector solutions for  the following coupled Kirchhoff-type system in $\R^3$:
\begin{equation}\label{11}
\begin{cases}
-\left(a_1+b_1\int_{\R^3}|\nabla u|^2dx\right)\Delta u+\lambda V(x)u=\frac{\alpha}{\alpha+\beta}|u|^{\alpha-2}u|v|^{\beta},&x\in\R^3,\\
-\left(a_2+b_2\int_{\R^3}|\nabla v|^2dx\right)\Delta v+\lambda W(x)v=\frac{\beta}{\alpha+\beta}|u|^{\alpha}|v|^{\beta-2}v,&x\in\R^3,\\
u,v\in \mathcal{D}^{1,2}(\R^3),
\end{cases}\tag{$\mathcal{K}_{\bf{b},\lambda}$}
\end{equation}
where $a_i>0$ are constants, $\lambda,b_i>0$ are  parameters for $i=1,2$, $\alpha>1$, $\beta>1$ satisfy $\alpha+\beta\le4$. We assume that $V(x)$ and $W(x)$ satisfy the following hypotheses:
\vskip4mm
\par\noindent
\begin{tabular}{@{}rp{14.80cm}}
	$(H_1)$ & $V(x),W(x)\in C(\R^3,[0,\infty))$. $\Omega_1=\textup{int } V^{-1}(0)$ and $\Omega_2=\textup{int } W^{-1}(0)$ have smooth boundaries with $\overline{\Omega}_1= V^{-1}(0)$, $\overline{\Omega}_2= W^{-1}(0)$, and ${\Omega}_1\cap\Omega_2\neq\emptyset$;\\
	$(H_2)$&There exists  $c>0$ such that the set $\mathcal{M}=\{x\in\R^3:V(x)W(x)\leq c^2\} $ has finite positive Lebesgue measure.   \\
\end{tabular}
\vskip4mm
\par
This kind of conditions, first proposed by Bartsch--Wang\cite{BW}, implies that $\lambda V$ and $\lambda W$ represent potential wells whose depths are controlled by $\lambda$, and has attracted the attention of several researchers, see \cite{BPW,SJT1,FSX,DYB,AN} and the
references therein. $\lambda V$ and $\lambda W$ are called steep potential wells if $\lambda$ is large enough. We note that condition $(H_2)$ implies that $\Omega_1$ and $\Omega_2$ have finite positive Lebesgue measures.
\par
In recent years, the following scalar Kirchhoff-type equation
\begin{equation}\label{1}
\begin{cases}
-\left(a+b\int_{\R^3}|\nabla u|^2dx\right)\Delta u+ V(x)u=f(x,u),\hspace{1ex}x\in\R^3,\\
u\in H^1(\R^3),
\end{cases}
\end{equation}
where $V:\R^3\to\R$ is a potential function, $f\in C(\R^3\times\R,\R)$ and $a,b>0$ are constants, has been widely studied by many researchers and arises in an interesting physical meaning. Indeed, if we replace $\R^3$ by a bounded domain $\Omega\subset\R^3$ and set $V(x)=0$ in \eqref{1}, then it becomes the following Kirchhoff Dirichlet problem
\begin{equation}\label{12}
\begin{cases}
-\left(a+b\int_{\Omega}|\nabla u|^2dx\right)\Delta u =f(x,u),&x\in\Omega,\\
u=0,&x\in\partial\Omega,
\end{cases}
\end{equation}
which is related to the stationary analogue of the equation
\begin{equation}\label{133}
u_{tt}-\left(a+b\int_{\Omega}|\nabla u|^2dx\right)\Delta u =f(x,u),x\in\Omega,
\end{equation}
presented by Kirchhoff in \cite{Kir} as an extension of the classical D'Alembert’s wave equations for free vibration of a stretched string or membrane, particularly, taking into account the changes in string length caused by oscillations, where $u$ denotes the displacement, $f(x,u)$ is the external force and $a$ is the initial tension while $b$ is related to the intrinsic properties of the string (such as Young's modulus). Similar nonlocal problems also appear in other fields as biological systems, where $u$ describes a process that depends on the average of itself, for example, population density. For more physical background, we refer the readers to the papers\cite{MB,AAS} and the references therein.  After J.L. Lions \cite{Lions} proposed a functional analysis approach to the equation \eqref{133},  the Kirchhoff
Dirichlet problem \eqref{12} received much attention, see \cite{TXH,kz,zzt,AF} and the references therein.
\par
Mathematically, equation \eqref{1} is often referred to as a nonlocal problem as the appearance of the term $(\int_{\R^3}|\nabla u|^2dx)\Delta u$, which implies that equation \eqref{1} is no longer a pointwise identity. This phenomenon causes some mathematical difficulties, which makes tmake the study of \eqref{1} particularly interesting. The qualitative analysis of nontrivial solutions for equation \eqref{1}, including the existence, multiplicity of positive solutions and of sign-changing solutions, has been extensively studied in the last decade, see \cite{Zou,DYB1,gn,TXH1} and the references therein.
\par
More recently, Sun and Wu \cite{SJT1} considered the following scalar case of \eqref{11}
\begin{equation}\label{14}
\begin{cases}
-\left(a+b\int_{\R^3}|\nabla u|^2dx\right)\Delta u+\lambda V(x)u=f(x,u),\hspace{1ex}x\in\R^N,\\
u\in H^1(\R^N),
\end{cases}
\end{equation}
where $N\geq 3$ and $f(x,s)$ is the potential $V$ satisfies the following conditions:
\vskip4mm
\par\noindent
\begin{tabular}{@{}rp{14.80cm}}
	$(V_1)$& $V\in C(\R^N,[0,\infty))$, $\Omega=\textup{int }V^{-1}(0)$ is nonempty and has smooth boundary with $\overline{\Omega}=V^{-1}(0)$;\\
	$(V_2)$& there exists $c>0$ such that the set $\{x\in\R^N|V(x)<c\}$ is nonempty and has finite Lebesgue measure.\\
\end{tabular}
\vskip4mm
\par
By assuming different conditions on the nonlinear term $f(x,u)$, they obtained the existence and the nonexistence of nontrivial solutions of \eqref{14} by using variational methods. Moreover, the authors also explored the asymptotic behavior of nontrivial solutions with the help of the vanishing lemma \cite{Lions1}.
Very recently, Zhang and Du  \cite{zhang} studied \eqref{14} when $N=3$ and $f(x,u)=|u|^{p-2}u$ with $2<p<4$ and proved the existence and asymptotic behavior of positive solutions by combining the truncation technique and the parameter-dependent compactness condition. For more results about Kirchhoff type problems like \eqref{14}, we refer the readers to \cite{DYB,SJT,SJT2} and the references therein.
\par
We observe that there exist extensive papers  in the study of the coupled elliptic systems on the whole space, we refer the readers to \cite{Lv,FSX,Lv1,SHX,cgf,Wu,zhou} and the references therein. Let us state some
known results. When $a_1=a_2=1$, $b_1=b_2=0$  and the three-dimensional space $\R^3$ is replaced by $\R^N$, the system \eqref{11} is reduced to the following semilinear Schr\"{o}dinger system:
\begin{equation}
\begin{cases}
-\Delta u+\lambda V(x)u=\frac{\alpha}{\alpha+\beta}|u|^{\alpha-2}u|v|^{\beta},&x\in\R^N,\\
-\Delta v+\lambda W(x)v=\frac{\beta}{\alpha+\beta}|u|^{\alpha}|v|^{\beta-2}v,&x\in\R^N,\\
u,v\in \mathcal{D}^{1,2}(\R^N),
\end{cases}
\end{equation}
where $N\geq 3$, $\lambda>0$ is a parameter, $\alpha, \beta>1$ and $\alpha+\beta<2^*=2N/(N-2)$. $V(x)$ and $W(x)$ satisfy $(H_1)$ and $(H_2)$.
In \cite{FSX}, Furtado, Silva and Xavier studied the existence and multiplicity of solutions  when the parameter $\lambda$ is large enough. Moreover, they also studied the asymptotic behavior of these solutions when $\lambda\to\infty$ via the Nehari manifold on the limit system.
Shi and Chen \cite{SHX} considered system \eqref{11} when $\lambda=1$ and the nonlinear terms are replaced by $\mu F_u(x,u,v)+|u|^{\tau-2}u$ and $\mu F_v(x,u,v)+|v|^4v$, respectively, where $4<\tau<6$ and $\mu>0$ is a parameter. The authors proved  the existence of ground state solutions when $V(x)=W(x)$ is a continuous and asymptotically periodic potential functional and the nonlinear term $F(x,u,v)$ is 4-superlinear at infinity and satisfies the following monotonicity condition, i.e.,
\vskip4mm
\par\noindent
\begin{tabular}{@{}rp{14.80cm}}
	$(F_1)$& $\frac{F(x,u,v)}{|(u,v)|^4}\to+\infty$ uniformly in $x\in\R^3$ as $(u,v)\to+\infty$ ;\\
	$(F_2)$& $u\to \frac{F_u(x,u,v)}{|u|^3}$ and $v\to \frac{F_v(x,u,v)}{|v|^3}$ are nondecreasing on $(-\infty,0)\cup(0,+\infty)$.\\
\end{tabular}
\vskip4mm
\par
L\"{u} and Xiao \cite{Lv} investigated the following coupled Kirchhoff-type systems
\begin{equation}
\begin{cases}
-\left(a+b\int_{\R^3}|\nabla u|^2dx\right)\Delta u+\lambda V(x)u=\frac{2\alpha}{\alpha+\beta}|u|^{\alpha-2}u|v|^{\beta},&x\in\R^3,\\
-\left(a+b\int_{\R^3}|\nabla v|^2dx\right)\Delta v+\lambda W(x)v=\frac{2\beta}{\alpha+\beta}|u|^{\alpha}|v|^{\beta-2}v,&x\in\R^3,
\end{cases}
\end{equation}
where $\alpha,\beta>2$ with $\alpha+\beta<2^*=6$ and the potentials $V$ and $W$ satisfy the following conditions:
\vskip4mm
\par\noindent
\begin{tabular}{@{}rp{14.80cm}}
	$(H'_1)$ & $V(x),W(x)\in C(\R^3,[0,\infty))$. $\Omega=\textup{int } V^{-1}(0)=\textup{int } W^{-1}(0)$ is nonempty with smooth boundary and $\overline{\Omega}=V^{-1}(0)=W^{-1}(0)$;\\
	$(H'_2)$&There exist  $M_1,M_2>0$ such that the sets $\{x\in\R^3:V(x)\leq M_1\}$ and $\{x\in\R^3:W(x)\leq M_2\}$ have finite positive Lebesgue measures.   \\
\end{tabular}
\vskip4mm
\par
By using Nehari manifold and Mountain Pass Theorem, they obtain the existence and multiplicity of nontrivial vector solutions when the parameter $\lambda$ is  large enough, but don't study the asymptotic behavior of nontrivial vector solutions. In \cite{Lv}, the assumption $4<\alpha+\beta$ is crucial to prove the boundedness of Palais-Smale sequences, which implies the 4-superlinear of the nonlinear term at infinity.
\par
However, to the best knowledge of us, the existence of positive vector solutions for \eqref{11} were established only in the case where $\alpha,\beta>2$ with $\alpha+\beta<6$, and there are no existence  results of positive vector solutions for \eqref{11} in the case where  $\alpha,\beta>1$ with $\alpha+\beta\le4$. We wish to point out that the one key difficulty in studying \eqref{11} comes from the fact that standard variational methods cannot be applied directly, even restricting the energy functional $\mathcal{J}_{\bf{b},\lambda}$(see \eqref{16} below) on the Nehari manifold, because Palais-Smale sequences cannot be bounded. Moreover, the fact that the nonlinear term $F(x,u,v)=|u|^\alpha|v|^\beta$ of \eqref{11} does not satisfy the condition $(F_2)$ prevents us from using Nehari manifold and fibering methods. There is another difficulty that the work space is a subspace of $\mathcal{D}^{1,2}(\R^3)\times\mathcal{D}^{1,2}(\R^3)$, but not of $H^{1}(\R^3)\times H^{1}(\R^3)$, which prevents us from applying Lions' vanishing lemma to study the asymptotic behavior of the solution, Therefore, the methods of \cite{zhang,SHX,Lv,SJT1} are not applicable for \eqref{11} when $2<\alpha+\beta\le4$.
\par
Motivated by the works described above, the purpose of this paper is to study the strongly coupled Kirchhoff-type system \eqref{11} in the case where $\alpha,\beta>1$ with $\alpha+\beta\le4$. More specifically, we shall study the existence of positive vector solutions for \eqref{11} in this case. Moreover, we also explore the asymptotic behavior of these solutions as $\bf{b}\to0$ and $\lambda\to \infty$.
\par
It is worth mentioning that the conditions $(H'_1)$ and $(H'_2)$ are stronger than $(H_1)$ and $(H_2)$, respectively. They are both two different equivalent forms of the conditions $(V_1)$ and $(V_2)$ if we let $V(x)=W(x)$. In \cite{SJT1},
the condition $(V_2)$ is crucial to prove the continuity of the work space embedded into the Sobolev space $H^1(\R^3)$, which can be applied directly in \cite{Lv} under the condition $(H'_2)$, but not in the present paper under the condition $(H_2)$. Hence, we don't expect to find positive vector solutions for \eqref{11} in $H^1(\R^3)$. However, by virtue of the strong coupling of the system and the condition $(H_2)$, we can apply variational methods to study \eqref{11} by considering the energy functional $\mathcal{J}_{\bf{b},\lambda}$ defined in an appropriate closed subspace of $\mathcal{D}^{1,2}(\R^3)\times\mathcal{D}^{1,2}(\R^3)$.
\par
To state our main results, we introduce the following closed subspaces of $\mathcal{D}^{1,2}(\R^3)$ associated with the potentials $V$ and $W$:\[E_{V}=\left\{u\in \mathcal{D}^{1,2}(\R^3):\int_{\R^3}V(x)u^2dx<+\infty \right\}\]
and
\[E_{W}=\left\{u\in \mathcal{D}^{1,2}(\R^3):\int_{\R^3}W(x)u^2dx<+\infty \right\}.\]
Thus the natural space in this paper is the Hilbert space $E=E_{V}\times E_{W}$ with the inner product and norm
\[\langle z,\zeta\rangle=\int_{\R^3}\left[\nabla u\nabla\varphi+\nabla v\nabla\psi+V(x)u\varphi+W(x)v\psi\right]dx,\hspace{2ex} ||z||=\langle z,z\rangle^{1/2}\]
for any $z=(u,v)\in E$ and $\zeta=(\varphi, \psi)\in E$. Throughout this paper, without loss of generality, we always assume that $a_1=a_2=1$ in \eqref{11}. Given $\lambda>0$, we let  $E_\lambda=(E,||\cdot||_{\lambda})$ be equipped with the inner product and norm
\[\langle z,\zeta\rangle_\lambda=\int_{\R^3}\left[\nabla u\nabla\varphi+\nabla v\nabla\psi+\lambda V(x)u\varphi+\lambda W(x)v\psi\right]dx,\hspace{2ex} ||z||_\lambda=\langle z,z\rangle_\lambda^{1/2}.\]
It is easy to see that $||z||\leq||z||_{\lambda}$ for $\lambda\geq1$.
\par
We shall  prove the existence of positive vector solutions of \eqref{11}, to do this we consider $\mathcal{J}_{{\bf{b}},\lambda}:E_\lambda\to\R$ given by
\begin{equation}\label{16}
\mathcal{J}_{\bf{b},\lambda}(z)=\frac{1}{2}||z||_\lambda^2+\frac{1}{4}(b_1|\nabla u|^4_2+b_2|\nabla v|^4_2)-\frac{1}{\alpha+\beta}\int_{ \R^3}|u^+|^\alpha|v^+|^\beta dx.
\end{equation}
In view of the conditions $(H_1)$ and $(H_2)$, the functional $\mathcal{J}_{\bf{b},\lambda}$ is well defined and of class $C^1$. Moreover, it is standard to see that \eqref{11} is variational and its solutions are the critical points of \eqref{16}.
\par
In our first result we study the existence of positive solutions of \eqref{11}. Here we observe that we call $z=(u,v)$ a positive vector function if the functions  $u$ and $v$ are positive a.e. in $\R^3$. We also observe that, if $z=(u,v)\in H_0^1(\Omega_1)\times H_0^1(\Omega_2)$ is a nontrivial vector solution of \eqref{19}, then by the condition $(H_1)$, $z$ is also a nontrivial vector solution of \eqref{11} for any $\lambda>0$. Therefore, we are interested in  positive vector solution of \eqref{11}, and obviously it does not lie in $H_0^1(\Omega_1)\times H_0^1(\Omega_2)$. Moreover, it is reasonable to regard $b_1$ and $b_2$ as parameters in \eqref{11}, since $b_1$ and $b_2$ are related to the intrinsic properties of the string, as already mentioned before.
\vskip4mm
\begin{theorem}\label{the1}
	Suppose that  $\alpha$, $\beta>1$ with $\alpha+\beta\le4$ and conditions $(H_1)$--$(H_2)$ hold. Then there exist four numbers $\Lambda^*>1$, $b_*>0$ and $\delta$, $T>0$ such that for any $\lambda>\Lambda^*$ and $0<b_1+b_2<b_*$, the system \eqref{11} possesses at least one positive vector solution $z_{\bf{b},\lambda}=(u_{\bf{b},\lambda},  v_{\bf{b},\lambda})\in E_\lambda$  satisfying
	\begin{equation}\label{13}
	\delta\le||z_{\bf{b},\lambda}||_\lambda\le T,\text{ for any } b_1,b_2 \text{ and } \lambda.
	\end{equation}
\end{theorem}
\begin{remark}
	\rm{}\eqref{11} has no nontrivial vector solution in $H_0^1(\Omega_1)\times H_0^1(\Omega_2)$ when $b_1+b_2$ is large enough. Indeed, if $z=(u,v)\in H_0^1(\Omega_1)\times H_0^1(\Omega_2)$ is a nontrivial vector solution of \eqref{11}, then similar to the proof of Lemma \ref{L21}, we have
	\[\begin{aligned}
	0=&|\nabla u|_{2,\Omega_1}^2 + |\nabla v|_{2,\Omega_2}^2+b_1|\nabla u|_{2,\Omega_1}^4+b_2|\nabla v|_{2,\Omega_2}^4-\int_{\Omega_1\cap\Omega_2}|u|^\alpha|v|^\beta dx\\
	\geq&|\nabla u|_{2,\Omega_1}^2 + |\nabla v|_{2,\Omega_2}^2+b_1|\nabla u|_{2,\Omega_1}^4+b_2|\nabla v|_{2,\Omega_2}^4-C\left(|\nabla u|_{2,\Omega_1}^2 + |\nabla v|_{2,\Omega_2}^2\right)^{\frac{\alpha+\beta}{2}},\\
	\end{aligned}\]
	where $C=S^{-\frac{\alpha+\beta}{2}}|\Omega_1\cap\Omega_2|^{\frac{6-\alpha-\beta}{6}}>0$ and $|\cdot|$ is the Lebesgue measure. This is a contradiction for $\max\{b_1,b_2\}$ large enough.
\end{remark}
Now we give our main ideas for the proof of Theorem \ref{the1}. To overcome the obstacle of finding bounded Palais-Smale sequence for $\mathcal{J}_{\bf{b},\lambda}$, as in \cite{SJP}, we use a cut-off function $\xi\in C^\infty(\R_+,[0,1])$ satisfying
\begin{equation}\label{199}
\begin{cases}
\xi(t)=1,&t\in[0,1],\\
\xi(t)=0,&t\in[2,+\infty),\\
\xi'(t)\le0,&t\in(0,\infty),\\
||\xi'||_\infty\le2,
\end{cases}
\end{equation}
and  move to study the following truncated functional $\mathcal{J}_{\bf{b},\lambda}^T$: $E_\lambda\to\R$ defined by
\begin{equation}\label{110}
\mathcal{J}_{\bf{b},\lambda}^T(z)=\frac{1}{2}||z||_\lambda^2+\frac{1}{4}\xi\left(\frac{||z||_\lambda^2}{T^2}\right)(b_1|\nabla u|^4_2+b_2|\nabla v|^4_2)-\frac{1}{\alpha+\beta}\int_{ \R^3}|u^+|^\alpha|v^+|^\beta dx,
\end{equation}
for any $T>0$. Firstly, by the Mountain Pass Theorem without (PS) condition \cite{mw},, we can obtain a  Palais-Smale sequence $\{z_n\}=\{(u_n,v_n)\}$ of $\mathcal{J}_{\bf{b},\lambda}^T$ at the mountain pass level $c_{\bf{b},\lambda}^T$ for small $b_1+b_2>0$. Secondly, we have an important upper bound for $c_{\bf{b},\lambda}^T$ that is independent of $T$, $b_1,b_2$ and $\lambda$. By choosing an appropriate $T>0$, we can deduce that, up to a subsequence, $||z_n||_\lambda\le T$ for all $n\in\mathbb{N}$ by restricting $b_1+b_2>0$ small enough, and so $\{z_n\}$ is a bounded Palais-Smale sequence of $\mathcal{J}_{\bf{b},\lambda}$, i.e.,
\par\noindent
$$\sup\limits_{n\in\mathbb{N}}||z_n||_{\lambda}\le T,\hspace{1ex}\mathcal{J}_{\bf{b},\lambda}(z_n)\to c^{T}_{\bf{b},\lambda}\hspace{1ex}\text{and}\hspace{1ex}||\mathcal{J}_{\bf{b},\lambda}'(z_n)||_{E^*_\lambda}\to 0,\hspace{1ex}\text{as}\hspace{1ex}n\to\infty,$$
where $E^*_\lambda$ is the dual space of $E_\lambda$. Finally, to recover the compactness, inspired by \cite{BW1}, we establish the parameter-dependent compactness condition to prove $\mathcal{J}_{\bf{b},\lambda}$ satisfies $\textup{(PS)}_{c^{T}_{\bf{b},\lambda}}$ condition for $\lambda>0$ large enough.
\par
In our next results, we study the asymptotic behavior of the positive vector solutions obtained by Theorem \ref{the1} as $\bf{b}\to\bf{0}$ and $\lambda\to\infty$. It is worth emphasizing there are two usual methods to explore the asymptotic behavior of solutions, one is with the help of the Lions' vanishing lemma \cite{Lions,mw}, see e.g. \cite{BPW,SJT,zhang,zhao}, the other is via the Nehari manifold on the limit system, see \cite{FSX}. As already mentioned before, the Nehari method is not applicable in our case. Moreover, the work space $E_\lambda$ is a subspace of $\mathcal{D}^{1,2}(\R^3)\times\mathcal{D}^{1,2}(\R^3)$, but not of $H^{1}(\R^3)\times H^{1}(\R^3)$, which prevents us from applying Lions' vanishing lemma. By Lemma \ref{L21}, the key inequalities \eqref{46}  and \eqref{38} are obtained, allowing us to overcome this difficulty.
\begin{theorem}\label{the2}
	Let $z_{\bf{b},\lambda}$ be the positive vector solutions of \eqref{11} obtained by Theorem \ref{the1}. Then for any $b_1+b_2\in (0,b_*)$ fixed, $z_{\bf{b},\lambda}\to z_{\bf{b}}$ in $E$ as $\lambda\to\infty$, where $z_{\bf{b}}=(u_{\bf{b}},v_{\bf{b}})\in H_0^1(\Omega_1)\times H_0^1(\Omega_2)$ is a positive vector solution of
	\begin{equation}\label{19}
	\begin{cases}
	-\left(1+b_1\int_{\Omega_1}|\nabla u|^2dx\right)\Delta u=\frac{\alpha}{\alpha+\beta}|u|^{\alpha-2}u|v|^{\beta},&x\in\Omega_1,\\
	-\left(1+b_2\int_{\Omega_2}|\nabla v|^2dx\right)\Delta v=\frac{\beta}{\alpha+\beta}|u|^{\alpha}|v|^{\beta-2}v,&x\in\Omega_2,\\
	u\in H_0^1(\Omega_1),v\in H^1_0(\Omega_2).
	\end{cases}\tag{$\mathcal{K}_{\bf{b},\infty}$}
	\end{equation}
\end{theorem}
\begin{theorem}\label{the3}
	Let $z_{\bf{b},\lambda}$ be the positive vector   solutions of \eqref{11} obtained by Theorem \ref{the1}. Then for any $\lambda\in (\Lambda^*, \infty )$ fixed, $z_{\bf{b},\lambda}\to z_\lambda$ in $E_\lambda$ as $\bf{b}\to\bf{0}$, where $z_\lambda\in E_\lambda$ is a positive vector solution of
	\begin{equation}\label{3res}
	\begin{cases}
	-\Delta u+\lambda V(x)u=\frac{\alpha}{\alpha+\beta}|u|^{\alpha-2}u|v|^{\beta},&x\in\R^3,\\
	-\Delta v+\lambda W(x)v=\frac{\beta}{\alpha+\beta}|u|^{\alpha}|v|^{\beta-2}v,&x\in\R^3,\\
	u,v\in \mathcal{D}^{1,2}(\R^3).
	\end{cases}\tag{$\mathcal{K}_{\bf{0},\lambda}$}
	\end{equation}
\end{theorem}
\begin{theorem}\label{the4}
	Let $z_{\bf{b},\lambda}$ be the positive vector   solutions of \eqref{11} obtained by Theorem \ref{the1}. Then $z_{\bf{b},\lambda}\to z_{\bf{0}}$ in $E$ as  $\bf{b}\to\bf{0}$ and $\lambda\to\infty$, where $z_{\bf{0}}=(u_{\bf{0}},v_{\bf{0}})\in H_0^1(\Omega_1)\times H_0^1(\Omega_2)$ is a positive vector solution of
	\begin{equation}\label{4res}
	\begin{cases}
	-\Delta u=\frac{\alpha}{\alpha+\beta}|u|^{\alpha-2}u|v|^{\beta},&x\in\Omega_1,\\
	-\Delta v=\frac{\beta}{\alpha+\beta}|u|^{\alpha}|v|^{\beta-2}v,&x\in\Omega_2,\\
	u\in H_0^1(\Omega_1),v\in H^1_0(\Omega_2).
	\end{cases}\tag{$\mathcal{K}_{0,\infty}$}
	\end{equation}
\end{theorem}
\begin{remark}
	\text{}\\
	\begin{tabular}{@{}rp{15.05cm}}
		\rm{(\romannumeral1)}&\rm{}Let $b_1(\text{or } b_2)>0$ be a small fixed parameter, similar to the proof of Theorem \ref{the2}, we can also obtain the asymptotic behavior of these solutions obtained by Theorem \ref{the1} and  as $b_2(\text{or } b_1)\to0$ and  $\lambda\to\infty$.\\
		\rm{(\romannumeral2)}&\rm{}Let $b_1(\text{or } b_2)>0$ be a small fixed parameter and $\lambda>0$ be a large fixed-paramete, similar to the proof of Theorem \ref{the3}, we can also obtain the asymptotic behavior of these solutions obtained by Theorem \ref{the1} as $b_2(\text{or } b_1)\to0$.\\
	\end{tabular}
\end{remark}
\par
The remainder of this paper is organized as follows. In Section \ref{sec2}, we give some preliminary results which are crucial throughout the paper. we prove Theorem \ref{the1} in Section \ref{sec3}. The proofs of Theorems \ref{the2}, \ref{the3} and \ref{the4} will be given in Section \ref{sec4}.
\par
Throughout the paper, we use the following notations:\par
$\bullet$ $H^1(\R^3)$ is the usual Sobolev space equipped with the inner product abd norm
$$(u,v)=\int_{\R^3}(\nabla u\nabla v)+uv)dx;\hspace{2ex}||u||^2=\int_{ \R^3}(|\nabla u|^2+u^2)dx.$$\par
$\bullet$ $\mathcal{D}^{1,2}(\R^3) $ is the completion of $C^{\infty}_0(\R^3)$ with respect to the norm
$$||u||_{\mathcal{D}^{1,2}}^2=|\nabla u|_2^2=\int_{ \R^3}|\nabla u|^2dx.$$\par
$\bullet$ $L^r(\Omega)$, $1\le r\le\infty$, $\Omega\subset\R^3$, denotes a Lebesgue space,  the norm in $L^r(\R^3)$ is denoted by\par \hspace{1.25ex} $|u|_{r,\Omega}$,  where $\Omega$ is a proper subset of $\R^3$, by $|u|_r$ where $\Omega=\R^3$.\par
$\bullet$ $S$ is the best constant for the embedding of  $\mathcal{D}^{1,2}(\R^3)\hookrightarrow$ $ L^6(\R^3)$, i.e.,
\begin{equation}\label{S}
|u|_6\le S^{-\frac{1}{2}}||u||_{\mathcal{D}^{1,2}}\quad \forall u\in \mathcal{D}^{1,2}(\R^3).
\end{equation}\par
$\bullet$ For a mensurable function $u$, we denote by $u^+$ and $u^-$ its positive and negative parts\par
\hspace{1.25ex} respectively, given by
\begin{equation*}
u^+=\max\{u,0\},\hspace{1ex} u^-=\max\{-u,0\}.\end{equation*}\par
$\bullet$ $\rightarrow$ and $\rightharpoonup$ denote the strong and weak convergence in the related function space respectively.\par
$\bullet$ $o(1)$ denotes any quantity which tends to zero when $n\to\infty$.\par
$\bullet$ Let $E$ be a real  Banach space and $I\in C^1(E,\R)$. We call that $I$ satisfies the Palais-Smale\par \hspace{1.25ex} condition at $c$ ($\textup{(PS)}_{c}$ condition for short), if any sequence  $\{u_n\}\subset E$ satisfying  $I(u_n)\to c$\par \hspace{1.25ex}  and $||I'(u_n)||_{E^*}\to0$ possesses a convergent subsequence.
\vskip4mm
{\section{Preliminary results}\label{sec2}}
\setcounter{equation}{0}
\vskip2mm
In this section, we give some preliminary lemmas for the proof of Theorem \ref{the1}. It is notable that the following lemma allows us to show the (PS) condition of $\mathcal{J}_{{\bf{b}},\lambda}$ in another way.
\begin{lemma}\label{L21}
	Suppose that  $\alpha$, $\beta>1$ with $\alpha+\beta\le4$ and $(H_1)$--$(H_2)$ hold. Then for any $z=(u,v)\in E_{\lambda}$ and $\lambda>0$, there exists a constant $\hat{c}>0$, independent of $\lambda$, such that
	\begin{equation*}
	\int_{ \R^3}|u|^\alpha|v|^\beta dx \leq \hat{c} ||z||^{\alpha+\beta}.
	\end{equation*}
\end{lemma}
\begin{proof}
	We borrow an idea from \cite[Lemmas 2.1-2.2]{FSX}. It follows from  $(H_2)$ that
	\begin{equation}\label{21}
	\begin{aligned}
	\int_{\R^3}|uv|dx&\le \int_{\mathcal{M}}|uv|dx+\frac{1}{ c}\int_{{\mathcal{M}}^{c}}( V(x)u^2)^{\frac{1}{2}}( W(x)v^2)^{\frac{1}{2}}dx\\
	&\le |u|_6|v|_6|\mathcal{M}|^{\frac{2}{3}}+\frac{1}{c}
	\left(\int_{{\mathcal{M}}^{c}} V(x)u^2dx\right)^{\frac{1}{2}}\left(\int_{{\mathcal{M}}^{c}}W(x)v^2dx\right)^{\frac{1}{2}}\\
	&\le S^{-1}|\nabla u|_2|\nabla v|_2|\mathcal{M}|^{\frac{2}{3}}+\frac{1}{2 c}\int_{\R^3}V(x)u^2+ W(x)v^2dx\\
	&\le \frac{1}{2}\max\left\{S^{-1}|\mathcal{M}|^{\frac{2}{3}}, \frac{1}{ c} \right\}||z||^2,
	\end{aligned}
	\end{equation}
	where we have used the H\"{o}lder inequality and \eqref{S}.
	Taking $r=\frac{6}{8-\alpha-\beta}$, then we can obtain that
	\begin{equation}\label{22}
	\frac{\alpha-1}{6}+\frac{\beta-1}{6}+\frac{1}{r}=1.
	\end{equation}
	Since $1<r\le\frac{3}{2}$, there exists $\theta=\frac{\alpha+\beta-2}{16-2(\alpha+\beta)}\in (0,\frac{1}{4}]$ such that $r=2\theta+1$. We can deduce from H\"{o}lder inequality and \eqref{S} that
	\begin{equation}\label{23}
	\begin{aligned}
	\int_{\R^3}|uv|^rdx&\le \left(\int_{\R^3}|uv|^3dx\right)^{\theta}\left(\int_{\R^3}|uv|dx\right)^{1-\theta}\\
	&\le \left[\frac{1}{2}\left(|u|_6^6+|v|^6_6\right) \right]^\theta\left(\int_{\R^3}|uv|dx\right)^{1-\theta}\\
	&\le \left[\frac{1}{2}S^{-3}\left(|\nabla u|_2^6+|\nabla v|^6_2\right) \right]^\theta\left(\int_{\R^3}|uv|dx\right)^{1-\theta}\\
	&\le S^{-3\theta}||z||^{6\theta}\left(\int_{\R^3}|uv|dx\right)^{1-\theta}.
	\end{aligned}
	\end{equation}
	From \eqref{21}--\eqref{23}, we can derive
	\begin{equation}\label{24}
	\begin{aligned}
	\int_{\R^3}|u|^\alpha |v|^\beta dx&=\int_{\R^3}|u|^{\alpha-1} |v|^{\beta-1}|uv| dx\\
	&\le
	\left(\int_{\R^3}|u|^6dx\right)^\frac{\alpha-1}{6}\left(\int_{\R^3}|v|^6dx\right)^\frac{\beta-1}{6}\left(\int_{\R^3}|uv|^r\right)^\frac{1}{r}\\
	&\le S^{-\frac{\alpha+\beta-2}{2}}||z||^{\alpha+\beta-2}\left(\int_{\R^3}|uv|^r\right)^\frac{1}{r}\\
	&\le S^{-\frac{9\theta}{2\theta+1}}||z||^{\alpha+\beta-2+\frac{6\theta}{r}}\left(\int_{\R^3}|uv|dx\right)^{\frac{1-\theta}{r}}\\
	&\le \left(\frac{1}{2}\max\left\{S^{-1}|\mathcal{M}|^{\frac{2}{3}}, \frac{1}{c} \right\}\right)^{\frac{1-\theta}{2\theta+1}}S^{-\frac{9\theta}{2\theta+1}}||z||^{\alpha+\beta}\\
	&= \hat{c}||z||^{\alpha+\beta},
	\end{aligned}
	\end{equation}
	where $\hat{c}=\left(\frac{1}{2}\max\left\{S^{-1}|\mathcal{M}|^{\frac{2}{3}}, \frac{1}{c} \right\}\right)^{\frac{1-\theta}{2\theta+1}}S^{-\frac{9\theta}{2\theta+1}}$ is independent of $\lambda$. The proof is complete.
\end{proof}
The following lemma shows that the the coupled term $\int_{ \R^3}|u|^\alpha|v|^\beta dx$ has BL-splitting property, which is an another version of Br\'{e}zis-Lieb lemma \cite{BL}.
\begin{lemma}(\label{BL}\cite[Lemma 4.2]{FSX})
	Let $\{z_n\}=\{(u_n,v_n)\}\subset E_\lambda$ be such that $(u_n, v_n)\rightharpoonup(u,v)$ in $E_\lambda$. Then
	$$\lim\limits_{n\to\infty}\int_{ \R^3}(|u_n|^\alpha|v_n|^\beta-|u_n-u|^\alpha|v_n-v|^\beta)dx=\int_{ \R^3}|u|^\alpha|v|^\beta dx.$$
\end{lemma}
\par
To overcome the obstacle of finding bounded Palais-Smale sequence for $\mathcal{J}_{\bf{b},\lambda}$, we first define a cut-off function $\xi\in C^\infty(\R_+,[0,1])$ (see \eqref{199} above). For any $T>0$, we consider the truncated functional $\mathcal{J}_{\bf{b},\lambda}^T:E_\lambda\to\R$
\begin{equation}
\mathcal{J}_{\bf{b},\lambda}^T(z)=\frac{1}{2}||z||_\lambda^2+\frac{1}{4}\xi\left(\frac{||z||_\lambda^2}{T^2}\right)\left(b_1|\nabla u|_2^4+b_2|\nabla v|_2^4\right)-\frac{1}{\alpha+\beta}\int_{ \R^3}|u^+|^\alpha|v^+|^\beta dx.
\end{equation}
It is easy to see that $\mathcal{J}_{\bf{b},\lambda}^T$ is of class $C^1$ with derivative
\begin{equation}\label{26}
\begin{aligned}
\langle (\mathcal{J}_{\bf{b},\lambda}^T)'(z),\zeta\rangle=&
\langle z,\zeta\rangle_\lambda+\frac{1}{2T^2}\xi'\left(\frac{||z||_\lambda^2}{T^2}\right)\left(b_1|\nabla u|_2^4+b_2|\nabla v|_2^4\right)\langle z,\zeta\rangle_\lambda\\
&+\xi\left(\frac{||z||_\lambda^2}{T^2}\right)\left(b_1|\nabla u|_2^2\int_{ \R^3}\nabla u\nabla \varphi dx +b_2|\nabla v|_2^2\int_{ \R^3}\nabla v\nabla\psi dx \right)\\
&-\frac{\alpha}{\alpha+\beta}\int_{\R^3}|u^+|^{\alpha-2}u^+\varphi|v^+|^\beta dx-\frac{\beta}{\alpha+\beta}\int_{\R^3}|u^+|^{\alpha}|v^+|^{\beta-2}v^+\psi dx,
\end{aligned}
\end{equation}
for any $z=(u,v)\in E_\lambda$ and $\zeta=(\varphi, \psi)\in E_\lambda$.
With this penalization, for a properly chosen $T>0$ and $b_1+b_2$ small enough, we are able to find a Palais-Smale sequence $\{z_n\}$ of $\mathcal{J}_{\bf{b},\lambda}^T$ satisfying $||z_n||_\lambda\le T$ and so $\{z_n\}$ is also a Palais-Smale sequence of $\mathcal{J}_{\bf{b},\lambda}$ satisfying $||z_n||_\lambda\le T$.
\par
The following Lemmas \ref{L233}--\ref{L244} imply that $\mathcal{J}_{\bf{b},\lambda}^T$ has the mountain pass geometry.
\begin{lemma}\label{L233}
	Suppose that $\alpha>1$, $\beta>1$ with $\alpha+\beta\le 4$ and $(H_1)$--$(H_2)$ hold. Then there exist $\rho, \eta>0$ such that for any $b_1$, $b_2$, $T>0$ and $\lambda\geq1$, there holds
	\[\inf\{\mathcal{J}^T_{\bf{b},\lambda}:z\in E_\lambda\hspace{1ex} \textup{with}\hspace{1ex} ||z||_{\lambda}=\rho\}\geq\eta.\]
\end{lemma}
\begin{proof}
	It follows from Lemma \ref{L21} that
	\[\mathcal{J}_{\bf{b},\lambda}^T(z)\geq\frac{1}{2}||z||_\lambda^2-\frac{\hat{c}}{\alpha+\beta}||z||_\lambda^{\alpha+\beta}=||z||_\lambda^2\left(\frac{1}{2}-\frac{\hat{c}}{\alpha+\beta}||z||_\lambda^{\alpha+\beta-2}\right),\]
	where $\hat{c}$ is given in Lemma \ref{L21} and independent of $b_1$, $b_2$, $T$ and $\lambda.$ Consequently, by choosing $\rho=\left(\frac{\alpha+\beta}{4\hat{c}}\right)^{\frac{1}{\alpha+\beta-2}}>0$ and $\eta=\frac{1}{4}\rho^2$, we arrive at the desired result.
\end{proof}
\begin{lemma}\label{L244}
	Suppose that $\alpha$, $\beta>1$ with $\alpha+\beta\le4$ and $(H_1)$--$(H_2)$ hold. Then there exist $b^*>0$ and $e=(e_0,e_0)\in C_0^\infty(\Omega_1)\times C_0^\infty(\Omega_2)$ such that for any $T$, $\lambda>0$ and $b_1+b_2\in (0,b^*)$, we have $\mathcal{J}_{\bf{b},\lambda}^T(e)<0$ with $||e||_\lambda>\rho$.
\end{lemma}
\begin{proof}
	We define the functional $\mathcal{I}_\lambda:E_\lambda\to\R$ by
	\[\mathcal{I}_\lambda(z)=\frac{1}{2}||z||_\lambda^2-\frac{1}{\alpha+\beta}\int_{ \R^3}|u^+|^\alpha|v^+|^\beta dx.\]
	Let $e_1=e_2\in C_0^\infty(\Omega_1\cap\Omega_2)$, then we have $V(x)e_1=W(x)e_1=0$ on $\R^3$ and $||(e_1,e_1)||_\lambda=\sqrt{2}|\nabla e_1|_2$. Hence, by $\alpha+\beta>2$, we can see that
	\[\lim\limits_{t\to\infty}\mathcal{I}_{\lambda}(t(e_1,e_1))=\lim\limits_{t\to\infty}\left(t^2\int\limits_{\Omega_1\cap\Omega_2}|\nabla e_1|^2dx-\frac{t^{\alpha+\beta}}{\alpha+\beta}\int\limits_{\Omega_1\cap\Omega_2}|e_1|^{\alpha+\beta} dx\right)=-\infty.\]
	So we further have that there exists $e_0\in C_0^\infty(\Omega_1\cap\Omega_2)$ with $|\nabla e_0|_2>\frac{\rho}{\sqrt{2}}$ and $e_0^+\not\equiv0$, independent of $b_1$, $b_2$, $T$ and $\lambda$, such that $\mathcal{I}_{\lambda}(e)\le -1$. Since
	\[\mathcal{J}_{\bf{b},\lambda}^T(e)=\mathcal{I}_{\lambda}(e)+\frac{b_1+b_2}{4}\xi\left(\frac{||e||_\lambda^2}{T^2}\right)|\nabla e_0|_2^4\le -1+\frac{b_1+b_2}{4}|\nabla e_0|_2^4,\]
	there exists  $b^*>0$, independent of $\lambda$ and $T$, such that $\mathcal{J}_{\bf{b},\lambda}^T(e)<0$ for any $\lambda, T>0$ and $b_1+b_2\in(0,b^*)$. The proof is complete.
\end{proof}
By Lemmas \ref{L233}--\ref{L244} and the Mountain Pass Theorem without (PS) condition \cite{mw}, we can obtain that for any $\lambda\geq 1$, $T>0$ and $b_1+b_2\in(0,b^*)$, there is a $\textup{(PS)}_{c^{T}_{\bf{b},\lambda}}$ sequence $\{z_n\}=\{(u_n,v_n)\}\subset E_\lambda$ such that
\begin{equation}\label{27}
\mathcal{J}_{\bf{b},\lambda}^T\to c^{T}_{\bf{b},\lambda}\geq\eta>0\hspace{1ex}\text{and}\hspace{1ex}||(\mathcal{J}_{\bf{b},\lambda}^T)'(z_n)||_{E_\lambda^*}\to 0,\hspace{1ex}\text{as}\hspace{1ex}n\to\infty,
\end{equation}
where
\[c^{T}_{\bf{b},\lambda}:=\inf\limits_{\gamma \in \Gamma }\max\limits_{t\in[0,1]} \mathcal{J}_{\bf{b},\lambda}^T(\gamma(t)),\]
and
\[\Gamma=\{\gamma\in C([0,1],E_\lambda);\gamma(0)=0,\gamma(1)=e\}.\]
To find a bounded Palais-Smale sequence,  we have an important upper bound for $c^{T}_{\bf{b},\lambda}$, which is the keystone of the truncation technique.
\begin{lemma}\label{L255}
	Suppose that $\alpha$, $\beta>1$ with $\alpha+\beta\le4$ and the conditions $(H_1)$--$(H_2)$ hold. Then for for any $\lambda\geq 1$, $T>0$ and $b_1+b_2\in(0,b^*)$, there exists $D>0$, independent of $b_1$, $b_2$, $T$ and $\lambda$, such tha $c^{T}_{\bf{b},\lambda}\le D$.
\end{lemma}
\begin{proof}
	It follows from $e_0\in C_0^\infty(\Omega_1\cap\Omega_2)$ that
	\[\mathcal{J}_{\bf{b},\lambda}^T(t(e_0,e_0))\le t^2\int\limits_{\Omega_1\cap\Omega_2}|\nabla e_0|^2dx+\frac{b^*t^4}{4}\left(\int\limits_{\Omega_1\cap\Omega_2}|\nabla e_0|^2dx\right)^2-\frac{t^{\alpha+\beta}}{\alpha+\beta}\int\limits_{\Omega_1\cap\Omega_2}|e_0^+|^{\alpha+\beta} dx.\]
	Hence, by the definition of $c^{T}_{\bf{b},\lambda}$, there exists a constant $D>0$ such that
	\[c^{T}_{\bf{b},\lambda}\le \max\limits_{t\in[0,1]}\mathcal{J}_{\bf{b},\lambda}^T(te)\le D.\]
	The proof is complete.
\end{proof}
\par
The following key lemma shows that for a properly  chosen $T> 0$, up to a subsequence, $\{z_n\}$ is $\textup{(PS)}_{c^{T}_{\bf{b},\lambda}}$ sequence of $\mathcal{J}_{\bf{b},\lambda}^T$ satisfying $||z_n||_\lambda\le T$, and so $\{z_n\}$ is also a bounded $\textup{(PS)}_{c^{T}_{\bf{b},\lambda}}$ sequence of  $\mathcal{J}_{\bf{b},\lambda}$, i.e.,
\begin{equation}\label{277}
\mathcal{J}_{\bf{b},\lambda}\to c^{T}_{\bf{b},\lambda}\geq\eta>0\hspace{1ex}\text{and}\hspace{1ex}||\mathcal{J}_{\bf{b},\lambda}'(z_n)||_{E_\lambda^*}\to 0,\hspace{1ex}\text{as}\hspace{1ex}n\to\infty.
\end{equation}
\begin{lemma}\label{L26}
	Suppose that $\alpha$, $\beta>1$ with $\alpha+\beta\le4$, the conditions $(H_1)$--$(H_2)$ hold, $\{z_n\}=\{(u_n,v_n)\}\subset E_\lambda$ is a $\textup{(PS)}_{c^{T}_{\bf{b},\lambda}}$ sequence of  $\mathcal{J}^T_{\bf{b},\lambda}$, and let $T=\sqrt{\frac{2(\alpha+\beta)(D+1)}{\alpha+\beta-2}}$. Then there exists $b_*\in(0,b^*)$ such that for any $\lambda\geq 1$ and $b_1+b_2\in(0,b_*)$, up to a subsequence, $||z_n||_\lambda\le T$.
\end{lemma}
\begin{proof}
	Arguing indirectly, suppose that $||z_n||_{\lambda}>T$ up to a subsequence. We next divide the proof into two separate cases.
	\par
	Case 1: $T<||z_n||_{\lambda}\le\sqrt{2}T$. In this case, by \eqref{26}--\eqref{27} and $\xi'(t)\le0$ for all $t>0$, we have
	\begin{equation}\label{29}
	\begin{aligned}
	c^{T}_{\bf{b},\lambda}=&\lim\limits_{n\to\infty}\left(\mathcal{J}_{\bf{b},\lambda}^T(z_n)-\frac{1}{\alpha+\beta}\langle (\mathcal{J}_{\bf{b},\lambda}^T)'(z_n),z_n\rangle\right)\\
	=&\lim\limits_{n\to\infty}\left((\frac{1}{2}-\frac{1}{\alpha+\beta})||z_n||^2_{\lambda}-(\frac{1}{\alpha+\beta}-\frac{1}{4})\xi\left(\frac{||z||_\lambda^2}{T^2}\right)\left(b_1|\nabla u_n|_2^4+b_2|\nabla v_n|_2^4\right)\right.\\
	&\phantom{=\;\;}
	\left.-\frac{1}{2(\alpha+\beta)T^2}\xi'\left(\frac{||z||_\lambda^2}{T^2}\right)\left(b_1|\nabla u_n|_2^4+b_2|\nabla v_n|_2^4\right)||z_n||^2_{\lambda}\right)\\
	\geq& \liminf\limits_{n\to\infty}\left((\frac{1}{2}-\frac{1}{\alpha+\beta})||z_n||^2_{\lambda}-(\frac{1}{\alpha+\beta}-\frac{1}{4})\left(b_1+b_2\right)||z_n||^4_{\lambda}\right)\\
	\geq& (D+1)-\frac{4(4-\alpha-\beta)(b_1+b_2)}{(\alpha+\beta-2)^2}(D+1)^2,
	\end{aligned}
	\end{equation}
	which contradicts the fact that  $c^{T}_{\bf{b},\lambda}\le D$ by choosing $b_*$ small enough.
	\par
	Case 2: $||z_n||_{\lambda}>\sqrt{2}T$. In this case, by the similar computation to \eqref{29} and using the fact that
	$\xi\left(\frac{||z||_\lambda^2}{T^2}\right)=\xi'\left(\frac{||z||_\lambda^2}{T^2}\right)=0$, we get
	\begin{equation}
	\begin{aligned}
	c^{T}_{\bf{b},\lambda}=&\lim\limits_{n\to\infty}\left(\mathcal{J}_{\bf{b},\lambda}^T(z_n)-\frac{1}{\alpha+\beta}\langle (\mathcal{J}_{\bf{b},\lambda}^T)'(z_n),z_n\rangle\right)\\
	=&\lim\limits_{n\to\infty}\left((\frac{1}{2}-\frac{1}{\alpha+\beta})||z_n||^2_{\lambda}-(\frac{1}{\alpha+\beta}-\frac{1}{4})\xi\left(\frac{||z||_\lambda^2}{T^2}\right)\left(b_1|\nabla u_n|_2^4+b_2|\nabla v_n|_2^4\right)\notag\right.\\
	&\phantom{=\;\;}
	\left.-\frac{1}{2(\alpha+\beta)T^2}\xi'\left(\frac{||z||_\lambda^2}{T^2}\right)\left(b_1|\nabla u_n|_2^4+b_2|\nabla v_n|_2^4\right)||z_n||^2_{\lambda}\right)\\
	\geq& 2(D+1),
	\end{aligned}
	\end{equation}
	which is again a contradiction. Hence $||z_n||_{\lambda}\le T$ and the proof is complete.
\end{proof}
\vskip2mm
{\section{Existence of positive vector solutions }\label{sec3}}
\setcounter{equation}{0}
\vskip2mm
To find the critical points of $\mathcal{J}_{\bf{b},\lambda}$, we first  establish the following
parameter-dependent compactness condition for the functional $\mathcal{J}_{\bf{b},\lambda}$.
\vskip2mm
\begin{proposition}\label{L28}
	Suppose that $\alpha>1$, $\beta>1$ with $\alpha+\beta\le4$ and the conditions $(V_1)$--$(V_2)$ hold, and let $T=\sqrt{\frac{2(\alpha+\beta)(D+1)}{\alpha+\beta-2}}$. Then there exists $\Lambda^*> 1$ such that  $\mathcal{J}_{\bf{b},\lambda}$ satisfies $\textup{(PS)}_{c^{T}_{\bf{b},\lambda}}$ condition in $E_\lambda$ for any $b_1+b_2\in(0,b_*)$ and $\lambda\in(\Lambda^*,\infty)$.
\end{proposition}
\begin{proof}
	Let $\{z_n\}=\{(u_n,v_n)\}\subset E$ be a $\textup{(PS)}_{c^{T}_{\bf{b},\lambda}}$ sequence for $\mathcal{J}_{\bf{b},\lambda}$. From Lemma \ref{L26}, we know that $\{z_n\}$ is bounded in $E_\lambda$, i.e.,
	\[||z_n||_\lambda\le T=\sqrt{\frac{2(\alpha+\beta)(D+1)}{\alpha+\beta-2}}.\]
	Then there exist a subsequence $\{z_n\}$ and $z=(u,v)$ in $E_\lambda$ such that
	\begin{equation}\label{28}
	\begin{cases}
	(u_n,v_n)\rightharpoonup(u,v),& \hspace{1ex}\textup{in} \hspace{1ex}E_\lambda;\\
	(u_n,v_n)\to(u,v),& \hspace{1ex}\textup{a.e. in} \hspace{1ex}\R^3;\\
	(u_n,v_n)\to(u,v),& \hspace{1ex}\textup{in} \hspace{1ex}L_{loc}^{s_1}(\R^3)\times L_{loc}^{s_2}(\R^3),\hspace{1ex} 2\le s_1, s_2<6.\\
	\end{cases}
	\end{equation}
	Now we prove that $z_n\to z$ in $E_\lambda$. Let $\tilde{z}_n=z_n-z$ with $\tilde{u}_n=u_n-u$ and $\tilde{v}_n=v_n-v$. By \eqref{28}, we have
	\[||z||_\lambda\le \liminf\limits_{n\to\infty}||z_n||_\lambda\le T,\]
	and then
	\begin{equation}\label{211}
	||\tilde{z}_n||_\lambda=||z_n-z||_\lambda\le 2T.
	\end{equation}
	It follows from $(H_2)$ that
	\[\begin{aligned}
	\int_{\R^3}|\tilde{u}_n\tilde{v}_n|dx&=\int_{\mathcal{M}^c}|\tilde{u}_n\tilde{v}_n|dx+\int_{\mathcal{M}}|\tilde{u}_n\tilde{v}_n|dx\\
	&\le\frac{1}{2\lambda c}\int_{\mathcal{M}^c}\lambda V(x)\tilde{u}_n^2+\lambda W(x)\tilde{v}_n^2dx+ |\tilde{u}_n|_{5,\mathcal{M}}|\tilde{v}_n|_{5,\mathcal{M}}|\mathcal{M}|^{\frac{3}{5}}\\
	&\le \frac{1}{2\lambda c}||\tilde{z}_n||_\lambda^2+o(1),
	\end{aligned}\]
	this together with Lemma \ref{L21} implies that
	\begin{equation}\label{2133}
	\begin{aligned}
	\int_{\R^3}|\tilde{u}_n|^\alpha|\tilde{v}_n|^\beta dx
	&\le S^{-\frac{9\theta}{2\theta+1}}||\tilde{z}_n||_\lambda^{\alpha+\beta-2+\frac{6\theta}{2\theta+1}}\left(\int_{\R^3}|\tilde{u}_n\tilde{v}_n|dx\right)^{\frac{1-\theta}{2\theta+1}}\\
	&\le \left(\frac{1}{2\lambda c}\right)^{\frac{1-\theta}{2\theta+1}}S^{-\frac{9\theta}{2\theta+1}}||\tilde{z}_n||_\lambda^{\alpha+\beta}+o(1),
	\end{aligned}
	\end{equation}
	where $\theta=\frac{\alpha+\beta-2}{16-2(\alpha+\beta)}\in (0,\frac{1}{4}]$.
	Moreover, there exist $A$, $B\in\R$ such that
	\[\int_{ \R^3}|\nabla u_n|^2dx\to A^2,\hspace{2ex}\int_{ \R^3}|\nabla v_n|^2dx\to B^2\]
	and
	\[\int_{ \R^3}|\nabla u|^2dx\le A^2,\hspace{2ex}\int_{ \R^3}|\nabla v|^2dx\le B^2.\]
	As in \cite{FSX}, it follows from \eqref{28} and the Lebesgue Dominated Convergence Theorem that
	\begin{equation}\label{212}
	\lim\limits_{n\to\infty}\int_{ \R^3}|u_n|^{\alpha-2}u_n\varphi|v_n|^\beta dx=\int_{ \R^3}|u|^{\alpha-2}u\varphi|v|^\beta dx,\textup{ for any } \varphi\in C_0^\infty(\R^3)
	\end{equation}
	and
	\begin{equation}\label{213}
	\lim\limits_{n\to\infty}\int_{ \R^3}|u_n|^{\alpha}|v_n|^{\beta-2}v_n\psi dx=\int_{ \R^3}|u|^{\alpha}|v|^{\beta-2}v\psi dx,\textup{ for any } \psi\in C_0^\infty(\R^3).
	\end{equation}
	Then by $ \mathcal{J}_{\bf{b},\lambda}'(z_n)\to0$ and \eqref{212}--\eqref{213} , we can get that for any  $\zeta=(\varphi,\psi)\in C_0^\infty(\R^3)\times C_0^\infty(\R^3)$,
	\begin{equation}\label{214}
	\begin{aligned}
	&\langle z,\zeta\rangle_\lambda-\frac{1}{\alpha+\beta}\left(\alpha\int_{\R^3}|u^+|^{\alpha-2}u^+\varphi|v^+|^\beta dx+{\beta}\int_{\R^3}|u^+|^{\alpha}|v^+|^{\beta-2}v^+\psi dx\right) \\
	&+b_1A^2\int_{ \R^3}\nabla u\nabla \varphi dx+b_2B^2\int_{ \R^3}\nabla v\nabla\psi dx=0.
	\end{aligned}
	\end{equation}
	Taking $\zeta=z$ in \eqref{214}, we have
	\begin{equation}\label{215}
	\begin{aligned}
	||z||^2_\lambda
	+b_1A^2|\nabla u|_2^2+b_2B^2|\nabla v|_2^2-\int_{\R^3}|u^+|^{\alpha}|v^+|^\beta dx =0.
	\end{aligned}
	\end{equation}
	Note that
	\begin{equation}
	\begin{aligned}
	||z_n||^2_\lambda
	+b_1|\nabla u_n|_2^4+b_2|\nabla v_n|_2^4 dx-\int_{\R^3}|u_n^+|^{\alpha}|v_n^+|^\beta dx =o(1).
	\end{aligned}
	\end{equation}
	By the Br\'{e}zis--Lieb Lemma, we have that
	\begin{equation}
	||\tilde{z}_n||_\lambda^2=||z_n||_\lambda^2-||z||_\lambda^2+o(1),
	\end{equation}
	\begin{equation}\label{218}
	|\nabla \tilde{u}_n|_2^2=|\nabla u_n|_2^2-|\nabla u|_2^2+o(1),\hspace{2ex} |\nabla \tilde{v}_n|_2^2=|\nabla v_n|_2^2-|\nabla v|_2^2+o(1).
	\end{equation}
	It follows from \eqref{215}--\eqref{218} and Lemma \ref{BL} that
	\begin{equation}\label{221}
	\begin{aligned}
	o(1)=&||z_n||^2_\lambda
	+b_1|\nabla u_n|_2^4+b_2|\nabla v_n|_2^4 dx-\int_{\R^3}|u_n^+|^{\alpha}|v_n^+|^\beta dx \\
	&- ||z||^2_\lambda
	-b_1A^2|\nabla u|_2^2-b_2B^2|\nabla v|_2^2+\int_{\R^3}|u^+|^{\alpha}|v^+|^\beta dx \\
	=&||\tilde{z}_n||^2_\lambda +b_1A^2|\nabla \tilde{u}_n|_2^2+b_2B^2|\nabla \tilde{v}_n|_2^2-\int_{\R^3}|\tilde{u}_n^+|^{\alpha}|\tilde{v}_n^+|^\beta dx.
	\end{aligned}
	\end{equation}
	Then by \eqref{24}, \eqref{211}--\eqref{2133} and \eqref{221}, we get that
	\begin{equation}
	\begin{aligned}
	o(1)=&||\tilde{z}_n||^2_\lambda-\int_{\R^3}|\tilde{u}_n^+|^{\alpha}|\tilde{v}_n^+|^\beta dx +b_1A^2|\nabla \tilde{u}_n|_2^2+b_2B^2|\nabla \tilde{v}_n|_2^2\\
	\geq& ||\tilde{z}_n||^2_\lambda-\int_{\R^3}|\tilde{u}_n|^{\alpha}|\tilde{v}_n|^\beta dx\\
	=&||\tilde{z}_n||^2_\lambda-\left(\int_{\R^3}|\tilde{u}_n|^{\alpha}|\tilde{v}_n|^\beta dx\right)^{\frac{\alpha+\beta-2}{\alpha+\beta}}
	\left(\int_{\R^3}|\tilde{u}_n|^{\alpha}|\tilde{v}_n|^\beta dx\right)^{\frac{2}{\alpha+\beta}}\\
	\geq &||\tilde{z}_n||^2_\lambda-(\hat{c})^{\frac{\alpha+\beta-2}{\alpha+\beta}}(2T)^{\alpha+\beta-2}\left(\frac{1}{2\lambda c}\right)^{\frac{1-\theta}{8\theta+1}}S^{-\frac{9\theta}{8\theta+1}}||\tilde{z}_n||^2_\lambda,
	\end{aligned}
	\end{equation}
	which implies that there exists $\Lambda^*> 1$ such that $\tilde{z}_n\to 0$ in $E_\lambda$ for all $\lambda>\Lambda^*$ and $0<b_1+b_2<b_*$. Consequently, this completes the proof.
\end{proof}
\begin{proof}[\rm{}\textbf{Proof of Theorem \ref{the1}}]
	Let $T=\sqrt{\frac{2(\alpha+\beta)(D+1)}{\alpha+\beta-2}}$.  By Lemmas \ref{L233}--\ref{L244}, we see that $\mathcal{J}_{\bf{b},\lambda}^T$ possesses the Mountain Pass geometry. Using a version of the Mountain Pass Theorem without (PS) condition \cite{mw}, we can obtain that for any $\lambda\geq 1$, and $0<b_1+b_2<b^*$, $\mathcal{J}_{\bf{b},\lambda}^T$ has a $\textup{(PS)}_{c^{T}_{\bf{b},\lambda}}$ sequence $\{z_n\}=\{(u_n,v_n)\}\subset E_\lambda$. By Lemmas \ref{L255}--\ref{L26}, we can get that there exists $b_*\in(0,b^*)$ such that for any $\lambda\geq 1$ and $b_1+b_2\in(0,b_*)$, up to a  subsequence, $\{z_n\}$ is a $\textup{(PS)}_{c^{T}_{\bf{b},\lambda}}$ sequence of $\mathcal{J}_{\bf{b},\lambda}$ satisfying
	\[\sup\limits_{n\in\mathbb{N}}||z_n||_{\lambda}\le T,\hspace{1ex}\mathcal{J}_{\bf{b},\lambda}(z_n)\to c^{T}_{\bf{b},\lambda}\hspace{1ex}\text{and}\hspace{1ex}||\mathcal{J}_{\bf{b},\lambda}'(z_n)||_{E^*_\lambda}\to 0,\hspace{1ex}\text{as}\hspace{1ex}n\to\infty.\]
	By Proposition \ref{L28}, there exists $\Lambda^*> 1$ such that  $\mathcal{J}_{\bf{b},\lambda}$ satisfies $\textup{(PS)}_{c^{T}_{\bf{b},\lambda}}$ condition in $E_\lambda$ for any $b_1+b_2\in(0,b_*)$ and $\lambda\in(\Lambda^*,\infty)$. Thus, we may assume that $z_n\to z_{\bf{b},\lambda}$ in $E_\lambda$, and then
	\[||z_{\bf{b},\lambda}||_{\lambda}\le T,\hspace{1ex}\mathcal{J}_{\bf{b},\lambda}(z_{\bf{b},\lambda})= c^{T}_{\bf{b},\lambda}\geq\eta>0\hspace{1ex}\text{and}\hspace{1ex}\mathcal{J}_{\bf{b},\lambda}'(z_{\bf{b},\lambda})= 0.\]
	Since
	\[\begin{aligned}
	0=&\langle \mathcal{J}_{\bf{b},\lambda}'(z_{\bf{b},\lambda}),(u^-_{\bf{b},\lambda},v^-_{\bf{b},\lambda})\rangle\\
	=&||z^-_{\bf{b},\lambda}||^2_\lambda+b_1|\nabla u|^2_2\cdot|\nabla u^-|^2_2+b_2|\nabla v|^2_2\cdot|\nabla v^-|^2_2\\
	\geq&  ||z^-_{\bf{b},\lambda}||^2_\lambda,
	\end{aligned}\]
	we get that $u_{\bf{b},\lambda},v_{\bf{b},\lambda}\geq 0$ in $\R^3$. By using the strong maximum principle in each equation of \eqref{11}, we have that $u_{\bf{b},\lambda},v_{\bf{b},\lambda}> 0$ in $\R^3$. Moreover, note that $z_{\bf{b},\lambda}\neq \bf{0}$ and
	\begin{equation}\label{3133}
	\begin{aligned}
	0=&\langle \mathcal{J}_{\bf{b},\lambda}'(z_{\bf{b},\lambda}),z_{\bf{b},\lambda}\rangle\\
	=&||z_{\bf{b},\lambda}||^2_\lambda
	+b_1|\nabla u_{\bf{b},\lambda}|_2^4+b_2|\nabla v_{\bf{b},\lambda}|_2^4 -\int_{\R^3}|u_{\bf{b},\lambda}|^{\alpha}|u_{\bf{b},\lambda}|^\beta dx\\
	\geq &||z_{\bf{b},\lambda}||^2_\lambda-\hat{c}||z_{\bf{b},\lambda}||^{\alpha+\beta},
	\end{aligned}
	\end{equation}
	then, there exists $\delta >0$, independent of $b_1,b_2$ and $\lambda$, such that $||z_{\bf{b},\lambda}||_\lambda>\delta$ for any $b_1,b_2$ and $\lambda$. The proof is complete.
\end{proof}
\vskip2mm
{\section{Asymptotic behavior of positive vector solutions }\label{sec4}}
\setcounter{equation}{0}
\vskip2mm
In this section, we study the asymptotic behavior of positive vector solutions for \eqref{11} and and give the proofs of Theorem \ref{the2}, \ref{the3} and \ref{the4}.
\par
\vskip2mm
\begin{proof}[\rm{}\textbf{Proof of Theorem \ref{the2}}]
	For any sequence $\lambda_n\to\infty$, let $b_1+b_2\in(0,b_*)$ be fixed and $z_n:=z_{\textbf{b},\lambda_n}$ be the  positive vector  solution of $\mathcal{K}_{\textbf{b},\lambda_n}$ obtained in Theorem \ref{the1}. By \eqref{13}, we can get that for any $n\in\mathbb{N}$,
	\begin{equation}\label{31}
	0<\delta\le ||z_n||_{\lambda_n}\le T.
	\end{equation}
	Then there exist a subsequence $\{z_n\}=\{(u_n,v_n)\}$ and $z_{\bf{b}}=(u_{\bf{b}},v_{\bf{b}})$ in $E$ such that
	\begin{equation}\label{322}
	\begin{cases}
	(u_n,v_n)\rightharpoonup(u_{\bf{b}},v_{\bf{b}}),& \hspace{1ex}\textup{in} \hspace{1ex}E;\\
	(u_n,v_n)\to(u_{\bf{b}},v_{\bf{b}}),& \hspace{1ex}\textup{a.e. in} \hspace{1ex}\R^3;\\
	(u_n,v_n)\to(u_{\bf{b}},v_{\bf{b}}),& \hspace{1ex}\textup{in} \hspace{1ex}L_{loc}^{s_1}(\R^3)\times L_{loc}^{s_2}(\R^3),\hspace{1ex} 2\le s_1, s_2<6.\\
	\end{cases}
	\end{equation}
	We first prove $u_{\bf{b}}=0$ in $\Omega_1^c$. For $m\in\mathbb{N}$, set $C_m=\left\{ x\in\R^3: V(x)\geq\frac{1}{m} \right\}$, it follows from \eqref{31} that
	\[\int_{C_m}u_n^2dx\le\frac{m}{\lambda_n}\int_{C_m}\lambda_nV(x)u_n^2dx\le\frac{m}{\lambda_n}||z_n||_{\lambda_n}^2\]
	and thus
	\[\int_{C_m}u_{\bf{b}}^2dx\le \liminf\limits_{n\to \infty }\int_{C_m}u_n^2dx=0,\]
	which implies that $u_{\bf{b}}\equiv 0$ in $\Omega_1^c=\bigcup_{m=1}^{n} C_m$. Since $\Omega_1$ has smooth boundary, we get that $u_{\bf{b}}\in H_0^1(\Omega_1)$. Similarly, we have $v_{\bf{b}}\in H_0^1(\Omega_2)$.
	\par
	For any $\varphi\in C_0^\infty(\Omega_1)$, using $(\varphi,0)$ as a test function, we have
	\begin{equation}
	\left(1+b_1|\nabla u_n|^2_2\right)\int_{\R^3}\nabla u_n\nabla\varphi dx=\frac{\alpha}{\alpha+\beta}\int_{\R^3}|u_n^+|^{\alpha-2}u_n^+\varphi|v_n^+|^\beta dx,
	\end{equation}
	Since $\varphi$ has compact support, using the same argument as that in \eqref{212} and taking limit in the above equality, we get
	\begin{equation}\label{34}
	\left(1+b_1|\nabla u_n|^2_2\right)\int\limits_{\Omega_1\cup\Omega_2}\nabla u_{\bf{b}}\nabla\varphi dx=\frac{\alpha}{\alpha+\beta}\int\limits_{\Omega_1\cup\Omega_2}|u_{\bf{b}}^+|^{\alpha-2}u_{\bf{b}}^+\varphi|v_{\bf{b}}^+|^\beta dx,\hspace{1ex} \forall \varphi\in C_0^\infty(\Omega_1).
	\end{equation}
	Similarly, for any $\psi\in C_0^\infty(\Omega_2)$, we have
	\begin{equation}\label{35}
	\left(1+b_2|\nabla v_n|^2_2\right)\int\limits_{\Omega_1\cup\Omega_2}\nabla v_{\bf{b}}\nabla\psi dx=\frac{\beta}{\alpha+\beta}\int\limits_{\Omega_1\cup\Omega_2}|u_{\bf{b}}^+|^{\alpha}|v_{\bf{b}}^+|^{\beta-2}v_{\bf{b}}^+\psi dx,\hspace{1ex} \forall \psi\in C_0^\infty(\Omega_2).
	\end{equation}
	\par
	Next we show that $z_n\to z_{\bf{b}}$ in $E$.
	Denote $\tilde{u}_n=u_n-u_{\bf{b}}$ and $\tilde{v}_n=v_n-v_{\bf{b}}$, then $\tilde{z}_n=z_n-z_{\bf{b}}$. Similiar to \eqref{2133}, we have
	\begin{equation}\label{46}
	\begin{aligned}
	\int_{\R^3}|\tilde{u}_n|^\alpha|\tilde{v}_n|^\beta dx
	&\le S^{-\frac{9\theta}{2\theta+1}}||\tilde{z}_n||_{\lambda_n}^{\alpha+\beta-2+\frac{6\theta}{2\theta+1}}\left(\int_{\R^3}|\tilde{u}_n\tilde{v}_n|dx\right)^{\frac{1-\theta}{2\theta+1}}\\
	&\le \left(\frac{1}{2\lambda_n c}\right)^{\frac{1-\theta}{2\theta+1}}S^{-\frac{9\theta}{2\theta+1}}||\tilde{z}_n||_{\lambda_n}^{\alpha+\beta}+o(1)=o(1),
	\end{aligned}
	\end{equation}
	this together with Lemma \ref{BL} implies that
	\begin{equation}\label{38}
	\lim\limits_{n\to\infty}\int_{\R^3}|u_n^+|^{\alpha}|v_n^+|^{\beta} dx=\int_{\R^3}|u_{\bf{b}}^+|^{\alpha}|v_{\bf{b}}^+|^{\beta} dx=\int\limits_{\Omega_1\cup\Omega_2}|u_{\bf{b}}^+|^{\alpha}|v_{\bf{b}}^+|^{\beta} dx.
	\end{equation}
	By \eqref{212}--\eqref{213} and
	\[\langle \mathcal{J}_{b,\lambda_n}'(z_n),z_n\rangle=\langle \mathcal{J}_{b,\lambda_n}'(z_n),(u_{\bf{b}},0)\rangle=\langle \mathcal{J}_{b,\lambda_n}'(z_n),(0,v_{\bf{b}})\rangle=o(1),\]
	we have
	\begin{equation}
	||z_n||_{\lambda_n}^2+b_1|\nabla u_n|^4_2+b_2|\nabla v_n|^4_2=\int_{\R^3}|u_n^+|^{\alpha}|v_n^+|^{\beta} dx+o(1),
	\end{equation}
	\begin{equation}
	||u_{\bf{b}}||^2+b_1|\nabla u_n|^2_2|\nabla u_{\bf{b}}|^2_2=\frac{\alpha}{\alpha+\beta}\int_{\R^3}|u_{\bf{b}}^+|^{\alpha}|v_{\bf{b}}^+|^{\beta} dx+o(1),
	\end{equation}
	\begin{equation}\label{311}
	||v_{\bf{b}}||^2+b_2|\nabla v_n|^2_2|\nabla v_{\bf{b}}|^2_2=\frac{\beta}{\alpha+\beta}\int_{\R^3}|u_{\bf{b}}^+|^{\alpha}|v_{\bf{b}}^+|^{\beta} dx+o(1).
	\end{equation}
	Up to a subsequence, we may assume that there exists $l_i>0(i=1,2,3)$ such that
	\[||z_n||_{\lambda_n}^2\to l_1,\hspace{1ex}|\nabla u_n|^2_2\to l_2,\hspace{1ex}|\nabla v_n|^2_2\to l_3.\]
	It follows from \eqref{38}--\eqref{311} that
	\[\begin{aligned}
	l_1+b_1l_2^2+b_2l_3^2&=||u_{\bf{b}}||^2+b_1l_2|\nabla u_{\bf{b}}|^2_2+||v_{\bf{b}}||^2+b_2l_3|\nabla v_{\bf{b}}|^2_2\\
	&\le ||z_{\bf{b}}||^2+b_1l_2^2+b_2l_3^2.
	\end{aligned}\]
	Then, $l_1\le ||z_{\bf{b}}||^2$. By the weakly lower semi-continuity of $||\cdot||$, we have
	\[||z_{\bf{b}}||^2\le \liminf\limits_{n\to\infty}||z_n||^2\le \limsup\limits_{n\to\infty}||z_n||^2\le \lim\limits_{n\to\infty}||z_n||^2_{\lambda_n}=l_1\le||z_{\bf{b}}||^2,\]
	which implies that $z_n\to z_{\bf{b}}$ in $E$.
	\par
	Finally, we are going to show that $z_{\bf{b}}$ is a positive vector solution of $\mathcal{K}_{\bf{b},\infty}$. It follows from \eqref{34}--\eqref{35} that for any $\varphi\in C_0^\infty(\Omega_1)$ and $\psi\in C_0^\infty(\Omega_2)$,
	\begin{equation}
	\left(1+b_1\int_{\Omega_1}|\nabla u_{\bf{b}}|^2 dx\right)\int_{\Omega_1}\nabla u_{\bf{b}}\nabla\varphi dx=\frac{\alpha}{\alpha+\beta}\int_{\Omega_1}|u_{\bf{b}}^+|^{\alpha-2}u_{\bf{b}}^+\varphi|v_{\bf{b}}^+|^\beta dx
	\end{equation}
	and
	\begin{equation}
	\left(1+b_2\int_{\Omega_2}|\nabla v_{\bf{b}}|^2 dx\right)\int_{\Omega_2}\nabla v_{\bf{b}}\nabla\psi dx=\frac{\beta}{\alpha+\beta}\int_{\Omega_2}|u_{\bf{b}}^+|^{\alpha}|v_{\bf{b}}^+|^{\beta-2}v_{\bf{b}}^+\psi dx.
	\end{equation}
	Thus, $z_{\bf{b}}=(u_{\bf{b}},v_{\bf{b}})$ is a nonnegative vector solution of $\mathcal{K}_{\bf{b},\infty}$ by the density of $C_0^\infty(\Omega_i)$ in $H_0^1(\Omega_i)(i=1,2)$. By \eqref{31} and $z_n\to z_{\bf{b}}$ in $E$, we have
	\[||z_{\bf{b}}||=\lim\limits_{n\to\infty}||z_n||_{\lambda_n}\geq \delta>0,\]
	and so $u_{\bf{b}},v_{\bf{b}}\neq0$.
	We can deduce that $u_{\bf{b}}>0$ in $\Omega_1$ and $v_{\bf{b}}>0$ in $\Omega_2$ by the strong maximum principle. The proof is complete.
\end{proof}
\begin{proof}[\rm{}\textbf{Proof of Theorem \ref{the3}}]
	For any sequence $b_n\to0$, let $\lambda\in(\Lambda^*,\infty)$ be fixed and $z_n:=z_{\textbf{b}_n,\lambda}$ be the  positive vector  solution of $\mathcal{K}_{\textbf{b}_n,\lambda}$ obtained in Theorem \ref{the1}. By \eqref{13}, we can get that for any $n\in\mathbb{N}$,
	\begin{equation}\label{313}
	0<\delta\le ||z_n||_{\lambda}\le T,
	\end{equation}
	i.e., $\{z_n\}$ is bounded in $E_\lambda$. Then there exist a subsequence $\{z_n\}$ and $z_\lambda=(u_\lambda,v_\lambda)$ in $E_\lambda$ such that $z_n\rightharpoonup z_\lambda$ in $E_\lambda$. By a similar argument as used in the proof of Proposition \ref{L28}, we get that $z_n\to z_\lambda$ in $E_\lambda$.
	\par
	Now, it is enough to show that $z_\lambda$ is a positive vector solution of \ref{3res}. Note that $\mathcal{J}_{b_n,\lambda}'(z_n)\to0$, then for any $\zeta=(\varphi,\psi)\in E_\lambda$, using $(\varphi,0)$ as a test function, we have
	\begin{equation}\label{314}
	\int_{\R^3}(\nabla u_\lambda\nabla\varphi+\lambda V(x)u_\lambda\varphi) dx=\frac{\alpha}{\alpha+\beta}\int_{\R^3}|u_\lambda^+|^{\alpha-2}u_\lambda^+\varphi|v_\lambda^+|^\beta dx.
	\end{equation}
	Similarly, using $(0,\psi)$ as a test function, we have
	\begin{equation}\label{315}
	\int_{\R^3}(\nabla v_\lambda\nabla\psi+\lambda W(x)v_\lambda\psi) dx=\frac{\beta}{\alpha+\beta}\int_{\R^3}|u_n^+|^{\alpha}|v_n^+|^{\beta-2}v_\lambda^+\psi dx.
	\end{equation}
	Thus, $z_\lambda=(u_\lambda,v_\lambda)$ is a nonnegative vector solution of \ref{3res}. By \eqref{313} and $z_n\to z_\lambda$ in $E_\lambda$, we have
	\[||z_\lambda||_\lambda=\lim\limits_{n\to\infty}||z_n||_{\lambda}\geq \delta>0,\]
	and so $u_\lambda,v_\lambda\neq0$.
	We can deduce that $u_\lambda>0$ in $\R^3$ and $v_\lambda>0$ in $\R^3$ by the strong maximum principle. The proof is complete.
\end{proof}
\begin{proof}[\rm{}\textbf{Proof of Theorem \ref{the4}}] The proof of Theorem \ref{the4} is essentially same as the proof of Theorem \ref{the2}, so we omit it here.
\end{proof}

\par\noindent
\textbf{Acknowledgments}
\vskip2mm
This work is  supported by National Natural Science Foundation of China (No. 12071486).

\vskip6mm


\section{Introduction}\label{sec:intro}
Environmental modeling is a common application in field robotics, where a robot is used to create a model of the environment, such as a spatial field.
One application is in ocean sampling, where one can create models of sea-surface temperature, or algae distributions.
Informative sampling is a proven technique for spatial field modeling~\cite{Krause2007,Singh2007,Low2008}.
In informative sampling a model is created of the field using, for example, Gaussian Process (GP) regression.
Information-theoretic metrics, such as entropy or mutual information, are then used to decide where the robot should sample next.
When sampling is run on-line, incorporating data sampled by the robot, this is called adaptive sampling.

In this work, we investigate how humans select waypoints, when given the same information as a robot sampling a field, and we compare the modeling performance between the manual and autonomous waypoint selection methods.
We construct a graphical user interface where humans can iteratively set waypoints for a robot, and where they get to see what data is revealed along the way.
The humans are given the instruction to guide the robot to best model the underlying field.
Simulations are run for robotic adaptive sampling, with one random sampling approach and two informative sampling approaches.
Results show that for scenarios that correspond to the robot's on-board model assumptions, the robot performs better than the average human, though similar to the best human.
For scenarios that do not correspond to the modeling assumptions, the robot performs similar to the average human, and random sampling performs on average as well as sampling based on GP models.

\section{Related Work}\label{sec:rel_work}

Borji and Itti~\cite{Borji2013} performed a similar study where they compared human search behavior with various search strategies, including Bayesian Optimization methods and Gaussian Processes.
Their study was more extensive in terms of number of search strategies considered and the number of participants recruited.
However, in their study, the task was simplified to seeking a maximum of a 1D function. As a stopping criteria, participants judged when they had found the maximum.
The study results showed that humans are surprisingly good at this task, and that Bayesian Optimization methods based on Gaussian Processes came closest to the human performance. Furthermore, they concluded that the participants focus on exploring uncertain regions, in active tasks.
This supports our study into comparing adaptive sampling by humans to sampling with robots that use Gaussian Process models.

Chien et al.~\cite{Chien2010} compared manual (human) path planning with automated path planning for search and rescue scenarios.
They used a deterministic roadmap planner, maximizing information gain, and making local decisions.
This automated approach was evaluated against data from human subject experiments. Results showed that for coverage achieved and regions explored, there was no significant difference between human and automated path planning.
The automated path planning did lead to a higher count for victims found, though the manual path planning required less teleoperation.
Overall they concluded that a human path planner could be replaced with an automated path planner without loss in performance.
We are interested to see if their results hold in our environmental modeling scenarios, given the differences in modeling approach and application domain.

Other studies have looked at deducing human strategies and human-like behavior to create planners for robots.
These studies actively try to imitate human or animal behavior.
Henry et al.~\cite{Henry2010} used Gaussian Process regression to estimate the density and flow of crowds of people. They then employed maximum entropy inverse reinforcement learning to learn cost weights for the features, i.e. the density and flow of crowds, to understand how to navigate through crowded spaces.
Sexton and Ren~\cite{Sexton2016} developed a crowdsourcing game to learn human search strategies. They made the assumption that human actions stem from Bayesian Optimization (Gaussian Processes), and used maximum likelihood estimation to learn the model covariance function.

In contrast to these kinds of studies, we are not attempting to learn strategies from our human subject experiment.
Instead, we are comparing performance of the human and adaptive sampling strategies over time, and we gather explicit feedback from the participants about their strategies through a questionnaire.
We wish to better understand how well humans perform at adaptive sampling if given the same information as a robot, what strategies they use, and which approaches are able to learn a good model more quickly.
The results of this project could help identify the potential benefits of using autonomous sampling techniques. This is especially useful for applications where the manual setting of waypoints by a human operator is infeasible or undesirable due to communication constraints or because the human operator has to control too many robots at the same time.

\pagestyle{plain}

\vspace{.2cm}
\section{Gaussian Process Regression}\label{sec:theory}

Gaussian Process (GP) regression is a standard method for spatial field modeling~\cite{Rasmussen2006}.
We briefly recap the theory here:
A GP model is specified by its prior mean and covariance function, also known as the kernel.
A standard approach is to use a zero mean prior and the isotropic squared exponential kernel.
The squared exponential kernel is given by~\cite{Rasmussen2006}:
\begin{equation}
  k({\bf x},{\bf x}' ) = \sigma_{f}^2 \exp \left\{ - \frac{1}{2 l^2} | {\bf x} - {\bf x}' |^2 \right\}
\end{equation}
where ${\bf x}$ and ${\bf x}'$ are locations of a training sample, $\sigma_{f}^2$ is the signal variance (or amplitude), and $l$ is the kernel's length scale.
$\sigma_{f}^2$ and $l$ are the GP's hyperparameters.
The hyperparameters are typically estimated from the data, by using, for example, maximum likelihood estimation~\cite{Rasmussen2006}.

To make predictions based on a GP model, one calculates the predictive mean and variance for test locations:
\begin{equation}
  {\bf y_* ({\bf X_*}, {\bf X}, {\bf y})} = K({\bf X_*}, {\bf X}) K({\bf X}, {\bf X})^{-1}  {\bf y}
\end{equation}
\begin{dmath}
  {\bf \Sigma_*} ({\bf X_*}, {\bf X}) = K({\bf X_*}, {\bf X_*}) - \\ K({\bf X}_*, {\bf X}) K({\bf X}, {\bf X})^{-1} K({\bf X}, {\bf X}_*)
\end{dmath}
where ${\bf X}$ are the training locations, ${\bf y}$ are the training samples, and ${\bf X_*}$ are the test locations.
The entropy for every test location can then be calculated using the standard equation for entropy:
\begin{equation}\label{eq:entropy}
  H = \frac{1}{2} \ln ( 2 \pi \exp \sigma^2 )
\end{equation}

For our robotic adaptive sampling approaches, we test three  different approaches for waypoint selection: one random and two informative sampling approaches.
The first informative approach uses the GP and selects waypoint locations based on maximum entropy.
Note that the entropy does not depend on the measured values, and only depends on the locations of the training samples.
Therefore, we also test another selection method, where we combine the entropy with the predictive mean:
\begin{equation}\label{eq:entropyplusmean}
  {\bf G} = \frac{{\bf H}}{\max({\bf H})} + \alpha \frac{\bm{\mu_{*}}}{\max(\bm{\mu_{*}})}
\end{equation}
where ${\bf G}$ is the new acquisition function, ${\bf H}$ is the entropy for all test locations, $\bm{\mu_{*}}$ is the predictive mean for all test locations, and $\alpha$ is an exploration-exploitation parameter which determines the importance of the predictive mean in the waypoint selection.
We set $\alpha = 0.25$.

\vspace{.2cm}
\section{Experimental Set-Up}\label{sec:setup}

We use the Matlab GPML libraries~\cite{gpml} both for running GP regression on the sampled data, as well as for generating simulation scenarios.

\subsection{Scenarios}\label{sub:scenarios}

Figure~\ref{fig:scenarios} shows the 12 scenarios created for this study.
The first six scenarios are random samples from a Gaussian Process model, whose hyperparameters were set based on averages from prior field data to $l = \exp \{-7.81\}$ and $\sigma_f^2 = \exp \{1.68\}$.
In doing so, we obtain the same smoothness parameters as observed fields.
The latter six scenarios are Gaussian mixture models (GMMs), where the number of Gaussians, their locations and their orientations are all randomly sampled.
Random noise is added up to $5\%$ of the maximum data value.
Data values are ranged from $0$ to $20$ as a proxy for Chlorophyll measurements ($\mu g/L$), which indicates algae abundance.
For all scenarios, data are simulated over a space of $400$x$200\,m$, along a grid with $10\,m$ grid spacing.

\begin{figure}[!b]
  \centering
  \includegraphics[width=\columnwidth]{scenarios.png}
  \caption{Simulated fields: the first six scenarios are randomly sampled from a GP, the latter six are randomly generated Gaussian mixture models. The colorscale ranges from dark blue, low values to yellow, high values.}
  \label{fig:scenarios}
\end{figure}

\subsection{Humans}

For the human subject experiment, we created a graphical user interface (GUI) using python.
The GUI loads an image of the scenario, and creates a black overlay over the scenario.
The user then iteratively selects waypoints, after which the underlying data is revealed between the previously selected and currently selected waypoints.
To simplify data collection and viewing, the black overlay is a grid of squares of the same size as the simulated data grid.
This means that sampled data can easily be indexed from the generated data file.
Furthermore, the revealed squares are large enough for humans to easily observe.

\begin{figure*}[!t]
  \centering
  \begin{minipage}{\textwidth}
    \centering
    \includegraphics[width=.9\textwidth]{blank_screen_screenshot.png}
    \vspace{.1cm}
  \end{minipage}
  \begin{minipage}{\textwidth}
    \centering
    \includegraphics[width=.9\textwidth]{human_path.png}
  \end{minipage}
  \caption{GUI used by human subjects. The scenario is hidden by a black overlay (top),  which is removed as the human iteratively selects waypoints and creates a path for the robot to traverse (bottom).}
  \label{fig:gui}
\end{figure*}

Figure~\ref{fig:gui} shows the created GUI both at the initial state, and after a user has selected waypoints up to the maximum revealed number of squares, i.e. the path length.
The GUI includes a description of the task, such that the instructions are the same for all participants.
It also includes a colorbar to indicate the meaning of the colors, and it shows how many squares can still be revealed.
The total number of squares revealed is limited to 190, which corresponds to 1900-2660 meters of sampling, depending on whether the path is completely horizontal or vertical, or diagonal.
This length was chosen because it is approximately 40\% of the time required to run both a north-south and east-west $20\,m$ grid spaced lawnmower survey over the area.
Without a time limit, every model would end up being equally good.
Though the path length was limited, participants were not restricted in time regarding how long they took to decide where to place the waypoints.
Similarly, computation time for the robot was not considered, though it is on the order of seconds on an Intel i7 quad-core machine with 32GB RAM.

We recruited 11 human participants for this study: 7 men and 4 women, aged between early $20$s and mid $40$s.
Nine participants were computer science and robotics students and professors, one participant works in oceanography and marine biology, and one participant works in student affairs.
None of the participants had run this study before, nor had they seen the scenarios.
All participants ran all scenarios in the order given in Figure~\ref{fig:scenarios}.
Participants were asked to follow the instructions on the screen.
No feedback was given about their performance, throughout the experiment.
Some participants asked questions about what their strategy should be, and were told that they should choose waypoints so as to  collect data that would create the best model of the underlying field.

All participants were also given a survey after running the experiment, which included the following questions:
\begin{itemize}
\item Describe your strategy when choosing waypoints, if you had one.
\item Did the data revealed along the path change the way you chose waypoints?
\end{itemize}

\subsection{Robot}

For the robot experiments,
three approaches were tested: (1) random waypoints, (2) GP model with entropy optimization, (3) GP model optimizing for entropy plus predictive mean.
Equations~\ref{eq:entropy}-\ref{eq:entropyplusmean} showed the calculations for approaches (2) and (3).
Given that there is no sensor or sampling noise, a single run was done for each approach.
The robots created GP models whose kernel is a combination of the isotropic squared exponential covariance function and a white noise kernel, such that there are three hyperparameters: the length scale $l$, signal variance or amplitude $\sigma_f^2$ and signal noise $\sigma_n^2$.
The log of the hyperparameters are initialized to $-7.5, 0.5$ and $1.0$, respectively, based on prior simulations~\cite{Kemna2017}.
The length scale, which is set in longitude/latitude degrees, roughly corresponds to a length of $55\,m$ at a latitude of 34.086.
The hyperparameters are re-estimated when every batch of samples is received and a new waypoint needs to be calculated.

\section{Results}\label{sec:results}

For evaluation, we use two methods of reconstructing the fields using the sampled data.
The first method is a spline-based method, V4, which is one of Matlab's standard interpolation methods, which also allows for extrapolation.
The second method is GP regression, using the GPML Matlab libraries~\cite{gpml}, with a standard isotropic squared-exponential covariance function.
Hyperparameters are estimated on the collected data for every model.
We evaluate performance in terms of the root mean squared error (RMSE) between the ground truth (simulated scenarios) and the learned models.

\begin{figure*}[!t]
\centering
\includegraphics[width=.63\textwidth]{happ_v42.png}
\caption{Robot versus human RMSE boxplots (median, 25th, 75th percentiles) for each field, using a spline-based method (V4) for field interpolation.}
\label{fig:happ_v4}
\end{figure*}

\begin{figure*}[!t]
\centering
\includegraphics[width=.63\textwidth]{happ_gp2.png}
\caption{Robot versus human RMSE boxplots for each field, using GP regression for field interpolation.}
\label{fig:happ_gp}
\end{figure*}

Figures~\ref{fig:happ_v4} and~\ref{fig:happ_gp} show the results of the experiments, with human performance as boxplots (median, 25th, 75th percentiles), and robot performance as points.
For the V4 interpolation, shown in Figure~\ref{fig:happ_v4}, we see that the robot was able to perform better than the average human in most cases, though the best human was able to outperform the robot in all cases, except scenarios 10 and 11.
For the GP interpolation, we see the RMSE greatly declines for the initial 6 scenarios.
This result is logical because these were generated from a GP model.
For the first six scenarios, where the robot optimizes on a GP model, it always performs well when using entropy optimization and approximately as good as the best human.
For the latter six scenarios, the robot performance is similar to the average human performance for all methods.

Table~\ref{tab:avg_rmse} shows the modeling error averaged over all scenarios, for all sampling and modeling methods.
This table shows that models created using GP regression were on average better than models created with V4 interpolation, due to the great improvement for scenarios $1$-$6$ with GP modeling.
Overall, for the robot performance, the best approach was the one that used a GP model with entropy for optimization.

{\setlength{\tabcolsep}{1em}
\bgroup
\def\arraystretch{1.2
\begin{table}[!t]
\centering
\begin{tabular}{|l|l|}
  \hline
  {\bf Method} & {\bf Average RMSE}\\
  \hline
   {\bf GP model} & \\
   Human & 0.048\\
   GP + entropy & {\bf 0.029} \\
   GP + entropy + mean & 0.037\\
   Random & 0.056\\
  \hline
  {\bf V4 model} &\\
   Human & 0.064\\
   GP + entropy & {\bf 0.042}\\
   GP + entropy + mean & 0.048\\
   Random & 0.083\\
  \hline
\end{tabular}
\caption{Average RMSE over all scenarios, for each method of sampling and modeling.}
\label{tab:avg_rmse}
\end{table}
\egroup
}

We further evaluate the modeling performance over time, from the start of each run until all samples are collected.
We evaluate the timed performance for the GP interpolation, given the overall better performance compared to V4 interpolation.
We divide the 190 samples into batches of 20 samples (10 for the last group), and create a GP model at every batch.
We calculate the RMSE between the GP model and the ground truth.
Because of the difference in the generation method for the scenarios, we split this evaluation between fields $1$-$6$ and fields $7$-$12$.
Figures~\ref{fig:rmse_time1} and~\ref{fig:rmse_time2} show the RMSE boxplots (median, 25th and 75th percentiles) over time, for the two groups of fields.

Figure~\ref{fig:rmse_time1} shows that overall the performance of both humans and the robot using the GP with entropy criterion are the most consistent for the six scenarios that were sampled from a GP.
For the initial 20 samples (number of samples axis label $n = 1$), performance is quite similar across all approaches, and none are much better than random.
After 40-60 samples have been gathered ($n = 2$-$3$), the GP with entropy criterion is on average better than the humans at improving the model.
This advantage disappears somewhat over time, though it has a consistently good model across all scenarios towards the end of all simulations.
The random approach is overall the worst in consistency and performance.
The GP with entropy plus mean has more variance in performance across the scenarios, and is less consistent overall, though occasionally better than the other approaches.

\begin{figure*}[!t]
\centering
  \includegraphics[width=.7\textwidth]{rmse_timesteps_fields1-6.png}\\
\caption{RMSE boxplots (median, 25th, 75th percentiles) evaluated over time, 20 samples per time step, for the first six scenarios, averaged over the GP sampled scenarios.}
\label{fig:rmse_time1}
\end{figure*}

\begin{figure*}[!t]
\centering
  \includegraphics[width=.7\textwidth]{rmse_timesteps_fields7-12.png}\\
\caption{RMSE boxplots evaluated over time, 20 samples per time step, for the last six scenarios, averaged over the GMM scenarios. \vspace{-.2cm}}
\label{fig:rmse_time2}
\end{figure*}

Figure~\ref{fig:rmse_time2}, showing RMSE performance over time for the GMM scenarios, shows variable performance across scenarios for all evaluated methods.
Human performance is a lot more consistent for the first $20$-$40$ samples ($n = 1$-$2$), and stays more consistent till the end.
During initial sampling, i.e. $20$-$80$ samples ($n = 1$-$4$), the robotic sampling techniques are on average slightly better at collecting informative samples for the model.
Overall, the difference between the robot approaches is not significant, and the GP-based approaches do not perform better than random sampling initially, though on average they are slightly better after the final $20$-$30$ samples ($160$-$190$, $n = 9$-$10$).

\subsection*{Survey results}
Ten of the eleven participants filled out the survey.
Table~\ref{tab:strategies} lists the most commonly used strategies.
The most commonly employed strategy was to try and cover the whole area. The next most common strategies were to seek out maxima and to sample                                                                                boundary regions rather than inside the extrema.
Two participants are counted for two strategies, because both started every scenario with coverage sampling and then refined the sampling afterward with a different strategy: one in combination with finding maxima and the other in combination with finding boundaries.
Out of the 10 participants who filled out the survey, 8 participants said they used the information they were given, i.e. the points revealed.
The other two participants reported they did not use this information partially or completely. These participants optimized waypoint selection for coverage.

{\setlength{\tabcolsep}{1em}
\bgroup
\def\arraystretch{1.2
\begin{table}[!t]
\centering
\begin{tabular}{|l|l|}
\hline
{\bf Strategy} & {\bf Number of participants}\\
\hline
Coverage & 6/10\\
Find maxima & 3/10\\
Find boundaries & 3/10\\
\hline
\end{tabular}
\caption{Strategies employed by participants.\vspace{-.4cm}}
\label{tab:strategies}
\end{table}
\egroup
}

\section{Discussion \& Future Work}\label{sec:disc}

The results of this study show that when the robot's modeling assumptions match the characteristics of the scenario, the robot performance is better than the average human.
However, in most cases the robot was unable to outperform the best human.
When the scenario does not fit the modeling assumptions well, the robot achieves similar performance to the average human.
The results evaluating RMSE over time confirmed this observation, and also showed that the robot and human performance for the GMM scenarios were not significantly better than random sampling.
This is surprising and interesting, because these scenarios are commonly used as simple algae bloom or underwater plume models.
The results suggest that the use of GP models with entropy-based optimization is not the most effective for these scenarios.
It would be interesting to further explore the GMM scenarios, to see if different modeling and path planning methods could further improve robot performance.

The current study had eleven participants with predominantly computer science and robotics backgrounds.
It would be interesting to run this experiment with a wider range of participants, including subject experts, for example oceanographers, biologists, or people who often process field data, as well as more participants that do not have computer science or robotics backgrounds.
If such a wide pool of participants were recruited, we could also investigate the differences in performance between different backgrounds.

\section{Conclusions}
Adaptive informative sampling approaches are commonly used in robotics for spatial field modeling.
In this paper, we evaluate the modeling performance between an adaptive sampling algorithm and human participants.
We developed a GUI where the human participants could iteratively select waypoints and reveal data, to obtain similar amounts of information that a robot might collect.
Our goal was to better understand how well humans can choose waypoints for adaptive sampling if given the same amount of information, and which strategies they employ.
Results showed that for scenarios that were generated under different assumptions than the model, the GP models with entropy-based criteria and the humans performed no better than random sampling.
When the modeling assumptions corresponded well to the characteristics of the actual field, the robot performance using the GP model and entropy selection criterion, was the most consistent and the best on average.
However, some humans were still able to outperform the robot.
The most commonly used strategy by humans was to optimize for coverage, followed by optimizing for maxima and exploring boundary zones.

\section*{Acknowledgment}

The authors would like to thank everyone who participated in the experiment!

\bibliographystyle{IEEEtran}