## Abstract

The electrically induced steady-state solutions of a gas bubble in a dielectric liquid under the action of a steady electric field are considered using the leaky dielectric model. Representing the shape deformation by a sum of spherical harmonics, it is shown that for a given parameter set there exists a critical value of the ratio of the electric to surfaces stresses beyond which no steady states exist, thus implying bubble instability and possible fragmentation. Previous studies imply that bubble instability can only be achieved if either the dielectric constant or the conductivity of the gaseous contents of the bubble is large. We show that on accounting for compressibility of the bubble, no such restriction applies for bubble instability. Below these critical values, multiple steady states are found. It is shown that a more approximate model, which assumes that the bubble is a prolate ellipsoid, can be used to represent the results for weak electric fields, but cannot be used for the prediction of the critical value of the strength of the electric field beyond which no steady-state exists.

## 1. Introduction

The determination of the possible steady states of a dynamical system under external force plays an integral part in accessing the stability and even existence of given phenomena. Here we examine the steady states of a gas bubble under a temporally and spatially uniform electric field. Introducing the dimensionless parameter *W*, which represents the ratio of electric to surface tension stresses, we show that as *W* increases, the non-existence of any steady states has important implications on bubble behaviour. In particular, we address here the practical issue of determining for what applied external fields a bubble can exist in a stable state or undergo disintegration.

Garton & Krasucki (1964) showed that for drops (bubbles) to become unstable, the dielectric constant of the liquid (gas) has to be much larger than that of the surrounding fluid. Their experiments also showed that beyond this critical value, conical points develop at the poles of the drops with either a fluid jet or a fine mist being ejected from the tip of these conical points. Employing a boundary integral formulation, Sherwood (1988) showed through numerical simulation that if the conductivity of a drop is sufficiently large, it may break up in a different manner, by deforming into separate blobs. Grigor’ev *et al.* (1999) concluded that as the gas permittivity is normally much smaller than that of the liquid, only the second mechanism identified by Sherwood (1988) could possibly be applicable to bubbles, with the conditions being met by assuming a good conductivity inside the bubble (and hence a constant interior potential) as the result of a large mobility of charge being present on the interface between the two phases. (An earlier study by Cheng & Chaddock (1984) using a free energy argument was unable to find critical conditions beyond which bubbles would become unstable. They attributed this result to the fact that the prolate spheroid shape they assumed for the bubble was too approximate, but still expected the bubble to disintegrate if sufficient elongation was incurred due to the application of the electric field.) Under their chosen assumptions, Grigor’ev *et al.* (1999) obtained the appropriate pressures at the pole and at the equator of a bubble deformed into a prolate spheroid by an applied electric field, using the approach of Taylor (1964). They then demonstrated that beyond critical values of the ratio of the electric to surface stresses, no steady-state solutions exist, thus implying bubble instability. For values less than this, two solution branches were found. Although compressibility inside the bubble was accounted for within their model by assuming bubble pressures satisfy the ideal gas law, the assumption of good conductivity inside the bubble is key to the predicted instability in Grigor’ev *et al.* (1999).

In this work, we revisit the problem of the stability of a gas bubble in a liquid, without requiring the conductivity of the bubble to be large. We shall allow for both phases to be conducting (most, but not all of the results presented are for non-conducting bubbles) and permit the inducement of charge on the interface between these phases. Thus a leaky dielectric model is employed (see Saville 1997 for a review); a perfect dielectric corresponds to a leaky dielectric with zero conductivity. Two approaches are presented: in the first, the bubble deformation is represented by a sum of spherical harmonics, whereas in the second a dipole approximation in a prolate spheroidal geometry is employed. The latter approach has parallels with a theoretical analysis conducted by Sherwood (1988) (based on earlier work by Stratton 1958) for the case wherein both fluids are perfect dielectrics and in the absence of charge. It is found that in this limiting case, our prolate spheroidal analysis is consistent with the theoretical result derived by Sherwood (1988). Our prolate spheroidal analysis for compressible, non-conducting bubbles predicts one solution branch and bubble stability for all values of *W* considered.

If instead the bubble deformation is not assumed to be ellipsoidal, but is expressed in terms of spherical harmonics (such that the interfacial conditions are satisfied more accurately on the entire bubble surface, not just the poles and equator), a critical value of *W* is determined beyond which no steady states are found. This result is found to be valid even if the conductivity inside the bubble is zero. Further, setting the conductivities in both phases to zero, i.e. the dielectric limit, was found not to have any apparent effect on our observations. For values of *W* less than this critical value, multiple steady states are found. Shaw *et al.* (in press) applied domain perturbation analysis within a spherical harmonic framework and represented the bubble deformation by the summed effect of a substantial number of Legendre polynomials, taking the analysis to second order in the small interaction terms. Their work predicted the existence of two main solution branches for the steady states, confirming the previous work of Spelt & Matar (2006). The origin of the multiple steady states predicted by the latter work was a consequence of assuming the contents of the bubble were compressible. Spelt & Matar (2006) considered the violent collapse of bubbles under electric fields in the context of sonoluminescence studies and thus calculated the pressure inside the bubble using van der Waals’ model. The same model is used both in the work of Shaw *et al.* (in press) and here. Intriguingly, Shaw *et al.* (in press) also found indication of a third solution branch. For , it was found that the predicted shape mode amplitudes for this branch were beyond the validity of the perturbation analysis, but as *W* was increased, predicted amplitudes become permissible. However, the perturbation approach applied by Shaw *et al.* (in press) was constrained to small values of *W*, and the valid *W* values for the observed third branch were near the limits placed on this parameter. In this paper, we relax the small values of *W* assumption and derive appropriate equations to second order in the perturbation analysis, exploring the third branch in more detail, as well as whether there is a maximum strength of the electric field for stable bubbles to exist, for a range of parameter values. We note that a van der Waals’ model is also used in our prolate spheroidal analysis, but as noted above, only one solution branch for non-conducting bubbles is found in this case.

## 2. Problem statement and method of solution

A gas bubble is assumed to exist in a dielectric liquid and is subject to a temporally and spatially uniform electric field **E**_{0} orientated along the *z*-axis of a Cartesian coordinate system. (Note, the orientation of the coordinate system is arbitrary in the absence of an electric field.) Steady-state solutions will be sought of the governing equations for the fully transient problem, by setting all time derivatives to zero. Although the translation of the bubble at constant velocity does represent a possible steady-state solution, attention is restricted here to stationary bubbles; in a separate analysis, this latter possibility was allowed for, but no non-zero velocity solutions were found.

The resultant electric fields in the liquid and inside the bubble are assumed to be solenoidal and irrotational, permitting an electric potential-based approach. Employing the subscript ‘in’ to represent quantities within the bubble, we have that (**E**,**E**_{in})=(∇*ψ*(**x**),∇*ψ*_{in}(**x**)), where the respective electric potentials *ψ* and *ψ*_{in} satisfy appropriate Laplace equations and the position coordinate **x** is measured relative to the bubble centre. The electrostatic boundary conditions (Saville 1997; Spelt & Matar 2006) at the interface *S* between the gas and liquid phases are the jump condition for the normal component of the electric field, the continuity of the tangential component of the electric field and the balance of free charge at the bubble interface (Spelt & Matar 2006), and are given, respectively, by
2.1
2.2
and
2.3
Within these boundary conditions, and denote, respectively, the relative permittivities of the liquid and gas phases, the permittivity of a vacuum, *κ* and *κ*_{in} the respective conductivities of the liquid and gas phases, the outward-pointing unit normal vector to the bubble surface and *q* the interfacial charge density.

The remaining conditions applied at the bubble surface are imposed on the normal and tangential components of the total stress. For the latter boundary condition, it is assumed that the tangential stress is zero. The normal stress boundary condition is
2.4
where *σ* is the interfacial surface tension, **M** the Maxwell stress tensor and the background equilibrium pressure in the liquid. The pressure inside the bubble *p*_{g} is modelled by van der Waals’ equation
2.5
The term *V*_{0} is the volume of a spherical bubble in the absence of the electric field, whereas the corresponding pressure inside the spherical bubble is denoted by *p*_{g0}; *γ* is the ratio of specific heats, *b* the van der Waals’ excluded volume (where, for most cases we use *b*=*V*_{0}/8.86^{3}, which is representative for argon bubbles; Hilgenfeldt *et al.* 1996) and *V* the bubble volume on application of the electric field.

In terms of the spherical polar coordinate system **x**=(*r*,*θ*,*Φ*) defined relative to the bubble centre, for small, axisymmetric deformation the bubble surface *S* is defined as
2.6
where the spherical radius of the bubble is represented by *R*, and *P*_{k}(*μ*) is the Legendre polynomial of degree *k*. The amplitudes *a*_{k} of the respective Legendre polynomials are scaled by the small parameter *ϵ* (≪1), to ensure small deformation. These amplitudes are often termed ‘shape mode’ amplitudes.

Thus, the general bubble volume *V* is determined by evaluating the volume integral , which from equation (2.6) gives
2.7
where *b*_{1}=3*b*/4*π* and *R*_{0} is the bubble radius in the absence of the electric field.

Assuming the bubble is uncharged before the electric field is applied, the electric potentials in the liquid and inside the bubble are (Spelt & Matar 2006) given, respectively, by 2.8 and 2.9 The free charge induced on the bubble surface by the application of the electric field is Fourier decomposed as 2.10

The determination of the coefficients *A*, *B*, *A*_{k}, *B*_{k}, occurring in equations (2.8) and (2.9), together with the subsequent generation of the differential equations the spherical bubble mode, the shape mode amplitudes and the charge density components satisfy, is achieved by applying the surface boundary conditions. To aid the analysis, terms are non-dimensionalized in the manner of Spelt & Matar (2006), resulting in the non-dimensional parameters
2.11
where *W* is the dielectric Weber number, *Π*_{0} a non-dimensional pressure, *λ* denotes the ratio of dielectric constants and *K* and *K*_{in} are dimensionless conductivities.

Within this work, we include terms to *O*(*ϵ*^{2}) in a perturbation analysis similar to the unsteady work of Doinikov (2004) for *W*=0. We allow *W* to be *O*(1) and derive appropriate equations to second order in the perturbation parameter. This contrasts with the analysis of Shaw *et al.* (in press), who restricted attention to cases where *W*=*O*(*ϵ*). (Also note that our analysis includes surface tension contributions to *O*(*ϵ*^{2}), whereas Doinikov’s corresponding contributions are to *O*(*ϵ*) only; these terms are found to be significant in the work considered below.) Thus, many more terms enter into the analysis at *O*(*ϵ*^{2}) through the action of the electric field. Consequently, the pertinent equations we determine for the steady-state problem have all been generated using the symbolic manipulation package MAPLE. The procedure by which this was done is given below. Note, in the following when we state that orthogonality is exploited we mean that the resultant expressions are multiplied by and integrated with respect to *θ*, and thus in the process we exploit the orthogonality properties of Legendre polynomials. Also, we truncate all the expansions at the finite number *k*=*N*. The impact of this truncation on the steady-state analysis will be accessed.

The first step of the procedure is to substitute the truncated expansions for the electric potentials (2.8) and (2.9), together with the equation for the bubble surface (2.6) into the normal and tangential components of the electric field interfacial boundary condition. In the case of the former boundary condition, a truncated form of equation (2.10) also requires substitution. Exploiting orthogonality and eliminating the coefficients *B* and *B*_{k} gives suitable expressions for the coefficients *A* and *A*_{k}. Next, the truncated expansions (2.6) and (2.8)–(2.10) are substituted into the unsteady balance of free charge boundary condition, orthogonality is exploited and equations for the interfacial charge density components are obtained. The coefficients *A* and *A*_{k} are then eliminated from this set of equations using the results from the first step.

Application of the normal component of the total stress boundary condition yields the equations governing the time evolution of the spherical oscillations and shape mode amplitudes. At this stage, the coefficients *D*_{j}, *A* and *A*_{k} are eliminated from this set of equations. Finally, the steady-state terms within this resultant set of equations, together with the appropriate corresponding terms in the interfacial charge density equations, are extracted yielding our pertinent steady-state system of equations.

To ensure the correctness of this automated procedure, we checked both the full equations and the corresponding steady-state equations against those presented by Shaw *et al.* (in press) for the case of *W*=*O*(*ϵ*) (which in itself was validated against the inviscid, incompressible free oscillation work of Shaw (2006)) and complete agreement was found. Further comparisons with previous work are made in the following sections.

## 3. Results

All of the results presented in this section are in terms of non-dimensional variables. A typical result for the spherical radius and the aspect ratio *χ* of a bubble as a function of *W* is shown in figure 1, whereas the corresponding variation in the amplitude of the leading order shape mode is displayed in figure 2*a*. Along one solution branch, at low values of *W*, the radius decreases, whereas the amplitude of the *P*_{2} shape mode increases with increasing *W*, for the starting conditions (*W*,*R*,{*ϵ**a*_{k}})=(0,1,0), *k*∈[2,*N*]. Along the other solution branch, the bubble radius also decreases as *W* increases, but in this case the *P*_{2} shape mode amplitude decreases. Inspection of the results reveals that as , the lower branch is consistent with the asymptotic expressions presented by Spelt & Matar (2006). Also the lower solution branch corresponds to that identified by the same name (but for different parameter sets and *W*≤0.3) by both Spelt & Matar (2006) and Shaw *et al.* (in press). The upper branch corresponds to what Shaw *et al.* (in press) term the ‘third’ solution branch. This branch was not found by Spelt & Matar (2006), and therefore must be attributed to the inclusion of terms that are nonlinear in the shape mode amplitudes. The ‘upper’ branch in these two previous works for the case considered in the present figure corresponds to unacceptably large shape deformations and is therefore not shown here.

As noted in §1, Spelt & Matar (2006) predicted the existence of two main steady-state solution branches, despite their work only being nonlinear in *R*. The source of this nonlinearity was identified as the term *p*_{g}, i.e. the pressure inside the bubble, which is given here by equation (2.7). This term, in turn, makes the normal stress balance a nonlinear equation in *R*, and thus a mechanism for the multiplicity of steady-state solutions is provided by the compressibility of the gas inside the bubble. The same is also true for the present analysis. Mathematically, this nonlinearity in the normal stress balance, and thus possible multiple solutions, still holds for *W*=0. In figure 1, the upper branch is seen to extend to *W*=0, thus predicting two solutions for this value of *W*: one solution corresponds to a sphere and the second to a deformed spheroid of larger volume. From the work of Kushch *et al.* (2002), it is expected that spherical harmonics will only converge with respect to the value of *N* if the aspect ratio of the bubble is sufficiently small. We shall therefore treat the upper branch seen in figure 1 with caution and provide evidence of convergence of the most significant part of the results below, including the turning point. A second, non-spherical solution at *W*=0 cannot outrightly be considered incorrect because of the spherical symmetry of the problem at *W*=0: the orientation of the coordinate system is arbitrary in the absence of the electric field. Intuitively, however, a non-spherical solution at *W*=0 may appear to be incorrect, despite the fact that it would be an admissible second solution of a nonlinear equation, and its origin must therefore briefly be investigated here to ensure the analysis is correct. Its mathematical origin can be found as follows. Consider the non-dimensional steady-state volume oscillation equation for zero electric fields
3.1
where *Π*_{0}−2 is the non-dimensional pressure in the liquid and . As the second solution corresponds to a larger bubble volume, both the and 2/*R* terms decrease in magnitude, but, generally, by different amounts as a consequence of the different dependencies on *R*. The inclusion of the *O*(*ϵ*^{2}) contributions provides a means by which a balance can be returned in equation (3.1). Therefore, it is a combination of the bubble compressibility and the inclusion of higher-order shape mode interactions that results in the prediction of a second *W*=0 solution. We shall investigate the convergence of the results in more detail below.

The most important trend in the results shown in figure 1, and the main subject of this paper, is that the identified branches are seen to turn over at a critical value of *W*, which is denoted henceforth by *W*_{c}. No other solution was found for the cases presented here beyond this turning point. As noted above, a separate solution branch is known to exist (Spelt & Matar 2006) but, on that branch, the shape modes are usually very large, especially when *Π*_{0}=*O*(1) or larger, and in the present case (as we have already observed) would not be admissible. We also see from figure 1*b* that the shape distortion becomes significant even on the presumably stable branch when *W*_{c} is approached from below. We conclude therefore that no steady state is expected to be observed for *W*>*W*_{c}.

Sample bubble shapes are shown in figure 3 for *W*=1.7 on the lower and upper branches, respectively, and at *W*=*W*_{c}. It is seen that on the upper solution branch, the curvature in the (*r*,*z*)-plane is relatively small, and one may speculate that a small perturbation may cause the curvature in the (*r*,*θ*)-plane to lead to rupture of the bubble through surface tension. Note, for this figure, (*r*,*θ*,*z*) is a cylindrical polar coordinate system.

The convergence behaviour of the results is investigated in figure 2*a* where the impact of the chosen value of the truncation parameter *N* on the resultant steady-state predictions for the shape mode amplitude *ϵ**a*_{2} as a function of *W* is studied. We expect this to be the dominant shape distortion mode for our chosen starting conditions. For values of *W* less than *W*_{c}, the lower branch shows little dependence on *N*, but the upper branch and the critical point *W*_{c} clearly do. Convergence of the critical value *W*_{c} is illustrated. Evidently, it is important to truncate at a rather large value of *N*, particularly in the context of obtaining the critical value for *W* beyond which the bubble is predicted to be unstable. However, even values of *N* appear to monotonically converge to approximately the same limit as that for odd values, but from opposing directions, and we note for future reference that the results in figure 2*a* for *N*=10 and 11 give a reasonably accurate estimate of the value of *W*_{c}. The solution along the upper branch does not obey this rather fortunate convergence behaviour, and very large values of *N* are expected to be required in order to obtain an accurate estimate there. In figure 2*b*, larger values of *N* could be used, because for *λ*=0 no charge resulted. The results on both sides of the critical region are seen to be very close, and to converge.

Note, Dr Thiele (2009, personal communication), employing continuation methods (AUTO2000), has checked our solutions for the case of *N*=4 in figure 2*a* and found exact agreement.

Further details of the convergence of the bubble radius, the amplitude of the first six even shape modes and the corresponding interfacial charge density components, at the critical value of *W*, are given in table 1 for the case of figure 2*a*. In table 2, the convergence behaviour of the coefficients of the electric field potential expansions, evaluated at the bubble equator (i.e. *μ*=0) and for criticality, are given. As the tangential electric field boundary condition (2.2) is *ψ*=*ψ*_{in} on the bubble surface, only the coefficients *A* and *A*_{k} are shown.

We have tested the sensitivity of the convergence properties to the manner in which the series in each equation are truncated. In the first method, as we have already noted, the expansions (2.6) and (2.8)–(2.10) are truncated at *k*=*N* and the subsequent analysis retains terms to *O*(*ϵ*^{2}). In an alternative approach, the cut-off value for the series expansions is increased, generating a higher number of shape mode equations, but to the same order in *ϵ*. We then set all shape mode amplitudes higher than the original truncation value *N* to zero. The reason for this is to test whether the nonlinear interaction of any Legendre modes, where at least one is larger than *N*, modifies the coefficients of the shape mode terms less than or equal to *N*. It was found that with this alternative manner of truncation, the convergence appears to be monotonic, i.e. as *N* is increased, the value of *W*_{c} found in the results presented here is approached from above, with, e.g. the result for *N*=8 being similar to that for *N*=9, etc. In what follows, then, we continue to use the first truncation technique.

It is found that for the chosen starting conditions, all the calculated steady-state solutions have zero odd shape mode amplitudes and even interfacial charge density components. The even shape mode amplitudes *ϵ**a*_{k} for *k*>2, the leading order charge density component , together with the higher-order odd charge density components corresponding to figures 1 and 2*a* are shown in figure 4. The amplitudes of the given shape modes and charge density components are seen to decrease in magnitude as the modal number increases. Both sets of distributions reflect the *ϵ**a*_{2} against *W* variations as opposed to the *R* variations in that the turning point is singular. For the spherical bubble radius, the slope ∂*R*/∂*W* appears to be continuous (cf. figure 1).

This difference in the nature of the turning point was observed for all cases studied. To summarize the results, the value of *W*=*W*_{c} at the turning point is shown in figure 5 as a function of *Π*_{0} and *λ*. Given the convergence behaviour of the results in figure 2, we present results for *N*=10 and 11 in this parametric study. As *W* is assumed to be *O*(1) in the analysis, we limit the results to cases where *W*_{c} is not very large. Note, however, figure 5*a* indicates that at sufficiently low values of *Π*_{0}, the critical value of *W* increases rapidly, i.e. bubbles are relatively stable when the imposed pressure is low compared to surface tension.

The results for *W*_{c} as a function of *Π*_{0} can be well represented by empirical fits of the form , where *E* and *F* are constants, but the exponent *k* is found not to be universal, and to depend, e.g., on the value of *λ*. Overall, we can conclude that the turning point is indeed observed in a wide range of cases. Additional tests were performed to check whether a different value of *γ* (=1.4) or *K*_{in} (=0.5) would lead to markedly different behaviour. In both cases, similar results with a turning point were found. Again we emphasize the importance of this critical point as beyond our results indicate that the bubble is unstable and thus will either split or disintegrate.

To determine whether there is a degree of compressibility below which *W*_{c} becomes very large, we investigated further the effect of the value of *γ* on the critical value, a summary of results being shown in table 3. With increasing *γ*, the value of *W*_{c} is seen to decrease, whereas *R* tends towards one, implying a decrease in compressibility.

Finally, we accessed whether the value chosen for the excluded volume parameter in van der Waals’ equation has any major influence on the resultant observations. In terms of the non-dimensional parameter , a number of values for in the range 2–20 were considered. Only small variations were observed.

## 4. A simplified model for the turning point

We have found the simplest case for which the model still predicts a turning point to be for *γ*=1, *λ*=0, *b*=0 with *N*=4. There is no charge distribution in this case because *λ*=0. The relative simplicity of this case allows us to investigate the conditions required for a turning point in more detail, as follows.

The model reduces in this case to the following set of three equations, corresponding to equations for *R*, *ϵ**a*_{2} and *ϵ**a*_{4}, respectively:
4.1
4.2
and
4.3

Further reduction of these equations by setting *ϵ**a*_{4}=0 and omitting the corresponding equation (4.3) does not give a turning point. However, setting *R*=1 and omitting equation (4.1), a turning point is still found. The resulting solutions for the latter case, for *ϵ**a*_{2}, are shown in figure 6.

To investigate this reduced model further, we have eliminated *ϵ**a*_{4} from equations (4.2) and (4.3) with *R*=1. This gives two solutions: inspection showed that one of these returns unphysical results (with *ϵ**a*_{2}≫1), whereas the other corresponds to the solution found numerically. The result of using the appropriate solution branch for *ϵ**a*_{4} is an equation of the form *f*(*a*_{2};*W*)=0. A turning point is now expected if ∂*f*/∂*a*_{2}=0 (we note here that ∂*f*/∂*W* is regular). This condition was evaluated analytically (it is not reproduced here because it is too long), and solved numerically, both with MAPLE. The resulting values for *ϵ**a*_{2} required for a turning point are indicated by the dashed line in figure 6 (two other branches were found following from this requirement, both of which returned only negative values for the present system, and are therefore not shown). It is seen that the dashed line intersects the solution curve for *ϵ**a*_{2} of equations (4.2) and (4.3) (*R*=1) at the turning point.

We therefore conclude that the turning point arises due to the mutual interaction of the shape modes, predominantly *ϵ**a*_{2} and *ϵ**a*_{4}, and that the volume of the bubble does not play a significant role in this. For larger values of *N*, no doubt the actual location of the turning point is somewhat modified, but its existence persists.

Finally, to note that we have verified that carrying the analysis to *O*(*ϵ*^{3}) gives similar results, and the turning point is still present. In this case, more branches appear for solutions of ∂*f*/∂*a*_{2}=0, but importantly one of these still intersects the solution of the system corresponding to equations (4.2) and (4.3).

## 5. Approximate model for small deformation

We now investigate the extent to which a simplified model can be used to represent the results obtained with the involved model used in the previous section. In this simplified analysis, the bubble shape is prescribed as a prolate spheroid, i.e. there is no other deformation. Essentially, we assume we have a dipole due to the electric field rather than a dipole due to a flow as in the work of Kushch *et al.* (2002). Such an analysis is certainly expected to capture the results for modest values of *W* and would therefore be practically useful. A limitation will be, however, the fact that the turning point was seen to be the result of a rather specific interaction between shape modes in the previous section; a model wherein the shape is assumed may struggle to predict such a turning point. We present the analysis in dimensional form, but again the results are for the corresponding equations in non-dimensional form. Define the prolate spheroidal coordinate system relative to the bubble centre **x**=(*ξ*,*η*,*Φ*). For an electric field oriented along the Cartesian *z*-axis, as in the previous section, the respective electric potentials with regard to this new coordinate system are
5.1
where *d* is the distance between the foci of the bubble and *Q*_{1}(*ξ*) is the Legendre polynomial of the second kind of degree 1. Applying the boundary conditions on the continuity of the tangential component of the electric field and the induced interfacial charge density yields
5.2
where the subscript *S* denotes evaluations on the bubble surface. For completeness, we provide the corresponding result for the interfacial charge density, which follows from the jump condition of the normal component of the electric field across the interface
5.3
where and . Note, in this section the relative permittivities of the liquid and gas phases are denoted by ϵ and ϵ_{in} respectively, while ϵ_{0} denotes the permittivity of a vacuum. For the Legendre polynomial of the second kind of degree 1 and its derivative, in these results we employ the following expressions ():
5.4
and
To eliminate *d*, we make use of the volume of the prolate spheroid. If *a* denotes the semi-major axis, and *b* denotes the semi-minor axis, then (Abramowitz & Stegun 1972)
5.5
where the aspect ratio .

These results are a generalization of those given by Sherwood (1988) in that an interfacial charge density is accounted for here. Application of the normal and tangential components of the electric field boundary conditions in the absence of any charge recovers the results given by Sherwood (1988).

The volume and aspect ratio of the bubble are determined from the balance of normal stresses across the interface. Contributions to this balance arise from surface tension, the Maxwell stress **M** and the difference in pressure in the two fluids. Hence, on *ξ*=*ξ*_{S},
5.6
where the pressure inside the bubble *p*_{g} is still modelled by van der Waals’ equation, and is the background ambient pressure. As the bubble shape is prescribed in this analysis, this condition is only approximately satisfied. We enforce equation (5.6) at the poles and equator of the bubble.

As , then the curvature on *ξ*=*ξ*_{S} is given by
5.7
where *h*_{ξ}, *h*_{η} and *h*_{Ψ} are the coordinate scale factors. The net Maxwell stress is
5.8
Therefore, as for *η*=0, collection of all contributions results in the following balance of normal stress at the interface at *η*=0:
5.9
and, as for *η*=1,
5.10

These two equations can be combined to give 5.11

Equations (5.9) and (5.10) (when converted into dimensionless form) were solved numerically for the unknown bubble volume *V* and aspect ratio χ. A comparison with the results obtained using the approach based on spherical harmonics is presented in figure 1 (the radius *R* corresponds to (3*V*/4*π*)^{(1/3)}). The results are seen to agree very well up to modest values of *W* with the lower branches in figure 1*a*,*b*. However, because the bubble shape is prescribed, and the normal stress interfacial condition is only satisfied approximately, not all features of the results obtained from the full analysis are recovered. Hence, the upper branch in these figures is not found (not even at small values of *W*), and it is predicted that a steady-state bubble shape does exist at large values of *W*.

## 6. Conclusions

The steady-state deformation of a gas bubble in a dielectric liquid subjected to a spatially uniform electric field has been considered in terms of a sum of Legendre polynomials. The results of our analysis indicate clearly that there is a critical value of a dielectric Weber number *W* beyond which no steady states exist and thus bubble disintegration is expected. For values just less than this, our bubbles will exist in higher distorted steady-state shapes with the smoothing effect of surface tension balanced by the electrical forcing. It is shown that approximate analyses based on either a few spherical harmonic modes, or an approximate analysis wherein the bubble is assumed to a prolate ellipsoid, are not capable of capturing the critical dielectric Weber number result. On the other hand, such models are useful, as they represent the results of the full analysis rather well for the lower solution branch at low values of *W*.

## Acknowledgements

The authors would like to acknowledge financial support from EPSRC under the grant EP/D5037 1X/1. We also thank Dr Uwe Thiele of Loughborough University for checking the *N*=4 results corresponding to those in figure 2*a*, using continuation methods.

## Footnotes

- Received March 31, 2009.
- Accepted June 30, 2009.

- © 2009 The Royal Society