Effective rheology of two-phase flow in a capillary fiber bundle model

TL;DR

This paper analyzes how two immiscible fluids flow through a bundle of capillary tubes with varying diameters, revealing a transition from linear to non-linear flow regimes and deriving scaling laws for flow rate near critical thresholds.

Contribution

The study provides an analytical and numerical investigation of the transition from linear to non-linear flow in a capillary bundle model with random threshold distributions, including new scaling laws.

Findings

Flow rate transitions from linear to non-linear with an exponent of 2.

Introduction of a lower cutoff in threshold distribution alters the non-linear exponent.

Flow rate scales as (|ΔP|-Pm)^{3/2} near the threshold Pm.

Abstract

We investigate the effective rheology of two-phase flow in a bundle of parallel capillary tubes carrying two immiscible fluids under an external pressure drop. The diameter of each tube varies along its length and the corresponding capillary threshold pressures are considered to be distributed randomly according to a uniform probability distribution. We demonstrate through analytical calculations that a transition from a linear Darcy regime to a non-linear behavior occurs while decreasing the pressure drop $\Delta P$, where the total flow rate $\langle Q \rangle$ varies with $\Delta P$ with an exponent $2$. This exponent for the non-linear regime changes when a lower cut-off $P_m$ is introduced in the threshold distribution. We demonstrate analytically that, in the limit where $\Delta P$ approaches $P_m$, the flow rate scales as $\langle Q \rangle \sim (|\Delta P|-P_m)^{3/2}$. We have also provided some numerical results in support to our analytical findings.

Equations

$$A_p=Na.$$

$$q=-\frac{a^2}{8\pi {\mu }_{\mathrm{av}}L}\Theta (|\Delta P|-P_c)\left[\Delta P-P_c\right],$$

$${\mu }_{\mathrm{av}}=S_w{\mu }_w+S_n{\mu }_n,$$

$$\overline{q}(P_c)=$$

$$-\frac{a^2}{8\pi {\mu }_{\mathrm{av}}L}\mathrm{sgn}(\Delta P)\Theta (|\Delta P|-P_c)\sqrt{\Delta P^2-P_c^2},$$

$$\Pi(P_{c})=\left\{\begin{array}[]{ll}0&\mbox{, $P_{c}\leq 0$\;,}\\ \frac{P_{c}}{P_{M}}&\mbox{, $0<P_{c}\leq P_{M}$\>,}\\ 1&\mbox{, $P_{c}>P_{M}$\;.}\\ \end{array}\right.$$

$$\Pi (P_c(k))=\frac{k}{N+1},$$

$$P_c(k)=P_M\frac{k}{N+1}.$$

$$⟨Q⟩=\sum _{k=1}^{(N+1)\min \left(\frac{|\Delta P|}{P_M},1\right)}\overline{q}\left(P_c(k)\right)$$

$$\frac{⟨Q⟩}{N}$$

$$\frac{⟨Q⟩}{N}=-\frac{a^2}{32{\mu }_{\mathrm{av}}L}\left|\frac{\Delta P}{P_M}\right|\Delta P$$

$$\frac{⟨Q⟩}{N}=-\frac{a^2{P}_M\mathrm{sgn}(\Delta P)}{16\pi {\mu }_{\mathrm{av}}L}$$

$$[\sqrt{{\left(\frac{|\Delta P|}{P_M}\right)}^2-1}$$

$$+{\left(\frac{|\Delta P|}{P_M}\right)}^2\arcsin \left(\frac{P_M}{|\Delta P|}\right)],$$

$$\frac{⟨Q⟩}{N}=-\frac{a^2}{8\pi {\mu }_{\mathrm{av}}L}\Delta P.$$

$$k_c=\frac{|\Delta P|}{P_M}(N+1).$$

$$l_{\perp }=\frac{N+1}{k_c}=\frac{P_M}{|\Delta P|},$$

$$\Pi(P_{c})=\left\{\begin{array}[]{ll}0&\mbox{, $P_{c}\leq 0$\;,}\\ \left(\frac{P_{c}}{P_{M}}\right)^{\alpha}&\mbox{, $0<P_{c}\leq P_{M}$\>,}\\ 1&\mbox{, $P_{c}>P_{M}$\;,}\\ \end{array}\right.$$

$$P_c(k)=P_M{\left(\frac{k}{N+1}\right)}^{1/\alpha },$$

$$\frac{⟨Q⟩}{N}$$

$$\frac{⟨Q⟩}{N}=-\frac{a^2\alpha }{32\sqrt{\pi }{\mu }_{\mathrm{av}}L}\frac{\Gamma \left(\frac{\alpha }{2}\right)}{\Gamma \left(\frac{3+\alpha }{2}\right)}{\left(\frac{|\Delta P|}{P_M}\right)}^{\alpha }\Delta P,$$

$$\Pi(P_{c})=\left\{\begin{array}[]{ll}0&\mbox{, $P_{c}\leq P_{m}$\;,}\\ \frac{P_{c}-P_{m}}{P_{M}-P_{m}}&\mbox{, $P_{m}<P_{c}\leq P_{M}$\>,}\\ 1&\mbox{, $P_{c}>P_{M}$\;.}\\ \end{array}\right.$$

$$P_c(k)=P_m+(P_M-P_m)\frac{k}{N+1}.$$

$$\frac{⟨Q⟩}{N}$$

$$\frac{⟨Q⟩}{N}=-\frac{a^2\mathrm{sgn}(\Delta P)}{3\sqrt{2}\pi {\mu }_{\mathrm{av}}L}\frac{\sqrt{P_m}}{(P_M-P_m)}(|\Delta P|-P_m{)}^{3/2},$$

$$\Pi(P_{c})=\left\{\begin{array}[]{ll}0&\mbox{, $P_{c}\leq P_{m}$\;,}\\ \left(\frac{P_{c}-P_{m}}{P_{M}-P_{m}}\right)^{\alpha}&\mbox{, $P_{m}<P_{c}\leq P_{M}$\>,}\\ 1&\mbox{, $P_{c}>P_{M}$\;,}\\ \end{array}\right.$$

$$\Pi(P_{c})=\left\{\begin{array}[]{ll}0&\mbox{, $P_{c}\leq P_{m}$\;,}\\ \frac{\log\left(\frac{P_{c}}{P_{m}}\right)}{\log\left(\frac{P_{c}}{P_{M}}\right)}&\mbox{, $P_{m}<P_{c}\leq P_{M}$\>,}\\ 1&\mbox{, $P_{c}>P_{M}$\;,}\\ \end{array}\right.$$

$$\Pi(P_{c})=\left\{\begin{array}[]{ll}0&\mbox{, $P_{c}\leq P_{m}$\;,}\\ 1-e^{-P_{c}-P_{m}/P_{d}}&\mbox{, $P_{m}<P_{c}$\>,}\\ \end{array}\right.$$

Full text

Understanding the hydrodynamic properties of simultaneous flow of two or more immiscible fluids is essential due its relevance to a wide variety of different systems in industrial, geophysical and medical sectors b72 ; d92 . Different applications, such as bubble generation in microfluidics, blood flow in capillary vessels, catalyst supports used in the automotive industry, transport in fuel cells, oil recovery, ground water management and ${\rm CO}_{2}$ sequestration, deal with the flow of bubble trains in different types of systems, ranging from single capillaries to more complex porous media. The underlying physical mechanisms in multiphase flow are controlled by a number of factors, such as the capillary forces at the interfaces, viscosity contrast between the fluids, wettability and geometry of the system, which make the flow properties different compared to single phase flow. When one immiscible fluid invades a porous medium filled with another fluid, different types of transient flow patterns, namely viscous fingering cw85 ; mfj85 , stable displacement ltz88 and capillary fingering lz85 are observed while tuning the physical parameters lmtsm04 . These transient flow patterns were modeled by invasion percolation ww83 and diffusion limited aggregation (DLA) models ws81 . When steady state sets in after the initial instabilities, the flow properties in are characterized by relations between the global quantities such as flow rate, pressure drop and fluid saturation v18 ; hsbkgv18 . It has been observed theoretically and experimentally that, in the regime where capillary forces compete with the viscous forces, the two-phase flow rate of Newtonian fluids in the steady state no longer obeys the linear Darcy law d56 ; w86 but varies as a power law with the applied pressure drop tkrlmtf09 ; tlkrfm09 ; rcs11 ; sbdk17 . Tallakstad et al. tkrlmtf09 ; tlkrfm09 experimentally measured the exponent of the power law to be close to two ($=1/0.54$) in a two-dimensional system and followed this observation up with arguments why the exponent should be two. Rassi et al. rcs11 found a value for the exponent varying between $2.2$ ($=1/0.45$) and $3.0$ ($=1/0.33$) in a three-dimensional system. Sinha et al. sbdk17 considered a similar system to that which had been studied by Rassi et al. finding an exponent $2.17\pm 0.24$ ($=1/(0.46\pm 0.05)$). The reason behind the discrepancy between the results of Rassi et al. and those of Sinha et al. is the possibility of a non-zero threshold pressure that observed in the later study, under which there would be no flow, which was assumed to be zero in the former study. The reciprocals in the brackets are provided in order to compare the exponent values reported in the literature tkrlmtf09 ; tlkrfm09 ; rcs11 ; sbdk17 with those we present here in this article, as we express our results as $\langle Q\rangle$ as a power law in $\Delta P$, whereas the cited papers write $\Delta P$ as a power law in $\langle Q\rangle$.

This power law behavior is in contrast to the assumption of linearity in the relation between flow rate and pressure drop that is generally assumed in the relative permeability approach dominating reservoir simulations wb36 .

For a single capillary tube with varying diameter, Sinha et al. shbk13 showed that the average volumetric flow rate $\overline{q}$ in the steady state has a non-linear square-root type relationship with the pressure drop $\Delta P$ as $\overline{q}\sim\sqrt{\Delta P^{2}-P_{c}^{2}}$. This was shown analytically by integrating the instantaneous linear two-phase flow equation over the whole capillary tube. Here $P_{c}$ is the threshold pressure difference below which there is no flow. It appears due to the capillary barriers at the interfaces at the narrow pore throats. Extending this non-linear relationship to a network of disordered pores, the relationship between the steady-state flow rate and an excess pressure drop leads to a quadratic relationship in the capillary dominated regime sh12 . The quadratic relationship for the pore network, both in two and in three dimensions, was obtained analytically by mean-field calculations and numerically with pore network modeling sbdk17 ; sh12 .

While increasing the pressure drop, the capillary forces become negligible compared to the viscous forces. This leads to a crossover from the non-linear regime to a linear Darcy regime for both the single capillary and for the pore network. Such non-linear quadratic relationship at low flow rate and a crossover to a linear regime at high flow rate was also observed in case of the single-phase flow of Bingham viscoplastic fluid in porous media rh87 ; ct15 . A Bingham fluid is a yield threshold fluid which behaves like a solid below the threshold and flows like a Newtonian fluid above it. The origin of the quadratic regime for the Bingham fluid flowing in a porous media can be understood intuitively in this way: the flow starts when one connected channel appears in the system just above a threshold pressure and the flow rate varies linearly with the excess pressure drop; while increasing the applied pressure drop further, more number of connected flow channels start to appear enhancing the overall flow rate more rapidly than the applied pressure drop leading to the quadratic relationship. Finally, when all possible flow paths become active, the flow become Newtonian following the linear Darcy law. Note that, in general, the rheology of the Bingham fluid is linear above the yield threshold. It is the disorder in the yield thresholds due to the porous medium that creates the quadratic regime.

The argument presented by Tallakstad et al. tkrlmtf09 ; tlkrfm09 focused on the successive opening of fluid channels when the pressure drop across the system was increased. When $|\Delta P|$ is small, the flow will occur along isolated channels. The volumetric flow rate in such a channel will be proportional to $|\Delta P|/L$. Between the channels there will be fluid clusters held in place by capillary forces, say of the order $p_{t}$. There is a pressure gradient $|\Delta P|/L$ in the flow direction. A given cluster of length $l_{\parallel}$ will be stuck if $p_{t}>l_{\parallel}|\Delta P|/L$. The largest stuck cluster will then have a size $l_{\parallel,\max}=Lp_{t}/|\Delta P|$. If we now assume that this length, $l_{\parallel,\max}$ is same as the distance between the channels where there is flow, $l_{\perp}$, then the total flow rate must be equal to the number of channels, which is proportional to $1/l_{\perp}$, multiplied by the flow rate in each channel. Hence, we have $Q\propto(1/l_{\perp})\ |\Delta P|\propto|\Delta P|^{2}$. Though this argument provides the same behavior as the one based on the mean field calculation sh12 for two-dimensional networks, a difference appears in three dimensions. Following the same arguments, it leads to the flow rate varying as the pressure drop to the third power as long as the isolated channels remain one-dimensional strings rather than two-dimensional sheets in three dimensions.

We present in this article a capillary fiber bundle model s53 ; s74 , which is a system of $N$ parallel capillary tubes, disconnected from each other, each carrying an independent bubble trail of two immiscible fluids under an external pressure drop. In a porous medium, a typical pore constitutes of two wide pore bodies at the ends and a narrow pore throat in the middle. When an interface moves along the pore, the capillary pressure associated with the interface becomes position dependent due to the change in the radius of curvature. This introduces an overall threshold pressure that depends on the position of all the interfaces shbk13 . One can simplify such shape of the pore by a sinusoidal type and a long capillary tube with varying radius can be seen as a series of many pores. In the capillary bundle model, the diameter of each tube varies along the axis identically and the disorder in the threshold appear due to the different interface positions in different tubes. This model is essentially the only model for immiscible two-phase flow which is analytically tractable. We calculate the total average flow rate as a function of the applied pressure drop and study the effect of disorder in the threshold distribution. We point out that, here we do not address the question of the relation between the fluid distributions in the capillaries and the respective threshold distributions. Our aim with this model is to investigate how the range of the disorder in the threshold distribution controls the effective flow properties. This provides insight into the non-linearities in steady-state two-phase flow. We will see that the exponent for the non-linear regime depends on the lower cut-off of the threshold distribution as well as on the behavior of the distribution near the cut-off. The possibility to study analytically for this model how the competition between viscous and capillary forces renders the Darcy relation non-linear, is a new and useful discovery.

The capillary fiber bundle model is a hydrodynamic analog of the fiber bundle model used in fracture mechanics to model mechanical failure under stress hhp15 . The fiber bundle model is an ideal example of a disordered system in statistical mechanics driven by threshold activated dynamics. It is a simple, yet very rich model to understand failure events in mechanical systems. In its simplest form it is analytically tractable. In more complex versions of the model, analytical calculations go hand in hand with numerical simulations.

Figure 1 illustrates a bundle with $N=5$ parallel capillaries. Each capillary tube has a length $L$ and an average inner area $a$. For each capillary, the diameter varies along the long axis identically. Each capillary is filled with a bubble train of wetting and non-wetting fluids. Due to the varying diameter, the capillary forces at the interfaces vary as the bubble train moves along the tubes. We assume that the wetting fluid does not form films along the pore walls so that the fluids do not pass each other. The lengths of the wetting and non-wetting fluids in each tube is $L_{w}$ and $L_{n}$ respectively such that the volume of the wetting fluid in each tube is $L_{w}a$ and the volume of the non-wetting fluid is $L_{n}a$, where $a$ is the average cross-sectional area of the capillary tubes. Hence, the saturations are $S_{w}=L_{w}/L$ and $S_{n}=L_{n}/L$ for each capillary tube. The cross-sectional pore area of the capillary fiber bundle is

[TABLE]

Though each tube contains the same amount of each fluid, it has its own division of the fluids into bubbles. We average over the ensemble of capillary tubes which constitute the capillary fiber bundle by averaging over the fluids in each tube that passes at a given instance through the imaginary cut shown in the figure. We will obtain the same averages if we consider a single capillary tube, averaging over a time interval the fluid passing the imaginary cut shbk13 ; sshb16 . This shows that the model is ergodic.

The volumetric flow rate in a capillary tube is given by shbk13

[TABLE]

where $\Delta P$ is the pressure drop across the capillary tube, $P_{c}$ the sum of all the capillary forces along the capillary tube due to the interfaces and

[TABLE]

is the effective viscosity. $\Theta(|\Delta P|-P_{c})$ is the Heaviside function which is zero for negative arguments and one for positive arguments. The function ${\rm sgn}(\Delta P)$ is the sign of the argument.

Sinha et al. shbk13 showed that the time average when the pressure difference across the tube is kept fixed is given by

[TABLE]

Suppose now that the thresholds $P_{c}$ are distributed uniformly between zero and a maximum value $P_{M}$. The cumulative threshold probability is then

[TABLE]

We have $N$ capillary tubes. Using order statistics, we may order the $N$ averaged threshold values,

[TABLE]

where $1\leq k\leq N$. Hence,

[TABLE]

The average volumetric flow rate through the capillary fiber bundle for $|\Delta P|>0$ is then

[TABLE]

We assume the limit $N\to\infty$ turning the sum into an integral,

[TABLE]

This integral is doable and we find

[TABLE]

when $|\Delta P|\leq P_{M}$ and

[TABLE]

when $|\Delta P|>P_{M}$. In the limit $|\Delta P|\gg P_{M}$, Equation (Effective rheology of two-phase flow in a capillary fiber bundle model) gives

[TABLE]

Hence, the Darcy relation for a tube is recovered.

We see that this picture is consistent with that central to the arguments of Tallakstad et al. tkrlmtf09 ; tlkrfm09 leading to the quadratic dependence of $Q$ on $\Delta P$. From Equation (7) we deduce that a number $k_{c}$ of the capillary tubes are active, where

[TABLE]

The typical distance between active capillary tubes in units of the distance between the tubes is then given by

[TABLE]

in accordance with the argument of Tallakstad et al.

How stable is the square law $Q\propto|\Delta P|^{2}$? That is, how much does it hinge on the choice of cumulative threshold probability $\Pi(P_{c})$. So far we have only considered the one given in Equation (5). Let us now generalize it to

[TABLE]

where $\alpha>0$. The average ordered threshold are then given by

[TABLE]

and when combined with the expression for $\langle Q\rangle$, Equation (8) in the limit $N\to\infty$, we find

[TABLE]

Since we are interested in the behavior for $|\Delta P|\to 0$, we do this integral under the assumption that $|\Delta P|<P_{M}$ finding

[TABLE]

where the $\Gamma$ function for real positive $z$ is defined as, $\Gamma(z)=\displaystyle\int_{-\infty}^{\infty}t^{z-1}e^{-t}dt$. When $\alpha=1$, we recover Equation (10).

Equation (10) shows the behavior observed experimentally in References tkrlmtf09 and tlkrfm09 . With Equation (18), we have just shown that $\langle Q\rangle/N\sim|\Delta P|^{\gamma}$ as $|\Delta P|\to 0$, where $\gamma$ depends on the threshold distribution, i.e. on $\alpha$ in Equation (15). Does this imply that there is no universality; that the experimentally observed behavior is due to the presence of a very specific threshold distribution?

As we now argue, there is universality. We note that the threshold distribution $p(P_{c})=d\Pi(P_{c})/dP_{c}$ behaves as $p(P_{c})\propto P_{c}^{\alpha-1}$. Hence, if $\alpha>1$, the distribution vanishes as $P_{c}\to 0$, whereas it diverges for $\alpha<1$. Hence, the behavior of the distribution is vastly different for these two cases, and this causes $\gamma$ to depend on $\alpha$. However, for $\alpha=1$, the distribution reaches a constant, non-zero value for $P_{c}\to 0$. Any threshold distribution with this behavior for small $P_{c}$, i.e., $p(P_{c})$ reaching a non-zero value and $dp(P_{c})/dP_{c}\to 0$ in the limit $P_{c}\to 0$ will give rise to the square power law seen in Equation (10). Such distributions are ubiquitous, and $\gamma=2$ is universal over this class of distributions.

We now consider $\alpha=1$ again, but introduce a minimum threshold $P_{m}$ so that the cumulative threshold probability is given by

[TABLE]

Equation (6) yields in this case the ordered threshold sequence

[TABLE]

Equation (8) now becomes in the limit $N\to\infty$

[TABLE]

when we assume $P_{m}\leq|\Delta P|\leq P_{M}$. We find

[TABLE]

We find to lowest order in $(|\Delta P|-P_{m})$, that this expression behaves as

[TABLE]

as $|\Delta P|\to P_{m}$.

We now turn to numerical simulations and observe that the numerical results are in good a agreement with the analytical findings. The numerical simulations also allow us to explore the regions which are analytically challenging. Results are shown in Figure 2 for a bundle containing $N=10^{5}$ capillary tubes and averaged over $10^{4}$ configurations. In Figure 2(a), we show the behavior of the volumetric flow rate $\langle Q\rangle$ as a function of increasing pressure drop $\Delta P$ for uniform threshold distributions with $P_{m}=0$ and $P_{m}>0$, given by Equations (5) and (19) respectively. The results show that, for each threshold distribution, the relationship is linear for high $\Delta P$ obeying the Darcy law as predicted by Equation (12). For small pressure drops, $\langle Q\rangle$ follows a power law in $\Delta P$ with an exponent $2$ when there is no lower cut-off in the threshold distribution, i.e. $P_{m}=0$. This is predicted in Equation (10). When a lower cut-off is introduced in the threshold distribution ($P_{m}>0$), this exponent shifts from $2$ to $3/2$ as predicted in Equation (23). These exponents in the non-linear regime are not sensitive to the span of the distribution as shown in Figure 2 (a). An insight to a more generalized picture is presented in Figure 2(b) and (c) for a threshold distribution given by a generalization of Equation (15) with an introduction of a lower cut-off $P_{m}$,

[TABLE]

With this distribution of thresholds, the exponent $\gamma$ in the non-linear region shows a continuous variation with $\alpha$ as $\gamma=\alpha+1$ for $P_{m}=0$. Such variation is given in Equation (18) and matches well with the numerical results as shown in see Figure 2 (d). In presence of a lower cut-off $P_{m}>0$, $\gamma$ varies as $(\alpha+1/2)$ instead of $(\alpha+1)$ irrespective of the position of the lower cut-off. An analytical treatment for a general $\alpha$ value with $P_{m}>0$ is rather challenging. Nevertheless, our numerical result matches with the analytical study (see Equation 23) in the limit $\alpha=1$.

Equation (18) predicts an exponent $\gamma=\alpha+1$. A simple argument, related to that given by Roux and Herrmann rh87 , goes as follows: The number of active capillary tubes is proportional to $(|\Delta P|-P_{m})^{\alpha}$. This behavior is observed in the insets in Figures 2 (b) and (c). The flow rate in an active capillary is proportional to $(|\Delta P|-P_{m})^{1/2}$. Hence, the total flow rate should be $\langle Q\rangle\propto(|\Delta P|-P_{m})^{\alpha+1/2}$. It is accidental that this argument works out for $P_{m}>0$ (Figure 2 (c)), as it does not when $P_{m}=0$, where $\gamma=\alpha+1$. For the argument to function, the distribution of active capillaries and the flow rate in each capillary should be uncorrelated. They are not.

We find the same behavior with respect to the cut-off: An exponent 3/2 for the cumulative threshold probability

[TABLE]

where $P_{m}=10^{-\beta}$ and $P_{M}=10^{\beta}$ and $\beta$ ranging from $0.5$ to $1.5$. The same goes for the cumulative threshold probability

[TABLE]

where we have set $P_{m}=0.1$ and $P_{d}=1$. In both of these cases, the probability density at $P_{c}=P_{m}$ is finite.

We have presented an analytical study supported by numerical simulations of steady-state two-phase flow in a system of parallel capillary tubes. Considering a uniform distribution for the threshold pressures for the capillaries, we have calculated the average flow rate as a function of the applied pressure drop. When the thresholds are distributed according to a uniform distribution between zero and a maximum value — or more generally, the threshold distribution approaches a non-zero value in the limit of zero thresholds — we obtain a quadratic relationship between the flow rate and the applied pressure drop when the applied pressure drop is below the maximum threshold pressure, and the linear Darcy relationship for higher pressure drops. This crossover between a quadratic non-linear and linear flow regimes is in agreement with many existing results of two-phase flow in porous media which shows that this simple model can capture effective two-phase flow properties of more complex porous media. When a lower cut-off is introduced in the threshold distribution, the quadratic relationship changes, and the flow rate varies with an excess pressure drop with an exponent $3/2$ as the pressure drop approaches to the lowest threshold pressure.

The difference between the capillary fiber bundle model and a porous medium is that in the latter, the fluids meet and mix at the nodes of the pore network. This is an essential mechanism that leads to the non-linear Darcy law is a power law with an exponent two as seen in the experiments, the numerical simulations and the mean-field calculations. However, it remains a mystery how the mixing at the nodes leads to this universality.

The authors thank Dick Bedeaux, Carl Fredrik Berg, Eirik G. Flekkøy, Signe Kjelstrup, Knut Jørgen Måløy, Per Arne Slotte and Ole Torsæter for interesting discussions. This work was partly supported by the Research Council of Norway through its Centres of Excellence funding scheme, project number 262644. SS was supported by the National Natural Science Foundation of China under grant number 11750110430.