## Abstract

The formation of pores in anodic aluminium oxide films is treated with a model equation. The model treats the oxide layer as a thin viscous liquid in two dimensions. Surface tension on the top boundary, electrostriction due to the external electric field and mass flow through the bottom boundary due to oxide formation are all included. Viscous flow is treated with the creeping flow assumption. The model equation is solved numerically using a Fourier spectral method in space and Adams–Bashforth/Adams–Moulton methods in time. Initial conditions include sinusoidal shapes as well as random shapes. The results show that pores form at the trough of the initial sinusoidal shape. Random shapes get smoothed before forming pore structures with spacing different than predicted by linear theory.

## 1. Introduction

Anodic oxide films grown on passivating metals at high driving force may form a hexagonal array of pores tens of nanometers in diameter that extend from the film surface to within a few nanometers of the metal–oxide interface. This type of coating accounts for an important process technology for protection of aluminium from abrasion and corrosion. The earliest models of pore formation were based on local dissolution of the oxide accelerated by the electrochemical field driving film growth [1–3]. The emergence of a pore array in anodic aluminium oxide (AAO) has also been attributed to mechanical stresses originating with mismatches in the specific volume between the metal substrate and the oxide film [4] or to ionic migration at the interface [5]. Stability analyses based on these driving forces indicate that either mechanism can lead to instability [3–8].^{1}

More recently, tracer studies have shown that the oxide layer ‘flows’ during anodizing [10,11]. A number of investigators have attributed pore formation to this oxide flow, and have indicated that the flow is driven by stresses in the film that arise from electrostriction or from the larger volume of the oxide in comparison with that of the metal [12–21]. Generally, these previous investigations posit that oxide flows from the thin layer at the pore bottoms outwards to the space between pores.

Among the stress-driven flow models, there are several broad views of the mechanism causing flow. In one view, the flow is driven by the electrostriction force alone [12–16]. Under this mechanism, fluctuations in the local film thickness are amplified as the current and hence the field and electrostriction force are greater where the film is thinner. In a previous paper, we applied linear stability analysis to this model [16]. The results show that an interval of wavelengths is unstable, and this interval has a maximum representing the fastest growing mode. The wavelength for this mode depends on the film thickness at which the instability occurs, and so the fastest mode is not unambiguously predicted by the model. However, Vega *et al.* [18] performed experiments in which the initial film thickness was controlled at low field strength before application of a large destabilizing potential. Their experimental pore spacings are in good agreement with the wavelength of the most unstable mode from our previous work [16]. Also, in this previous work, we proposed a mechanism for re-stabilization of the thin film under the pore bottom. As the film is thinned locally under the electrostriction instability, the rate of production of oxide increases until it reaches a balance with the electrostriction force and the viscous resistance to flow. The model predicted a ratio of pore diameter to barrier layer thickness of 0.2, substantially different from the observed value of about 1.2 [1]. Evidently, the forces acting at the pore bottom are more complicated than our simple linear model.

A second view of pore formation holds that stresses induced by the volume of oxide formation contribute to the observed instability [19–21]. Houser & Hebert [19] considered the steady-state form of the film with pores and compared it with a conduction model for the current distribution. They found that two elements are required to obtain a consistent result: a high field conduction law and the presence of a space charge in the oxide. They also treated the viscous flow of oxide in the film [19,20]. They did not specify the origin of the force driving flow, but they gave a quantitative assessment of various candidates including electrostriction.

Recent theoretical and experimental work by Hebert and co-workers [21–24] provides support for other mechanisms: instability based on oxide dissolution or the efficiency of oxide formation rather than oxide flow. In their paper on instability, they show that depending on the formation efficiency, film formation may be stable, may lead to instability with disordered pore formation or lead to instability with ordered pore formation [22]. In experimental work on stress during oxide growth, they have identified several sources of stress, both tensile and compressive [23,24]. Based on a current interruption technique, they find that the electrostriction force, which vanishes upon interruption of the current, is smaller than residual stresses that relax more slowly through diffusional processes. They also found that the electrostriction force is accurately modelled by the method used below [23]. Here, we focus on electrostriction, guided by the fact that our previous stability results with electrostriction produce approximately the experiment pore spacing. We ignore the other mechanisms in the present work, and plan to treat these in a future effort.

The initial linear phase of pore formation with electrostriction appears to be adequately described by linear stability analysis. However, the development of the pores in the nonlinear regime and their re-stabilization into a quasi-steady geometry is not as well understood. The objective of this work is to examine the nonlinear phase in detail and provide insight into the formation of pore geometry and spacing. The approach is to treat a nonlinear evolution equation that approximately models the pore development. The evolution equation assumes the oxide is a thin viscous two-dimensional layer that is subject to the effects of surface tension at the interface between oxide and acid, electrostriction due to the applied voltage and mass flux through the bottom due to oxide formation. Similar evolution equations that model the dynamics of thin liquid films have been treated previously, as reviewed by Oron *et al.* [25], including such effects as heat transfer and evaporation, which are not important here. Thin liquid films with evaporation and melting can also develop pores very similar to some of the results given here. These pores can even cause ultimate rupture of the thin film [25]. Oron *et al*. also discuss van der Waals forces, which are important when films are 100–1000 ångström, significantly thinner than the films being considered here; thus van der Waals forces are neglected. All of these previous studies include the effect of surface tension, but have not included electrostriction and oxide formation at the bottom of the film.

The results given below show that an initially sinusoidal interface shape that is unstable according to linear theory develops distinct pores at the trough of the initial oscillation. Random initial interface shapes are smoothed at first by surface tension, then develop a sequence of pores that do not match the equivalent linear evolution. Section 2 provides the development of the mathematical model, including the models for the squeezing force, surface tension and mass flux due to dissolution. Section 3 discusses rescaling and the numerical methods, and provides results for both sinusoidal initial shapes and random initial shapes. Finally, §4 summarizes and provides overall conclusions.

## 2. The mathematical model

The configuration is shown schematically in figure 1. An aluminium substrate is submerged in an acid solution. An electric potential (voltage) is applied at the top of the acid (*ϕ*_{1}) and at the bottom of the aluminium (*ϕ*_{2}), creating an electric field through the experiment. This ultimately causes a film of oxide material to form on the surface of the aluminium. After an initialization of the process, the oxide is presumed to exist as a thin layer. There are then two interfaces, the aluminium–oxide interface and the acid–oxide interface, both of which may assume complex shapes. The results given below will assume that the oxide is ‘thin’, and for such thin layers, the aluminium–oxide interface may be assumed to remain flat. The complexity of both interfaces is then measured by the thickness of the oxide layer, given by *z*=*η*(*x*,*y*,*t*). Here, *x*,*y* are the Cartesian coordinates in the plane tangent to the substrate, while *z* is normal to the substrate.

### (a) The electric field

The oxide material has a much higher electrical resistance than either the acid or the aluminium. This fact suggests the idealization used here where both the solution and the aluminium are taken to be perfect conductors, while the oxide is a perfect insulator. When an electric field is applied to a conductor, free charge will move through the conductor and accumulate at the boundary of the conductor forming a charged layer. In the chosen configuration, this process will continue until the field induced by the charged layers is exactly equal in strength to the externally applied field but with opposite sense resulting in a zero net electric field. This is assumed to happen in both the acid and the aluminium. As the oxide is assumed to be a perfect insulator, the end result is a charged layer on the top and bottom of the oxide with opposite sense. These charged layers act to mechanically squeeze the oxide film.

In the region outside of the oxide layer, the total electric field is zero and therefore the electric potential *ϕ* is constant, *ϕ*_{1} (constant) in the acid layer and *ϕ*_{2} (also constant) in the aluminium layer. In the oxide layer, the potential is determined using Laplace’s equation,
*ϕ*=*ϕ*_{1},*ϕ*_{2} on top and bottom. The charge distribution is then determined by evaluating
*q* is the charge per unit area, *ϵ*_{0} is the permitivity of free space and ** n** is a unit normal vector.

For complicated displacements of the interface between oxide and acid, a solution to Laplace’s equation must be achieved in order to precisely determine the charge distribution. However, for thin layers (relative to the transverse dimension), Laplace’s equation may be simplified by rescaling
*λ* is the transverse dimension (defined later) and *η*_{0} is the mean initial (constant) thickness of the oxide. Laplace’s equation becomes

### (b) The squeezing force

Now, each interface is considered to be thin charged sheets, with charge density *q*_{1} for the top sheet and *q*_{2} for the bottom sheet (later set *q*_{1}=*q*_{2}=*q*). The electric field generated by a differential patch of the bottom sheet is
*S* is the surface area of the differential patch, *R* is the distance from the differential element and *e*_{R} is a base vector. The total electric field due to the bottom sheet is
*q*_{1} is assumed constant and the integral is performed over an unbounded tangential distance, resulting in
** k** is a unit vector normal to the substrate. Thus, the electric field is perpendicular to the flat aluminium substrate, since the substrate is flat, infinite in the tangential direction, and has constant charge density

*q*

_{2}.

The attracting force between the charged sheets can be obtained by choosing the field point to be the top of the oxide layer. Coulomb’s law gives
*S* is now the area of a patch of the top surface. Define the squeezing stress *τ* to be the force per unit area. Using the above,

### (c) Surface tension

The charged sheet on top of the layer will act as a membrane with the extra feature of the applied Coulomb stress *τ*. At this juncture, we choose to allow deformations only in a single tangential direction, making the problem two dimensional.

The tension in a curved membrane results in a force normal to the membrane that is given by
*τ*:
*z*=*η*, where *σ*_{nn} is the normal stress at the membrane in the oxide and *σ*_{0} is the normal stress in the acid. The value of *σ*_{0} will be taken as constant, implying that the acid is motionless. Note that *τ* is normal to the substrate, while *τ* normal to the membrane. For small deformations, this normal stress boundary condition simplifies to
*z*=*η*_{0}, where

The oxide is presumed to deform plastically as a result of the squeezing force, as discussed in [16]. To include this permanent deformation, the oxide is modelled as a viscous liquid, discussed further in the next section. The stress in the oxide is thus given by
*w* is the normal velocity and *p* is the pressure in the oxide. Furthermore, ∂*w*/∂*z* may also be neglected, as indicated in the next section. Finally,
*z*=*η*_{0}.

A force balance tangent to the membrane gives
*z*=*η*, where *σ*_{nt} is the shear stress in the oxide. For thin layers,
*z*=*η*_{0}.

### (d) Fluid flow relations

It is presumed here that the oxide layer is squeezed sufficiently to result in plastic deformation of the oxide. This deformation is modelled here as the flow of a very viscous fluid. Viscous flow with constant properties is governed by the Navier–Stokes equations. For two-dimensional flow, the stream function *ψ* may be employed, defined by *u*=*ψ*_{z}, *w*=−*ψ*_{x}, allowing the equations to be reduced to
*A*,*B*,*C* and *D* could depend on the tangential distance *x*. On the flat bottom, the no-slip condition requires
*z*=0, resulting in *B*=0. The tangential velocity is then

The volume flux past any tangential position *Q* is
*C* in (2.8) gives
*p*_{z}≈0.

The tangential pressure gradient may be determined then by taking the tangential derivative of the pressure boundary condition, (2.4)

### (e) Mass flux

Dissolution is treated with a mass flux through the flat bottom,
*z*=0, where *ρ* is the density and *m* is the mass flux. This condition determines *A* in (2.6), giving
*m* is proportional to the inverse of the thickness,
*M* is constant. The final amplitude equation is
*F*_{0} and *M* are zero, then this equation becomes a well-known equation for the evolution of a thin viscous film with a free surface. Furthermore, the second term has the same form as a typical diffusion term, except with the opposite sign. This term now represents negative diffusion.

## 3. Results

The amplitude equation (2.15) is treated here with a numerical technique. First, all quantities are made dimensionless. Choose
*L* is the length of the domain and *λ* is a length-scale associated with the initial conditions. Note that the length scale used for *x* is different from that for *t*. This choice allows multiple oscillations to be included in the initial conditions without changing the scaling. The dimensionless equations are
*x*∈(0,1) and

The amplitude equation (3.1) is integrated in time using the second-order Adams–Moulton method for linear terms and the second-order Adams–Bashforth method [26] for nonlinear terms:
*j* is a positive integer indicating the time step.

The side boundaries are taken to be periodic and *η* and *Q* are expanded in Fourier series:
*k* is an integer. This Fourier spectral method has an error that decays faster than any algebraic method (finite difference), as discussed in Boyd [27]. Spatial derivatives are treated by taking derivatives of the Fourier functions directly, which allows the equation to be rearranged into an explicit form:
*Q*_{k} at each time step are evaluated pseudospectrally, where values of *η* and its derivatives are determined at equally spaced gridpoints using an inverse fast Fourier transform (FFT), then *Q* is determined at these same gridpoints and finally transformed to obtain *Q*_{k} with a forward FFT. Details of the pseudospectral method are also given by Boyd [27].

Accumulation of energy at the highest modes is avoided by filtering at each time step using the sequence of filters by Vandeven [28]. This family of filters maintains the spectral accuracy of the Fourier transform. They are increasingly sharp with an increase in the parameter *p*: the value *p*=15 was used for all calculations.

The value *N*=256 is used for all cases, corresponding to 512 Fourier coefficients, and the same number of gridpoints for the pseudospectral evaluation of the nonlinear terms. Higher resolution was found to produce the same results and therefore unnecessary.

### (a) Linear stability

A small amplitude disturbance is unstable for some parameter intervals, as previously reported [16]. For the above model equation, choose
*η*^{′} is a small perturbation. Inserting this into the above equation gives
*c*_{M} is assumed constant, it follows that the average layer thickness increases linearly. The disturbance *η*^{′} is governed by
*η*_{0} is the original layer thickness. Note that the layer is assumed to be growing slowly in this stability analysis such that the average thickness *η*_{0} for short time periods. This assumption requires the growth of the instability to be much faster than the growth of the layer thickness, which is generally true when the layer is unstable. The solution to (3.3) is
*σ*>0, or

If the wavelength of the critical wavenumber is chosen to equal the domain length *L*, then *k*_{cr}=2*π*. Waves longer than *L* are unstable; however, such waves are not included, and this represents the cut-off for instability. Hence, instability exists when
*C*_{T} to be less than this critical value.

### (b) Sinusoidal initial conditions

Choose the initial shape to be a sinusoid given by
*β* must be the inverse of an even integer for this initial shape to be periodic, for example
*C*_{T}, *C*_{M}, *α*) are listed in table 2. With these parameter values, the integer *n*=5 corresponds approximately to the fastest growing sinusoidal mode. Other parameter values have a different value of *n* for the most unstable mode, but otherwise show very similar behaviour.

Figure 2 shows the evolution of the thickness for this case. There are three panels in figure 2*a*–*c*, each panel providing a profile of *η*(*x*) for a fixed time, shown with the solid line. Figure 2*a* is the initial condition, while figure 2*b,c* is later times. The dashed line in each panel is the profile *η*(*x*) for the linear solution, matching the evolution of the *n*=5 mode with the same parameter values. Note that *C*_{M} is set to zero for the evolution of the linear mode.

The solid line in figure 2*b* shows that the ‘troughs’ of the initial condition have plunged to a very small thickness, forming structures that will be called ‘pores’, analogous to the pores that are seen in experiments. Comparing the solid and dashed lines in figure 2*b*, it can be seen that the nonlinear theory predicts that the pores form much faster than the linear theory. Note that the thickness at the pore bottom is small but not zero.

The solid line in figure 2*c* shows that the pore bottoms have formed new dendritic structures at the location of the bottom of the pores. The thickness of the new structure is growing monotonically with time and the width of the original pore is increasing to accommodate the new structure. Overall, the number of dendritic structures appears to be doubling by this mechanism.

The solid line in figure 2*c* also shows that the tops of the original dendrites have begun to ‘flatten’. This flattening appears to be the beginning of a tip-splitting phenomena; however, tip-splitting cannot be confirmed with the present model. The difficulty is that the flattening and splitting process is complex, and is accompanied by the growth of small-scale structures. The small-scale structures have proved to be difficult to resolve in practice, and result finally in failure of the simulation. The small-scale structures only form at the top of dendritic shapes, and are unsteady with a spectrum of time scales including very fast motions. Only very drastic filtering can control these structures with the present model; however, this filtering also effects the larger scales to an unacceptable level. Hence no reliable results are reported here beyond the time when these small-scale structures form. A more sophisticated model is probably necessary to continue a simulation for longer times.

The amplitude of the initial condition *A* has the rather large value of 0.5 in the above case. Extensive testing with smaller values of *A* showed that the amplitude grows in a linear manner while the average thickness of the layer increases. However, eventually the pores appear as above, and the overall pattern of evolution is the same as above.

Figure 3 shows the evolution for the same case as figure 2, except that the initial condition has *n*=2. This longer mode is still unstable, but the growth rate is much slower than the *n*=5 case. Figure 3*a* again shows the initial condition, now just a single oscillation across the domain. Figure 3*b* shows that the trough of this initial oscillation has formed a more complex structure. The length of this new structure is approximately the same as the wavelength of the fastest growing mode. The initial condition does not have any significant energy at the wavelength of the fastest mode, but as the longer mode grows, it is distorted nonlinearly, and the harmonic that matches the fastest mode grows and forms the structures present in figure 3. Interestingly, the ‘trough’ of the initial condition does not form a ‘pore’, mimicking the secondary dendrite formation of the previous case. Also, note that the crest of the initial condition does not show any significant smaller structures, but will eventually be polluted by the fastest growing mode.

### (c) Random initial conditions

The initial film thickness in a physical experiment will not be sinusoidal, and here we consider a more general initial shape. The initial shape is determined with a random number generator.

Figure 4 shows the same case as before, with parameter values listed in table 2. Figure 4*a* shows the complex initial condition. Figure 4*b* shows a later time, and as before, the dashed line is the layer thickness for linear evolution, while the solid line is for nonlinear evolution. Both linear and nonlinear cases use the same initial condition. Figure 4*b* shows that the high-frequency oscillations in the initial condition have been smoothed for both linear and nonlinear cases. This smoothing is due to surface tension. The fastest-growing mode from linear stability is the same mode as before, thus has the same five ‘troughs’ or ‘pores’ as figure 2. The dashed line in figure 4*b* shows that five minima have formed, presumably the beginning of this same fastest-growing mode. The nonlinear case in figure 4*b* only shows four minima.

Figure 4*c* shows that the linear solution has grown such that the pores have penetrated the bottom (there is no mass flux through the bottom of the linear solution to prohibit this). However, there are still five minima. The solid line in figure 4*c* shows a different pattern, with only four minima. Thus, an initial condition along with nonlinear effects can result in a different pore spacing.

Note in figure 4*c* that the linear solution is much larger than the nonlinear, suggesting a different growth rate. The growth rates are compared in figure 5, which shows a time history of the maximum and minimum of the layer thickness for both linear and nonlinear cases. The upper dashed line, which shows the maximum thickness for the linear case, shows the expected exponential growth of the linear mode. However, the upper solid line in figure 5, which corresponds to the maximum for the nonlinear case, appears nearly straight, suggesting linear rather than exponential growth. The lower solid and dashed lines show the minimum thickness, thus providing a measure of the thickness of the pore bottom. Note in figure 5 that the minimum thickness increases at first, as a result of the increase in the average layer thickness. Once the pore begins to form, however, the minimum thickness decreases very rapidly, much faster than the linear case. The pore also develops faster than the previous sinusoidal case. Sometimes this results in the pore bottom penetrating the layer, as in this case.

Figure 6 shows the Fourier spectra for the same three time values of figure 4. Note the log scale for the ordinate. Figure 6*a* shows the spectrum for the initial condition, which is nearly uniform across all scales, similar to white noise. Figure 6*b* shows that at a later time the higher frequencies contain very little energy, due to the effects of surface tension. It is also evident in figure 6*b* that a maximum exists at *k*=5, the most unstable mode. Finally, figure 6*c* shows the spectrum at a much later time, and energy at high frequency has reappeared, due to the effects of nonlinearity. Also in figure 6*c*, the most unstable mode is now nearly an order of magnitude larger than all other frequencies.

Figure 7 shows a profile of the magnitude of each the three terms in the evolution equation (3.1) at a single time value. The individual terms are
*f*_{1}, which can be traced to the effects of the squeezing stress, has consistently the largest magnitude, while *f*_{2}, the surface tension term, is the smallest. However, the three terms have the same order-of-magnitude, thus all three effects are important in the evolution of the film. These trends are very similar at other time values.

Figure 8 is another case with random initial shape and with the same parameter values as in figure 4. Furthermore, the initial shapes, while different, have the same mean and standard deviation. Figure 8*b* shows that the linear solution has developed six minima instead of the five that are expected. Apparently, after surface tension has smoothed the initial shape, there is not enough energy in the fastest growing mode for it to appear significantly, and instead, a nearby mode appears.

The nonlinear solution is much different, as seen in figure 8*c*. The solid line only shows three minima, and only one has formed a pore, although the others will probably continue to form two more. Note that this nonlinear evolution for this case is much slower than the previous case. This can be seen in figure 9, which again shows the time history of maximum and minimum of layer thickness. The maximum thickness again grows approximately linearly, until flattening of the top of the dendrite begins. Simultaneously, the pore formation is delayed by more than twice the time scale of the linear pore formation. Note that, in all cases, the rapid pore formation and the ‘flattening’ of an adjacent dendrite occur together. This fact is probably related to the mass flux through the bottom, which increases dramatically during pore formation. This increased mass gets pushed to either side of the pore, resulting in a local increase in thickness on both sides of the pore. We speculate that it is this local thickening that is responsible for the small-scale oscillations discussed above.

Figure 10 shows the spectra of the film thickness for the case of figure 8. Figure 10*a* shows again that the initial spectrum has significant energy at all frequencies, similar to white noise. Note here, however, that the value of |*η*_{5}|, the most unstable mode, is significantly smaller in size than nearby frequencies. Figure 10*b* shows the intermediate time where surface tension has eliminated higher frequencies. Here, the frequency of the most unstable mode corresponds to the largest value of |*η*_{k}|, but this value is not dominant. Figure 10*c* shows that, at a later time, higher frequencies are again gaining energy. In this case, however, several lower frequencies have larger values of |*η*_{k}| than the most unstable mode |*η*_{5}|. Thus nonlinear effects have caused nearby lower frequencies to dominate. This feature appears when the initial spectrum has low energy at the most unstable frequency, as here, giving the nearby frequencies the opportunity to gain energy at a rate that is apparently faster than the linear growth of the most unstable mode.

Other values of the parameters gave similar results, with the primary difference being the time for pore and dendrite development. For example, small values of *C*_{M} show the same basic pore development, only the pores appear much sooner, as shown in figures 11 and 12. Figure 11 shows the time evolution of the maximum and minimum film thicknesses for several values of *C*_{M}. Clearly, the minimum thickness, which is the thickness at a pore bottom, approaches zero much sooner with small values of *C*_{M}. Other than the time, the evolution is approximately the same. This is shown in figure 12 with a single profile of *η* for three different cases. The time value for each case was chosen to coincide with the time when the minimum thickness was near zero. Note that all figure 12*a*–*c* shows that a single pore has developed, and on either side of the pore the film has become thicker. While these three profiles are not identical, it does indicate that approximately the same pore development occurs with different *C*_{M} values, with the primary difference being the speed of development.

## 4. Summary and conclusion

Previous laboratory experiments [1] have shown that aluminium oxide films will develop pores when grown with high driving voltages. Other experiments suggest that this pore formation may be due to significant plastic deformation of the oxide that occurs as a result of electrostriction. Thus, the results here have assumed that the oxide is a thin layer experiencing electrostriction and modelled as a very viscous fluid. Recent stability analysis [16] of this configuration with the inclusion of surface tension at the top of the layer has shown that an instability exists and the fastest growing mode approximately matches the scale of the observed pores. Nonlinear pore growth is treated here with an approximate evolution equation. Along with nonlinear effects, mass flow through the bottom of the layer has been added to model oxide formation.

The nonlinear evolution equation is solved numerically. Results with a sinusoidal initial layer thickness tuned to match the most unstable linear mode show that the troughs rapidly decrease in thickness and form pores, with pore spacing matching the wavelength of the initial mode. A new dendritic mode then appears at the pore bottom, which grows rapidly, doubling the number of pores and dendrites. Results with a much longer initial wavelength show the smaller most unstable wavelength grows in the trough of the oscillation, forming pores only within the trough. An initial shape consisting of a broad spectrum is first smoothed by surface tension, followed by pore growth. The final spacing of the pores depends on the initial interface shape, and can be close to the spacing predicted by linear theory. However, a different random initial shape that has the same mean and standard deviation can produce a much different pore spacing if the initial energy in the most unstable mode is small. The growth rate of both the pore and dendrite thickness can be much different from linear growth rates.

The driving force for pore formation is assumed here to be electrostriction only. Recent experiments [22] indicate that other stresses may also be important, and the magnitude of these other effects may even be stronger than electrostriction. The present results show that electrostriction by itself can lead to pores, similar to experiments. However, it is possible that these other mechanisms will also lead to such pore formation. Other mechanisms will be included in future studies.

## Ethics

The work did not involve any active collection of experimental data, only mathematical theory and computer simulations of material behaviour.

## Authors' contributions

J.M. carried out the modelling and numerical aspects of the work, participated in the design of the study and drafted the manuscript. D.B. conceived the study, participated in design of the study, participated in the modelling aspects of the work, and helped draft the manuscript. Both authors gave their final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

No funding has been received for this article.

## Footnotes

- Received December 21, 2016.
- Accepted April 28, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.