Scale-invariant relaxation dynamics in two-component Bose-Einstein condensates with large particle-number imbalance

TL;DR

This paper investigates the scale-invariant relaxation dynamics in segregating two-component Bose-Einstein condensates with large particle-number imbalance, revealing the role of droplet random walk and deriving scaling laws through simulations and analytical models.

Contribution

It introduces a novel random walk model for droplets that explains the scale-invariant relaxation and predicts crossover phenomena in two-component Bose-Einstein condensates.

Findings

Emergence of dynamical scaling in correlation functions.

Derivation of 1/3 and 1/2 power laws for scaling exponents.

Identification of a crossover in the scaling behavior.

Abstract

We theoretically study the scale-invariant relaxation dynamics in segregating two-component Bose-Einstein condensates with large particle-number imbalance, and uncover that random walk of droplet for the minor component plays a fundamental role in the relaxation process. Our numerical simulations based on the binary Gross-Pitaevskii model reveal the emergence of the dynamical scaling during the relaxation, which is a hallmark of scale-invariant dynamics, in a correlation function for the minor condensate. Tracking exponents characterizing the dynamical scaling in time, we find out a crossover phenomenon that features the change in power exponents of the correlation length. To understand the fundamental mechanism inherent in the scale-invariant relaxation dynamics, we construct a random walk model for droplets. Employing the model, we analytically derive the $1/3$ and $1/2$ power laws and predict the crossover of the scaling. These exponents are in reasonable agreement with the values obtained in the numerical calculations. We also discuss a possible experimental setup for observing the scale-invariant dynamics.

Tables

Table 1. Table 1: Dependence of Λ Λ \Lambda and L ′ / ξ superscript 𝐿 ′ 𝜉 L^{\prime}/\xi in Fig.  3 on the imbalance parameter R 𝑅 R . The third column shows Λ Λ \Lambda for R = 0.25 𝑅 0.25 R=0.25 and 0.35 0.35 0.35 divided by that for R = 0.5 𝑅 0.5 R=0.5 . Despite the naive estimation of Λ Λ \Lambda , the values in the third column are reasonably close to those in the fourth one.

$R$		$\Lambda$		$\Lambda \left(R\right)/\Lambda \left(0.5\right)$		$L^{\prime }/\xi$
0.5		10.1		1		1
0.35		11.4		1.13		1.10
0.25		13.0		1.29		1.14

Equations

$$i\hslash \frac{\partial }{\partial t}{\psi }_1(\mathbf{r},t)$$

$$i\hslash \frac{\partial }{\partial t}{\psi }_2(\mathbf{r},t)$$

$$N_m={\int }_{\mathcal{S}}|{\psi }_m(\mathbf{r},t){|}^{2}d\mathbf{r}$$

$${\psi }_1(\mathbf{r},0)$$

$${\psi }_2(\mathbf{r},0)=\frac{1}{\sqrt{S}}\sum _{\mathbf{k}\in \{\mathbf{k}|{ϵ}_{\mathbf{k}}\le \mu ,\mathbf{k}\ne 0\}}{\alpha }_{2,\mathbf{k}}{e}^{i\mathbf{k}\cdot \mathbf{r}},$$

$$u_{\mathbf{k}}=\sqrt{\frac{{ϵ}_{\mathbf{k}}+\mu }{2E_{\mathbf{k}}}+\frac{1}{2}},$$

$$v_{\mathbf{k}}=-\sqrt{\frac{{ϵ}_{\mathbf{k}}+\mu }{2E_{\mathbf{k}}}-\frac{1}{2}},$$

$$E_{\mathbf{k}}=\sqrt{{ϵ}_{\mathbf{k}}({ϵ}_{\mathbf{k}}+2\mu )},$$

$$\displaystyle W(\{\alpha_{j,\bm{k}}\})=A\prod_{\bm{k}}{\rm exp}\Bigl{(}-2|\alpha_{j,\bm{k}}|^{2}\Bigl{)}$$

$$\Omega (t)=\frac{1}{\tau }\theta ({\tau }_{\mathrm{hold}}-t)\theta (t)$$

$$\displaystyle C(r,t)=\biggl{\langle}\frac{1}{S}\int_{\mathcal{S}}\psi_{2}^{*}({\bm{r}}+{\bm{x}},t)\psi_{2}({\bm{x}},t)d\bm{x}\biggl{\rangle}_{\Omega},$$

$$C(L(t),t)=0.5C(0,t).$$

$$C(r,t)=s^{-\gamma }C(r{s}^{\beta },ts)$$

$$C(r,t)=t^{\gamma }G(r{t}^{-\beta })$$

$$L(t)$$

$$C(0,t)$$

$$\frac{{\lambda }_{\mathrm{in}}}{\xi }=2\pi \sqrt{\frac{2\mu }{-\mu +\sqrt{g_0^2({\rho }_1-{\rho }_2{)}^2+4{\rho }_1{\rho }_2{g}_1^2}}}$$

$$(\beta ,\gamma )\simeq \{\begin{array}{ll}(0.31,0)& (t<t^{\ast }),\\ (0.44,0)& (t^{\ast }<t)\end{array}$$

$$f_{\nu }=\frac{1}{\rho }\sum _{m,n=1,2}{\psi }_m^{\ast }({\sigma }_{\nu }{)}_{mn}{\psi }_n,$$

$$\rho =|{\psi }_1{|}^2+|{\psi }_2{|}^2$$

$$\frac{\partial }{\partial t}\rho f_z+\nabla \cdot {\mathbf{J}}_z=0,$$

$$\displaystyle{\bm{J}}_{z}=\rho\biggl{[}f_{z}\bm{v}+\frac{\hbar}{2M}\Big{(}f_{y}{\bm{\nabla}}f_{x}-f_{x}{\bm{\nabla}}f_{y}\Big{)}\biggl{]},$$

$$\displaystyle\bm{v}=\frac{\hbar}{2M\rho i}\sum_{m=1}^{2}\Bigl{(}\psi_{m}^{*}\bm{\nabla}\psi_{m}-\psi_{m}\bm{\nabla}\psi_{m}^{*}\Bigl{)}.$$

$$\mathbf{X}(t)=\frac{{\int }_{\mathcal{S}}\mathbf{r}D(\mathbf{r},t)d\mathbf{r}}{{\int }_{\mathcal{S}}D(\mathbf{r},t)d\mathbf{r}},$$

$$\displaystyle D(\bm{r},t)=\frac{1}{2}\Bigl{(}1-f_{z}(\bm{r},t)\Bigl{)}.$$

$$S_{\mathrm{dp}}\simeq {\int }_{\mathcal{S}}D(\mathbf{r},t)d\mathbf{r}.$$

$${\int }_{\mathcal{S}}{f}_z(\mathbf{r},t)d\mathbf{r}=\mathrm{Const}.$$

$$\frac{d}{dt}{X}_{\nu }(t)$$

$$\frac{d}{dt}{X}_{\nu }(t)\simeq \frac{1}{2{\rho }_{\mathrm{tot}}{S}_{\mathrm{dp}}}{\int }_{{\mathcal{S}}_{\mathrm{wall}}}{r}_{\nu }\nabla \cdot {\mathbf{J}}_z(\mathbf{r},t)d\mathbf{r}.$$

$$\frac{{\mathbf{J}}_z}{{\rho }_{\mathrm{tot}}{C}_{\mathrm{ch}}}=\frac{f_z\mathbf{v}}{C_{\mathrm{ch}}}+f_y\bar{\nabla }{f}_x-f_x\bar{\nabla }{f}_y$$

$$\displaystyle\biggl{\langle}\frac{|f_{z}\bm{v}|}{C_{\rm ch}}\biggl{\rangle}_{\rm wall}$$

Full text

I Introduction

Ultracold atomic gases have been one of the ideal playgrounds for studying nonequilibrium phenomena in quantum systems because of the high tunability of physical parameters and long enough time scales for observing the dynamics in experiments Bloch et al. (2008); Polkovnikov et al. (2011); Eisert et al. (2015). In fact, a variety of fundamental issues such as thermalization Kinoshita et al. (2006); Gring et al. (2012); Trotzky et al. (2012), phase transition Sadler et al. (2006); Weiler et al. (2008); Nicklas et al. (2015a); Lamporesi et al. (2013); Navon et al. (2015); Braun et al. (2015); Clark et al. (2016); Anquez et al. (2016), and turbulence Henn et al. (2009); Kwon et al. (2014); Kang et al. (2017); Navon et al. (2016); Mossman et al. (2018) have been investigated. Among them, universal aspects of nonequilibrium phenomena have been attracting a lot of attention. A typical example is scale-invariant relaxation dynamics, features of which are captured by a correlation function $C(r,t)$ that exhibits dynamical scaling $C(r,t)=s^{-\gamma}C(rs^{\beta},ts)$ with a constant $s$ Cardy (1996); Hohenberg and Halperin (1977); Onuki (2002); Bray (1994); Täuber (2014). The exponents $\beta$ and $\gamma$ classify universality of relaxation dynamics. The recent works in ultracold quantum gases reveal the existence of several universality classes in terms of coarsening dynamics Damle et al. (1996); De et al. (2014); Kudo and Kawaguchi (2011, 2013, 2015); Hofmann et al. (2014); Williamson and Blakie (2016a, b); Bourges and Blakie (2017); Symes and Blakie (2017); Williamson and Blakie (2017); Shitara et al. (2017); Takeuchi et al. (2015); Takeuchi (2016, 2018, 2019); De et al. (2014); Phuc et al. (2017); Fujimoto et al. (2018); Liu et al. (2016, 2018); Schmied et al. (2019a); Pietraszewicz et al. (2019) and non-thermal fixed points (NTFPs)Nowak et al. (2011, 2012); Schole et al. (2012); Karl et al. (2013a, b); Karl and Gasenzer (2017); Schmied et al. (2019b); Fujimoto et al. (2019); Erne et al. (2018); Prüfer et al. (2018). In the context of coarsening dynamics, segregating binary Bose-Einstein condensates (BECs) are shown to have $\beta=2/3$ and belong to the universality class for classical binary liquids in a inertial regime Bray (1994) when the particle numbers in the two condensates are equal Kudo and Kawaguchi (2013); Hofmann et al. (2014); Williamson and Blakie (2016a).

This kind of universal relaxation with particle-number balance has been originally studied in classical binary liquid at a critical quench for a long time Onuki (2002); Bray (1994). However, many systems with particle-number imbalance have given another essential testbed for comprehending scale-invariant relaxation dynamics. In fact, Ostwald ripening, which is an important phenomenon in such a system, has been actively studied. In this phenomenon, many small droplets of the minor component are nucleated and subsequently grow or shrink in time. One of the established theories is originally developed by Lifshitz and Slyozov Lifshitz and Slyozov (1961) and Wagner Wagner (1961), who uncovered that diffusion processes (evaporation and concentration) are of importance in the dynamics and that the averaged droplet size is governed by a $\beta=1/3$ power law. The original idea for the theory has been improved and extended to other systems by many researchers Rogers and Desai (1989); Farrell and Valls (1991); Yao et al. (1993); Gorentsveig et al. (1997); Huse (1986), and thus Ostwald ripening has played a fundamental role in developing a basis for the universal dynamics.

Against its backdrop, we bring up an interesting question; “Does a kind of Ostwald ripening emerge even in quantum gases?” Usually, ultracold atomic gases are well separated from the environment, and are regarded as isolated quantum systems. Thus, the total energy is conserved in contrast to the corresponding classical systems, and the diffusion process mentioned above is generally prohibited. This implies that droplets cannot shrink nor grow through the Ostwald ripening mechanism. On the other hand, coarsening dynamics actually proceeds even in the absence of diffusive processes in a binary BEC with particle-number balance. We expect that similar coarsening dynamics also occurs in a particle-number imbalanced system. It is then natural to ask what kind of universal phase-ordering kinetics emerges.

In this paper, we theoretically find that the particle-number imbalance greatly affects the universal aspects of the relaxation dynamics in two-dimensional (2D) segregating binary BECs, and demonstrates that the relaxation mechanism is essentially different from Ostwald ripening in classical systems. Our theoretical model is a two-component Gross-Pitaevskii (GP) equation, and we firstly perform the numerical simulations starting from initial states generated by the truncated Wigner method Sinatra et al. (2002); Blakie et al. (2008); Polkovnikov (2010). Then, we reveal that many droplets of the minor condensate are formed and that the correlation function of the wavefunction for the minor component shows dynamical scaling. Furthermore, we numerically find a crossover of the power law: the exponent $\beta$ is close to $1/3$ in the early stage of the dynamics and approaches to $1/2$ in the late time. To unveil a mechanism behind the scale-invariant dynamics, we employ an idea that random collisions between droplets promote the relaxation dynamics, and then analytically derive $1/3$ and $1/2$ power laws for the correlation length, which well explain the numerical results. Therefore, we conclude that the droplet random walk is the origin of the power laws, and argue that the mechanism is quite different from Ostwald ripening in classical systems in which the diffusion processes are crucial. Here, it is noteworthy that there are a few theoretical works related to a percolation problem in binary BECs with the imbalance Takeuchi et al. (2015); Takeuchi (2016); Shitara et al. (2017), but the several things mentioned above is not yet to be considered.

The rest of paper is organized as follows. In Sec. II, we describe our theoretical tools such as the binary GP model, a tuning protocol for the particle-number imbalance, and the dynamical scaling in a correlation function. Section III shows our numerical results, and we find that the particle-number imbalance gives a significant impact on the scale-invariant relaxation. To understand the universal aspects of the relaxation, in Sec. IV, we analytically obtain two power exponents for the growth laws of the correlation length, and unveil the fundamental relaxation mechanism. In Sec. V, we explain the relation between the universal relaxation dynamics found here and previous literature, and discuss an experimental setup for observing our theoretical prediction. Finally, we summarize our results in Sec. VI.

II Theoretical formulation

II.1 Model

We consider a two-component BEC close to zero temperature in a 2D uniform system. The system is well described by two macroscopic wavefunctions $\psi_{m}(\bm{r},t)~{}(m=1,2)$ obeying the following binary Gross-Pitaevskii (GP) equation Pethick and Smith (2008); Pitaevskii and Stringari (2016):

[TABLE]

with a particle mass $M$, intra- and inter-species interaction couplings $g_{0}>0$ and $g_{1}>0$, and a Rabi coupling $\Omega(t)$. The ground state with $\Omega(t)=0$ is immiscible (miscible) for $g_{1}>g_{0}$ ($g_{1}<g_{0}$), in which the two condensates are separated into two spatially distinct regions (are overlapped and spread in the whole region). In this work, we assume an immiscible binary system, i.e., $g_{1}>g_{0}$. As explained in Sec. II.3, the Rabi coupling $\Omega(t)$ is used to adjust a particle-number imbalance parameter $R=N_{2}/N_{1}$. Here, the particle number $N_{m}$ for the component $m$ is defined by

[TABLE]

with the whole system space $\mathcal{S}$.

II.2 How to prepare initial states

Our system is a 2D square box $L\times L$ with the periodic boundary condition. An initial wavefunction of Eqs. (1) and (2) is prepared by employing the truncated Wigner method with a Bogoliubov vacuum Sinatra et al. (2002); Blakie et al. (2008); Polkovnikov (2010). We consider a single condensed state of the component $1$ as the vacuum. Then, we generate the initial state by using the following expressions with random variables $\alpha_{j,\bm{k}}~{}(j=1,2)$:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\rho_{\rm tot}=N_{0}/S$ is the number density of the condensed particles with $N_{0}$ being the number of condensed particles and $S=L^{2}$ the system area, $\mu=g_{0}\rho_{\rm tot}$ is the chemical potential, and $\epsilon_{\bm{k}}=\hbar^{2}k^{2}/2M$ is the single-particle kinetic energy. The complex random numbers $\alpha_{j,\bm{k}}$ are sampled by the Wigner function $W(\{\alpha_{j,\bm{k}}\})$ corresponding to the Bogoliubov vacuum:

[TABLE]

with the normalization factor $A$.

II.3 Protocol for tuning the particle-number imbalance

We study how the particle-number imbalance affects the dynamical scaling behavior of a segregating binary mixture in the immiscible parameter regime. For this purpose, we use a protocol to tune $R$ as shown in Fig. 1(a), where the Rabi coupling is used to adjust the imbalance parameter $R$. Our work assumes the coupling $\Omega(t)$ is given by

[TABLE]

with the holding time $\tau_{\rm hold}$, the characteristic time $\tau=\hbar/\mu$, and the step function $\theta(x)$. Under the tuning protocol, the total particle number $N_{\rm tot}=N_{1}+N_{2}$ is conserved in time, and we can freely tune the parameter $R$ as shown in Fig. 1(b), where we show numerical results of Eqs. (1) and (2) on how the imbalance parameter $R$ changes in time. This method has already been implemented in several experiments for studying non-equilibrium phenomena in binary BECs Hamner et al. (2013); Nicklas et al. (2015b); Bersano et al. (2018). In this work, considering only $R<1$ in which the component $2$ is minority, we study relaxation dynamics, namely phase separation kinetics after the protocol.

II.4 Correlation function and dynamical scaling

Scale-invariant relaxation dynamics is well known to be captured by several correlation functions Cardy (1996); Hohenberg and Halperin (1977); Onuki (2002); Bray (1994); Täuber (2014). To quantify what extent our system exhibits scale-invariant properties, we introduce a correlation function for $\psi_{2}({\bm{r}},t)$:

[TABLE]

where $\langle\cdots\rangle_{\Omega}$ denotes an average over the angle of $\bm{r}$. The correlation function $C(r,t)$ takes only real values because of the angle average (see Appendix A for the proof). We define a correlation length $L(t)$ as

[TABLE]

The length $L(t)$ is a characteristic scale for the minor condensate and corresponds to a typical size of droplets generated after the tuning protocol as described in Sec. III.

Here, we comment on a density correlation function for $\rho_{2}(\bm{r},t)$. This correlation function includes information for a domain distribution around a droplet because $\rho_{2}(\bm{r},t)$ always gives positive contribution. On the other hand, in Eq. (11), the phase of the macroscopic wavefunction $\psi_{2}(\bm{r},t)$ on droplets are almost independent, and thus effects of the surrounding domains can be canceled due to the integral over space. Thus, Eq. (11) is suitable for investigating growth of the droplet size rather than using the density correlation function.

When the system has scale-invariant dynamics, the correlation function obeys the dynamical scaling:

[TABLE]

with a variable $s$. Substituting $s=t^{-1}$ in Eq. (13), we obtain

[TABLE]

with a function $G(x)=C(x,1)$. This means that the correlation length $L(t)$ and the on-site correlation $C(0,t)$ have power laws:

[TABLE]

In our system without the Rabi coupling, the particle number for each component is independently conserved in time, and thus the exponent $\gamma$ is zero. This fact is numerically confirmed in the next section.

This scaling law has been discussed in the context of critical dynamics close to second-order phase transitions in equilibrium states over several decades. However, several theoretical studies have recently found the emergence of such dynamical scalings in isolated quantum systems even without critical points, and the universality are investigated from the perspective of NTFPs Berges (2015); Nowak et al. (2011, 2012); Schole et al. (2012); Karl et al. (2013a, b); Karl and Gasenzer (2017); Schmied et al. (2019b); Berges and Sexty (2011); Piñeiro Orioli et al. (2015); Berges and Wallisch (2017); Mazeliauskas and Berges (2019); Fujimoto et al. (2019). So far, the two experiments Erne et al. (2018); Prüfer et al. (2018) report observations of the NTFP scenario in scalar and spinor Bose gases.

III Numerical results

We numerically study the relaxation dynamics with large particle-number imbalance by tuning several values of $R$ smaller than unity. This section describes scale-invariant aspects of the relaxation using the density distributions and the correlation functions. The parameters used in our numerics are given as follows: $\rho_{\rm tot}=15/\xi^{2}$, $g_{1}/g_{0}=1.5$ and $L=256\xi$ with the coherence length $\xi=\hbar/\sqrt{2M\mu}$. Using the pseudo-spectral method with the Fourth-order Runge-Kutta time evolution, we numerically solve Eqs. (1) and (2). The numerical resolutions for space and time are $\xi$ and $0.001\hbar/\mu$, respectively.

Figures 1(c) and 1(d) show time evolution of $\rho_{2}(\bm{r},t)=|\psi_{2}(\bm{r},t)|^{2}$ for $R=0.25$ and $0.5$ after the tuning protocol of Fig. 1(b). We find that many droplets are generated just after adjusting a target value of $R$. The droplet-nucleation dynamics can be well understood by dynamical instability predicted by the Bogoliubov theory (linear analysis for small fluctuations), from which the most unstable wavelength is given by

[TABLE]

with $\rho_{1}=\rho_{\rm tot}/(1+R)$ and $\rho_{2}=\rho_{\rm tot}R/(1+R)$. The derivation of Eq. (17) is given in Appendix B. Our calculations have $\lambda_{\rm in}/\xi\sim 15.2$ and $13.2$ for $R=0.25$ and $0.5$, respectively, which are consistent with the distributions in the left panels of Figs. 1 (c) and 1(d). In the late stage of the dynamics, the linear size of each droplet obviously increases in time, and tiny droplet configurations are kept for a longer time when the value of $R$ is smaller. However, as time goes by, such structures are gradually lost, and minor domains become large with fine surface fluctuations. Also, one finds a density hole in the domains, which is denoted by a solid arrow in the distributions at $t/\tau=800$ and $1600$ of Fig. 1(d). This structure is a quantized vortex for the component 2, and there the component 1 occupies the vortex core. Because of topological stability, the composite vortex survives until when it encounters its anti-vortex or goes outside of the domain. However, such a small structure is expected to hardly affect the domain growth dynamics.

To see the scale-invariant aspects of the relaxation dynamics, we calculate the correlation function (11), and find the dynamical scaling of Eq. (14) as shown in Fig. 2. The main panels of Fig. 2 are the plots of the raw numerical data, which clearly exhibit monotonous growth of the correlation lengths in time for both cases. The insets show the same data in the normalized ordinate and abscissa by the on-site correlations $C(0,t)$ and the correlation lengths $L(t)$ respectively, clearly demonstrating that the numerical data at the different times can collapse to a single curve. This behavior is consistent with the dynamical scaling of Eq. (14), which is a hallmark of scale-invariant relaxation dynamics.

The dynamical scaling is characterized by the exponents $\beta$ and $\gamma$. To determine $\beta$, we first plot the time evolution of $L(t)$. The inset of Fig. 3 shows the time evolution of $L(t)$ for $R=0.25$, $0.35$, and $0.5$. These curves seemingly show different time dependences, but they can be expressed by a single function when we divide the ordinate and abscissa for each $R$ by constants $L^{\prime}$ and $\tau^{\prime}$, respectively, as shown in the main panel of Fig. 3. This implies that the relaxation dynamics for the different imbalance parameters is identical to each other. From the figure, we find that the exponent $\beta$ is close to $0.31$ in the early time and approaches to $0.44$ in the late. These exponents are estimated by fitting a power law function $y=ax^{b}$ to the data with the Marquardt-Levenberg method in Gunplot. Here, $a$ and $b$ are the fitting parameters. The fitting regions for the two exponents are determined by eye. On the other hand, as we mentioned in Sec. II.4, the exponent $\gamma$ is zero because the particle number $N_{2}$ is conserved in time after the tuning protocol and $C(0,t)$ is identical to $N_{2}/S$. Thus, we obtain the following scaling exponents:

[TABLE]

with the crossover time $t^{*}$. As described in the next section, the fundamental mechanism behind the crossover is that the droplets randomly move with the droplet-size-dependent velocity.

We also investigate the density correlation function in the relaxation dynamics, and find that the correlation function does not obey the dynamical scaling when $R$ is far from unity. We expect that this is attributed to existence of two length scales: one is a typical domain size and the other is an average distance between droplets. As mentioned in Sec. II.4, the density correlation function includes contributions from a density distribution around a droplet, which related to the averaged distance. Thus, the correlation function cannot be described by the dynamical scaling because Eq. (14) means that the system is characterized only by a single length scale.

Before closing this section, we should note that the interaction-parameters in the simulations are far from the transition points ($g_{0}=g_{1}$). In this sense our result of the relaxation dynamics supports the NTFP conjecture described in Sec. II.4.

IV Analytical derivation of the power laws for the correlation length

We analytically derive power laws of the correlation length $L(t)$ using the GP equation (1) and (2). The essential idea is that many droplets of the minor component randomly move with the size-dependent velocity $V(L_{\rm dp}(t))$, where $L_{\rm dp}(t)$ is the linear size of a droplet. The derivation is divided into two steps. Firstly, focusing on the dynamics of a single droplet, we show that the velocity obeys two power laws $V(L_{\rm dp})\propto L_{\rm dp}^{-n}~{}(n=1,2)$. Secondly, assuming random collisions between droplets, we construct an equation of motion for the averaged droplet size. As a result, we derive the power exponent $\beta=1/3$ and $1/2$ for the early and late stages of the dynamics, respectively. Our theoretical results well describe the numerical results presented in Fig. 3.

IV.1 Dynamics of a single droplet

We investigate dynamics of a single droplet by employing a continuity equation of spin density vectors defined by

[TABLE]

[TABLE]

with the Pauli matrix $(\sigma_{\nu})_{mn}~{}(\nu=x,y,z)$. The continuity equation for $f_{z}$ Lamacraft (2008); Kudo and Kawaguchi (2010); Yukawa and Ueda (2012); Xu and Han (2012) is given by

[TABLE]

[TABLE]

where ${\bm{J}}_{z}$ is the spin current for $f_{z}$ and $\bm{v}$ being the mass current defined by

[TABLE]

The derivation of Eqs. (21) – (23) is explained in Appendix C. It is worthy to note that the total density distribution $\rho$ is almost uniform because of a small repulsion interaction $|g_{0}-g_{1}|/(g_{0}+g_{1})=0.2$ between the different components. In the following, we assume $\rho\simeq\rho_{\rm tot}$.

Here, we assume that there is a single droplet of the minor component in the system. The centroid ${\bm{X}}(t)$ of the droplet is defined by

[TABLE]

[TABLE]

The function $D(x)$ describes a reshaped droplet configuration as shown in Fig. 4(a), which has a sharp edge with a width in the order of $\xi$.

We then calculate the time derivative of ${\bm{X}}(t)$. We start with a time derivative of the denominator of Eq. (24) from the perspective of stability of the single droplet. By definition, the denominator is almost the same as the droplet area $S_{\rm dp}$ because the droplet has the flat profile of $f_{z}$ except for the edge:

[TABLE]

Since our model is symmetric under rotation about the $z$-axis in the spin space, the net longitudinal magnetization should be conserved:

[TABLE]

Thus, we can safely assume that the area $S_{\rm dp}$ is independent of time. This means that the single droplet itself is quite stable and the area is conserved although the shape can change in time.

Utilizing Eqs. (21) and (26), we calculate the time derivative of $X_{\nu}(t)~{}(\nu=x,y)$:

[TABLE]

Here, note that the divergence $\bm{\nabla}\cdot\bm{J}_{z}(\bm{r},t)$ should have large values around the center of the domain wall because $D(\bm{r})$ rapidly changes over the edge as shown in Fig. 4(a). The typical length scale for this change is $\xi$, and the main contribution of the spatial integral to Eq. (LABEL:droplet4) is expected to come from $\mathcal{S}_{\rm wall}$ schematically shown in Fig. 4(b). Thus, we can replace the integral region $\mathcal{S}$ by $\mathcal{S}_{\rm wall}$:

[TABLE]

To analyze the integral (29), we evaluate the spin current $\bm{J}_{z}(\bm{r},t)$ of Eq. (22). We normalize it by the characteristic length $\xi$ and time $\tau$, and the current becomes

[TABLE]

with the non-dimensional gradient $\bar{\bm{\nabla}}=\xi\bm{\nabla}$ and the characteristic velocity $C_{\rm ch}=\xi/\tau$. To estimate the order of each term on the right-hand side of Eq. (30), we assume that the spatial gradients in Eq. (30) can be replaced by the size of the droplet $L_{\rm dp}$ because of the following reason. In the relaxation dynamics in Figs. 1(c) and (d), surfaces of droplets have various excitations whose wavelengths are smaller than $L_{\rm dp}$. However, we expect that such fine surface excitations are not essential in the droplet dynamics since the effect may be canceled in the integral of Eq. (29) over space. Thus, the main contribution to the integral is expected to come from the spatial structures of $f_{\nu}(\bm{r},t)$ and $v_{\nu}(\bm{r},t)$ on the scale of the droplet size $L_{\rm dp}$. By using this approximation, we can obtain the following naive condition to neglect the mass current term in Eq. (30):

[TABLE]

[TABLE]

where $\langle\cdots\rangle_{\rm wall}$ the spatial average inside the wall $\mathcal{S}_{\rm wall}$. The left (right)-hand side on inequality (31) comes from the first (second) term of the right-hand side on Eq. (30).

To obtain the values of $\Lambda$, we numerically calculate the absolute values of $f_{z}{\bm{v}}_{\nu}/C_{\rm ch}$ and $f_{\nu}$ and average them inside the walls. As shown in Fig. 5, $\sqrt{f_{x}^{2}+f_{y}^{2}}$ is about $0.953$ irrespective of $R$ while $|f_{z}\bm{v}/C_{\rm ch}|$ is about $0.07$, $0.08$, and $0.09$ for $R=0.25$, $0.35$, and $0.5$, respectively. Thus, our numerical simulations for $R=0.25$, $0.35$, and $0.5$ leads to $\Lambda=13.0$, $11.4$, and $10.1$, respectively. Combining all the results obtained so far, we derive an expression for the spin current depending on the droplet size $L_{\rm dp}$:

[TABLE]

Substituting Eq. (34) into Eq. (29), we obtain an expression for the centroid velocity $X_{\nu}(t)~{}(\nu=x,y)$ of the single droplet:

[TABLE]

This demonstrates that the main driving mechanisms for the droplet can change in accordance with the droplet size. As described in the next subsection, the change of the predominant term leads to the two power laws, and thus we can predict the crossover phenomena of the scaling in Fig. 3.

Finally, we estimate the $L_{\rm dp}$dependence of the integral in the right-hand side of Eq. (35). We notice that the spin density $f_{\nu}$ and the mass current $\bm{v}$ in $\mathcal{S}_{\rm wall}$ are expected to modulate on the scale of $L_{\rm dp}$, and thus the integral can be estimated as

[TABLE]

[TABLE]

Combining Eqs. (35), (36) and (37), we conclude

[TABLE]

Here, we use the relation $S_{\rm dp}\propto L_{\rm dp}^{2}$.

IV.2 Growth law of droplets

We construct the equation of motion for the correlation length $L(t)$ using Eq. (38), and derive the $1/3$ and $1/2$ power laws. Suppose that we have many droplets with the average size $L_{\rm dp}(t)$ and the mean distance between them is denoted by $L_{\rm md}(t)$. Then, the number $N_{\rm dp}(t)$ of the droplets is given by

[TABLE]

with a constant $A$ proportional to the system area $S$. Assuming that these droplets merge through random collisions between them, we expect that an average collision-time is about $L_{\rm md}(t)/V(L_{\rm dp}(t))$. Thus, the the simple rate equation for $N_{\rm dp}(t)$ is given by

[TABLE]

with a constant $B$. We note that $L_{\rm md}(t)$ is proportional to $L_{\rm dp}(t)$:

[TABLE]

This relation can be derived by considering the conservation laws for $N_{1}$ and $N_{2}$, the detail of which is described in the Appendix D. Substituting Eqs. (38), (39), and (41) into Eq. (40), we obtain

[TABLE]

Regarding that the droplet size $L_{\rm dp}(t)$ is proportional to the correlation length $L(t)$, we finally obtain the two power laws by solving Eq. (42):

[TABLE]

This shows that the correlation length monotonically grows in time, and thus we readily expect that the $1/3$ power law appears in the early stage of the dynamics and the $1/2$ power law gradually does in the late time.

Comparing Eq. (43) with the numerical results of Fig. 3, we notice that there is a small difference for the exponent in the early time, but the entire evolution of $L(t)$ is consistent with our theoretical prediction. To quantify the consistency, we compare the values of the relation between $L^{\prime}$ used in Fig, 3 and $\Lambda\xi$. The former is the rescaling factor of $L(t)$ which is defined such that the curves for $R=0.25$ and $0.35$ coincide with that for $R=0.5$. In other words, $L^{\prime}/\xi$’s for $R=0.25$ and $0.35$ give the ratio of the crossover lengths $\Lambda\xi$, at which the power law changes, to that for $R=0.5$. Thus, the values of $\Lambda\xi$ for $R=0.25$ and $0.35$ divided by that for $R=0.5$ can be close to $L^{\prime}/\xi$ because $\Lambda\xi$ is a measure of the crossover points. We summarize the numerical numbers in Table 1, which shows that the two quantities have the reasonable values despite the naive estimation of $\Lambda$ in Eq. (33).

V Discussion

We discuss four issues. Firstly, we consider stability of droplet structures, which is the assumption used in our analytical model of Sec. IV. Secondly, introducing previous works on relaxation dynamics with large imbalances in classical systems, we clarify the difference and argue that the mechanism of the power laws in our system is unconventional. Thirdly, the dependence of the relaxation dynamics on the imbalance parameter $R$ is discussed. Finally, we propose an experimental setup to observe our theoretical prediction such as crossover phenomena.

V.1 Stability of isolated droplets

Our analytical model described in Sec. IV has the important assumption that an isolated droplet is stable against fluctuations behind it. Here, we verify this by numerically investigating the isolated droplet dynamics.

V.1.1 Droplet dynamics in a realistic situation

We prepare an initial state with four isolated droplets using the imaginary time step of Eqs. (1) and (2) with $\Omega=0$. To stabilize the droplets with the sizes $L_{p}/\xi=8,10,12$, and $14$, we apply the following external potential $V_{p}$ only to the component $1$:

[TABLE]

where the bulk chemical potential $\mu$ is identical to that in Figs. 1(c) and (d). Here, the region $\mathcal{A}_{\rm dp}$ is given by $\mathcal{A}_{1}\cup\mathcal{A}_{2}\cup\mathcal{A}_{3}\cup\mathcal{A}_{4}$ with

[TABLE]

where the potential position $(x_{j},y_{j})$ and the radius $R_{j}$. In our simulation, we use $(x_{1}/\xi,y_{1}/\xi,R_{1}/\xi)=(64,64,4)$, $(x_{2}/\xi,y_{2}/\xi,R_{2}/\xi)=(192,64,5)$, $(x_{3}/\xi,y_{3}/\xi,R_{3}/\xi)=(64,192,6)$, and $(x_{4}/\xi,y_{4}/\xi,R_{4}/\xi)=(192,192,7)$. As a result, we obtain a macroscopic wavefunction $\psi_{m}^{\rm(dp)}(\bm{r})$ with the four droplets.

We add initial noises into $\psi_{m}^{\rm(dp)}(\bm{r})$ to mimic droplet motions in the quench dynamics of Figs. 1(c) and (d), and then numerically solve Eqs. (1) and (2) under $\Omega=0$ and $V_{p}=0$. The initial wavefunction used here is

[TABLE]

[TABLE]

where $\beta_{1,\bm{k}}$ and $\beta_{2,\bm{k}}$ are sampled from the probability distribution functions defined by

[TABLE]

[TABLE]

Here $B_{1}$ and $B_{2}$ are normalization constants, and widths $b_{1}$ and $b_{2}$ of the distributions are $2.7$ and $0.6$, which give almost same amplitude of the density fluctuations after the domain nucleation in Fig. 1. The form of these noises are identical to Eqs. (LABEL:GP5_1) and (5), except for the values of $b_{1}$ and $b_{2}$.

Figures 6(a)–(d) show time evolution for spatial distributions of $\rho_{2}(\bm{r},t)$ and $\rho(\bm{r},t)$, from which we find that all the isolated droplets are stable. We quantitatively estimate the stability by calculating a particle number contained in each droplet, which is measured by

[TABLE]

with $\mathcal{B}_{j}(t)=\{\bm{r}|\rho_{2}(\bm{r},t)>\rho_{\rm tot}/2\}\cap\{\bm{r}|x_{1,j}\leq x<x_{2,j},y_{1,j}\leq y<y_{2,j}\}$. Here, the coordinates used here are given by $(x_{1,1},x_{2,1},y_{1,1},y_{2,1})=(0,128\xi,0,128\xi)$, $(x_{1,2},x_{2,2},y_{1,2},y_{2,2})=(129\xi,255\xi,0,128\xi)$, $(x_{1,3},x_{2,3},y_{1,3},y_{2,3})=(0,128\xi,129\xi,255\xi)$, and $(x_{1,4},x_{2,4},y_{1,4},y_{2,4})=(129\xi,255\xi,129\xi,255\xi)$. As shown in Fig. 6(e), we find that the particle numbers in the droplets are approximately conserved in the long time dynamics. Here, we note the early stage of the dynamics, in which there are small decreases of the droplet particle numbers. We expect that this is caused by instabilities arising from our initial droplets, which are externally stabilized by applying the potential $V_{p}(\bm{r})$. Thus the initially prepared droplets are slightly different from the genuine configurations, and emit several partilces to reduce the sizes.

In the simulation of the quench dynamic described in Sec. III, we indeed confirm that the stable droplets collide each other, and the typical size grows in time as shown in Fig. 7.

V.1.2 Droplet stability mechanism

We now give an essential insight for the droplet stability, and argue that fluctuations of the component $2$ can be emitted from or absorbed into droplets, and that the in- and out-flow balance supported by the two processes makes droplets stable.

To investigate this droplet stability mechanism, we perform a numerical simulation starting from the following initial state without any noise terms in the component $2$:

[TABLE]

[TABLE]

Figures 8(a)–(d) shows time evolution of $\rho_{2}(\bm{r},t)$ and $\rho(\bm{r},t)$, and we find that the droplets become smaller in time. Correspondingly, the domain particle numbers within $\mathcal{B}_{j}$ exhibit clear decreases as in Fig. 8(e). Thus, the domains show diffusive behavior by emitting their particles.

Considering the above numerical results with the two different initial states [Eqs. (LABEL:noise2_1), (47), (LABEL:noise3_1), and (52)], we expect that the droplets in Fig. 6 are stabilized by balance between emission and absorption of the fluctuations of $\rho_{2}(\bm{r},t)$. This description is consistent with our results so far. In the simulation of Fig. 8, the component-$2$ particles outside the droplets do not exit, and thus the absorption process is almost absent before fluctuations emitted from other droplets spread over space. As a result, the droplets continue to emit their particles by the fluctuations of the component $1$, and the sizes decrease in time. On the other hand, in the simulation of Fig. 6, there exist the component-$2$ particles around the droplets as noises, both the emission and absorption processes can work well, and thus the isolated droplets become stable due to the balance. We also perform a numerical simulation with large noises in the component $2$, where we use $b_{2}=1.5$ larger than $0.6$ in Eq. (47) and other parameters are same as those in Fig. 6. If our argument is correct, the fluctuations of the component $2$ are excessively absorbed into the droplets, and their sizes should grow over time. As shown in Fig. 9, we in fact confirm the droplet growth, which supports the droplet stability mechanism mentioned above.

In classical dissipative systems, fluctuations around droplets vanish due to a diffusion process, and this kind of stability mechanism is not valid. However, because our isolated system conserves the total energy, the fluctuations survive for a long time and stabilize the droplets.

Finally, we comment on possible scenarios of droplet collapse. If strength of the fluctuations is strong enough, droplets can collapse. However, in our quench dynamics, this process may be rare because this occurs only when a droplet size is comparable to the coherence length $\xi$. In fact, as shown in Fig. 6, the droplet with $L_{p}=8\xi$, which is comparable to the smallest domain size in the quench dynamics, is still stable. As another mechanism, quantum fluctuations may break droplets. Several 1D spin system actually finds the domain melting Halimeh et al. (2014). However, compared to the 1D system, quantum fluctuations should not be dominant in our case because the system considered in this work is a weakly interacting BEC. Thus, we do not expect that such domain melting often occurs in 2D binary BECs, although our numerical method with the truncated Wigner approximation-like noises does not fully include the effect of quantum fluctuations.

V.2 $1/3$ power law in conventional classical systems and the relation with our results

Relaxation dynamics considered here has been investigated for a long time in terms of Ostwald ripening and off-critical quench in binary mixtures Bray (1994); Onuki (2002); Lifshitz and Slyozov (1961); Wagner (1961); Rogers and Desai (1989); Farrell and Valls (1991); Yao et al. (1993); Gorentsveig et al. (1997); Huse (1986), in which the correlation length has been found to obey the $1/3$ power law. The exponent is identical to our analytical result (43) in the early time. In this section, we argue that the fundamental mechanism behind the $1/3$ power in the binary BECs is essentially different from that in classical systems.

Let us briefly review a physical origin of the $1/3$ in classical systems. One standard explanation for Ostwald ripening is given by the Lifshitz-Slyozov-Wagner (LSW) theory Lifshitz and Slyozov (1961); Wagner (1961), which predicted that small droplets shrink and vanishes while large ones grow in time because of diffusion processes of concentration fields around the droplets. Based on this evaporation and concentration, the LSW theory analytically derives the $1/3$ power law, and many numerical simulations for several models have confirmed it Rogers and Desai (1989); Farrell and Valls (1991); Yao et al. (1993); Gorentsveig et al. (1997); Huse (1986).

However, the classical mechanism of the $1/3$ power law does not emerge in binary BECs because of the energy conversation in Eqs. (1) and (2). Although in the LSW theory diffusion processes of fields play a pivotal role in relaxation of classical mixtures, the conservation law basically prohibits such process in the atomic BEC. Our model, actually, does not have diffusive terms breaking the energy conservation. Due to the conservation law, a droplet itself is quite stable and cannot shrink in our relaxation dynamics. In fact, this behavior is numerically confirmed in Sec. V.1. As an exceptional process, too small droplets comparable to the coherence length $\xi$ can collapse by emitting the partilces when the background fluctuations are quite strong. However, this kind of collapse is rare as mentioned in Sec. V.1. Thus, we conclude that the $1/3$ power in the binary BECs is fundamentally different from the conventional LSW mechanism.

We comment on the previous works Binder and Stauffer (1974); Siggia (1979); Furukawa (1985a, b); San Miguel et al. (1985); Shimizu and Tanaka (2015), where random motion of droplets in classical systems is considered and various power exponents $\beta$ are analytically obtained. This analytical approach is similar to our calculation in Sec. IV, but most of the models used in previous literature have diffusive terms breaking the energy conservation and furthermore the droplets are not driven by the spatial gradient of spin and velocity. Thus, the fundamental relaxation mechanisms seem to be different from that in our systems. It is interesting to consider whether or not such droplet dynamics in classical systems still emerges in ultracold quantum gases, and it remains as a future work.

V.3 Dependence of the dynamics on $R$

We discuss how the relaxation dynamics changes depending on the particle-number imbalance parameter $R$. Our analytical model in Sec. IV is based on the stable droplet structures and their random collisions. It is obvious that such description should break down when $R$ is close to unity since large scale domains have surface excitations with various length-scales. However, our numerical simulations with $R\sim 1$ tend to show the power law behavior similar to Eq. (43). In this current work, we cannot sufficiently understand the behavior, and we leave it in the future work.

V.4 Possible experimental setups

We discuss experimental possibility for observing our theoretical predictions. As described in the previous section, the universal scale-invariant dynamics is induced by random collisions between tiny droplet structures, so that we need to prepare a situation in which a droplet size is much smaller than a system size just after the tuning protocol. A diameter of a system in several experiments Kang et al. (2019); Kim et al. (2017) is about $200~{}{\rm\mu m}$ at the most, so that an averaged droplet linear-size is required to be smaller than at least $20~{}{\rm\mu m}$. The typical size is roughly estimated as $\lambda_{\rm in}/2$ in Eq. (17). We plot its dependence on $g_{1}/g_{0}$ and $R$ in Fig. 10 using the typical coherence length $\xi=0.5~{}{\rm\mu m}$. We finds that the size can be smaller than $20~{}{\rm\mu m}$ if the ratio $g_{1}/g_{0}$ of the interaction couplings is larger than $1.1$.

We consider two internal states $|F=1,m=-1\rangle$ and $|F=2,m=-2\rangle$ of ${}^{87}{\rm Rb}$ as a candidate to realize the suitable system mentioned above because the binary condensate has a long lifetime. This kind of system has in fact been used to study long-time dynamics Yan et al. (2011). The interaction couplings, however, are almost the same and the ground state is a miscible state. When we can tune the interaction coupling between the different components by means of the Feshbach resonance, our theoretical prediction should be observed.

VI Conclusion

We have theoretically studied the relaxation dynamics in binary BECs with large particle-number imbalances, and uncovered the emergence of the scale-invariant relaxation dynamics. The significance of our finding is that the relaxation dynamics exhibits the dynamical scaling of the correlation function and the scaling law have two different power exponents depending on the elapsed time. We have analytically derived the $1/3$ and $1/2$ power laws for the growth of the correlation length, predicting the crossover in time and then numerically confirm it. As described in the derivation, the origin of the power laws is random collisions between droplets, whose centroids fluctuate with size-dependent velocities. The size dependence is determined by which is the main contribution to the spin current, the spatial variation of the spin density vector or the mass current. The former (latter) is dominant for a small (large) droplet size and leads to $1/3~{}(1/2)$ power law. Hence, the crossover from the $1/3$ to $1/2$ power law occurs as the droplet size grows. Our prediction can be experimentally realized by utilizing the spin polarized mixture of $F=1$ and $2$ hyperfine states of ${}^{87}{\rm Rb}$ atoms with a controlled inter-species interaction by means of the Feshbach resonance.

Acknowledgements.

This work was supported by JST-CREST (Grant No. JPMJCR16F2) and JSPS KAKENHI Grant Nos. JP15K17726, JP18K03538, JP19H01824, and JP19K14628.

Appendix A Imaginary part of the correlation function

As mentioned in Sec. II.4, the imaginary part of $C(r,t)$ is zero. Here, we give the proof in the following. The wavefunction can be expressed in the wave number space:

[TABLE]

with the Fourier component $\bar{\psi}_{2}({\bm{k}},t)$. Substituting Eq. (A-1) into Eq. (11), we obtain

[TABLE]

with the Bessel function $J_{n}(x)~{}(n=0,1,2,\cdots$) of the first kind. Here, we calculate the angle-averaged part as

[TABLE]

Thus, from Eq. (A-2), we find that the imaginary part of $C(r,t)$ becomes zero.

Appendix B Initial instability in the Bogoliubov theory

We derive the expression of Eq. (17) for the most unstable wavelength and describe the strength of the instability. Applying the conventional Bogoliubov theory to Eqs. (1) and (2) with $\Omega=0$ and a fixed particle-number imbalance $R$, we derive the Bogoliubov excitation energy Pethick and Smith (2008):

[TABLE]

where

[TABLE]

[TABLE]

with $\rho_{1}=\rho_{\rm tot}/(1+R)$ and $\rho_{2}=\rho_{\rm tot}R/(1+R)$. The wavenumber $k_{\rm in}$ determined by $\epsilon(k_{\rm in})=A/2$ gives the largest imaginary part of $E(k)$ due to the quadratic form with respect to $\epsilon(k)$ in Eq. (B-1) when the phase-separation condition $g_{0}<g_{1}$ is satisfied. The analytical expression for $k_{\rm in}$ becomes

[TABLE]

which leads to Eq. (17) through $\lambda_{\rm in}=2\pi/k_{\rm in}$.

Next, we give an expression for the most largest imaginary part of $E(k)$. Substituting Eq. (B-4) into Eq. (B-1), we obtain

[TABLE]

When the parameters $R$, $g_{0}$, and $g_{1}$ satisfy $g_{0}<g_{1}$ and $0<R<1$, the amplitude $|{{\rm Im}E(k_{\rm in})}|$ is a monotonically increasing function of $R$. This readily confirmed by

[TABLE]

Thus, we find that the initial instability becomes stronger when $R$ is larger. As described in the caption of Fig. 5(a), our numerical simulations are consistent with this dependence.

Appendix C Derivation of the continuity equation for $f_{z}$

This appendix derives Eqs. (21) and (22) from the GP equations (1) and (2) with $\Omega=0$. The time derivative of $f_{z}$ in the GP model readily leads to

[TABLE]

[TABLE]

To rewrite Eq. (C-2) in terms of the spin density vector $f_{\nu}$, we use the following notation for the macroscopic wavefunctions:

[TABLE]

where $\Phi$ is the overall phase and $\theta$ and $\phi$ are the polar and azimuthal angles of the spin density vector, i.e., $\bm{f}$ is given by

[TABLE]

Substituting Eq. (C-3) into Eqs. (C-2) and (23), we obtaion

[TABLE]

By calculating the spatial gradient of $f_{x}$ and $f_{y}$, the following equality is obtained:

[TABLE]

We substitute Eqs. (C-8) and (C-9) into Eq. (C-7), and then obtain

[TABLE]

Appendix D Relation between $L_{\rm dp}(t)$ and $L_{\rm md}(t)$

In Sec. IV.2, we have introduced two length scales $L_{\rm md}(t)$ and $L_{\rm dp}(t)$ as the mean-droplet distance and the typical linear size of droplets, respectively. In this Appendix, we show that they are related to each other due to the particle number conservation. Suppose that we have $N_{\rm dp}(t)$ droplets with the typical linear size $L_{\rm dp}(t)$. Then, the imbalance parameter $R$ is estimated as

[TABLE]

with a constant $\eta\sim\mathcal{O}(1)$. The mean distance $L_{\rm md}(t)$ between the droplets can be expressed by

[TABLE]

with a constant $\zeta\sim\mathcal{O}(1)$. Then, combining the above equations, we obtain

[TABLE]