Let $z=x+iy \in \mathbb{H}:=\{z= x+ i y\in\mathbb{C}: y>0\}$ and $
\theta (s;z)=\sum_{(m,n)\in\mathbb{Z}^2 } e^{-s \frac{\pi }{y }|mz+n|^2}$ be the theta function associated with the lattice $\Lambda ={\mathbb Z}\oplus z{\mathbb Z}$. In this paper we consider the following pair of minimization problems
\begin{equation}\aligned
\min_{ \mathbb{H} } \theta (2;\frac{z+1}{2})+\rho\theta (1;z),\;\;\rho\in[0,\infty), \\
\min_{ \mathbb{H} } \theta (1; \frac{z+1}{2})+\rho\theta (2; z),\;\;\rho\in[0,\infty),
\endaligned\end{equation}
where the parameter $\rho\in[0,\infty)$ represents the competition of two intertwining lattices, the particular selection of the parameters $s=1,2$ is determined by the physical model, which can be generalized by our strategy and method proposed here.
We find that as $\rho$ varies the optimal lattices admit a novel pattern: they move from rectangular (the ratio of long and short sides changes from $\sqrt3$ to 1 continuously), square, rhombus (the angle changes from $\pi/2$ to $\pi/3$ continuously) to hexagonal continuously; geometrically, up to an invariant group(a subgroup of the classical modular group), they move continuously on a special curve ; furthermore, there exists a closed interval of $\rho$ such that
the optimal lattices is always a square lattice. This is the first, novel and also the complete result on the minimizer problem for theta functions with parameter $\rho$.
This is in sharp contrast to optimal lattice shapes for single theta function ($\rho=\infty$ case), for which the hexagonal lattice prevails.
As a consequence, we give a partial and positive answer to optimal lattice arrangements of vortices in competing systems of Bose-Einstein condensates as conjectured (and numerically and experimentally verified) by Mueller-Ho \cite{Mue2002}, this is the first progress on Mueller-Ho conjecture.