\section{III-V Ternary Alloy Parameter Calculations}\label{app:iii_v_params}  %
Following the standard linear and quadratic interpolation schemes for III-V
ternary alloys \AxBxC{x}{1-x}, with composition parameter $x$ and in terms of experimentally measured values of their binary constituents, AB and BC, we have the lattice constant $a$, energy gap $E$, and effective mass at the $\Gamma$ point $m^{\Gamma\ast}$ as~\cite{Adachi2017}
    \begin{eqnarray}
        a_{\text{\AxBxC{x}{1-x}}}  &= x a\textsubscript{AC} + (1 - x)
        a_{\text{BC}}, \label{eq:a_AxB1xC} \\
        E_{\text{\AxBxC{x}{1-x}}} &= x E_{\text{AC}} + (1 - x)
        E_{\text{BC}} + x (1 - x) E_{\text{AB}},
        \label{eq:E_AxB1xC} \\
        m^{\Gamma\ast}_{\AxBxC{x}{1-x}} &= x
        m^{\Gamma}_{\text{AC}}
        + (1 - x) m^{\Gamma}_{\text{BC}} + x (1 - x)
        m^{\Gamma}_{\text{AB}}. \label{eq:mgamma_AxB1xC}
    \end{eqnarray}
Similarly, the hole effective masses follow from a quadratic interpolation
scheme of the AB, AC binary components as computed from a spherical band
approximation of the valence band edge~\cite{Adachi2017}
    \begin{eqnarray}
      m_{\mathrm{p,dos}} = \left(m_{\mathrm{lh}}^{3/2} + m_{\mathrm{hh}}^{3/2}\right)^{2/3},
      \label{eq:mp_dos} \\
      m_{\mathrm{p, c}} = \frac{m_{\mathrm{lh}}^{5/2} + m_{\mathrm{hh}}^{5/2}}{m_{\mathrm{p, dos}}},
      \label{eq:mp_c} \\
      m_{\mathrm{p, c, \AxBxC{x}{1-x}}} =
      x m_{\mathrm{p, c, AC}}
      + (1-x) m_{\mathrm{p, c, BC}}, \label{eq:mpc_AxB1xC} \\
      m_{\mathrm{p, dos, \AxBxC{x}{1-x}}} =
      x m_{\mathrm{p, dos, AC}}
      + (1-x) m_{\mathrm{p, dos, BC}}. \label{eq:mpdos_AxB1xC}
    \end{eqnarray}
We recognize that the spherical band approximation may not apply to the III-V materials in our study, but it gives an estimate for density of states and conduction band effective masses that are inputs to the COMSOL Semiconductor Module materials models.

To estimate the conduction band offsets between the \InxAlxAs{x}{1-x} and \InxGaxAs{x}{1-x} layers, we followed another interpolation scheme that computes the absolute
conduction band edges $E_{\mathrm{c}}$ using experimentally measured parameters of InAs, AlAs, and GaAs\cite{Krijn1991}
   \begin{eqnarray}
       E_{\mathrm{c}} &= E\textsubscript{v,avg} + \frac{\Delta_0}{3} +  E_{\mathrm{g}} + \Delta E_{\mathrm{c}}^{\mathrm{hy}}
       \label{eq:Ec} \\
       \Delta E_{\mathrm{c}} &= E_{\mathrm{c}}^B - E_{\mathrm{c}}^A, \label{eq:dEc}
   \end{eqnarray}
where $E_{\mathrm{v,avg}}$ is the average valence band edge, $\Delta_0$ is the
spin-orbit splitting in the absence of strain, $E_{\mathrm{g}}$ is the band gap energy, and
$\Delta E_{\mathrm{c}}^{\mathrm{hy}}$ is the shift of the conduction band edge due to
hydrostatic strain. These parameters are calculated from the following
expressions with coefficients $C_{ij}$ read-off from Table 3 compiled by Krijn~\cite{Krijn1991}
    \begin{eqnarray}
        E\textsubscript{v, avg} &= \sum_{i=1}^2
        C_{i0} (E\textsubscript{v,avg})x^i, \label{eq:Evavg} \\
        \Delta_0 &= \sum_{i=1}^2 C_{i0} (\Delta_0) x^i,\label{eq:delta0} \\
        \Delta E_{\mathrm{c}}^{\mathrm{hy}} &=
        \frac{\Delta a(x)}{a(x)} \sum_{i=0}^1
        C_{i0}(\Delta E_{\mathrm{c}}^{\mathrm{hy}}) x^i, \label{eq:dEchy}\\
        \Delta a(x) &= a_0 - a(x). \label{eq:delta_a}
    \end{eqnarray}
\section{Charge--Charge Interaction Matrix Element Derivation}\label{app:charge_matrix}
    \begin{figure}[!ht]
        \centering

        \includegraphics[width=0.6\textwidth]{tunable_cap_2deg_phases_parasitic.eps}
        \caption{Coupler circuit models. (a) The 2DEG coupler compact representation with a single gate and a pair of transmon qubits compared with (b) the simplified circuit used in the derivation of the charge-charge interaction matrix in the main text. (c) Parasitic capacitance circuit model of the coupler and two transmons; capacitances taken from figure~\ref{fig:cmatrix}.}
        \label{fig:circuit_phases}
    \end{figure}
Starting from the two transmon circuit coupled by a voltage-controlled Josephson junction (our 2DEG coupler) in figure~\ref{fig:circuit_phases} (b), with phases $\varphi_1, \varphi_2, \varphi_3$ referring to the left, right, and coupling junctions, respectively, we have~\cite{Orlando1999}
    \begin{eqnarray}
        \varphi_1 - \varphi_2 + \varphi_3 = -\frac{2\pi}{\Phi_0} \Phi_{\mathrm{ext}}, \label{eq:phase_loop_rule}
    \end{eqnarray}
where $\Phi_{\mathrm{ext}}$ is the flux threading the loop formed by the three Josephson junctions as in a typical flux qubit circuit. For a finite $\Phi_{\mathrm{ext}}$, the potential energy $U$ is given by
    \begin{eqnarray}
        U(\mathbf{\varphi}) &= \sum_j E_{\mathrm{J}_j}(1 - \cos\varphi_j) \nonumber \\
        &= E_{\mathrm{J}_1}(1 - \cos\varphi_1) + E_{\mathrm{J}_2}(1 - \cos\varphi_2) \nonumber \\
        &+ E_{\mathrm{J}_3}(1 - \cos(\varphi_2 - \varphi_1 - 2\pi\Phi_{\mathrm{ext}}/\Phi_0)), \label{eq:U_phi}
    \end{eqnarray}
with the signs on the phases following figure~\ref{fig:circuit_phases} (b), preserving the conventions chosen in.~\cite{Orlando1999} We order the phases in a single column vector as
    \begin{eqnarray}
        \mathbf{\varphi} =
        \begin{pmatrix}
        \varphi_1 \\
        \varphi_2
        \end{pmatrix}. \label{eq:def_phi_vector}
    \end{eqnarray}
Setting $\Phi_{\mathrm{ext}}$ = 0, we compute the kinetic energy $T$ by using the Josephson equation relating the voltages at nodes with $k=\left\{1,\ 2\right\}$, $V_k = (\Phi_0/2\pi)\dot{\varphi}_k$ and the definition of $T$ in terms of $\dot{\varphi}_k$
    \begin{eqnarray}
        T &= \frac{1}{2}\left(C_1 V_1^2 + C_2 V_2^2 + C_3 V_3^2\right) \nonumber \\
          &= \frac{1}{2}\left(\frac{\Phi_0}{2\pi}\right)^2\left(C_1\dot{\varphi}^2_1 + C_2\dot{\varphi}^2_2
          + C_3\left(\dot{\varphi}_2 - \dot{\varphi}_1\right)^2\right) \nonumber \\
          &= \frac{1}{2}\left(\frac{\Phi_0}{2\pi}\right)^2 \dot{\mathbf{\varphi}}^T \mathbf{C} \dot{\mathbf{\varphi}}, \label{eq:T_phi_dot}
    \end{eqnarray}
and reading off the capacitance matrix
    \begin{eqnarray}
        \mathbf{C} =
        \begin{pmatrix}
        C_1 + C_3 & -C_3 \\
        -C_3 & C_2 + C_3
        \end{pmatrix}. \label{eq:def_c_matrix_qq}
    \end{eqnarray}
Relating the total capacitances (both the intrinsic junction and external capacitance, commonly referred to as $C_{\Sigma}$~\cite{Schuster2007}) shunting the junctions, $C_1,\ C_2$, to the anharmonicities extracted from the EPR calculations, we have, from the asymptotic expressions derived by Koch~\etal~\cite{Koch2007}
    \begin{eqnarray}
        C_k &= \frac{e^2}{2E_{\mathrm{C}}} \simeq -\frac{e^2}{2\alpha_k} \label{eq:alpha_to_CJ}
    \end{eqnarray}
and we take $C_3$ = $C_{12}(V_{\mathrm{g}})$, the gate voltage-dependent capacitance across the 2DEG coupler. The classical Lagrangian $\mathcal{L}$ and Hamiltonian $\mathcal{H}$ associated with the kinetic and potential energies above, then read~\cite{Orlando1999}
    \begin{eqnarray}
        \mathcal{L(\pmb{\varphi}, \pmb{\dot{\varphi}})} &= T - U \nonumber \\
        &= \frac{1}{2} \left(\frac{\Phi_0}{2\pi}\right)^2 \dot{\pmb{\varphi}}^T \mathbf{C} \dot{\pmb{\varphi}}
        - \sum_j E_{J_j}(1 - \cos\varphi_j) \label{eq:def_classical_lagrangian} \\
        \mathcal{H} &= \mathbf{P}^T \pmb{\dot{\varphi}} - \mathcal{L} \nonumber \\
        &= \frac{1}{2}\mathbf{Q}^T \mathbf{C}^{-1} \mathbf{Q} + U(\pmb{\varphi}) \label{eq:def_classical_hamiltonian} \\
        P_j &= \frac{\partial \mathcal{L}}{\partial \dot{\varphi_j}}
        = \left(\frac{\Phi_0}{2\pi}\right)^2 \sum_kC_{jk}\pmb{\dot{\varphi}}_k, \ \mathbf{Q} = \frac{2\pi}{\Phi_0} \mathbf{P} \label{eq:momenta}
    \end{eqnarray}
We take the form of the quantized Hamiltonian to be the same as the classical one in (\ref{eq:def_classical_hamiltonian}) with classical variables promoted to operators, and identify the charge-charge matrix elements as $e^2[\mathbf{C}^{-1}]_{ij}/2$. Similarly, we write the Lagrangian and identify the capacitance matrix corresponding to the parasitic capacitance model given by the circuit in Fig.~\ref{fig:circuit_phases} (c) as
    \begin{eqnarray}
        \mathcal{L} &= \frac{1}{2}
        \left( \frac{\Phi_0}{2\pi} \right)^2
        \left[
          (C_1 + C_{11}) \dot{\varphi}^2_1
        + (C_2 + C_{22}) \dot{\varphi}^2_2 \right. \nonumber \\
        &+ \left. C_{33} \dot{\varphi}^2_3
        + C_{12}(\dot{\varphi}_2 - \dot{\varphi}_1)^2\right. \nonumber \\
        &+ \left. C_{13}(\dot{\varphi}_3 - \dot{\varphi}_1)^2
         + C_{23}(\dot{\varphi}_2 - \dot{\varphi}_3)^2
        \right] - U(\varphi) \label{eq:lagrangian_parasitic}
    \end{eqnarray}

        \begin{eqnarray}
            \mathbf{C} &=
            \begin{pmatrix}

            \widetilde{C}_{11} & -C_{12} & -C_{13} \\
            -C_{12} & \widetilde{C}_{22} & -C_{23} \\
            -C_{13} & -C_{23} & \widetilde{C}_{33}
            \end{pmatrix}, \label{eq:cmatrix_parasitic}
        \end{eqnarray}

where $\widetilde{C}_{11}=C_1 + C_{11} + C_{12} + C_{13}$, $\widetilde{C}_{22}=C_{2} + C_{22} + C_{12} + C_{23}$, and $\widetilde{C}_{33}=C_{13} + C_{23} + C_{33}$.

\subsection{Description and Realization of the Tuning Element}
We propose a voltage controlled capacitive coupling element between neighboring superconducting qubits. Interaction Hamiltonian realized  is of the form $H_{\mathrm{int}, 12}\propto Q_1 Q_2$, where $Q_i$ refers to the charge operator associated with the $i$-th qubit.

The device is composed of a pair of fixed air-gap capacitors and a single, effective capacitance tuned by a collection of metal-oxide gates that deplete a two dimensional electron gas (2DEG) in an InAs/InGaAs semiconductor heterostructure.\cite{Mayer2019} Air-gap capacitors isolate the tuning elements, namely the gates, from the qubits, minimizing stray electric fields on the neighboring qubits. This is an improvement over existing inductive couplers, where stray magnetic fields can lead to classical cross-talk between qubits.~\cite{Bialczak2011} We suspect that the 2DEG coupler will introduce more charge noise than the inductive couplers, though transmons qubits are first order insensitive to charge noise.~\cite{Koch2007}

\begin{figure}[!ht]
    \centering
    \includesvg[width=0.5\columnwidth]{tunable_cap_2deg.svg}
    \caption{Schematic of two qubits and the 2DEG coupler.}
    \label{fig:capacitor_qubit_schematic}
\end{figure}
We recognize the recent work by~\citeauthor{Casparis2019} that uses a similar coupling mechanism in their superconducting field effect transistor design to realize a tunable frequency superconducting resonator bus.\cite{Casparis2019} Compared with our implementation, the tunable bus design has a higher circuit complexity and lower on/off ratio, although we expect the experimental on/off ratio to differ from the reported value in figure~\ref{fig:cmatrix}.
\begin{figure}[!ht]
    \centering
    \includesvg[width=0.7\columnwidth]{simplified_three_gate_device_stack_scaled.svg}
    \caption{Schematic of the device as modeled in COMSOL Multiphysics. Regions labeled $n^+$ refer to highly doped $n$-type InAs, $\epsilon_0$ refers to air-gapped regions, and dark gray regions are source (drain) terminals. Gold electrodes are used for the gate terminals. We have scaled the vertical dimensions a factor of four relative to the horizontal dimensions to enhance the 2DEG region.}
    \label{fig:device_stack}
\end{figure}

\subsection{Two Qubit Coupler}
In figure~\ref{fig:two_qubit_coupler} we have a microwave circuit model of two transmon qubits coupled by a lumped impedance $Z_{\mathrm{JJ}_{\mathrm{c}}}(\omega) = 1/\left(1/R + i\omega C\right)$, where $R$ and $C$ take the values of $R_{12}$ and $C_{12}$ in either the fully conducting or fully depleted limits of the 2DEG coupler. In the conducting limit, where some current can flow across the coupler and act like a Josephson junction, one might consider adding an inductance to the coupler lumped element model. Taking $R_{12}$ to be the normal resistance of a Josephson junction and computing the junction inductance with the Ambegaokar-Baratoff formula~\cite{Ambegaokar1963}, the junction inductances would be very small, resulting in high coupler mode frequencies, outside of the frequency band of the finite element electromagnetic field solver, Ansys HFSS. For this reason and expected small modifications to qubit-qubit interactions, we omit these inductances in our model and use Ansys HFSS to compute the lowest electromagnetic eigenmodes of the device with the two transmon qubits, indexed by $j$, defined as parallel $LC$ lumped elements, $Z_{\mathrm{JJ}_{\mathrm{q}},j} = 1/\left(1/(i\omega L_{\mathrm{q}, j}) + i\omega C_{\mathrm{q},j}\right).$ In the following section, we use these eigenmode solutions to estimate the Hamiltonian matrix elements corresponding to qubit-qubit mode and qubit-coupler mode coupling strengths. We will differentiate between this modal coupling and direct capacitive coupling in the final part of this section, where we calculate the direct charge-charge interaction matrix elements.
\begin{figure}[!ht]
    \centering

    \includegraphics[width=0.6\textwidth]{qubit_coupler_mask_with_insets.eps}
    \caption{False color geometry of two transmon qubits with the capacitive coupler in between used in the HFSS simulations. Lumped impedances defined in the gold regions of the insets, represent the linear response of the Josephson junctions and capacitive coupling element in the HFSS model.}
    \label{fig:two_qubit_coupler}
\end{figure}
\subsection{Energy Participation Ratios and Quantization}
To extract the coupling matrix elements between the qubits in our microwave device layout, we employ the energy participation ratio (EPR) method developed by Minev~\cite{Minev2019}. This approach goes beyond the larger family of black box quantization methods~\cite{Nigg2012, Solgun2014}, where the Hamiltonian describes a collection of Josephson junction-based qubits interacting with any number of harmonic modes separates into linear and nonlinear terms.

One can relate the modal decomposition of the classical electromagnetic response, e.g. impedance, admittance, or electromagnetic energies, with the linear and nonlinear parts of the Hamiltonian. Additional inputs describing the Josephson junction energy scales, $E_{\mathrm{J}}$ and $E_{\mathrm{C}}$, related to the inductive and capacitive energies of the junction, complete the model. The total Hamiltonian, accounting for $M$ modes, in the dispersive regime and under the rotating wave approximation, reads
    \begin{eqnarray}
        H &= H_{\mathrm{lin}} + H_{\mathrm{nl}} \label{eq:H_disp} \\
        H_{\mathrm{lin}} / \hbar &= \sum_{m=1}^M\omega_m a_m^{\dagger}a_m \label{eq:Hlin} \\
        H_{\mathrm{nl}} / \hbar &= -\sum_{m=1}^M\left(\Delta_ma_m^{\dagger}a_m  + \frac{1}{2}\alpha_m
        a_m^{\dagger\ 2}a_m^2 \right) \nonumber \\
        &+ \frac{1}{2}\sum_{m\neq n} \chi_{mn} a_m^{\dagger} a_m a_n^{\dagger}a_n, \label{eq:Hnl}
    \end{eqnarray}
where the Lamb shifts $\Delta_m$, cross-Kerr coefficients $\chi_{mn}$, and anharmonicities $\alpha_m$ are given by~\cite{Minev2019}
    \begin{eqnarray}
        \Delta_m &= \frac{1}{2}\sum_{n=1}^{M}\chi_{mn} \label{eq:lamb_shift} \\
        \chi_{mn} &= -\sum_{j\in J}\frac{1}{2} \frac{\hbar\omega_m\omega_n}{4E_{\mathrm{J}_j}} \label{eq:chi_mn} \\
        \alpha_m &= \frac{1}{2} \chi_{mm}. \label{eq:alpha_m}
    \end{eqnarray}

Table~\ref{tab:epr_chi_w_alpha} gives a summary of the cross- and self-Kerr (anharmonicities) coefficients along with the three lowest electromagnetic modes of the device as extracted with the pyEPR Python package~\cite{Minev2019, Minev2020}. The eigenfrequencies and quality factors in the first two columns of table~\ref{tab:epr_chi_w_alpha} follow from the HFSS eigenmode solutions.
    \begin{table}[!ht]
        \centering
        \begin{tabular}{p{0.2\linewidth}p{0.2\linewidth}%
                        p{0.12\linewidth}p{0.12\linewidth}%
                        p{0.12\linewidth}}
          \hline\hline
          Qubit Index & $\omega/2\pi$ [GHz] & Q & \multicolumn{2}{c}{$\chi/2\pi$ [MHz]} \\
          \hline

          1 (d) &  6.0228 & 1.7E7  & 129   & 1.02  \\
          2 (d) &  8.6135 & 5.1E9  & 1.02  & 0.002 \\
          \hline
          1 (c) &  6.0228 & 4.5E8  & 223  & 67.1  \\
          2 (c) &  8.6135 & 1.3E9  & 67.1  & 223  \\

          \hline\hline
        \end{tabular}
        \caption{Extracted matrix elements from energy participation ratio calculations for values of $R_{12}$, $C_{12}$ in the fully depleted (d) and fully conducting (c) limits of the 2DEG. Diagonal entries of $\chi$ are scaled by 1/2 to denote the anharmonicities from (\ref{eq:alpha_m}).}
        \label{tab:epr_chi_w_alpha}
    \end{table}
    \begin{figure}[!ht]
        \centering

        \includegraphics[width=0.6\textwidth]{hfss_two_qubit_coupler_modes_220528.eps}
        \caption{Electric field magnitude (dB scale to enhance color contrast) for the first two eigenmode solutions computed with HFSS. (a) 6.0228 GHz and (b) 8.6135 GHz qubits in the fully depleted limit of the coupler.}

        \label{fig:hfss_eigenmodes}
    \end{figure}
\subsection{Extraction of the Exchange Interaction}
To compute the charge-charge interaction strength between the transmon qubits in our HFSS model, we consider the capacitance matrix associated with a persistent current or flux qubit following the derivation by Orlando~\cite{Orlando1999}. For details on the derivation of the capacitance matrix, see \ref{app:charge_matrix}. The Hamiltonian for the coupled two transmons, written in terms of the Josephson junction phases $\varphi_j$ and node charges $q_j$, is given by
    \begin{eqnarray}
        H &= \frac{1}{2}\mathbf{Q}^T\mathbf{C}^{-1}\mathbf{Q} + U(\mathbf{\varphi}) \label{eq:def_H_charge_phase}\\
        U(\mathbf{\varphi}) &= \sum_jE_{\mathrm{J}_j}\left(1 - \cos\varphi_j\right) \label{eq:def_U_phase} \\
        \mathbf{C} &=
        \begin{pmatrix}
        C_1 + C_3 & -C_3 \\
        -C_3 & C_2 + C_3
        \end{pmatrix}. \label{eq:def_c_matrix_qq_main}
    \end{eqnarray}

In (\ref{eq:def_H_charge_phase}), the charge-charge matrix elements are one half the entries of the inverse of the capacitance matrix. We numerically inverted $\mathbf{C}$ in (\ref{eq:def_c_matrix_qq_main}) using values for $C_1$, $C_2$ obtained from (\ref{eq:alpha_to_CJ}) and $C_3=C_{12}(V_{\mathrm{g}})$ in the depleting and conducting limits. Table~\ref{tab:exchange_matrix_elements} summarizes these results, recovering a on/off of more than one order of magnitude, $r_{\mathrm{int}}\approx$14.
\begin{table}[!ht]
    \centering
    \begin{tabular}{p{0.3\linewidth}p{0.25\linewidth}p{0.25\linewidth}}
    \hline\hline
    Qubit Index & \multicolumn{2}{c}{$e^2[\mathbf{C}^{-1}]_{ij}$ [MHz]} \\
    \hline
    1 (d)   &  \phantom{-}198  & -2.22 \\
    2 (d)   &  -2.22 & \phantom{-}198 \\
    \hline
    1 (c)   & \phantom{-}161 & -31.6  \\
    2 (c)   & -31.6 & \phantom{-}161   \\
    \hline\hline
    \end{tabular}
    \caption{Simplified capacitance model charge-charge interaction matrix elements in the (d) depleting and (c) conducting limits of the 2DEG coupler. The ratio of the matrix elements in the conducting and depleted limits gives an on/off interaction ratio of $r_{\mathrm{int}}\approx 14$.}
    \label{tab:exchange_matrix_elements}
\end{table}
\begin{table}[!ht]
    \centering
    \begin{tabular}{p{0.3\linewidth}p{0.2\linewidth}p{0.2\linewidth}p{0.2\linewidth}}
    \hline\hline
    Node Index & \multicolumn{3}{c}{$e^2[\mathbf{C}^{-1}]_{ij}$ [MHz]} \\
    \hline
    1 (d)   &  \phantom{-}195  & -2.30 & -49.2   \\
    2 (d)   &  -2.30 & \phantom{-}195  & -49.2   \\
    3 (d)   &  -49.2 & -49.2 & \phantom{-}1.44E4 \\
    \hline
    1 (c)   & \phantom{-}125 & -24.7  & -37.5   \\
    2 (c)   & -24.7 & \phantom{-}125  & -37.5   \\
    3 (c)   & -37.5 & -37.5 & \phantom{-}1.17E4 \\
    \hline\hline
    \end{tabular}
    \caption{Parasitic capacitance charge-charge interaction matrix elements in the (d) depleting and (c) conducting limits of the 2DEG coupler. The ratio of the matrix elements in the conducting and depleted limits gives an on/off interaction ratio of $r_{\mathrm{int}}\approx 11$.}
    \label{tab:exchange_matrix_elements_parasitic}
\end{table}

We emphasize here that the off-diagonal charge-charge interaction matrix elements give a more accurate description of the coupling between the qubits mediated by the 2DEG coupler than the EPR calculations of the cross-Kerr coefficients.

A more detailed treatment of the coupler including the parasitic capacitances from the Maxwell capacitance matrices in figure~\ref{fig:cmatrix} gives a three-by-three coupling capacitance matrix in the Lagrangian of the form in (\ref{eq:cmatrix_parasitic}). Table~\ref{tab:exchange_matrix_elements_parasitic} summarizes the modified charge-charge interaction matrix elements in the full parasitic capacitance model.
\subsection{Estimation of Coupler Coherence Limit}
To estimate the total coherence limit of our coupler in the two qubit device in figure~\ref{fig:hfss_eigenmodes}, we found that a back-of-the-envelope calculation of the energy stored in the coupler surface, in either qubit mode, is on the order of 10\textsuperscript{-3}, resulting in a coherence limit of a few to tens of $\mu$s, when considering the loss to be dominated by the gate dielectric and the top InGaAs layer. This is consistent with previous studies of transmon qubits whose participations are near unity in their given mode~\cite{Nigg2012,Minev2020} and gives us further confidence that a coupler of a similar geometry could support fast parametric operations with moderate coherence.

\section{Introduction}
\input{intro.tex}

\section{Methods and Modeling}
\subsection{Conceptual Design}\label{sec:conceptual_design}
\input{conceptual_design.tex}

\subsection{Classical Modeling}\label{sec:coupler_modeling}
\input{modeling.tex}

\section{Results}
\input{results.tex}

\section{Discussion}\label{sec:discussion}
\input{discussion.tex}

\section{Conclusion}\label{sec:conclusion}
\input{conclusion.tex}

\section*{Acknowledgments}
We acknowledge funding from NSF Grant PHY-1653820 and Army Research Office Grant Nos. W911NF-18-1-0125 and W911NF-18-1-0067. N. M. acknowledges funding from the Graduate Fellowships for STEM Diversity and W. M. S. acknowledges funding from the Army Research Office QuaCGR Fellowship. We would like to thank Bradley Llyod, Chien Liu, Paul Niyonkuru, Alan Phillips, John Rose, Meenakshi Singh, and David Young for many insightful discussions.

\subsubsection{Semiconductor 2DEG Calculations}
To estimate the capacitance of the 2DEG coupler, we compute the electron concentrations in the active region of the device (InGaAs/InAs/InGaAs layers) using the COMSOL Multiphysics Semiconductor Module~\cite{SemiconductorModuleUsersGuide}. Equilibrium solutions to the drift-diffusion equations with Fermi-Dirac statistics serve to identify regions of high depletion under the gate(s) when applying negative voltages on the order of a few volts, overcoming the work function of the aluminum gate contact.

We use a layer structure typical of gatemon qubits as in figure~\ref{fig:device_stack} and refer to this structure as the ``device stack''~\cite{Shabani2016, Wickramasinghe2018, Strickland2020}. To model the device stack in COMSOL, we specified the following electronic properties of the semiconductor materials and the dielectric constant of the gate oxide: electron and hole effective conduction band masses $m^{\ast}_{\mathrm{n(p),c}}$, low-field mobilities $\mu^{\mathrm{lf}}_{\mathrm{n(p)}}$, band gap energies $E\textsubscript{g}$, conduction band offsets $\Delta E\textsubscript{c}$ between neighboring semiconductors, dielectric constants $\varepsilon\textsubscript{r}$, and effective densities of states for the conduction and valence bands $N\textsubscript{c(v)}$. Taking the electron affinity $\chi$ for InAs as given by the COMSOL material library,  we calculated the remaining affinities using Anderson's affinity rule and the conduction band offsets of each material~\cite{Streetman2015Ch5}. Table~\ref{tab:mat_table} gives a summary of the material parameters used in these semiconductor simulations; see \ref{app:iii_v_params} for detailed calculations of the energy gaps, effective masses, and conduction band offsets for \InxGaxAs{x}{1-x} and \InxAlxAs{x}{1-x} as functions of the composition parameter (x).
\begin{figure}[!ht]
    \centering

    \includegraphics[width=0.5\textwidth]{full_stack_iiiv_220521.eps}
    \caption{Schematic of the 2DEG coupler as modeled in COMSOL Multiphysics based on~\cite{Mayer2019}. An aluminum contact deposited on~\AlxOy{x}{y}~defines the gate terminal. We abbreviate the fixed composition ternary III-V alloys In\textsubscript{0.81}Ga\textsubscript{0.19}As and In\textsubscript{0.81}Al\textsubscript{0.19}As, as InGaAs and InAlAs, respectively. Not shown or modeled is the superlattice graded buffer layer between InP and InAlAs~\cite{Shabani2016, Wickramasinghe2018, Strickland2020}.}
    \label{fig:device_stack}
\end{figure}
\begin{table}[!ht]
        \centering
        \begin{tabular}{p{0.25\linewidth}p{0.18\linewidth}p{0.18\linewidth}p{0.18\linewidth}p{0.18\linewidth}}
        \hline
        \                &  InAs  & InGaAs & InAlAs & InP\\
        \hline\
        $E\textsubscript{g}$ [eV]                      &  0.354           & 0.473            & 0.752  & 1.344 \\
        $\Delta E\textsubscript{c}$ [eV]               &  -               & 0.200            & 0.201  & 0.12\\
        $\varepsilon_r$                 & 15.15            & 14.03            & 13.13  & 12.9\\
        $N\textsubscript{c}$ [cm\textsuperscript{-3}]             & 6.6E18 & 1.4E17  & 2.1E17 & \\
        $N\textsubscript{v}$ [cm\textsuperscript{-3}]             & 8.73E16 & 6.4E18 & 7.8E18 & \\
        $\chi$ [eV]                     & 4.9              & 4.7              & 4.5 & 4.38 \\
        $\mu^{\mathrm{lf}}_{\mathrm{n(p)}}$ [cm\textsuperscript{2}V\textsuperscript{-1}s\textsuperscript{-1}]& 14.4E3 (500)   & 14.4E3 (450)     & 14.E3 (384) & 5.4E3 (200)\\
        $m^{\ast}_{\mathrm{n(p),c}}$ [$m_0$]     & 0.023 (1.00)     & 0.03 (0.25)      & 0.04 (0.31) & 0.08 (0.60)\\
        \hline\hline
        \end{tabular}
        \caption{Materials parameters used in the COMSOL Semiconductor Module calculations. InGaAs and InAlAs abbreviate In\textsubscript{0.81}Ga\textsubscript{0.19}As and In\textsubscript{0.81}Al\textsubscript{0.19}As. $m_0$ corresponds to the rest mass of an electron (0.511 MeV $c^{-2}$). Out-of-plane effective electron and hole masses of InAs are set to $m_0$ in the model to simulate 2DEG confinement in the xy-plane. Electron mobilities for InAs, InGaAs, and InAlAs are all set to the same value as extracted from measurements of a similar device at millikelvin temperatures~\cite{Wickramasinghe2018}. Values not in parenthesis (in parenthesis) correspond to electron (hole) properties.}
        \label{tab:mat_table}
    \end{table}

We specify the geometry in figure~\ref{fig:device_stack} using the native COMSOL CAD editor to define domains (surfaces or planes) and boundaries (lines or edges), solving for the electron density in the domains and on the boundaries. Electronic properties assigned to each domain follow from table~\ref{tab:mat_table}. We model the terminals (source -- 1, drain -- 2, gate -- 3, as in figure~\ref{fig:capacitor_qubit_schematic}) as Terminal boundary conditions with voltages $V_1$, $V_2$, $V_3$ and contact work functions $\Phi\textsubscript{c,1}$, $\Phi\textsubscript{c,2}$, $\Phi\textsubscript{c,3}$ = 4 V~\cite{SemiconductorModuleUsersGuide}.

We selected the density gradient discretization scheme~\cite{Ancona2011} in COMSOL to approximate the quantum confinement effects in the 2DEG more efficiently than a self-consistent Schr\"{o}dinger-Poisson equation calculation. The density gradients modify the equilibrium electron ($n$) and hole ($p$) concentrations by~\cite{SemiconductorModuleUsersGuide}
    \begin{eqnarray}
        n = N\textsubscript{c} F_{1/2}\left(\frac{E\textsubscript{fn} - E\textsubscript{c} + q V^{\mathrm{DG}}_{\mathrm{n}}}{k\textsubscript{B}T}\right) \label{eq:n_dg} \\
        p = N\textsubscript{v} F_{1/2}\left(\frac{E\textsubscript{v} - E\textsubscript{fp} + q V_{\mathrm{p}}^{\mathrm{DG}}}{k\textsubscript{B}T}\right) \label{eq:p_dg} \\
        N_{\mathrm{c(v)}} = \left(\frac{2 m_{\mathrm{n(p)}}^{\ast} \pi k_{\mathrm{B}} T}{h^2}\right)^{3/2}, \label{eq:N_cv}
    \end{eqnarray}
where $E_{\mathrm{c(v)}}$ is a given material's conduction (valence) band edge, $E_{\mathrm{fn(p)}}$ are the electron (hole) quasi-Fermi level energies, $F_{1/2}(\eta)$ is the Fermi-Dirac integral~\cite{Kim2019}, $k_{\mathrm{B}}$ is Boltzmann's constant, $T$ is the temperature of the device (approximate temperature of the mixing chamber stage of typical dilution refrigerators $\sim$10 mK), and $q$ is the charge of an electron or hole. The quantum potentials $V_{\mathrm{n(p)}}^{\mathrm{DG}}$ are defined in terms of the density gradients by~\cite{SemiconductorModuleUsersGuide}
    \begin{eqnarray}
        \nabla\cdot \left(\mathbf{b}_{\mathrm{n}} \nabla \sqrt{n}\right) &= \frac{1}{2}\sqrt{n}\ V_{\mathrm{n}}^{\mathrm{DG}} \label{eq:grad_n} \\
        \nabla\cdot \left(\mathbf{b}_{\mathrm{p}} \nabla \sqrt{p}\right) &= \frac{1}{2}\sqrt{p}\ V_{\mathrm{p}}^{\mathrm{DG}}, \label{eq:grad_p}
    \end{eqnarray}
with the density gradient tensors $\mathbf{b}_{\mathrm{n(p)}}$ for electrons (holes) expressed in terms of the effective mass tensors $\mathbf{m}_{\mathrm{n(p)}}^{\ast}$
    \begin{eqnarray}
        \mathbf{b}_{\mathrm{n}} =
        \frac{\hbar^2}{12q}\left[\mathbf{m}_{\mathrm{n}} \right]^{-1} \label{eq:b_n} \\
        \mathbf{b}_{\mathrm{p}} =
        \frac{\hbar^2}{12q}\left[\mathbf{m}_{\mathrm{p}}\right]^{-1}. \label{eq:b_p}
    \end{eqnarray}
Note the distinction between the scalar effective masses $m_{\mathrm{n(p)}}^{\ast}$ and, the effective mass tensors $\mathbf{m}_{\mathrm{n(p)}}$. Anisotropy in the effective mass tensors emulates the quantum confinement effects in the 2DEG, constraining electron movement to one plane.

For the remaining materials, InAlAs, \AlxOy{2}{3}, and air, we used the Electric Charge Conservation interface, including the following constitutive relations for each dielectric in terms of its electric permittivity tensor $\pmb{\varepsilon}$~\cite{SemiconductorModuleUsersGuide}
    \begin{eqnarray}
       \mathbf{D} = \varepsilon_0\pmb{\varepsilon}:\mathbf{E}, \label{eq:D_epsE}
    \end{eqnarray}
where $\mathbf{D}$ is the electric displacement field, $\varepsilon_0$ is the permittivity of free space, $\mathbf{E}$ is the electric field, and $\pmb{\varepsilon}:\mathbf{E}$ is a tensor contraction (matrix vector product) between $\pmb{\varepsilon}$ and $\mathbf{E}$. Modeling these regions as pure dielectrics reduces the size of the system of equations relative to a drift-diffusion calculation applied to the materials that behave as perfect insulators, e.g. large dielectric layers (InP and InAlAs) and oxide layers.
    \begin{figure}[!ht]
        \centering

        \includegraphics[width=0.5\textwidth]{n_gated_0Vg_n4Vg_one_220527.eps}
        \caption{Electron concentrations [cm\textsuperscript{-3}] on a base-10 logarithmic scale with source-drain bias $V_{\mathrm{sd}}$ = +5 mV - (-5 mV) = 10 mV for (a) the fully conducting $V_{\mathrm{g}}$ = 0 V and (b) fully depleted $V_{\mathrm{g}}$ = -4 V operating points. We do not solve for $n$ in the regions where we applied the Electric Charge Conservation equations, i.e. in the \AlxOy{x}{y}, InAlAs, and InP regions, yet the electric fields respect the boundary conditions set by those layers. The vertical axis is scaled by a factor of four relative to the horizontal to improve the visibility of the active region.}
        \label{fig:n_conducting_depleted}
    \end{figure}
\subsubsection{Electrostatic Capacitance Matrix Calculations}
Following the semiconductor calculations described in the previous section, we developed a simplified model to extract the capacitance and conductance matrices. This model replaces regions with high electron concentration with perfect electric conductors and regions with low electron concentration with perfect dielectrics. The COMSOL Electrostatics Interface defines the Maxwell Capacitance Matrix $\mathbf{C}$ in terms of the charges and voltages on each node of an $N$-terminal network as~\cite{ACDCModuleUsersGuide}
    \begin{eqnarray}
        \begin{pmatrix}
        Q_1 \\
        Q_2 \\
        \vdots \\
        Q_N
        \end{pmatrix}
        =
        \begin{pmatrix}
        C_{11} & C_{12} & \hdots & C_{1N} \\
        C_{21} & C_{22} & \hdots & C_{2N} \\
        \vdots & \vdots & \      & \vdots \\
        C_{N1} & C_{N2} & \hdots & C_{NN}
        \end{pmatrix}
        \begin{pmatrix}
        V_1 \\
        V_2 \\
        \vdots \\
        V_N
        \end{pmatrix}, \label{eq:def_maxwell_cmatrix}
    \end{eqnarray}
where $Q_k$ and $V_k$ are the charges and voltages on the $k$-th terminal. The charges and voltages follow from solutions to Poisson's equation, and the charge continuity equation. In two dimensions these equations combine to give~\cite{ACDCModuleUsersGuide}
    \begin{eqnarray}
       -\nabla\cdot d \left(\varepsilon_0 \nabla V - \mathbf{P}\right) = \rho, \label{eq:electrostatics_equation}
    \end{eqnarray}
where $V$ is the electric potential, $\mathbf{P}$ is the polarization vector, $\rho$ is the space charge density, and $d$ is the out-of-plane thickness. This model gives an approximation of the electron densities from the two-dimensional semiconductor interface calculations in the previous section, with $d$ = 5 $\mu$m.

We set the boundaries of domains representing perfect electric conductors as Floating Potential boundary conditions and define terminals as ohmic contacts on the corresponding dielectrics in figure~\ref{fig:device_stack}~\cite{ACDCModuleUsersGuide}.

In the fully conducting and depleted limits ($V_{\mathrm{g}}$ = 0 V and $V_{\mathrm{g}}$ = -4 V), we compute the capacitance matrices in figures~\ref{fig:cmatrix} (a) and (b). The matrix elements of interest, $C_{12}=C_{21}$ represent the effective capacitance between the source and drain terminals. These terminals form capacitive contacts with any pair of qubits. The tuning ratio $r$, or on/off contrast of the 2DEG coupler is given by $r$ = $C_{12, \mathrm{conducting}}$ / $C_{12, \mathrm{depleted}}$ $\approx$ 29 from the electrostatic solutions in figures~\ref{fig:cmatrix} (a) and (c).
\begin{figure}[!ht]
    \centering

    \includegraphics[width=0.6\textwidth]{csd_one_gate_220527.eps}
    \caption{Capacitance matrix computed with the COMSOL Electrostatics Interface for (a) the fully conducting and (c) the fully depleted limit of the 2DEG. Capacitance matrix extracted from the admittance matrix computed with the COMSOL Electric Currents Interface for (c) the fully conducting and (d) the fully depleted limit of the 2DEG. Vertical and horizontal axes labels correspond to terminal indices in figure~\ref{fig:capacitor_qubit_schematic}. The tuning ratio between the conducting and depleted 2DEG is given by $r = C_{12,\mathrm{conducting}} / C_{12,\mathrm{depleted}}>20$.}
    \label{fig:cmatrix}
\end{figure}
\subsubsection{Electric Currents Admittance Matrix Calculations}
To extract the conductance matrix and verify the capacitance matrix of the device under high frequency excitation agrees with the electrostatic result, we use the AC/DC Electric Currents Interface to compute the admittance matrix $\mathbf{Y}$ defined in terms of the $N$ terminal voltages $V_k$ and currents $I_k$~\cite{ACDCModuleUsersGuide}
    \begin{eqnarray}
        \begin{pmatrix}
        I_1 \\
        I_2 \\
        \vdots \\
        I_N
        \end{pmatrix}.
        =
        \begin{pmatrix}
        Y_{11} & Y_{12} & \hdots & Y_{1N} \\
        Y_{21} & Y_{22} & \hdots & Y_{2N} \\
        \vdots & \vdots & \      & \vdots \\
        Y_{N1} & Y_{N2} & \hdots & Y_{NN}
        \end{pmatrix}
        \begin{pmatrix}
        V_1 \\
        V_2 \\
        \vdots \\
        V_N
        \end{pmatrix}
        \label{eq:def_ymatrix}
    \end{eqnarray}
In the frequency domain, the voltages and currents become phasors of the form $\widetilde{V}_ke^{i\omega t}$ and $\widetilde{I}_ke^{i\omega t}$, with the admittance matrix given by
    \begin{eqnarray}
        \mathbf{Y} = \mathbf{G} + i\omega \mathbf{C}, \label{eq:ymatrix_freq_domain}
    \end{eqnarray}
where $\mathbf{G}$ and $\mathbf{C}$ are the conductance and capacitance matrices, $i=\sqrt{-1}$, and $\omega$ is the angular frequency~\cite{ACDCModuleUsersGuide}. Both matrices are symmetric for our linear, three-terminal device in figure~\ref{fig:cmatrix}.

The Electric Currents Interface in COMSOL computes the currents and voltages, in (\ref{eq:def_ymatrix}) by solving the time harmonic equations~\cite{ACDCModuleUsersGuide}
    \begin{eqnarray}
        \nabla\cdot\left(\pmb{\sigma}:\mathbf{E} + \mathbf{J}_e\right) + i\omega\rho = 0\label{eq:continuity_eq}, \\
       \nabla \cdot \mathbf{D} = \rho, \label{eq:divD}
    \end{eqnarray}
where $\mathbf{E}$ is the electric field, $\mathbf{J}_e$ is the externally generated electric current density, $\pmb{\sigma}$ is the electrical conductivity tensor, and $\rho$ is the electric space charge density. We compute the currents at each terminal, given voltage source excitations, as a function of frequency $\omega$ in the band of 1-10 GHz relevant to superconducting qubit and resonator frequencies, and extracted the conductance matrix as the real, frequency-independent part of $\mathbf{Y}$ and the capacitance matrix as the derivative of the imaginary part of $\mathbf{Y}$ from (\ref{eq:ymatrix_freq_domain})
    \begin{eqnarray}
        G_{ij} = \mathrm{Re}\{Y_{ij}(\omega)\} \label{eq:G_ij} \\
        C_{ij} = \mathrm{Im}\left\{\frac{\rmd Y_{ij}(\omega)}{\rmd\omega}\right\}, \label{eq:C_ij}
    \end{eqnarray}
where $\rmd Y_{ij}(\omega) / \rmd\omega$ is a constant in our case, as we omit the junction inductance $L_{\mathrm{J}_0}$ leading to discontinuities at resonance frequencies proportional to $(C_{ij}L_{\mathrm{J}_0})^{-1/2}$~\cite{Nigg2012}.

\subsection{Coupler Loss Estimates}
We give bounds on the losses introduced by the 2DEG coupler from experimental measurements of the high participating gate dielectrics and InGaAs upper layer, along with the other layers in the device stack. In table~\ref{tab:pratios}, we compute the electric field participation ratios $p_j$ following the procedure developed by~\cite{Wenner2011, Wang2015}. The relaxation time $T_1$ at a given angular frequency $\omega$, as a function of the dielectric material properties and geometric factors, reads~\cite{Wang2015}
    \begin{eqnarray}
        T_1^{-1} &= \frac{\omega}{Q} = \omega\sum_{j}\frac{p_j}{Q_j} + \Gamma_0 \label{eq:T1_participation} \\
        Q_j^{-1} &= \tan\delta_j \label{eq:Qj_tan_delta_j} \\
        p_j &= W_{\mathrm{e}}^{-1}t_{\mathrm{oxide}}\varepsilon_{1,j}
        \int_{S_j}\left|\mathbf{E}\right|^2 dS \label{eq:pj} \\
        W_{\mathrm{e}} &= \int_{V} |\mathbf{E}|^2 dV \label{eq:We}
    \end{eqnarray}
where $W_e$ is the electric field energy density stored in the volume of the entire geometry $V$, $Q_j$ are the quality factors, $\tan\delta_j$ are the loss tangents, $\varepsilon_{1,j}$ are the real parts of the dielectric function, and $t_{\mathrm{oxide}}$ is the thickness of the participating lossy surface, assumed to be 3 nm for all materials~\cite{Wang2015}. The participation ratios give the fraction of the electrical energy stored in a given surface $S_j$ relative to the total electrical energy stored in the entire volume  of the device. The last term in (\ref{eq:T1_participation}), $\Gamma_0$, includes all other loss mechanisms contributing to $T_1$ besides dielectric loss~\cite{Wang2015}. Note, these participation ratios differ from those in subsequent calculations involving \textit{energy} participation ratios referenced to a given mode rather than a particular surface.

Other sources of loss relevant to III-V semiconductor materials, but not considered in this study, include piezoelectricity~\cite{Scigliuzzo2020}, non-equilibrium quasiparticles~\cite{Nguyen2022}, cosmic ray muon flux~\cite{Wilen2021}, and, to a lesser extent, stray magnetic fields~\cite{Kringhoj2021}.
\begin{table}[ht]
        \centering
        \begin{tabular}{lllll}
            \hline\hline
            Depleted & $t_j$ [nm] & $p_{j,\mathrm{norm}}$ & $\tan\delta_j^{\ast}$ & $T_1$ [$\mu$s] \\
            \hline
            InGaAs (Top)               & 10    & 4.19E-2 & 4.1E-4   & 1.85E+0  \\
            InAs                       & 4     & 1.03E-2 & 4.1E-4   & 7.53E+0  \\
            InGaAs (Bottom)            & 4     & 7.85E-3 & 4.1E-4   & 9.89E+0  \\
            InAlAs                     & 100   & 2.92E-2 & 4.1E-4   & 2.67E+0  \\
            \AlxOy{2}{3}~\cite{McRae2020b} & 50    & 9.04E-1 & 5E-3 & 6.87E-3 \\
            InP                        & 3.5E3 & 6.97E-3 & 4.1E-4   & 8.95E+0  \\
            Total                      & -     & 1       & 7.3E-3   & 6.81E-3  \\
            \hline
            Conducting & $t_j$ [nm] & $p_{j,\mathrm{norm}}$ & $\tan\delta_j^{\ast}$ & $T_1$ [$\mu$s] \\
            \hline
            InGaAs (Top)                & 10    & 1.01E-8 & 4.1E-4     & 7.69E+6 \\
            InAs                        & 4     & 3.73E-9 & 4.1E-4     & 2.08E+7 \\
            InGaAs (Bottom)             & 4     & 4.03E-9 & 4.1E-4     & 1.93E+7 \\
            InAlAs                      & 100   & 1.10E-9 & 4.1E-4     & 7.06E+7 \\
            \AlxOy{2}{3}~\cite{McRae2020b}  & 50    & 9.9999E-1 & 5E-3 & 6.24E-3 \\
            InP                         & 3.5E3 & 1.41E-5 & 4.1E-4     & 4.43E+3 \\
            Total                       & -     & 1       & 7.3E-3     & 6.24E-1 \\
            \hline\hline
        \end{tabular}
        \caption{Participation ratios $p_j$, dielectric loss tangents $\tan\delta_j$, layer thicknesses $t_j$, and estimated dielectric-loss-limited $T_{1,j}$. All $T_{1,j}$ times are referenced to a qubit frequency of $\omega/2\pi$ = 5 GHz and $\tan\delta_j^{\ast}$ indicates that in the absence of reliable loss tangent data for the individual InAs, InGaAs, InAlAs, and InP layers, we used the low power loss extracted from measurements of an Al patterned CPW resonator on the full III-V stack modeled in this work and measured at 100 mK as an estimate~\cite{Phan2022}.}
        \label{tab:pratios}
\end{table}
\subsection{Integration with Circuit QED}\label{sec:int_cqed}
\subsubsection{Two Qubit Coupler}
In figure~\ref{fig:two_qubit_coupler} we have a microwave circuit model of two transmon qubits coupled by a lumped impedance $Z_{\mathrm{JJ}_{\mathrm{c}}}(\omega) = 1/\left(1/R + i\omega C\right)$, where $R$ and $C$ take the values of $R_{12}$ and $C_{12}$ in either the fully conducting or fully depleted limits of the 2DEG coupler. In the conducting limit, where some current can flow across the coupler and act like a Josephson junction, one might consider adding an inductance to the coupler lumped element model. Taking $R_{12}$ to be the normal resistance of a Josephson junction and computing the junction inductance with the Ambegaokar-Baratoff formula~\cite{Ambegaokar1963}, the junction inductances would be very small, on the order of a few aH to tens of fH, resulting in high coupler mode frequencies, far outside of the frequency band of the finite element electromagnetic field solver, Ansys HFSS. For this reason and expected small modifications to qubit-qubit interactions, we omit these inductances in our model and use Ansys HFSS to compute the lowest electromagnetic eigenmodes of the device with the two transmon qubits, indexed by $j$, defined as parallel $LC$ lumped elements, $Z_{\mathrm{JJ}_{\mathrm{q}},j} = 1/\left(1/(i\omega L_{\mathrm{q}, j}) + i\omega C_{\mathrm{q},j}\right).$ In the following section, we use these eigenmode solutions to estimate the Hamiltonian matrix elements corresponding to qubit-qubit mode and qubit-coupler mode coupling strengths. We will differentiate between this modal coupling from direct capacitive coupling in the final part of this section, where we calculate the direct charge-charge interaction matrix elements.
\begin{figure}[!ht]
    \centering
    \includegraphics[width=0.6\textwidth]{qubit_coupler_mask_with_insets.eps}
    \caption{False color geometry of two transmon qubits with the capacitive coupler in between used in the HFSS simulations. Lumped impedances defined in the gold regions of the insets, represent the linear response of the Josephson junctions and capacitive coupling element in the HFSS model.}
    \label{fig:two_qubit_coupler}
\end{figure}
\subsubsection{Energy Participation Ratios and Quantization}
To extract the coupling matrix elements between the qubits in our microwave device layout, we employ the energy participation ratio (EPR) method developed by Minev~\cite{Minev2019}. This approach goes beyond the larger family of black box quantization methods~\cite{Nigg2012, Solgun2014}, where the Hamiltonian describes a collection of Josephson junction-based qubits interacting with any number of harmonic modes separates into linear and nonlinear terms.

One can relate the modal decomposition of the classical electromagnetic response, e.g. impedance, admittance, or electromagnetic energies, with the linear parts of the Hamiltonian. Additional inputs describing the Josephson junction energy scales, $E_{\mathrm{J}}$ and $E_{\mathrm{C}}$, related to the inductive and capacitive energies of the junction, account for the nonlinear terms. The total Hamiltonian, accounting for $M$ modes, in the dispersive regime and under the rotating wave approximation, reads
    \begin{eqnarray}
        H &= H_{\mathrm{lin}} + H_{\mathrm{nl}} \label{eq:H_disp} \\
        H_{\mathrm{lin}} / \hbar &= \sum_{m=1}^M\omega_m a_m^{\dagger}a_m \label{eq:Hlin} \\
        H_{\mathrm{nl}} / \hbar &= -\sum_{m=1}^M\left(\Delta_ma_m^{\dagger}a_m  + \frac{1}{2}\alpha_m
        a_m^{\dagger\ 2}a_m^2 \right) \nonumber \\
        &+ \frac{1}{2}\sum_{m\neq n} \chi_{mn} a_m^{\dagger} a_m a_n^{\dagger}a_n, \label{eq:Hnl}
    \end{eqnarray}
where the Lamb shifts $\Delta_m$, cross-Kerr coefficients $\chi_{mn}$, and anharmonicities $\alpha_m$ are given by~\cite{Minev2019}
    \begin{eqnarray}
        \Delta_m &= \frac{1}{2}\sum_{n=1}^{M}\chi_{mn} \label{eq:lamb_shift} \\
        \chi_{mn} &= -\sum_{j\in J}\frac{1}{2} \frac{\hbar\omega_m\omega_n}{4E_{\mathrm{J}_j}} \label{eq:chi_mn} \\
        \alpha_m &= \frac{1}{2} \chi_{mm}. \label{eq:alpha_m}
    \end{eqnarray}

Table~\ref{tab:epr_chi_w_alpha} gives a summary of the cross- and self-Kerr (anharmonicities) coefficients along with the three lowest electromagnetic modes of the device as extracted with the pyEPR Python package~\cite{Minev2019, Minev2020}. The eigenfrequencies and quality factors in the first two columns of table~\ref{tab:epr_chi_w_alpha} follow from the HFSS eigenmode solutions.
    \begin{table}[!ht]
        \centering
        \begin{tabular}{p{0.2\linewidth}p{0.2\linewidth}%
                        p{0.12\linewidth}p{0.12\linewidth}%
                        p{0.12\linewidth}}
          \hline\hline
          Qubit Index & $\omega/2\pi$ [GHz] & Q & \multicolumn{2}{c}{$\chi/2\pi$ [MHz]} \\
          \hline

          1 (d) &  6.0228 & 1.7E7  & 129   & 1.02  \\
          2 (d) &  8.6135 & 5.1E9  & 1.02  & 0.002 \\
          \hline
          1 (c) &  6.0228 & 4.5E8  & 223  & 67.1  \\
          2 (c) &  8.6135 & 1.3E9  & 67.1  & 223  \\

          \hline\hline
        \end{tabular}
        \caption{Extracted matrix elements from energy participation ratio calculations for values of $R_{12}$, $C_{12}$ in the fully depleted (d) and fully conducting (c) limits of the 2DEG. Diagonal entries of $\chi$ are scaled by 1/2 to denote the anharmonicities from (\ref{eq:alpha_m}).}
        \label{tab:epr_chi_w_alpha}
    \end{table}
    \begin{figure}[!ht]
        \centering

        \includegraphics[width=0.5\textwidth]{hfss_two_qubit_coupler_modes_220528.eps}
        \caption{Electric field magnitude (dB scale to enhance color contrast) for the first two eigenmode solutions computed with HFSS. (a) 6.0228 GHz and (b) 8.6135 GHz qubits in the fully depleted limit of the coupler.}

        \label{fig:hfss_eigenmodes}
    \end{figure}
\subsubsection{Extraction of the Exchange Interaction}
To compute the charge-charge interaction strength between the transmon qubits in our HFSS model, we consider the capacitance matrix associated with a persistent current or flux qubit following the derivation by Orlando~\cite{Orlando1999}. For details on the derivation of the capacitance matrix, see \ref{app:charge_matrix}. The Hamiltonian for the coupled two transmons, written in terms of the Josephson junction phases $\varphi_j$ and node charges $q_j$, is given by
    \begin{eqnarray}
        H &= \frac{1}{2}\mathbf{Q}^T\mathbf{C}^{-1}\mathbf{Q} + U(\mathbf{\varphi}) \label{eq:def_H_charge_phase}\\
        U(\mathbf{\varphi}) &= \sum_jE_{\mathrm{J}_j}\left(1 - \cos\varphi_j\right) \label{eq:def_U_phase} \\
        \mathbf{C} &=
        \begin{pmatrix}
        C_1 + C_3 & -C_3 \\
        -C_3 & C_2 + C_3
        \end{pmatrix}. \label{eq:def_c_matrix_qq_main}
    \end{eqnarray}

In (\ref{eq:def_H_charge_phase}), the charge-charge matrix elements are one half the entries of the inverse of the capacitance matrix. We numerically inverted $\mathbf{C}$ in (\ref{eq:def_c_matrix_qq_main}) using values for $C_1$, $C_2$ obtained from (\ref{eq:alpha_to_CJ}) and $C_3=C_{12}(V_{\mathrm{g}})$ in the depleting and conducting limits. Table~\ref{tab:exchange_matrix_elements} summarizes these results, recovering an on/off interation ratio of more than one order of magnitude, $r_{\mathrm{int}}\approx$14.
\begin{table}[!ht]
    \centering
    \begin{tabular}{p{0.3\linewidth}p{0.25\linewidth}p{0.25\linewidth}}
    \hline\hline
    Qubit Index & \multicolumn{2}{c}{$e^2[\mathbf{C}^{-1}]_{ij}$ [MHz]} \\
    \hline
    1 (d)   &  \phantom{-}198  & -2.22 \\
    2 (d)   &  -2.22 & \phantom{-}198 \\
    \hline
    1 (c)   & \phantom{-}161 & -31.6  \\
    2 (c)   & -31.6 & \phantom{-}161   \\
    \hline\hline
    \end{tabular}
    \caption{Simplified capacitance model charge-charge interaction matrix elements in the (d) depleting and (c) conducting limits of the 2DEG coupler.}
    \label{tab:exchange_matrix_elements}
\end{table}
\begin{table}[!ht]
    \centering
    \begin{tabular}{p{0.3\linewidth}p{0.2\linewidth}p{0.2\linewidth}p{0.2\linewidth}}
    \hline\hline
    Node Index & \multicolumn{3}{c}{$e^2[\mathbf{C}^{-1}]_{ij}$ [MHz]} \\
    \hline
    1 (d)   &  \phantom{-}195  & -2.30 & -49.2   \\
    2 (d)   &  -2.30 & \phantom{-}195  & -49.2   \\
    3 (d)   &  -49.2 & -49.2 & \phantom{-}1.44E4 \\
    \hline
    1 (c)   & \phantom{-}125 & -24.7  & -37.5   \\
    2 (c)   & -24.7 & \phantom{-}125  & -37.5   \\
    3 (c)   & -37.5 & -37.5 & \phantom{-}1.17E4 \\
    \hline\hline
    \end{tabular}
    \caption{Parasitic capacitance charge-charge interaction matrix elements in the (d) depleting and (c) conducting limits of the 2DEG coupler. }
    \label{tab:exchange_matrix_elements_parasitic}
\end{table}

We emphasize here that the off-diagonal charge-charge interaction matrix elements give a more accurate description of the coupling between the qubits mediated by the 2DEG coupler than the EPR calculations of the cross-Kerr coefficients.

A more detailed treatment of the coupler including the parasitic capacitances from the Maxwell capacitance matrices in figure~\ref{fig:cmatrix} gives a three-by-three coupling capacitance matrix in the Lagrangian of the form in (\ref{eq:cmatrix_parasitic}). Table~\ref{tab:exchange_matrix_elements_parasitic} summarizes the modified charge-charge interaction matrix elements in the full parasitic capacitance model. In the parasitic capacitance model, we find an on/off interaction ratio of $r_{\mathrm{int}}\approx 11$.
\subsection{Estimation of Coupler Coherence Limit}
To estimate the total coherence limit of our coupler in the two qubit device in figure~\ref{fig:hfss_eigenmodes}, a back-of-the-envelope calculation of the energy stored in the coupler surface, in either qubit mode, gives an electric field participation ratio on the order of 10\textsuperscript{-3}, resulting in a coherence limit of a few to tens of $\mu$s, when considering the loss to be dominated by the gate dielectric and the top InGaAs layer. This is consistent with previous studies of transmon qubits whose participations are near unity in their given mode~\cite{Nigg2012,Minev2020} and gives us further confidence that a coupler of a similar geometry could support fast parametric operations with moderate coherence.


\section{Introduction}\label{Sec 1}
When an acoustic field encounters inhomogeneity, it exerts acoustic radiation force on it. Here by inhomogeneity, we mean non-uniform or discontinuous variation of physical properties in a system such as particles/cells suspended in fluid, emulsions, co-flowing streams of miscible or immiscible fluids, and fluid subjected to a temperature gradient. The acoustic forces acting on inhomogeneity are extensively studied in microscale flows, and this field is known as 'microscale acoustofluidics'\citep{Friend2011Jun}. Over the last two decades, acoustofluidics has found a wide range of applications in biological\citep{Ahmed2016Mar,Iranmanesh2015Sep,Collins2015Nov,Christakou2013Apr,Lakshmanan2020Sep}, chemical\citep{Shi2009Dec,Xie2020Dec}, and medical\citep{Li2015Apr,Lu2019Feb,Zhang2020Mar} sciences.

Recently, the relocation and stabilization of inhomogeneous co-flowing fluid streams in microchannels has gained the attention of the research community which is evident from the following works. Through silicon-glass microchannel experiments, \citet{Deshmukh2014Jul} could relocate high-impedance sodium chloride solution to node (center) and low-impedance water to anti-node (sides). Also, they could stabilize high-impedance fluid at the center (and low-impedance fluid to the sides) against gravity stratification using acoustic forces. In addition to the above experiments on miscible fluids, \citet{Hemachandran2019Apr} demonstrated the relocation of immiscible fluids using acoustic fields by overcoming the interfacial tension forces. Followed by this, \citet{Karlsen2018Jan} showed that acoustic forces acting on stable inhomogeneous fluid configuration could effectively suppress the boundary-driven Rayleigh streaming in the bulk. The theoretical framework and understanding of the above non-linear acoustic forces on inhomogeneous fluids are provided by \citet{Rajendran2022Jun,Karlsen2016Sep,Karlsen2018Jan}. Other notable works on the practical applications of acoustic forces on co-flowing inhomogeneous fluids include iso-acoustic focusing of cells \citep{augustsson2016iso}, acoustic focusing of sub-micron particles \citep{VanAssche2020Feb,Gautam2018Oct}, tweezing and patterning of inhomogeneous fluids in a microchannel \citep{Karlsen2017Mar,Baudoin2020Aug}, rapid mixing of fluids using an alternating multinode method \citep{Pothuri2019Dec}, reversible stream-droplet transition in a microfluidic co-flowing immiscible system \citep{Hemachandran2021Sep}. Despite the above recent advancements and practical importance, the criterion at which the inhomogeneous fluid system becomes unstable or stable under acoustic fields has not been clearly established. This paper aims to establish the stability criterion of inhomogeneous co-flowing fluids subjected to standing acoustic wave fields.

In this work, using linear stability analysis, we derive the dispersion relation that governs the stability of inhomogeneous fluids (with and without interfacial tension) under acoustic body force. We study the various parameters such as the initial arrangement of fluids, the position of the interface with respect to the node, acoustic energy density, the height of the channel, and surface tension to establish the necessary and sufficient conditions for relocation and stability. For fluids with interfacial tension, a non-dimensional number called acoustic Bond number is obtained theoretically which characterizes stable and unstable (relocation) regime. Also, we deduce a relation between critical acoustic energy density and the height of the channel which paves a way for relocating fluids with higher interfacial tension $(\mathcal{O}(10^1$ to $10^2) mN/m)$ in a microchannel. Furthermore, numerical simulations are carried out using generalized acoustic body force which agrees well with the derived theoretical stability criterion.

\section{Physics of the problem}\label{Sec 2}
The hydrodynamics of the inhomogeneous fluids involved in this study is governed by the mass-continuity and momentum equations \citep{Landau1987Aug},
\begin{subequations}
\label{eq 1}
    \begin{gather}
    \label{eq 1a}
        \frac{\partial\rho}{\partial t}+\boldsymbol{\nabla\cdot}(\rho \textbf{\emph{V}})=0, \\
    \label{eq 1b}
        \rho \frac{D\textbf{\emph{V}}}{Dt}= -\boldsymbol{\nabla} P +\eta\nabla^2 \textbf{\emph{V}} +\beta\eta\boldsymbol{\nabla}(\boldsymbol{\nabla\cdot}\textbf{\emph{V}}) + \textbf{\emph{f}}_{ac}.
    \end{gather}
\end{subequations}
where $\rho$ represents density, $\textbf{\emph{V}}$ represents the velocity vector field, $P$ represents the pressure field, $\eta$ is the dynamic viscosity of the fluid,  $\beta=(\xi/\eta)+(1/3)$, $\xi$ is the bulk viscosity, and $D/Dt$ denotes the material derivative ($D/Dt=\partial_t+\textbf{V}\boldsymbol{\cdot\nabla}$). Here the body force $\textbf{\emph{f}}_{ac}$ is only due to acoustics. The gravitational body force is neglected since it is dominated by acoustic force in microscale flows. The acoustic body force $\textbf{\emph{f}}_{ac}$ is given as \citep{Rajendran2022Jun}
\begin{align}
\label{eq 2}
     \textbf{\emph{f}}_{ac} = -\boldsymbol{\nabla\cdot}\langle\rho_0\textbf{\emph{v}}_1 \textbf{\emph{v}}_1\rangle = & \left( \frac{1}{2}\boldsymbol{\nabla}\left(\kappa_0\langle|p_1|^2\rangle - \rho_0 \langle|\textbf{\emph{v}}_1|^2\right\rangle)\right) + \biggl[\langle\textbf{\emph{v}}_1\times\boldsymbol{\nabla}\times(\rho_0\textbf{\emph{v}}_1)\rangle\biggr] \nonumber \\
      & + \left(-\frac{1}{2}\langle|p_1|^2\rangle\boldsymbol{\nabla}\kappa_0-\frac{1}{2}\langle\textbf{\emph{v}}_1^2\rangle\boldsymbol{\nabla}\rho_0\right)  \\
      \label{eq 3}
            =& \ (\textbf{\emph{f}}_{1})+[\textbf{\emph{f}}_{2}]+(\textbf{\emph{f}}_{3}).
\end{align}
where $p_1$ and $\textbf{\emph{v}}_1$ denote the first-order (fast time scale) pressure and velocity fields due to acoustic waves (see Appendix \ref{Sec A1}) and $\langle...\rangle$ is the time average in one oscillation period (the time average of two first-order fields $\langle\textbf{\emph{u}}_1\textbf{\emph{v}}_1\rangle$ is defined as $\frac{1}{2}$Real($\textbf{\emph{u}}_1^{\star}\textbf{\emph{v}}_1$), where $\star$ denotes complex conjugation). The terms $\rho_0$ and $\kappa_0$ denote the zeroth-order (background) density and compressibility of the fluid.
In (\ref{eq 2}), the first term is a conservative or gradient term that induces pressure and not fluid flow, the second term is only dominant at boundary layers and is responsible for boundary-driven Rayleigh streaming and the third term is responsible for relocation and stabilization of inhomogeneous fluids. Hence, only the relevant third term ($\textbf{\emph{f}}_{3}$ in (\ref{eq 3})) is considered for theoretical analysis. For the standing acoustic wave applied along the X-direction, $\langle|p_1|^2\rangle=p_a^2\sin^2(k_wx)$, $\langle|\textbf{\emph{v}}_1^2|\rangle=p_a^2/(\rho_0^2c_0^2)\cos^2(k_wx)$, $k_w = 2\upi/\lambda_w$ denotes the wavenumber ($\lambda_w$ denotes the wavelength) and $p_a$ denotes the pressure amplitude. Then the relocation force $f_{rl}$ can be well approximated in terms of impedance gradient as \citep{Rajendran2022Jun}
\begin{equation}
\label{eq 4}
    \textbf{\emph{f}}_{3}=\textbf{\emph{f}}_{rl}=-E_{ac}\cos({2k_wx})\boldsymbol{\nabla}\hat Z.
\end{equation}
where $E_{ac}=p_a^2/(4\rho_{avg}c_{avg}^2)=(v_a^2\rho_{avg})/4$ is the acoustic energy density, $ Z = \rho_0 c_0$ denotes impedance, $c_0$ denotes background (zeroth-order) speed of sound in a medium, $\hat Z=Z/z_{avg}$, $\hat c_0 = c_0/ c_{avg}$ and $\hat \rho_0 = \rho_0/ \rho_{avg}$, where the subscript 'avg' denotes the respective average quantities of fluid A and B.
\begin{figure}
    \center
    \includegraphics[width=1\linewidth]{figure1.png}        \caption{\label{figure 1} Inhomogeneous fluids (of different impedance $Z_A$ and $Z_B$) separated by a plane interface subjected to a standing acoustic half-wave. In the absence of interfacial tension, the fluid system \textbf{(a)} is in an unstable equilibrium, \textbf{(b)} is in a stable equilibrium, \textbf{(c)} is in neutral equilibrium and \textbf{(d)} is in non-equilibrium state. In the presence of interfacial tension, the fluid system \textbf{(a)} is in a conditionally stable equilibrium, \textbf{(b)} is in stable equilibrium, \textbf{(c)} is in stable equilibrium and \textbf{(d)} is in a conditionally stable equilibrium. Note: \textbf{(d)} is not in the scope of this work.}
\end{figure}

\subsection{Stability analysis of inhomogeneous fluids in the absence of interfacial tension\label{Sec 2.1}}
A two-dimensional fluid domain subjected to a standing acoustic half-wave in the X-direction, with two fluids separated by a sharp vertical interface as shown in figure \ref{figure 1}.a-c is considered for the stability analysis. Before beginning the analysis, it is necessary to understand the equilibrium of the system in the absence of interfacial tension. In a completely enclosed domain, a fluid initially at rest ($\textbf{\emph{V}}=0$) will remain at rest (or equilibrium) if the body force can be completely absorbed in pressure, $\textbf{\emph{f}}_{rl}=\boldsymbol{\nabla}P$ from (\ref{eq 1b}) and (\ref{eq 4}). By taking the curl of the above relation, the condition for equilibrium is given as $\boldsymbol{\nabla} \times \textbf{\emph{f}}_{rl} = 0$. Thus,
\begin{equation}
\label{eq 5}
   -\frac{E_{ac}}{z_{avg}}\left[\frac{\partial}{\partial x}\left( \cos(2k_wx)\frac{\partial z}{\partial y}\right)-\frac{\partial}{\partial y}\left(\cos(2k_wx)\frac{\partial z}{\partial x}\right)\right]  =  0.
\end{equation}
It is clear from (\ref{eq 5}) that the given fluid configuration will be in an equilibrium state, only if the direction of the acoustic standing wave is normal to the fluid-fluid interface (the direction of the acoustic standing wave is parallel to the direction of the impedance gradient) as shown in figure \ref{figure 1}($a$-$c$). Since $\boldsymbol{\nabla} \times \textbf{\emph{f}}_{rl} \neq 0$ for the configuration shown in figure \ref{figure 1}$(d)$, it is not in equilibrium and tends to relocate to the stable configuration without any perturbations. The stability nature of the equilibrium configurations is analysed by imposing infinitesimal perturbations on the interface.
Now we proceed to show that in the absence of interfacial tension, the configuration shown in figure \ref{figure 1}$(a)$ is in unstable equilibrium (perturbations grow), figure \ref{figure 1}$(b)$ is in stable equilibrium (perturbations decay), and figure \ref{figure 1}($c$) is in neutral equilibrium (perturbations neither grow nor decay).

The effect of viscosity is neglected in the stability analysis, as it governs only the timescale of the phenomenon and does not contribute to the stability criterion. Although the physical properties are non-uniform in an inhomogeneous system, the fluid particles considered in the flow field have constant density $\rho$, speed of sound $c$, and impedance $Z$. Thus, the material derivative of all properties is zero, which includes the incompressibility condition $(D\rho/Dt=\partial \rho/\partial t+\textbf{\emph{V}}\boldsymbol{\cdot\nabla}\rho = 0)$. By combining the incompressibility condition with (\ref{eq 1a}) and neglecting the viscosity, the governing equations (\ref{eq 1}) reduce to
\begin{subequations}
\label{eq 6}
\begin{gather}
    \label{eq 6a}
    \frac{\partial U}{\partial x}+\frac{\partial V}{\partial y}=0, \\
    \label{eq 6b}
    \rho\frac{DU}{Dt}=-\frac{\partial P}{\partial x}-\frac{E_{ac}\cos({2k_wx})}{z_{avg}}\frac{\partial Z}{\partial x}, \\
    \label{eq 6c}
    \rho\frac{DV}{Dt}=-\frac{\partial P}{\partial y}-\frac{E_{ac}\cos({2k_wx})}{z_{avg}}\frac{\partial Z}{\partial y},
\end{gather}
where $U, V$ are the X-component and Y-component of the velocity field $\textbf{\emph{V}}$. Since the body force term is a function of impedance, the below impedance relation is required for the closure.
\begin{equation}
    \label{eq 6d}
    \frac{DZ}{Dt}=\frac{\partial Z}{\partial t}+U\frac{\partial Z}{\partial x}+V\frac{\partial Z}{\partial y} = 0.
\end{equation}
\end{subequations}
Now, the flow fields are decomposed into an unperturbed zeroth-order stationary state and infinitesimal perturbations as $U=u_0+\delta u,\ V=v_0+\delta v,\ P=p_0+\delta p,\ \rho =\rho_0+\delta\rho$ and $Z=z_0+\delta z$. In this study, the variation of acoustic impedance is considered only in the X-direction (figure \ref{figure 1}($a$-$c$)), $z_0=z_0(x)$. At the stationary state $(u_0=v_0=0)$, the unperturbed zeroth-order equations become, $\frac{\partial p_0}{\partial x} = -\frac{E_{ac}\cos({2k_wx})}{z_{avg}}\frac{\partial z_0}{\partial x}$ from \ref{eq 6b}, $\frac{\partial p_0}{\partial x}=0$ from  \ref{eq 6c} and $\frac{\partial z_0}{\partial t}=0$ from \ref{eq 6d}. Using the above zeroth order relations and  neglecting the second-order terms in (\ref{eq 6}), the first-order perturbation equations governing the stability becomes
\begin{subequations}
\label{eq 7}
    \begin{gather}
    \label{eq 7a}
        \frac{\partial \delta u}{\partial x}+\frac{\partial \delta v}{\partial y}=0,\\
    \label{eq 7b}
        \rho_0\frac{\partial \delta u}{\partial t}= -\frac{\partial \delta p}{\partial x} - \frac{E_{ac}\cos({2k_wx})}{z_{avg}}\frac{\partial \delta z}{\partial x},\\
    \label{eq 7c}
        \rho_0\frac{\partial \delta v}{\partial t} = -\frac{\partial \delta p}{\partial y} - \frac{E_{ac}\cos({2k_wx})}{z_{avg}}\frac{\partial \delta z}{\partial y}, \\
    \label{eq 7d}
        \frac{\partial(\delta z)}{\partial t}=-\delta u\frac{\partial z_0}{\partial x}.
    \end{gather}
\end{subequations}
Analysing the disturbances into normal modes, the amplitude of the disturbances $\delta u,\ \delta v,\ \delta \rho,\ \delta p,$ and $\delta z$ takes the following form
\begin{equation}
    \label{eq 8}
    A(x,y,t)=A(x)exp(ik_yy+nt),
\end{equation}
where $k_y$ is the wavenumber considered along the Y-direction. Applying the above amplitude relations in the form (\ref{eq 8}) in (\ref{eq 7}),
\begin{subequations}
\label{eq 9}
    \begin{gather}
    \label{eq 9a}
        \frac{\partial\delta u}{\partial x}+ik_y\delta v=0, \\
    \label{eq 9b}
        \rho_0 n\delta u=-\frac{\partial \delta p}{\partial x}-\frac{E_{ac}\cos{\left(2k_wx\right)}}{z_{avg}}\frac{\partial \delta z}{\partial x}, \\
    \label{eq 9c}
        \rho_0 n\delta v=-ik_y\delta p -ik_y\frac{E_{ac}\cos{\left(2k_wx\right)}}{z_{avg}}\delta z,\\
\label{eq 9d}
    n\delta z=-\delta u\frac{\partial z_0}{\partial x}.
\end{gather}
\end{subequations}
The partial notation is dropped since the only derivatives in (\ref{eq 9}) are with respect to the $x$ coordinate. Multiplying by $ik_y$ throughout (\ref{eq 9c}) and combining with (\ref{eq 9a}) and (\ref{eq 9d}), we obtain,
\begin{equation}
\label{eq 10}
\delta p = -\rho_0 \frac{n}{{k_y}^2}\frac{d\delta u}{dx}+E_{ac}\frac{\cos(2k_wx)}{z_{avg}}\frac{\delta u}{n}\frac{dz_0}{dx}.
\end{equation}
substituting (\ref{eq 9d}) and (\ref{eq 10}) in (\ref{eq 9b}) results in,
\begin{equation}
\label{eq 11}
     \frac{d}{dx}\left(\rho_0 \frac{d\delta u}{dx} \right) - \rho_0 k_y^2 \delta u = -E_{ac}\frac{2k_w\delta u}{z_{avg}}\frac{k_y^2}{n^2}\frac{dz_0}{dx} \sin(2k_wx).
\end{equation}
Considering two uniform fluids of different impedance $Z_A$ and $Z_B$ separated by interfaces positioned at  $x_s$,
\begin{subequations}
\label{eq 12}
    \begin{equation}
    \label{eq 12a}
    z_0=z_A+(z_B-z_A)H(x-x_s),
    \end{equation}
    \begin{equation}
    \label{eq 12b}
    \frac{dz_0}{dx}=(z_B-z_A)\delta (x-x_s),
    \end{equation}
\end{subequations}
where $H(x-x_s)$ is the Heaviside step function at $x=x_s$ and  $\delta(x-x_s)$ is the Dirac’s $\delta$-function at $x=x_s$. Substituting (\ref{eq 12b}) in (\ref{eq 11}),
\begin{equation}
\label{eq 13}
    \frac{d}{dx}\left(\rho_0 \frac{d\delta u}{dx}\right) -\rho_0k_y^2\delta u= -E_{ac}\frac{2k_w\delta u}{z_{avg}}\frac{k_y^2}{n^2}\sin(2k_wx)(z_B-z_A)\delta(x-x_s).
\end{equation}
Equation (\ref{eq 13}) is the governing differential equation for the stability of inhomogeneous fluids (without interfacial tension). For a uniform region on either side of the interface(s) where there are no discontinuities in the impedance, the governing equation (\ref{eq 13}) reduces to
\begin{equation}
    \label{eq 14}
    \frac{d^2\delta u}{dx^2}-k_y^2\delta u=0.
\end{equation}
The solution of (\ref{eq 14}) is of the form $\delta u=C_1e^{k_y(x-x_s)}+C_2e^{-k_y(x-x_s)}$ where $C_1,C_2$ are constants. Since $\delta u$ must vanish at the boundaries, we can write the solution as,
\begin{subequations}
\label{eq 15}
\begin{equation}
    \label{eq 15a}
    \delta u_B = Ce^{k_y(x-x_s)}\quad\quad(x<x_s),
\end{equation}
\begin{equation}
    \label{eq 15b}
    \delta u_A = Ce^{-k_y(x-x_s)}\quad\quad(x>x_s),
\end{equation}
\end{subequations}
where the constant $C$ in (\ref{eq 15}) is chosen to ensure continuity in velocity across the interfaces. For the solution at the interface $(x=x_s)$, we integrate (\ref{eq 13}) along infinitesimal distance $(dx\approx 0)$, the second term in the left-hand side of the equation is zero and the remaining terms are,
\begin{equation}
\label{eq 16}
\Delta\left(\rho_0 \frac{d\delta u_s}{dx}\right)= -E_{ac}\frac{2k_w\delta u_s}{z_{avg}}\frac{k_y^2}{n^2}(z_B-z_A)\int{ \left(\sin(2k_wx)\delta(x-x_s)\right)dx},
\end{equation}
where $\delta u_s$ is the value of $\delta u$ at $x=x_s$. Using,  (\ref{eq 15}) and the Dirac delta identity $\int f(x)\delta(x-a)dx = f(a)$ to solve for eigenvalue $n$ in (\ref{eq 16}).
\begin{equation}
    \label{eq 17}
   \rho_A (- k_y \delta u_s) - \rho_B (k_y \delta u_s) = -E_{ac}\frac{2k_w\delta u_s}{z_{avg}}\frac{k_y^2}{n^2}(z_B-z_A) \sin (2 k_w x_s),
\end{equation}
Rearranging (\ref{eq 17}), the dispersion relation $n$ for the stability problem becomes
\begin{equation}
\label{eq 18}
    n = \sqrt{\frac{k_y}{\rho_A+\rho_B}\phi E_{ac}(z_B-z_A)\sin(2k_wx_s)}.
\end{equation}
where $\phi = 2k_w/z_{avg}$. The dispersion relation (\ref{eq 18}) establishes the acoustic stability criterion when inhomogeneous fluids (without interfacial tension) in a microchannel are subjected to a standing acoustic wave. If the eigenvalue $n$ is imaginary in (\ref{eq 18}), then the configuration is in a stable equilibrium and the configuration is in an unstable equilibrium when the eigenvalue $n$ is real. For a standing acoustic half-wave, in (\ref{eq 18}), the values of $\frac{k_{y}}{\rho_{B}+\rho_{A}}, \phi$ and $E_{a c}$ are always positive. Thus, the sign of $z_{B}-z_{A}$ (initial configuration of the fluids) and $\sin (2 k_{w} x_s)$ (relative location of the interface with respect to the standing acoustic wave) decide the nature of the eigenvalue in (\ref{eq 18}). $z_{B}-z_{A}$ is positive when high-impedance fluid is present to the right of the interface, and negative when high-impedance fluid is present to the left of the interface. $\sin (2 k_{w} x_s)$ has a negative value to the left of the node ($x_s$ is negative), a positive value to the right of the node ($x_s$ is positive) and zero when the interface coincides with the node (centre of the microchannel) or anti-node (sides of the microchannel) $(x_s = 0)$.

As per the above arguments, the inhomogeneous system in figure \ref{figure 1}($a$) is in an unstable equilibrium as eigenvalue $n$ is real, and the system in figure \ref{figure 1}($b$) is in a stable equilibrium as eigenvalue $n$ is imaginary. It can be concluded from the above discussion and figures \ref{figure 2}($a$-$i$) and \ref{figure 2}($a$-$ii$) that, a system is said to be acoustically stable (unstable) if the initial configuration of the fluids is in such a way that the low (high) impedance fluid is present at the anti-node(s) and the high (low) impedance fluid is present at the node(s). This conclusion is consistent with the demonstration of acoustic relocation of fluids within a microchannel by \citet{Deshmukh2014Jul}. For the case where the interface coincides with the node, $\sin (2 k_{w} x_s)=0$. Thus, the system is in a neutral equilibrium ($n=0$) as shown in figures \ref{figure 1}($c$) and \ref{figure 2}($a$-$iii$). The above analysis can be easily extended to an inhomogeneous system consisting of multiple interfaces. In this case, the eigenvalues evaluated at the fluid interfaces govern the nature of the system. Figure \ref{figure 2}($b$) shows the stability of two interface systems that are widely used in acoustofluidic applications. It can be seen from figure \ref{figure 2}($b$-$i$) that when high impedance fluid is at the sides (anti-nodes), the eigenvalue at both the interfaces $(IF_1$ and $IF_2)$ is real and hence the system is in unstable equilibrium. The system is in stable equilibrium in figure \ref{figure 2}($b$-$ii$), as the eigenvalue at both interfaces is imaginary.
\begin{figure}
    \center
    \includegraphics[width=1\linewidth]{figure2.png}
    \caption{\label{figure 2} Different inhomogeneous fluid configurations commonly used in microfluidics and their equilibrium nature. \textbf{a)} Single interface configurations; \textbf{b)} Double interface configurations. Equations (\ref{eq 1}) and (\ref{eq 2}) are used to calculate $n$ for fluids without interfacial tension and with interfacial tension.}
\end{figure}

\subsection{Stability analysis of inhomogeneous fluids in the presence of interfacial tension\label{Sec 2.2}}

Proceeding to solve for immiscible fluids, the effect of surface tension must be accounted for. The discontinuity in impedance occurring in the interfaces ($x_s$) is modelled by including the interfacial tension effects in the X momentum equation (\ref{eq 9b}) as, \citep{Chandrasekhar1961}
\begin{equation}
\label{eq 19}
    \rho_0 n\delta u = -\frac{\partial \delta p}{\partial x} - \frac{E_{ac}\cos{\left(2k_wx\right)}}{z_{avg}}\frac{\partial \delta z}{\partial x} - k_y^2\sum_{s}T\delta x_s\delta(x-x_s).
\end{equation}
where $T$ is the interfacial tension and $\delta x_s$ denotes the perturbation of the interfaces and $\frac{d}{dt}\delta x_s = \delta u_s \implies \delta x_s = \frac{\delta u_s}{n}$. The governing differential equation for stability between inhomogeneous fluids with interfacial tension is obtained similar to the case without interfacial tension, as in the previous $\S$ \ref{Sec 2.1},
\begin{multline}
\label{eq 20}
       \frac{d}{dx}\left(\rho_0 \frac{d\delta u}{dx} \right) - \rho_0 k_y^2 \delta u = - E_{ac}\frac{2k_w\delta u}{z_{avg}}\frac{k_y^2}{n^2}\sin(2k_wx)(z_B-z_A)\delta(x-x_s)\\+\frac{k_y^2}{n^2}\sum_{s}{k_y^2\left(T \delta u_s\right)\delta(x-x_s)},
\end{multline}
Integrating (\ref{eq 20}) across an infinitesimal distance ($dx \approx 0$) and solving for the dispersion relation $n$,
\begin{equation}
\label{eq 21}
    n = \sqrt{\frac{k_y}{\rho_1+\rho_2}\left(\phi E_{ac}(z_B-z_A)\sin(2k_wx_s)-k_y^2T\right)}.
\end{equation}
Equation (\ref{eq 21}) establishes the acoustic stability criterion when fluids with interfacial tension are subjected to a standing acoustic wave. It can be seen from (\ref{eq 21}) that the interfacial tension $(T)$ and wavenumber of the perturbation $(k_y)$ play a role in the stability of immiscible fluids.

In the presence of interfacial tension ($T>0$), the fluid system shown in figure \ref{figure 1}($b$) is always stable, as the negative sign of $(z_{B}-z_{A}) \sin (2k_{w}x_s)$ results in an imaginary eigenvalue $n$ in (\ref{eq 21}). Whereas, for the fluid system shown in figure \ref{figure 1}($a$), the sign of  $(z_{B}-z_{A}) \sin (2k_{w}x_s)$ is positive in (\ref{eq 21}). Thus, the system is conditionally stable, and the stability is determined by the relative magnitudes of $\phi E_{a c}(z_{B}-z_{A}) \sin (2 k_{w} x_s)$ and $k_y^{2} T$. The fluid system (figure \ref{figure 1}($a$) becomes unstable ($n$ is real) if the acoustic force density $F_{rl}$ ($\phi E_{a c}(z_{B}-z_{A}) \sin (2 k_{w} x_s)$) dominates (or is greater than) the interfacial force density $F_{int}$ ($k_y^{2} T$) and becomes stable ($n$ is imaginary) if the interfacial force density dominates the acoustic force density. For the case where the interface coincides with the node ($\sin (2 k_{w} x_s)=0$ and eigenvalue $n$ is imaginary) and the system is in a stable equilibrium, as shown in figures \ref{figure 1}($c$) and \ref{figure 2}($a$-$iii$).

Now, for conditionally stable configuration, we proceed to find the minimum energy density required to relocate the fluid systems in figures \ref{figure 1}($a$) and \ref{figure 2}($b$-$i$) with interfacial tension ($T>0$). Since the interface height, $h$ is finite, this leads to the quantization of the possible modes $k_y=k_{h_n}=n\pi/h$. The minimum (critical) acoustic energy density $(E_{cr})$ required to relocate the fluid system is decided by the first conceivable mode, $k_{h_1}=k_h=\pi/h$ and the critical acoustic energy density is obtained by limiting the eigenvalue $n$ to zero in (\ref{eq 21}). Thus,
\begin{equation}
\label{eq 22}
    E_{cr} = \frac{k_h^2Tz_{avg}}{sin(2k_wx_s)2k_w(z_B-z_A)}.
\end{equation}
If the applied energy density $E_{ac}$ is less than the critical energy density $E_{cr}$ ($E_{a c}<E_{c r}$), the interfacial tension succeeds in stabilizing a potentially unstable configuration. The same system becomes unstable and eventually relocates to a stable configuration when $E_{a c} > E_{c r}$. The above discussions on the equilibrium nature of different inhomogeneous fluid configurations (with and without interfacial tension) are clearly summarised in figure \ref{figure 2}.
\section{Numerical results and discussion}\label{Sec 3}
In this section ($\S$ \ref{Sec 3}), we numerically analyze the stability of inhomogeneous fluids (with and without interfacial tension) under acoustic fields and compare them with the results obtained by the theoretical analysis in the previous section ($\S$ \ref{Sec 2}). At first, we study the stability and relocation using acoustic relocation force $\textbf{\emph{f}}_{rl}$ (\ref{eq 4}) where the acoustic energy density is assumed to be constant (as the variation of first-order pressure and velocity are not considered). We further extend the numerical analysis using the generalized acoustic body force $\textbf{\emph{f}}_{ac}$ (\ref{eq 2}) where the first-order pressure and velocity vary during relocation (thus $E_{ac}$ varies) \citep{Rajendran2022Jun}.

The numerical analysis is carried out on a two-dimensional fluid domain of height $h=160\ \mu m$ and width $w= 360\ \mu m$ in COMSOL Multiphysics 6.0. For this study, the fluids mineral oil $(Z=1.23\ MPa\ s/m)$ and silicone oil $(Z=0.961\ MPa\ s/m)$ are used. A mesh refinement procedure, similar to those employed by \citet{Rajendran2022Jun} is used to confirm that the numerical findings are not affected by grid size. Three different fluid configurations are considered for the study, namely,
\begin{itemize}
    \item High-Low-High (HLH) configuration where the high impedance fluid is present at the anti-nodes (sides) and the low impedance fluid is present at the node (center) as shown in figure \ref{figure 3}($a$).
    \item Low-High-Low (LHL) configuration where the low impedance fluid is present at the anti-nodes (sides) and the high impedance fluid is present at the node (center) as shown in figure \ref{figure 3}($b$).
    \item High-Low (HL) configuration where the high impedance fluid occupies the domain to the left of the center of the microchannel and the low impedance fluid occupies the domain to the right of the center of the microchannel as shown in figure \ref{figure 3}($c$).
\end{itemize}

For the sake of brevity, the configurations shown in figures \ref{figure 2}($a$-$i$) (or \ref{figure 1}) and \ref{figure 2}($a$-$ii$) are not discussed explicitly as their stability and relocation are captured by HLH and LHL configurations. The Low-High (LH) configuration is also not discussed, as it would be analogous to the HL configuration. For all the analyses, the initial interface is perturbed and modelled as $x_s(y)= A_0\cos\left(\frac{2\pi}{h}y+\frac{h}{2}\right)$, where $A_0=0.01h$ is the perturbation amplitude.

\begin{figure}
\center
    \includegraphics[width=1\linewidth]{figure3.png}
    \caption{\label{figure 3} Stabilization and relocation of inhomogeneous fluids using simplified body force (\ref{eq 4}) with constant $E_{ac}$ - \textbf{(a)} High-Low-High (HLH) configuration, \textbf{(b)} Low-High-Low (LHL) configuration,
    \textbf{(c)} High-Low (HL) configuration}
\end{figure}

\subsection{Numerical analysis of stability using constant acoustic energy density}\label{Section 3.1}

For the numerical simulations shown in \ref{figure 3}, we employ equation (\ref{eq 4}) as body force and assumed $E_{ac}$ to be constant throughout the relocation process. The boundary condition for the analysis is no slip at the walls and the pressure is constrained at a point (bottom left corner of the channel). In the absence of interfacial tension $(T=0\ mN/m)$, it is observed that for any $E_{ac}>0$, the HLH configuration undergoes relocation to a stable LHL configuration as in figure \ref{figure 3}($a$) (the simulation is shown for $E_{ac}$ $=80$ $J/m^3$). In this case, the magnitude of $E_{ac}$ only influences the timescale of the relocation process by competing with the viscosity. While, in the presence of interfacial tension, $T=1\ mN/m$, the HLH fluid configuration remained stable for all energy densities below $88\ J/m^3$, and relocation is observed for all energy densities above $89\ J/m^3$. These simulations are in close agreement with the critical acoustic energy density $E_{cr}=88.78\ J/m^3$ predicted by  (\ref{eq 22}) for mineral-silicon oil combination. Simulation results of other fluid combinations shown in figure \ref{figure 5} also agree with (\ref{eq 22}). When the applied $E_{ac}$ is just above $E_{cr}$, the fluids take a much longer time to relocate. Thus for convenience, the simulation is  shown for $E_{ac}=120\ J/m^3$ in figure \ref{figure 3}($a$).

For LHL configuration with and without interfacial tension ($T\ge0$), for any $E_{ac}>0$, the relocation of fluid is not observed, and the system remained stable as shown in figure \ref{figure 3}($b$) (the simulation is shown for $E_{ac}$ $=120$ $J/m^3$). In the HL configuration, the node of the standing acoustic half-wave coincides with the fluid-fluid interface. Here for fluids with interfacial tension, relocation is not observed for any $E_{ac}>0$, and the fluid system remained stable (figure \ref{figure 3}($c$)). Whereas, for fluids without interfacial tension, the HL configuration is observed to be in neutral equilibrium (figure \ref{figure 3}($c$)).
These simulation results of unstable, stable, and neutral equilibrium of inhomogeneous fluids (figure \ref{figure 3}) are in agreement with the stability criteria (from (\ref{eq 18}) and (\ref{eq 21})) that we established theoretically in $\S$ \ref{Sec 2}.

\subsection{Numerical analysis of stability using generalized body force $\textbf{f}_{ac}$\label{Section 3.2}}

\begin{figure}
    \center
    \includegraphics[width=1\linewidth]{figure4.png}
    \caption{\label{figure 4}
    Stabilization and relocation of inhomogeneous fluids using generalized body force $f_{ac}$ (\ref{eq 2}) along with the first-order pressure field ($|p_1|=\sqrt{Real(p_1^{\star}p_1)}$) for different fluid configurations. 1-D actuation is imposed on \textbf{(a-c)}, and 2-D actuation is imposed on \textbf{(d)}. \textbf{(a)} HLH configuration remained stable up to $E_{ac}=85.58\ J/m^3$ ($p_a = 0.67\ MPa$, $d_0$ = $0.21$ $nm$, $\nu=$  $1.73$ $M Hz$). \textbf{(b)} HLH configuration undergoes relocation above $E_{ac}=86.22$ $J/m^3$ ($p_a = 0.69\ MPa$, $d_0$ = $0.22$ $nm$, $\nu=$  $1.73$ $M Hz$). Significant  variation in $|p_1|$ during relocation is observed. \textbf{(c)} HL configuration where the fluid interface coincides with the node remained in stable equilibrium even at $E_{ac}=2334\ J/m^3$ ($p_a = 3.54\ MPa$, $d_0$ = $20$ $nm$, $\nu=$  $1.73$ $M Hz$). \textbf{(d)} Relocation of HL configuration due to 2-D wall actuation $(p_a=1.37\ MPa,\  d_0=20\ nm,\ \nu=2.1\ MHz)$
  }
\end{figure}

Thus far, in the theoretical stability analysis ($\S$ \ref{Sec 2}) as well as in the numerical simulations ($\S$ \ref{Section 3.1}), a simplified equation $\textbf{\emph{f}}_{rl}$ (\ref{eq 4}) is employed as a body force with the assumption of constant $E_{ac}$ (the amplitudes of first-order fields $p_a$ and $v_a$ do not vary during relocation). In this section, the generalized acoustic body force $\textbf{\emph{f}}_{ac}$ (\ref{eq 2}) is employed and the first-order fields required to calculate the above $\textbf{\emph{f}}_{ac}$ are obtained from the wave equations (frequency domain - see appendix \ref{Sec A1}) by actuating the channel walls at a frequency $\nu$ with a wall displacement $d_0$. There are two reasons for using generalized acoustic body force $\textbf{\emph{f}}_{ac}$: 1. To show the relocation predicted by $\textbf{\emph{f}}_{rl}$ and $\textbf{\emph{f}}_{ac}$ is approximately the same. When we use much simpler $\textbf{\emph{f}}_{rl}$ instead of the complex $\textbf{\emph{f}}_{ac}$, the first-order field equations are not required to be solved which will significantly reduce the computation time for simulation of relocation of inhomogeneous fluids. 2. To explain the previous microchannel experiments in immiscible fluid relocation \citep{Hemachandran2019Apr}.

For one-directional (1-D) standing half-wave simulations, the sidewalls are actuated in phase at a displacement $d_0$ at a frequency $\nu$. In laminar flow equations, the boundary conditions used are no-slip at all walls, and the pressure is constrained at a point (bottom left corner of the channel). To disregard the effect of streaming, the first-order acoustic fields (see Appendix \ref{Sec A1}) are allowed to slip in the frequency domain.

Figure \ref{figure 6}($a$), shows the HLH configuration subjected to 1-D standing half-wave by actuating sidewalls at a displacement $d_0$ of $0.21$ $nm$ and a frequency $\nu$ of $1.73\ MHz$. In this case, it is observed that the resulting pressure amplitude $P_a$ of $0.67\ MPa$ ($E_{ac}$ = $85.58\ J/m^3$), could not relocate the fluids in the HLH configuration and thus  remains stable. Whereas, when the displacement is increased to $0.22$ $nm$, the resulting pressure amplitude of $P_a=0.69\ MPa$ ($E_{ac}$ = $86.22\ J/m^3$), could relocate the HLH configuration to a stable equilibrium as shown in figure \ref{figure 6}($b$).
From the above discussion, the critical acoustic energy density is found to be  $E_{cr}$ = $85.9\ \pm 0.32$ $J/m^3$. This value of $E_{cr}$ obtained through generalized body force $\textbf{\emph{f}}_{ac}$ is in close agreement (deviation of $3.24\%$) with the simplified relocation force $f_{rl}$ employed to derive stability criterion ($\S$ \ref{Sec 2}).

In the case of a 1-D standing half-wave, when the interface of the fluid coincides with the pressure node ($x_s = 0$), for any $E_{ac}$, relocation is not observed using both $\textbf{\emph{f}}_{ac}$ and $\textbf{\emph{f}}_{rl}$ (figures \ref{figure 3}($c$) and \ref{figure 4}($c$) as predicted by the stability criteria (\ref{eq 21}). However, \citet{Hemachandran2019Apr} through experiments demonstrated the relocation of fluids irrespective of the location of the vertical interface $x_s$. In their experiments, the frequency employed ($2.1$ $MHz$) is far from the 1-D resonant half-wave frequency ($\nu$ = $1.6$ $MHz$ $\approx$ $c_{avg}/2w$). In our previous work \citep{Rajendran2022Jun}, we have shown that the above relocation is due to standing two-directional (2-D) acoustic wave (frequency $f$ $=$ $2.1\ MHz$ between $c_{a v g} / 2 w$ and $c_{a v g} / 2h$) as shown in figure \ref{figure 4}($d$). From figure \ref{figure 4}($d$) it is clear that the pressure node is not vertical but inclined with respect to the fluid interface owing to the 2-D actuation (all four walls are actuated at $d_0$).
The above 2-D relocation can be clearly explained by the fact that if the fluid interface and node are not perpendicular to each other, then
$\boldsymbol{\nabla} \times \textbf{\emph{f}}_{rl}$ $\neq  0$. This implies when a sufficient energy density is applied, the fluid system \ref{figure 4}($d$) will not be in equilibrium and relocation begins without imposing any perturbations unlike the other relocation discussed in this work.

\subsection{Characterization of stable and unstable (relocation) regime\label{Section 3.3}}
When the 1-D acoustic standing wave is imposed on fluids with interfacial tension, the configurations (figures \ref{figure 1}($b$), \ref{figure 2}($b$-$i$), \ref{figure 3}($a$) having high impedance fluid at the anti-node and the low impedance fluid at the node, become conditionally stable. From (\ref{eq 21}), it is evident that the stability of the above inhomogeneous fluid configurations is governed by the ratio of $F_{rl}$ and $F_{int}$, which is called as acoustic Bond number ($Bo_{a}$), given by
\begin{equation}
\label{eq 3.1}
Bo_{a}=\frac{F_{rl}}{F_{int}}=\frac{\phi E_{a c} \Delta Z \sin \left(2 k_{w} x_s\right)}{k_{h}^{2} T}
\end{equation}
The $Bo_a$ that separates the stable and unstable region is called critical acoustic Bond number $Bo_{a,cr}$. From (\ref{eq 21})
\begin{equation}
    \label{eq 3.2}
    Bo_{a,cr}=1
\end{equation}
For $Bo_a>Bo_{a,cr}$, the above configurations become unstable (relocation occurs), and for $Bo_a<Bo_{a,cr}$ the configurations remain stable. Figure \ref{figure 5} shows the simulation results of different immiscible fluid combinations. The relocation and non-relocation regimes predicted by the simulations are in line with (\ref{eq 3.2}). It must also be noted that the fluids with higher interfacial tension require a higher acoustic force for relocation.

\begin{figure}
    \center
    \includegraphics[width=0.7\linewidth]{figure5.png}
    \caption{\label{figure 5} Characterization of relocation and non-relocation regimes of immiscible fluids using acoustic Bond number $Bo_{a}$}
\end{figure}

\subsection{Effect of the height of the channel on stability\label{Section 3.4}}

\begin{figure}
\center
    \includegraphics[width=1\linewidth]{figure6.png}
    \caption{\label{figure 6} Effect of channel height on stability, \textbf{(a)} channel height $h= 80\ \mu m$ - no relocation is observed as applied energy density $(E_{ac} = 120\ J/m^3)$ is less than the critical energy density $(E_{cr} = 384\ J/m^3)$,  \textbf{(b)} channel height $ h = 160\ \mu m$ - relocation is observed as applied energy density $(E_{ac} = 120\ J/m^3)$ is high than the critical energy density $(E_{cr} = 88.78\ J/m^3)$.
    This demonstrates that the interfacial tension force weakens as the height of the channel increases and thus higher the height of the channel, the lower the $E_{ac}$ required for relocation as $(E_{cr}\propto1/h^2)$.}
\end{figure}

The height of the channel $h$ plays a critical role in the stability of immiscible fluids. For a given $E_{ac}$, the increase in $h$  weakens the stabilizing effect of interfacial tension force, as analysed theoretically in $\S$ \ref{Sec 2}.
From (\ref{eq 22}), it can be inferred that the critical acoustic energy density is inversely proportional to the square of the channel height $(E_{cr} \propto 1/h^2)$. In figure \ref{figure 6}($a$), for a microchannel of height $h=80\ \mu m$ consisting of mineral-silicone oil with interfacial tension $T=1\ mN/m$, the fluid system is stable as the applied $E_{ac}\ (120\ J/m^3$) is lower than the critical energy density ($E_{cr}=384\ J/m^3$). Whereas for $h=160\ \mu m$ and keeping the remaining parameters same, fluid relocation is observed as applied $E_{ac}\ (120\ J/m^3)$ is higher than the critical energy density $(E_{cr}=88.78\ J/m^3)$.

The above discussion on the effect of channel height on acoustic relocation has high relevance in practical applications. To relocate fluids with high interfacial tension of $\mathcal{O}(10^1-10^2)$ $mN/m$, in commonly used acoustofluidic channels of height ranging from $100\ \mu m$ to $200\ \mu m$, the required $E_{ac}$ becomes $\approx \mathcal{O}(10^4)\ J/m^3$, which is very high compared to the $E_{ac}$ employed in typical acoustofluidic experiments ($\mathcal{O}(10^2-10^3)$ $J/m^3$). The equation (\ref{eq 22}) tells that the above problem can be solved by increasing the channel height as $E_{ac} \propto 1/h^2$. Hence, the depth (height) of the channel is a crucial aspect to be considered during the fabrication of an acoustofluidic microchannel for handling high interfacial tension fluids.

\section{Conclusion}
We have theoretically established the stability criteria for inhomogeneous fluids subjected to standing acoustic fields, which is consistent with the previous experimental investigations on miscible \citep{Deshmukh2014Jul}, \citep{Karlsen2016Sep}, and immiscible fluids \citep{Hemachandran2019Apr}. Numerical simulations on the same were carried out using simplified and generalized body force to understand the various parameters that contribute towards stability and relocation of fluids. However, the effect of boundary layer-driven streaming on relocation is neglected in this work, which will be addressed in an upcoming paper. The insights gained from this study can have potential applications in inhomogeneous fluid handling and particle manipulation in the field of acoustofluidics.