Air/water interfacial waves with a droplet at the tip of their crest

In nature, it is common to observe water wave crests with a droplet at their tip. This fascinating configuration remains unexplained from the physical point of view. The present study explores such a unique local configuration numerically. Solitary waves that propagate at the interface between two layers of irrotational fluids are considered. Extending the work of Guan et al. [ “ A local model for the limiting configuration of interfacial solitary waves, ” J. Fluid Mech. 921 , A9 (2021)], the density ratio has been decreased to a very small value equal to 0.001, which is close to the air/water density ratio at sea level (0.0013). A highly accurate solution for the limiting configuration of solitary waves is obtained by solving the irrotational Euler equations using the boundary integral method and Newton ’ s iterations. It is confirmed that the limiting configuration consisting of a droplet sitting on a crest with a 120 (cid:1) angle exists for very small density ratios. This limiting configuration obviously does not exist for surface waves with a void on the top, thus stressing the crucial role played by the air. The droplet is stationary in a frame of reference moving with the wave and experiences intense shear at its tip. From the energy point of view, the formation of a crest with a droplet is accompanied by a remarkable drop of kinetic and potential energies of water in the vicinity of the crest. Furthermore, we present a simple set of scaling relations for the fall of the droplet.


I. INTRODUCTION
Water waves with a droplet on their crest can be observed on the open ocean, on the seashore, and on lakes. Figure 1 shows some examples observed on the seashore in south Dublin, Ireland, and one may also refer to the field observation in Ref. 1 (Fig. 5). This common, yet unexplained, phenomenon may be related to wave breaking as these droplets quickly bend and splash on the sea surface, thus generating foam and spume. At the end of their experimental study, Villermaux and Pomeau 1 commented that the formation of droplets on the crest may be relevant to the superfall of wave crests. Superfall means here that the fluid can experience accelerations that are larger than the acceleration due to gravity only. However, their experiments were carried out with a different model than that of water waves, i.e., free fall of water filling a vertical tube with expanding cross-section and open ends. They did observe a concentrated "nipple" on the top of the essentially flat base solution and attributed it to the instability of the interface due to a pressure gradient during the superfall.
In fact, the droplet-type crest could be one of the limiting configurations of water waves. The investigation of the limiting configuration of wave crests goes back to the conjecture made by Stokes one and a half centuries ago about the crest profile of periodic traveling waves when they reach their maximal amplitude. They are known as the Stokes highest waves. After many years of discussion from asymptotic, numerical, and analytical perspectives, [2][3][4] it is now widely accepted that the limiting configuration of the Stokes highest wave is a stagnation point at the crest with an enclosed angle of 120 . At the same time, the limiting configuration of solitary waves has also been investigated. It has been shown that extreme solitary waves have the same limiting configuration as the periodic waves. [5][6][7][8] Indeed, the limiting configuration at the very tip is localized and independent of wavelength and water depth. 2 Comparing this limiting configuration obtained with a surface wave model to the droplet-type crest, it looks like the surface wave model can hardly reproduce and explain the formation of a droplet on the crest. One possible reason is that the Physics of Fluids ARTICLE scitation.org/journal/phf influence of air which has been neglected so far turns out to be more important than originally thought. In this paper, we take air into account and investigate the limiting configuration of interfacial waves between two layers of homogeneous fluids. The aforementioned sharp crest of 120 is no longer possible. Otherwise, the fluid velocity in the upper layer near the crest would become infinite. 9 A wave profile with a vertical tangent was initially regarded as the configuration of the highest-possible interfacial wave. 10 Later, the existence of overhanging configurations was found numerically. 11,12 Meiron and Saffman 9 asserted that these overhanging solutions would continue to exist until the interface becomes selfintersecting. Finally, a stagnant fluid "bubble" on the top of a 120 angle was proposed by Grimshaw and Pullin 13 as the limiting configuration of interfacial waves. Recently, a highly accurate numerical investigation of the limiting configuration of interfacial waves was conducted (see Ref. 14, Fig. 7), in which the smallest density ratio used was 0.01. This configuration being localized, it is also valid for solitary interfacial waves. 15,16 In fact, the limiting configuration of a stagnant fluid bubble on a wave crest is highly similar to the droplet-type crest. This similarity led us to wonder about their connections.
The density ratio q could be the most crucial number for the localized configuration of the crest under the precondition that the depth allows the wave to develop freely. However, the aforementioned investigations are restricted to moderate values of q. Although interfacial waves with q ¼ 10 À4 were examined in Ref. 17, the overhanging profile with the smallest value of q is only given for q ¼ 0:15. Recently, Guan et al. 18 revisited the limiting configuration of solitary waves and decreased the value of q to 0.01. To overcome the rapid shrinking of the local limiting configuration when decreasing q even further, a simplified model was proposed. In that study, it was found that the formation of the limiting configuration is always accompanied by a decrease in the amplitude. It seems to be consistent with the concept that nipples are caused by the superfall of the crest, see Ref. 1. Considering that the actual density ratio of air to water is of the order of Oð10 À3 Þ, that is at least one order of magnitude smaller than the values used in the existing studies. For such a low density ratio, the existence of a limiting configuration with a stagnant fluid bubble on the crest is still an open question. If this limiting configuration does exist, the actual solution has not been presented yet, and the possible relationship with the droplet-type crests has not been discussed yet.
In the present study, solitary interfacial waves are investigated. We follow Ref. 18 and decrease the value of q to the air-water density ratio ($ 10 À3 ). It turns out that the limiting configuration with a stagnant droplet on the wave crest does exist for such a low density ratio. The actual configuration, the velocity field and the shear along the interface are presented. The variation in the kinetic and potential energies of the fluid in the crest region along the solution branch is discussed. It is shown that the decrease in the wave amplitude is not a

ARTICLE
scitation.org/journal/phf necessary condition for reaching the limiting configuration. Instead, the remarkable drop of mechanical energy inside the crest is. Furthermore, a simple set of scaling relations is introduced for estimating the splash of the main droplet. Our results may be helpful for understanding the common, yet unexplained, droplet-type crest phenomenon.

II. FORMULATION OF THE PROBLEM
We consider a two-dimensional air layer of density q 2 and thickness h 2 over a water layer of density q 1 and thickness h 1 , as illustrated in Fig. 2. The top and bottom boundary conditions are solid horizontal walls. The stratification is overall stable as q 1 > q 2 . An interfacial solitary wave travels at speed c along the x-axis from the left (x ! À1) to the right (x ! 1), and it is symmetric with respect to the vertical y-axis. The system of coordinates moves with the wave, and the level y ¼ 0 is the undisturbed level of the interface. Gravity acts in the negative y direction. The quantities h 1 , c, and q 1 are used as length, velocity, and density units, respectively, to nondimensionalize the system.
Under the assumption of irrotational flow, the motion of water and air can be described by the Laplace equation (1) / j is the nondimensional velocity potential, and X j is the fluid domain shown in Fig. 2. Here and hereafter, the subscripts j ¼ 1 and 2 correspond to the water layer and the air layer, respectively. The two kinematic and the dynamic boundary conditions on the interface read @/ j @y À @/ j @x @g @x and qjr/ 2 j 2 À jr/ 1 j 2 þ 2ðq À 1Þ g Here, q ¼ q 2 =q 1 represents the air/water density ratio, g is the nondimensional elevation of the interface, F ¼ c= ffiffiffiffiffiffi ffi gh 1 p is the Froude number, and g is the acceleration due to gravity.
The two fluid layers are bounded by the horizontal upper and bottom walls. The corresponding boundary conditions can be expressed as The depth ratio h is defined as h 2 =h 1 .
In the frame of reference moving with the wave speed, the boundary conditions at infinity for the solitary wave are The solitary waves are characterized by three dimensionless parameters: h, q, and F. The impact of h is insignificant when it is large so that the wave can develop freely. According to Kataoka (Ref. 17 ,  Table I), the threshold value is of the order of h > 5 for small q. In this study, h is fixed at 80. This implies that when q will be further fixed at 10 À3 , a value close to the air/water density ratio, there will be only one governing parameter: the Froude number F.

III. NUMERICAL METHOD
The numerical method closely follows the numerical method used in Ref. 18. Some details are repeated for the sake of completeness, and the differences are emphasized. The interface is parameterized by the arc length s as x ¼ XðsÞ and y ¼ gðsÞ. It gives the arc length equation X 0 ðsÞ 2 þ g 0 ðsÞ 2 ¼ 1: Let s ¼ 0 at x ¼ 0. There are two conditions for X and g Here, A is the nondimensional wave amplitude. Following Refs. 18 and 19, the Laplace equation (1) can be reformulated using the Cauchy integral formula. Let u and v denote the horizontal and vertical velocities, respectively. Introducing W j ðzÞ ¼ u j ðzÞ À iv j ðzÞ þ 1, which is an analytic function of z ¼ x þ iy, the Cauchy integral formula gives for z 0 on the boundary of the integration domain X j . We use the Schwarz reflection principle to satisfy the boundary conditions on the upper and bottom walls (4). The integration domain in (8) for water can be X 1 plus its reflection about the bottom wall, which is bounded by g and g R;1 . Similarly, for the air, its integration domain X 2 is bounded by g and g R;2 . g R;j are expressed as Based on the arc length parameterization, the integration in (8) carried out in the original plus reflected domains can be rewritten as

Physics of Fluids
According to (5), W j ! 0 when x ! 1. Thus, the integral above can be calculated in s 2 ðÀL; LÞ with large L. Using the symmetry of the interface with respect to the y-axis and taking the real part, (10) can be finally simplified as Here, X 6 ¼ XðsÞ6Xðs 0 Þ and j ¼ 1, 2.
Meanwhile, the kinematic boundary condition (2) can also be written in the parameterized form v j ðsÞX 0 ðsÞ À u j ðsÞg 0 ðsÞ ¼ 0; j ¼ 1; 2: The integration domain s 2 ð0; LÞ is divided into N À 1 elements by seeds s k , where k ¼ 1; 2; 3; …; N. The length of the kth element ds k is equal to s kþ1 À s k , and the middle point of each element r k is equal to ðs kþ1 þ s k Þ=2, where k ¼ 1; 2; 3; …; N À 1. In Ref. 18, uniform meshes were used. However, as the size of the droplet-type crest decreases significantly with the decrease in q, nonuniform meshes are adopted in the integral domain in the present paper. More specifically, uniform finer meshes for the crest region, uniform larger meshes in the tail, and meshes with fixed stretching ratio ds k =ds kÀ1 ¼ 1:004 in between the crest and tail are used.
The midpoint rule is applied for the integral equations (11) to avoid the singularities. It gives 2N À 2 algebraic equations. The boundary conditions of the interfaces (3), (6), and (12) are discretized on the seeds, which gives 2N, N, and N algebraic equations, respectively. Using the boundary condition (7), the number of unknown variables can be reduced to 6N, which are x 0 ðs k Þ; g 0 ðs k Þ; u j ðs k Þ; v j ðs k Þ; k ¼ 1; 2; 3; …; N. Together with the conditions g 0 ð0Þ ¼ 0 and u 1 ðLÞ ¼ À1, we are ready to search the branch of steady solutions in the F À A plane.
Lagrange second-order interpolation is adopted for the derivatives. Iterations are carried out by Newton's method with a maximum residual error less than 10 À9 . When the first turning point on the branches appears, Lagrange third-order interpolation and a maximum residual error threshold of 10 À11 are adopted.
The solution branches are explored by the following procedures: (i) a Gaussian profile with infinitesimal amplitude is produced as an initial guess of the solitary wave, and then, the first solution, a point on the branch, can be obtained by iteration; (ii) the former solution is used as the initial guess for the next step with a slight change of A (or F); thus, the next new solution can be obtained by iteration; (iii) repeat step (ii) until the droplet-type configuration develops; and (iv) the solution branch in the F À A plane is then built by collecting all the solutions. Figure 3 shows that the mesh size at the crest ds 1 has a significant impact on the crest profile, whereas the domain size L mainly affects the profile of the tail of the wave and the wave amplitude. Using 1900 elements with ds 1 ¼ 1:0 Â 10 À3 and L ¼ 9.2 was found to be a good balance between computational cost and precision.  Fig. 5. It is shown that the angle of the crest decreases and converges to 120 along the branch. Figure 5

ARTICLE
scitation.org/journal/phf confirms that the droplet-type crest does exist for the small density ratio q ¼ 0:001 and that it develops after the second turning point S6. Figure 4 also shows that the wave amplitude increases monotonically along the solution branch. However, for the solution branches with larger values of the density ratio q (see Fig. 14 for q ¼ 0.1 and illustrations for other density ratios in Ref. 18), the amplitudes always steepen along the solution branch to reach the highest wave and then decrease to form the droplet-type crest. Moreover, Villermaux and Pomeau 1 conjectured that the nipples on the wave crest may be caused by its superfall. It may give an illusion that the formation of the droplet-type crest is always accompanied by the decrease in A. However, the case q ¼ 0:001 considered in the present study indicates that a decrease in A is not a necessary condition for the formation of the limiting configuration.
To extend the solution branch past S8, much finer meshes are needed as the length of the neck region of the droplet crest will decrease further. In Ref. 18, a simplified model was introduced: the crest was replaced by two straight solid walls intersecting at the origin with an angle of 120 and a closed fluid bubble with a flat bottom on the top of the angle. It was found that the profiles for different values of the density ratio q are geometrically similar. The limiting configuration obtained with the local model is also presented for q ¼ 0:001 in Fig. 5(b). Comparing the S8 profile with the local model, it is seen that the neck region may shrink further, leading to selfintersection after S8 and the local model may still be a good representation of the true solution.
Next, we explore the velocity fields. The velocity components in the moving frame at any location z i in X j , j ¼ 1, 2, are calculated again using the Cauchy integral formula, yielding Let us take S4 and S8 as solutions that are representative of the classical and droplet-type crests, respectively. Their velocity fields with the wave and interfacial shear velocity are plotted in Fig. 6. The components of the interfacial shear velocity (Du; Dv) are (u 1 À u 2 ; v 1 À v 2 ). Figures 6(a) and 6(c) show that the most intense shear occurs at the tips because of the significant slowdown of the water (in the moving frame) and speedup of air around the crest, and the shear velocity can be more than twice the wave speed. For the droplet-type situation, Fig. 6(b) indicates that the droplet at the wave tip is almost stationary. The most severe shear still occurs at the crest. However, it has been further strengthened to be more than five times larger than the wave speed, see Fig. 6(d). It is interesting that there exists a segment on the interface without shear, denoted by S b , which corresponds to the segment labeled by the blue line at the bottom of the droplet, see Fig. 6(b). In fact, the velocity components at segment S b vanish as well.
The ratio between the speed of the water particle at the tip and the wave speed, u 1t þ 1, along the solution branch is plotted in Fig. 7(a). It can be seen that the speed of water at the tip approaches the wave speed with an obvious speedup after point S4 along the solution branch. According to the kinematic wave breaking criterion (e.g., Ref. 20), it suggests that breaking won't occur for the droplet-type crest if 1.0 is regarded as the threshold ratio between the particle speed at the tip and the wave speed. Note that the value of 1.0 was revisited by Barthelemy et al. 21 and brought down to 0.85, in which case breaking would occur. Figure 7(b) shows the square of the shear velocity at the tip, along the solution branch. It seems that there is a simple exponential relation FIG. 4. Solution branch in the F À A plane for q ¼ 0:001 and h ¼ 80. The points labeled S1-S8 correspond to typical solutions selected to observe the wave profile and the shear at the interface.

ARTICLE
scitation.org/journal/phf with a ¼ 2:1 for small and moderate amplitude and a ¼ 42:5 for extreme high amplitude.

B. Crest mass and energies
The kinetic and potential energies in a domain X 0 around the crest, called the "crest domain" hereafter, are of interest as the limiting configuration is localized. The kinetic energy E T;j , the potential energy E P;j , and the mass M j of water and air in that domain are calculated as where j ¼ 1, 2. It is noteworthy that E T;j is calculated in the still frame of reference, and all the energies are nondimensionalized by q 1 gh 3 1 . The mass and energy can be influenced by the shape and size of the crest domain. In practice, the crest domain is a square domain of length l, and the wave tip is located at the height 5=6l and the horizontal center as marked by the domain X 0 in Fig. 2 7. Ratio between the speed of the water particle at the tip of the crest and wave speed (a) and shear velocity at the very tip of the wave (b). u jt ¼ u j (s ¼ 0), j ¼ 1, 2, q ¼ 0.001 and h ¼ 80. With increasing wave amplitude along the solution branch shown in Fig. 4, the wave crest steepens from an initially flat interface and finally develops into a droplet-type crest. This means that part of the water in the crest domain is replaced by air, resulting in the decrease in M 1 . Figure 8 shows that up to 38% of water has been lost along the solution branch. The sharp decrease starts at S4, the first turning point in the solution branch.

Physics of Fluids
In general, the kinetic and potential energies of water and air are divided into two stages along the solution branch, i.e., a gentle variation stage followed by a sharp rise/drop stage, see Figs. 9 and 10. However, the turning points are not synchronized exactly. In fact, the size of the crest domain l can also influence the location of the turning point slightly. Figure 9 indicates that E T;1 and E P;1 drop significantly during the formation of the droplet-type crest, although the wave amplitude does not decrease.
With the increase in the velocity of water at the tip, the kinetic energy of water in the crest domain first increases. The sharp drop of E T;1 in Fig. 9(a) between S6 and S8 is mainly due to the remarkable decrease in the water mass in the crest domain, see Fig. 8. This implies that the variation in the kinetic energy of water is a result of the competition between an increase in the velocity and a decrease in the mass in the crest domain. The sharp drop of E P;1 in Fig. 9(b) is also due to the decrease in M 1 .
The monotonic increase in E T;2 and E P;2 is straightforward, because the velocity of air at the tip keeps increasing up to more than five times larger than the wave speed [see Fig. 6(d)], and up to 38% of water is monotonically replaced by air (see Fig. 8) along the solution branch. The insignificant difference of E T;2 between S7 and S8 is due to the flow of air by the side of the neck, i.e., the region between the profiles of S7 and S8 [see Fig. 5(b)] is characterized by small velocities.
The total energy of water and air, E ¼ RE T;j þ RE P;j , in the crest domain ("crest energy" hereafter) is plotted vs F in Fig. 11. The maximum value of the crest energy as a function of F is marked with a red cross. Figures 11(a)-11(d) show that the location of the maximum of the crest energy is sensitive to the size of the crest domain and shifts from S4 to S6 when l decreases. This is because the energy of the smaller crest domain tends to reflect more the energy variation of the very tip. The droplet-typed crest formed at the very tip after S6 is accompanied by a significant energy drop. Overall, the formation of a droplet-type crest won't occur before that maximum.
Solitary waves can be regarded as waves with infinite wavelength. Therefore, the wavelength of any disturbance will be shorter than the solitary wave. The stability of the basic wave to disturbances with a smaller wavelength is called superharmonic instability. 22 The exchange of superharmonic instability occurs at stationary values of the total wave energy relative to the wave speed, as first found numerically by Tanaka 23 for periodic gravity waves in deep water. It was then analytically proved 24 and extended to surface solitary waves. 25 However, the situation is subtle for waves in water of finite depth. Indeed, Kataoka 26 showed that the exchange of stability for surface gravity waves occurs when the total wave energy becomes stationary as a function of the wave speed for fixed mean surface height and not for fixed Bernoulli's constant. Kataoka 17 showed that the energy criterion is also valid for interfacial waves. The application of the criterion for exchange of superharmonic instability to the branch of solutions computed in the present work is left for future work that will be fully dedicated to the stability of the present solitary waves.
We now present the results on the total energy in the whole computation domain as they are interesting in themselves, even without mentioning stability. This energy, which we call the "global energy" hereafter, is calculated as  Physics of Fluids ARTICLE scitation.org/journal/phf boundary condition and (6), and nondimensionalizing by q 1 gh 3 1 , the energies are expressed as and Here, n is the normal to the interface. The global energy along the branch with fixed Bernoulli's constant is presented in Fig. 12. The global energy is insensitive to L as long as L is large enough. Indeed, uj x¼L ¼ vj x¼L $ 0, so there is no contribution to the kinetic energy. As for the potential energy, it is on the undisturbed free surface, and the constant term has been dropped during integration. The first stationary point ðdE G =dF ¼ 0Þ occurs before S4, and the second stationary point occurs around S6.
Comparing Fig. 12 with Fig. 11, we conclude that the energy will first reach its maximum globally and then locally at the tip along the solution branch.

C. Fall of the droplet
So far, we considered only stationary solutions in the frame of reference moving with the wave. In this section, we imagine that we consider time-dependent solutions. As time increases, the droplet will fall on the front side of the wave under gravity, considering the weak support from the neck, and finally, it will splash on the interface forming foam resulting in intense energy transfer and

ARTICLE
scitation.org/journal/phf dissipation. At the same time, the wave will experience stretching as shown by the red profile in the sketch of Fig. 13. The stretched water body may turn into a thin film. This thin film can encounter significant aerodynamic effects and air entrainment and may further break into some micro-droplets released to the air. There exist two different scenarios: one is that the droplet can fall a distance approximately equal to A, designated as uninterrupted fall, and the other one is that the droplet splashes on the slope profile of the wave itself with a fall that is much smaller than A, designated as interrupted fall. During the short period of the fall of the droplet, we introduce the following assumptions: (i) the limiting configuration of the wave consists of a droplet on the top, a main water body of the wave which is confined by two flat interfaces intersecting with an angle of 120 , and a neck region in between connecting them, as illustrated by the black dash line in Fig. 13; (ii) the droplet falls freely under gravity, and the aerodynamic force is much smaller than the gravity force; and (iii) the main water body of the wave propagates at a speed much slower than the wave speed with insignificant deformation. In fact, the last two assumptions imply that intense stretching would occur at the neck region. Based on the assumption of free fall, the timescale for the main droplet to fall a distance s v is ð2s v h 1 =gÞ 1=2 , where s v is normalized by h 1 and 0 < s v < A. If cð2s v h 1 =gÞ 1=2 is always larger than s v h 1 tan ðp=3Þ, it goes to the uninterrupted fall. Otherwise, it goes to the interrupted fall. That gives a criterion identifying these two scenarios, which reads The criterion (21) has been plotted in the F À A plane with the solution branches of q ¼ 0:001 and 0.1, as shown in Fig. 14. It suggests that the main droplet will experience an uninterrupted fall for the case of q ¼ 0:001, whereas an interrupted fall of the main droplet can occur for a larger density ratio as indicated by the solution branch of q ¼ 0:1.
Using the fact that the droplet falls a distance A without interruption, the length scale of the droplet and its impacting velocity and angle can be further estimated. As discussed in Fig. 6 Fig. 13). Accordingly, the dynamic boundary condition (3) can be simplified, which gives the vertical coordinates of S a and S b as 1=2F 2 þ 1=2qu 2 2t F 2 =ð1 À qÞ and 1=2F 2 , respectively. Their difference gives the following dimensional length scale l s for the droplet: This relationship for the length scale is valid as long as the droplettype solution exists, whatever the value of the density ratio is. However, measuring the speed u 2t is usually not easy. We may estimate u 2t using the exponential relation in (14), u 2 2t $ A a , with a ¼ 42:5, according to Fig. 7(b).
Based on the free fall timescale ð2Ah 1 =gÞ 1=2 and on the wave speed, we can also obtain the scales of the dimensional impact velocity V s and of the angle h s . They read

Physics of Fluids
The combination of (22) and (23) may lead to a quantitative analysis of the splashing event using the corresponding splashing models.

V. CONCLUSION
The local limiting crest configuration of progressive solitary interfacial waves confined by two walls parallel to the undisturbed interface is fully determined by three dimensionless numbers: the depth ratio, the density ratio, and the Froude number. In this study, the density ratio is set equal to the small value of 0.001 that corresponds roughly to that of the air-water situation at sea level; the depth ratio is fixed at 80-it is large enough to allow the solitary waves to develop freely. The irrotational flows in both layers are numerically solved using the boundary integral equation method and Newton's iteration method. The solution branch in the amplitude-Froude number plane is obtained, and a highly accurate solution of the limiting configuration with the maximum residual error less than 10 À11 is presented.
Along the solution branch, there exist two different stages. In the first stage, the wave has a non-overhanging crest and the energy increases. Then, a limiting configuration characterized by a droplet sitting on the wave crest with an angle of 120 can be observed. The formation of the droplet-type crest is accompanied by a remarkable drop of both kinetic and potential energies of the water in the crest. This limiting configuration experiences intense shear at the tip, and the shear velocity can be more than five times larger than the wave speed. Furthermore, a simple set of scaling relations is also presented for the fall of the droplet, including the length scale of the droplet and the scales of its impacting velocity and angle. Our results confirm that the droplet-type crest exists even for the low air-water density ratio, without the presence of surface tension. It may be helpful for understanding the common, yet unexplained, phenomena of water waves with a droplet on their crest. The stability analysis is left for future work.