4.1 Invariant solutions in an electrified film

This section examines the oscillatory phenomenon in an electrified film. Previous numerical simulations have shown that an oscillation between big droplets can occur when subjected to an electric field (Wray et al. Reference Wray, Papageorgiou and Matar2013a ; Ding et al. Reference Ding, Xie, Wong and Liu2014). Stability analysis of a uniform film, $S=1$ , has shown that for $\unicode[STIX]{x1D6FD}<\exp (1)$ the electric field can enhance the linear instability. For $\unicode[STIX]{x1D6FD}>\exp (1)$ , however, the electric field impedes the linear instability, but is destabilizing with respect to finite-amplitude perturbations. The flow then proceeds to the formation of an oscillatory state or spike-like singular state, dependent on the value of $\unicode[STIX]{x1D6FD}$ (Ding et al. Reference Ding, Xie, Wong and Liu2014). In this paper, we are interested in the oscillatory state and therefore fix the radius of the outer electrode at $\unicode[STIX]{x1D6FD}=\exp (1.1)$ . In the following discussion, we fix the small parameter $\unicode[STIX]{x1D716}$ , and the radius of the fibre $\unicode[STIX]{x1D6FC}$ at $(\unicode[STIX]{x1D716},\unicode[STIX]{x1D6FC})=(0.2,0.25)$ following previous studies, (e.g. Craster & Matar Reference Craster and Matar2006; Liu et al. Reference Liu, Chen and Ding2018). To understand the oscillation dynamics between big droplets, we consider a periodic computational domain that is large enough for two big droplets. To determine an appropriate domain size $L$ , we refer to the linear stability analysis by Ding et al. (Reference Ding, Xie, Wong and Liu2014), where the wavelength of the most unstable mode in the electrified film is

(4.1) $$\begin{eqnarray}\frac{2\sqrt{2}\,\unicode[STIX]{x1D716}\unicode[STIX]{x03C0}}{\displaystyle \sqrt{1+\frac{E_{b}(1-\ln (\unicode[STIX]{x1D6FD}))}{\ln ^{3}(\unicode[STIX]{x1D6FD})}}}\sim 2.\end{eqnarray}$$

We therefore start with domain sizes of $L\sim 4$ for $0\leqslant E_{b}\leqslant 4$ to allow two big droplets. Figure 3 shows the type of solution observed in simulations for the range $4\leqslant L\leqslant 5$ . For $L>6$ , we expect three or more big droplets to form in the domain.

Figure 3.

Phase diagram of solution regimes in the $L$ versus $E_{b}$ plane. The boundaries separating different solution regimes are obtained starting from RPOs and tracking where they undergo transition to spikes or TWs.

Figure 4.

(a) The profiles of the ‘two different/identical droplets’ (type I/type II) illustrated by ‘blue/red’ lines for $E_{b}=1$ , $L=4$ . (b–d) The wave speed $c$ versus the electric number $E_{b}$ for $L=4$ , 4.5 and 5, respectively. The type-I TW branch for $L=5$ stops at $E_{b}=2.76$ after it turns back, as the singular structure appears. Curves for TWs are obtained by a direct Newton–Raphson iteration using the code of Ding et al. (Reference Ding, Xie, Wong and Liu2014). Curves for RPOs in (b–d) are obtained from the Newton–Krylov method implemented in this work. Stable/unstable TW branches are solid/dashed lines. The Newton–Krylov method is validated by reproducing the TW solution as marked by diamonds in (b). The bullets indicate where the type-II TW connects with type-I TW.

Initially, to find stable states in simulations, we perturb the trivial solution $S=1$ by harmonic waves

(4.2) $$\begin{eqnarray}S(z,0)=1+0.1\cos \left(n\frac{2\unicode[STIX]{x03C0}}{L}z\right),\end{eqnarray}$$

with $n=1$ or $n=2$ and with $E_{b}$ small. For $L=4$ , we find that the system evolves into one of two stable TW states. For $n=2$ , two identical droplets form, hereafter called type I. Their spatial period is $L/2$ , but they are stable to small disturbances of $n=1$ type. Starting with $n=1$ in (4.2), two slightly different droplets form, hereafter referred to as type II. The two forms of TW solution are illustrated in figure 4(a). Starting from these non-oscillatory solutions at $L=4$ , $E_{b}=0$ , pseudo-arclength continuation with the code of Ding et al. (Reference Ding, Xie, Wong and Liu2014) has been used to track the TW solutions and the wave speed $c$ against the electric Weber number $E_{b}$ , shown in figure 4(b–d). For the case of type-II TW, the $c$ versus $E_{b}$ curve is shown in figure 4(b), which folds back at $E_{b}\approx 1.94$ for $L=4$ , while the curve stops at around $E_{b}=4.3$ for $L=4.5$ . Thus, for both $L=4$ and $L=4.5$ this indicates that the type-II TWs do not exist when the electric field exceeds some critical value. For $L=5$ neither type of solution exists for $E_{b}\gtrsim 4$ .

Figure 5.

(a) Evolution of the energy norm $\mathscr{E}=\langle S^{2}/2\rangle$ at different electric number subject to the initial condition $S(z,0)=1+0.1\cos (\unicode[STIX]{x03C0}z/2)$ where $z\in [0,4]$ . For $E_{b}=0,1.9$ , steady type-II TWs are obtained (see that $\mathscr{E}$ becomes invariant). (b) The spike-like structure plotted at $t\approx 8.5$ at $E_{b}=4$ , $L=4$ .

Numerical simulations of TW solutions show that they lose stability to an oscillatory bound state when $E_{b}$ exceeds a critical value (see figure 5 a). The oscillation suggests the existence of RPOs, which will be discussed presently. For larger $E_{b}$ , we observe that two droplets start to merge, and the minimal gap between the liquid interface and outer electrode reduces quickly – Ding et al. (Reference Ding, Xie, Wong and Liu2014) showed that the outer electrode attracts the interface, and the liquid interface becomes unstable. It is quickly attracted to the outer electrode when the gap is small, leading to the formation of the spike-like structure. (An example of the spike-like structure can be found in figure 5 b.) It is natural to examine whether the oscillatory solutions are linked to the TW solutions through oscillatory bifurcations.

The two leading eigenvalues (either two real or a complex pair) of the type-I TWs are plotted against $E_{b}$ in figure 6 for (a,b) $L=4$ , (c,d) $L=4.5$ and (e,f) $L=5$ . The eigenvalues for the type-II TW for these parameters are shown in figure 7.

Figure 6.

Eigenvalue of the type-I TWs versus the electric number $E_{b}$ at (a,b)  $L=4$ , (c,d)  $L=4.5$ and (e,f)  $L=5$ . The red solid (blue dashed) lines represent the stationary (oscillatory) mode. For $L=5$ , the stability of the lower branch of type-I TW is investigated in (e,f).

Figure 7.

(a) Leading eigenvalue for type-II TWs versus the electric number $E_{b}$ at $L=4.5$ and 5. The leading mode is oscillatory. (b) The spectra of the eigenvalue problem of type-II TW at $E_{b}=3$ , $L=4.5$ ; two unstable eigenvalues were detected.

At $L=4$ the type-I TW becomes unstable when $E_{b}>0.66$ , but the unstable mode is not oscillatory until $E_{b}\gtrsim 1.94$ . Numerical simulation shows that type-I TWs lose stability and evolve into type-II TWs for $0.66<E_{b}<1.94$ , and the type-II TWs are stable over this range. For $E_{b}>1.94$ the dynamics develop into an RPO.

For $L=4.5$ , figure 6(c,d) shows that the oscillatory eigenmode is always stable, and the unstable mode is always real. In this case, RPOs do not originate from the instability of type-I TWs, but instead come from the instability of type-II TWs when $E_{b}>2.6$ , as demonstrated by figure 7(a). Figure 4(c) shows that the RPO branch connects with the type-II TW branch.

Figure 8.

The periods of RPOs versus the electric Weber number at different domain sizes $L=4$ , 4.5 and 5.

For $L=5$ , type-I TWs are unstable to a real mode when $E_{b}<0.92$ or $E_{b}>3.92$ , while the mode is unstable and oscillatory for $2.23<E_{b}<3.92$ . It should be stressed that the type-II TW is a stable attractor when $2.23<E_{b}<2.52$ at $L=5$ and numerical simulation confirms that type-II TWs rather than RPOs are observed in this range. For type-II TWs, figure 7(a,b) shows that the type-II TWs become unstable to an oscillatory mode for $E_{b}>2.52$ . Hence, for $L=5$ , RPOs may be born from the instability of either type-I or type-II TWs for $E_{b}>2.52$ .

For the $E_{b}$ values where there are oscillations, time series of the energy norm $\mathscr{E}$ are found to be periodic, implying the existence of RPOs. Using a snapshot of $S(z,t_{0})$ at some $t_{0}$ from the numerical simulations as the initial guess for our Newton–Krylov code, an exact RPO solution can be obtained. Typical profiles of the RPOs at $L=4$ are shown in figure 9(a,b). A spatio-temporal diagram of the film radius $h$ in figure 9(c) illustrates the oscillating behaviour for two different sized droplets. A more clear picture of the oscillation can be seen in our supplementary movie (available at https://doi.org/10.1017/jfm.2019.450). The dynamics of the RPOs is projected into the ${\mathcal{I}}-{\mathcal{D}}$ versus ${\mathcal{D}}$ plane in figure 10(a). TWs are fixed points on the horizon axis ${\mathcal{I}}={\mathcal{D}}$ , while RPOs are closed loops that cross this axis. For $L=4$ , using a parametric continuation method, RPOs are found for $E_{b}$ up to $E_{b}\approx 4$ . Further increase in $E_{b}$ leads to the formation of the spike-like structure. When the length of the computational domain is increased for fixed $E_{b}=3$ , RPOs are found up to around $L=4.75$ , shown in figure 10(b). Further increasing $L$ also leads to formation of the spike-like structure.

Figure 9.

(a,b) The interfacial profiles at times $t$ , $t+T/2$ and $t+T$ . The electric number and time period in (a) and (b) are respectively $(E_{b}=2,T=26.6521,l=0.4520)$ and $(E_{b}=3,T=10.2997,l=-0.9633)$ . (c) The space–temporal diagram of $h(x,t)$ in a period at $E_{b}=4$ . Here, the coordinate $\unicode[STIX]{x1D709}=(2\unicode[STIX]{x03C0}/L)z$ and $L=4$ .

Figure 10.

RPOs in the ${\mathcal{I}}-{\mathcal{D}}$ versus ${\mathcal{D}}$ plane at (a)  $L=4$ and (b)  $E_{b}=3$ . Diamonds are steady TW solutions.

The period of the RPOs, $T$ , is plotted against the electric Weber number $E_{b}$ in figure 8(a). The electric field has dual effects on the oscillation. For $L=4,5$ , the oscillation becomes faster ( $T$ decreases) as $E_{b}$ increases initially. However, beyond some critical point, the oscillation becomes slower ( $T$ increases). For $L=4.5$ , the dynamics is a bit more complicated. The oscillation becomes slower ( $T$ increases) as $E_{b}$ increases from the bifurcation point $E_{b}\approx 2.6$ . Then, the oscillation becomes faster ( $T$ decreases) as $E_{b}>2.78$ . However, as $E_{b}$ increases further, the oscillation becomes slower ( $T$ increases). Referring to figure 4(b–d), the RPO runs faster/slower when the oscillation is slower/faster. The RPO is slower than the type-I TW, which might be faster than type-II TW (see figures 4 c and 15 b).

For $L=4$ , as $E_{b}$ is reduced towards $1.94$ , we find that $T\rightarrow \infty$ and our Newton–GMRES code stops converging, due to the breakdown of the linearization (3.7) for large $T$ . We are able to get to sufficiently large $T$ , however, to identify that the RPO connects to the type-I TW. The divergence of $T$ is in agreement with the linear stability analysis where $\unicode[STIX]{x1D706}_{i}\rightarrow 0$ as $E_{b}\rightarrow 1.94^{+}$ (figure 6 b). For $L=4.5$ , a Hopf bifurcation occurs at $E_{b}\approx 2.6$ for the type-II TW (figure 7 a). The eigenvalue at the bifurcation point $E_{b}=2.6$ predicts a period $T=2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D706}_{i}=2\unicode[STIX]{x03C0}/0.34\approx 18.48$ . The RPO branch is quite contorted as it approaches the type-II branch, but at the last points computed on the branch the period of the RPO is ${\approx}18.70$ and decreasing, in good agreement with the linear stability calculation. For $L=5$ , we have been unable to identify a connection between the RPO and the TWs (see figure 4 d), as far as we are able to reliably perform the calculations.

To summarize the results in this section, the solution regimes are constructed in figure 3, and, in particular, the oscillatory bound states are represented by RPOs. We have demonstrated that the two large droplets cannot survive in a steady bound state but can do so in an oscillatory bound state when $E_{b}$ exceeds a critical value (e.g.  $E_{b}=2$ and see figure 5 a). This indicates that we may expect to see steady travelling waves, relative periodic solutions or spike-like structures in experimental study, depending on the strength of the electric field and droplet spacing. When the electric field is weak, we should observe the formation of a steady liquid viscous bead train due to the Rayleigh–Plateau mechanism. The steady wave train will lose stability and an oscillating state will emerge as the strength of the electric field increases. The spike-like structure should appear when the electric field is strong.

4.2 Invariant solutions in a non-electrified film

In the absence of an electric field, it is observed that droplets can move in either a stable or an unstable bound state (Duprat et al. Reference Duprat, Ruyer-Quil, Kalliadasis and Giorgiutti-Dauphiné2007), but this oscillation is different from that discussed in § 4.1 – here, experimental and numerical studies have shown that the oscillation is caused by the interaction between a large droplet and a much smaller droplet (Craster & Matar Reference Craster and Matar2006; Duprat et al. Reference Duprat, Ruyer-Quil, Kalliadasis and Giorgiutti-Dauphiné2007). The large droplet eats small droplets, and can then break up into small droplets. Such a process is cyclical, and we will demonstrate that it can be represented by RPOs.

Figure 11.

The space–time diagram for type-II TW (a), and type-I TW (b). The bright (dark) shading indicates elevated (depressed) regions. The dependent parameters are $Bo=\unicode[STIX]{x1D716}=0.2$ , $\unicode[STIX]{x1D6FC}=0.25$ and $L=12$ .

Space–time diagrams at $L=12$ are shown in figure 11. In fact, for $L=12$ , other types of solutions do exist, e.g. a one-droplet solution, a three-droplet solution, and so on. The multi-droplet solutions can be stable or unstable bound states. Here we only consider two big droplets in the computational domain so as to isolate the oscillation dynamics. The leading eigenvalues of type-I and type-II TWs are tracked against the domain size $L$ , shown in figure 12. Results in figure 12(a) show that type-II TWs are stable when $L\lesssim 6.7$ . Type-I TWs are unstable when $4.4\lesssim L\lesssim 5.2$ or $L\gtrsim 6.3$ . However, type-I TWs do not lose stability to RPOs when $4.4\lesssim L\lesssim 5.2$ , as $\unicode[STIX]{x1D706}_{i}=0$ . Figure 12(a,b) demonstrates that the most unstable modes are oscillatory when $L\gtrsim 6.7$ for both type-I and type-II TWs, suggesting the existence of RPOs. Interestingly, the sawtooth $\unicode[STIX]{x1D706}_{i}{-}L$ curve originates from the mode switching in the system: the growth rate $\unicode[STIX]{x1D706}_{r}$ of higher-frequency modes increases faster than that of the low-frequency modes, leading to the switching. Therefore, the frequency of the most unstable mode jumps as $L$ increases, leading to the sawtooth structure of the $\unicode[STIX]{x1D706}_{i}{-}L$ curve.

Figure 12.

(a,b) The leading eigenvalue versus the domain size for type-I and type-II TWs. The dependent parameters are $Bo=\unicode[STIX]{x1D716}=0.2$ and $\unicode[STIX]{x1D6FC}=0.25$ .

Figure 13.

The space–time diagrams for the RPOs plotted in three periods at $L=12$ : (a) RPO1 with period $T=1.9122$ ; (b) RPO2 with period $T=1.3501$ ; and (c) RPO3 with period $T=1.0344$ . The dependent parameters are $Bo=\unicode[STIX]{x1D716}=0.2$ and $\unicode[STIX]{x1D6FC}=0.25$ .

Numerical simulations were performed in order to obtain initial guesses for the RPOs. We fixed $L=12$ and obtained three converged exact RPOs, namely RPO1, RPO2 and RPO3, whose space–time diagrams are shown in figure 13. We expect that it is possible to find more RPOs of longer period, but we have not pursued this in the present study. The profiles of the three RPOs at some instants in time are shown in figure 14(b). For RPO1, between the two main big droplets, the film is unstable and develops into a ripple-like structure. For RPO2, the unstable film between the two main droplets evolves into two small droplets. For RPO3, there are three small droplets between the two main droplets. As there are more small droplets between the two main droplets, the mean distance between the big droplet and its nearby small droplet decreases. Therefore, it takes less time for the big droplet to consume the nearby small droplet and the oscillation becomes faster, i.e. the periods decrease from RPO1 to RPO3 (see figure 15 a). We observe that the big droplet eats the front small droplet or ripple and itself is unstable at the rear region, which breaks up into new ripples or small droplets. The breakup of the big droplet at its rear region is due to the Plateau–Rayleigh mechanism. Such a cyclic process (also confirmed in figure 14 a) is self-sustaining, which gives the birth of exact RPOs in the non-electrified film.

Figure 14.

(a) Three exact RPOs in the ${\mathcal{I}}$ – ${\mathcal{D}}$ plane at $L=12$ . (b) The instant profiles of the liquid films for the RPOs. The dependent parameters are $Bo=\unicode[STIX]{x1D716}=0.2$ and $\unicode[STIX]{x1D6FC}=0.25$ . Three movies for these RPOs are shown in the supplementary material. Each movie shows three periods of the RPO.

Figure 15.

(a) Period $T$ versus the domain size $L$ for RPOs. Solutions of different regimes are illustrated. When $L<6.85$ , only TWs are detected. RPO1 is detected when $L>6.85$ ; RPO2 is detected when $L>8.7$ ; and RPO3 is detected when $L>10.7$ . (b) The wave speed versus the domain size $L$ . The dependent parameters are $Bo=\unicode[STIX]{x1D716}=0.2$ and $\unicode[STIX]{x1D6FC}=0.25$ . PD stands for period-doubling bifurcation.

It is interesting to investigate how the period and wave speed of the RPO depend on the domain size. Figure 15(a) shows that the period $T$ becomes longer as the domain size $L$ increases. We speculate that this is because the distance between the big droplet and its front small droplet (the closest one to the big droplet) becomes larger, and therefore longer time is required for the coalescence process. The domain sizes of the three observed different flow regimes in Kliakhandler et al. (Reference Kliakhandler, Davis and Bankhoff2001) are within the current computation, and we do not investigate longer domains in the present study.

Different solution regimes are shown in figure 15(a). RPO1 is not detected when $L<6.85$ . The curve for RPO1 in figure 15(b) connects with the type-I TW solution at $L\approx 6.85$ . We have verified that the spatial form for RPO1 is very similar to the type-I TW at $L=6.85$ , where the ripple between the two big droplets is very small. This indicates that RPO1 is born from the instability of the type-I TW, and explains why we do not detect RPO1 below $L=6.8$ . RPO2 and RPO3 were not detected for $L\lesssim 8.7$ and $10.7$ , respectively. A possible reason why they do not exist in smaller domains may be due to insufficient space for varying numbers of small beads or ripples. The wave speed $\bar{c}$ reduces for a stronger oscillation (RPO1 is faster than RPO2, and RPO2 is faster than RPO3). As the domain size $L$ increases, the wave speed of RPOs increases (see figures 15 b and 16 c). This is because the average height of the large droplet increases and therefore a faster wave (Craster & Matar Reference Craster and Matar2006). When the domain size is large, although these RPOs can be tracked using numerical continuation, they lose stability and the system is typically chaotic; hence we expect a host of RPOs to describe the droplet coalescence and breakup processes discussed by Craster & Matar (Reference Craster and Matar2006). A selection of the shortest orbits, however, is expected to be sufficient to capture most of the dynamics (Budanur et al. Reference Budanur, Short, Farazmand and Willis2017).

Figure 16.

(a) A snapshot of exact RPO. (b) Regime ‘c’ adapted from Kliakhandler et al. (Reference Kliakhandler, Davis and Bankhoff2001). The axis in (a) is rescaled to $\unicode[STIX]{x1D701}=z/\unicode[STIX]{x1D716}$ so as to compare with the experimental observation. (c) Wave speed of the RPO versus domain size, in which there is only one big droplet. A movie for the RPO in the supplementary material demonstrates the unsteady flow. The dependent parameters $\unicode[STIX]{x1D716}=0.178$ , $\unicode[STIX]{x1D6FC}=0.3262$ and $L=8.185$ for regime ‘c’ were taken from Craster & Matar (Reference Craster and Matar2006).

Finally, we revisit the experimental study of the flow regime ‘c’ in Kliakhandler et al. (Reference Kliakhandler, Davis and Bankhoff2001), where there is a large droplet with some small ripples at its front. In contrast to the ‘steady’ flow regime ‘b’, the flow regime ‘c’ is an essentially unsteady state due to the droplet coalescence and breakup process (Craster & Matar Reference Craster and Matar2006). Although the TW solution excellently captures the dynamics of flow regime ‘b’, it is questionable whether or not we can compare a ‘steady’ TW with flow regime ‘c’. Craster & Matar (Reference Craster and Matar2006) observed that the large droplet oscillates in the experiment and proposed that flow regime ‘c’ is a superposition of a steady TW and small perturbations. However, their numerical experiment does not compare well with the experimental study, not only in the wave speed but also in the wave profile. In Ruyer-Quil et al. (Reference Ruyer-Quil, Treveleyan, Giorgiutti-Dauphiné and Kalliadasis2008), the streamwise dissipation effect was considered, which leads to a smaller wave speed ( $c=0.94$ ) than in Craster & Matar (Reference Craster and Matar2006). Novbari & Oron (Reference Novbari and Oron2009) used an energy integral model and predicted a bit higher phase speed than Craster & Matar (Reference Craster and Matar2006), $c=1.48$ , and the streamwise dissipation effect was neglected by them. (The wave speeds by Ruyer-Quil et al. (Reference Ruyer-Quil, Treveleyan, Giorgiutti-Dauphiné and Kalliadasis2008) and Novbari & Oron (Reference Novbari and Oron2009) are rescaled by $W=\unicode[STIX]{x1D70C}h_{0}^{2}g/\unicode[STIX]{x1D707}$ here.) Indeed, Novbari & Oron (Reference Novbari and Oron2009) have realized that flow regime ‘c’ is not a TW because it was unsteady, although they compared it with a TW. All these TWs, however, do not capture a correct film profile: there are no small beads coexisting with the large droplet. Within our set of periodic orbits, we have found profiles of the film very much like experimental observation (figure 16). It should be stressed here that the small droplets/ripples in the front of the big droplet in figure 16 are captured by the RPO, which are not seen in the TW. Kliakhandler et al. (Reference Kliakhandler, Davis and Bankhoff2001) simulated flow regime ‘c’ using their thick film model and showed that these small droplets can survive with big droplets in an unsteady state. The large droplet catches up with and consumes the small droplets in the front, and the rear region of the large droplet breaks up into small droplets. Such a process recurs in the system, i.e. the profile in figure 16(a) appears periodically, which is strikingly similar to that in flow regime ‘c’. Furthermore, the period of the obtained RPO in figure 16(a) is $T=1.7156$ . However, no experimental study on RPOs in this film flow was reported, and the period of flow regime ‘c’ in experiments remains unknown. The wave speed of this RPO $\bar{c}=1.34$ is a little smaller than the wave speed of the TW, $c=1.36$ (Craster & Matar Reference Craster and Matar2006), and is closer to the experimental observation $c=0.93$ . Hence, the RPO captures the dynamical scenario, i.e. oscillation, the wave speed of flow regime ‘c’, and the film profile better than the steady TWs, which provides a better understanding of this unsteady state than the TW solution. The mean wave speed predicted by the RPO, however, is still approximately 30 % higher than the measured one. Note that the wave speed of flow regime ‘a’ is correctly predicted by Craster & Matar (Reference Craster and Matar2006). It can be inferred that the 30 % difference is not caused by dropping the inertia term because the Reynolds number in flow regime ‘c’ is much smaller than in flow regime ‘a’. As is shown by Ruyer-Quil et al. (Reference Ruyer-Quil, Treveleyan, Giorgiutti-Dauphiné and Kalliadasis2008), the streamwise viscous dissipation effect can significantly improve the prediction of the wave speed of steady TWs. This suggests that the streamwise diffusion is important when the Reynolds number is small, which is the reason why the predicted wave speed by the asymptotic model is higher than the experimental value.