Amelia Jiménez-Sánchez

Shadi Albarqouni

Diana Mateus

Capsule Networks against Medical Imaging Data Challenges

Keywords: capsule networks, small datasets, class imbalance

A key component to the success of deep learning is the availability of massive amounts of training data. Building and annotating large datasets for solving medical image classification problems is today a bottleneck for many applications. Recently, capsule networks were proposed to deal with shortcomings of Convolutional Neural Networks (ConvNets). In this work, we compare the behavior of capsule networks against ConvNets under typical datasets constraints of medical image analysis, namely, small amounts of annotated data and class-imbalance. We evaluate our experiments on MNIST, Fashion-MNIST and medical (histological and retina images) publicly available datasets. Our results suggest that capsule networks can be trained with less amount of data for the same or better performance and are more robust to an imbalanced class distribution, which makes our approach very promising for the medical imaging community.

Introduction

Currently, numerous state of the art solutions for medical image analysis tasks such as computer-aided detection or diagnosis rely on Convolutional Neural Networks (ConvNets) [START_REF] Litjens | A survey on deep learning in medical image analysis[END_REF]. The popularity of ConvNets relies on their capability to learn meaningful and hierarchical image representations directly from examples, resulting in a feature extraction approach that is flexible, general and capable of encoding complex patterns. However, their success depends on the availability of very-large databases representative of the full-variations of the input source. This is a problem when dealing with medical images as their collection and labeling are confronted with both data privacy issues and the need for time-consuming expert annotations. Furthermore, we have poor control of the class distributions in medical databases, i.e. there is often an imbalance problem. Although strategies like transfer learning [START_REF] Zhou | Using convolutional neural networks and transfer learning for bone age classification[END_REF], data augmentation [START_REF] Vasconcelos | Increasing deep learning melanoma classification by classical and expert knowledge based image transforms[END_REF] or crowdsourcing [START_REF] Albarqouni | Aggnet: Deep learning from crowds for mitosis detection in breast cancer histology images[END_REF] have been proposed, data collection and annotations is for many medical applications still a bottleneck [START_REF] Cardoso | Intravascular Imaging and Computer Assisted Stenting, and Large-Scale Annotation of Biomedical Data and Expert Label Synthesis[END_REF].

ConvNets' requirement for big amounts of data is commonly justified by a large number of network parameters to train under a non-convex optimization [START_REF]MICCAI Grand Challenge Tumor Proliferation Assessment Challenge (TUPAC16)[END_REF] shows the difference between the sigmoid and squashing functions. Eq. ( 2) is a weighted sum of the inputs (ConvNets use bias). In CapsNets, c ij are the coupling coefficients. In [START_REF] Cardoso | Intravascular Imaging and Computer Assisted Stenting, and Large-Scale Annotation of Biomedical Data and Expert Label Synthesis[END_REF], ûj|i is the transformed input to the j-th capsule/neuron. In CapsNets, the input from the i-th capsule is transformed with the weights W ij . While in ConvNets, the raw input from the previous neuron is used. scheme. We argue, however, that part of these data requirements is there to cope with their poor modeling of spatial invariance. As it is known, purely convolutional networks are not natively spatially invariant. Instead, they rely on pooling layers to achieve translation invariance, and on data-augmentation to handle rotation invariance. With pooling, the convolution filters learn the distinctive features of the object of interest irrespective of their location. Thereby losing the spatial relationship among features which might be essential to determine their class (e.g. the presence of plane parts in an image does not ensure that it contains a plane).

Recently, capsule networks [START_REF] Sabour | Dynamic routing between capsules[END_REF] were introduced as an alternative deep learning architecture and training approach to model the spatial/viewpoint variability of an object in the image. Inspired by computer graphics, capsule networks not only learn good weights for feature extraction and image classification but also learn how to infer pose parameters from the image. Poses are modeled as multidimensional vectors whose entries parametrize spatial variations such as rotation, thickness, skewness, etc. As an example, a capsule network learns to determine whether a plane is in the image, but also if the plane is located to the left or right or if it is rotated. This is known as equivariance and it is a property of human one-shot learning type of vision.

In this paper, we experimentally demonstrate that the equivariance properties of CapsNets reduce the strong data requirements, and are therefore very promising for medical image analysis. Focusing on computer-aided diagnosis (classification) tasks, we address the problems of the limited amount of annotated data and imbalance of class distributions. To ensure the validity of our claims, we perform a large number of controlled experiments on two vision (MNIST and Fashion-MNIST) and two medical datasets that targets: mitosis detection (TUPAC16) and diabetic retinopathy detection (DIARETDB1). To the best of our knowledge, this is the first study to address data challenges in the medical image analysis community with Capsule Networks.

Methods

In the following, we focus on the image classification problem characteristic of computer-aided diagnosis systems. Our objective is to study the behavior of Capsule Networks (CapsNets) [START_REF] Sabour | Dynamic routing between capsules[END_REF] in comparison to standard Convolutional Networks (ConvNets) under typical constraints of biomedical image databases, such as a limited amount of labeled data and class imbalance. We discuss the technical advantages that make CapsNets better suited to deal with the above-mentioned challenges and experimentally demonstrate their improved performance.

Capsule vs Convolutional Networks

Similar to ConvNet approaches, CapsNets build a hierarchical image representation by passing an image through multiple layers of the network. However, as opposed to the tendency towards deeper models, the original CapsNet is formed with only two layers: a first primary caps layer, capturing low-level cues, followed by a specialized secondary caps, capable of predicting both the presence and pose of an object in the image. The main technical differences of CapsNets w.r.t. ConvNets are: i) Convolutions are only performed as the first operation of the primary caps layer, leading as usual to a series of feature channels.

ii) Instead of applying a non-linearity to the scalar outputs of the convolution filters, CapsNets build tensors by grouping multiple feature channels (see the grid in Fig. 1). The non-linearity, a squashing function, becomes also a multidimensional operation, that takes the j -th vector s j and restricts its range to the [0,1] interval to model probabilities while preserving the vector orientation. The result of the squashing function is a vector v j , whose magnitude can be then interpreted as the probability of the presence of a capsule's entity, while the direction encodes its pose. v j is then the output of the capsule j.

iii) The weights W ij connecting the i primary capsule to the the j -th secondary capsule are an affine transformation. These transformations allow learning part/whole relationships, instead of detecting independent features by filtering at different scales portions of the image.

iv) The transformation weights W ij are not optimized with the regular backpropagation but with a routing-by-agreement algorithm. The principal idea of the algorithm is that a lower level capsule will send its input to the higher level capsule that agrees better with its input, this way is possible to establish the connection between lower-and higher-level information (refer to [START_REF] Sabour | Dynamic routing between capsules[END_REF] for details).

v) Finally, the output of a ConvNet is typically a softmax layer with crossentropy loss:

L ce = -x g l (x) log(p l (x)).
Instead, for every secondary capsule, CapsNet computes the margin loss for class k:

L k = T k max(0, m + -||v k ||) 2 + λ (1 -T k ) max(0, ||v k || -m -) 2 , (1)
where the one-hot encoded labels T k are 1 iff an entity of class k is present and m + = 0.9 and m -= 0.1, i.e. if an entity of class k is present, its probability is expected to be above 0.9 (||v k || > 0.9), and if it is absent ||v k || < 0.1. Since the threshold is not set as 0.5, the marginal loss forces the distances of the positive instances to be close to each other, resulting in a more robust classifier. The weight λ = 0.5. As regularization method, CapsNet uses a decoder branch composed of two fully connected layers of 512 and 1024 filters respectively. The loss of this branch is the mean square error between the input image x and its reconstruction x both of size N × M ,

L MSE = 1 N • M N n=1 M m=1 (x(n, m) -x(n, m)) 2 ) (2)
The final loss, is a weighted average of the margin loss and the reconstruction loss

L total = N k k=1 L k + α L M SE .

Medical Data Challenges

It is frequent for medical image datasets to be small and highly imbalanced. Particularly, for rare disorders or volumetric segmentation, healthy samples are the majority against the abnormal ones. The cost of miss-predictions in the minority class is higher than in the majority one since high-risk patients tend to be in the minority class. There are two common strategies to cope with such scenarios: i) increase the number of data samples and balance the class distribution, and ii) use weights to penalize stronger miss-predictions of the minority class.

We propose here to rely on the equivariance property of CapsNets to exploit the structural redundancy in the images and thereby reduce the number of images needed for training. For example, in Fig. 1, we can see a fundus image in which diabetic retinopathy is present. There are different patterns present in the image that could lead to a positive diagnosis. Particularly, one can find soft and hard exudates or hemorrhages. While a ConvNet would tend to detect the presence of any of these features to make a decision, CapsNet routing algorithm is instead designed to learn to find relations between features. Redundant features are collected by the routing algorithm instead of replicated in several parts of the network to cope with invariance. We claim that the above advantages directly affect the number of data samples needed to train the networks. To demonstrate

N k caps 16D 1 × 1 512 ch 1 × 1 1024 ch 8.2M
Table 1: Details of each of the architectures. For convolution, we specify the size of the kernel and the number of output channels. In the case of pooling, the size of the kernel. And for capsule layers, first, the number of capsules and, in the second row, the number of dimensions of each capsule.

our hypothesis we have carefully designed a systematic and large set of experiments comparing a traditional ConvNet: LeNet [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] and a standard ConvNet:

Baseline from [START_REF] Sabour | Dynamic routing between capsules[END_REF], against a Capsule Network [START_REF] Sabour | Dynamic routing between capsules[END_REF]. We focus on comparing their performance with regard to the medical data challenges to answer the following questions:

• How do networks behave under decreasing amounts of training data?

• Is there a change in their response to class-imbalance?

• Is there any benefit from data augmentation as a complementary strategy?

To study the generalization of our claims, our designed experiments are evaluated on four publicly available datasets for two vision and two medical applications: i) Handwritten Digit Recognition (MNIST), ii) Clothes Classification (FASHION MNIST), iii) Mitosis detection, a sub-task of mitosis counting, which is the standard way of assessing tumor proliferation in breast cancer images (TU-PAC16 challenge [START_REF]MICCAI Grand Challenge Tumor Proliferation Assessment Challenge (TUPAC16)[END_REF]), and iv) Diabetic Retinopathy, an eye disease, that due to diabetes could end up in eye blindness over time. It is detected by a retinal screening test (DIARETDB1 dataset). Next, we provide some implementation details of the compared methods.

Architectures Since research of capsules is still in its infancy, we pick the first ConvNet, LeNet [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF] for a comparison. Though this network has not many parameters (approx. 60K), it is important to notice the presence of pooling layers which reduce the number of parameters and lose the spatial relationship among features. For a fairer comparison, we pick another ConvNet with similar complexity to CapsNet, in terms of training time, that has no pooling layers, which we name hereafter Baseline and was also used for comparison in [START_REF] Sabour | Dynamic routing between capsules[END_REF].

LeNet has two convolutional layers of 6 and 16 filters. Kernels are of size 5x5 and stride 1. Both are followed by a ReLU and pooling of size 2x2. Next, there are two fully connected layers with 120 and 84 filters. Baseline is composed of three convolutional layers of 256, 256, 128 channels, with 5x5 kernel and stride of 1. Followed by two fully connected layers of size 382, 192 and dropout. In both cases, the last layer is connected to a softmax layer with cross-entropy loss. For CapsNet [START_REF] Sabour | Dynamic routing between capsules[END_REF], we consider two convolutional layers of 256 filters with kernel size of 9x9 and stride of 1. Followed by two capsule layers of 8 and 16 dimensions, respectively, as depicted in Fig. 1. For each of the 16-dimensional vectors that we have per class, we compute the margin loss like [START_REF] Sabour | Dynamic routing between capsules[END_REF] and attach a decoder to reconstruct the input image. Details are summarized in Table 1.

Implementation. The networks were trained on a Linux-based system, with 32 GB RAM, Intel(R) Core(TM) CPU @ 3.70 GHz and 32 GB GeForce GTX 1080 graphics card. All models were implemented using Googles Machine Learning library TensorFlow4 . The convolutional layers are initialized with Xavier weights [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF]. All the models were trained in an end to end fashion, with Adam optimization algorithm [6], using grayscale images of size 28×28. The batch size was set to 128. For MNIST and Fashion-MNIST, we use the same learning rate and weight for the reconstruction loss as [START_REF] Sabour | Dynamic routing between capsules[END_REF], while for AMIDA and DIARETDB1 we reduced both by 10. If not otherwise stated, the models were trained for 50 epochs. The reported results were tested at minimum validation loss.

Experimental Validation

Our systematic experimental validation compares the performance of LeNet, a Baseline ConvNet and CapsNet with regard to the three mentioned datachallenges, namely the limited amount of training data, the class-imbalance, and the utility of data-augmentation. We trained in total 432 networks, using 3 different architectures, under 9 different data conditions, for 4 repetitions, and for 4 publicly available datasets. The two first datasets are the well known MNIST [START_REF] Lecun | MNIST handwritten digit database[END_REF] and Fashion-MNIST [START_REF] Xiao | Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms[END_REF], with 10 classes and, 60K and 10K images for training and test respectively.

For mitosis detection, we use the histological images of the first auxiliary dataset from the TUPAC16 challenge [START_REF]MICCAI Grand Challenge Tumor Proliferation Assessment Challenge (TUPAC16)[END_REF]. There are a total of 73 breast cancer images, of 2K × 2K pixels each, and with the annotated location coordinates of the mitotic figures. Images are normalized using color deconvolution [START_REF] Vahadane | Structure-preserving color normalization and sparse stain separation for histological images[END_REF] and only the hematoxylin channel is kept. We extract patches of size 100 × 100 pixels that are downsampled to 28 × 28, leading to about 60K and 8K images for training and test respectively. The two classes are approximately class-wise balanced after sampling.

For the diabetic retinopathy detection, we consider DIARETDB1 dataset [START_REF] Kalesnykiene | Diaretdb1 diabetic retinopathy database and evaluation protocol 6[END_REF]. It consists of 89 color fundus images of size 1.1K × 1.5K pixels, of which 84 contain at least mild signs of the diabetic retinopathy, and 5 are considered as normal. Ground truth is provided as masks. We enhance the contrast of the fundus images by applying contrast limited adaptive histogram equalization (CLAHE) on the lab color space and keep only the green channel. We extract patches of 200 × 200 pixels that are resized to 28 × 28. This results in about 50K and 3K images for training and test respectively. They are approximately class-wise balanced after sampling.

Limited amount of training data

We compare the performance of the two networks for the different classification tasks when the original amount of training data is reduced to 50%, 10%, 5%, and 1% while keeping the original class distribution. We run each of the models for the same number of iterations that are required to train 50 full epochs using all the training data. Early-stop is applied if the validation loss does not improve in the last 20 epochs.

The results are shown in Table 2a. For almost all scenarios CapsNet performs better than LeNet and Baseline. We can observe in Figure 2 how for MNIST the gap is higher for a small amount of data and is reduced when more data is included. LeNet with 5% of the data has a similar performance to CapsNet, and better than Baseline, with 1% of the data for DIARETDB1. We attribute this behavior to the structures that are present in this type of images. All the experiments validated the significance test with a p-value < 0.05, except for those on the TUPAC16 dataset, we presume this is associated to the CapsNet limitations that we present in Section 4.

Class-imbalance

For the medical datasets, we simulate class imbalance by reducing to 20% one of the two classes. Initially, we reduce abnormal class and, afterward, the healthy class. For the other two datasets, we decrease two classes at the same time. For MNIST, we first consider reducing the classes "0" and "1" and secondly, the classes "2" and "8". Similar for Fashion-MNIST, we reduce the classes "Tshirt/top" and "Trouser", and in the second scenario, "Pullover" and "Shirt".

In Table 2b results are reported. Again, CapsNet surpasses the performance of ConvNets for all cases, except for Fashion-MNIST where the f1-scores are similar. At least one of the imbalance cases verified the significance test for all datasets.

Data augmentation

In the last series of experiments, we compare the performance of the three networks using data augmentation, a common technique to increase the amount of training data and balance class distributions. The original dataset is augmented with ±10 degrees rotations, with a translation of ±30 pixels for medical datasets, and with flips (horizontal for Fashion-MNIST and, both horizontal and vertical for TUPAC16 and DIARETDB1). MNIST and Fashion-MNIST are augmented by 5%, for the other two datasets we consider the no augmented version to be 50% (TUPAC16) and 90% (DIARETDB1) smaller.

The performances in Table 2c show that, CapsNet without data augmentation achieves a similar (TUPAC16, MNIST, Fashion-MNIST) or even better (DI-ARETDB1) performance than ConvNets using data augmentation. All results are significant, the only Baseline for MNIST is comparable to the performance of CapsNet. These results confirm the benefits of equivariance over invariance.

Conclusion

In this work, we experimentally demonstrate the effectiveness of using Cap-sNet to improve CADx classification performance under medical data challenges. In particular, we demonstrate the increased generalization ability of CapsNets vs. ConvNets when dealing with the limited amount of data and class-imbalance. The performance improvement is a result of CapsNets equivariance modeling, that is, its ability to learn pose parameters along with filter weights. Together with the routing-by-agreement algorithm, this paradigm change requires to see fewer viewpoints of the object of interest, and therefore fewer images, in order to learn the discriminative features to classify them. We have also reported limitations to this otherwise general improvement of CapsNets over ConvNets, their improvement in performance is significant but has a limit that we observed for the more complex TUPAC dataset at 1% (5.5K training samples).

Classification tasks where the global spatial structure plays a role can better exploit the advantages of CapsNets (DIARETDB1).

One of the disadvantages of routing-by-agreement is that is slower than regular backpropagation, CapsNet with 8.2M parameters take about the same training time per epoch than Baseline with 35.4M (a ResNet-50 has 25.6M parameters). These architectures lack purposed layers, e.g. batch normalization, that could help to ease the convergence. Depending on the number of classes, Cap-sNet and Baseline need between 1-3 minutes per epoch, while LeNet runs in 1-2 seconds.

Also, when visualizing the images reconstructed through the encoder-decoder branch (Fig. 3), we observe that they are blurry, especially for medical datasets with complex backgrounds. The fully-connected layers of this branch seem to be good enough to regularize the parameter optimization but lose a lot of information. Our future work includes replacing these layers with deconvolutions to get a better insight into the learned latent space.

We recommend the use of capsule networks for medical datasets where the structure is important and patterns appear in different parts of the input images, as it is for retina. Our results confirm that they perform better than standard ConvNets for the limited amount of data, at least of the order of 10k. Another potential application would be the detection of rare diseases or segmentation due to the high performance under class-imbalance. cial support through the "la Caixa" INPhINIT Fellowship Grant for Doctoral studies at Spanish Research Centres of Excellence, "la Caixa" Banking Foundation, Barcelona, Spain. The authors would like to thank Nvidia for the GPU donation and Aurélien Geron for his tutorial and code on Capsule Networks.

Fig. 1 :

 1 Fig. 1: Comparison of the flow and connections of ConvNets vs. CapsNets. Eq. (1)shows the difference between the sigmoid and squashing functions. Eq. (2) is a weighted sum of the inputs (ConvNets use bias). In CapsNets, c ij are the coupling coefficients. In (3), ûj|i is the transformed input to the j-th capsule/neuron. In CapsNets, the input from the i-th capsule is transformed with the weights W ij . While in ConvNets, the raw input from the previous neuron is used.

Fig. 2 :

 2 Fig. 2: Mean F 1 -score and standard deviation (4 runs) for different amounts of training data. Solid line: CapsNet, dotted line: Baseline, and dashed line: LeNet.

Fig. 3 :

 3 Fig. 3: Test input images and their reconstructions.

  Mean F1-score reported for different class-imbalance scenarios.

	Training Data		1%			5%			10%		50%
		LeNet Base. CapsNet LeNet Base. CapsNet LeNet Base. CapsNet LeNet Base. CapsNet
	TUPAC16	0.822 0.784	0.809	0.872 0.835	0.872	0.890 0.852 0.898 0.908 0.903 0.923
	DIARETDB1 0.870 0.847 0.875 0.877 0.852 0.893 0.883 0.863 0.907 0.895 0.854 0.908
	Fashion-M.	0.759 0.749 0.772 0.841 0.817 0.846 0.856 0.847 0.866 0.885 0.889 0.896
	MNIST	0.909 0.916 0.943 0.961 0.966 0.975 0.975 0.978 0.985 0.987 0.989 0.992
		(a) Mean F1-score using different amounts of training data.
	Scenario		Balanced			Imbalanced 1		Imbalanced 2
		LeNet	Baseline CapsNet	LeNet	Baseline CapsNet	LeNet	Baseline CapsNet
	TUPAC16	0.914	0.913	0.932	0.881	0.813	0.892	0.905	0.874	0.909
	DIARETDB1	0.895	0.863	0.899	0.869	0.839	0.887	0.889	0.874	0.898
	Fashion-M.	0.899	0.911	0.910	0.890	0.902	0.889	0.871	0.881	0.863
	MNIST	0.989	0.991	0.991	0.988	0.989	0.993	0.985	0.987	0.992
	(b) Data Augmentation		No				Yes
				LeNet	Baseline CapsNet	LeNet	Baseline CapsNet
	TUPAC16		0.904	0.892	0.914	0.914	0.913	0.932
	DIARETDB1		0.883	0.864	0.895	0.892	0.863	0.899
	Fashion-MNIST		0.899	0.911	0.910	0.902	0.911	0.913
	MNIST			0.989	0.991	0.991	0.990	0.993	0.994
		(c) Mean F1-score with and wihtout data augmentation.

Table 2 :

 2 F-1 scores under different data-challenges.

https://www.tensorflow.org/

Acknowledgment. This work has received funding from the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No. 713673. Amelia Jiménez-Sánchez has received finan-


S L Gavrilyuk
email: sergey.gavrilyuk@univ.u-3mrs.fr

H Gouin

V M Teshukov
email: teshukov@hydro.nsc.ru

Bubble effect on Kelvin-Helmholtz' instability

We derive boundary conditions at interfaces (contact discontinuities) for a class of Lagrangian models describing, in particular, bubbly flows. We use these conditions to study Kelvin-Helmholtz' instability which develops in the flow of two superposed layers of a pure incompressible fluid and a fluid containing gas bubbles, co-flowing with different velocities. We show that the presence of bubbles in one layer stabilizes the flow in some intervals of wave lengths.

Introduction

Many mathematical models of fluid mechanics are derived through the approximation of the solution of boundary value problems for Euler equations. For example, equations for bubbly flows are derived as an approximation of the solution of a complex free-boundary problem describing the motion of the mixture of water and gas bubbles [START_REF] Iordanski | On the equations of motion of the liquid containing gas bubbles[END_REF], [START_REF] Kogarko | On the model of cavitating liquid[END_REF], Wijngaarden (1968)). Green-Naghdi's model describing wave motions of a liquid layer of finite depth in a shallow water approximation with account of dispersion effects is another example [START_REF] Green | On the theory of water waves[END_REF], [START_REF] Green | A derivation of equations for wave propagation in water of variable depth[END_REF]). Due to the fact that averaging procedures and asymptotic expansions have been used in the derivation, it is not obvious to decide what boundary conditions are natural for these systems of equations. In this paper, a Hamiltonian formulation of the problem in (t, x)-space is proposed which allows one to find boundary conditions for a general class of models. This class includes a model of bubbly fluid and dispersive shallow water. We use the boundary conditions to study the stability of two co-flowing layers of bubbly and pure fluids. In absence of gravity and capillarity, it is well known that Kelvin-Helmholtz' instability develops for any wave lengths in the flow of two superposed layers of pure incompressible fluids (see, for example, [START_REF] References Drazin | Hydrodynamic Stability[END_REF]). We prove that the presence of bubbles in one layer can produce a stabilizing effect: the flow becomes stable with respect to perturbations of some wave lengths.

2 Governing equations and conditions on moving interfaces

Variation of Hamilton's action

Here we calculate the variation of Hamilton's action for a special class of Lagrangians. Usually this procedure is used only for derivation of governing equations. We obtain below not only governing equations but also boundary conditions at inner contact surfaces.

Let us consider the Lagrangian of the form:

L = L(J, ∂J ∂z , z) (2.1)
where z = t x ≡ z i , (i = 0, 1, 2, 3), t is the time, x is the space variable, J = ρ ρu , ρ is the fluid density and u is the velocity field. Let us calculate the variation of Hamilton's action in the case when the 4-D momentum J verifies the equation of continuity:

Div J ≡ ∂ρ ∂t + div(ρu) = 0 (2.2)
Hamilton's action is defined by

a = Ω Ldz (2.3)
where Ω is a material domain. Considering a smooth one-parameter family of virtual motions

z = Φ(Z, ε), Φ(Z, 0) = ϕ(Z)
(Z stands for the Lagrangian coordinates, ε is a small parameter at the vicinity of zero and z = ϕ(Z) is the real motion), we define the virtual displacements ζ(Z) and the Lagrangian variations δJ(Z) by the formulae:

ζ(Z) = ∂Φ(Z, ε) ∂ε ε=0 , δJ(Z) = ∂J(Z,ε) ∂ε ε=0 (2.4)
Due to the fact that Z = ϕ -1 (z), we consider the variations as functions of Eulerian coordinates and use the same notations ζ(z) and δJ(z) in Eulerian variables. Hamilton's principle assumes that ζ(z) = 0 on the boundary ∂Ω of Ω.

In the following, the transposition is denoted by T . For any vectors a and b we use the notation a T b for their scalar product a • b and ab T for their tensor product a⊗ b. The divergence of the second order tensor A is a covector defined by Div (Ah) = Div (A) h

where h is any constant vector field. Div and div, Grad and ∇ are respectively divergence and gradient operators in the 4-D and 3-D space. The identity matrix and the zero matrix of dimension n are denoted by I n and O n .

In calculations we shall use the equality

δJ = ∂ζ ∂z -(Div ζ) I 4 J (2.5)
which was proved in Gavrilyuk & Gouin (1999). Variation of the Hamilton action is

δa = da dε ε=0 = Ω δL + L tr ∂ζ ∂z dz (2.6) Since δL = ∂L ∂J δJ+tr     ∂L ∂ ∂J ∂z δ ∂J ∂z     + ∂L ∂z ζ and δ ∂J ∂z = ∂δJ ∂z - ∂J ∂z ∂ζ ∂z
we get from (2.5) and (2.6):

δa = Ω ∂L ∂J ∂ζ ∂z -(Div ζ) I 4 J+tr A T ∂δJ ∂z - ∂J ∂z ∂ζ ∂z + ∂L ∂z ζ + L Div ζ dz with A T = ∂L ∂ ∂J ∂z or A T j i = ∂L ∂J i ∂z j (2.7)
For the sake of simplicity the measure of integration will not be indicated. Since for any linear transformation A and vector field v

Div(Av) = (Div A)v+tr(A ∂v ∂z ) we get δa = Ω tr J ∂L ∂J -A T ∂J ∂z ∂ζ ∂z + L - ∂L ∂J J Div ζ+tr A T ∂δJ ∂z + ∂L ∂z ζ = = Ω tr J ∂L ∂J -A T ∂J ∂z ∂ζ ∂z + L - ∂L ∂J J Div ζ+ +Div A T δJ -Div A T δJ + ∂L ∂z ζ ,
and finally

δa = Ω tr J ∂L ∂J -A T ∂J ∂z ∂ζ ∂z + L - ∂L ∂J J Div ζ + Div A T δJ - -Div A T ∂ζ ∂z -(Div ζ) I 4 J + ∂L ∂z ζ = = Ω tr J ∂L ∂J -A T ∂J ∂z -JDiv A T ∂ζ ∂z + + L - ∂L ∂J J + DivA T J Div ζ + Div A T δJ + ∂L ∂z ζ = = Ω Div J ∂L ∂J -A T ∂J ∂z -J Div A T ζ + A T δJ - -Div J ∂L ∂J -A T ∂J ∂z -J Div A T ζ+ +Div L - ∂L ∂J J+ DivA T J ζ -Grad L - ∂L ∂J J+ DivA T J T ζ+ ∂L ∂z ζ Let us denote δL δJ ≡ ∂L ∂J -Div A T = K T . Then δa = Ω ∂L ∂z -Div J K T -A T ∂J ∂z + L -K T J I 4 ζ+ (2.8) + ∂Ω N T J K T -A T ∂J ∂z + L -K T J I 4 ζ + N T A T δJ
with N T = (-D n , n T ), D n denotes the surface velocity and n is the space unit normal vector. Virtual displacements vanish at the boundary ∂Ω and the surface integral is zero. The volume integral yields the equations of motion in conservative form as in Gavrilyuk & Gouin (1999).

In the case when fluid tensorial quantities are discontinuous at the inner interface Σ, expression (2.8) becomes

δa = Ω ∂L ∂z -Div J K T -A T ∂J ∂z + L -K T J I 4 ζ+ (2.9) + Σ N T J K T -A T ∂J ∂z + L -K T J I 4 ζ + N T A T δJ
where the jump through Σ is denoted by [ ].

Governing equations and boundary conditions

We explicit the governing equations and the inner boundary conditions for the Lagrangian

L = 1 2 ρ |u| 2 -W ρ,
• ρ , where

• ( ) = d dt = ∂ ( ) ∂t + u T ∇ ( ) (2.10)
Such a Lagrangian appears in the study of wave propagation in both shallow water flows with dispersion and bubbly flows (a complete discussion of these models is given in [START_REF] Gavrilyuk | Generalized vorticity for bubbly liquid and dispersive shallow water equations[END_REF]). We get

∂J ∂z =      ∂ρ ∂t ∂ρ ∂x ∂j ∂t ∂j ∂x      , ∂L ∂J = - |u| 2 2 - ∂W ∂ρ + 1 ρ ∂W ∂ • ρ (∇ρ) T u, u T - 1 ρ ∂W ∂ • ρ (∇ρ) T Matrix (2.7) becomes A T = ∂L ∂ ∂J ∂z =       - ∂W ∂ • ρ 0 T - ∂W ∂ • ρ u O 3       .
We get

Div(A T ) = - ∂ ∂t ∂W ∂ • ρ -div ∂W ∂ • ρ u , 0 T , (2.11) K T = - |u| 2 2 - δW δρ + 1 ρ ∂W ∂ • ρ (∇ρ) T u, u T - 1 ρ ∂W ∂ • ρ (∇ρ) T , N T J = ρ n T u -D n , N T A T = - ∂W ∂ • ρ n T u-D n , 0 T .
We study the case when Σ is a contact surface and, consequently, N T J = 0, N T A T = 0 T . We see that the surface integral in (2.9) vanishes if the boundary conditions :

L -K T J ≡ [p] = 0 (2.12) are fulfilled with p ≡ L -K T J.
Vanishing the volume integral in (2.9) and using relations (2.11), we obtain the governing equations in the form

∂ρ ∂t + div(ρu) = 0, (2.13) ∂ρu ∂t + div ρuu T +pI = 0, p ≡ ρ ∂W ∂ρ - ∂ ∂t ∂W ∂ • ρ -div ∂W ∂ • ρ u -W = ρ δW δρ -W.
A complete set of boundary conditions at the contact interfaces for the class of models considered is:

n T u-D n = 0, [p] = 0 (2.14)
Equations (2.13) have been obtained earlier in [START_REF] Gavrilyuk | Media with equations of state that depend on derivatives[END_REF] by means of the classical method of Lagrange multipliers. Notice that this method did not give jump conditions (2.12). System (2.13) is reminiscent of Euler equations for compressible fluids. The term p in equations of motion (2.13) and boundary conditions (2.14) stands for the pressure. Nevertheless, p is not a function of density as in the case of barotropic fluids: it depends also on material time derivatives of the density. Thus, p is not usual thermodynamic pressure.

In general, it implies that the boundary conditions are not just consequences of conservation laws. For example, for Korteweg-type fluids where the stress tensor depends on the density gradient, the boundary conditions contain normal derivatives of the density at interfaces [START_REF] Seppecher | The limit conditions for a fluid described by the second gradient theory: the case of capillarity[END_REF], Gouin & Gavrilyuk (1999)). As we have shown here, the boundary conditions at inner interfaces in bubbly fluids do not depend on normal derivatives of the density: they involve only tangential derivatives.

If the flow domain is shared between two domains consisting of pure ideal fluid and bubbly fluid, all previous calculations are still valid and boundary conditions (2.14) can be extended to such a case.

Exact statement of the problem

Now, we propose to use conditions (2.14) to study Kelvin-Helmholtz' instability of a parallel flow between two rigid walls of a layer of bubbly fluid in contact with a layer of pure incompressible fluid.

Consider the 2-D flow of two superposed finite layers of inviscid fluids in the channel: -∞ < x < +∞, 0 < y < H 0 . The first one is governed by the Euler equations :

div(u 1 ) = 0, ∂u 1 ∂t + u T 1 ∇ u 1 + ∇p 1 ρ 10 = 0, -∞ < x < ∞, 0 < y < h(t, x)
(2.15) Here ρ 10 is the constant density of pure fluid, u 1 = (u 1 , v 1 ) T is the velocity field and p 1 is the pressure. In the domain h(t, x) < y < H 0 governing equations (2.13) for the second fluid are

∂ρ 2 ∂t + div(ρ 2 u 2 ) = 0, (2.16) ∂u 2 ∂t + u T 2 ∇ u 2 + ∇p 2 ρ 2 = 0, p 2 = ρ 2 ∂W ∂ρ 2 - ∂ ∂t ∂W ∂ • ρ 2 -div ∂W ∂ • ρ 2 u 2 -W, W = W (ρ 2 , • ρ 2 )
Here ρ 2 , u 2 = (u 2 , v 2 ) T and p 2 are the average density, velocity and mixture pressure, respectively; W = W (ρ 2 ,

• ρ 2 ) is a given potential. For the flow of compressible bubbles having the same radius, potential W has the form (see, for example, [START_REF] Gavrilyuk | Large amplitude oscillations and their thermodynamics for continua "with memory[END_REF] and [START_REF] Gavrilyuk | Generalized vorticity for bubbly liquid and dispersive shallow water equations[END_REF]):

W (ρ 2 , • ρ 2 ) = ρ 2 c g ε g (ρ g ) -2πnρ l R 3 • R 2 (2.17)
where ε g is the internal energy of the gas in bubbles, R is the bubble radius, ρ g is the gas density, ρ l = const is the density of carrying phase, c g = const is the bubble mass concentration, n is the number of bubbles per unit mass. Potential W is the difference between the internal energy of gas and the kinetic energy of fluid due to radial bubble oscillations. The bubble radius and the bubble density are functions of the average density ρ 2 :

4 3 πR 3 = 1 n 1 ρ 2 - 1 -c g ρ l , ρ g = c g 1 ρ 2 - 1 -c g ρ l -1
It can be shown that governing equations (2.16) with potential (2.17) coincide exactly with classical equations of bubbly fluids [START_REF] Kogarko | On the model of cavitating liquid[END_REF] and van Wijngaarden ( 1968)). In particular, the role of "equation of state"

p 2 = ρ 2 δW δρ 2 -W, W = W (ρ 2 , • ρ 2 ) (2.18)
plays the Rayleigh-Lamb equation which governs the radial oscillations of a spherical bubble:

R •• R + 3 2 • R 2 = 1 ρ l (p g -p 2 ), p g = ρ 2 g dε g dρ g (2.19)
The equivalence between (2.19) and (2.18) with W given by (2.17) has been proved, for example, in [START_REF] Gavrilyuk | Large amplitude oscillations and their thermodynamics for continua "with memory[END_REF]. Dissipation-free model (2.16)-(2.17) assumes that the sliding between components is negligible. Moreover, this model is valid only when the volume fraction of bubbles is very small. Notice that a detailed description of the bubble interaction has been recently done by [START_REF] Russo | Kinetic theory of bubbly flow I: Collisionless case[END_REF] and [START_REF] Herrero | On the motion of dispersed balls in a potential flow: a kinetic description of the added mass effect[END_REF] in the case of rigid bubbles, and [START_REF] Smereka | A Vlasov equation for pressure wave propagation in bubbly fluids[END_REF]), Teshukov & Gavrilyuk (2002) in the case of compressible bubbles, by using a kinetic approach. At the rigid walls y = 0 and y = H 0 , the vertical components of the velocity are equal to zero

v 1 (t, x, 0) = v 2 (t, x, H 0 ) = 0.
In accordance with the previous Section, we prescribe the following boundary conditions at the contact interface y = h(t, x):

h t + u 1 h x = v 1 , h t + u 2 h x = v 2 , p 1 = p 2 (2.20)
3 Linear stability problem

Linearization

Consider the following main parallel flow of two fluids:

u 1 = u 10 = const, v 1 = 0, p 1 = p 0 = const, 0 < y < h 0 , h 0 = const (3.1) u 2 = u 20 = const, v 2 = 0, p 2 = p 0 = const, ρ 2 = ρ 20 = const, h 0 < y < H 0 (3.
2) Small perturbations denoted by ′ :

u 1 = u 10 + u ′ 1 , v 1 = v ′ 1 , p 1 = p 0 + p ′ 1 , ρ 2 = ρ 20 + ρ ′ 2 , u 2 = u 20 + u ′ 2 , v 2 = v ′ 2 , p 2 = p 0 + p ′ 2
satisfy the linearized system of equations

u ′ 1x + v ′ 1y = 0, ρ 10 D 1 u ′ 1 + p ′ 1x = 0, ρ 10 D 1 v ′ 1 + p ′ 1y = 0, 0 ≤ y ≤ h 0 (3.3) D 2 ρ ′ 2 +ρ 20 (u ′ 2x +v ′ 2y ) = 0, ρ 20 D 2 u ′ 2 +p ′ 2x = 0, ρ 20 D 2 v ′ +p ′ 2y = 0, h 0 ≤ y ≤ H 0 .
(3.4) We use the notations

D i = ∂ ∂t + u 0i ∂ ∂x , i = 1, 2 In (3.4) p ′ 2 = a 2 ρ ′ 2 + b 2 D 2 2 ρ ′ 2 with a 2 = ρ 20 ∂ 2 W ∂ρ 2 2 (ρ 20 , 0), b 2 = -ρ 20 ∂ 2 W ∂ • ρ 2 2 (ρ 20 , 0).
Here a is the equilibrium sound velocity, and b is a characteristic wave length depending on the bubble radius and gas volume fraction. We suppose that

∂ 2 W ∂ρ 2 ∂ • ρ 2 (ρ 20 , 0) = 0
This condition is obviously fulfilled for potential (2.17). The coefficients a and b calculated in equilibrium R = R 0 and p 2 = p 0 are :

a 2 = - dp g dτ (τ 0 ) 1 N 0 ρ 20 , b 2 = 1 4πN 0 R 0 ρ l ρ 20 (3.5)
For dilute mixtures we can simplify these expressions by replacing ρ 20 by ρ l :

a 2 = - dp g dτ (τ 0 ) 1 N 0 ρ l , b 2 = 1 4πN 0 R 0 = R 2 0 3α 0
Here N 0 = ρ 20 n is the number of bubbles per unit volume, R = R 0 is the equilibrium bubble radius, α 0 = τ 0 N 0 is the volume fraction of gas, τ 0 = 4 3 πR 3 0 is the bubble volume and p g = p g (τ ) is the gas pressure in bubbles expressed as a function of bubble volume. The equilibrium sound speed a is usually small with respect to the gas sound velocity. Expressions (3.5) can be obtained directly after linearization of Rayleigh-Lamb's equation for bubbles (2.19).

Boundary conditions (2.20) at y = h 0 are

h ′ t + u 10 h ′ x = v ′ 1 , h ′ t + u 20 h ′ x = v ′ 2 , p ′ 1 = p ′ 2 (3.6)
Let us consider the normal modes of the linear problem (3.3), (3.4), (3.6):

u ′ 1 = U 1 (y) exp(i(kx -ωt)), v ′ 1 = ikV 1 (y) exp(i(kx -ωt)), (3.7) p ′ 1 = P 1 (y) exp(i(kx -ωt)), h ′ = H exp(i(kx -ωt)) ρ ′ 2 = R 2 (y) exp(i(kx -ωt)), u ′ 2 = U 2 (y) exp(i(kx -ωt)), v ′ 2 = ikV 2 (y) exp(i(kx -ωt)), p ′ 2 = P 2 (y) exp(i(kx -ωt))
By substituting into the linearized system we get the system of equations for unknown amplitudes

U 1 + V 1y = 0, ρ 10 (u 10 -c)U 1 + P 1 = 0, (3.8) -ρ 10 k 2 (u 10 -c)V 1 + P 1y = 0, (u 20 -c)R 2 + ρ 20 (U 2 + V 2y ) = 0, ρ 20 (u 20 -c)U 2 + P 2 = 0, -ρ 20 k 2 (u 20 -c)V 2 + P 2y = 0, P 2 = a 2 -b 2 k 2 (u 20 -c) 2 R 2
with the following boundary conditions

V 1 (0) = 0, V 2 (H 0 ) = 0, P 1 (H 0 ) = P 2 (H 0 ), V 1 (h 0 ) u 10 -c = V 2 (h 0 ) u 20 -c = H
Here c = ω k is the phase velocity.

From equations (3.8) we obtain the eigenvalue problem for pressure amplitudes:

P 1yy -k 2 P 1 = 0, 0 < y < h 0 , (3.9) P 2yy -λ 2 k 2 P 2 = 0, h 0 < y < H 0 , P 1 (h 0 ) = P 2 (h 0 ), P 1y (h 0 ) ρ 10 (u 10 -c) 2 = P 2y (h 0 ) ρ 20 (u 20 -c) 2 , P 1y (0) = P 2y (H 0 ) = 0 where λ 2 = 1 - (u 20 -c) 2 a 2 -b 2 k 2 (u 20 -c) 2 .
It follows from (3.9) that the eigenvalues c are solutions of the equation

th(kh 0 ) ρ 10 (u 10 -c) 2 + λ th(λk (H 0 -h 0 )) ρ 20 (u 20 -c) 2 = 0 (3.10)
Next we assume that both layers are thin: k (H 0 -h 0 ) ≪ 1, kh 0 ≪ 1, and dispersion relation (3.10) is simplified into

(H 0 -h 0 ) ρ 20 1 (u 20 -c) 2 - 1 a 2 -b 2 k 2 (u 20 -c) 2 + h 0 ρ 10 (u 10 -c) 2 = 0 (3.11)

Study of dispersion relation

In case u 20 = u 10 , equation (3.11) reduces to the quadratic equation

(1 + A)a 2 -(1 + A) b 2 k 2 + 1 (u 20 -c) 2 = 0, with A = h 0 ρ 20 (H 0 -h 0 )ρ 10
.

which has only real roots. It means that the flows with equal velocities are stable.

Let us consider the general case when u 20 = u 10 . Equation (3.11) can be rewritten as a polynomial of fourth degree. The stability needs that all four roots are real. To study this problem, it is convenient to rewrite the dispersion relation (3.11) in the form F (z) + A ≡ 0 (3.12)

with

F (z) = (1 -N z) 2 1 - d 2 z 2 -1 , (3.13) z = ad u 20 -c , d = 1 bk , N = M d , M = u 20 -u 10 a
Here M is similar to the Mach number, d is the dimensionless length of perturbation wave. Obviously, z is real if and only if c is real. The derivative of (3.13) is

F ′ (z) = -2N d 2 (1 -N z)ϕ(z) (3.14)
where

ϕ(z) = 1 d 2 + 1 -N -1 z (z 2 -1) 2 . (3.15)
Four cases must be considered

1 • 0 < N < 1, d 2 ≤ N -2 -1 (3.16) 2 • 0 < N < 1, d 2 > N -2 -1 (3.17) 10 3 • N > 1, d 2 > 16N 2 27 8N 2 -6 1 - 3N -2 4 + 8N 2 -9 (3.18) 4 • N > 1, d 2 ≤ 16N 2 27 8N 2 -6 1 - 3N -2 4 + 8N 2 -9 (3.19)
Each case is illustrated in Figures 1234. For greater convenience, the roots z i of (3.12) are shown as the intersection points of the two graphs η = F (z) and η = -A.

The results are the following:

In case 2 the dispersion relation defined by (3.12)-(3.13) has four real roots for A not large.

In case 3 the dispersion relation has four real roots for "intermediate" values of A.

In other cases the dispersion relation has two real roots and two complex roots.

Precise formulations are given in the following Propositions.

Proposition 1 If inequalities (3.16) are satisfied, equation (3.12) has two real roots and two complex roots for any positive A (Figure 1).

Proposition 2 If inequalities (3.17) are satisfied, the equation ϕ(z) = 0 defined by (3.15) has one root z 01 on interval I 2 . Equation (3.12) has four real roots for 0 < A ≤ -F (z 01 ) (Figure 2a) and two real and two complex roots for A > -F (z 01 ) (Figure 2b).

Proposition 3 If inequalities (3.18) are satisfied, the equation ϕ(z) = 0 defined by (3.15) has two roots z 02 and z 03 on interval I 2 . Equation (3.12) has four real roots for -F (z 02 ) < A ≤ -F (z 03 ) (Figure 3b), two real and two complex roots for A < -F (z 02 ) (Figure 3a) and A > -F (z 03 ) (Figure 3c).

Proposition 4 If inequalities (3.19) are satisfied, equation (3.12) has two real roots and two complex roots for any positive A (Figure 4).

The proofs are given in Appendix.

It is convenient to rewrite conditions (3.16)-(3.19) in terms of b, k and M :

1 • : 0 < M < 1, k < √ 1 -M 2 bM (3.20) 2 • : 0 < M < 1, √ 1 -M 2 bM < k < 1 bM or M > 1, k < 1 bM (3.21) 3 • : 1 bM < k < 1 bd 0 (M ) (3.22) 4 • : 1 bd 0 (M ) < k (3.23)
Here d 0 (M ) is the root of the equation f (d, M ) = 0, where

f (d, M ) = d 2 - 16M 2 27d 2 8 M d 2 -6 1 - 3d 2 4M 2 + 8 M d 2 -9 (0 < d 0 (M ) < M ). This equation has a unique solution d = d 0 (M ) in a domain 0 < d < M because of the properties f d (d, M ) > 0, f (d, M ) → -∞ as d → 0 and f (M, M ) = M 2 > 0.
Inequalities (3.20)-(3.23) are illustrated in Figure 5. For any given flow parameters M and b, one can find a sufficiently large wave number k satisfying inequality (3.23). It means that the flow is always unstable with respect to perturbations with sufficiently short wave lengths (as it follows from Proposition 4, the dispersion relation has complex roots). The flow with subsonic relative velocity is stabilized in an intermediate interval of wave lengths if the pure liquid layer is thin enough with respect to the bubbly layer. In the case of supersonic relative velocity a similar stabilization is observed for long waves (see Proposition 2 and inequalities (3.21)). For another intermediate interval of wave lengths (3.22) the stabilization of perturbations is attained for layer depths of the same order (see Proposition 3).

Conclusion and discussion

We have derived from Hamilton's principle of stationary action governing equations and boundary conditions at the contact interfaces in bubbly fluids. It has been shown that the dynamic condition on the interface reduces to the continuity of the average pressure. By using the boundary conditions derived, we have studied the Kelvin-Helmholtz instability of two superposed layers of a pure incompressible fluid and a bubbly fluid. We have shown that in contrast to the case of two incompressible fluids when the instability develops for any length of perturbations (if the gravity, capillarity or compressibility are not taken into account), the presence of bubbles can stabilize the flow in some range of perturbation wave lengths.

The stabilizing effect is due to the following reason. The development of classical Kelvin-Helmholtz instability leads to the appearance of wave-like bulges at the interface between two fluids. The presence of bubbles in a fluid (i.e. new interfaces) permits one to transform a part of the energy responsible for the bulge formation into the energy of radial oscillations of bubbles.

The viscosity effect on the flow stabilization is an important issue. The method of viscous potential flows developed by [START_REF] Joseph | Breakup of a liquid drop suddenly exposed to a high-speed airstream[END_REF] in the analysis of Rayleigh-Taylor and applied by [START_REF] Funada | Viscous potential flow analysis of Kelvin-Helmholtz instability in a channel[END_REF] in the analysis of Kelvin-Helmholtz, could be used here. Acknowledgement Teshukov, V.M. & Gavrilyuk, S.L. 2002 Kinetic model for the motion of compressible bubbles in a perfect fluid, European J. Mech. B/Fluids 21, 469-491.

van Wijngaarden, L. 1968 On the equations of motion for mixtures of liquid and gas bubbles. J. Fluid Mech. 33, 465-474.

Appendix

Real roots of equation (3.12) belong to the intervals

I 1 : -(1 + d 2 ) 1/2 < z < -1, I 2 : 1 < z < (1 + d 2 ) 1/2
Assuming that u 20 > u 10 and the wave number k is positive, we see that the relative Mach number M and the constant d are positive. Equation (3.12) has always at least two real roots for any positive values of d, A and N because on the real axis F (z) tends to -∞ as |z| tends to 1 + 0, and F (z) = 0 for |z| = (1 + d 2 ) 1/2 . Additional two real roots can appear in subsets of real axis where F ′ (z) changes sign. Notice, that the interval I 1 contains exactly one root of (3.12) because, obviously, F ′ (z) < 0 for z < -1 (see (3.14)). We prove below that, depending on the flow parameters, equation (3.12) has one or three real roots on interval I 2 .

First, let us show that ϕ

′ (z) > 0 on I 2 if 0 < N < 1. If √ 3/2 < N < 1, then the derivative ϕ ′ (z) = 3z 2 -4N z + 1 N (z 2 -1) 3
vanishes at points z 1 < z 2 , where

z 2,1 = 2N 3 1 ± 1 - 3N -2 4 . (6.1)
The inequality ϕ ′ (z) > 0 is fulfilled because in this case the roots z i satisfy the inequalities 0 < z 1 < z 2 < 1, and the interval where ϕ ′ (z) < 0 has no common points with I 2 . For 0 < N < √ 3/2 the roots z i of the equation ϕ ′ (z) = 0 are complex and ϕ ′ (z) > 0 on I 2 . Proof of Proposition 1. It follows from inequalities (3.16) that 1-N z > 0 for z ∈ I 2 . Let us show that F ′ (z) > 0 on this interval. It was proved above that ϕ

′ (z) > 0 on I 2 if 0 < N < 1. Taking into account that ϕ( d 2 + 1) = √ d 2 + 1 N d 4 N d 2 + 1 -1 < 0 (6.2)
we obtain inequalities ϕ(z) < 0, F ′ (z) > 0 on interval (1, √ d 2 + 1). Then the equation (3.12) has only one root on I 2 . Proposition 1 is proved.

Proof of Proposition 2. In this case N -1 ∈ I 2 , F ′ (z) > 0 for 1 < z < N -1 and F (z) has a local maximum at z = N -1 , where F ′ (N -1 ) = 0,

F (N -1 ) = 0 . Taking into account that ϕ ′ (z) > 0 , ϕ( √ d 2 + 1) > 0 and ϕ(N -1 ) = (d 2 + 1)N 2 -1 d 2 (N 2 -1)
< 0, we see that there exists the unique point z 01 on I 2 , z 01 ∈ (N -1 , √ d 2 + 1), such that ϕ(z 01 ) = 0, F ′ (z 01 ) = 0. At this point F (z) takes a local negative minimum: F (z 01 ) < 0. It is obvious that equation (3.12) has three roots on I 2 only if 0 < A < -F (z 01 ) (see Figure 2a). Proposition 2 is proved.

Proof of Proposition 3. When inequalities (3.18) are satisfied, the function ϕ(z) is positive in a neighbourhood of the end points of I 2 (ϕ(z) → ∞ as z → 1 + 0, ϕ( √ d 2 + 1) > 0) and 1 -N z < 0 . The function F (z) can be non monotone only in the case when ϕ(z) takes negative values and, consequently, ϕ ′ (z) vanishes at some point of I 2 . It follows from (6.1) that the roots of the equation ϕ ′ (z) = 0 are such that z 1 < 1, z 2 > 1 for N > 1. It is easy to verify that a stronger inequality z 2 > N is valid for N > 1. Hence only the root z 2 of the equation ϕ ′ (z) = 0 can belong to I 2 . The inequalities

a) d 2 > z 2 2 -1, b) ϕ(z 2 ) = 1 d 2 + 1 -N -1 z 2 (z 2 2 -1) 2 < 0 (6.3)
provide the inclusion z 2 ∈ I 2 and the existence of the roots z 02 , z 03 of the equations ϕ(z

02 ) = F ′ (z 02 ) = 0, ϕ(z 03 ) = F ′ (z 03 ) = 0 satisfying inequalities 1 < z 02 < z 03 < √ d 2 + 1.
The function F (z) has a local maximum at z 02 and a local minimum at z 02 . It means that equation (3.12) has three roots on I 2 if and only if -F (z 03 ) < A < -F (z 02 ) (see Figure 3a). For A > -F (z 02 ) and 0 < A < -F (z 03 ) this equation has only one root (see Figures 3b,3c). Notice that inequality (6.3a) is a consequence of (6.3b) because for N > 1 (z 2 2 -1) 2 N -1 z 2 -1 > z 2 2 -1 (6.4)

Using the identity

(z 2 2 -1) 2 N -1 z 2 -1 = 16N 2 27 8N 2 -6 1 - 3N -2 4 + 8N 2 -9
one can show that inequality (6.3b) is equivalent to the second inequality in (3.18). Proposition 3 is proved.

Proof of Proposition 4. In accordance with (6.4) two cases are possible:

a) (z 2 2 -1) 2 N -1 z 2 -1 > d 2 > z 2 2 -1, b) (z 2 2 -1) 2 N -1 z 2 -1 > z 2 2 -1 > d 2 (6.5)
Inequalities (6.5a) are equivalent to the inclusion z 2 ∈ I 2 (see (6.3)) and the positiveness of minimal value of ϕ(z) on I 2 . Using the inequalities ϕ(z) > 0, 1 -N z < 0 and (3.14) we show that F ′ (z) > 0. Consequently, equation (3.12) has only one root on I 2 .

If inequality (6.5b) is fulfilled, then z 2 / ∈ I 2 and ϕ ′ (z) does not change sign on I 2 . The function ϕ(z) is positive on I 2 , and F (z) is a monotone function. We see that equation (3.12) also has only one root on I 2 . Proposition 4 is proved.

Figure 1 :

 1 Figure 1: The variation of F as a function of z for case 1. The dispersion relation F (z) + A = 0 has only two real roots z i , i = 1, 2 for any positive A. Graph is drawn for N = 0.25 and d = 2 √ 2.

Figure 2 :

 2 Figure 2: The variation of F as a function of z for case 2. The dispersion relation F (z) + A = 0 has four real roots z i , i = 1, 2, 3, 4 for not large values of A and two real roots z i , i = 1, 2 for large values of A. Graph is drawn for N = 3 and d = √ 150.

Figure 3 :

 3 Figure 3: The variation of F as a function of z for case 3. The dispersion relation F (z) + A = 0 has two real roots z i , i = 1, 2 for A small or large, and four real roots z i , i = 1, 2, 3, 4 for A intermediate. Graph is drawn for N = 1.1 and d = 4.45.

Figure 4 : 5 Figure 5 :

 455 Figure 4: The variation of F as a function of z for case 4. The dispersion relation F (z) + A = 0 has only two real roots z i , i = 1, 2 for any positive A. Graph is drawn for N = 2 and d = 10.5