## 3 The Variational Method  
To gain an analytical insight into the condensate dynamics, we use the variational approximation (VA) with the Gaussian ansatz as a trial wave function for the solution of Equation (6) in the absence of an external trap [23]:  
\[\psi(\rho,t)=\frac{1}{R(t)\sqrt{\pi}}\exp\left(-\frac{\rho^{2}}{2R(t)^{2}}+i \beta(t)\rho^{2}\right), \tag{8}\]  
Here the variables \(R(t)\) and \(\beta(t)\) are real-valued radius and radial chirp of the self-trapped state, while the amplitude is determined by the normalization condition (2). The Lagrangian density generating Equation (6) in the absence of the trap (\(d(t)=0\)) is given by  
\[\mathcal{L} = \frac{i}{2}\bigg{(}\psi\frac{\partial\psi^{*}}{\partial t}-\psi^{ *}\frac{\partial\psi}{\partial t}\bigg{)}+\frac{|\
abla_{\rho}\psi|^{2}}{2}+ \frac{2\pi a(t)N}{\sqrt{2\pi}d_{z}}|\psi|^{4}+\frac{\chi(t)}{3\sqrt{3}\pi d_{ z}^{2}}|\psi|^{6} \tag{9}\] \[\qquad+\frac{2\pi a_{dd}N}{\sqrt{2\pi}d_{z}}|\psi|^{2}\int\frac{ d\mathbf{k}_{p}}{(2\pi)^{2}}e^{i\mathbf{k}_{p}\cdot\rho}n(\mathbf{k}_{p},t)h_{2D}( \sigma)\bigg{)}.\]  
where \(\psi^{*}\) denotes the complex conjugate of the wave function \(\psi\). The trial wave function (8) is substituted into the Lagrangian density (9) and the effective Lagrangian is calculated by integrating the Lagrangian density. Then, the following Euler-Lagrangian equations for the variational parameters can be derived from the effective Lagrangian as  
\[\frac{dR(t)}{dt} = 2R(t)\beta(t), \tag{10}\] \[\frac{d\beta(t)}{dt} = \frac{1}{2R(t)^{4}}+\frac{Na(t)}{\sqrt{2\pi}d_{z}R(t)^{4}}-\frac {Na_{dd}\Lambda(\xi(t))}{2\sqrt{2\pi}d_{z}R(t)^{4}}+\frac{2N^{2}\chi(t)}{9\sqrt {3}d_{z}^{2}\pi^{3}R(t)^{6}}-2\beta(t)^{2}. \tag{11}\]  
By combining Equations (10) and (11), we obtain the following equation for the evolution of the width \(R(t)\) as  
\[\frac{d^{2}R(t)}{dt^{2}}= \frac{1}{R(t)^{3}}+\frac{\sqrt{2}Na(t)}{\sqrt{\pi}d_{z}R(t)^{3}}- \frac{Na_{dd}\Lambda(\xi(t))}{\sqrt{2\pi}d_{z}R(t)^{3}}+\frac{4N^{2}\chi(t)}{9 \sqrt{3}d_{z}^{2}\pi^{3}R(t)^{5}} \tag{12}\]  
where \(\Lambda(\xi(t))=2-7\xi(t)^{2}-4\xi(t)^{4}+9\xi(t)^{4}d(\xi(t))/(1-\xi(t)^{2})^ {2}\), \(\xi(t)=R(t)/d_{z}\). We intend to consider the effects of periodic modulation of the contact interaction in the form of  
\[a(t)=\epsilon_{0}+\epsilon_{1}\sin\left(\Omega t\right),\ \chi(t)=\chi_{0}+\chi_{1}\sin \left(\Omega t\right) \tag{13}\]  
on the stability of DBEC where \(\epsilon_{0}\), \(\chi_{0}\) are constant parts of the two-body and three-body interactions, respectively, and \(\epsilon_{1}\), \(\chi_{1}\) are the amplitudes of their temporal modulation. The Kapitza averaging method can be employed to treat the oscillatory parts [30]. With the explicitly included oscillating nonlinearity, we get the following equation for the evolution of the radial width (cf. Ref. [23]):  
\[\frac{d^{2}R(t)}{dt^{2}}=\frac{1}{R(t)^{3}}+\frac{\sqrt{2}N[\epsilon_{0}+ \epsilon_{1}\sin\left(\Omega t\right)]}{\sqrt{\pi}d_{z}R(t)^{3}}-\frac{Na_{dd} \Lambda(\xi(t))}{\sqrt{2\pi}d_{z}R(t)^{3}}+\frac{4N^{2}[\chi_{0}+\chi_{1}\sin \left(\Omega t\right)]}{9\sqrt{3}d_{z}^{2}\pi^{3}R(t)^{5}} \tag{14}\]  
Now, \(R(t)\) may be separated into a slowly varying part \(R_{0}(t)\) and a rapidly varying one \(R_{1}(t),R(t)=R_{0}(t)+R_{1}(t)\).