1. Introduction
In laboratory experiments, forced shear deformation of nominally uniform, partially molten rocks causes melt segregation into high-porosity bands oriented at low angle (
$15^{\circ }$
–
$20^{\circ }$
) to the shear plane (Holtzman et al.
Reference Holtzman, Groebner, Zimmerman, Ginsberg and Kohlstedt2003; Holtzman & Kohlstedt Reference Holtzman and Kohlstedt2007; King, Zimmerman & Kohlstedt Reference King, Zimmerman and Kohlstedt2010). Stevenson (Reference Stevenson1989) predicted the emergence of such bands in a self-reinforcing feedback arising from the porosity weakening of the
$\text{crystal}+\text{magma}$
aggregate, but the angle predicted by this theory was
$45^{\circ }$
(Spiegelman Reference Spiegelman2003), much higher than observed. The low angle of high-porosity bands is widely thought to provide an additional constraint on the rheology of the aggregate, but it has proven challenging to understand. Katz, Spiegelman & Holtzman (Reference Katz, Spiegelman and Holtzman2006) found that non-Newtonian viscosity with a high sensitivity to stress could reproduce the low angle of bands, but King et al. (Reference King, Zimmerman and Kohlstedt2010) subsequently showed that the viscosity of experiments that produce low-angle bands is actually close to Newtonian.
Theory by Takei & Holtzman (Reference Takei and Holtzman2009a ,Reference Takei and Holtzman b ) of anisotropic viscosity under diffusion creep of a partially molten aggregate represents a possible solution. This theory is motivated by observations of the coherent alignment of melt pockets between solid grains under a deviatoric stress (e.g. Daines & Kohlstedt Reference Daines and Kohlstedt1997; Zimmerman et al. Reference Zimmerman, Zhang, Kohlstedt and Karato1999) and of the enhancement of diffusion creep by melt at grain boundaries and triple junctions (e.g. Cooper, Kohlstedt & Chyung Reference Cooper, Kohlstedt and Chyung1989). The melt provides fast pathways for diffusional transport of solid constituents around grains; the alignment of melt with respect to the principal stress directions hypothetically results in anisotropic viscosity of the aggregate (Takei & Holtzman Reference Takei and Holtzman2009a ).
Analysis of the theory of anisotropic viscosity by Takei & Holtzman (Reference Takei and Holtzman2009b ), Katz & Takei (Reference Katz and Takei2013), Takei & Katz (Reference Takei and Katz2013) and Allwright & Katz (Reference Allwright and Katz2014) shows that it introduces qualitatively different behaviour from previous models with isotropic (and even power-law) viscosity. Shear and normal components of stress and strain rate are coupled under viscous anisotropy. As a result of this coupling, a gradient in shear stress becomes a driving force for melt segregation that is not present in the isotropic system. Under Poiseuille flow, melt segregates towards higher-stress regions; under torsional flow, compressive hoop stresses drive the solid outwards and the magma inwards. The mechanics of this ‘base-state’ melt segregation are explained in detail by Takei & Katz (Reference Takei and Katz2013). An experimental test of radial melt segregation in torsional flow by Qi et al. (Reference Qi, Kohlstedt, Katz and Takei2015) shows striking consistency with predictions.
Furthermore, theoretical work has demonstrated that there is a connection between the strength of anisotropy and the angle of high-porosity bands that emerge by unstable growth. This was shown with linearised stability analysis (Takei & Holtzman Reference Takei and Holtzman2009b ; Takei & Katz Reference Takei and Katz2013) and numerical simulations (Butler Reference Butler2012; Katz & Takei Reference Katz and Takei2013), where the strength and orientation of anisotropy are assumed to be known and are imposed a priori. In those static anisotropy calculations, high-porosity bands emerge at low angles to the shear plane only when viscous anisotropy is at or near saturation. This is a rather restrictive condition that may be incompatible with the robust appearance and consistently low angle of bands in experiments (Holtzman & Kohlstedt Reference Holtzman and Kohlstedt2007). However, in numerical simulations that allow anisotropy strength and direction to vary dynamically in space and time (Katz & Takei Reference Katz and Takei2013), band angles are significantly lowered and appear to be less sensitive to the mean strength of anisotropy. These findings raise several basic unanswered questions. Why do the mechanics of viscous anisotropy give rise to low-angle bands? Why is dynamic anisotropy more effective in this regard than static anisotropy? What are the general conditions under which low-angle high-porosity bands should form?
The present paper addresses these questions through a combination of linearised stability analysis and physical reasoning. The crucial, enabling advance is to perform the analysis in a coordinate system that is rotated to align with the porosity bands (rather than with the plane of shear). This drastically simplifies the expressions for growth rate under static anisotropy (Takei & Katz Reference Takei and Katz2013), making them readily interpretable in physical terms. Moreover, it allows us to extend the analysis to dynamic anisotropy in a form that exposes the physical differences from static anisotropy. Finally, the same coordinate rotation clarifies the physical reason for low angles under isotropic non-Newtonian viscosity.
The paper is organised as follows. In the next section, we briefly discuss the non-dimensionalised governing equations and present an anisotropic viscous constitutive model for the two-phase partially molten aggregate. The full nonlinear system is solved numerically in § 3 for static and dynamic cases, to elucidate the questions listed above. The coordinate rotation is introduced and the linearised stability analysis is developed in § 4. In particular, § 4.3 develops an expression for the growth rate of porosity perturbations under the fully dynamic model of § 2. This expression is challenging to understand and so we subsequently consider it under reducing assumptions of static anisotropy (§ 5.1), which includes the simplest case of Newtonian isotropic model. We build on this to explain the full complexity in §§ 5.2 and 5.3. We conclude with a summary and discussion of the results in terms of the motivating questions.
2. Governing and constitutive equations
In the theory of magma–mantle interaction, the macroscopic behaviour of a two-phase aggregate is treated within the framework of continuum mechanics (e.g. Drew Reference Drew1983; McKenzie Reference McKenzie1984). This theory is concerned with the evolution of macroscopic fields, including the volume fraction of melt or porosity
${\it\phi}$
, the velocity of the solid phase
$\boldsymbol{V}$
, the liquid pressure
$P$
(compression positive) and the bulk or phase-averaged stress tensor
${\it\sigma}_{ij}=(1-{\it\phi}){\it\sigma}_{ij}^{S}-{\it\phi}P{\it\delta}_{ij}$
, where
${\it\sigma}_{ij}^{S}$
is the stress tensor of the solid phase (tension positive). Further details of the two-phase flow theory were previously presented (e.g. Rudge, Bercovici & Spiegelman Reference Rudge, Bercovici and Spiegelman2011; Takei & Katz Reference Takei and Katz2013) and are not repeated here.
We proceed directly to the non-dimensional governing equations,











Takei & Holtzman (Reference Takei and Holtzman2009b
) and Takei & Katz (Reference Takei and Katz2013) proposed a model of anisotropic viscosity caused by stress-induced microstructural anisotropy. In partially molten rocks, the melt phase is contained within a permeable network of tubules between grains. The solid matrix is formed by a contiguous skeleton of solid grains. The area of grain-to-grain contact is known as the contiguity. Contiguity is the microstructural variable that determines the macroscopic (i.e. continuum) mechanical properties of the matrix (Takei Reference Takei1998; Takei & Holtzman Reference Takei and Holtzman2009a
). Although the equilibrium microstructure developed under hydrostatic stress has isotropic contiguity, deviations from the equilibrium microstructure have been observed in experimentally deformed, partially molten samples (e.g. Daines & Kohlstedt Reference Daines and Kohlstedt1997; Takei Reference Takei2010). Based on these observations, we infer that, under a differential stress, the grain-to-grain contacts with normals that are parallel to the maximum tensile stress (
${\it\tau}_{3}$
) are reduced in area. Similarly, the areas of those with normals parallel to the maximum compressive stress (
${\it\tau}_{1}$
) are increased. Using a microstructure-based model of aggregate viscosity (Takei & Holtzman Reference Takei and Holtzman2009a
) and a coordinate transformation (Takei & Katz Reference Takei and Katz2013), the constitutive law and the viscosity tensor are














The factor in front of the matrix represents the normalised shear viscosity
${\it\eta}({\it\phi})/{\it\eta}({\it\phi}_{0})$
; it decreases exponentially with increasing melt fraction
${\it\phi}$
, and so
${\it\lambda}$
is called the porosity-weakening factor. We take
${\it\lambda}=27$
based on the experimental results (e.g. Mei et al.
Reference Mei, Bai, Hiraga and Kohlstedt2002). The parameter
$r_{{\it\xi}}$
represents the bulk-to-shear viscosity ratio,
$r_{{\it\xi}}={\it\xi}/{\it\eta}$
, which is assumed to be constant
$(=5/3)$
based on theoretical results by Takei & Holtzman (Reference Takei and Holtzman2009a
) (although see Simpson, Spiegelman & Weinstein Reference Simpson, Spiegelman and Weinstein2010a
,Reference Simpson, Spiegelman and Weinstein
b
). Parameters
${\it\alpha}$
,
${\it\beta}$
and
${\it\Theta}$
represent the magnitude and direction of microstructural anisotropy:
${\it\alpha}$
and
${\it\beta}$
quantify the amplitude of contiguity reduction and increase, respectively, in the directions of principal stress
${\it\tau}_{3}$
and
${\it\tau}_{1}$
;
${\it\Theta}$
represents the angle that the most tensile stress (
${\it\tau}_{3}$
) direction makes with the
$X$
axis of the coordinate system. Using the local differential stress
${\it\tau}_{ij}$
,
${\it\Theta}$
is given by

and
${\it\alpha}$
and
${\it\beta}$
are modelled as
















For simplicity in previous linearised analyses (Takei & Katz Reference Takei and Katz2013; Allwright & Katz Reference Allwright and Katz2014), parameters
${\it\alpha}$
,
${\it\beta}$
and
${\it\Theta}$
were fixed to their initial values. We call this simplifying assumption the static anisotropy model. In contrast, the complete model with stress-dependent direction and magnitude is called the dynamic anisotropy model. Katz & Takei (Reference Katz and Takei2013) discovered a remarkable difference between static and dynamic anisotropy; this difference motivates the present study and is demonstrated in the next section.
3. Numerical solutions
Numerical solutions of (2.1) and (2.2) highlight the difference between results obtained for static and dynamic anisotropy. The solutions are computed with a finite-volume method on a fully staggered grid that is periodic in the
$X$
direction; in this section, the
$X$
axis is taken parallel to the initial flow direction. No-slip impermeable boundary conditions enforce a constant displacement rate of
$\pm {\textstyle \frac{1}{2}}\hat{\boldsymbol{X}}$
on the top and bottom boundaries, respectively. A semi-implicit Crank–Nicolson scheme is used to discretise time, and the hyperbolic equation for porosity evolution is solved separately from the elliptic system in a Picard loop with two iterations at each time step. The solutions are obtained in the context of the Portable, Extensible Toolkit for Scientific Computation (PETSc; Balay et al.
Reference Balay, Buschelman, Gropp, Kaushik, Knepley, McInnes, Smith and Zhang2001, Reference Balay, Buschelman, Gropp, Kaushik, Knepley, McInnes, Smith and Zhang2004; Katz et al.
Reference Katz, Knepley, Smith, Spiegelman and Coon2007). Full details and references are provided by Katz & Takei (Reference Katz and Takei2013).

Figure 1. Comparison of numerical solutions to (2.1) and (2.2) with static and dynamic anisotropy. For both calculations,
$r_{{\it\beta}}=1,\,R=1,\,r_{{\it\xi}}=5/3,\,{\it\phi}_{0}=0.05,\,{\it\epsilon}|{\it\phi}_{1}(\boldsymbol{X})|\leqslant 0.005$
and the domain is discretised into
$600\times 300$
square cells. (a) Porosity field at a strain of 0.75 for a simulation with
${\it\alpha}=1.8$
and
${\it\Theta}=45^{\circ }$
throughout the domain. (b) Porosity field at a strain of 1.25 for a simulation with anisotropy calculated according to (2.3) and (2.4) with
${\it\tau}_{sat}=1,\,{\it\tau}_{offset}=1.5$
. (c) Spectral power binned by wavefront angle
${\it\theta}$
to the shear plane (after Katz et al.
Reference Katz, Spiegelman and Holtzman2006) for the porosity fields shown in (a,b). Each spectrum is normalised by its maximum power. (d,e) Two-dimensional histograms derived from the simulation with dynamic anisotropy at a strain of 1.25 (after Katz & Takei Reference Katz and Takei2013, figure 12). Red dashed lines have a slope given by the ratio of perturbation quantities
${\it\alpha}_{1}/{\it\phi}_{1}$
and
${\it\Theta}_{1}/{\it\phi}_{1}$
from the stability analysis in § 4.
Figure 1 compares solutions with fixed and dynamic anisotropy. In figure 1(a), anisotropy parameters are prescribed as
${\it\alpha}={\it\beta}=1.8$
,
${\it\Theta}=45^{\circ }$
; in figure 1(b), these parameters are calculated cellwise using (2.3) and (2.4), with
$r_{{\it\beta}}=1$
. Both calculations have
$R=1$
(compaction length equal to domain height) and are initialised with the same porosity field,
${\it\phi}(\boldsymbol{X},t=0)={\it\phi}_{0}+{\it\epsilon}{\it\phi}_{1}(\boldsymbol{X})$
. Here
${\it\phi}_{0}=0.05$
and
${\it\epsilon}=0.005$
; and
${\it\phi}_{1}(\boldsymbol{X})$
is a smooth random field with unit amplitude, generated by filtering grid-scale white noise to remove variation at wavelengths below 15 grid cells. Because the growth rate of porosity perturbations differs for fixed and dynamic anisotropy, the simulations are shown at different values of the average simple shear strain
${\it\gamma}$
.
The different orientation of high-porosity features is evident in figure 1(a) and (b): dynamic anisotropy is associated with lower angles. This is quantified by the power spectrum in figure 1(c), where the power from a two-dimensional fast Fourier transform of the porosity field is binned according to the angle between the wavefront and the shear plane (Katz et al.
Reference Katz, Spiegelman and Holtzman2006). Dynamic anisotropy produces a peak at
${\sim}10^{\circ }$
whereas static anisotropy produces a peak at
${\sim}23^{\circ }$
. There is also a high-angle (
${\sim}80^{\circ }$
) peak for static anisotropy (corresponding to features visible in figure 1(a)) that does not survive at large strain. Figure 1(d) and (e) show the covariation of
${\it\alpha}$
and
${\it\Theta}$
with
${\it\phi}$
in figure 1(b); black dotted lines indicate mean values. These means are closely matched with the parameter values used in the fixed anisotropy simulation. It is therefore clear that the difference in the dominant band angle (figure 1
c) arises from the coupling between stress and the variations in
${\it\alpha}$
and
${\it\Theta}$
. What is unclear, however, is the physical explanation for this difference and, indeed, why viscous anisotropy gives rise to bands at angles less than
$45^{\circ }$
to the shear plane at all. We clarify these points below.
4. Linearised analysis with a perturbation-oriented coordinate system
Let
${\it\epsilon}\ll {\it\phi}_{0}$
be the initial amplitude of a porosity perturbation. We express the problem variables as a series expansion about the base state in which the porosity is uniform and equal to
${\it\phi}_{0}$
. We truncate the series after the first-order terms,

The first term of (4.1) with index
$0$
represents a simple shear flow and its associated anisotropy, which is the base-state solution of order one (
${\it\epsilon}^{0}$
), corresponding to the uniform porosity
${\it\phi}_{0}$
. The second term of (4.1) with index 1 represents the perturbation of order
${\it\epsilon}^{1}$
caused by
${\it\epsilon}{\it\phi}_{1}$
. By substituting (4.1) into (2.1), using
$\boldsymbol{{\rm\nabla}}\!\boldsymbol{\cdot }\boldsymbol{V}^{(0)}=0$
, and balancing terms at the order of
${\it\epsilon}^{1}$
, we derive the governing equations for the perturbations as




Following previous studies, the porosity perturbations
${\it\phi}_{1}$
take the form of a plane wave oriented at a given angle to the base-state shear plane. Past workers chose to align the coordinate system with the base-state shear plane, such that the base-state strain-rate tensor has a simple form (e.g. Spiegelman Reference Spiegelman2003). Although the coordinate system was so aligned in the numerical models above, in this section the coordinates are rotated such that the
$Y$
axis is parallel to the wavevector of the initial perturbation, as shown in figure 2. With this choice,
${\it\theta}$
again represents the angle between the perturbation wavefronts and the base-state shear plane. However, in the rotated coordinate system,
${\it\Theta}$
depends on both the direction of
${\it\tau}_{3}$
and the orientation of the bands. The base-state direction of maximum tensile stress
${\it\tau}_{3}$
makes an angle
${\rm\pi}/4$
to the shear plane (2.3) and so, for a coordinate rotation by
${\it\theta}$
, we have
${\it\Theta}_{0}={\rm\pi}/4+{\it\theta}$
(figure 2
a).

Figure 2. Schematic diagrams of the coordinate axes and porosity perturbation. (a) The coordinate system
$(X,Y)$
for the linearised analysis is taken such that the
$X$
axis is parallel to the initial perturbation wavefronts. The shear plane of the base-state simple shear flow is then rotated by an angle
${\it\theta}$
. (b) The base-state normal stress
${\it\tau}_{YY}^{(0)}$
is oriented parallel to the initial perturbation wavevector. (c) The base-state shear stress
${\it\tau}_{XY}^{(0)}$
is parallel to the wavefronts.
In the following part of this section, we give an outline of the linearised approach, which shows that the new coordinate system reduces the complexity of the analysis and exposes the physical mechanisms of perturbation growth. This enables us to clarify the mechanics leading to low-angle bands in complicated problems such as under dynamic anisotropy. We first consider the base-state simple shear flow at the order of
${\it\epsilon}^{0}$
(§ 4.1) and then the linearised governing equations at the order of
${\it\epsilon}^{1}$
(§ 4.2). Finally, in § 4.3, we obtain the growth rate of porosity perturbations
${\it\phi}_{1}$
for the most general case of dynamic anisotropy. The result obtained is used in § 5 to clarify the mechanisms of low-angle band formation.
4.1. Base-state simple shear flow
Using the angle
${\it\theta}$
between the initial perturbation wavefronts (aligned with the
$X$
direction) and the base-state shear plane (figure 2
a), the components of the base-state strain-rate tensor in the rotated coordinate system are

As shown in figure 2(b,c),
${\it\tau}_{YY}^{(0)}$
and
${\it\tau}_{XY}^{(0)}$
represent, respectively, the base-state tensile and shear stresses normal and parallel to the perturbation wavefronts, which play important roles in understanding the growth of these perturbations. Noting that
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(0)}=-\dot{\unicode[STIX]{x1D626}}_{XX}^{(0)}$
, these components are given by

with
${\it\Theta}_{0}={\rm\pi}/4+{\it\theta}$
. In figure 3,
${\it\tau}_{YY}^{(0)}$
and
${\it\tau}_{XY}^{(0)}$
are plotted as a function of angle
${\it\theta}$
. An understanding of their systematics is needed to interpret the results of the stability analysis.

Figure 3. Normal stress
${\it\tau}_{YY}^{(0)}$
(tension positive) and shear stress
${\it\tau}_{XY}^{(0)}$
produced by the base-state simple shear flow as functions of angle
${\it\theta}$
between the shear plane and the
$X$
axis (figure 2
a). We use
$r_{{\it\xi}}=5/3$
here and throughout (Takei & Katz Reference Takei and Katz2013). Each panel is computed with a different set
$({\it\alpha}_{0},{\it\beta}_{0})$
as labelled above each panel.
Considering first the solid curves representing normal stress, we see that, for
${\it\alpha}_{0}={\it\beta}_{0}=0$
(isotropic, figure 3
a),
${\it\tau}_{YY}^{(0)}$
follows the expected pattern of
$\sin 2{\it\theta}$
. It is tensile for
${\it\theta}<90^{\circ }$
and compressive for
${\it\theta}>90^{\circ }$
(Spiegelman Reference Spiegelman2003). However, as
${\it\alpha}_{0}$
increases (for
${\it\beta}_{0}=0$
, figure 3
a–c),
${\it\tau}_{YY}^{(0)}$
becomes negative (compressive) at all angles. The mechanism for this change is twofold. First,
${\it\alpha}_{0}$
decreases the viscous resistance to extension in the
${\it\tau}_{3}$
direction and reduces the maximum tensile stress. This is because increasing
${\it\alpha}_{0}$
reduces
$(\unicode[STIX]{x1D60A}_{YYYY}^{(0)}-\unicode[STIX]{x1D60A}_{YYXX}^{(0)})$
for angles near
${\it\theta}=45^{\circ }$
(figure 4
a, black solid line). Superimposed on this is a compressive stress around
${\it\theta}=0^{\circ }$
and
$90^{\circ }$
that emerges as a consequence of shear strain rate coupled to normal stress via the
$\unicode[STIX]{x1D60A}_{YYXY}^{(0)}$
viscosity (figure 4
a, grey solid line). The product
$\unicode[STIX]{x1D60A}_{YYXY}^{(0)}\dot{\unicode[STIX]{x1D626}}_{XY}^{(0)}$
is negative for all
${\it\theta}$
, motivating us to name this coupling ‘shear strain-induced compression’. This non-trivial result comes from the fact that the stress-induced softening occurs in the tensile (
${\it\tau}_{3}$
) direction, as schematically illustrated in Takei & Katz (Reference Takei and Katz2013, figure 7b).
The effect of increasing
${\it\beta}_{0}$
on
${\it\tau}_{YY}^{(0)}$
is shown down the columns in figure 3. Similar to
${\it\alpha}_{0}$
,
${\it\beta}_{0}$
couples the shear strain rate to compressive normal stress for angles near
${\it\theta}=0^{\circ }$
and
$90^{\circ }$
via the
$\unicode[STIX]{x1D60A}_{YYXY}^{(0)}$
viscosity. In contrast to
${\it\alpha}_{0}$
, however,
${\it\beta}_{0}$
strengthens the aggregate in the
${\it\tau}_{1}$
direction and increases the normal stress amplitude near
${\it\theta}=135^{\circ }$
(figures 3 and 4
b). As a result, the sign change of normal stress caused by
${\it\beta}_{0}$
occurs in a limited range of
${\it\theta}\lesssim 90^{\circ }$
and
${\it\theta}\gtrsim 0^{\circ }$
(figure 3
g).
The shear stress curves in figure 3 (dashed lines) also change with increasing
${\it\alpha}_{0}$
and/or
${\it\beta}_{0}$
. For zero anisotropy in figure 3(a), the shear stress follows
$\cos 2{\it\theta}$
, as expected for coordinate rotation only. Anisotropy does not change the mean of
${\it\tau}_{XY}^{(0)}({\it\theta})$
, as required by symmetry of the stress tensor. Increasing
${\it\alpha}_{0}$
(or
${\it\beta}_{0}$
) has an overall weakening (or strengthening) effect, changing only the amplitude of
${\it\tau}_{XY}^{(0)}$
. This is in contrast to the effect of anisotropy on normal stress, which has a strong dependence on angle
${\it\theta}$
.

Figure 4. Viscosity components from (4.4) as functions of angle
${\it\theta}$
between the shear plane and the
$X$
axis. The curves are computed for the anisotropy parameters (a)
${\it\alpha}_{0}=1$
and
${\it\beta}_{0}=0$
, and (b)
${\it\alpha}_{0}=0$
and
${\it\beta}_{0}=1$
. The thin black lines represent each component for the isotropic case (
${\it\alpha}_{0}={\it\beta}_{0}=0$
).
4.2. Growth of porosity perturbations
For simplicity in this linearised analysis, we consider the case of liquid viscosity
${\it\eta}_{L}=0$
, giving a non-dimensional compaction length
$R\rightarrow \infty$
. In this limit, liquid segregation over any length scale occurs at vanishingly small pressure gradients. Therefore, the pressure gradient terms in (4.2c
) are negligible. Pressure then drops out of the problem and we no longer need to solve (4.2b
) (which has become indeterminate anyway!). This is equivalent to considering only the subset of perturbations with wavelengths much smaller than the dimensional compaction length (e.g. Katz et al.
Reference Katz, Spiegelman and Holtzman2006).
Using the initial wavenumber vector
$\boldsymbol{K}=(0,\boldsymbol{K})$
, the porosity perturbation at time
$t$
is

which accounts for rotation of the wavevector due to advection by the base-state flow (Spiegelman Reference Spiegelman2003). At
$t=0$
, by choice of the coordinate system, perturbations are uniform in the
$X$
direction. Therefore, partial derivatives of the first-order quantities with respect to
$X$
are zero. Using
$\dot{\unicode[STIX]{x1D626}}_{XX}^{(0)}=-\dot{\unicode[STIX]{x1D626}}_{YY}^{(0)}$
,
$\dot{\unicode[STIX]{x1D626}}_{XX}^{(1)}=0$
and
$\boldsymbol{{\rm\nabla}}P_{1}=\mathbf{0}$
, equations (4.2) become




The right-hand sides of (4.6b,c
) represent
${\it\tau}_{YY}^{(1)}$
and
${\it\tau}_{XY}^{(1)}$
, respectively. Since pressure gradients are negligible, these stresses must be spatially uniform for the system to be in balance. Therefore the first-order product of viscosity and strain rate must sum to zero; viscosity reduction associated with porosity perturbations (within the first square brackets on the right-hand side) is compensated by the strain-rate perturbations (within the second square brackets).
To facilitate the physical interpretation of (4.6b,c ), these equations are re-expressed as

with equivalent (‘forcing’) stresses





where
$\widetilde{\unicode[STIX]{x1D60A}}_{YYYY}^{(0)}$
and
$\widetilde{\unicode[STIX]{x1D60A}}_{YYXY}^{(0)}$
are the compliances defined by











4.3. General solution
Equations (4.6) are solved here to obtain an explicit expression for
${\dot{s}}$
for the full model of dynamic anisotropy. The first-order viscosity tensor is written in terms of the porosity and anisotropy perturbations
${\it\phi}_{1}$
,
${\it\alpha}_{1}$
,
${\it\beta}_{1}$
and
${\it\Theta}_{1}$
. The anisotropy perturbations are then expressed in terms of the porosity and strain-rate perturbations
${\it\phi}_{1}$
,
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
and
$\dot{\unicode[STIX]{x1D626}}_{XY}^{(1)}$
. These calculations are sketched in appendix B. The components of the first-order viscosity tensor are then substituted into (4.6b,c
), which are manipulated to solve for
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
and
$\dot{\unicode[STIX]{x1D626}}_{XY}^{(1)}$
as functions of
${\it\phi}_{1}$
. The normal strain rate
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
obtained by this approach is substituted into (4.6a
) to give an expression for the growth rate of perturbations,

with the compliances defined by (4.10) and dynamic factors



The constant coefficients
$p$
and
$q$
that express the sensitivity of the growth rate to dynamic anisotropy perturbations are









5. Physical interpretation in various limits
The growth rate in (4.11) is a general result for the full model presented in § 2 above (with the sole assumption of
$R\rightarrow \infty$
). To build up a physical understanding of this equation, we return to the simpler case of static anisotropy, which includes the simplest case of the Newtonian isotropic model. The static anisotropy model has previously been studied by Takei & Katz (Reference Takei and Katz2013) using linearised analysis. However, the mathematical complexity of their results precluded a detailed mechanical interpretation. A reconsideration using the perturbation-oriented coordinate system enables a physical understanding of the instability mechanism and the rheological control on the dominant band angle. These are needed to understand the more complicated, dynamic model. To facilitate this (in §§ 5.1 and 5.2), we make the simplifying assumption that
${\it\beta}=0$
, i.e. that there is no contiguity increase in the
${\it\tau}_{1}$
direction. The effect of non-zero
${\it\beta}$
is discussed in § 5.3, where we show that its role is minor compared to that of
${\it\alpha}$
, the contiguity decrease in the
${\it\tau}_{3}$
direction.
5.1. Static anisotropy
When
${\it\alpha}_{1}={\it\beta}_{1}={\it\Theta}_{1}=0$
, the mechanical equilibrium conditions (4.6b,c
) are written as









It is evident from (4.6a
) that the normal strain rate
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
causes an increase in the amplitude of porosity perturbations; the shear strain rate
$\dot{\unicode[STIX]{x1D626}}_{XY}^{(1)}$
does not cause the porosity to change. When the viscosity is anisotropic, both the normal and the shear stress drive
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
and hence contribute to perturbation growth. This does not occur under isotropic viscosity.
5.1.1. Instability mechanism in the isotropic system
For an isotropic aggregate (
${\it\alpha}_{0}={\it\beta}_{0}=0$
),
$\unicode[STIX]{x1D60A}_{YYXY}^{(0)}=\unicode[STIX]{x1D60A}_{XYYY}^{(0)}=0$
in (4.7) and (5.1), and the compliance
$\widetilde{\unicode[STIX]{x1D60A}}_{YYXY}^{(0)}$
that couples shear stress to normal strain rate is zero. In this case,
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
is driven only by the base-state normal stress
${\it\tau}_{YY}^{(0)}$
. The growth rate
${\dot{s}}$
is given by

Perturbations are unstable under tensile stress (
${\it\tau}_{YY}^{(0)}>0$
) normal to the perturbation wavefronts, stable under compressive normal stress (
${\it\tau}_{YY}^{(0)}<0$
) and unaffected by shear stress
${\it\tau}_{XY}^{(0)}$
. We therefore term this the tensile stress-induced instability, or tensile instability. When band angle
${\it\theta}$
relative to the simple shear flow is
$45^{\circ }$
, the tensile stress
${\it\tau}_{YY}^{(0)}$
attains its maximum (figure 3
a) and hence the growth rate
${\dot{s}}$
is also at a maximum, as shown in figure 5(a) and by Spiegelman (Reference Spiegelman2003). The occurrence of the tensile instability in a porosity-weakening two-phase aggregate was first predicted by Stevenson (Reference Stevenson1989).
5.1.2. Two instability mechanisms in the anisotropic system
For an anisotropic aggregate (
${\it\alpha}_{0}>0$
and/or
${\it\beta}_{0}>0$
), there is a coupling between shear and normal components via
$\unicode[STIX]{x1D60A}_{YYXY}^{(0)}=\unicode[STIX]{x1D60A}_{XYYY}^{(0)}\neq 0$
. In this case,
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
is forced by both normal stress across perturbations and shear stress along perturbations. The growth rate is

using the compliances given by (4.10). The first term on the right-hand side of (5.3) represents the tensile instability, generalised to the anisotropic aggregate. The second term represents a shear stress-induced instability that does not occur in the isotropic system. The total growth rate
${\dot{s}}$
versus band angle
${\it\theta}$
is plotted in figure 5(a–c) for various anisotropy amplitudes
${\it\alpha}_{0}$
(thick lines). Consistent with previous work, as
${\it\alpha}_{0}$
is increased, the single growth-rate peak splits into two peaks at low and high angles to the shear plane (figure 5
c). Because the lower-angle peak dominates the higher-angle peak after a finite time (Katz et al.
Reference Katz, Spiegelman and Holtzman2006; Takei & Katz Reference Takei and Katz2013), this result means a significant lowering of the dominant band angle by the viscous anisotropy – if the magnitude of anisotropy
${\it\alpha}$
is sufficiently close to saturation (
${\it\alpha}\simeq 2$
).

Figure 5. Characteristics of the static anisotropy model as a function of the angle between porosity perturbations and the shear plane. Each column is for a different value of
${\it\alpha}_{0}$
. In all cases,
${\it\beta}_{0}=0$
,
${\it\lambda}=27$
and
${\it\phi}_{0}=0.05$
. (a–c) Growth rate
${\dot{s}}$
of perturbations
${\it\phi}_{1}$
from (5.3). The heavy line represents the total growth rate. Light lines show the growth rate decomposed into two parts: the contribution from the tensile instability (first term of (5.3)) and that from the shear instability (second term). (d–f) Compliances
$\widetilde{\unicode[STIX]{x1D60A}}_{YYYY}^{(0)}$
and
$\widetilde{\unicode[STIX]{x1D60A}}_{YYXY}^{(0)}$
in (5.3).
5.1.3. How viscous anisotropy causes lowering of band angle
Although the effect of viscous anisotropy is evident from the total growth rate shown in figure 5(a–c), it is not immediately obvious why the dominant band angle is lowered by viscous anisotropy. The physical mechanism can be understood by considering the tensile and shear components of the growth rate independently (first and second terms of (5.3), respectively). In figure 5(a–c), these two growth rates are plotted separately for various values of
${\it\alpha}_{0}$
(thin solid curve for tensile instability; thin dashed curve for shear instability). Comparison of figure 5(a–c) reveals that the peak split occurs through (i) stabilisation of the tensile instability and (ii) emergence of the shear instability with increasing magnitude of anisotropy
${\it\alpha}$
. We consider each of these in turn.
To understand why viscous anisotropy stabilises the tensile instability, we return to the systematics of the base-state stress (§ 4.1). Comparison of the three columns of figure 3 shows that, as
${\it\alpha}_{0}$
increases, the tensile stress
${\it\tau}_{YY}^{(0)}$
decreases in amplitude and becomes compressive at all angles. With
${\it\tau}_{YY}^{(0)}\leqslant 0$
, the first term in (5.3) is always less than or equal to zero, and hence stable.
To understand why viscous anisotropy destabilises the shear mechanism, we consider the coupling between the shear stress that drives the instability and the normal strain rate that is responsible for its growth. As shown by (4.9) with
${\it\tau}_{YY}^{(f)}={\it\lambda}{\it\phi}_{1}{\it\tau}_{YY}^{(0)}$
and
${\it\tau}_{XY}^{(f)}={\it\lambda}{\it\phi}_{1}{\it\tau}_{XY}^{(0)}$
, the shear stress
${\it\tau}_{XY}^{(0)}$
is coupled to the normal strain rate
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
via
$\widetilde{\unicode[STIX]{x1D60A}}_{YYXY}^{(0)}$
. The angular dependence of
$\widetilde{\unicode[STIX]{x1D60A}}_{YYXY}^{(0)}$
is shown by dashed curves in figure 5(d–f). If
$\widetilde{\unicode[STIX]{x1D60A}}_{YYXY}^{(0)}{\it\tau}_{XY}^{(0)}$
is positive, then
${\dot{s}}$
is positive (or
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
is in phase with
${\it\phi}_{1}$
) and the shear mechanism contributes to unstable growth of porosity perturbations. In fact, this product is positive (or zero) for all
${\it\theta}$
, enabling us to name this coupling ‘shear stress-induced expansion’. This non-trivial result comes from the assumed microstructural behaviour: that stress-induced softening occurs in the tensile (
${\it\tau}_{3}$
) direction, as illustrated in Takei & Katz (Reference Takei and Katz2013, figure 7a). So the porosity perturbation grows because of the shear mechanism, for which the low angle is favourable.
5.1.4. Summary of static anisotropy model
As a recap and summary, note that, under isotropic viscosity, the growth of bands at
$45^{\circ }$
to the shear plane is caused by a tensile instability (Stevenson Reference Stevenson1989; Spiegelman Reference Spiegelman2003). In contrast, under anisotropic viscosity, the peak growth rate of bands is controlled by a distinct shear instability. Although the peak growth rate of the shear instability occurs at
${\it\theta}<15^{\circ }$
, stabilisation at these low angles by the tensile mechanism acts to give a maximum in the combined growth rate at
${\it\theta}=15^{\circ }$
.
The comparison between isotropic and anisotropic systems developed above is summarised in the first two rows of table 1. The tensile instability is separated into porosity weakening,
${\it\lambda}$
, which is fundamental to all models, and the tensile stress across bands,
${\it\tau}_{YY}^{(0)}$
, which affects both isotropic and anisotropic cases. A shift of
${\it\tau}_{YY}^{(0)}$
to more negative, compressive values (represented by ▵) stabilises the tensile instability. In contrast, the difference in shear instability can be simply shown by its existence or non-existence (○ or —). It is the leading-order terms
${\it\tau}_{YY}^{(0)}$
and
$\widetilde{\unicode[STIX]{x1D60A}}_{YYXY}^{(0)}$
that are responsible for these differences.
Katz et al. (Reference Katz, Spiegelman and Holtzman2006) extended the analysis of isotropic viscosity to include a power-law dependence of viscosity on strain rate (or equivalently on stress). They showed that strain-rate weakening viscosity leads to lowering of band angle. In the discussion section, we compare the angle-lowering mechanism of viscous anisotropy to that of the power-law viscosity. This is enabled by a reanalysis of the power-law model using the rotated coordinate system.
Table 1. Summary of band formation models: ○, exists; —, does not exist; ▵, is modified.

5.2. Dynamic anisotropy
We return to the full expression for the growth rate of bands, (4.11), to develop a physical understanding of why dynamic anisotropy lowers band angles, as observed in the numerical solutions (figure 1). To do so we take
${\it\alpha}_{0}=1$
and again make the simplifying assumption that
${\it\beta}=0$
(though see § 5.3).
The perturbations in
${\it\Theta}$
and
${\it\alpha}$
under dynamic anisotropy are obtained by linearisation of (2.3) and (2.4) with respect to the stress perturbation
${\it\tau}_{ij}^{(1)}$
. The expansion is conducted around the base-state values
${\it\Theta}_{0}$
and
${\it\alpha}_{0}$
. In appendix B, we show that the sensitivity of
${\it\alpha}$
to variations in deviatoric stress is given by the parameter

This parameter allows us to write
${\it\alpha}_{1}=\breve{{\it\alpha}}{\rm\Delta}{\it\tau}^{(1)}$
, and hence to see that static anisotropy corresponds to the case where
$\breve{{\it\alpha}}=0$
. The situation for
${\it\Theta}$
is more complicated because there is no single parameter that controls its sensitivity to deviatoric stress; variations of
${\it\Theta}$
can either be fully considered or fully neglected. Fortunately, numerical and analytical results show that these variations (
${\it\Theta}_{1}$
) play an insignificant role in the understanding of band angles, and hence we need consider only the magnitude of anisotropy
${\it\alpha}$
. This is achieved by looking at the dependence of key quantities (especially
${\dot{s}}$
) on
$\breve{{\it\alpha}}$
.

Figure 6. Characteristics of the dynamic anisotropy model for various values of
$\breve{{\it\alpha}}$
as a function of the angle
${\it\theta}$
between porosity perturbations and the shear plane. In all panels,
${\it\alpha}_{0}=1$
,
${\it\beta}_{0}=0$
,
${\it\lambda}=27$
and
${\it\phi}_{0}=0.05$
. (a) Full growth rate
${\dot{s}}$
from (4.11). (b) Growth rate from (4.11) split into the tensile-instability term (solid lines) and the shear-instability term (dashed lines). (c) Band-normal forcing stress (solid lines; (4.8a
)) and band-parallel perturbation stress (dashed lines; (B 5)). Both are divided by
${\it\lambda}{\it\phi}_{1}$
. Calculation details are in appendix B. (d,e) The anisotropy perturbations
${\it\alpha}_{1}$
and
${\it\Theta}_{1}$
divided by
${\it\lambda}{\it\phi}_{1}$
, calculated according to (B 8a,b
).
The growth rate of porosity perturbations
${\dot{s}}$
is shown in figure 6(a) for
${\it\alpha}_{0}=1$
and for values of
$\breve{{\it\alpha}}$
ranging from 0 to 2. Although
$\breve{{\it\alpha}}=0$
does not exclude linearised variations in
${\it\Theta}$
, comparison with the
${\dot{s}}$
curve in figure 5(b) confirms that variations in
${\it\Theta}$
are ineffectual; with
$\breve{{\it\alpha}}=0$
the expected band angle is
$45^{\circ }$
. For increasing
$\breve{{\it\alpha}}$
, the growth-rate peak again splits into peaks at low and high angles. It is important to note that the mean value of
${\it\alpha}$
(
${\it\alpha}_{0}$
) is not changed in this exercise – only the amplitude of variations about that mean. Consistent with the numerical results of figure 1, dynamic variations in the magnitude of anisotropy can sharply reduce band angles, even at moderate
${\it\alpha}_{0}$
for which the static anisotropy model predicts a high band angle (
$45^{\circ }$
).
Figure 6(b) breaks the full growth rate into two parts, each associated with one of the terms of (4.11). Dashed lines, representing the shear instability, are almost unaffected by
$\breve{{\it\alpha}}$
. In contrast, the tensile instability is strongly stabilised with increasing
$\breve{{\it\alpha}}$
. This stabilisation causes the peak of the full growth rate in figure 6(a) to split into low- and high-angle peaks. To understand why dynamic anisotropy promotes low band angles, it is therefore sufficient to understand why it stabilises the tensile instability.
The tensile instability is driven by
${\it\tau}_{YY}^{(f)}$
, as discussed in § 5.1. This represents the normal stress (tension positive) that arises when viscosity perturbations interact with the base-state strain rate. The detailed form of the forcing stress for the dynamic anisotropy model is given in (B 4). Figure 6(c) shows that the forcing normal stress
${\it\tau}_{YY}^{(f)}$
varies significantly with
$\breve{{\it\alpha}}$
(whereas the forcing shear stress, not shown, is almost unaffected by
$\breve{{\it\alpha}}$
). The system compliances, which are leading-order quantities, are not affected by dynamic anisotropy. Therefore, it is the variation of
${\it\tau}_{YY}^{(f)}$
that is responsible for stabilisation of the tensile instability under dynamic anisotropy.
To develop a physical understanding of the detailed dependence of
${\it\tau}_{YY}^{(f)}$
on
${\it\phi}_{1}$
,
${\it\alpha}_{1}$
and
${\it\Theta}_{1}$
(B 4a
), focus attention on
${\it\theta}=45^{\circ }$
, as this is the dominant band angle when
$\breve{{\it\alpha}}=0$
. For bands at
$45^{\circ }$
, the solid curves in figure 6(c) show that the forcing stress goes from a positive perturbation (in phase with
${\it\phi}_{1}$
) to a negative perturbation (in antiphase with
${\it\phi}_{1}$
) with increasing
$\breve{{\it\alpha}}$
– hence the forcing stress
${\it\tau}_{YY}^{(f)}$
in the high-porosity bands goes from tensile to compressive. This change is due to an increase in the magnitude of anisotropy perturbation
${\it\alpha}_{1}=\breve{{\it\alpha}}{\rm\Delta}{\it\tau}^{(1)}$
, shown in figure 6(d). Since
${\it\tau}_{XY}^{(1)}={\it\tau}_{YY}^{(1)}=0$
, the deviatoric stress perturbation
${\rm\Delta}{\it\tau}^{(1)}$
is entirely due to the band-parallel normal stress perturbation
${\it\tau}_{XX}^{(1)}$
(according to (B 5)), which is shown by dashed curves in figure 6(c). Because
${\it\tau}_{XX}^{(0)}<0$
,
${\it\tau}_{XX}^{(1)}>0$
signifies a magnitude reduction of
${\it\tau}_{XX}$
in the high-porosity bands; the largest change occurs for bands at
${\it\theta}=45^{\circ }$
. As sensitivity
$\breve{{\it\alpha}}$
to deviatoric stress increases,
${\it\alpha}_{1}$
becomes more negative (figure 6
d). Negative values of
${\it\alpha}_{1}$
(in antiphase with
${\it\phi}_{1}$
) mean that high-porosity bands have lower deviatoric stress and weaker anisotropy than the low-porosity inter-band regions. This is consistent with numerical results in figure 1(d).
Figure 5(d–f) shows that
${\it\alpha}_{0}$
increases the normal compliance
$\widetilde{\unicode[STIX]{x1D60A}}_{YYYY}^{(0)}$
at angles between zero and 90
$^{\circ }$
. A negative perturbation to
${\it\alpha}_{0}$
therefore makes the high-porosity bands in this range of angles less compliant to tensile stress and the low-porosity inter-bands more compliant. Overall, then, the perturbation in anisotropy amplitude
${\it\alpha}_{1}$
tends to cancel the direct effect of the porosity perturbation
${\it\phi}_{1}$
on the normal compliance, and hence
${\it\alpha}_{1}$
works to stabilise the tensile instability.
The comparison between the static and dynamic anisotropy models developed in this section is summarised in table 1. These two models are identical at leading order but different at first order. Therefore, stabilisation of the tensile mechanism due to more compressive base-state stress (
${\it\tau}_{YY}^{(0)}$
) and destabilisation of the shear mechanism due to shear stress-induced expansion (
$\widetilde{\unicode[STIX]{x1D60A}}_{YYXY}^{(0)}$
) occur in both the static and the dynamic anisotropy models. These two cases differ, however, in that further stabilisation of the tensile mechanism occurs due to the dynamic variation of anisotropy magnitude (
${\it\alpha}_{1}$
). This effect hardens the band regions and weakens the inter-band regions under dynamic anisotropy. This additional factor (○ in table 1) significantly lowers the band angle.
It is interesting to note that dynamic perturbations to the angle of anisotropy
${\it\Theta}_{1}$
are not an important control on band angle. Figure 6(e) shows that they are not affected by
$\breve{{\it\alpha}}$
. More importantly,
${\it\Theta}_{1}$
is always zero for bands oriented at
${\it\theta}=45^{\circ }$
. This indicates that the stabilisation of the tensile instability and the lowering of band angle under dynamic anisotropy cannot be attributed to
${\it\Theta}_{1}$
. In numerical simulations (figure 1
e), the variations of
${\it\Theta}$
do not contribute to the lowering of band angle that is observed in figure 1(c), though they are well explained by the stability analysis at
${\it\theta}\simeq 10^{\circ }$
(red dashed line).
5.3. The effect of contiguity increase in the
${\it\tau}_{1}$
direction
Until now, we have neglected
${\it\beta}$
and focused on the effects of
${\it\alpha}$
, which quantifies contiguity decrease in the direction of maximum tension. Non-zero
${\it\alpha}$
represents a weakening in the
${\it\tau}_{3}$
direction that (i) reduces the magnitude of tensile stress and leads to (ii) shear strain-induced compression and (iii) shear stress-induced extension (Takei & Katz Reference Takei and Katz2013). We have shown that the tensile mechanism is stabilised around
${\it\theta}=45^{\circ }$
by the first of these and is stabilised around
${\it\theta}=0^{\circ }$
and
$90^{\circ }$
by the second; we have also shown that the shear mechanism is destabilised by (iii). The parameter
${\it\beta}$
quantifies the contiguity increase in the direction of maximum compression. Even if
${\it\alpha}$
is zero, a non-zero
${\it\beta}$
creates viscous anisotropy (see (2.2b
)), causing the couplings (ii) and (iii). However, figure 4 shows that
${\it\beta}_{0}$
does not cause the weakening (i). It is this weakening, by
${\it\alpha}$
only, that is responsible for splitting the growth-rate peak in both static and dynamic models. On this basis, we expect the effect of
${\it\beta}_{0}$
to be small. This is indeed the case: as shown below,
${\it\beta}$
alone does not cause a lowering of band angle, but it can affect the lowering by
${\it\alpha}$
.

Figure 7. The effect of
${\it\beta}$
on the growth rate of porosity perturbations for static (a,b) and dynamic (c–e) anisotropy. (a) Full growth rate
${\dot{s}}$
from (5.3) for
${\it\alpha}_{0}=2$
and various values of
${\it\beta}_{0}$
. A single curve for
${\it\alpha}_{0}=0$
and
${\it\beta}_{0}=2$
is also shown. (b) Cases with
${\it\alpha}_{0}=2$
decomposed into tensile and shear parts. Line grey scale has the same meaning as in (a); there is no curve for
${\it\alpha}_{0}=0$
. (c) Full growth rate
${\dot{s}}$
from (4.11) for
${\it\alpha}_{0}=\breve{{\it\alpha}}=1$
and various values of
${\it\beta}_{0}=\breve{{\it\beta}}$
. A single curve for
${\it\alpha}_{0}=\breve{{\it\alpha}}=0$
and
${\it\beta}_{0}=\breve{{\it\beta}}=1$
is also plotted. (d) Cases with
${\it\alpha}_{0}=\breve{{\it\alpha}}=1$
decomposed into tensile and shear parts. (e) The band-parallel normal stress perturbation
${\it\tau}_{XX}^{(1)}$
.
Figure 7(a,b) illustrates the effect of
${\it\beta}$
under static anisotropy. Figure 7(a) shows that, under static anisotropy,
${\dot{s}}$
is split into high- and low-angle peaks for any value of
${\it\beta}_{0}=0$
when
${\it\alpha}_{0}=2$
(solid curves), whereas it is peaked at
$45^{\circ }$
for any value of
${\it\beta}_{0}$
when
${\it\alpha}_{0}=0$
(
${\it\beta}_{0}=2$
shown by dash-dotted curve). For
${\it\alpha}_{0}=2$
, increasing
${\it\beta}_{0}$
causes a modest shift to more compressive
${\it\tau}_{YY}^{(0)}$
at
${\it\theta}\sim 0^{\circ }$
and
${\sim}90^{\circ }$
and a modest increase in the amplitude of shear stress
${\it\tau}_{XY}^{(0)}$
(figure 3
c,f,i). Therefore, as figure 7(b) shows,
${\it\beta}_{0}$
causes stabilisation of the tensile instability and destabilisation of the shear instability in equal measure. These two effects compensate each other and the solid growth-rate curves in figure 7(a) are thus all very similar to that for
${\it\beta}_{0}=0$
.
Figure 7(c–e) shows how
${\it\beta}$
affects dynamic anisotropy. Figure 7(c) shows that, for
${\dot{s}}$
in the dynamic anisotropy model, a two-peaked growth rate occurs for
${\it\alpha}_{0}=\breve{{\it\alpha}}=1$
and
${\it\beta}_{0}=\breve{{\it\beta}}=0$
(light grey curve) but does not for
${\it\alpha}_{0}=\breve{{\it\alpha}}=0$
and
${\it\beta}_{0}=\breve{{\it\beta}}=1$
(dash-dotted curve). In the former case of non-zero
${\it\alpha}$
with a double peak, increasing
${\it\beta}_{0}$
enhances the stabilisation of the tensile mechanism at
$45^{\circ }$
and deepens the valley between low- and high-angle peaks of
${\dot{s}}$
(figure 7
c,d). This occurs because
${\it\tau}_{XX}^{(1)}$
is enhanced by the overall strengthening effect of
${\it\beta}_{0}$
(figure 7
e). The very low band angles that emerge in the numerical simulation with dynamic anisotropy and
$r_{{\it\beta}}=1$
(figure 1) are therefore a consequence of both the dynamic effect of
${\it\alpha}$
and the enhancement by
${\it\beta}$
.
6. Summary and discussion
We have developed and analysed a model of coupled magma–mantle dynamics with anisotropic viscosity. The anisotropy is controlled by the orientation of principal stresses and the amount of deviatoric stress. The model presented here introduces small modifications on that of Takei & Katz (Reference Takei and Katz2013). In particular, the parameter
${\it\beta}$
models an increase in contiguity of grains in the direction of maximum compressive stress and the parameter
${\it\tau}_{offset}$
allows for a finer control on the magnitude of anisotropy and its sensitivity to stress (for
$r_{{\it\beta}}=\text{const.}$
). This description of viscous anisotropy is physically consistent with experiments and relatively simple, so its analysis should clarify the mechanics of rocks for which the assumptions hold. Existing experimental data, however, are not enough to quantitatively constrain all parameter values. The parameter studies performed here aim to understand the underlying physics.
It is known from previous theoretical work that anisotropic viscosity lowers the angle of emergent high-porosity bands. Numerical solutions (figure 1) compare uniform anisotropy imposed a priori with anisotropy that varies according to local conditions of stress. They show that dynamic anisotropy leads to lowering of band angle as compared with uniform anisotropy, where the mean magnitude and angle from the dynamic case are used in the static case. Moreover, dynamic anisotropy produces low-angle bands even when its mean values would not do so if applied uniformly and held constant. The physical reasons for this have not previously been clear. Indeed, the question of why anisotropic viscosity lowers band angle at all has not previously been addressed.
Static viscous anisotropy, in which viscous resistance to extension in the most tensile direction is decreased, predicts low angles of high-porosity bands for two reasons: (i) it suppresses the mode of instability in which tension causes extension across high-porosity bands; and (ii) it creates a mode of instability in which shear stress causes extension across high-porosity bands. The tensile instability has a peak perturbation growth rate in the maximum tensile direction (
${\it\theta}=45^{\circ }$
). When this instability is suppressed by static anisotropy, the peak growth rate shifts to the smaller angles that are favoured by the emergent shear instability. Although the growth of the lowest-angle bands are enhanced by the shear instability, perturbations parallel to the shear plane (
${\it\theta}=0^{\circ }$
) are stable because of the compressive stress created by the base-state flow. Therefore, a low but finite angle of high-porosity bands is predicted by this model. Allowing for an increase in contiguity and viscosity in the direction of maximum compression has counterbalancing effects that leave predicted band angles almost unchanged.
Dynamic viscous anisotropy, in which the anisotropy parameters are allowed to vary with the local orientation and magnitude of deviatoric stress, tends to further lower band angles. It does so because it suppresses the tensile instability around
${\it\theta}=45^{\circ }$
via the following dynamic effect. Lower deviatoric stress in viscously weak bands gives lower anisotropy there, which makes them less compliant to tensile stress across them. Enhanced anisotropy in the interleaved lower-porosity regions makes those regions more compliant. This effect overcompensates the compliance variations directly due to porosity weakening; it favours melt segregation from the bands into the inter-bands. Allowing for an increase in contiguity and viscosity in the direction of maximum compressional stress increases the contrast in band-parallel compressional stress (and deviatoric stress) between bands and inter-bands. This enhances the contrast in anisotropy and further suppresses the tensile instability. Dynamic anisotropy makes almost no modification to the shear instability.
The additional effects of dynamic anisotropy and the anisotropic increase of contiguity are important because they make more robust the prediction of low band angles. Under static anisotropy, the mean magnitude of anisotropy must be quite high to produce low-angle bands; moderate levels are insufficient. In contrast, under dynamic anisotropy with contiguity increase in the direction of maximum compression, moderate levels of mean anisotropy efficiently produce low-angle bands. This helps to support the hypothesis that low-angle bands in experiments are due to anisotropic viscosity because it expands the parameter space in which the theoretical predictions should hold.
These conclusions were reached by use of stability analysis in a coordinate system that is rotated with respect to the plane of simple shear; in particular, the coordinate system is aligned with the wavefronts of the harmonic perturbations. This rotation leads to simpler expressions for the growth rate of perturbations: the tensile and shear modes appear as distinct terms that are amenable to physical understanding. For this reason, our analysis represents a framework in which to test and understand the family of rheologies that potentially produce low-angle bands in shearing flows. This includes variants of isotropic and anisotropic viscosity, but also potentially of dilatational granular rheology, damage, or composite rheologies (e.g. Rudge & Bercovici Reference Rudge and Bercovici2015).
An application of the rotated coordinate system to the isotropic power-law creep model with stress exponent
$n$
(Katz et al.
Reference Katz, Spiegelman and Holtzman2006) is presented in appendix A. As with all other models considered here, porosity perturbations reduce viscosity in the bands, resulting in the enhancement of the normal and shear strain rates,
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
and
$\dot{\unicode[STIX]{x1D626}}_{XY}^{(1)}$
. In this model, however, the enhanced strain rates further reduce the viscosity, which feeds back to further enhance the strain rates. The importance of this non-Newtonian feedback relative to the porosity-weakening feedback is roughly approximated by
$n-1$
. Both normal and shear strain rates contribute to the non-Newtonian feedback; the relative importance of the shear component increases with increasing
$r_{{\it\xi}}$
. Therefore, if
$n$
and
$r_{{\it\xi}}$
are sufficiently large, shear strain rate is the key weakening factor and the growth rate reaches a maximum at a substantially lowered angle. However, in contrast to anisotropic viscosity, strain-rate weakening viscosity does not give rise to the shear instability – it merely lowers the most favourable angle for the tensile instability (comparison in table 1). Although the details differ, both models predict an important role for shear stress in the lowering mechanism; both predict a low but finite angle with localised shear strain in the higher-porosity bands.
The model of viscous anisotropy used here seems promising as an explanation for laboratory experiments on the deformation of partially molten rocks. Although its detailed form must be considered tentative, we are not aware of another theory that reproduces the low-angle bands found in experiments (Holtzman et al. Reference Holtzman, Groebner, Zimmerman, Ginsberg and Kohlstedt2003) while respecting the measured stress dependence of creep viscosity (King et al. Reference King, Zimmerman and Kohlstedt2010). Furthermore, radially inward migration of magma in experiments employing torsional deformation (King, Holtzman & Kohlstedt Reference King, Holtzman and Kohlstedt2011; Qi et al. Reference Qi, Kohlstedt, Katz and Takei2015) may be direct evidence of base-state segregation, a feature that arises naturally from viscous anisotropy (Takei & Katz Reference Takei and Katz2013) but may be impossible to reconcile with isotropic viscosity. Although the present study focuses on the angle of bands, the growth rate of bands is also affected by static and dynamic anisotropy. The growth rate is lowered by static anisotropy and further lowered by dynamic anisotropy (figures 5 and 6). This can also be discerned in the different total strain and different ranges of porosity in figure 1(a,b). Therefore, a quantitative comparison between the measured and predicted growth rate becomes important for further refining and testing the theory.
In the present theory,
${\it\alpha}$
,
${\it\beta}$
and
${\it\Theta}$
are assumed to depend on stress, based on the experimental results by Daines & Kohlstedt (Reference Daines and Kohlstedt1997) and Takei (Reference Takei2010). Although this assumption is considered to be valid at small strain, possible evolution of these parameters with increasing strain has to be investigated to model the system at large strains. Indeed, for more than 200 % strain under simple shear, Zimmerman et al. (Reference Zimmerman, Zhang, Kohlstedt and Karato1999) observed that the long axis of melt pockets is predominantly oriented at an angle of
$20^{\circ }$
from
${\it\tau}_{1}$
; this is difficult to explain by stress alone. It should be noted, however, that microstructural analysis in laboratory studies has been performed in terms of shape and orientation of melt pockets. An analysis in terms of observed contiguity is more appropriate for comparison with and incorporation into the model. Numerical simulation using dynamic anisotropy and an empirically justified evolution equation for contiguity will be important in future work.
Acknowledgements
The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Seventh Framework Programme (FP7/2007–2013)/ERC grant agreement 279925. R.F.K. visited the Earthquake Research Institute of the University of Tokyo with support from the International Research Promotion Office; he is grateful for support by the Leverhulme Trust. Numerical simulations were performed at the Advanced Research Computing facility of the University of Oxford. The authors are grateful for stimulating discussions with M. Spiegelman, D. L. Kohlstedt and C. Qi, and for helpful and encouraging reviews by S. Butler and two anonymous referees.
Appendix A. Power-law creep model by Katz et al. (2006)
The model of band formation under power-law viscosity by Katz et al. (Reference Katz, Spiegelman and Holtzman2006) is formulated by (2.1) and the viscous constitutive relations

where only six of the 16 components of the two-dimensional version are shown due to the symmetry of
$\unicode[STIX]{x1D60A}_{ijkl}$
. The normalised shear viscosity
${\it\eta}({\it\phi},\dot{\unicode[STIX]{x1D626}}_{II})$
depends on porosity and the second invariant of the strain-rate tensor,
$\dot{\unicode[STIX]{x1D626}}_{II}=\sqrt{\dot{\unicode[STIX]{x1D626}}_{ij}\dot{\unicode[STIX]{x1D626}}_{ij}/2}$
, as

(Katz et al.
Reference Katz, Spiegelman and Holtzman2006; Takei & Holtzman Reference Takei and Holtzman2009b
). Equation (A 2) represents a power-law viscosity that, to represent deformation by dislocation creep, has an exponent
$n\approx 3.5$
(e.g. Karato & Wu Reference Karato and Wu1993) (
$n=1$
corresponds to Newtonian viscosity). Under dislocation creep, the strain rate is highly sensitive to the stress because dislocation velocity and density both increase with increasing stress. Hence the model of Katz et al. (Reference Katz, Spiegelman and Holtzman2006) incorporates strain-rate weakening in addition to the porosity weakening.
Katz et al. (Reference Katz, Spiegelman and Holtzman2006) demonstrated that strain-rate weakening viscosity works to lower the band angle. Although the mechanism of this lowering is briefly discussed in their paper, further analysis of their model using the perturbation-oriented coordinate system is helpful to understand their explanation and to compare it with the mechanism of viscous anisotropy. For consistency with the foregoing development,
${\it\eta}_{L}=0$
is assumed here. We can expand (A 2) into base-state and perturbation terms as

where
${\it\eta}_{0}=(\dot{\unicode[STIX]{x1D626}}_{II}^{(0)})^{(1-n)/n}=2^{(n-1)/n}$
and we have used
$\dot{\unicode[STIX]{x1D626}}_{II}^{(0)}=1/2$
and
$\dot{\unicode[STIX]{x1D626}}_{XX}^{(1)}=0$
. Combining (A 1) and (A 3) with stress balance equations (4.6b
) and (4.6c
), we obtain




The cross-coupling terms in (A 4) arise because normal
$\dot{\unicode[STIX]{x1D626}}_{YY}$
and shear
$\dot{\unicode[STIX]{x1D626}}_{XY}$
components both affect
$\dot{\unicode[STIX]{x1D626}}_{II}$
and hence modify the viscosity (by way of an increase in dislocation density). Two feedback mechanisms are thus at work, causing growth of porosity perturbations. The first of these is a direct effect: when
${\it\lambda}>0$
, high-porosity bands are weaker by virtue of their higher porosity. The second is indirect: porosity-weakened bands have a larger strain rate that, when
$n>1$
, further weakens them through the nonlinear viscosity. The relative importance of the second mechanism to the first one increases with increasing
$n-1$
.

Figure 8. Characteristics of the power-law isotropic viscosity model as a function of the angle between porosity perturbations and the shear plane. In both panels,
$n=5$
. (a) Normal strain-rate perturbation, which shows a double peak. (b) Forcing normal stress
${\it\tau}_{YY}^{(f)}$
due to porosity weakening (solid curve), strain-rate weakening associated with
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
(dashed curve) and strain-rate weakening associated with
$\dot{\unicode[STIX]{x1D626}}_{XY}^{(1)}$
(dotted curve). The three curves in (b) sum to the curve in (a).
Solving (A 4) for
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
and
$\dot{\unicode[STIX]{x1D626}}_{XY}^{(1)}$
and using (4.6a
) and (A 1), the growth rate is

Takei & Holtzman (Reference Takei and Holtzman2009b
) obtained the identical result for
${\dot{s}}$
and showed that a single peak splits into two at large
$n$
. Here, to understand the mechanism of the split, the forcing stress associated with tension,
${\it\tau}_{YY}^{(f)}$
, is plotted in figure 8(b) for each of the three terms on the right-hand side of (A 4a
). Although the tensile forcing stress due to the porosity and normal strain-rate perturbations are maximum at
${\it\theta}=45^{\circ }$
(solid and dashed curves), that due to the enhanced shear strain rate
$\dot{\unicode[STIX]{x1D626}}_{XY}^{(1)}$
has peaks at
${\it\theta}\simeq 10^{\circ }$
and
$80^{\circ }$
(dotted curve). The sum of these three curves determines the profile of
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
in figure 8(a) and hence determines the growth rate. Figure 8(b) confirms that the weakening of viscosity by the enhanced shear strain rate
$\dot{\unicode[STIX]{x1D626}}_{XY}^{(1)}$
is the main cause of the peak split of the growth rate.
Appendix B. Calculation of the dynamic anisotropy growth rate
To solve the first-order equations of force balance (4.6), we need an expansion of the viscosity tensor into its base-state and perturbation components. Equation (2.2b
) gives
$\unicode[STIX]{x1D60A}_{ijkl}$
as a function of
${\it\phi}$
,
${\it\alpha}$
,
${\it\beta}$
and
${\it\Theta}$
. Under the dynamic anisotropy model, it is necessary to account for non-zero perturbations
${\it\alpha}_{1}$
,
${\it\beta}_{1}$
and
${\it\Theta}_{1}$
. In that case,
$\unicode[STIX]{x1D60A}_{ijkl}^{(0)}$
and
$\unicode[STIX]{x1D60A}_{ijkl}^{(1)}$
are calculated as



where
${\rm\Delta}{\it\tau}^{(0)}=2(1-({\it\alpha}_{0}-{\it\beta}_{0})/4)$
. Then from (B 1b
), the stress perturbation is written as

Using (B 3), (2.2b
), (4.3) and
${\it\Theta}_{0}={\rm\pi}/4+{\it\theta}$
, the mechanical equilibrium conditions
${\it\tau}_{XY}^{(1)}={\it\tau}_{YY}^{(1)}=0$
from (4.6) are written as






Microstructural anisotropy is determined by deviatoric stress. From the total differentials of (2.3) and (2.4a
), and from
${\it\tau}_{XY}^{(1)}={\it\tau}_{YY}^{(1)}=0$
,
${\it\alpha}_{1}$
and
${\it\Theta}_{1}$
are related to
${\it\tau}_{XX}^{(1)}$
as



Then we use the expression (B 5) for
${\it\tau}_{XX}^{(1)}$
and (B 6) to obtain



Equations (B 8) can be substituted into the stress-balance equations (B 4), giving a system in which the only first-order quantities are
${\it\phi}_{1}$
,
$\dot{\unicode[STIX]{x1D626}}_{XY}^{(1)}$
and
$\dot{\unicode[STIX]{x1D626}}_{YY}^{(1)}$
.